Add support for non-cartesian manifolds in codim problems by ee-nn · Pull Request #1264 · festim-dev/FESTIM

ee-nn · 2026-09-12T05:33:12Z

Description

v2.2-rc.2 released co-dimensional subdomains but only for cartesian coordinates. This adds functionality for cylindrical and cartesian coordinates.

Related Issues

Fixes #1263

Type of Change

🐛 Bug fix (non-breaking change which fixes an issue)
✨ New feature (non-breaking change which adds functionality)
💥 Breaking change (fix or feature that would cause existing functionality to not work as expected)
🔨 Code refactoring (no functional changes, no API changes)
📝 Documentation update
✅ Test update (adding missing tests or correcting existing tests)
🔧 Build/CI configuration change

Testing

All existing tests pass locally (pytest)
I have added new tests that prove my fix is effective or that my feature works

Code Quality Checklist

My code follows the code style of this project (Ruff formatted: ruff format .)
My code passes linting checks (ruff check .)
I have performed a self-review of my own code
I have commented my code, particularly in hard-to-understand areas

Documentation

I have updated the documentation accordingly (if applicable)
I have added docstrings to new functions/classes following the project conventions

Screenshots/Examples

Here is a small mwe:

from mpi4py import MPI

import dolfinx
import numpy as np

import festim as F


def run(name, mesh):
    bulk = F.VolumeSubdomain(id=1, material=F.Material(D_0=1, E_D=0))
    surface = F.VolumeSubdomain(  # the co-dim part 
        id=2,
        dim=mesh.vdim - 1,
        locator=lambda x: np.isclose(x[0], 1.0),
        material=F.Material(D_0=0.1, E_D=0),
    )
    c_bulk = F.Species("bulk", subdomains=[bulk])
    c_surface = F.Species("surface", subdomains=[surface])
    dependencies = {"cb": c_bulk, "cs": c_surface}

    # Hydrogen initially in the bulk transfers to the initially empty surface.
    # The pipe starts with an axial gradient so its surface profile is nonuniform.
    model = F.HydrogenTransportProblemDiscontinuous(
        mesh=mesh,
        subdomains=[bulk, surface],
        species=[c_bulk, c_surface],
        temperature=500,
        initial_conditions=[
            F.InitialConcentration(
                value=(lambda x: 1 + x[1]) if name == "pipe" else 1.0,
                species=c_bulk,
                volume=bulk,
            ),
        ],
        boundary_conditions=[
            F.ParticleFluxBC(
                subdomain=surface,
                species=c_bulk,
                value=lambda cb, cs: cs - cb,  # assume k = 1
                species_dependent_value=dependencies,
            ),
        ],
        sources=[
            F.ParticleSource(
                volume=surface,
                species=c_surface,
                value=lambda cb, cs: cb - cs,
                species_dependent_value=dependencies,
            ),
        ],
        # A codim manifold carries its own equation, so use a volume quantity.
        exports=[
            F.AverageVolume(
                field=c_surface,
                volume=surface,
                filename=f"{name}_surface.csv",
            ),
        ],
        settings=F.Settings(
            transient=True,
            stepsize=0.01,
            final_time=0.2,
            atol=1e-10,
            rtol=1e-10,
        ),
    )
    model.initialise()
    model.run()


# Pipe: inner radius 0.5, outer radius 1, length 2; outer wall has dim=1.
pipe = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, [np.array([0.5, 0.0]), np.array([1.0, 2.0])], [8, 24]
)
run("pipe", F.Mesh(pipe, coordinate_system="cylindrical"))

# Ball: radius 1; outer spherical surface has dim=0 in the radial mesh.
ball = dolfinx.mesh.create_interval(MPI.COMM_WORLD, 32, [0.0, 1.0])
run("ball", F.Mesh(ball, coordinate_system="spherical"))

…omains

codecov · 2026-09-12T05:40:18Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.03%. Comparing base (3e37dab) to head (822f47f).

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1264      +/-   ##
==========================================
+ Coverage   95.86%   96.03%   +0.16%     
==========================================
  Files          57       57              
  Lines        5036     5065      +29     
==========================================
+ Hits         4828     4864      +36     
+ Misses        208      201       -7

Flag	Coverage Δ
dolfinx-v0.10.0	`93.08% <100.00%> (+0.17%)`	⬆️
dolfinx-v0.11.0	`95.89% <98.41%> (+0.14%)`	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

RemDelaporteMathurin · 2026-09-14T12:50:41Z

+    def integration_weight(self, mesh):
+        """Physical measure per unit mesh measure, assuming full angular coverage.
+
+        Use the integration domain's coordinates, including when it is a manifold
+        submesh. Cylindrical 1D quantities are per unit axial length.
+        """
+        if self == CoordinateSystem.CARTESIAN:
+            return 1
+        r = spatial_coordinate(mesh)[0]
+        if self == CoordinateSystem.CYLINDRICAL:
+            return 2 * math.pi * r
+        return 4 * math.pi * r**2


Good idea to store this in the CoordinateSystem class.

RemDelaporteMathurin · 2026-09-14T12:51:17Z

    def define_meshtags(self, surface_subdomains, volume_subdomains, interfaces=None):
        # check if all borders are defined
-        self.check_borders(volume_subdomains)
+        self.check_borders([v for v in volume_subdomains if v.codim(self.vdim) == 0])


