From a7222a1b3445544e34ba898840777a31ddb9be77 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 21 Oct 2026 10:44:08 +0000 Subject: [PATCH 1/2] Support checkpointing of RT/Nedelec on submeshes --- src/io4dolfinx/submesh.py | 14 +++++-- tests/test_submesh.py | 83 +++++++++++++++++++++++---------------- 2 files changed, 59 insertions(+), 38 deletions(-) diff --git a/src/io4dolfinx/submesh.py b/src/io4dolfinx/submesh.py index e17117e..43c76ca 100644 --- a/src/io4dolfinx/submesh.py +++ b/src/io4dolfinx/submesh.py @@ -354,9 +354,10 @@ def transfer_submesh_function( routed by. From :attr:`SubmeshCheckpoint.post_code`. Raises: - NotImplementedError: For H(div)/H(curl) spaces on a manifold submesh, - where DOLFINx's interpolation would return wrong values without - reporting an error. + NotImplementedError: For H(div)/H(curl) spaces on a manifold submesh + with a DOLFINx lacking ``FiniteElement.physical_base_value_size``, + whose interpolation would return wrong values without reporting an + error. """ V_dest = u_dest.function_space dest_mesh = V_dest.mesh @@ -365,7 +366,12 @@ def transfer_submesh_function( assert isinstance(comm, MPI.Intracomm) tdim = dest_mesh.topology.dim - if tdim < dest_mesh.geometry.dim and V_dest.element.needs_dof_transformations: + supports_manifold_transform = hasattr(dolfinx.fem.FiniteElement, "physical_base_value_size") + manifold_transformations = ( + tdim < dest_mesh.geometry.dim and V_dest.element.needs_dof_transformations + ) + + if manifold_transformations and not supports_manifold_transform: raise NotImplementedError( "Transferring an H(div)/H(curl) function onto a submesh with" f" tdim ({tdim}) < gdim ({dest_mesh.geometry.dim}) is not supported." diff --git a/tests/test_submesh.py b/tests/test_submesh.py index c2add2c..f978df7 100644 --- a/tests/test_submesh.py +++ b/tests/test_submesh.py @@ -178,19 +178,34 @@ def f(x): assert _max_error(u_sub, reference_sub) < 1e-13 -@pytest.mark.parametrize("family", ["N1curl", "RT"]) -def test_manifold_hcurl_roundtrips_but_transfer_raises(tmp_path, backend, family): - """H(div)/H(curl) on a manifold submesh: storable, but not reattachable. - - Writing and reading the standalone submesh is exact for these spaces, and - that is checked here so the capability is not lost by accident. +def _tangential_reference(x, family: str): + """A field tangent to the facet submesh ``x = 0`` that lies in ``(family, 1)``. - The transfer, though, is refused. DOLFINx cannot reconcile the reference and - physical value sizes of these families when ``tdim < gdim``, and its - interpolation does not report the problem -- it returns values that are - simply wrong (the L2 norm of an N1curl field measured 9.096 before and 4.904 - after). A loud refusal is the only safe behaviour. + On that plane RT1 contains ``a + c (y, z)`` and N1curl1 ``a + c (-z, y)``. + ``x`` is either the point array or ``ufl.SpatialCoordinate``. """ + zero = 0.0 * x[0] + if family == "RT": + return zero, 0.3 + 0.5 * x[1], -0.7 + 0.5 * x[2] + return zero, 0.3 - 0.5 * x[2], -0.7 + 0.5 * x[1] + + +def _l2_error(u, family: str) -> float: + """L2 distance between ``u`` and :func:`_tangential_reference`.""" + mesh = u.function_space.mesh + e = u - ufl.as_vector(_tangential_reference(ufl.SpatialCoordinate(mesh), family)) + local = dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(e, e) * ufl.dx)) + return np.sqrt(mesh.comm.allreduce(local, MPI.SUM)) + + +@pytest.mark.skipif( + not hasattr(dolfinx.fem.FiniteElement, "physical_base_value_size"), + reason="DOLFINx cannot interpolate into H(div)/H(curl) on a manifold", +) +@pytest.mark.parametrize("family", ["N1curl", "RT"]) +@pytest.mark.parametrize("degree", [1, 2]) +def test_manifold_hdiv_hcurl(tmp_path, backend, family, degree): + """H(div)/H(curl) on a manifold submesh: stored, read and reattached exactly.""" comm = MPI.COMM_WORLD folder = comm.bcast(tmp_path, root=0) parent_file = folder / f"parent{SUFFIX[backend]}" @@ -198,18 +213,9 @@ def test_manifold_hcurl_roundtrips_but_transfer_raises(tmp_path, backend, family mesh, submesh, cell_map = _make_submesh(comm, codim=1) # Spaces mapped by the contravariant Piola map need these on a manifold. fenicsx_compat.create_facet_orientations(submesh.topology) - V = dolfinx.fem.functionspace(submesh, (family, 1)) - u = dolfinx.fem.Function(V, name="u") - # These spaces cannot be interpolated into on a manifold, so fill by global - # dof index instead: reproducible no matter how the mesh is partitioned. - imap = V.dofmap.index_map - num_owned = imap.size_local * V.dofmap.index_map_bs - global_index = imap.local_to_global(np.arange(imap.size_local, dtype=np.int32)) - u.x.array[:num_owned] = np.sin(0.37 * global_index.astype(np.float64) + 0.11) - u.x.scatter_forward() - norm_written = comm.allreduce( - dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u, u) * ufl.dx)), MPI.SUM - ) + u = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, (family, degree)), name="u") + u.interpolate(lambda x: np.vstack(_tangential_reference(x, family))) + assert _l2_error(u, family) < 1e-13 io4dolfinx.write_mesh(parent_file, mesh, backend=backend) io4dolfinx.write_submesh( @@ -218,28 +224,37 @@ def test_manifold_hcurl_roundtrips_but_transfer_raises(tmp_path, backend, family io4dolfinx.write_function( parent_file, u, time=0.0, mode=io4dolfinx.FileMode.append, mesh_name="wall", backend=backend ) - del mesh, submesh, cell_map + del mesh, submesh, cell_map, u parent = io4dolfinx.read_mesh(parent_file, comm, backend=backend) stored = io4dolfinx.read_mesh(parent_file, comm, mesh_name="wall", backend=backend) fenicsx_compat.create_facet_orientations(stored.topology) - u_stored = dolfinx.fem.Function(dolfinx.fem.functionspace(stored, (family, 1)), name="u") + u_stored = dolfinx.fem.Function(dolfinx.fem.functionspace(stored, (family, degree)), name="u") io4dolfinx.read_function( parent_file, u_stored, time=0.0, name="u", mesh_name="wall", backend=backend ) - - # Steps 1-3 are exact even here. - norm_read = comm.allreduce( - dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.inner(u_stored, u_stored) * ufl.dx)), - MPI.SUM, - ) - assert np.isclose(norm_read, norm_written, rtol=0, atol=1e-12) + assert _l2_error(u_stored, family) < 1e-13 checkpoint = io4dolfinx.read_submesh(parent_file, parent, mesh_name="wall", backend=backend) fenicsx_compat.create_facet_orientations(checkpoint.submesh.topology) - u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(checkpoint.submesh, (family, 1))) + u_sub = dolfinx.fem.Function(dolfinx.fem.functionspace(checkpoint.submesh, (family, degree))) + io4dolfinx.transfer_submesh_function(u_stored, u_sub, checkpoint.post_code) + assert _l2_error(u_sub, family) < 1e-13 + + +@pytest.mark.skipif( + hasattr(dolfinx.fem.FiniteElement, "physical_base_value_size"), + reason="DOLFINx supports H(div)/H(curl) on a manifold", +) +@pytest.mark.parametrize("family", ["N1curl", "RT"]) +def test_manifold_hdiv_hcurl_transfer_raises(family): + """Older DOLFINx interpolates these spaces wrongly on a manifold, so refuse.""" + _, submesh, _ = _make_submesh(MPI.COMM_WORLD, codim=1) + fenicsx_compat.create_facet_orientations(submesh.topology) + u = dolfinx.fem.Function(dolfinx.fem.functionspace(submesh, (family, 1))) + post_code = np.zeros(0, dtype=np.int64) with pytest.raises(NotImplementedError, match="not supported"): - io4dolfinx.transfer_submesh_function(u_stored, u_sub, checkpoint.post_code) + io4dolfinx.transfer_submesh_function(u, u, post_code) @pytest.mark.parametrize("codim", [0, 1]) From 18897b57da6a3ea0d0661fd7dbe77411c2cc0284 Mon Sep 17 00:00:00 2001 From: Joergen Schartum Dokken Date: Wed, 21 Oct 2026 10:47:44 +0000 Subject: [PATCH 2/2] Update note about the limitations of checkpointing on submeshes, as it is all resolved on main --- docs/submesh.py | 39 +++++++++++++++++++++++---------------- src/io4dolfinx/submesh.py | 7 ++++--- 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/docs/submesh.py b/docs/submesh.py index 2a2ac47..48fac76 100644 --- a/docs/submesh.py +++ b/docs/submesh.py @@ -153,24 +153,31 @@ print(f"Area of the submesh assembled over the parent: {area:.3f}") # - -# ## What is not supported +# ## H(div) and H(curl) on a manifold # -# H(div) and H(curl) on a submesh of co-dimension greater than zero can be -# **written and read** exactly -- steps 1-3 above -- but cannot be transferred to -# a re-derived submesh. `transfer_submesh_function` raises `NotImplementedError` -# for them. +# On a submesh of co-dimension greater than zero, H(div) and H(curl) spaces +# depend on two DOLFINx changes that are on the main branch but in no release +# yet (0.11 and earlier lack them): # -# The reason is worth knowing. A finite element space on a manifold does not -# carry its value shape in *physical* space, so DOLFINx cannot reconcile the -# reference and physical value sizes of these families when `tdim < gdim` -# ([FEniCS/dolfinx#3619](https://github.com/FEniCS/dolfinx/issues/3619), open at -# the time of writing). Elsewhere that surfaces as an outright error: -# *"Interpolation: elements have different value dimensions"*, but in this -# particular path it does not: `interpolate_nonmatching` returns successfully -# with values that are wrong. The guard exists so that this -# shows up as an error rather than as a plausible-looking result, and it can be -# lifted once that issue is fixed. +# - [FEniCS/dolfinx#4511](https://github.com/FEniCS/dolfinx/pull/4511) gives a +# finite element its value size in *physical* space, which differs from the +# reference one when `tdim < gdim` +# ([FEniCS/dolfinx#3619](https://github.com/FEniCS/dolfinx/issues/3619)). +# - [FEniCS/dolfinx#4625](https://github.com/FEniCS/dolfinx/pull/4625) signs the +# facets on which neighbouring cells disagree on the normal, without which a +# contravariant Piola space such as Raviart-Thomas is not conforming. Call +# `create_facet_orientations` on the topology of every manifold mesh before +# creating such a space on it -- the original submesh, the stored one and the +# re-derived one -- or use `fenicsx_compat.create_facet_orientations`, which +# does nothing on DOLFINx versions without it. # -# There is no sound workaround. The stored submesh is an independent mesh, +# With both, writing, reading and transferring are exact for these spaces too. +# +# On older DOLFINx they can still be **written and read** exactly, as under +# *Writing* and *Reading* above, but not transferred to a re-derived submesh: +# `transfer_submesh_function` raises `NotImplementedError`. There, +# `interpolate_nonmatching` returns successfully with values that are wrong, +# and the guard makes that an error rather than a plausible-looking result. +# There is no sound workaround: the stored submesh is an independent mesh, # derived from nothing, so a projection -- or any other form assembled between it # and the re-derived submesh -- would have no basis. diff --git a/src/io4dolfinx/submesh.py b/src/io4dolfinx/submesh.py index 43c76ca..572c8cf 100644 --- a/src/io4dolfinx/submesh.py +++ b/src/io4dolfinx/submesh.py @@ -375,9 +375,10 @@ def transfer_submesh_function( raise NotImplementedError( "Transferring an H(div)/H(curl) function onto a submesh with" f" tdim ({tdim}) < gdim ({dest_mesh.geometry.dim}) is not supported." - " DOLFINx cannot reconcile the reference and physical value sizes of" - " these families on a manifold (FEniCS/dolfinx#3619," - " https://github.com/FEniCS/dolfinx/issues/3619), and the interpolation" + " This DOLFINx cannot reconcile the reference and physical value sizes" + " of these families on a manifold (FEniCS/dolfinx#3619," + " https://github.com/FEniCS/dolfinx/issues/3619, fixed on the main" + " branch by FEniCS/dolfinx#4511), and the interpolation" " this transfer relies on does not report that: it returns silently" " with values that are wrong (measured: the L2 norm of an N1curl field" " on a facet submesh fell from 9.096 to 4.904). Writing and reading the"