Replace scifem.compat and scifem.ufl_compat with fenicsx-compat

README.md

@@ -63,6 +63,8 @@ We aim to support the last two stable releases of DOLFINx. Sometimes features ar
 
 Some features might depend on functionality added to DOLFINx, which means that they won't be backwards compatible.
 
+Workarounds for API differences between DOLFINx versions live in [`fenicsx-compat`](https://github.com/scientificcomputing/fenicsx-compat), which `scifem` depends on. It requires DOLFINx 0.10 or newer and Python 3.12 or newer.
+
 # Installation
 
 The package is partly written in C++ and relies on `dolfinx`. Users are encouraged to install `scifem` with `pip` in an environment where `dolfinx` is already installed or with `conda`.

docs/api.rst

@@ -134,13 +134,3 @@ XDMF output
 
 .. automodule:: scifem.xdmf
     :members:
-
-
-Compat functions
-################
-
-.. automodule:: scifem.compat
-    :members:
-
-.. automodule:: scifem.ufl_compat
-    :members:

examples/periodic_poisson.py

@@ -48,8 +48,8 @@
 # :class: tip dropdown
 #
 # On the `main` branch of DOLFINx, ghost mode is supplied directly to {py:func}`dolfinx.mesh.create_mesh`,
-# rather than through the partitioner. Use {py:func}`scifem.compat.create_partitioner` to get
-# a partitioner that works on all versions.
+# rather than through the partitioner. Use {py:func}`fenicsx_compat.create_mesh` to supply a
+# ghost mode in a way that works on all versions.
 # ```
 
 N = 25

examples/transfer_tags_to_submesh.py

@@ -12,7 +12,8 @@
 import os
 import sys
 import pyvista
-from scifem import transfer_meshtags_to_submesh, compat
+import fenicsx_compat
+from scifem import transfer_meshtags_to_submesh
 from mpi4py import MPI
 import gmsh
 import dolfinx
@@ -82,7 +83,7 @@ def plot_mesh(mesh: dolfinx.mesh.Mesh, values=None, figure=None):
     plotter = pyvista.Plotter()
     V_linear = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
     linear_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(V_linear))
-    if compat.cmap(mesh).degree > 1:
+    if fenicsx_compat.cmap(mesh).degree > 1:
         ugrid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
         if values is not None:
             ugrid.cell_data["Marker"] = values

pyproject.toml

@@ -13,13 +13,19 @@ name = "scifem"
 version = "0.27.0.dev0"
 description = "Scientific tools for finite element methods"
 readme = "README.md"
-requires-python = ">=3.10"
+requires-python = ">=3.12"
 license = "MIT"
 authors = [
     { email = "dokken@simula.no", name = "Jørgen S. Dokken" },
     { email = "henriknf@simula.no", name = "Henrik N.T. Finsberg" },
 ]
-dependencies = ["fenics-dolfinx", "numpy", "packaging", "scipy", "typing-extensions"]
+dependencies = [
+    "fenics-dolfinx>=0.10.0",
+    "fenicsx-compat>=0.2.0",
+    "numpy",
+    "scipy",
+    "typing-extensions",
+]
 
 [project.urls]
 repository = "https://github.com/scientificcomputing/scifem.git"
@@ -38,7 +44,7 @@ docs = [
     "scifem[gmsh]",
 ]
 dev = ["ruff", "mypy", "bump-my-version", "pre-commit"]
-test = ["pytest", "petsc4py", "h5py", "scifem[biomed]"]
+test = ["pytest", "petsc4py", "h5py", "packaging", "scifem[biomed]"]
 biomed = ["nibabel"]
 gmsh = ["gmsh"]
 all = ["scifem[docs,dev,test,adios2,h5py,biomed,gmsh]"]