Good catch. Is there a test catching this bug?

I don't believe so, I can add one

Added a small direct unit test for this just to be safe, although several of the system tests added do fail if this bug is present

RemDelaporteMathurin · 2026-09-14T12:52:54Z

            "Lagrange",
            function_space.mesh.topology.cell_name(),
-            1,
+            0 if function_space.mesh.topology.dim == 0 else 1,


advection in 0D?

I don't even think basix allows ("Lagrange", 0), is this tested?

Sorry I glossed over this initially! You're right I don't think this makes sense, perhaps we can replace this with a warning/error, similar to your comment further down.

I did double-check the syntax and basix.ufl.element("Lagrange", "point", 0) does work. "point" should already be what is set from dolfinx.mesh.create_submesh.

RemDelaporteMathurin · 2026-09-14T12:56:15Z

+def spatial_coordinate(mesh):
+    """Coordinates on an integration mesh, including a zero-dimensional submesh.
+
+    FFCx 0.11 cannot compile SpatialCoordinate on a vertex cell (its coordinate
+    table is identically one). A P0 coefficient stores the exact coordinates at
+    each point, also when a point subdomain contains several disconnected points.
+    """
+    domain = mesh if isinstance(mesh, ufl.Mesh) else mesh.ufl_domain()
+    # UFL exposed this as a method before the release bundled with DOLFINx 0.11.
+    tdim = domain.topological_dimension
+    if callable(tdim):
+        tdim = tdim()
+    if tdim > 0:
+        return ufl.SpatialCoordinate(domain)
+    if isinstance(mesh, ufl.Mesh):
+        mesh = dolfinx.mesh.Mesh(domain.ufl_cargo(), domain)
+    gdim = mesh.geometry.dim
+    V = fem.functionspace(mesh, ("P", 0, (gdim,)))
+    coordinate = fem.Function(V)
+    coordinate.interpolate(lambda x: x[:gdim])
+    coordinate.x.scatter_forward()
+    return coordinate


So this is only to catch the situation where you have a vertex cell (a point mesh)?

In that case, could we not have a try/except instead?

try: return ufl.SpatialCoordinate(domain) except [INSERT ERROR WITH 0D MESH]: tdim = .... V = .... coordinate = ...

Yes we could do something like that instead, however after a bit of testing I realize that we'd need to wrap ufl.SpatialCoordinate with dolfinx.fem.form() because it only fails at form compilation

I left it as-is for now since trying to return ufl.SpatialCoordinate(domain) won't directly trip a try/except, but let me know if you think it would still be good to have a try/except block here

RemDelaporteMathurin · 2026-09-14T12:58:02Z

+        # A radial surface is a point in a 1D mesh. Basix only supports P0 there.
+        if subdomain.submesh.topology.dim == 0:
+            element_degree = 0


Do we agree that this is not necessarily related to non-cartesian manifolds but 0D manifolds in general yes?

Yes that's my understanding as well

RemDelaporteMathurin · 2026-09-14T13:00:33Z

+            if subdomain.submesh.topology.dim == 0:
+                # A point has no tangential direction in which a species can drift.
+                continue


Should we raise a warning or maybe even an error here?

Yeah I agree, I think a warning might be good

…M into 1263-non-cartesian-manifolds

ee-nn · 2026-09-18T15:46:12Z

@RemDelaporteMathurin Let me know if this looks okay now, thanks for all the feedback so far! Once it's merged it should make #1271 easier

RemDelaporteMathurin · 2026-09-18T16:01:20Z

@ee-nn i'm traveling today so I won't have time to review it again. @hhy2022 Could you base this other PR off this branch instead?

hhy2022 · 2026-09-18T16:09:37Z

Sure! I will base my PR on this branch when I start working on #1271. @ee-nn @RemDelaporteMathurin

ee-nn added 2 commits September 11, 2026 23:04

add cylindrical and spherical coordinate functionality for codim subd…

3aa7b4d

…omains

update docs

9ffcfaf

ee-nn self-assigned this Sep 12, 2026

ee-nn requested review from RemDelaporteMathurin and jhdark September 13, 2026 15:16

RemDelaporteMathurin reviewed Sep 14, 2026

View reviewed changes

ee-nn and others added 5 commits September 15, 2026 21:33

Add unit test for mesh_1d border check

2ba11b0

add warnings for drift on 0D domain

8bbc471

Merge branch 'main' into 1263-non-cartesian-manifolds

0e782e0

add test for warnings

238389c

Merge branch '1263-non-cartesian-manifolds' of github.com:ee-nn/FESTI…

822f47f

…M into 1263-non-cartesian-manifolds

ee-nn requested a review from RemDelaporteMathurin September 16, 2026 17:11

This was referenced Sep 16, 2026

Support cylindrical meshes in the gas enclosure #1271

Open

Derived quantities for cylindrical and spherical coordinate systems #1023

Closed

Conversation

ee-nn commented Sep 12, 2026

Description

Related Issues

Type of Change

Testing

Code Quality Checklist

Documentation

Screenshots/Examples

Uh oh!

codecov Bot commented Sep 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ee-nn commented Sep 18, 2026

Uh oh!

RemDelaporteMathurin commented Sep 18, 2026

Uh oh!

hhy2022 commented Sep 18, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

codecov Bot commented Sep 12, 2026 •

edited

Loading