diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 62a2f2d7..30aecaa6 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -32,6 +32,7 @@ class Concatenation(pybamm.Symbol):
         raise NotImplementedError
 
     def get_children_domains(self, children):
+        # combine domains from children
         domain = []
         for child in children:
             child_domain = child.domain
@@ -39,6 +40,11 @@ class Concatenation(pybamm.Symbol):
                 domain += child_domain
             else:
                 raise pybamm.DomainError("""domain of children must be disjoint""")
+
+        # ensure domain is sorted according to KNOWN_DOMAINS
+        domain_dict = {d: pybamm.KNOWN_DOMAINS.index(d) for d in domain}
+        domain = sorted(domain_dict, key=domain_dict.__getitem__)
+
         # Simplify domain if concatenation spans the whole cell
         if domain == ["negative electrode", "separator", "positive electrode"]:
             domain = ["whole cell"]
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 6f6acfc6..2a6c62ba 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -61,12 +61,23 @@ class Symbol(anytree.NodeMixin):
         except TypeError:
             raise TypeError("Domain: argument domain is not iterable")
         else:
-            for d in domain:
-                assert d in pybamm.KNOWN_DOMAINS, ValueError(
+            # check that domains are all known domains
+            try:
+                indicies = [pybamm.KNOWN_DOMAINS.index(d) for d in domain]
+            except ValueError:
+                raise ValueError(
                     """domain "{}" is not in known domains ({})""".format(
-                        d, str(pybamm.KNOWN_DOMAINS)
-                    )
+                        domain, str(pybamm.KNOWN_DOMAINS))
                 )
+
+            # check that domains are sorted correctly
+            is_sorted = all(a <= b for a, b in zip(indicies, indicies[1:]))
+            if not is_sorted:
+                raise ValueError(
+                    """domain "{}" is not sorted according to known domains ({})"""
+                    .format(domain, str(pybamm.KNOWN_DOMAINS))
+                )
+
             self._domain = domain
 
     @property
