Repository: tinosulzer/PyBaMM

Issue:
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.

Discussion:

Patch:
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


Test patch:
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")