From b5698d3bfca9ae3e747b4a1c440dee5f9c257621 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Fri, 20 Feb 2026 20:02:13 +0000 Subject: [PATCH 1/3] Make it possible to reconstruct mesh with new coordinate element degree. --- src/io4dolfinx/__init__.py | 2 ++ src/io4dolfinx/utils.py | 62 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+) diff --git a/src/io4dolfinx/__init__.py b/src/io4dolfinx/__init__.py index 79fb409..1d89ff1 100644 --- a/src/io4dolfinx/__init__.py +++ b/src/io4dolfinx/__init__.py @@ -31,6 +31,7 @@ read_point_data, ) from .snapshot import snapshot_checkpoint +from .utils import reconstruct_mesh meta = metadata("io4dolfinx") __version__ = meta["Version"] @@ -62,4 +63,5 @@ "get_backend", "write_cell_data", "write_point_data", + "reconstruct_mesh", ] diff --git a/src/io4dolfinx/utils.py b/src/io4dolfinx/utils.py index 1144389..4efb763 100644 --- a/src/io4dolfinx/utils.py +++ b/src/io4dolfinx/utils.py @@ -14,9 +14,11 @@ from mpi4py import MPI +import basix.ufl import dolfinx import numpy as np import numpy.typing as npt +import ufl from packaging.version import Version __all__ = [ @@ -27,6 +29,7 @@ "unroll_dofmap", "compute_insert_position", "unroll_insert_position", + "reconstruct_mesh", ] @@ -196,3 +199,62 @@ def compute_dofmap_pos( local_cell[indicator] = cell_indicator[markers].reshape(-1) dof_pos[indicator] = local_indices[markers].reshape(-1) return local_cell, dof_pos + + +def reconstruct_mesh(mesh: dolfinx.mesh.Mesh, coordinate_element_degree: int) -> dolfinx.mesh.Mesh: + """ + Make a copy of a mesh and potentially change the element of the coordinate element. + + Note: + The topology is shared with the original mesh but the geometry is reconstructed. + + Args: + mesh: Mesh to reconstruct + coordinate_element_degree: Degree to use for coordinate element + + Returns: + The new mesh + + """ + # Extract cell properties + ud = mesh.ufl_domain() + assert ud is not None + c_el = ud.ufl_coordinate_element() + family = c_el.family_name + lvar = c_el.lagrange_variant + ct = c_el.cell_type + + # Create new UFL element + new_c_el = basix.ufl.element( + family, + ct, + coordinate_element_degree, + shape=(mesh.geometry.dim,), + lagrange_variant=lvar, + dtype=mesh.geometry.x.dtype, + ) + + # Extract new node coordinates + V_tmp = dolfinx.fem.functionspace(mesh, new_c_el) + gdim = mesh.geometry.dim + x = V_tmp.tabulate_dof_coordinates()[:, :gdim] + + # Create new geoemtry + geom_imap = V_tmp.dofmap.index_map + geom_dofmap = V_tmp.dofmap.list + num_nodes_local = geom_imap.size_local + geom_imap.num_ghosts + original_input_indices = geom_imap.local_to_global(np.arange(num_nodes_local, dtype=np.int32)) + coordinate_element = dolfinx.fem.coordinate_element( + mesh.topology.cell_type, coordinate_element_degree, lvar, dtype=mesh.geometry.x.dtype + ) + # Could use create_geometry here when things are fixed + geom = dolfinx.mesh.Geometry( + mesh.geometry._cpp_object.__class__( + geom_imap, geom_dofmap, coordinate_element._cpp_object, x, original_input_indices + ) + ) + + # Create new mesh + new_top = mesh.topology + cpp_mesh = mesh._cpp_object.__class__(mesh.comm, new_top._cpp_object, geom._cpp_object) + return dolfinx.mesh.Mesh(cpp_mesh, ufl.Mesh(new_c_el)) From a6a56b131d298d7379a977ed2dba12e993d8da14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Mon, 23 Feb 2026 11:41:27 +0100 Subject: [PATCH 2/3] Apply suggestions from code review Co-authored-by: Henrik Finsberg --- src/io4dolfinx/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io4dolfinx/utils.py b/src/io4dolfinx/utils.py index 4efb763..a6deae6 100644 --- a/src/io4dolfinx/utils.py +++ b/src/io4dolfinx/utils.py @@ -249,12 +249,12 @@ def reconstruct_mesh(mesh: dolfinx.mesh.Mesh, coordinate_element_degree: int) -> ) # Could use create_geometry here when things are fixed geom = dolfinx.mesh.Geometry( - mesh.geometry._cpp_object.__class__( + type(mesh.geometry._cpp_object)( geom_imap, geom_dofmap, coordinate_element._cpp_object, x, original_input_indices ) ) # Create new mesh new_top = mesh.topology - cpp_mesh = mesh._cpp_object.__class__(mesh.comm, new_top._cpp_object, geom._cpp_object) + cpp_mesh = type(mesh._cpp_object)(mesh.comm, new_top._cpp_object, geom._cpp_object) return dolfinx.mesh.Mesh(cpp_mesh, ufl.Mesh(new_c_el)) From 7d4784541913aee3a4e2a8c21373ecfd81693829 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Mon, 23 Feb 2026 21:30:45 +0000 Subject: [PATCH 3/3] Add test --- tests/test_mesh.py | 56 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 tests/test_mesh.py diff --git a/tests/test_mesh.py b/tests/test_mesh.py new file mode 100644 index 0000000..c09d38c --- /dev/null +++ b/tests/test_mesh.py @@ -0,0 +1,56 @@ +from mpi4py import MPI + +import dolfinx +import numpy as np +import pytest +import ufl + +from io4dolfinx import reconstruct_mesh + + +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("degree", [2, 3]) +@pytest.mark.parametrize("R", [0.1, 1, 10]) +def test_curve_mesh(degree, dtype, R): + N = 8 + mesh = dolfinx.mesh.create_rectangle( + MPI.COMM_WORLD, + [[-1, -1], [1, 1]], + [N, N], + diagonal=dolfinx.mesh.DiagonalType.crossed, + dtype=dtype, + ) + org_area = dolfinx.fem.form(1 * ufl.dx(domain=mesh), dtype=dtype) + + curved_mesh = reconstruct_mesh(mesh, degree) + + def transform(x): + u = R * x[0] * np.sqrt(1 - (x[1] ** 2 / (2))) + v = R * x[1] * np.sqrt(1 - (x[0] ** 2 / (2))) + return np.asarray([u, v]) + + curved_mesh.geometry.x[:, : curved_mesh.geometry.dim] = transform(curved_mesh.geometry.x.T).T + + area = dolfinx.fem.form(1 * ufl.dx(domain=curved_mesh), dtype=dtype) + circumference = dolfinx.fem.form(1 * ufl.ds(domain=curved_mesh), dtype=dtype) + + computed_area = curved_mesh.comm.allreduce(dolfinx.fem.assemble_scalar(area), op=MPI.SUM) + computed_circumference = curved_mesh.comm.allreduce( + dolfinx.fem.assemble_scalar(circumference), op=MPI.SUM + ) + + tol = 10 * np.finfo(dtype).eps + assert np.isclose(computed_area, np.pi * R**2, atol=tol) + assert np.isclose(computed_circumference, 2 * np.pi * R, atol=tol) + + linear_mesh = reconstruct_mesh(curved_mesh, 1) + linear_area = dolfinx.fem.form(1 * ufl.dx(domain=linear_mesh), dtype=dtype) + + recovered_area = linear_mesh.comm.allreduce( + dolfinx.fem.assemble_scalar(linear_area), op=MPI.SUM + ) + + # Curve original mesh + mesh.geometry.x[:, : mesh.geometry.dim] = transform(mesh.geometry.x.T).T + ref_area = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(org_area), op=MPI.SUM) + assert np.isclose(recovered_area, ref_area, atol=tol)