@@ -79,7 +90,8 @@ class Symbol(anytree.NodeMixin):
         which would then mess with loop-checking in the anytree module
         """
         return hash(
-            (self.__class__, self.name) + tuple([child.id for child in self.children])
+            (self.__class__, self.name) +
+            tuple([child.id for child in self.children])
         )
 
     def render(self):
diff --git a/setup.py b/setup.py
index 3376d2d4..cbaa5739 100644
--- a/setup.py
+++ b/setup.py
@@ -19,6 +19,7 @@ setup(
         "numpy>=1.13",
         "scipy>=1.0",
         "pandas>=0.23",
+        "anytree>=2.4.3",
         # Note: Matplotlib is loaded for debug plots, but to ensure pints runs
         # on systems without an attached display, it should never be imported
         # outside of plot() methods.
@@ -26,7 +27,10 @@ setup(
         "matplotlib>=2.0",
     ],
     extras_require={
-        "docs": ["sphinx>=1.5"],  # For doc generation
+        "docs": [
+            "sphinx>=1.5",
+            "guzzle-sphinx-theme",
+        ],  # For doc generation
         "dev": [
             "flake8>=3",  # For code style checking
             "jupyter",  # For documentation and testing

Constrain combined domains
As suggested by @martinjrobins in #83 : 

> If you write:
> 
> a = Scalar(3, domain=["negative electrode"]
> b = Scalar(4, domain=["positive electrode"]
> c = Scalar(5, domain=["separator"]
> eq = Concatenate(a,b,c)
> the domain of eq will be ["negative electrode", "positive electrode", "separator"]. But with
> 
> a = Scalar(3, domain=["negative electrode"]
> b = Scalar(4, domain=["positive electrode"]
> c = Scalar(5, domain=["separator"]
> eq = Concatenate(a,c,b)
> the domain will be identified as ["whole cell"]
> 
> Plus in discretisation the resultant Vector is different between the two concatenations.
>
> We could either impose an ordering on combined domains, or it could be the responsibility of the model creator. The former is preferable I think (might want to open another issue for this).

I agree that imposing a constraint is preferable. Best place is probably in the `domain` setter in `Symbol`, so that it also catches domains that have been created in the wrong order. 
Should raise a `DomainError` if domains are given in the wrong order.
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index dee1112f..adfe5b30 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -298,16 +298,16 @@ class TestDiscretise(unittest.TestCase):
         a_vec = disc.scalar_to_vector(a)
         np.testing.assert_allclose(a_vec.evaluate(), expected_vector)
 
-        a_vec = disc.scalar_to_vector(a, ["whole cell", "negative electrode"])
+        a_vec = disc.scalar_to_vector(a, ["negative electrode", "whole cell"])
         expected_vector = np.concatenate(
             [
-                5 * np.ones_like(mesh["whole cell"].nodes),
                 5 * np.ones_like(mesh["negative electrode"].nodes),
+                5 * np.ones_like(mesh["whole cell"].nodes),
             ]
         )
         np.testing.assert_allclose(a_vec.evaluate(), expected_vector)
 
-        a = pybamm.Scalar(5, domain=["whole cell", "negative electrode"])
+        a = pybamm.Scalar(5, domain=["negative electrode", "whole cell"])
         a_vec = disc.scalar_to_vector(a)
         np.testing.assert_allclose(a_vec.evaluate(), expected_vector)
 
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 4769a8c0..a815b4dd 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -29,14 +29,23 @@ class TestConcatenations(unittest.TestCase):
             conc.domain,
             ["negative electrode", "separator", "positive electrode", "test"],
         )
+
         # Whole cell concatenations should simplify
         conc = pybamm.Concatenation(a, b)
         self.assertEqual(conc.domain, ["whole cell"])
+
         # Can't concatenate nodes with overlapping domains
         d = pybamm.Symbol("d", domain=["separator"])
         with self.assertRaises(pybamm.DomainError):
             pybamm.Concatenation(a, b, d)
 
+        # ensure concatenated domains are sorted correctly
+        conc = pybamm.Concatenation(c, a, b)
+        self.assertEqual(
+            conc.domain,
+            ["negative electrode", "separator", "positive electrode", "test"],
+        )
+
     def test_numpy_concatenation_vectors(self):
         # with entries
         y = np.linspace(0, 1, 15)
diff --git a/tests/test_expression_tree/test_symbol.py b/tests/test_expression_tree/test_symbol.py
index b787d7f7..0724ab60 100644
--- a/tests/test_expression_tree/test_symbol.py
+++ b/tests/test_expression_tree/test_symbol.py
@@ -14,6 +14,19 @@ class TestSymbol(unittest.TestCase):
         self.assertEqual(sym.name, "a symbol")
         self.assertEqual(str(sym), "a symbol")
 
+    def test_symbol_domains(self):
+        a = pybamm.Symbol("a", domain=pybamm.KNOWN_DOMAINS[0])
+        self.assertEqual(a.domain, [pybamm.KNOWN_DOMAINS[0]])
+        a = pybamm.Symbol("a", domain=pybamm.KNOWN_DOMAINS[:2])
+        self.assertEqual(a.domain, pybamm.KNOWN_DOMAINS[:2])
+        with self.assertRaises(TypeError):
+            a = pybamm.Symbol("a", domain=1)
+        with self.assertRaises(ValueError):
+            a = pybamm.Symbol("a", domain=["unknown domain"])
+        with self.assertRaises(ValueError):
+            a = pybamm.Symbol("a", domain=[pybamm.KNOWN_DOMAINS[1],
+                                           pybamm.KNOWN_DOMAINS[0]])
+
     def test_symbol_methods(self):
         a = pybamm.Symbol("a")
         b = pybamm.Symbol("b")

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
3.6
requirements.txt requirements-dev.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs==22.2.0
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.18.1
entrypoints==0.4
flake8==3.9.2
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.8.3
iniconfig==1.1.1
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.6.1
mistune==0.8.4
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging==21.3
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy==1.0.0
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py==1.11.0
-e git+https://github.com/tinosulzer/PyBaMM.git@0d54ab4c834bce7ca0f068904f5f9561fca4bba1#egg=pybamm
pycodestyle==2.7.0
pycparser==2.21
pyflakes==2.3.1
Pygments==2.14.0
pyparsing==3.1.4
pyrsistent==0.18.0
pytest==7.0.1
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==5.3.0
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
tomli==1.2.3
tornado==6.1
traitlets==4.3.3
typing_extensions==4.1.1
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp==3.6.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - attrs==22.2.0
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.18.1
      - entrypoints==0.4
      - flake8==3.9.2
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.8.3
      - iniconfig==1.1.1
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.6.1
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - pluggy==1.0.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - py==1.11.0
      - pycodestyle==2.7.0
      - pycparser==2.21
      - pyflakes==2.3.1
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyrsistent==0.18.0
      - pytest==7.0.1
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==5.3.0
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tomli==1.2.3
      - tornado==6.1
      - traitlets==4.3.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
      - zipp==3.6.0
prefix: /opt/conda/envs/PyBaMM

tests/test_expression_tree/test_concatenations.py::TestConcatenations::test_concatenation_domains tests/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_domains
tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_concatenation tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_concatenation_of_scalars tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_core_NotImplementedErrors tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_discretise_slicing tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_discretise_spatial_operator tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_complex_expression tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_dict tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_initial_conditions tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_model tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_symbol_base tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_scalar_to_vector tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_vector_of_ones tests/test_expression_tree/test_concatenations.py::TestConcatenations::test_base_concatenation tests/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_concatenation_vector_scalar tests/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_concatenation_vectors tests/test_expression_tree/test_symbol.py::TestSymbol::test_multiple_symbols tests/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_evaluation tests/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_init tests/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_methods tests/test_expression_tree/test_symbol.py::TestSymbol::test_symbol_repr
setup.py pybamm/expression_tree/symbol.py pybamm/expression_tree/concatenations.py
setup.py pybamm/expression_tree/symbol.py pybamm/expression_tree/concatenations.py