diff --git a/README.md b/README.md index 3a6b397..05ab1e4 100644 --- a/README.md +++ b/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`. diff --git a/docs/api.rst b/docs/api.rst index e9db16c..7db0ece 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -134,13 +134,3 @@ XDMF output .. automodule:: scifem.xdmf :members: - - -Compat functions -################ - -.. automodule:: scifem.compat - :members: - -.. automodule:: scifem.ufl_compat - :members: diff --git a/examples/periodic_poisson.py b/examples/periodic_poisson.py index f553f61..40cf907 100644 --- a/examples/periodic_poisson.py +++ b/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 diff --git a/examples/transfer_tags_to_submesh.py b/examples/transfer_tags_to_submesh.py index a081f1e..35f7808 100644 --- a/examples/transfer_tags_to_submesh.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index ca85f85..96d792b 100644 --- a/pyproject.toml +++ b/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. diff --git a/src/scifem/__init__.py b/src/scifem/__init__.py index 11df3d7..c591fe8 100644 --- a/src/scifem/__init__.py +++ b/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", diff --git a/src/scifem/assembly.py b/src/scifem/assembly.py index 9f4b747..af65ddc 100644 --- a/src/scifem/assembly.py +++ b/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. diff --git a/src/scifem/bcs.py b/src/scifem/bcs.py index 592708d..264d375 100644 --- a/src/scifem/bcs.py +++ b/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 diff --git a/src/scifem/biomedical.py b/src/scifem/biomedical.py index bade023..a7e6155 100644 --- a/src/scifem/biomedical.py +++ b/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()) diff --git a/src/scifem/compat.py b/src/scifem/compat.py deleted file mode 100644 index a2f4209..0000000 --- a/src/scifem/compat.py +++ /dev/null @@ -1,278 +0,0 @@ -"""Layer for small backward compatibility wrappers for DOLFINx""" - -import numpy.typing as npt -import numpy as np -import dolfinx -import inspect - -__all__ = [ - "create_partitioner", - "create_cell_permutations", - "compute_integration_domains", - "get_facet_permutations", - "create_cpp_finite_element", -] - - -def create_partitioner( - ghost_mode: dolfinx.mesh.GhostMode = dolfinx.mesh.GhostMode.shared_facet, - max_facet_to_cell_links: int = 2, -): - """Create a partitioner across the supported DOLFINx versions. - - The constructor has changed shape more than once and none of the forms is - introspectable, so they are told apart by the `TypeError` the call itself raises: - - 1. communicator and ghost mode only, everything else through setters; - 2. the same with `max_facet_to_cell_links` passed positionally; - 3. as (2) without the communicator. - - Args: - comm: The communicator of the new partitioner. - ghost_mode: The ghosting mode to use. - max_facet_to_cell_links: Maximum number of cells that can share a facet. Used by - forms (2) and (3) only, which take it directly; form (1) does not accept it. - - Returns: - The default partitioner - """ - if not hasattr(dolfinx.mesh, "create_cell_partitioner"): - partitioner = dolfinx.graph.partitioner() - else: - sig = inspect.signature(dolfinx.mesh.create_cell_partitioner) - part_kwargs = {} - if "max_facet_to_cell_links" in sig.parameters: - part_kwargs["max_facet_to_cell_links"] = max_facet_to_cell_links - - partitioner = dolfinx.mesh.create_cell_partitioner(ghost_mode, **part_kwargs) - return partitioner - - -def cmap(mesh: dolfinx.mesh.Mesh) -> dolfinx.fem.CoordinateElement: - """Get the basix Cmap for the mesh.""" - if hasattr(mesh.geometry, "cmaps"): - return mesh.geometry.cmaps[0] - if callable(mesh.geometry.cmap): - return mesh.geometry.cmap() - else: - return mesh.geometry.cmap - - -def dofmap(mesh: dolfinx.mesh.Mesh) -> npt.NDArray[np.int32]: - """Get the dofmap for the geometry.""" - if hasattr(mesh.geometry, "dofmaps"): - return mesh.geometry.dofmaps[0] - if callable(mesh.geometry.dofmap): - return mesh.geometry.dofmap() - else: - return mesh.geometry.dofmap - - -def form_map(form: dolfinx.fem.Form) -> tuple[dolfinx.common.IndexMap, int]: - try: - return ( - form.function_spaces[0].dofmaps(0).index_map, - form.function_spaces[0].dofmaps(0).index_map_bs, - ) - - except TypeError: - return ( - form.function_spaces[0].dofmaps[0].index_map, - form.function_spaces[0].dofmaps[0].index_map_bs, - ) - - -def index_map(comm, size_local, ghosts, owners, tag: int | None = None): - if dolfinx.common.IndexMap != dolfinx.cpp.common.IndexMap: - assert tag is not None, "Tag must be provided for dolfinx.common.index_map" - return dolfinx.common.index_map(comm, size_local, ghosts=(ghosts, owners), tag=tag) - else: - try: - return dolfinx.common.IndexMap(comm, size_local, ghosts, owners) - except TypeError: - assert tag is not None, "Tag must be provided for dolfinx.common.IndexMap" - return dolfinx.common.IndexMap(comm, size_local, ghosts, owners, tag=tag) - - -def extract_cpp_object(obj): - if hasattr(obj, "_cpp_object"): - return obj._cpp_object - else: - return obj - - -def topology(comm, cell_type, tdim, vertex_map, cell_map, c_to_v, v_to_v, original_cell_index): - """Construct a ``dolfinx.cpp.mesh.Topology`` across the supported DOLFINx versions. - - The constructor has changed shape more than once and none of the forms is - introspectable, so they are told apart by the `TypeError` the call itself raises: - - 1. communicator and cell type only, everything else through setters; - 2. the same with the maps, dofmap and original cell index passed positionally; - 3. as (2) without the communicator. - - Args: - comm: The communicator of the new topology. - cell_type: Its cell type. - tdim: Its topological dimension, for the index map and connectivity it is set at. - vertex_map: Index map for dimension 0. - cell_map: Index map for dimension `tdim`. - c_to_v: Cell-to-vertex connectivity, in the local numbering of `vertex_map`. - v_to_v: Vertex-to-vertex connectivity, i.e. the identity over `vertex_map`. Used - by form (1) only, which cannot derive it. - original_cell_index: Input global index of each cell. Used by forms (2) and (3), - which take it directly; form (1) does not accept it. - - Returns: - The topology, with both index maps and both connectivities set. - """ - try: - topology = dolfinx.cpp.mesh.Topology(comm, cell_type) - topology.set_index_map(0, vertex_map) - topology.set_index_map(tdim, cell_map) - topology.set_connectivity(v_to_v, 0, 0) - topology.set_connectivity(c_to_v, tdim, 0) - return topology - except TypeError: - pass - - args = ( - cell_type, - extract_cpp_object(vertex_map), - extract_cpp_object(cell_map), - extract_cpp_object(c_to_v), - original_cell_index, - ) - try: - return dolfinx.cpp.mesh.Topology(comm, *args) - except TypeError: - return dolfinx.cpp.mesh.Topology(*args) - - -def index_to_dest_ranks( - imap: dolfinx.common.IndexMap | dolfinx.cpp.common.IndexMap, tag: int -) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]: - """Get the destination ranks for each index in an index map, across DOLFINx versions. - - Args: - imap: The index map to query. - tag: MPI tag for the consensus exchange, the same on every process and not in use - by another exchange in flight. - - Returns: - (Sharing-rank data, offsets). Ranks sharing local - index `i` occupy `[offsets[i], offsets[i + 1])`. For an owned - index, these are the ranks that ghost it; for a ghost index, its - owner and the other ranks that ghost it. The caller is excluded. - """ - signature_inputs = inspect.signature(imap.index_to_dest_ranks).parameters - if len(signature_inputs) == 2: - return imap.index_to_dest_ranks(tag) # type: ignore[call-arg] - elif len(signature_inputs) == 0: - return imap.index_to_dest_ranks() - else: - raise RuntimeError( - f"Unexpected signature for index_to_dest_ranks: {signature_inputs}. " - "Please report this to the scifem developers." - ) - - -def ghosting_ranks(index_map, tag: int): - """The ranks that ghost each owned index of `index_map`, flat and with offsets. - - Compat: `IndexMap.index_to_dest_ranks` hands back the two arrays directly in DOLFINx - 0.12 and wraps them in an adjacency list in 0.11. - - Args: - index_map: The map to read the ghosting ranks of. - tag: MPI tag for the consensus exchange, the same on every process and not in use - by another exchange in flight. - - Returns: - ``(ranks, offsets)``: the ghosting ranks of owned index ``i`` are - ``ranks[offsets[i]:offsets[i + 1]]``. Collective. - """ - dest = index_to_dest_ranks(index_map, tag) - if isinstance(dest, tuple): - ranks, offsets = dest - else: - ranks, offsets = (dest.array, dest.offsets) # type: ignore - return np.asarray(ranks, dtype=np.int32), np.asarray(offsets, dtype=np.int64) - - -def create_cell_permutations(topology: dolfinx.mesh.Topology): - """Compute the packed per-cell permutation info, across DOLFINx versions. - - Args: - topology: The topology to compute the permutation info of. - """ - if hasattr(topology, "create_cell_permutations"): - topology.create_cell_permutations() - else: # DOLFINx < 0.10, where it takes no dimension - topology.create_entity_permutations() # type: ignore[call-arg] - - -def compute_integration_domains( - integral_type: dolfinx.fem.IntegralType, - topology: dolfinx.mesh.Topology, - entities: npt.NDArray[np.int32], -) -> npt.NDArray[np.int32]: - """:py:func:`dolfinx.fem.compute_integration_domains` across DOLFINx versions. - - Older versions also take the dimension of ``entities``, which follows from - ``integral_type``: the cells for a cell integral, the facets otherwise. - - Args: - integral_type: The type of integral the entities are for. - topology: The topology of the mesh the entities belong to. - entities: The entities, local to the process. - - Returns: - The integration entities, flattened, as returned by DOLFINx. - """ - try: - return dolfinx.fem.compute_integration_domains(integral_type, topology, entities) - except TypeError: - dim = topology.dim if integral_type == dolfinx.fem.IntegralType.cell else topology.dim - 1 - return dolfinx.fem.compute_integration_domains(integral_type, topology, entities, dim) - - -def get_facet_permutations(topology: dolfinx.mesh.Topology) -> npt.NDArray[np.uint8]: - """The permutation of every facet of every cell, across DOLFINx versions. - - Each value encodes how a facet is oriented as seen from a cell, relative to a low-to-high - ordering of its global vertex indices, as FFCx uses for ``quadrature_permutation``. - - Args: - topology: The topology to compute the facet permutations of. - - Returns: - The permutations, shape ``(num_cells, num_facets_per_cell)``, ghost cells included. - """ - fdim = topology.dim - 1 - num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(topology.cell_type, fdim) - if hasattr(topology, "create_cell_permutations"): - topology.create_entity_permutations(fdim) - permutations = topology.get_entity_permutations(fdim) - else: # DOLFINx < 0.10 - topology.create_entity_permutations() # type: ignore[call-arg] - permutations = topology.get_facet_permutations() - return np.asarray(permutations).reshape(-1, num_facets_per_cell) - - -def create_cpp_finite_element(constructor, element, gdim: int, block_shape: tuple[int, ...]): - """Create a blocked C++ finite element across the supported DOLFINx versions. - - FEniCS/dolfinx#4511 added the geometric dimension to the constructor. Versions from before - ``block_shape`` existed raise ``TypeError``, for the caller to fall back on ``block_size``. - - Args: - constructor: ``dolfinx.cpp.fem.FiniteElement_float32`` or ``_float64``. - element: The C++ Basix element. - gdim: Geometric dimension of the mesh. - block_shape: Block shape of the element. - """ - try: - return constructor(element, gdim=gdim, block_shape=block_shape, symmetric=False) - except TypeError: - return constructor(element, block_shape=block_shape, symmetric=False) diff --git a/src/scifem/eval.py b/src/scifem/eval.py index 582b3fd..512cd84 100644 --- a/src/scifem/eval.py +++ b/src/scifem/eval.py @@ -10,7 +10,7 @@ import ufl from scipy.optimize import minimize from ufl.algorithms.signature import compute_expression_signature -from . import compat +import fenicsx_compat T = typing.TypeVar("T", int, float) @@ -121,7 +121,7 @@ def find_cell_extrema( jit_options: dict | None = None, method: str | None = None, options: dict | None = None, - tol: float | None = None, + tol: float | np.floating | None = None, ) -> tuple[npt.NDArray[np.floating], np.floating]: """ Find the extrema of a :py:class:`ufl.core.expr.Expr` within a cell. @@ -148,7 +148,7 @@ def find_cell_extrema( uds = ufl.domain.extract_domains(u) if uds is None or len(uds) != 1: raise ValueError(f"Expression {u:s} must be defined on a single domain") - ud = uds[0] + ud = typing.cast(ufl.Mesh, uds[0]) mesh = dolfinx.mesh.Mesh(ud.ufl_cargo(), ud) if x0 is None: @@ -157,9 +157,9 @@ def find_cell_extrema( x_ref = x0 _cell = np.array([cell], dtype=np.int32) - mesh_nodes = mesh.geometry.x[compat.dofmap(mesh)[cell], : mesh.geometry.dim] + mesh_nodes = mesh.geometry.x[fenicsx_compat.dofmap(mesh)[cell], : mesh.geometry.dim] _x_p = np.zeros(3) - cmap = compat.cmap(mesh) + cmap = fenicsx_compat.cmap(mesh) def eval_J(x_ref): # Evaluating basis functions through :py:func:`dolfinx.fem.Function.eval` @@ -312,7 +312,7 @@ def compute_extrema( uds = ufl.domain.extract_domains(u) if uds is None or len(uds) != 1: raise ValueError(f"Expression {u:s} must be defined on a single domain") - ud = uds[0] + ud = typing.cast(ufl.Mesh, uds[0]) mesh = dolfinx.mesh.Mesh(ud.ufl_cargo(), ud) if kind is builtins.max: @@ -348,7 +348,7 @@ def compute_extrema( with tempfile.TemporaryDirectory(ignore_cleanup_errors=False, delete=True) as cache_dir: jit_options["cache_dir"] = cache_dir # Find local extrema - local_min = np.inf + local_min: float | np.floating = np.inf local_X = None for cell in candidate_cells: X_c, extrema = find_cell_extrema( diff --git a/src/scifem/geometry.py b/src/scifem/geometry.py index a06f64b..c464a71 100644 --- a/src/scifem/geometry.py +++ b/src/scifem/geometry.py @@ -3,7 +3,7 @@ import numpy.typing as npt import warnings from scifem._scifem import closest_point_projection_float64, closest_point_projection_float32 -from . import compat +import fenicsx_compat def project_onto_simplex( @@ -159,10 +159,14 @@ def _closest_point_projection( roundoff_tol = 100 * eps # Extract scalar element of mesh - element = mesh.ufl_domain().ufl_coordinate_element().sub_elements[0] + ufl_domain = mesh.ufl_domain() + assert ufl_domain is not None + element = ufl_domain.ufl_coordinate_element().sub_elements[0] tdim = mesh.topology.dim # Get the coordinates of the nodes for the specified cell - node_coords = mesh.geometry.x[compat.dofmap(mesh)[cells]][:, :, : mesh.geometry.dim] + sub_xdofs = fenicsx_compat.dofmap(mesh)[cells] + # Output form is (num_cells, num_xnodes_per_cell, gdim) + node_coords = mesh.geometry.x[sub_xdofs][:, :, : mesh.geometry.dim] target_points = target_points.reshape(-1, mesh.geometry.dim) # Constraints and Bounds @@ -219,7 +223,8 @@ def project(x): # for stretching of the reference space jac_norm = np.linalg.norm(tangents) - scaled_tol_grad = tol_grad * max(jac_norm, 1.0) + # mypy infers max(np.floating, float) as a comparison protocol, not a number + scaled_tol_grad = tol_grad * max(jac_norm, 1.0) # type: ignore[operator] if np.linalg.norm(g) < scaled_tol_grad: break diff --git a/src/scifem/interpolation.py b/src/scifem/interpolation.py index 311ae28..c167df2 100644 --- a/src/scifem/interpolation.py +++ b/src/scifem/interpolation.py @@ -2,16 +2,22 @@ from mpi4py import MPI -from packaging.version import Version import basix import dolfinx import ufl import numpy as np import numpy.typing as npt import scipy.linalg -from .bcs import build_quadrature_permutations, pull_back_to_reference_facet -from .compat import create_cell_permutations, get_facet_permutations -from .ufl_compat import apply_pullback_inverse +from fenicsx_compat import ( + apply_pullback_inverse, + before, + cell_permutation_info, + facet_permutations, + interpolate, + interpolate_to_submesh_entity_maps, + permute_facet_quadrature, +) +from .bcs import pull_back_to_reference_facet from .utils import group_by_key, unroll_dofmap from .mesh import _EntityMap @@ -94,10 +100,7 @@ def prepare_interpolation_data( raise RuntimeError("Only codim-1 interpolation matrices can be defined") # Extract quadrature points for expression (in Q space) - try: - q_points = Q.element.interpolation_points() - except TypeError: - q_points = Q.element.interpolation_points + q_points = Q.element.interpolation_points # Compile expression num_points = q_points.shape[0] @@ -123,7 +126,7 @@ def prepare_interpolation_data( except (RuntimeError, ValueError): Q_vs = 1 # If we do not have a basix element, assume value size is 1 assert isinstance(Q.ufl_element().pullback, ufl.pullback.IdentityPullback) - pull_back = lambda x: None + pull_back = lambda x: None # type: ignore[assignment] # Unused for identity pullback assert Q.element.interpolation_ident im = None @@ -147,8 +150,13 @@ def prepare_interpolation_data( for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): for q in range(Q.dofmap.bs): + # The scalar type of a dolfinx.fem.Expression is not inferred statically, + # so its evaluations are typed float32 regardless of the element's type new_array[:, q * Q_vs : (q + 1) * Q_vs, i] = pull_back( - array_evaluated[:, :, q * Q_vs : (q + 1) * Q_vs, i], jacs, detJs, Ks + array_evaluated[:, :, q * Q_vs : (q + 1) * Q_vs, i], # type: ignore[arg-type] + jacs, # type: ignore[arg-type] + detJs, # type: ignore[arg-type] + Ks, # type: ignore[arg-type] ).reshape(num_cells * num_points, Q_vs) new_array = new_array.reshape( num_cells, num_points, Q.dofmap.bs * Q_vs, V.dofmap.bs * V.dofmap.dof_layout.num_dofs @@ -187,6 +195,7 @@ def prepare_interpolation_data( for c in range(num_cells): for i in range(V.dofmap.bs * V.dofmap.dof_layout.num_dofs): tmp_array = np.zeros((int(num_points), Q.dofmap.bs * Q_vs), dtype=np.float64) + assert im is not None # Only None for an identity interpolation for p in range(num_points): tmp_array[p] = new_array[c, p, :, i] if Q.dofmap.bs == 1: @@ -199,8 +208,7 @@ def prepare_interpolation_data( if Q.element.needs_dof_transformations: # Apply dof transformation to each column (using Piola maps) - create_cell_permutations(mesh.topology) - cell_perm = mesh.topology.get_cell_permutation_info()[:num_cells] + cell_perm = cell_permutation_info(mesh.topology)[:num_cells] permuted_matrix = interpolated_matrix.flatten().copy() Q.element.Tt_inv_apply( @@ -247,7 +255,7 @@ def interpolation_matrix( def scatter( A: dolfinx.la.MatrixCSR, num_cells: int, - dofs_visited: npt.NDArray[np.int32], + dofs_visited: npt.NDArray[np.int8], num_rows_local: int, array_evaluated: npt.NDArray[np.inexact], dofmap0: npt.NDArray[np.int32], @@ -317,7 +325,7 @@ def petsc_interpolation_matrix( def scatter( A: PETSc.Mat, num_cells: int, - dofs_visited: npt.NDArray[np.int32], + dofs_visited: npt.NDArray[np.int8], num_rows_local: int, array_evaluated: npt.NDArray[np.inexact], dofmap0: npt.NDArray[np.int32], @@ -330,7 +338,13 @@ def scatter( A_local = array_evaluated[i].reshape(len(rows), len(cols)) row_filter = (dofs_visited[rows] == 1) | (rows >= num_rows_local) A_local[row_filter] = 0 - A.setValuesLocal(rows, cols, A_local, addv=PETSc.InsertMode.ADD_VALUES) + # petsc4py annotates the indices as Sequence[int], numpy arrays are accepted + A.setValuesLocal( + rows, # type: ignore[arg-type] + cols, # type: ignore[arg-type] + A_local, + addv=PETSc.InsertMode.ADD_VALUES, + ) dofs_visited[rows] = 1 row_dofmap = unroll_dofmap(Q.dofmap.list, Q.dofmap.bs) # (num_cells, num_rows) @@ -380,8 +394,7 @@ def compute_entity_closure_permutations( mesh = V.mesh element = V.element.basix_element entity_types = basix.cell.subentity_types(element.cell_type)[dim] - create_cell_permutations(mesh.topology) - cell_info = mesh.topology.get_cell_permutation_info() + cell_info = cell_permutation_info(mesh.topology) entities = np.asarray(entities, dtype=np.int32).reshape(-1, 2) cells, local_entities = entities.T if size is None: @@ -466,9 +479,7 @@ def __init__( integration_entities: npt.NDArray[np.int32], entity_maps: list[_EntityMap] | None = None, ): - if Version(dolfinx.__version__) < Version("0.10.0"): - raise RuntimeError("interpolate_to_submesh requires dolfinx version 0.10.0 or higher") - ufl_domains = ufl.domain.extract_domains(expr) + ufl_domains = typing.cast(tuple[ufl.Mesh, ...], ufl.domain.extract_domains(expr)) max_tdim_pos = np.argmax([domain.ufl_cargo().topology.dim for domain in ufl_domains]) self._mesh = dolfinx.mesh.Mesh( ufl_domains[max_tdim_pos].ufl_cargo(), ufl_domains[max_tdim_pos] @@ -479,19 +490,17 @@ def __init__( submesh = V_surface.mesh ip = V_surface.element.interpolation_points - try: - self._expression = dolfinx.fem.Expression(expr, ip, entity_maps=entity_maps) - except TypeError: - self._expression = dolfinx.fem.Expression(expr, ip) + # fenicsx_compat annotates the first argument as a Function, but any expression is accepted + self._expression = interpolate_to_submesh_entity_maps(expr, ip, entity_maps) # type: ignore[arg-type] self._mesh.topology.create_connectivity(self._mesh.topology.dim, submesh.topology.dim) self._mesh.topology.create_connectivity(submesh.topology.dim, self._mesh.topology.dim) - create_cell_permutations(self._mesh.topology) - create_cell_permutations(submesh.topology) + cell_permutation_info(self._mesh.topology) + cell_permutation_info(submesh.topology) # Before the introduction of https://github.com/FEniCS/dolfinx/pull/4140 # one needed to permute the data according to the facet permutations, # one permutation per evaluation point. self._permutations: npt.NDArray[np.int32] | None = None - if Version(dolfinx.__version__) < Version("0.11.0.dev0"): + if before("0.11.0.dev0"): V_vol = expr.function_space # type: ignore[attr-defined] self._permutations = compute_entity_closure_permutations( V_vol, V_vol.mesh.topology.dim - 1, integration_entities, ip.shape[0] @@ -525,16 +534,11 @@ def apply(self, u_surface: dolfinx.fem.Function): data = data.swapaxes(0, 1) if self._expression.value_size == 1: - shaped_data = data.flatten() + shaped_data: npt.NDArray[np.inexact] = data.flatten() else: shaped_data = data.reshape(self._expression.value_size, -1) - if hasattr(u_surface._cpp_object, "interpolate_f"): - interpolate_func = u_surface._cpp_object.interpolate_f - else: - interpolate_func = u_surface._cpp_object.interpolate - - interpolate_func(shaped_data, self._submesh_facets) + interpolate(u_surface._cpp_object, shaped_data, self._submesh_facets) u_surface.x.scatter_forward() @@ -735,14 +739,13 @@ def _setup_piola( entities = np.asarray(integration_entities, dtype=np.int32).reshape(-1, 2) cells, local_facets = entities.T - create_cell_permutations(mesh.topology) - cell_info = mesh.topology.get_cell_permutation_info()[cells] - permutations = get_facet_permutations(mesh.topology)[cells, local_facets] + cell_info = cell_permutation_info(mesh.topology)[cells] + permutations = facet_permutations(mesh.topology)[cells, local_facets] groups = group_by_key(permutations.astype(np.int64) * num_facets_per_cell + local_facets) v = ufl.TestFunction(V_surface) expr = apply_pullback_inverse(V_volume.ufl_element().pullback, v, mesh.ufl_domain()) - point_sets = build_quadrature_permutations(facet_type, reference_points) + point_sets = permute_facet_quadrature(facet_type, reference_points) # Compiled on each process for its own permutations only expressions = { perm: dolfinx.fem.Expression( diff --git a/src/scifem/mesh.py b/src/scifem/mesh.py index f2301fc..6b06a62 100644 --- a/src/scifem/mesh.py +++ b/src/scifem/mesh.py @@ -7,19 +7,16 @@ from typing_extensions import deprecated from . import _scifem # type: ignore -from .compat import compute_integration_domains, create_cpp_finite_element import collections import dolfinx import typing import numpy as np import numpy.typing as npt from typing import Protocol -from packaging.version import Version import basix +import fenicsx_compat import ufl -from . import compat - __all__ = [ "create_entity_markers", "transfer_meshtags_to_submesh", @@ -48,9 +45,12 @@ class _EntityMap(Protocol): """Protocol for EntityMap-like objects.""" - sub_topology: dolfinx.mesh.Topology - dim: int - _cpp_object: object + @property + def sub_topology(self) -> dolfinx.mesh.Topology: ... + @property + def dim(self) -> int: ... + @property + def _cpp_object(self) -> object: ... def sub_topology_to_topology( self, entities: npt.NDArray[np.int32], inverse: bool @@ -73,13 +73,15 @@ def get_entity_map( Mapped indices of entities. """ try: - sub_top = entity_map.sub_topology + em = typing.cast(_EntityMap, entity_map) + sub_top = em.sub_topology assert isinstance(sub_top, dolfinx.mesh.Topology) - sub_map = sub_top.index_map(entity_map.dim) + sub_map = sub_top.index_map(em.dim) indices = np.arange(sub_map.size_local + sub_map.num_ghosts, dtype=np.int32) - return entity_map.sub_topology_to_topology(indices, inverse=False) + return em.sub_topology_to_topology(indices, inverse=False) except AttributeError: - return entity_map + # Deprecated array-based entity map (older DOLFINx) + return typing.cast(npt.NDArray[np.int32], entity_map) def create_entity_markers( @@ -207,8 +209,6 @@ def reverse_mark_entities( Returns: Local indices marked on any process sharing this entity """ - if hasattr(dolfinx.common, "index_map"): - entity_map = entity_map comm_vec = dolfinx.la.vector(entity_map, dtype=np.int32) comm_vec.array[:] = 0 comm_vec.array[entities] = 1 @@ -249,11 +249,11 @@ def extract_submesh( submesh, cell_map, vertex_map, node_map = dolfinx.mesh.create_submesh(mesh, edim, entities) # Transfer cell markers - if hasattr(dolfinx.mesh, "transfer_meshtags_to_submesh"): - new_et = dolfinx.mesh.transfer_meshtags_to_submesh( - entity_tag, submesh, vertex_to_parent=vertex_map, cell_to_parent=cell_map + try: + new_et = fenicsx_compat.transfer_meshtags_to_submesh( + entity_tag, submesh, cell_map, vertex_map ) - else: + except NotImplementedError: # DOLFINx 0.10 has no native version new_et, _ = transfer_meshtags_to_submesh( entity_tag, submesh, vertex_to_parent=vertex_map, cell_to_parent=cell_map ) @@ -276,11 +276,7 @@ def find_interface( Returns: The facets shared between the two domains. """ - # Wrapping C++ objects to be back compatible with older versions of dolfinx - if not isinstance(cell_tags.topology, dolfinx.cpp.mesh.Topology): - topology = cell_tags.topology - else: - topology = dolfinx.mesh.Topology(cell_tags.topology) + topology = fenicsx_compat.meshtags_topology(cell_tags) assert topology.dim == cell_tags.dim tdim = topology.dim @@ -352,7 +348,7 @@ def compute_subdomain_exterior_facets( sub_facets = dolfinx.mesh.exterior_facet_indices(sub_mesh.topology) # Map exterior facet to (submesh_cell, local_facet_index) tuples - integration_entities = compute_integration_domains( + integration_entities = dolfinx.fem.compute_integration_domains( dolfinx.fem.IntegralType.exterior_facet, sub_mesh.topology, sub_facets ).reshape(-1, 2) submap_array = get_entity_map(cell_map) @@ -393,13 +389,6 @@ def compute_interface_data( Returns: The integration data. """ - # Future compatibility check - integration_args: tuple[int] | tuple - if Version("0.10.0") <= Version(dolfinx.__version__): - integration_args = () - else: - fdim = cell_tags.dim - 1 - integration_args = (fdim,) if include_ghosts: if len(facet_indices) == 0: return np.empty((0, 4), dtype=np.int32) @@ -432,15 +421,10 @@ def compute_interface_data( assert np.unique(_row).shape[0] == len(_row) idata = np.vstack([cells, local_pos]).T.reshape(-1, 4) else: - if not isinstance(cell_tags.topology, dolfinx.cpp.mesh.Topology): - _cpp_topology = cell_tags.topology._cpp_object - else: - _cpp_topology = cell_tags.topology idata = dolfinx.cpp.fem.compute_integration_domains( dolfinx.fem.IntegralType.interior_facet, - _cpp_topology, + fenicsx_compat.meshtags_topology(cell_tags)._cpp_object, facet_indices, - *integration_args, ) ordered_idata = idata.reshape(-1, 4).copy() switch = cell_tags.values[ordered_idata[:, 0]] > cell_tags.values[ordered_idata[:, 2]] @@ -462,7 +446,7 @@ def create_geometry_function_space( """ geom_imap = mesh.geometry.index_map() - geom_dofmap = compat.dofmap(mesh) + geom_dofmap = fenicsx_compat.dofmap(mesh) ufl_domain = mesh.ufl_domain() assert ufl_domain is not None sub_el = ufl_domain.ufl_coordinate_element().sub_elements[0] @@ -480,6 +464,12 @@ def create_geometry_function_space( ufl_el = basix.ufl.blocked_element(sub_el, shape=(N,)) value_shape = (N,) + _fe_constructor: ( + type[dolfinx.cpp.fem.FiniteElement_float32] | type[dolfinx.cpp.fem.FiniteElement_float64] + ) + _fem_constructor: ( + type[dolfinx.cpp.fem.FunctionSpace_float32] | type[dolfinx.cpp.fem.FunctionSpace_float64] + ) if ufl_el.dtype == np.float32: _fe_constructor = dolfinx.cpp.fem.FiniteElement_float32 _fem_constructor = dolfinx.cpp.fem.FunctionSpace_float32 @@ -488,24 +478,20 @@ def create_geometry_function_space( _fem_constructor = dolfinx.cpp.fem.FunctionSpace_float64 else: raise RuntimeError(f"Unsupported type {ufl_el.dtype}") - try: - cpp_el = create_cpp_finite_element( - _fe_constructor, ufl_el.basix_element._e, mesh.geometry.dim, value_shape - ) - except TypeError: - cpp_el = _fe_constructor(ufl_el.basix_element._e, block_size=N, symmetric=False) + cpp_el = fenicsx_compat.finite_element_ctor_kwargs( + _fe_constructor, + ufl_el.basix_element._e, + gdim=mesh.geometry.dim, + value_shape=value_shape, + ) dof_layout = dolfinx.cpp.fem.create_element_dof_layout(cpp_el, []) - if not isinstance(geom_imap, dolfinx.cpp.common.IndexMap): - _cpp_im = geom_imap._cpp_object - else: - _cpp_im = geom_imap - cpp_dofmap = dolfinx.cpp.fem.DofMap(dof_layout, _cpp_im, N, adj_list, N) + cpp_dofmap = dolfinx.cpp.fem.DofMap( + dof_layout, fenicsx_compat.unwrap_index_map(geom_imap), N, adj_list, N + ) - # Create function space - try: - cpp_space = _fem_constructor(mesh._cpp_object, cpp_el, cpp_dofmap) - except TypeError: - cpp_space = _fem_constructor(mesh._cpp_object, cpp_el, cpp_dofmap, value_shape=value_shape) + # Create function space. The constructor matches the mesh dtype at runtime, + # which the float32/float64 stubs cannot express. + cpp_space = _fem_constructor(mesh._cpp_object, cpp_el, cpp_dofmap) # type: ignore[arg-type] return dolfinx.fem.FunctionSpace(mesh, ufl_el, cpp_space) @@ -551,9 +537,13 @@ def copy(mesh: dolfinx.mesh.Mesh) -> dolfinx.mesh.Mesh: The new mesh """ - # Initialize C++ mesh + # Initialize C++ mesh (geometry dtype matches the mesh class at runtime) cpp_mesh = type(mesh._cpp_object)( - mesh.comm, mesh.topology._cpp_object, mesh.geometry._cpp_object + mesh.comm, + mesh.topology._cpp_object, + mesh.geometry._cpp_object, # type: ignore[arg-type] ) # Initialize Python wrapper with new symbolic numbering in the `ufl.Mesh` - return dolfinx.mesh.Mesh(cpp_mesh, ufl.Mesh(mesh.ufl_domain().ufl_coordinate_element())) + ufl_domain = mesh.ufl_domain() + assert ufl_domain is not None + return dolfinx.mesh.Mesh(cpp_mesh, ufl.Mesh(ufl_domain.ufl_coordinate_element())) diff --git a/src/scifem/mpi_utils.py b/src/scifem/mpi_utils.py index b6dd49f..404d7ab 100644 --- a/src/scifem/mpi_utils.py +++ b/src/scifem/mpi_utils.py @@ -2,6 +2,7 @@ import numpy as np import numpy.typing as npt import dolfinx +import fenicsx_compat mpi_dtype = { np.float64: MPI.DOUBLE, @@ -199,14 +200,10 @@ def broadcast_marked_entities(mesh, dim, entities): def local_range(comm, num_indices): """The block of ``range(num_indices)`` this process is the post office for. - Uses ``dolfinx.common.local_range`` where it exists, so that the blocks match what the - rest of DOLFINx means by the same words; the fallback reproduces it. + The same blocks as ``dolfinx.common.local_range``, so that they match what the rest of + DOLFINx means by the same words. """ - if hasattr(dolfinx.common, "local_range"): - return tuple(dolfinx.common.local_range(comm.rank, int(num_indices), comm.size)) - per_rank, remainder = divmod(int(num_indices), comm.size) - low = comm.rank * per_rank + min(comm.rank, remainder) - return low, low + per_rank + (1 if comm.rank < remainder else 0) + return fenicsx_compat.local_range(comm.rank, int(num_indices), comm.size) def index_owner(comm, indices, num_indices): diff --git a/src/scifem/periodic/gmsh.py b/src/scifem/periodic/gmsh.py index 8c00edf..e5670ac 100644 --- a/src/scifem/periodic/gmsh.py +++ b/src/scifem/periodic/gmsh.py @@ -20,77 +20,17 @@ from __future__ import annotations -import typing -from mpi4py import MPI as _MPI -import dolfinx +import fenicsx_compat import numpy as np from .mesh import create_periodic_mesh_from_igi, DEFAULT_TAG_BASE from .utils import PeriodicNodes, resolve_to_roots -import inspect -from ..compat import create_partitioner __all__ = [ "extract_gmsh_periodic_nodes", "read_periodic_mesh_from_msh", - "model_to_mesh", ] -def model_to_mesh( - model, - comm, - rank: int = 0, - gdim: int = 3, - max_facet_to_cell_links: int = 2, - ghost_mode: dolfinx.mesh.GhostMode = dolfinx.mesh.GhostMode.shared_facet, - partitioner: typing.Callable[ - [_MPI.Comm, int, int, dolfinx.cpp.graph.AdjacencyList_int32], - dolfinx.cpp.graph.AdjacencyList_int32, - ] - | None = None, -): - """``model_to_mesh`` with shared-facet ghosting, wherever the version wants it told. - - On 0.12.0.dev0, ``model_to_mesh`` takes `ghost_mode` directly and threads it into whatever - partitioner is used -- the caller's, or its own default if `partitioner` is `None` -- - at call time. On 0.11 there is no such keyword, so `ghost_mode` can only reach the - mesh through a partitioner already built with it, which is constructed here unless - the caller supplied one. Getting this wrong does not fail here -- it fails much - later, when an interior facet integral finds an interprocess facet with only one - cell. - - Args: - model: An initialised, meshed ``gmsh.model``. - comm: The communicator to distribute the mesh over. - rank: The rank that reads the model. - gdim: Geometric dimension of the mesh. - max_facet_to_cell_links: Passed to `model_to_mesh` where it takes it directly, - and to the partitioner constructed here otherwise. - ghost_mode: The ghost mode the mesh must be built with; see - :py:func:`scifem.periodic.mesh.check_facet_ghosting` for why the rebuild - needs `shared_facet`. - partitioner: A caller-supplied partitioner, which takes priority over the one - built here from `ghost_mode`. On 0.11, supplying one means `ghost_mode` and - `max_facet_to_cell_links` are not applied, since there both only take effect - through the partitioner. - - Returns: - The mesh. - """ - params = inspect.signature(dolfinx.io.gmsh.model_to_mesh).parameters - kwargs: dict[str, typing.Any] = {"partitioner": partitioner} - if "ghost_mode" in params: - kwargs["ghost_mode"] = ghost_mode - if "max_facet_to_cell_links" in params: - kwargs["max_facet_to_cell_links"] = max_facet_to_cell_links - elif partitioner is None: - kwargs["partitioner"] = create_partitioner( - ghost_mode, max_facet_to_cell_links=max_facet_to_cell_links - ) - mesh_data = dolfinx.io.gmsh.model_to_mesh(model, comm, rank=rank, gdim=gdim, **kwargs) - return getattr(mesh_data, "mesh", mesh_data) - - def extract_gmsh_periodic_nodes( model, include_high_order: bool = False, tol: float | None = 1e-8 ) -> PeriodicNodes: @@ -192,10 +132,12 @@ def read_periodic_mesh_from_msh( comm: The communicator to distribute the mesh over. rank: The rank that reads the file. gdim: Geometric dimension of the mesh. - partitioner: Cell partitioner, passed through to ``model_to_mesh``. + partitioner: Cell partitioner, passed through to + :py:func:`fenicsx_compat.model_to_mesh`. tag_base: Passed through to :py:func:`create_periodic_mesh_from_igi`. - kwargs: Further arguments for ``model_to_mesh``, such as ``ghost_mode`` where the - installed DOLFINx takes it there. + kwargs: Further arguments for :py:func:`fenicsx_compat.model_to_mesh`, such as + ``ghost_mode``. Its default, ``shared_facet``, is what the rebuild needs; see + :py:func:`scifem.periodic.mesh.check_facet_ghosting`. Returns: ``(periodic_mesh, replaced_vertices, replacement_map)``, as @@ -218,17 +160,16 @@ def read_periodic_mesh_from_msh( gmsh.merge(str(filename)) pairs = extract_gmsh_periodic_nodes(gmsh.model) else: - pairs = PeriodicNodes(0) - mesh_data = model_to_mesh( - gmsh.model, comm, rank=rank, gdim=gdim, partitioner=partitioner, **kwargs + pairs = PeriodicNodes() + mesh_data = fenicsx_compat.model_to_mesh( + gmsh.model, comm, rank, gdim=gdim, partitioner=partitioner, **kwargs ) finally: if started_here and gmsh.isInitialized(): gmsh.finalize() - mesh = getattr(mesh_data, "mesh", mesh_data) return create_periodic_mesh_from_igi( - mesh, + mesh_data.mesh, pairs.replaced, pairs.partner, pairs.num_nodes_global, diff --git a/src/scifem/periodic/mesh.py b/src/scifem/periodic/mesh.py index 205b5a6..5cfc161 100644 --- a/src/scifem/periodic/mesh.py +++ b/src/scifem/periodic/mesh.py @@ -9,6 +9,13 @@ import numpy.typing as npt import dataclasses +from fenicsx_compat import ( + create_index_map, + create_sub_index_map, + index_to_dest_ranks, + unwrap_index_map, +) + from ..mpi_utils import ( broadcast_marked_entities, get_ownership, @@ -20,8 +27,6 @@ mpi_dtype, exchange_to_destinations, ) -from ..compat import index_map, topology as compat_topology -from ..compat import ghosting_ranks as _compat_ghosting_ranks from .utils import VertexCorrespondence, PeriodicNodes from .topological_search import periodic_correspondence_from_nodes from .geometrical_search import match_vertices_geometric @@ -121,22 +126,10 @@ def _reduced_vertex_map(mesh, indicator_vertices): keep_vertices = np.ones(num_vertices_local, dtype=np.bool_) keep_vertices[removed] = False reduced_vertices = np.flatnonzero(keep_vertices).astype(np.int32) - # Compat: 0.12 moved `create_sub_index_map` to `dolfinx.common` and made it report - # ownership changes instead of taking a flag. Once only that form is supported the - # branch goes away and this function is the four lines around the call. - if hasattr(dolfinx.common, "create_sub_index_map"): - sub_map, sub_to_parent, changed_owner = dolfinx.common.create_sub_index_map( - vertex_map, reduced_vertices - ) - if vertex_map.comm.allreduce(changed_owner, op=MPI.LOR): - raise RuntimeError( - "Vertex ownership has changed, which is not supported. " - "Please report this issue to the dolfinx developers." - ) - else: - sub_map, sub_to_parent = dolfinx.cpp.common.create_sub_index_map( - vertex_map, reduced_vertices, allow_owner_change=False - ) + # Vertex ownership must not change; this raises on every process if it would. + sub_map, sub_to_parent = create_sub_index_map( + vertex_map, reduced_vertices, allow_owner_change=False + ) parent_to_sub = np.full(num_vertices_local, -1, dtype=np.int32) parent_to_sub[sub_to_parent] = np.arange(sub_to_parent.size, dtype=np.int32) @@ -387,7 +380,7 @@ def _holders_of_partner_vertices(mesh, partner_vertex, dest_owner, tag: int): # --- from the owner to every holder. `index_to_dest_ranks` lists the ranks that ghost # each owned index; the owner holds it as well, and is appended as itself. - ghosting_ranks, ghosting_offsets = _compat_ghosting_ranks(vertex_map, tag) + ghosting_ranks, ghosting_offsets = index_to_dest_ranks(vertex_map, tag) local = vertex_map.global_to_local(np.ascontiguousarray(at_owner[:, 0])) assert (local != -1).all(), "A vertex was reported to a rank that does not own it" # Flatten the ragged per-vertex rank lists. Equivalent to: @@ -899,7 +892,9 @@ def _build_periodic_mesh( # Vertex map is temporary, as we need to extend it with additional ghosts # on the process taking over facets - tmp_vertex_map = index_map(comm, new_local_size, new_ghosts, new_owners, tag=tag_vertex_ghosts) + tmp_vertex_map = create_index_map( + comm, new_local_size, new_ghosts, new_owners, tag=tag_vertex_ghosts + ) tmp_vertex_ownership = get_ownership(tmp_vertex_map) # Create replacement map @@ -1058,8 +1053,8 @@ def _build_periodic_mesh( lost_geom_igi = mesh.geometry.input_global_indices[org_geom_dm_cells_losing_vertex] # Compute insertion position based on cell ownership - lost_src_ranks, num_send_lost_cells = np.unique(lost_dest_ranks, return_counts=True) - num_send_lost_cells = num_send_lost_cells.astype(np.int32) + lost_src_ranks, lost_counts = np.unique(lost_dest_ranks, return_counts=True) + num_send_lost_cells = lost_counts.astype(np.int32) lost_cell_insert_pos = compute_insert_position( lost_dest_ranks, lost_src_ranks, num_send_lost_cells ) @@ -1103,8 +1098,8 @@ def _build_periodic_mesh( lost_cells_to_gainer_comm = comm.Create_dist_graph( [comm.rank], [len(lost_src_ranks)], lost_src_ranks.tolist(), reorder=False ) - recv_lost_cells_ranks, sent_to_ranks, _ = lost_cells_to_gainer_comm.Get_dist_neighbors() - recv_lost_cells_ranks = np.asarray(recv_lost_cells_ranks, dtype=np.int32) + recv_ranks, sent_to_ranks, _ = lost_cells_to_gainer_comm.Get_dist_neighbors() + recv_lost_cells_ranks = np.asarray(recv_ranks, dtype=np.int32) assert np.array_equal(np.asarray(sent_to_ranks, dtype=np.int32), lost_src_ranks) num_recv_lost_cells = np.zeros(len(recv_lost_cells_ranks), dtype=np.int32) all_to_all(lost_cells_to_gainer_comm, num_send_lost_cells, num_recv_lost_cells) @@ -1311,34 +1306,31 @@ def _build_periodic_mesh( assert (all_owners != comm.rank).all(), "Ghosted vertices on owned process" # Create new cell and vertex map - new_cell_map = index_map( + new_cell_map = create_index_map( comm, cell_map.size_local, all_cell_ghosts, all_cell_owners, tag=tag_cell_map ) - new_vertex_map = index_map( + new_vertex_map = create_index_map( comm, tmp_vertex_map.size_local, all_ghosts, all_owners, tag=tag_vertex_map ) new_c_to_v = dolfinx.graph.adjacencylist( np.vstack([new_c, extra_dm, lost_cells_dofs_as_local.reshape(-1, num_vertices)]) ) - new_v_to_v = dolfinx.graph.adjacencylist( - np.arange(new_vertex_map.size_local + new_vertex_map.num_ghosts, dtype=np.int32) - ) assert (new_c_to_v.array < new_vertex_map.size_local + new_vertex_map.num_ghosts).all(), ( "Cell to vertex map is out of bounds" ) - topology = compat_topology( - comm, + topology = dolfinx.cpp.mesh.Topology( mesh.topology.cell_type, - mesh.topology.dim, - new_vertex_map, - new_cell_map, - new_c_to_v, - new_v_to_v, + unwrap_index_map(new_vertex_map), + unwrap_index_map(new_cell_map), + new_c_to_v._cpp_object, # type: ignore[arg-type] all_cell_oci, ) - c_el = dolfinx.fem.coordinate_element(mesh._ufl_domain.ufl_coordinate_element().basix_element) + # `coordinate_element` is a singledispatch; the basix-element overload takes no degree + c_el = dolfinx.fem.coordinate_element( # type: ignore[call-arg] + mesh._ufl_domain.ufl_coordinate_element().basix_element + ) assert ( (all_ghosts < tmp_vertex_map.local_range[0]) | (tmp_vertex_map.local_range[1] <= all_ghosts) @@ -1366,7 +1358,7 @@ def _build_periodic_mesh( extended_coords = np.vstack( [mesh.geometry.x, extra_node_coords, filtered_geometry_coords] ).astype(mesh.geometry.x.dtype)[:, : mesh.geometry.dim] - new_node_im = index_map( + new_node_im = create_index_map( comm, num_local_nodes, extended_geom_ghosts, @@ -1381,10 +1373,12 @@ def _build_periodic_mesh( geometry = dolfinx.mesh.create_geometry( new_node_im, extended_dofmap, c_el, extended_coords, extended_igi ) + # The geometry dtype matches the mesh class chosen at runtime + cpp_mesh: dolfinx.cpp.mesh.Mesh_float32 | dolfinx.cpp.mesh.Mesh_float64 if mesh.geometry.x.dtype == np.float64: - cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(comm, topology, geometry._cpp_object) + cpp_mesh = dolfinx.cpp.mesh.Mesh_float64(comm, topology, geometry._cpp_object) # type: ignore[arg-type] elif mesh.geometry.x.dtype == np.float32: - cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(comm, topology, geometry._cpp_object) + cpp_mesh = dolfinx.cpp.mesh.Mesh_float32(comm, topology, geometry._cpp_object) # type: ignore[arg-type] else: raise RuntimeError(f"Unsupported dtype for mesh {mesh.geometry.x.dtype}") diff --git a/src/scifem/petsc.py b/src/scifem/petsc.py index 3f51427..9494fac 100644 --- a/src/scifem/petsc.py +++ b/src/scifem/petsc.py @@ -1,144 +1,68 @@ -import dolfinx +from __future__ import annotations + import typing +import dolfinx +import fenicsx_compat +from fenicsx_compat import ghost_update, zero_petsc_vector + +if typing.TYPE_CHECKING: + from petsc4py import PETSc + __all__ = [ "zero_petsc_vector", "ghost_update", "apply_lifting_and_set_bc", ] -if dolfinx.has_petsc4py: - from petsc4py import PETSc - - try: - import dolfinx.la.petsc - - def zero_petsc_vector(b: PETSc.Vec) -> None: - """Zero a PETSc vector, including ghosts""" - dolfinx.la.petsc._zero_vector(b) - - def ghost_update( - x: PETSc.Vec, - insert_mode: PETSc.InsertMode, - scatter_mode: PETSc.ScatterMode, - ) -> None: - """Ghost update a vector""" - dolfinx.la.petsc._ghost_update(x, insert_mode, scatter_mode) - except ModuleNotFoundError: - - def zero_petsc_vector(b: PETSc.Vec) -> None: - """Zero a PETSc vector, including ghosts""" - - if b.getType() == PETSc.Vec.Type.NEST: - for b_sub in b.getNestSubVecs(): - with b_sub.localForm() as b_local: - b_local.set(0.0) - else: - with b.localForm() as b_loc: - b_loc.set(0) - - def ghost_update( - x: PETSc.Vec, insert_mode: PETSc.InsertMode, scatter_mode: PETSc.ScatterMode - ) -> None: - """Ghost update a vector""" - if x.getType() == PETSc.Vec.Type.NEST: - for x_sub in x.getNestSubVecs(): - x_sub.ghostUpdate(addv=insert_mode, mode=scatter_mode) - else: - x.ghostUpdate(addv=insert_mode, mode=scatter_mode) -else: - - def zero_petsc_vector(_b: "PETSc.Vec") -> None: # type: ignore - """Zero a PETSc vector, including ghosts""" - raise RuntimeError("petsc4py is not available. Cannot zero vector.") - - def ghost_update( # type: ignore - _x: "PETSc.Vec", # type: ignore - _insert_mode: "PETSc.InsertMode", # type: ignore - _scatter_mode: "PETSc.ScatterMode", # type: ignore - ) -> None: # type: ignore - """Ghost update a vector""" - raise RuntimeError("petsc4py is not available. Cannot ghost update vector.") - - -if dolfinx.has_petsc4py: - from petsc4py import PETSc - import dolfinx.fem.petsc - - def apply_lifting_and_set_bc( - b: PETSc.Vec, - a: typing.Union[ - typing.Iterable[dolfinx.fem.Form], typing.Iterable[typing.Iterable[dolfinx.fem.Form]] - ], - bcs: typing.Union[ - typing.Iterable[dolfinx.fem.DirichletBC], - typing.Iterable[typing.Iterable[dolfinx.fem.DirichletBC]], - ], - x: typing.Optional[PETSc.Vec] = None, - alpha: float = 1.0, - ): - """Apply lifting to a vector and set boundary conditions. - - Convenience function to apply lifting and set boundary conditions for multiple matrix types. - This modifies the vector `b` such that - - .. math:: - b = - \\begin{pmatrix} - b_{free} \\\\ - b_{bc} - \\end{pmatrix} = - \\begin{pmatrix} - b_{free} - \\alpha \\sum_{i=0}^n a[i] (u_{bc}[i] - x[i])\\\\ - u_{bc}[0]\\\\ - \\vdots\\\\ - u_{bc}[n] - \\end{pmatrix}. - - where :math:`b_{free}` is the free part of the vector, :math:`b_{bc}` is the part that has - boundary conditions applied, :math:`u_{bc}[i]` is the value of the ith boundary condition. - - Args: - b: The vector to apply lifting to. - a: Sequence of forms to apply lifting from. If the system is blocked or nested, - this is a nested list of forms. - bcs: The boundary conditions to apply. If the form is blocked or nested, this is a list, - while if it is a single form, this is a nested list. - x: Vector to subtract from the boundary conditions. Usually used in a Newton iteration. - alpha: The scaling factor for the boundary conditions. - """ - if hasattr(dolfinx.fem.petsc, "create_vector_nest"): - raise NotImplementedError( - "This function is only implemented for later versions of DOLFINx." - ) - - try: - bcs0 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(a, 0), bcs) - bcs1 = dolfinx.fem.bcs_by_block(dolfinx.fem.extract_function_spaces(a, 1), bcs) - except AssertionError: - bcs0 = bcs - bcs1 = bcs - except ValueError: # Changed error type post https://github.com/FEniCS/dolfinx/pull/4312 - bcs0 = bcs - bcs1 = bcs - - dolfinx.fem.petsc.apply_lifting(b, a, bcs=bcs1, x0=x, alpha=alpha) - ghost_update(b, PETSc.InsertMode.ADD_VALUES, PETSc.ScatterMode.REVERSE) - try: - dolfinx.fem.petsc.set_bc(b, bcs0, x0=x, alpha=alpha) - except AttributeError: - for _bcs in bcs0: - dolfinx.fem.petsc.set_bc(b, _bcs, x0=x, alpha=alpha) - except TypeError: - # Catch the case when bcs0 shouldn't be a nested list (e.g., for a single form) - assert len(bcs0) == 1, "bcs0 should be a single DirichletBC or a list of DirichletBCs." - dolfinx.fem.petsc.set_bc(b, bcs0[0], x0=x, alpha=alpha) - - ghost_update(b, PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD) - -else: - def apply_lifting_and_set_bc(_b, _a, _bcs, _x, _alpha): # type: ignore - raise RuntimeError( - "petsc4py is not available. Cannot apply lifting and set boundary conditions." +def apply_lifting_and_set_bc( + b: PETSc.Vec, + a: typing.Union[ + typing.Iterable[dolfinx.fem.Form], typing.Iterable[typing.Iterable[dolfinx.fem.Form]] + ], + bcs: typing.Union[ + typing.Iterable[dolfinx.fem.DirichletBC], + typing.Iterable[typing.Iterable[dolfinx.fem.DirichletBC]], + ], + x: typing.Optional[PETSc.Vec] = None, + alpha: float = 1.0, +): + """Apply lifting to a vector and set boundary conditions. + + Convenience function to apply lifting and set boundary conditions for multiple matrix types. + This modifies the vector `b` such that + + .. math:: + b = + \\begin{pmatrix} + b_{free} \\\\ + b_{bc} + \\end{pmatrix} = + \\begin{pmatrix} + b_{free} - \\alpha \\sum_{i=0}^n a[i] (u_{bc}[i] - x[i])\\\\ + u_{bc}[0]\\\\ + \\vdots\\\\ + u_{bc}[n] + \\end{pmatrix}. + + where :math:`b_{free}` is the free part of the vector, :math:`b_{bc}` is the part that has + boundary conditions applied, :math:`u_{bc}[i]` is the value of the ith boundary condition. + + Args: + b: The vector to apply lifting to. + a: Sequence of forms to apply lifting from. If the system is blocked or nested, + this is a nested list of forms. + bcs: The boundary conditions to apply. If the form is blocked or nested, this is a list, + while if it is a single form, this is a nested list. + x: Vector to subtract from the boundary conditions. Usually used in a Newton iteration. + alpha: The scaling factor for the boundary conditions. + """ + if isinstance(a, (list, tuple)) and len(a) == 1 and isinstance(a[0], dolfinx.fem.Form): + # A single form, which this function takes as ``[form]`` with ``[bcs]`` + form_bcs = list(bcs) + fenicsx_compat.apply_lifting_and_set_bc( + b, a[0], form_bcs[0] if form_bcs else [], x=x, alpha=alpha ) + return + fenicsx_compat.apply_lifting_and_set_bc(b, a, bcs, x=x, alpha=alpha) diff --git a/src/scifem/point_source.py b/src/scifem/point_source.py index 69b3017..1a2f037 100644 --- a/src/scifem/point_source.py +++ b/src/scifem/point_source.py @@ -5,7 +5,6 @@ from __future__ import annotations from typing import TYPE_CHECKING -from packaging.version import Version from mpi4py import MPI if TYPE_CHECKING: @@ -19,11 +18,11 @@ _HAS_PETSC = False import dolfinx +import fenicsx_compat import numpy as np import numpy.typing as npt import ufl -from . import compat from .utils import unroll_dofmap __all__ = ["PointSource"] @@ -56,10 +55,6 @@ def __init__( If the point source is outside the mesh, a ``ValueError`` will be raised. """ self._function_space = V - if V.dofmap.bs > 1 and dolfinx.__version__ == "0.8.0": - raise NotImplementedError( - "Block function spaces are not supported in dolfinx 0.8.0. Please upgrade dolfinx" - ) if points.ndim == 1: points = points.reshape(1, -3) elif points.ndim == 2 and points.shape[1] != 3: @@ -92,28 +87,13 @@ def recompute_sources(self, tol: float | None = None): # Determine what process owns a point and what cells it lies within mesh = self._function_space.mesh tol = self._tol if tol is None else tol - if dolfinx.__version__ == "0.8.0": - src_ranks, _, self._points, self._cells = ( - dolfinx.cpp.geometry.determine_point_ownership( - mesh._cpp_object, self._input_points, tol - ) - ) - self._points = np.array(self._points).reshape(-1, 3) - elif Version(dolfinx.__version__) >= Version("0.9.0.0"): - try: - collision_data = dolfinx.cpp.geometry.determine_point_ownership( - mesh._cpp_object, self._input_points, tol - ) - except TypeError: - collision_data = dolfinx.geometry.determine_point_ownership( - mesh, self._input_points, padding=tol - ) - - self._points = collision_data.dest_points - self._cells = collision_data.dest_cells - src_ranks = collision_data.src_owner - else: - raise NotImplementedError(f"Unsupported version of dolfinx: {dolfinx.__version__}") + # Points are float32 or float64; the constrained TypeVar cannot take the union + collision_data = dolfinx.geometry.determine_point_ownership( # type: ignore[type-var] + mesh, self._input_points, padding=tol + ) + self._points = collision_data.dest_points + self._cells = collision_data.dest_cells + src_ranks = collision_data.src_owner if -1 in src_ranks: raise ValueError("Point source is outside the mesh.") @@ -122,11 +102,11 @@ def compute_cell_contributions(self): mesh = self._function_space.mesh # Pull owning points back to reference cell mesh_nodes = mesh.geometry.x - cmap = compat.cmap(mesh) + cmap = fenicsx_compat.cmap(mesh) ref_x = np.zeros((len(self._cells), mesh.topology.dim), dtype=mesh.geometry.x.dtype) for i, (point, cell) in enumerate(zip(self._points, self._cells)): - geom_dofs = compat.dofmap(mesh)[cell] + geom_dofs = fenicsx_compat.dofmap(mesh)[cell] ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs]) # Create expression evaluating a trial function (i.e. just the basis function) @@ -140,29 +120,15 @@ def compute_cell_contributions(self): expr = dolfinx.fem.Expression(u, ref_x, comm=MPI.COMM_SELF) all_values = expr.eval(mesh, self._cells) # Diagonalize values (num_cells, num_points, num_dofs, bs) -> (num_cells, num_dofs) - if Version(dolfinx.__version__) <= Version("0.9.0"): - basis_values = np.empty((len(self._cells), num_dofs), dtype=all_values.dtype) + basis_values = np.zeros((len(self._cells), num_dofs), dtype=dolfinx.default_scalar_type) + if bs == 1: + # Values have shape (num_cells, num_points, num_dofs) for i in range(len(self._cells)): - basis_values[i] = sum( - [ - all_values[i, i * num_dofs * bs : (i + 1) * num_dofs * bs][ - j * num_dofs : (j + 1) * num_dofs - ] - for j in range(bs) - ] - ) + basis_values[i] = all_values[i, i, :] else: - basis_values = np.zeros( - (len(self._cells), num_dofs), dtype=dolfinx.default_scalar_type - ) - if bs == 1: - # Values have shape (num_cells, num_points, num_dofs) - for i in range(len(self._cells)): - basis_values[i] = all_values[i, i, :] - else: - # Values have shape (num_cells, num_points, bs, num_dofs*bs) - for i in range(len(self._cells)): - basis_values[i] = sum(all_values[i, i, j, :] for j in range(bs)) + # Values have shape (num_cells, num_points, bs, num_dofs*bs) + for i in range(len(self._cells)): + basis_values[i] = sum(all_values[i, i, j, :] for j in range(bs)) else: basis_values = np.zeros((0, num_dofs), dtype=dolfinx.default_scalar_type) self._basis_values = basis_values diff --git a/src/scifem/solvers.py b/src/scifem/solvers.py index 2970279..a913fe2 100644 --- a/src/scifem/solvers.py +++ b/src/scifem/solvers.py @@ -9,14 +9,10 @@ else: from typing_extensions import deprecated - -from packaging.version import parse as _v - import numpy as np import ufl import dolfinx - -from . import compat +import fenicsx_compat __all__ = ["NewtonSolver"] @@ -25,20 +21,14 @@ from petsc4py import PETSc import dolfinx.fem.petsc - # assemble_vector_block(scale=...) is renamed assemble_vector_block(alpha=...) - # in 0.9 - _alpha_kw: str = "alpha" - if _v(dolfinx.__version__) < _v("0.9"): - _alpha_kw = "scale" - @deprecated("NewtonSolver is deprecated in favor of `dolfinx.fem.petsc.NonlinearProblem`.") class NewtonSolver: max_iterations: int _bcs: list[dolfinx.fem.DirichletBC] A: PETSc.Mat b: PETSc.Vec - _J: dolfinx.fem.Form | ufl.form.Form - _F: dolfinx.fem.Form | ufl.form.Form + _J: list[list[dolfinx.fem.Form | None]] + _F: list[dolfinx.fem.Form] _dx: PETSc.Vec _error_on_convergence: bool _pre_solve_callback: Callable[["NewtonSolver"], None] | None @@ -46,8 +36,8 @@ class NewtonSolver: def __init__( self, - F: list[dolfinx.fem.Form], - J: list[list[dolfinx.fem.Form]], + F: list[ufl.Form | dolfinx.fem.Form], + J: list[list[ufl.Form | dolfinx.fem.Form | None]], w: list[dolfinx.fem.Function], bcs: list[dolfinx.fem.DirichletBC] | None = None, max_iterations: int = 5, @@ -68,9 +58,11 @@ def __init__( error_on_nonconvergence: Raise an error if the linear solver does not converge or Newton solver doesn't converge """ - # Compile forms if not compiled. Will throw error it requires entity maps - self._F = dolfinx.fem.form(F) - self._J = dolfinx.fem.form(J) + # Compile forms if not compiled. Will throw error it requires entity maps. + # `dolfinx.fem.form` passes compiled forms through, but is typed for UFL + # forms only, and as if every residual block could compile to `None`. + self._F = dolfinx.fem.form(F) # type: ignore[arg-type, assignment] + self._J = dolfinx.fem.form(J) # type: ignore[arg-type, assignment] # Store solution and accessible/modifiable properties self._pre_solve_callback = None @@ -81,24 +73,18 @@ def __init__( self.w = w # Create PETSc objects for block assembly - try: - self.b = dolfinx.fem.petsc.create_vector_block(self._F) - self.A = dolfinx.fem.petsc.create_matrix_block(self._J) - self.dx = dolfinx.fem.petsc.create_vector_block(self._F) - self.x = dolfinx.fem.petsc.create_vector_block(self._F) - except AttributeError: - self.b = dolfinx.fem.petsc.create_vector( - dolfinx.fem.extract_function_spaces(self._F), kind="mpi" - ) - self.A = dolfinx.fem.petsc.create_matrix(self._J, kind="mpi") - self.dx = dolfinx.fem.petsc.create_vector( - dolfinx.fem.extract_function_spaces(self._F), kind="mpi" - ) - self.x = dolfinx.fem.petsc.create_vector( - dolfinx.fem.extract_function_spaces(self._F), kind="mpi" - ) + self.b = dolfinx.fem.petsc.create_vector( + dolfinx.fem.extract_function_spaces(self._F), kind="mpi" + ) + self.A = dolfinx.fem.petsc.create_matrix(self._J, kind="mpi") + self.dx = dolfinx.fem.petsc.create_vector( + dolfinx.fem.extract_function_spaces(self._F), kind="mpi" + ) + self.x = dolfinx.fem.petsc.create_vector( + dolfinx.fem.extract_function_spaces(self._F), kind="mpi" + ) - self._maps = [compat.form_map(form) for form in self._F] + self._maps = [fenicsx_compat.form_map(form) for form in self._F] # Set PETSc options opts = PETSc.Options() @@ -107,7 +93,7 @@ def __init__( opts[k] = v # Define KSP solver - self._solver = PETSc.KSP().create(self.b.getComm().tompi4py()) + self._solver = PETSc.KSP().create(self.b.getComm()) self._solver.setOperators(self.A) self._solver.setFromOptions() @@ -144,55 +130,26 @@ def solve(self, atol=1e-6, rtol=1e-8, beta=1.0) -> int: self._pre_solve_callback(self) # Pack constants and coefficients - try: - constants_L = dolfinx.fem.pack_constants(self._F) - coeffs_L = dolfinx.fem.pack_coefficients(self._F) - constants_a = [ - dolfinx.fem.pack_constants(form) - if form is not None - else np.array([], dtype=self.x.array.dtype) - for form in self._J - ] - coeffs_a = [ - {} if form is None else dolfinx.fem.pack_coefficients(form) - for form in self._J - ] - except AttributeError: - # NOTE: DOLFINx 0.9 compatibility - # Pack constants and coefficients - constants_L = [ - form and dolfinx.cpp.fem.pack_constants(form._cpp_object) - for form in self._F - ] - coeffs_L = [ - dolfinx.cpp.fem.pack_coefficients(form._cpp_object) for form in self._F - ] - constants_a = [ - [ - dolfinx.cpp.fem.pack_constants(form._cpp_object) - if form is not None - else np.array([], dtype=self.x.array.dtype) - for form in forms - ] - for forms in self._J - ] - coeffs_a = [ - [ - {} - if form is None - else dolfinx.cpp.fem.pack_coefficients(form._cpp_object) - for form in forms - ] - for forms in self._J - ] + constants_L = dolfinx.fem.pack_constants(self._F) + coeffs_L = dolfinx.fem.pack_coefficients(self._F) + constants_a = [ + dolfinx.fem.pack_constants(form) + if form is not None + else np.array([], dtype=self.x.array.dtype) + for form in self._J + ] + coeffs_a = [ + {} if form is None else dolfinx.fem.pack_coefficients(form) for form in self._J + ] # Scatter previous solution `w` to `self.x`, the blocked version used for lifting - if hasattr(dolfinx.common, "index_map"): - _imaps = [s.function_space.dofmap.index_map._cpp_object for s in self.w] - else: - _imaps = [s.function_space.dofmap.index_map for s in self.w] + _imaps = [ + fenicsx_compat.unwrap_index_map(s.function_space.dofmap.index_map) + for s in self.w + ] dolfinx.cpp.la.petsc.scatter_local_vectors( self.x, - [si.x.petsc_vec.array_r for si in self.w], + # The binding is typed for the one scalar type PETSc was built with. + [si.x.petsc_vec.array_r for si in self.w], # type: ignore[misc] [ ( _imaps[i], @@ -206,67 +163,47 @@ def solve(self, atol=1e-6, rtol=1e-8, beta=1.0) -> int: # Assemble F(u_{i-1}) - J(u_D - u_{i-1}) and set du|_bc= u_D - u_{i-1} with self.b.localForm() as b_local: b_local.set(0.0) - try: - dolfinx.fem.petsc.assemble_vector_block( - self.b, - self._F, - self._J, - bcs=self.bcs, - x0=self.x, - coeffs_a=coeffs_a, - constants_a=constants_a, - coeffs_L=coeffs_L, - constants_L=constants_L, - # dolfinx 0.8 compatibility - # this is called 'scale' in 0.8, 'alpha' in 0.9 - **{_alpha_kw: -1.0}, - ) - except AttributeError: - dolfinx.fem.petsc.assemble_vector( - self.b, self._F, coeffs=coeffs_L, constants=constants_L - ) - bcs1 = dolfinx.fem.bcs_by_block( - dolfinx.fem.extract_function_spaces(self._J, 1), self.bcs - ) - dolfinx.fem.petsc.apply_lifting( - self.b, - self._J, - bcs=bcs1, - x0=self.x, - coeffs=coeffs_a, - constants=constants_a, - alpha=-1.0, - ) - self.b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - bcs0 = dolfinx.fem.bcs_by_block( - dolfinx.fem.extract_function_spaces(self._F), self.bcs - ) - dolfinx.fem.petsc.set_bc(self.b, bcs0, x0=self.x, alpha=-1.0) - self.b.ghostUpdate(PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD) + dolfinx.fem.petsc.assemble_vector( + self.b, self._F, coeffs=coeffs_L, constants=constants_L + ) + bcs1 = fenicsx_compat.bcs_by_block( + dolfinx.fem.extract_function_spaces(self._J, 1), self.bcs + ) + dolfinx.fem.petsc.apply_lifting( + self.b, + self._J, + bcs=bcs1, + # A blocked `x0` is the single block vector, which dolfinx's own + # lifting reads, although typed as one vector per block. + x0=self.x, # type: ignore[arg-type] + coeffs=coeffs_a, + # `pack_constants` is typed to return `None` blocks, which + # `apply_lifting` and `assemble_matrix` are typed to refuse. + constants=constants_a, # type: ignore[arg-type] + alpha=-1.0, + ) + self.b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + bcs0 = fenicsx_compat.bcs_by_block( + dolfinx.fem.extract_function_spaces(self._F), self.bcs + ) + dolfinx.fem.petsc.set_bc(self.b, bcs0, x0=self.x, alpha=-1.0) + self.b.ghostUpdate(PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD) # Assemble Jacobian self.A.zeroEntries() - try: - dolfinx.fem.petsc.assemble_matrix_block( - self.A, - self._J, - bcs=self.bcs, - constants=constants_a, - coeffs=coeffs_a, - ) - except AttributeError: - dolfinx.fem.petsc.assemble_matrix( - self.A, - self._J, - self.bcs, - coeffs=coeffs_a, - constants=constants_a, - ) + dolfinx.fem.petsc.assemble_matrix( + self.A, + self._J, + self.bcs, + coeffs=coeffs_a, + constants=constants_a, # type: ignore[arg-type] + ) self.A.assemble() self._solver.solve(self.b, self.dx) if self._error_on_convergence: - if (status := self._solver.getConvergedReason()) <= 0: + # petsc4py types the reason as its enum class; it is an `int`. + if (status := self._solver.getConvergedReason()) <= 0: # type: ignore[operator] raise RuntimeError(f"Linear solver did not converge, got reason: {status}") # Update solution @@ -276,7 +213,9 @@ def solve(self, atol=1e-6, rtol=1e-8, beta=1.0) -> int: s.function_space.dofmap.index_map.size_local * s.function_space.dofmap.index_map_bs ) - s.x.petsc_vec.array_w[:num_sub_dofs] -= ( + # numpy's in-place operators reject the real-or-complex union petsc4py + # types its arrays with; at runtime both are the one PETSc scalar type. + s.x.petsc_vec.array_w[:num_sub_dofs] -= ( # type: ignore[misc] beta * self.dx.array_r[offset_start : offset_start + num_sub_dofs] ) s.x.petsc_vec.ghostUpdate( diff --git a/src/scifem/spaces.py b/src/scifem/spaces.py index ee9f2b9..70c114d 100644 --- a/src/scifem/spaces.py +++ b/src/scifem/spaces.py @@ -4,10 +4,9 @@ import dolfinx import basix import numpy as np +import fenicsx_compat from . import _scifem # type: ignore -from .compat import create_cpp_finite_element from collections.abc import Sequence -from packaging.version import Version def create_real_functionspace( @@ -25,17 +24,14 @@ def create_real_functionspace( For scalar elements value shape is ``()``. """ - if Version(dolfinx.__version__) >= Version("0.11.0"): + if fenicsx_compat.at_least("0.11.0"): warn( "`create_real_functionspace()` is deprecated and will be removed in a future version." + "Please use `basix.ufl.create_real_element()` instead.", category=DeprecationWarning, stacklevel=2, ) - el = basix.ufl.real_element( - mesh.basix_cell(), value_shape=value_shape, dtype=mesh.geometry.x.dtype - ) - return dolfinx.fem.functionspace(mesh, el) + return fenicsx_compat.real_functionspace(mesh, value_shape) dtype = mesh.geometry.x.dtype ufl_e = basix.ufl.element( @@ -55,7 +51,7 @@ def create_space_of_simple_functions( mesh: dolfinx.mesh.Mesh, cell_tag: dolfinx.mesh.MeshTags, tags: Sequence[int | np.int32 | np.int64], - value_shape: tuple[()] | tuple[int | tuple[int]] | None = None, + value_shape: tuple[int, ...] | None = None, ) -> dolfinx.fem.FunctionSpace: """Create a space of simple functions. @@ -78,11 +74,7 @@ def create_space_of_simple_functions( """ value_shape = () if value_shape is None else value_shape - if not isinstance(cell_tag.topology, dolfinx.cpp.mesh.Topology): - _cpp_topology = cell_tag.topology._cpp_object - else: - _cpp_topology = cell_tag.topology - if mesh.topology._cpp_object != _cpp_topology: + if mesh.topology._cpp_object != fenicsx_compat.meshtags_topology(cell_tag)._cpp_object: raise ValueError("Topology of cell tag is not the mesh topology") if cell_tag.dim != mesh.topology.dim: raise ValueError( @@ -94,7 +86,8 @@ def create_space_of_simple_functions( # Sanity check cell tags if cell_map.size_local + cell_map.num_ghosts != len(cell_tag.indices): raise ValueError("Every cell in the cell tag must have a value assigned to it.") - if not np.isin(cell_tag.values, tags).all(): + # numpy stubs reject a sequence mixing Python and numpy integers + if not np.isin(cell_tag.values, tags).all(): # type: ignore[arg-type] raise ValueError("All values in the cell tag must be in the input tags.") # Determine the owner of the material degrees of freedom. @@ -110,16 +103,7 @@ def create_space_of_simple_functions( # Create index map describing the distribution of dofs across processes ghosts = np.arange(num_ghosts, dtype=np.int64) owners = np.full(num_ghosts, ghost_owner_rank, dtype=np.int32) - if Version(dolfinx.__version__) > Version("0.9.0"): - imap_kwargs = {"tag": 321} - else: - imap_kwargs = {} - if hasattr(dolfinx.common, "index_map"): - dof_imap = dolfinx.common.index_map( - mesh.comm, num_dofs, ghosts=(ghosts, owners), **imap_kwargs - ) - else: - dof_imap = dolfinx.common.IndexMap(mesh.comm, num_dofs, ghosts, owners, **imap_kwargs) + dof_imap = fenicsx_compat.create_index_map(mesh.comm, num_dofs, ghosts, owners, tag=321) # Create element dof layout (1 dof per cell, based of the DG-0 element) value_size = int(np.prod(value_shape)) @@ -133,34 +117,32 @@ def create_space_of_simple_functions( # Create dofmap mapping local dofs to position in `tags` input. adj_flattened = np.zeros_like(cell_tag.values) for i, tag in enumerate(tags): - adj_flattened[cell_tag.find(tag)] = i + # numpy integer tags are accepted at runtime; the stub only admits int | float + adj_flattened[cell_tag.find(tag)] = i # type: ignore[arg-type] # Create dofmap + # Cell tag values are int32 at runtime; the stub admits any MeshTags value dtype dofmap_adj = dolfinx.cpp.graph.AdjacencyList_int32( - adj_flattened, np.arange(len(cell_tag.values) + 1, dtype=np.int32) + adj_flattened, # type: ignore[arg-type] + np.arange(len(cell_tag.values) + 1, dtype=np.int32), + ) + cpp_dofmap = dolfinx.cpp.fem.DofMap( + e_layout, fenicsx_compat.unwrap_index_map(dof_imap), index_map_bs, dofmap_adj, bs ) - if not isinstance(dof_imap, dolfinx.cpp.common.IndexMap): - _cpp_im = dof_imap._cpp_object - else: - _cpp_im = dof_imap - cpp_dofmap = dolfinx.cpp.fem.DofMap(e_layout, _cpp_im, index_map_bs, dofmap_adj, bs) # Create function space - try: - cpp_el = create_cpp_finite_element( - dolfinx.cpp.fem.FiniteElement_float64, - el.basix_element._e, - mesh.geometry.dim, - value_shape, # type: ignore - ) - cpp_space = dolfinx.cpp.fem.FunctionSpace_float64(mesh._cpp_object, cpp_el, cpp_dofmap) - except TypeError: - cpp_el = dolfinx.cpp.fem.FiniteElement_float64( - el.basix_element._e, block_size=bs, symmetric=False - ) - cpp_space = dolfinx.cpp.fem.FunctionSpace_float64( - mesh._cpp_object, cpp_el, cpp_dofmap, value_shape=value_shape - ) + cpp_el = fenicsx_compat.finite_element_ctor_kwargs( + dolfinx.cpp.fem.FiniteElement_float64, + el.basix_element._e, + gdim=mesh.geometry.dim, + value_shape=value_shape, + ) + # Only float64 meshes are supported here + cpp_space = dolfinx.cpp.fem.FunctionSpace_float64( + mesh._cpp_object, # type: ignore[arg-type] + cpp_el, + cpp_dofmap, + ) _el = basix.ufl.element("P", mesh.basix_cell(), 0, shape=value_shape, discontinuous=True) V = dolfinx.fem.FunctionSpace(mesh, _el, cpp_space) diff --git a/src/scifem/ufl_compat.py b/src/scifem/ufl_compat.py deleted file mode 100644 index ce11487..0000000 --- a/src/scifem/ufl_compat.py +++ /dev/null @@ -1,59 +0,0 @@ -"""Layer for small backward compatibility wrappers for UFL""" - -import ufl - -__all__ = ["apply_pullback_inverse"] - - -def _inverse_identity(expr: ufl.core.expr.Expr, domain: ufl.Mesh) -> ufl.core.expr.Expr: - return expr - - -def _inverse_contravariant_piola(expr: ufl.core.expr.Expr, domain: ufl.Mesh) -> ufl.core.expr.Expr: - """``u = J u_ref / det J`` inverted: ``u_ref = det J K u``.""" - return ufl.JacobianDeterminant(domain) * ufl.dot(ufl.JacobianInverse(domain), expr) - - -def _inverse_covariant_piola(expr: ufl.core.expr.Expr, domain: ufl.Mesh) -> ufl.core.expr.Expr: - """``u = K^T u_ref`` inverted: ``u_ref = J^T u``.""" - return ufl.dot(ufl.transpose(ufl.Jacobian(domain)), expr) - - -_INVERSE_PULLBACKS = { - "IdentityPullback": _inverse_identity, - "ContravariantPiola": _inverse_contravariant_piola, - "CovariantPiola": _inverse_covariant_piola, -} - - -def apply_pullback_inverse( - pullback: ufl.pullback.AbstractPullback, expr: ufl.core.expr.Expr, domain: ufl.Mesh -) -> ufl.core.expr.Expr: - """Map ``expr`` from the physical cell to the reference cell, across UFL versions. - - Uses ``pullback.apply_inverse``, added in https://github.com/FEniCS/ufl/pull/511, when the - installed UFL has it, and otherwise a closed form for the identity, contravariant Piola and - covariant Piola pullbacks. - - Args: - pullback: The element's pullback. - expr: A physical-cell expression. - domain: The domain whose Jacobian relates the two cells. - - Returns: - ``expr`` pulled back to the reference cell. - - Raises: - NotImplementedError: If this UFL has no ``apply_inverse`` and there is no closed form - for ``pullback`` here. - """ - if hasattr(pullback, "apply_inverse"): - return pullback.apply_inverse(expr, domain) - name = type(pullback).__name__ - try: - return _INVERSE_PULLBACKS[name](expr, domain) - except KeyError: - raise NotImplementedError( - f"This UFL has no {name}.apply_inverse (added in FEniCS/ufl#511), and there is no " - "closed form for it in scifem. Upgrade UFL." - ) from None diff --git a/src/scifem/xdmf.py b/src/scifem/xdmf.py index 1a5a7e9..c9a26db 100644 --- a/src/scifem/xdmf.py +++ b/src/scifem/xdmf.py @@ -11,6 +11,7 @@ import warnings from mpi4py import MPI +import fenicsx_compat import numpy as np import numpy.typing as npt import dolfinx @@ -243,10 +244,7 @@ def write_hdf5_adios( """ import adios2 - def resolve_adios_scope(adios2): - return adios2 if not hasattr(adios2, "bindings") else adios2.bindings - - adios2 = resolve_adios_scope(adios2) + adios2 = fenicsx_compat.resolve_adios_scope(adios2) if len(functions) == 0: return @@ -313,7 +311,8 @@ def create_function_space_data(V: dolfinx.fem.FunctionSpace) -> FunctionSpaceDat points = V.tabulate_dof_coordinates() bs = V.dofmap.index_map_bs index_map = V.dofmap.index_map - comm = V.mesh.comm + # dolfinx annotates Mesh.comm as MPI.Comm, but a mesh communicator is an intracommunicator + comm = typing.cast(MPI.Intracomm, V.mesh.comm) return FunctionSpaceData( points=points, @@ -481,10 +480,7 @@ def _init_adios(self) -> None: logger.debug("Initializing ADIOS2") import adios2 - def resolve_adios_scope(adios2): - return adios2 if not hasattr(adios2, "bindings") else adios2.bindings - - adios2 = resolve_adios_scope(adios2) + adios2 = fenicsx_compat.resolve_adios_scope(adios2) # Create ADIOS2 reader self._adios = adios2.ADIOS(self._data.comm) diff --git a/tests/periodic_mesh/test_gmsh_periodic.py b/tests/periodic_mesh/test_gmsh_periodic.py index 9c78dd8..b5c5771 100644 --- a/tests/periodic_mesh/test_gmsh_periodic.py +++ b/tests/periodic_mesh/test_gmsh_periodic.py @@ -27,7 +27,7 @@ PeriodicNodes, read_periodic_mesh_from_msh, ) -from scifem.periodic.gmsh import model_to_mesh +import fenicsx_compat @pytest.fixture @@ -308,7 +308,7 @@ def periodic_square(comm, res=1.0 / 5, directions=("x", "y"), low_is_replaced=Fa else: pairs = PeriodicNodes() - mesh = model_to_mesh(gmsh.model, comm, rank=0, gdim=2) + mesh = fenicsx_compat.model_to_mesh(gmsh.model, comm, 0, gdim=2).mesh if started_here: gmsh.finalize() return mesh, pairs @@ -538,11 +538,7 @@ def test_read_periodic_mesh_from_msh_round_trip(tmp_path): gmsh.finalize() comm.Barrier() - ghost_mode = dolfinx.mesh.GhostMode.shared_facet - partitioner = scifem.compat.create_partitioner(ghost_mode=ghost_mode, max_facet_to_cell_links=2) - periodic_mesh, _, _ = read_periodic_mesh_from_msh( - filename, comm, gdim=2, partitioner=partitioner - ) + periodic_mesh, _, _ = read_periodic_mesh_from_msh(filename, comm, gdim=2) num_vertices, volume, bad, jump = torus_invariants(periodic_mesh) assert (num_vertices, bad) == (33, 0) @@ -618,7 +614,7 @@ def periodic_box(comm, res=1.0 / 4, low_is_replaced=True, order=1): else: pairs = PeriodicNodes() - mesh = model_to_mesh(gmsh.model, comm, rank=0, gdim=3) + mesh = fenicsx_compat.model_to_mesh(gmsh.model, comm, 0, gdim=3).mesh if started_here: gmsh.finalize() return mesh, pairs diff --git a/tests/periodic_mesh/test_periodic.py b/tests/periodic_mesh/test_periodic.py index 86493a4..260d9c1 100644 --- a/tests/periodic_mesh/test_periodic.py +++ b/tests/periodic_mesh/test_periodic.py @@ -15,7 +15,6 @@ """ from mpi4py import MPI -import inspect import basix.ufl import numpy as np @@ -23,8 +22,8 @@ import ufl import dolfinx +import fenicsx_compat -from scifem.compat import create_partitioner from scifem.periodic.geometrical_search import match_vertices_geometric from scifem.periodic.mesh import create_periodic_mesh from scifem.periodic.mesh import ( @@ -646,10 +645,6 @@ def test_a_facet_with_three_distinct_cells_is_not_a_collapse(): cannot ask. `shared_facet` is what brings the third one over. """ comm = MPI.COMM_WORLD - # `max_facet_to_cell_links` defaults to 2, which is DOLFINx refusing exactly this mesh - # unless asked. Probed by signature rather than by version, as elsewhere in the suite. - if "max_facet_to_cell_links" not in inspect.signature(dolfinx.mesh.create_mesh).parameters: - pytest.skip("this DOLFINx cannot be asked for more than two cells per facet") # The input goes in on one rank and is distributed from there; handing every rank the # same cells would build the mesh `comm.size` times over. @@ -674,19 +669,17 @@ def test_a_facet_with_three_distinct_cells_is_not_a_collapse(): points = np.zeros((0, 3), dtype=np.float64) cells = np.zeros((0, 3), dtype=np.int64) - # Ghosting is a `create_mesh` keyword on 0.12 and a partitioner on 0.11, the same split - # `model_to_mesh` carries in test_gmsh_periodic.py. Both places have to be told that a - # facet may carry three cells, or the third sheet is not ghosted and the check sees two. - partitioner = create_partitioner(dolfinx.mesh.GhostMode.shared_facet, max_facet_to_cell_links=3) - ghost_mode = dolfinx.mesh.GhostMode.shared_facet - kwargs = {} - if "ghost_mode" in inspect.signature(dolfinx.mesh.create_mesh).parameters: - kwargs["ghost_mode"] = ghost_mode - kwargs["partitioner"] = partitioner - + # `max_facet_to_cell_links` defaults to 2, which is DOLFINx refusing exactly this mesh + # unless asked. It has to reach the ghosting as well as the mesh, or the third sheet is + # not ghosted and the check sees two. domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))) - mesh = dolfinx.mesh.create_mesh( - comm, cells, domain, points, max_facet_to_cell_links=3, **kwargs + mesh = fenicsx_compat.create_mesh( + comm, + cells, + domain, + points, + ghost_mode=dolfinx.mesh.GhostMode.shared_facet, + max_facet_to_cell_links=3, ) assert mesh.topology.index_map(2).size_global == 6 diff --git a/tests/periodic_mesh/test_transfer_function.py b/tests/periodic_mesh/test_transfer_function.py index 4fc9a2e..ff8bc88 100644 --- a/tests/periodic_mesh/test_transfer_function.py +++ b/tests/periodic_mesh/test_transfer_function.py @@ -24,12 +24,11 @@ import pytest import dolfinx -import inspect +import fenicsx_compat from scifem.periodic.gmsh import extract_gmsh_periodic_nodes from scifem.periodic.mesh import create_periodic_mesh_from_igi from scifem.periodic.transfer import transfer_function_to_parent_mesh from scifem.periodic.utils import PeriodicNodes -from scifem.compat import create_partitioner def _periodic_scalar(x): @@ -75,20 +74,13 @@ def meshes(): else: pairs = PeriodicNodes(num_nodes_global=0) # The rebuild needs a ghost layer across every interprocess facet, which - # `model_to_mesh` does not give by default. - gmodel_to_mesh = inspect.signature(dolfinx.io.gmsh.model_to_mesh) - kwargs = {} - ghost_mode = dolfinx.mesh.GhostMode.shared_facet - if "ghost_mode" in gmodel_to_mesh.parameters: - kwargs["ghost_mode"] = ghost_mode - partitioner = create_partitioner(ghost_mode) - mesh_data = dolfinx.io.gmsh.model_to_mesh( - gmsh.model, comm, 0, gdim=2, partitioner=partitioner, **kwargs - ) + # `dolfinx.io.gmsh.model_to_mesh` does not give by default. + mesh = fenicsx_compat.model_to_mesh( + gmsh.model, comm, 0, gdim=2, ghost_mode=dolfinx.mesh.GhostMode.shared_facet + ).mesh finally: gmsh.finalize() - mesh = getattr(mesh_data, "mesh", mesh_data) periodic_mesh, _, _ = create_periodic_mesh_from_igi( mesh, pairs.replaced, pairs.partner, pairs.num_nodes_global, root=0 ) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index 89c5e49..e735039 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -3,7 +3,8 @@ import numpy as np import numpy.typing as npt from scipy.optimize import minimize -from scifem import closest_point_projection, compat +import fenicsx_compat +from scifem import closest_point_projection from scifem.geometry import project_onto_simplex, _closest_point_projection import ufl import basix.ufl @@ -32,10 +33,12 @@ def scipy_project_point_to_element( """ # Extract scalar element of mesh - element = mesh.ufl_domain().ufl_coordinate_element().sub_elements[0] + ufl_domain = mesh.ufl_domain() + assert ufl_domain is not None + element = ufl_domain.ufl_coordinate_element().sub_elements[0] # Get the coordinates of the nodes for the specified cell - node_coords = mesh.geometry.x[compat.dofmap(mesh)[cells]][:, :, : mesh.geometry.dim] + node_coords = mesh.geometry.x[fenicsx_compat.dofmap(mesh)[cells]][:, :, : mesh.geometry.dim] target_points = target_points.reshape(-1, mesh.geometry.dim) # cmap = mesh.geometry.cmap diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 3f144e0..efea843 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -1,7 +1,6 @@ from mpi4py import MPI from packaging.version import Version import dolfinx -import scifem.compat import scifem.interpolation import scifem.mesh import pytest @@ -372,7 +371,7 @@ def _exterior_facet_submesh(mesh): num_facets = submesh.topology.index_map(fdim).size_local submesh_facets = np.arange(num_facets, dtype=np.int32) parent_facets = scifem.mesh.get_entity_map(entity_map)[:num_facets] - entities = scifem.compat.compute_integration_domains( + entities = dolfinx.fem.compute_integration_domains( dolfinx.fem.IntegralType.exterior_facet, mesh.topology, parent_facets ).reshape(-1, 2) return submesh, submesh_facets, entities, entity_map diff --git a/tests/test_point_source.py b/tests/test_point_source.py index cb00a7c..e5ff2dd 100644 --- a/tests/test_point_source.py +++ b/tests/test_point_source.py @@ -4,7 +4,8 @@ import basix import numpy as np import ufl -from scifem import PointSource, compat +import fenicsx_compat +from scifem import PointSource def test_midpoint(): @@ -91,8 +92,10 @@ def test_outside(): if cells.offsets[-1] > 0: cell = cells.links(0)[0] - geom_dofs = compat.dofmap(mesh)[cell] - ref_x = compat.cmap(mesh).pull_back(point.reshape(-1, 3), mesh.geometry.x[geom_dofs]) + geom_dofs = fenicsx_compat.dofmap(mesh)[cell] + ref_x = fenicsx_compat.cmap(mesh).pull_back( + point.reshape(-1, 3), mesh.geometry.x[geom_dofs] + ) ref_values = V.ufl_element().tabulate(0, ref_x).flatten() b_nonzero = np.flatnonzero(b.x.array) dofs = V.dofmap.cell_dofs(cell)