import inspect
from collections.abc import Callable
from mpi4py import MPI
import basix.ufl
import dolfinx
import dolfinx.fem
import dolfinx.graph
import numpy as np
import numpy.typing as npt
import ufl
[docs]
def cmap(mesh: dolfinx.mesh.Mesh) -> dolfinx.fem.CoordinateElement:
"""Get the coordinate element for a mesh's geometry.
dolfinx 0.10 has only the `geometry.cmap` property. From 0.11,
`geometry.cmaps` exists and `cmap` emits a deprecation warning, so
prefer `cmaps[0]` where it exists. Dispatches on `hasattr`.
"""
if hasattr(mesh.geometry, "cmaps"):
return mesh.geometry.cmaps[0]
return mesh.geometry.cmap
[docs]
def dofmap(mesh: dolfinx.mesh.Mesh) -> npt.NDArray[np.int32]:
"""Get the geometry dofmap for a mesh, across the same `dofmaps` split as `cmap`."""
if hasattr(mesh.geometry, "dofmaps"):
return mesh.geometry.dofmaps[0]
return mesh.geometry.dofmap
[docs]
def cell_permutation_info(topology: dolfinx.mesh.Topology) -> npt.NDArray[np.uint32]:
"""Compute and return the packed per-cell permutation info, ghost cells included.
dolfinx 0.10 and 0.11 compute it, together with the facet
permutations, in a single argument-free
`topology.create_entity_permutations()`. dolfinx PR #3904 (post-0.11)
split that into `create_cell_permutations()` and a per-dimension
`create_entity_permutations(dim)`. Dispatches on `hasattr`.
"""
if hasattr(topology, "create_cell_permutations"):
topology.create_cell_permutations()
else:
topology.create_entity_permutations() # type: ignore[call-arg]
return topology.get_cell_permutation_info()
[docs]
def facet_permutations(topology: dolfinx.mesh.Topology) -> npt.NDArray[np.uint8]:
"""Compute and return the permutation of every facet of every cell.
Shape `(num_cells, num_facets_per_cell)`, ghost cells included. Each
value encodes the facet's orientation as seen from the cell, as used
for FFCx's `quadrature_permutation`.
Across the same dolfinx PR #3904 split as `cell_permutation_info`:
0.10 and 0.11 expose these via `get_facet_permutations()` after an
argument-free `create_entity_permutations()`; from #3904 they are
`get_entity_permutations(tdim - 1)` after
`create_entity_permutations(tdim - 1)`. Both return a flat array.
"""
fdim = topology.dim - 1
if hasattr(topology, "create_cell_permutations"):
topology.create_entity_permutations(fdim)
perms = topology.get_entity_permutations(fdim)
else:
topology.create_entity_permutations() # type: ignore[call-arg]
perms = topology.get_facet_permutations() # type: ignore[attr-defined]
num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(topology.cell_type, fdim)
return np.asarray(perms).reshape(-1, num_facets_per_cell)
[docs]
def create_facet_orientations(topology: dolfinx.mesh.Topology) -> None:
"""Compute the facet orientations that contravariant Piola spaces on a manifold need.
dolfinx PR #4625 (post-0.11) refuses to create a Raviart-Thomas or
other contravariant Piola space on a manifold until
`topology.create_facet_orientations()` has been called, and that
method does not exist on 0.10 and 0.11, which create such spaces
without it. Calls it where it exists and does nothing otherwise.
Dispatches on `hasattr`.
Call it before creating the space, on every rank: it is collective,
and changing the orientations afterwards invalidates functions on
existing spaces. Where it exists it raises `ValueError` unless the
topological dimension is 2, as the native method does.
"""
if hasattr(topology, "create_facet_orientations"):
topology.create_facet_orientations()
[docs]
def create_cell_partitioner(
ghost_mode: dolfinx.mesh.GhostMode,
max_facet_to_cell_links: int = 2,
) -> Callable:
"""Build a cell partitioner with ghost mode baked in, for dolfinx 0.10 and 0.11.
Not needed (and not called) on dolfinx generations where
`create_cell_partitioner` was removed (post dolfinx PR #4403) —
`create_mesh` takes `ghost_mode` directly on those. Raises
`NotImplementedError` if called anyway, so a caller doesn't silently
get an un-ghosted mesh.
"""
if not hasattr(dolfinx.mesh, "create_cell_partitioner"):
raise NotImplementedError(
"dolfinx.mesh.create_cell_partitioner was removed in this dolfinx version "
"(see dolfinx PR #4403); pass ghost_mode directly to "
"fenicsx_compat.mesh.create_mesh() instead of calling create_cell_partitioner()."
)
return dolfinx.mesh.create_cell_partitioner(ghost_mode, max_facet_to_cell_links)
[docs]
def create_mesh(
comm: MPI.Comm,
cells: npt.NDArray[np.int64],
e,
x: npt.NDArray[np.floating],
*,
partitioner: Callable | None = None,
ghost_mode: dolfinx.mesh.GhostMode = dolfinx.mesh.GhostMode.shared_facet,
max_facet_to_cell_links: int = 2,
num_threads: int = 1,
) -> dolfinx.mesh.Mesh:
"""Create a mesh from topology and geometry arrays, with a given ghost mode.
Post dolfinx PR #4403 (post-0.11), `ghost_mode` and `num_threads` are
`create_mesh`'s own kwargs and `create_cell_partitioner` is gone. On
0.10 and 0.11 the ghost mode must be baked into the partitioner via
`create_cell_partitioner`; `num_threads` is ignored there. Dispatches
on `inspect.signature`.
If `partitioner` is given explicitly it is always used as-is (never
silently replaced); `ghost_mode` is still passed to `create_mesh` on
generations that take it.
The default `ghost_mode` is `shared_facet`, while dolfinx's own
default is `none` on every version.
"""
kwargs: dict = {}
if "ghost_mode" in inspect.signature(dolfinx.mesh.create_mesh).parameters:
kwargs = {"ghost_mode": ghost_mode, "num_threads": num_threads}
if partitioner is None and comm.size > 1:
partitioner = dolfinx.graph.partitioner()
elif partitioner is None and comm.size > 1:
partitioner = create_cell_partitioner(ghost_mode, max_facet_to_cell_links)
return dolfinx.mesh.create_mesh(
comm,
cells,
e,
x,
partitioner=partitioner,
max_facet_to_cell_links=max_facet_to_cell_links,
**kwargs,
)
[docs]
def reconstruct_mesh(mesh: dolfinx.mesh.Mesh, coordinate_element_degree: int) -> dolfinx.mesh.Mesh:
"""Copy a mesh, changing its coordinate element degree.
The topology is shared with the original mesh; the geometry is
reconstructed. Delegates to the native `dolfinx.fem.interpolate_geometry`
from dolfinx>=0.11; for 0.10, reconstructs the geometry manually (no
native equivalent exists at that version).
"""
if hasattr(dolfinx.fem, "interpolate_geometry"):
new_cmap = dolfinx.fem.coordinate_element(
mesh.topology.cell_type,
coordinate_element_degree,
dtype=mesh.geometry.x.dtype,
variant=cmap(mesh).variant,
)
return dolfinx.fem.interpolate_geometry(mesh, new_cmap)
ud = mesh.ufl_domain()
assert ud is not None
c_el = ud.ufl_coordinate_element()
family = c_el.family_name
lvar = c_el.lagrange_variant
ct = c_el.cell_type
new_c_el = basix.ufl.element(
family,
ct,
coordinate_element_degree,
shape=(mesh.geometry.dim,),
lagrange_variant=lvar,
dtype=mesh.geometry.x.dtype,
)
V_tmp = dolfinx.fem.functionspace(mesh, new_c_el)
gdim = mesh.geometry.dim
x = V_tmp.tabulate_dof_coordinates()[:, :gdim]
geom_imap = V_tmp.dofmap.index_map
geom_dofmap = V_tmp.dofmap.list
num_nodes_local = geom_imap.size_local + geom_imap.num_ghosts
original_input_indices = geom_imap.local_to_global(np.arange(num_nodes_local, dtype=np.int32))
coordinate_element = dolfinx.fem.coordinate_element(
mesh.topology.cell_type, coordinate_element_degree, lvar, dtype=mesh.geometry.x.dtype
)
# This package is deliberately generic over the mesh's float32/float64
# scalar type, while dolfinx's Geometry/Mesh/CoordinateElement cpp
# bindings are separate per-scalar-type overloads (Geometry_float32 vs.
# Geometry_float64, etc.) -- mypy cannot narrow which concrete overload
# applies here, since both are valid depending on the input mesh's dtype.
geom = dolfinx.mesh.Geometry(
type(mesh.geometry._cpp_object)(
geom_imap, # type: ignore[arg-type]
geom_dofmap,
coordinate_element._cpp_object, # type: ignore[arg-type]
x,
original_input_indices,
)
)
new_top = mesh.topology
cpp_mesh = type(mesh._cpp_object)(
mesh.comm,
new_top._cpp_object,
geom._cpp_object, # type: ignore[arg-type]
)
return dolfinx.mesh.Mesh(cpp_mesh, ufl.Mesh(new_c_el))