From 2f7a2035c8bd2c76d2406bcbd7ddb157845cc680 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Apr 2025 21:24:22 +0100 Subject: [PATCH] Mesh: extend P1 cell_sizes to high-order meshes --- firedrake/mesh.py | 23 +++++++++++++++++-- .../regression/test_interpolate_zany.py | 13 ++++++++++- 2 files changed, 33 insertions(+), 3 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 558834c5b6..b9e34d2e22 100644 --- a/firedrake/mesh.py +++ b/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. diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index f814bf604d..37e44eb071 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/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)