Fix for MaximumSurface and MinimumSurface

festim/exports/derived_quantities/derived_quantities.py

@@ -1,6 +1,8 @@
 from festim import (
     MinimumVolume,
     MaximumVolume,
+    MinimumSurface,
+    MaximumSurface,
     DerivedQuantity,
 )
 import fenics as f
@@ -106,6 +108,7 @@ def make_header(self):
 
     def assign_measures_to_quantities(self, dx, ds):
         self.volume_markers = dx.subdomain_data()
+        self.surface_markers = ds.subdomain_data()
         for quantity in self:
             quantity.dx = dx
             quantity.ds = ds
@@ -129,6 +132,8 @@ def compute(self, t):
         for quantity in self:
             if isinstance(quantity, (MaximumVolume, MinimumVolume)):
                 value = quantity.compute(self.volume_markers)
+            elif isinstance(quantity, (MaximumSurface, MinimumSurface)):
+                value = quantity.compute(self.surface_markers)
             else:
                 value = quantity.compute()
 

festim/exports/derived_quantities/maximum_surface.py

@@ -1,5 +1,5 @@
 from festim import SurfaceQuantity
-import fenics as f
+from dolfin import MPI
 import numpy as np
 
 
@@ -46,16 +46,22 @@ def title(self):
     def compute(self, surface_markers):
         """Maximum of f over subdomains facets marked with self.surface"""
         V = self.function.function_space()
-
+        mesh = surface_markers.mesh()
         dm = V.dofmap()
 
