Mesh: extend P1 cell_sizes to high-order meshes

firedrake/mesh.py

@@ -2399,10 +2399,29 @@ 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
-        P1 = FunctionSpace(self, "Lagrange", 1)
-        return project(CellSize(self), P1)
+
+        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
 
     def clear_cell_sizes(self):
         """Reset the :attr:`cell_sizes` field on this mesh geometry.

tests/firedrake/regression/test_interpolate_zany.py

@@ -51,7 +51,7 @@ def expect(V, which):
         return a + b
 
 
-def test_interpolate(V, mesh, which, expect, tolerance):
+def test_interpolate_zany_into_cg(V, mesh, which, expect, tolerance):
     degree = V.ufl_element().degree()
     Vcg = FunctionSpace(mesh, "P", degree)
 
@@ -71,3 +71,14 @@ def test_interpolate(V, mesh, which, expect, tolerance):
         g.interpolate(a + b)
 
     assert numpy.allclose(norm(g - expect), 0, atol=tolerance)
+
+
+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)