Address review comments

Author: pbrubeck

firedrake/interpolation.py

@@ -1013,6 +1013,9 @@ def callable():
         # Make sure we have an expression of the right length i.e. a value for
         # each component in the value shape of each function space
         loops = []
+        if numpy.prod(expr.ufl_shape, dtype=int) != V.value_size:
+            raise RuntimeError('Expression of length %d required, got length %d'
+                               % (V.value_size, numpy.prod(expr.ufl_shape, dtype=int)))
 
         if len(V) == 1:
             loops.extend(_interpolator(V, tensor, expr, subset, arguments, access, bcs=bcs))

firedrake/mesh.py

@@ -2400,22 +2400,28 @@ 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.functionspace import FunctionSpace
         from firedrake.function import Function
+        from firedrake.functionspace import FunctionSpace
         from firedrake.projection import project
 
-        mesh = self
         if self.topological_dimension() == 0:
-            P0 = FunctionSpace(mesh, "DG", 0)
+            # On vertex-only meshes we define the cell sizes as 1
+            P0 = FunctionSpace(self, "DG", 0)
             return Function(P0).assign(1)
 
-        elif self.ufl_coordinate_element().degree() > 1:
-            V = self.coordinates.function_space().reconstruct(degree=1)
-            mesh = Mesh(Function(V).interpolate(self.coordinates))
+        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
 
+        # Project the CellSize into P1
         P1 = FunctionSpace(mesh, "Lagrange", 1)
         h = project(CellSize(mesh), P1)
-        if self is not mesh:
+
+        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