cleanup

Author: pbrubeck

firedrake/mesh.py

@@ -2400,29 +2400,10 @@ def cell_sizes(self):
 
         This is computed by the :math:`L^2` projection of the local mesh element size."""
         from firedrake.ufl_expr import CellSize
-        from firedrake.function import Function
         from firedrake.functionspace import FunctionSpace
         from firedrake.projection import project
-
-        if self.topological_dimension() == 0:
-            # On vertex-only meshes we define the cell sizes as 1
-            P0 = FunctionSpace(self, "DG", 0)
-            return Function(P0).assign(1)
-
-        if self.ufl_coordinate_element().degree() > 1:
-            # We need a P1 mesh, as CellSize is not defined on high-order meshes
-            VectorP1 = self.coordinates.function_space().reconstruct(degree=1)
-            mesh = Mesh(Function(VectorP1).interpolate(self.coordinates))
-        else:
-            mesh = self
-
-        P1 = FunctionSpace(mesh, "Lagrange", 1)
-        h = project(CellSize(mesh), P1)
-
-        if P1.mesh() is not self:
-            # Transfer the Function on the P1 mesh into the original mesh
-            h = Function(P1.reconstruct(mesh=self), val=h.dat)
-        return h
+        P1 = FunctionSpace(self, "Lagrange", 1)
+        return project(CellSize(self), P1)
 
     def clear_cell_sizes(self):
         """Reset the :attr:`cell_sizes` field on this mesh geometry.

tests/firedrake/regression/test_interpolate_zany.py

@@ -78,39 +78,51 @@ def vom(mesh):
     return VertexOnlyMesh(mesh, [(0.5, 0.5), (0.31, 0.72)])
 
 
-def test_interpolate_zany_into_vom(V, mesh, which, vom):
+@pytest.fixture
+def expr_at_vom(V, which, vom):
+    mesh = V.mesh()
     degree = V.ufl_element().degree()
     x, y = SpatialCoordinate(mesh)
     expr = (x + y)**degree
 
-    f = Function(V)
-    f.project(expr, solver_parameters={"ksp_type": "preonly",
-                                       "pc_type": "lu"})
-    fexpr = f
-    vexpr = TestFunction(V)
-    P0 = FunctionSpace(vom, "DG", 0)
     if which == "coefficient":
         P0 = FunctionSpace(vom, "DG", 0)
     elif which == "grad":
-        fexpr = grad(fexpr)
-        vexpr = grad(vexpr)
         expr = ufl.algorithms.expand_derivatives(grad(expr))
         P0 = VectorFunctionSpace(vom, "DG", 0)
 
-    expected = Function(P0)
-    point = Constant([0]*mesh.geometric_dimension())
+    fvom = Function(P0)
+    point = Constant([0] * mesh.geometric_dimension())
     expr_at_pt = ufl.replace(expr, {SpatialCoordinate(mesh): point})
     for i, pt in enumerate(vom.coordinates.dat.data_ro):
         point.assign(pt)
-        expected.dat.data[i] = numpy.asarray(expr_at_pt, dtype=float)
+        fvom.dat.data[i] = numpy.asarray(expr_at_pt, dtype=float)
+    return fvom
+
+
+def test_interpolate_zany_into_vom(V, mesh, which, expr_at_vom):
+    degree = V.ufl_element().degree()
+    x, y = SpatialCoordinate(mesh)
+    expr = (x + y)**degree
+
+    f = Function(V)
+    f.project(expr, solver_parameters={"ksp_type": "preonly",
+                                       "pc_type": "lu"})
+    fexpr = f
+    vexpr = TestFunction(V)
+    if which == "grad":
+        fexpr = grad(fexpr)
+        vexpr = grad(vexpr)
+
+    P0 = expr_at_vom.function_space()
 
     # Interpolate a Function into P0(vom)
     f_at_vom = assemble(Interpolate(fexpr, P0))
-    assert numpy.allclose(f_at_vom.dat.data_ro, expected.dat.data_ro)
+    assert numpy.allclose(f_at_vom.dat.data_ro, expr_at_vom.dat.data_ro)
 
     # Construct a Cofunction on P0(vom)*
     Fvom = Cofunction(P0.dual()).assign(1)
-    expected_action = assemble(action(Fvom, expected))
+    expected_action = assemble(action(Fvom, expr_at_vom))
 
     # Interpolate a Function into Fvom
     f_at_vom = assemble(Interpolate(fexpr, Fvom))
@@ -120,14 +132,3 @@ def test_interpolate_zany_into_vom(V, mesh, which, vom):
     expr_vom = assemble(Interpolate(vexpr, Fvom))
     f_at_vom = assemble(action(expr_vom, f))
     assert numpy.allclose(f_at_vom, expected_action)
-
-
-def test_high_order_mesh_cell_sizes():
-    msh1 = UnitSquareMesh(2, 2)
-    h1 = msh1.cell_sizes
-
-    P2 = msh1.coordinates.function_space().reconstruct(degree=2)
-    msh2 = Mesh(Function(P2).interpolate(msh1.coordinates))
-    h2 = msh2.cell_sizes
-
-    assert numpy.allclose(h1.dat.data, h2.dat.data)