Merge pull request #48 from pyroteus/45_mesh_copy

Copy mesh upon creating Mover

movement/monge_ampere.py

@@ -77,17 +77,17 @@ def __init__(self, mesh, monitor_function, **kwargs):
         super().__init__(mesh, monitor_function=monitor_function)
 
         # Create function spaces
-        self.P0 = firedrake.FunctionSpace(mesh, "DG", 0)
-        self.P1 = firedrake.FunctionSpace(mesh, "CG", 1)
-        self.P1_vec = firedrake.VectorFunctionSpace(mesh, "CG", 1)
-        self.P1_ten = firedrake.TensorFunctionSpace(mesh, "CG", 1)
+        self.P0 = firedrake.FunctionSpace(self.mesh, "DG", 0)
+        self.P1 = firedrake.FunctionSpace(self.mesh, "CG", 1)
+        self.P1_vec = firedrake.VectorFunctionSpace(self.mesh, "CG", 1)
+        self.P1_ten = firedrake.TensorFunctionSpace(self.mesh, "CG", 1)
 
         # Create objects used during the mesh movement
         self.theta = firedrake.Constant(0.0)
         self.monitor = firedrake.Function(self.P1, name="Monitor function")
         self.monitor.interpolate(self.monitor_function(self.mesh))
         self.volume = firedrake.Function(self.P0, name="Mesh volume")
-        self.volume.interpolate(ufl.CellVolume(mesh))
+        self.volume.interpolate(ufl.CellVolume(self.mesh))
         self.original_volume = firedrake.Function(self.volume)
         self.total_volume = firedrake.assemble(firedrake.Constant(1.0) * self.dx)
         self.L_P0 = firedrake.TestFunction(self.P0) * self.monitor * self.dx

movement/mover.py

@@ -12,28 +12,30 @@ class PrimeMover(object):
     """
 
     def __init__(self, mesh, monitor_function=None, **kwargs):
-        self.mesh = mesh
+        self.mesh = firedrake.Mesh(mesh.coordinates.copy(deepcopy=True))
         self.monitor_function = monitor_function
-        self.dim = mesh.topological_dimension()
+        self.dim = self.mesh.topological_dimension()
         if self.dim != 2:
             raise NotImplementedError(
                 f"Dimension {self.dim} has not been considered yet"
             )
-        self.gdim = mesh.geometric_dimension()
+        self.gdim = self.mesh.geometric_dimension()
         self.plex = self.mesh.topology_dm
         self.vertex_indices = self.plex.getDepthStratum(0)
         self.edge_indices = self.plex.getDepthStratum(1)
 
         # Measures
         degree = kwargs.get("quadrature_degree")
-        self.dx = firedrake.dx(domain=mesh, degree=degree)
-        self.ds = firedrake.ds(domain=mesh, degree=degree)
-        self.dS = firedrake.dS(domain=mesh, degree=degree)
+        self.dx = firedrake.dx(domain=self.mesh, degree=degree)
+        self.ds = firedrake.ds(domain=self.mesh, degree=degree)
+        self.dS = firedrake.dS(domain=self.mesh, degree=degree)
 
         # Mesh coordinate functions
-        self.coord_space = mesh.coordinates.function_space()
-        self._x = firedrake.Function(mesh.coordinates, name="Physical coordinates")
-        self.xi = firedrake.Function(mesh.coordinates, name="Computational coordinates")
+        self.coord_space = self.mesh.coordinates.function_space()
+        self._x = firedrake.Function(self.mesh.coordinates, name="Physical coordinates")
+        self.xi = firedrake.Function(
+            self.mesh.coordinates, name="Computational coordinates"
+        )
         self.v = firedrake.Function(self.coord_space, name="Mesh velocity")
 
     def _get_coordinate_section(self):

movement/spring.py

@@ -39,10 +39,10 @@ def __init__(self, mesh, **kwargs):
         :arg mesh: the physical mesh
         """
         super().__init__(mesh)
-        self.HDivTrace = firedrake.FunctionSpace(mesh, "HDiv Trace", 0)
-        self.HDivTrace_vec = firedrake.VectorFunctionSpace(mesh, "HDiv Trace", 0)
+        self.HDivTrace = firedrake.FunctionSpace(self.mesh, "HDiv Trace", 0)
+        self.HDivTrace_vec = firedrake.VectorFunctionSpace(self.mesh, "HDiv Trace", 0)
         self.f = firedrake.Function(self.mesh.coordinates.function_space())
-        self.displacement = np.zeros(mesh.num_vertices())
+        self.displacement = np.zeros(self.mesh.num_vertices())
 
     @property
     @PETSc.Log.EventDecorator("SpringMover_Base.facet_areas")

test/test_monge_ampere.py

@@ -26,7 +26,7 @@ def test_uniform_monitor(method, exports=False):
     mover = MongeAmpereMover(mesh, const_monitor, method=method)
     num_iterations = mover.move()
 
-    assert np.allclose(coords, mesh.coordinates.dat.data)
+    assert np.allclose(coords, mover.mesh.coordinates.dat.data)
     assert num_iterations == 0
 
 
@@ -42,14 +42,11 @@ def test_continue(method, exports=False):
     # Solve the problem to a weak tolerance
     mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=0.1)
     num_it_init = mover.move()
-    phi, sigma = mover.phi, mover.sigma
     if exports:
-        File("outputs/init.pvd").write(phi, sigma)
+        File("outputs/init.pvd").write(mover.phi, mover.sigma)
 
-    # Continue with a new Mover
-    mover = MongeAmpereMover(
-        mesh, ring_monitor, method=method, rtol=rtol, phi_init=phi, sigma_init=sigma
-    )
+    # Continue with a tighter tolerance
+    mover.rtol = rtol
     num_it_continue = mover.move()
     if exports:
         File("outputs/continue.pvd").write(mover.phi, mover.sigma)
@@ -83,14 +80,14 @@ def test_change_monitor(method, exports=False):
     mover.move()
     if exports:
         File("outputs/ring.pvd").write(mover.phi, mover.sigma)
-    assert not np.allclose(coords, mesh.coordinates.dat.data, atol=tol)
+    assert not np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol)
 
     # Adapt to a constant monitor
     mover.monitor_function = const_monitor
     mover.move()
     if exports:
         File("outputs/const.pvd").write(mover.phi, mover.sigma)
-    assert np.allclose(coords, mesh.coordinates.dat.data, atol=tol)
+    assert np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol)
 
 
 @pytest.mark.slow
@@ -113,10 +110,10 @@ def test_bcs(method, fix_boundary):
     mover.move()
 
     # Check boundary lengths are preserved
-    bnd_new = assemble(one * ds(domain=mesh))
+    bnd_new = assemble(one * ds(domain=mover.mesh))
     assert np.isclose(bnd, bnd_new)
 
     # Check boundaries are indeed fixed
     if fix_boundary:
-        bnd_coords_new = mesh.coordinates.dat.data[bnodes]
+        bnd_coords_new = mover.mesh.coordinates.dat.data[bnodes]
         assert np.allclose(bnd_coords, bnd_coords_new)