From 9506b9046c8b2fdbda1f029c06e72c0441bc5148 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Apr 2025 22:16:05 +0100 Subject: [PATCH 1/6] Interpolate physically-mapped elements into VertexOnlyMesh --- firedrake/assemble.py | 2 +- firedrake/interpolation.py | 2 +- firedrake/mesh.py | 18 +++++++++++++++-- tests/firedrake/meshes/test_meshes_volume.py | 9 +++++++++ .../regression/test_interpolate_zany.py | 16 ++++++++++++++- tsfc/driver.py | 6 ++++-- tsfc/fem.py | 3 ++- tsfc/kernel_interface/firedrake_loopy.py | 20 +++++++++---------- 8 files changed, 58 insertions(+), 18 deletions(-) diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 50287cc2f2..97104b05cc 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1114,7 +1114,7 @@ class ZeroFormAssembler(ParloopFormAssembler): Parameters ---------- - form : ufl.Form or slate.TensorBasehe + form : ufl.Form or slate.TensorBase 0-form. Notes diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 4222e9e325..05d3438318 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -1194,7 +1194,7 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None): co = target_mesh.cell_orientations() parloop_args.append(co.dat(op2.READ, co.cell_node_map())) if needs_cell_sizes: - cs = target_mesh.cell_sizes + cs = source_mesh.cell_sizes parloop_args.append(cs.dat(op2.READ, cs.cell_node_map())) for coefficient in coefficients: if isinstance(target_mesh.topology, firedrake.mesh.VertexOnlyMeshTopology): diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 0be58d6d95..6b01f1e14b 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2401,9 +2401,23 @@ 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.projection import project - P1 = FunctionSpace(self, "Lagrange", 1) - return project(CellSize(self), P1) + + mesh = self + if self.topological_dimension() == 0: + P0 = FunctionSpace(mesh, "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)) + + P1 = FunctionSpace(mesh, "Lagrange", 1) + h = project(CellSize(mesh), P1) + if self is not 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/meshes/test_meshes_volume.py b/tests/firedrake/meshes/test_meshes_volume.py index d355ca120e..af88e7a401 100644 --- a/tests/firedrake/meshes/test_meshes_volume.py +++ b/tests/firedrake/meshes/test_meshes_volume.py @@ -20,3 +20,12 @@ def test_meshes_volume_solidtorusmesh(): vol = assemble(Constant(1., domain=mesh) * dx) exact = (pi * r * r) * (2 * pi * R) assert abs(vol - exact) / exact < .0005 + + +def test_meshes_cell_sizes(): + msh = UnitDiskMesh(4) + msh.cell_sizes + + V = msh.coordinates.function_space().reconstruct(degree=2) + msh = Mesh(Function(V).interpolate(msh.coordinates)) + msh.cell_sizes diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index f814bf604d..58a56d0cdf 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,17 @@ def test_interpolate(V, mesh, which, expect, tolerance): g.interpolate(a + b) assert numpy.allclose(norm(g - expect), 0, atol=tolerance) + + +@pytest.fixture +def vom(mesh): + return VertexOnlyMesh(mesh, [(0.5, 0.5)]) + + +def test_interpolate_zany_into_vom(V, vom): + P0 = FunctionSpace(vom, "DG", 0) + vvom = TestFunction(P0) + Fvom = inner(1, vvom) *dx + + v = TestFunction(V) + assemble(Interpolate(v, Fvom)) diff --git a/tsfc/driver.py b/tsfc/driver.py index 9e518292bc..32a65f0542 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -208,11 +208,13 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, # Collect required coefficients and determine numbering coefficients = extract_coefficients(expression) orig_coefficients = extract_coefficients(orig_expression) - coefficient_numbers = tuple(orig_coefficients.index(c) for c in coefficients) + coefficient_numbers = tuple(map(orig_coefficients.index, coefficients)) builder.set_coefficient_numbers(coefficient_numbers) + elements = [f.ufl_element() for f in (*coefficients, *arguments)] + needs_external_coords = False - if has_type(expression, GeometricQuantity) or any(fem.needs_coordinate_mapping(c.ufl_element()) for c in coefficients): + if has_type(expression, GeometricQuantity) or any(map(fem.needs_coordinate_mapping, elements)): # Create a fake coordinate coefficient for a domain. coords_coefficient = ufl.Coefficient(ufl.FunctionSpace(domain, domain.ufl_coordinate_element())) builder.domain_coordinate[domain] = coords_coefficient diff --git a/tsfc/fem.py b/tsfc/fem.py index d86be550db..421f58906c 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -336,7 +336,8 @@ class GemPointContext(ContextBase): def basis_evaluation(self, finat_element, mt, entity_id): return finat_element.point_evaluation(mt.local_derivatives, self.point_expr, - (self.integration_dim, entity_id)) + (self.integration_dim, entity_id), + CoordinateMapping(mt, self)) class Translator(MultiFunction, ModifiedTerminalMixin, ufl2gem.Mixin): diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py index 6481e40bce..c767fc976e 100644 --- a/tsfc/kernel_interface/firedrake_loopy.py +++ b/tsfc/kernel_interface/firedrake_loopy.py @@ -104,6 +104,16 @@ def _coefficient(self, coefficient, name): self.coefficient_map[coefficient] = expr return expr + def set_coordinates(self, domain): + """Prepare the coordinate field. + + :arg domain: :class:`ufl.Domain` + """ + # Create a fake coordinate coefficient for a domain. + f = Coefficient(FunctionSpace(domain, domain.ufl_coordinate_element())) + self.domain_coordinate[domain] = f + self._coefficient(f, "coords") + def set_cell_sizes(self, domain): """Setup a fake coefficient for "cell sizes". @@ -304,16 +314,6 @@ def set_arguments(self, arguments): self.return_variables = return_variables self.argument_multiindices = argument_multiindices - def set_coordinates(self, domain): - """Prepare the coordinate field. - - :arg domain: :class:`ufl.Domain` - """ - # Create a fake coordinate coefficient for a domain. - f = Coefficient(FunctionSpace(domain, domain.ufl_coordinate_element())) - self.domain_coordinate[domain] = f - self._coefficient(f, "coords") - def set_coefficients(self, integral_data, form_data): """Prepare the coefficients of the form. From 4849fcd1d43c0570071c59ee40083b8f99e02cf5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Apr 2025 22:24:34 +0100 Subject: [PATCH 2/6] lint --- tests/firedrake/regression/test_interpolate_zany.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index 58a56d0cdf..9e971ad681 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/tests/firedrake/regression/test_interpolate_zany.py @@ -81,7 +81,7 @@ def vom(mesh): def test_interpolate_zany_into_vom(V, vom): P0 = FunctionSpace(vom, "DG", 0) vvom = TestFunction(P0) - Fvom = inner(1, vvom) *dx + Fvom = inner(1, vvom) * dx v = TestFunction(V) assemble(Interpolate(v, Fvom)) From de5f5ee5f9843ba1d9c96d31c1238b82f3482c0a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Apr 2025 10:52:02 +0100 Subject: [PATCH 3/6] More extensive tests --- firedrake/interpolation.py | 3 - tests/firedrake/meshes/test_meshes_volume.py | 9 --- .../regression/test_interpolate_zany.py | 58 +++++++++++++++++-- 3 files changed, 52 insertions(+), 18 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 05d3438318..24629cf52e 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -1013,9 +1013,6 @@ 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)) diff --git a/tests/firedrake/meshes/test_meshes_volume.py b/tests/firedrake/meshes/test_meshes_volume.py index af88e7a401..d355ca120e 100644 --- a/tests/firedrake/meshes/test_meshes_volume.py +++ b/tests/firedrake/meshes/test_meshes_volume.py @@ -20,12 +20,3 @@ def test_meshes_volume_solidtorusmesh(): vol = assemble(Constant(1., domain=mesh) * dx) exact = (pi * r * r) * (2 * pi * R) assert abs(vol - exact) / exact < .0005 - - -def test_meshes_cell_sizes(): - msh = UnitDiskMesh(4) - msh.cell_sizes - - V = msh.coordinates.function_space().reconstruct(degree=2) - msh = Mesh(Function(V).interpolate(msh.coordinates)) - msh.cell_sizes diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index 9e971ad681..39dbf58377 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/tests/firedrake/regression/test_interpolate_zany.py @@ -75,13 +75,59 @@ def test_interpolate_zany_into_cg(V, mesh, which, expect, tolerance): @pytest.fixture def vom(mesh): - return VertexOnlyMesh(mesh, [(0.5, 0.5)]) + return VertexOnlyMesh(mesh, [(0.5, 0.5), (0.31, 0.72)]) -def test_interpolate_zany_into_vom(V, vom): +def test_interpolate_zany_into_vom(V, mesh, which, 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) P0 = FunctionSpace(vom, "DG", 0) - vvom = TestFunction(P0) - Fvom = inner(1, vvom) * dx + 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()) + 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) + + # 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) + + # Construct a Cofunction on P0(vom)* + Fvom = Cofunction(P0.dual()).assign(1) + expected = assemble(action(Fvom, expected)) + + # Interpolate a Function into Fvom + f_at_vom = assemble(Interpolate(fexpr, Fvom)) + assert numpy.allclose(f_at_vom, expected) + + # Interpolate a TestFunction into Fvom + expr_vom = assemble(Interpolate(vexpr, Fvom)) + f_at_vom = assemble(action(expr_vom, f)) + assert numpy.allclose(f_at_vom, expected) + + +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 - v = TestFunction(V) - assemble(Interpolate(v, Fvom)) + assert numpy.allclose(h1.dat.data, h2.dat.data) From cefb39808b65e652b939a0088243aaa229e12d5d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Apr 2025 11:13:44 +0100 Subject: [PATCH 4/6] Address review comments --- firedrake/interpolation.py | 3 +++ firedrake/mesh.py | 20 +++++++++++++------- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 24629cf52e..05d3438318 100644 --- a/firedrake/interpolation.py +++ b/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)) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 6b01f1e14b..fecbbcf344 100644 --- a/firedrake/mesh.py +++ b/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 From 635b95ec4efa0951bbff71393e9fe0f3e3b714c0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Apr 2025 12:30:33 +0100 Subject: [PATCH 5/6] Apply suggestions from code review Co-authored-by: Connor Ward --- firedrake/mesh.py | 1 - tests/firedrake/regression/test_interpolate_zany.py | 6 +++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index fecbbcf344..93e2cd2f36 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2416,7 +2416,6 @@ def cell_sizes(self): else: mesh = self - # Project the CellSize into P1 P1 = FunctionSpace(mesh, "Lagrange", 1) h = project(CellSize(mesh), P1) diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index 39dbf58377..f2a8c3b90b 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/tests/firedrake/regression/test_interpolate_zany.py @@ -110,16 +110,16 @@ def test_interpolate_zany_into_vom(V, mesh, which, vom): # Construct a Cofunction on P0(vom)* Fvom = Cofunction(P0.dual()).assign(1) - expected = assemble(action(Fvom, expected)) + expected_action = assemble(action(Fvom, expected)) # Interpolate a Function into Fvom f_at_vom = assemble(Interpolate(fexpr, Fvom)) - assert numpy.allclose(f_at_vom, expected) + assert numpy.allclose(f_at_vom, expected_action) # Interpolate a TestFunction into Fvom expr_vom = assemble(Interpolate(vexpr, Fvom)) f_at_vom = assemble(action(expr_vom, f)) - assert numpy.allclose(f_at_vom, expected) + assert numpy.allclose(f_at_vom, expected_action) def test_high_order_mesh_cell_sizes(): From d33efa9ca7c51ef107caa4bad5ead5c809e2d48b Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 9 Apr 2025 21:44:05 +0100 Subject: [PATCH 6/6] cleanup --- firedrake/mesh.py | 23 +-------- .../regression/test_interpolate_zany.py | 51 ++++++++++--------- 2 files changed, 28 insertions(+), 46 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 93e2cd2f36..0be58d6d95 100644 --- a/firedrake/mesh.py +++ b/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. diff --git a/tests/firedrake/regression/test_interpolate_zany.py b/tests/firedrake/regression/test_interpolate_zany.py index f2a8c3b90b..b2054843cd 100644 --- a/tests/firedrake/regression/test_interpolate_zany.py +++ b/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)