-        subd_dofs = np.unique(
-            np.hstack(
-                [
-                    dm.cell_dofs(c.index())
-                    for c in f.SubsetIterator(surface_markers, self.surface)
-                ]
-            )
+        facets = surface_markers.where_equal(self.surface)
+        entity_closure_dofs = np.array(
+            dm.entity_closure_dofs(mesh, mesh.topology().dim() - 1, facets),
+            dtype=np.int32,
         )
+        local_dofs = entity_closure_dofs < len(self.function.vector().get_local())
+        if len(local_dofs) == 0:
+            local_max = -np.inf
+        else:
+            local_max = np.max(
+                self.function.vector().get_local()[entity_closure_dofs[local_dofs]]
+            )
+
+        global_max = MPI.max(mesh.mpi_comm(), local_max)
 
-        return np.max(self.function.vector().get_local()[subd_dofs])
+        return global_max

festim/exports/derived_quantities/minimum_surface.py

@@ -1,5 +1,5 @@
 from festim import SurfaceQuantity
-import fenics as f
+from dolfin import MPI
 import numpy as np
 
 
@@ -46,16 +46,22 @@ def title(self):
     def compute(self, surface_markers):
         """Minimum of f over subdomains facets marked with self.surface"""
         V = self.function.function_space()
-
+        mesh = surface_markers.mesh()
         dm = V.dofmap()
 
-        subd_dofs = np.unique(
-            np.hstack(
-                [
-                    dm.cell_dofs(c.index())
-                    for c in f.SubsetIterator(surface_markers, self.surface)
-                ]
-            )
+        facets = surface_markers.where_equal(self.surface)
+        entity_closure_dofs = np.array(
+            dm.entity_closure_dofs(mesh, mesh.topology().dim() - 1, facets),
+            dtype=np.int32,
         )
+        local_dofs = entity_closure_dofs < len(self.function.vector().get_local())
+        if len(local_dofs) == 0:
+            local_min = np.inf
+        else:
+            local_min = np.min(
+                self.function.vector().get_local()[entity_closure_dofs[local_dofs]]
+            )
+
+        global_min = MPI.max(mesh.mpi_comm(), local_min)
 
-        return np.min(self.function.vector().get_local()[subd_dofs])
+        return global_min

test/system/test_misc.py

@@ -905,3 +905,46 @@ def test_error_with_multiple_1d_domains_no_borders():
         match="borders attributes need to be set for multiple 1D domains",
     ):
         my_model.initialise()
+
+
+@pytest.mark.parametrize("mesh", [f.UnitIntervalMesh(10), f.UnitSquareMesh(10, 10)])
+def test_error_surface_quantities(mesh):
+    """Test to catch #983"""
+
+    model = F.Simulation()
+
+    model.materials = F.Material(id=1, D_0=1, E_D=0)
+
+    volume_markers = f.MeshFunction("size_t", mesh, mesh.topology().dim())
+    surface_markers = f.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
+
+    surface_markers.set_all(0)
+
+    left = f.CompiledSubDomain("near(x[0], 0) && on_boundary")
+    left.mark(surface_markers, 1)
+    right = f.CompiledSubDomain("near(x[0], 1) && on_boundary")
+    right.mark(surface_markers, 2)
+
+    model.mesh = F.Mesh(
+        mesh, volume_markers=volume_markers, surface_markers=surface_markers
+    )
+
+    model.materials = F.Material(id=1, D_0=1, E_D=1)
+
+    model.T = F.x + F.y * (mesh.topology().dim() - 1)
+
+    model.exports = F.DerivedQuantities(
+        [
+            F.MinimumSurface(field="T", surface=1),
+            F.MaximumSurface(field="T", surface=2),
+            F.MinimumSurface(field="T", surface=3),
+            F.MaximumSurface(field="T", surface=4),
+        ]
+    )
+
+    model.settings = F.Settings(
+        transient=False, absolute_tolerance=1e10, relative_tolerance=1e-10
+    )
+
+    model.initialise()
+    model.run()

test/unit/test_exports/test_derived_quantities/test_maximum_surface.py

@@ -2,6 +2,7 @@
 import fenics as f
 import numpy as np
 import pytest
+from dolfin import MPI
 
 
 @pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
@@ -19,17 +20,18 @@ def test_title(field, surface):
     assert my_max.title == f"Maximum {field} surface {surface} ({my_max.export_unit})"
 
 
-class TestCompute:
+def test_compute_maximum():
     """Test that the maximum surface export computes the correct value"""
 
-    mesh = f.UnitIntervalMesh(10)
+    mesh = f.UnitSquareMesh(10, 10)
     V = f.FunctionSpace(mesh, "P", 1)
 
-    c = f.interpolate(f.Expression("x[0]", degree=1), V)
+    c = f.interpolate(f.Expression("x[0] + x[1]", degree=1), V)
 
-    left = f.CompiledSubDomain("x[0] < 0.5")
-    surface_markers = f.MeshFunction("size_t", mesh, 1, 1)
-    left.mark(surface_markers, 2)
+    surface_markers = f.MeshFunction("size_t", mesh, 1, 0)
+
+    right = f.CompiledSubDomain("near(x[0], 1) && on_boundary")
+    right.mark(surface_markers, 1)
 
     dx = f.Measure("dx", domain=mesh, subdomain_data=surface_markers)
 
@@ -38,17 +40,8 @@ class TestCompute:
     my_max.function = c
     my_max.dx = dx
 
-    def test_minimum(self):
-        dm = self.V.dofmap()
-        subd_dofs = np.unique(
-            np.hstack(
-                [
-                    dm.cell_dofs(c.index())
-                    for c in f.SubsetIterator(self.surface_markers, self.surface)
-                ]
-            )
-        )
-        expected = np.max(self.c.vector().get_local()[subd_dofs])
-
-        produced = self.my_max.compute(self.surface_markers)
-        assert produced == expected
+    produced = my_max.compute(surface_markers)
+
+    expected = 2.0
+
+    assert produced == expected

test/unit/test_exports/test_derived_quantities/test_minimum_surface.py

@@ -2,6 +2,7 @@
 import fenics as f
 import numpy as np
 import pytest
+from dolfin import MPI
 
 
 @pytest.mark.parametrize("field,surface", [("solute", 1), ("T", 2)])
@@ -19,17 +20,18 @@ def test_title(field, surface):
     assert my_min.title == f"Minimum {field} surface {surface} ({my_min.export_unit})"
 
 
-class TestCompute:
+def test_compute_minimum():
     """Test that the minimum surface export computes the correct value"""
 
-    mesh = f.UnitIntervalMesh(10)
+    mesh = f.UnitSquareMesh(10, 10)
     V = f.FunctionSpace(mesh, "P", 1)
 
-    c = f.interpolate(f.Expression("x[0]", degree=1), V)
+    c = f.interpolate(f.Expression("x[0] + x[1]", degree=1), V)
 
-    left = f.CompiledSubDomain("x[0] < 0.5")
-    surface_markers = f.MeshFunction("size_t", mesh, 1, 1)
-    left.mark(surface_markers, 2)
+    surface_markers = f.MeshFunction("size_t", mesh, 1, 0)
+
+    right = f.CompiledSubDomain("near(x[0], 1) && on_boundary")
+    right.mark(surface_markers, 1)
 
     dx = f.Measure("dx", domain=mesh, subdomain_data=surface_markers)
 
@@ -38,17 +40,8 @@ class TestCompute:
     my_min.function = c
     my_min.dx = dx
 
-    def test_minimum(self):
-        dm = self.V.dofmap()
-        subd_dofs = np.unique(
-            np.hstack(
-                [
-                    dm.cell_dofs(c.index())
-                    for c in f.SubsetIterator(self.surface_markers, self.surface)
-                ]
-            )
-        )
-        expected = np.min(self.c.vector().get_local()[subd_dofs])
-
-        produced = self.my_min.compute(self.surface_markers)
-        assert produced == expected
+    produced = my_min.compute(surface_markers)
+
+    expected = 1.0
+
+    assert produced == expected