@@ -98,8 +104,8 @@ exclude = [
 # Same as Black.
 line-length = 100
 
-# Assume Python 3.10.
-target-version = "py310"
+# Assume Python 3.12.
+target-version = "py312"
 
 [tool.ruff.lint]
 # Enable pycodestyle (`E`) and Pyflakes (`F`) codes by default.

src/scifem/__init__.py

@@ -24,7 +24,6 @@
 from .interpolation import interpolation_matrix, prepare_interpolation_data
 from .geometry import closest_point_projection
 from scifem import periodic
-from scifem import compat
 
 meta = metadata("scifem")
 __version__ = meta["Version"]
@@ -44,7 +43,6 @@
     "closest_point_projection",
     "assemble_scalar",
     "create_space_of_simple_functions",
-    "compat",
     "compute_interface_data",
     "compute_subdomain_exterior_facets",
     "create_real_functionspace",

src/scifem/assembly.py

@@ -18,15 +18,14 @@ def _extract_dtype(expr: ufl.form.Form | ufl.core.expr.Expr | dolfinx.fem.Form)
 
     """
     try:
-        return expr.dtype
+        # Compiled forms, as well as Functions and Constants, carry a dtype
+        return expr.dtype  # type: ignore[union-attr]
     except AttributeError:
-        try:
-            geom_type = ufl.domain.extract_domains(expr)[0].ufl_cargo().geometry.x.dtype
-        except AttributeError:
-            # Legacy support
-            geom_type = expr.ufl_domain().ufl_cargo().geometry.x.dtype
-        coefficients = expr.coefficients()
-        constants = expr.constants()
+        form = typing.cast(ufl.form.Form, expr)
+        domain = typing.cast(ufl.Mesh, ufl.domain.extract_domains(form)[0])
+        geom_type = domain.ufl_cargo().geometry.x.dtype
+        coefficients = form.coefficients()
+        constants = form.constants()
         if (not coefficients) and (not constants):
             return geom_type
         else:
@@ -44,7 +43,7 @@ def _extract_dtype(expr: ufl.form.Form | ufl.core.expr.Expr | dolfinx.fem.Form)
 
 def assemble_scalar(
     J: ufl.form.Form | dolfinx.fem.Form,
-    entity_maps: typing.Optional[dict[dolfinx.mesh.Mesh, npt.NDArray[np.int32]]] = None,
+    entity_maps: typing.Sequence[dolfinx.mesh.EntityMap] | None = None,
 ) -> np.floating | np.complexfloating:
     """Assemble a scalar form and gather result across processes
 
@@ -56,7 +55,8 @@ def assemble_scalar(
         The accumulated value of the assembled form.
     """
     dtype = _extract_dtype(J)
-    compiled_form = dolfinx.fem.form(J, entity_maps=entity_maps, dtype=dtype)
+    # dolfinx.fem.form passes an already compiled form through unchanged
+    compiled_form = dolfinx.fem.form(J, entity_maps=entity_maps, dtype=dtype)  # type: ignore[arg-type]
 
     if (rank := compiled_form.rank) != 0:
         raise ValueError(f"Form must be a scalar form, got for of arity {rank}")
@@ -67,7 +67,7 @@ def assemble_scalar(
 def norm(
     expr: ufl.core.expr.Expr,
     norm_type: Literal["L2", "H1", "H10"],
-    entity_maps: typing.Optional[dict[dolfinx.mesh.Mesh, npt.NDArray[np.int32]]] = None,
+    entity_maps: typing.Sequence[dolfinx.mesh.EntityMap] | None = None,
 ) -> dolfinx.fem.Form:
     """
     Compile the norm of an UFL expression into a DOLFINx form.

src/scifem/bcs.py

@@ -5,59 +5,13 @@
 import ufl
 import numpy.typing as npt
 import numpy as np
-from packaging.version import Version
+import fenicsx_compat
 
-from .compat import get_facet_permutations
 from .utils import group_by_key
-from ffcx.ir.elementtables import (
-    permute_quadrature_interval,
-    permute_quadrature_triangle,
-    permute_quadrature_quadrilateral,
-)
 
 __all__ = ["interpolate_function_onto_facet_dofs"]
 
 
-def build_quadrature_permutations(facet_type, points):
-    if Version(dolfinx.__version__) < Version("0.11.0.dev0"):
-        # In older versions of dolfinx, the permutation is handled internally
-        # in the C++ code, so we can just return the original points
-        if facet_type == basix.CellType.interval:
-            num_permutations = 2
-        elif facet_type == basix.CellType.triangle:
-            num_permutations = 6
-        elif facet_type == basix.CellType.quadrilateral:
-            num_permutations = 8
-        else:
-            raise ValueError(f"Unsupported {facet_type=}")
-        return [points for _ in range(num_permutations)]
-    else:
-        if facet_type == basix.CellType.interval:
-            return [permute_quadrature_interval(points, ref) for ref in range(2)]
-        elif facet_type == basix.CellType.triangle:
-            perms = []
-            # FFCx order: rot is outer loop, ref is inner loop
-            for rot in range(3):
-                for ref in range(2):
-                    # Counteract the mapping with the inverse permutation
-                    rot_inv = (3 - rot) % 3 if ref == 0 else rot
-                    perms.append(permute_quadrature_triangle(points, ref, rot_inv))
-            return perms
-
-        elif facet_type == basix.CellType.quadrilateral:
-            perms = []
-            # FFCx order: rot is outer loop, ref is inner loop
-            for rot in range(4):
-                for ref in range(2):
-                    # Counteract the mapping with the inverse permutation
-                    rot_inv = (4 - rot) % 4 if ref == 0 else rot
-                    perms.append(permute_quadrature_quadrilateral(points, ref, rot_inv))
-            return perms
-
-        else:
-            raise ValueError(f"Unsupported {facet_type=}")
-
-
 def pull_back_to_reference_facet(
     cell_type: basix.CellType, points_per_facet: Sequence[npt.NDArray[np.floating]]
 ) -> npt.NDArray[np.floating]:
@@ -125,7 +79,7 @@ def interpolate_function_onto_facet_dofs(
     reference_facet_points = pull_back_to_reference_facet(
         domain.basix_cell(), interpolation_points[fdim]
     )
-    facet_points = build_quadrature_permutations(facet_type, reference_facet_points)
+    facet_points = fenicsx_compat.permute_facet_quadrature(facet_type, reference_facet_points)
     expressions = [dolfinx.fem.Expression(expr, points) for points in facet_points]
     points_per_entity = [sum(ip.shape[0] for ip in ips) for ips in interpolation_points]
     offsets = np.zeros(domain.topology.dim + 2, dtype=np.int32)
@@ -147,7 +101,7 @@ def interpolate_function_onto_facet_dofs(
     is_marked = np.zeros(num_facets_on_process, dtype=np.int8)
     is_marked[facets] = 1
     num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(domain.topology.cell_type, fdim)
-    facet_permutations = get_facet_permutations(domain.topology)
+    facet_permutations = fenicsx_compat.facet_permutations(domain.topology)
     # One evaluation per facet permutation, over every marked facet with that permutation
     local_facets = np.asarray(c_to_f.array, dtype=np.int32).reshape(-1, num_facets_per_cell)
     rows, slots = np.nonzero(is_marked[local_facets[all_connected_cells]])
@@ -157,21 +111,15 @@ def interpolate_function_onto_facet_dofs(
     values = values.reshape(len(all_connected_cells), offsets[-1], expr_value_size)
     for perm, selected in group_by_key(permutations):
         entities = np.column_stack((cells[selected], slots[selected])).astype(np.int32)
-        # Backwards compatibility
-        try:
-            evaluated = expressions[perm].eval(domain, entities)
-        except (AttributeError, AssertionError):
-            evaluated = expressions[perm].eval(domain, entities.flatten())
+        evaluated = fenicsx_compat.expression_eval(expressions[perm], domain, entities)
         positions = offsets[fdim] + num_points * slots[selected, None] + np.arange(num_points)
         values[rows[selected, None], positions] = evaluated.reshape(
             len(selected), num_points, expr_value_size
         )
 
     qh = dolfinx.fem.Function(Q)
-    if hasattr(qh._cpp_object, "interpolate_f"):
-        interpolate_func = qh._cpp_object.interpolate_f
-    else:
-        interpolate_func = qh._cpp_object.interpolate
-    interpolate_func(values.reshape(-1, expr_value_size).T.copy(), all_connected_cells)
+    fenicsx_compat.interpolate(
+        qh._cpp_object, values.reshape(-1, expr_value_size).T.copy(), all_connected_cells
+    )
     qh.x.scatter_forward()
     return qh

src/scifem/biomedical.py

@@ -112,7 +112,7 @@ def read_mri_data_to_function(
     cells: npt.NDArray[np.int32] | None = None,
     vox2ras_transform: npt.NDArray[np.floating] | None = None,
     degree: int = 0,
-    dtype: np.inexact = dolfinx.default_scalar_type,
+    dtype: type[np.inexact] = dolfinx.default_scalar_type,
     use_tkr: bool = False,
 ) -> dolfinx.fem.Function:
     """Read in MRI data over a set of cells in the mesh and attach it to an appropriate function.
@@ -143,7 +143,9 @@ def read_mri_data_to_function(
         raise ValueError("Degree must be a non-negative integer")
 
     if degree == 0:
-        el = basix.ufl.element("Lagrange", mesh.basix_cell(), 0, discontinuous=True)
+        el: basix.ufl._ElementBase = basix.ufl.element(
+            "Lagrange", mesh.basix_cell(), 0, discontinuous=True
+        )
     else:
         el = basix.ufl.quadrature_element(scheme="default", degree=degree, cell=mesh.basix_cell())