# Copyright (C) 2017-2024 Chris N. Richardson, Garth N. Wells and
# Jørgen S. Dokken
#
# This file is part of DOLFINx (https://www.fenicsproject.org)
#
# SPDX-License-Identifier: LGPL-3.0-or-later
"""Creation, refining and marking of meshes."""
from __future__ import annotations
import typing
import warnings
from collections.abc import Callable, Sequence
from mpi4py import MPI as _MPI
import numpy as np
import numpy.typing as npt
import basix
import basix.ufl
import ufl
from dolfinx import cpp as _cpp
from dolfinx import default_real_type
from dolfinx.common import IndexMap as _IndexMap
from dolfinx.cpp.mesh import (
CellType,
DiagonalType,
GhostMode,
build_dual_graph,
cell_dim,
to_string,
to_type,
)
from dolfinx.cpp.refinement import (
IdentityPartitionerPlaceholder,
RefinementOption,
mark_maximum,
)
from dolfinx.cpp.refinement import (
uniform_refine as _uniform_refine,
)
from dolfinx.fem import CoordinateElement as _CoordinateElement
from dolfinx.fem.element import _coordinate_element_from_basix
from dolfinx.graph import AdjacencyList
from dolfinx.typing import Real
__all__ = [
"CellType",
"DiagonalType",
"EntityMap",
"Geometry",
"GhostMode",
"IdentityPartitionerPlaceholder",
"Mesh",
"MeshTags",
"RefinementOption",
"Topology",
"build_dual_graph",
"cell_dim",
"compute_cell_centroids",
"compute_incident_entities",
"compute_midpoints",
"create_box",
"create_geometric_cell_partitioner",
"create_geometry",
"create_hybrid_cell_partitioner",
"create_interval",
"create_mesh",
"create_point_mesh",
"create_rectangle",
"create_submesh",
"create_unit_cube",
"create_unit_interval",
"create_unit_square",
"entities_to_geometry",
"exterior_facet_indices",
"locate_entities",
"locate_entities_boundary",
"mark_maximum",
"meshtags",
"meshtags_from_entities",
"refine",
"to_string",
"to_type",
"transfer_meshtag",
"transfer_meshtags_to_submesh",
"uniform_refine",
]
PartitioningFunc = typing.Callable[
[
_MPI.Comm,
int,
_cpp.graph.AdjacencyList_int64,
npt.NDArray[np.int32] | None,
npt.NDArray[np.int32] | None,
bool,
],
_cpp.graph.AdjacencyList_int32,
]
# float64 to match the C++ span<double> interface. Unlike
# PartitioningFunc/HybridPartitioningFunc, the return is a plain
# destination rank per node, not an AdjacencyList: a geometric
# partitioner never ghosts.
GeometricPartitioningFunc = typing.Callable[
[_MPI.Comm, int, npt.NDArray[np.float64], npt.NDArray[np.int32] | None],
npt.NDArray[np.int32],
]
HybridPartitioningFunc = typing.Callable[
[
_MPI.Comm,
int,
_cpp.graph.AdjacencyList_int64,
npt.NDArray[np.float64],
npt.NDArray[np.int32] | None,
npt.NDArray[np.int32] | None,
bool,
],
_cpp.graph.AdjacencyList_int32,
]
[docs]
def create_geometric_cell_partitioner(
part: GeometricPartitioningFunc,
) -> _cpp.graph.GeometricPartitioner:
"""Create a function to partition a mesh using cell positions.
``part`` is called by :func:`create_mesh` with each cell's centroid
(mean of its vertex positions), letting position-based partitioners
be used where :func:`create_mesh` otherwise expects graph edges.
The returned partitioner has no access to the dual graph, so it
never ghosts: the mesh always has ``GhostMode.none`` regardless of
what was requested. Use :func:`create_hybrid_cell_partitioner` if
the partitioner also needs the dual graph, e.g. to support ghosting.
``part`` also receives optional node weights, matching
:data:`PartitioningFunc`; :func:`create_mesh` supplies its
``cell_weights`` as the node weights. Not every geometric
partitioner honours these, e.g. ``partitioner_parmetis_geom`` raises
if given anything other than ``None``.
Args:
part: A custom geometric graph partitioning function. The
built-in geometric partitioners (``dolfinx.cpp.graph``'s
``partition_morton``, ``partition_hilbert`` and
``partitioner_parmetis_geom``) already return a
:class:`~dolfinx.cpp.graph.GeometricPartitioner` and do not
need to be wrapped.
Return:
Partitioning function, for use as :func:`create_mesh`'s
``partitioner`` argument.
"""
return _cpp.graph.create_geometric_cell_partitioner(part)
[docs]
def create_hybrid_cell_partitioner(
part: HybridPartitioningFunc,
) -> _cpp.graph.HybridPartitioner:
"""Create a function to partition a mesh using a hybrid partitioner.
Unlike :func:`create_geometric_cell_partitioner`, the dual graph is
always computed and passed to ``part`` along with cell positions,
regardless of ghost mode. Use this when the partitioner needs graph
edges for the partitioning decision itself, not only to determine
ghost cells.
The dual graph (its ``max_facet_to_cell_links`` and ``num_threads``)
is controlled at call time by :func:`create_mesh`'s matching
arguments, rather than fixed when the partitioner is created.
``part`` also receives optional node and edge weights, matching
:data:`PartitioningFunc`; :func:`create_mesh` supplies its
``cell_weights`` as the node weights and always passes ``None`` for
edge weights, since it has no way to compute them.
Args:
part: A custom hybrid graph partitioning function. ParMETIS's
built-in hybrid partitioner (``dolfinx.cpp.graph``'s
``partitioner_parmetis_hybrid``, i.e. ParMETIS `GeomKway`)
already returns a :class:`~dolfinx.cpp.graph.HybridPartitioner`
and does not need to be wrapped.
Return:
Partitioning function, for use as :func:`create_mesh`'s
``partitioner`` argument.
"""
return _cpp.graph.create_hybrid_cell_partitioner(part)
[docs]
def compute_cell_centroids(
comm: _MPI.Comm,
cell_types: Sequence[CellType],
cells: Sequence[npt.NDArray[np.int64]],
commg: _MPI.Comm,
x: npt.NDArray[np.floating],
) -> npt.NDArray[np.float64]:
"""Compute the centroid of each cell from its vertex positions.
:func:`create_mesh` performs this computation internally before
calling a partitioner created by
:func:`create_geometric_cell_partitioner` or
:func:`create_hybrid_cell_partitioner`. Use it directly to call such
a partitioner without going through :func:`create_mesh`.
Args:
comm: MPI communicator that ``cells`` is distributed over.
cell_types: Cell types in the mesh.
cells: Cells of each cell type, using global vertex indices, as
accepted by :func:`create_mesh`.
commg: MPI communicator that ``x`` is distributed over.
x: Geometry ('node') coordinates, distributed over ``commg``.
Returns:
Cell centroids, with shape ``(num_cells, gdim)``, with the cells
of each cell type concatenated in the order they appear in
``cells``.
"""
return _cpp.mesh.compute_cell_centroids(
comm, cell_types, cells, commg, np.ascontiguousarray(x, dtype=np.float64)
)
[docs]
class Topology:
"""Topology for a :class:`dolfinx.mesh.Mesh`."""
_cpp_object: _cpp.mesh.Topology
def __init__(self, topology: _cpp.mesh.Topology):
"""Initialize a topology from a C++ topology.
Args:
topology: The C++ topology object
Note:
Topology objects should usually be constructed with the
:func:`dolfinx.cpp.mesh.create_topology` and not this class
initializer.
"""
self._cpp_object = topology
[docs]
def cell_name(self) -> str:
"""String representation of the cell-type of the topology."""
return to_string(self._cpp_object.cell_type)
[docs]
def connectivity(self, d0: int, d1: int) -> _cpp.graph.AdjacencyList_int32:
"""Return connectivity.
Connectivity from entities of dimension ``d0`` to entities of
dimension ``d1``.
Args:
d0: Dimension of entity one is mapping from.
d1: Dimension of entity one is mapping to.
"""
if (conn := self._cpp_object.connectivity(d0, d1)) is not None:
return conn
else:
raise RuntimeError(
f"Connectivity between dimension {d0} and {d1} has not been computed.",
f"Please call `dolfinx.mesh.Topology.create_connectivity({d0}, {d1})` first.",
)
@property
def comm(self) -> _MPI.Comm:
"""Return the MPI communicator associated with the topology."""
return self._cpp_object.comm
[docs]
def create_connectivity(self, d0: int, d1: int) -> None:
"""Build entity connectivity ``d0 -> d1``.
Args:
d0: Dimension of entities connectivity is from.
d1: Dimension of entities connectivity is to.
"""
self._cpp_object.create_connectivity(d0, d1)
[docs]
def create_entities(self, dim: int, num_threads: int = 1) -> bool:
"""Create entities of given topological dimension.
Args:
dim: Topological dimension of entities to create.
num_threads: Number of CPU threads to use. Must be >= 1.
Returns:
``True` is entities are created, ``False`` is if entities
already existed.
"""
return self._cpp_object.create_entities(dim, num_threads)
[docs]
def create_entity_permutations(self, num_threads: int = 1) -> None:
"""Compute entity permutations and reflections.
Args:
num_threads: Number of CPU threads to use. Must be >= 1.
"""
self._cpp_object.create_entity_permutations(num_threads)
@property
def dim(self) -> int:
"""Return the topological dimension of the mesh."""
return self._cpp_object.dim
@property
def entity_types(self) -> list[list[CellType]]:
"""Entity types in the topology for all topological dimensions."""
return self._cpp_object.entity_types
[docs]
def get_cell_permutation_info(self) -> npt.NDArray[np.uint32]:
"""Get permutation information.
The returned data is used for packing coefficients and
assembling of tensors. The bits of each integer encodes the
number of reflections and permutations for each sub-entity of
the cell to be able to map it to the reference element.
"""
return self._cpp_object.get_cell_permutation_info()
[docs]
def get_facet_permutations(self) -> npt.NDArray[np.uint8]:
"""Get the permutation integer to apply to facets.
The bits of each integer describes the number of reflections and
rotations that has to be applied to a facet to map between a
facet in the mesh (relative to a cell) and the corresponding
facet on the reference element. The data has the shape
``(num_cells, num_facets_per_cell)``, flattened row-wise. The
number of cells include potential ghost cells.
Note:
The data can be unpacked with ``numpy.unpack_bits``.
"""
return self._cpp_object.get_facet_permutations()
[docs]
def index_map(self, dim: int) -> _cpp.common.IndexMap:
"""Index map for the parallel distribution of the mesh entities.
Args:
dim: Topological dimension.
Returns:
Index map for the entities of dimension ``dim``.
"""
return self._cpp_object.index_map(dim)
[docs]
def index_maps(self, dim: int) -> list[_cpp.common.IndexMap]:
"""Index maps for parallel distribution of the mesh entities.
Args:
dim: Topological dimension of the entities.
Returns:
Index maps for the entities of dimension ``dim``. May be
empty if not yet computed.
"""
return self._cpp_object.index_maps(dim)
[docs]
def interprocess_facets(self) -> npt.NDArray[np.int32]:
"""List of inter-process facets.
Empty if facet topology has not been computed.
"""
return self._cpp_object.interprocess_facets()
@property
def original_cell_index(self) -> npt.NDArray[np.int64]:
"""Original cell index."""
return self._cpp_object.original_cell_index
@property
def cell_type(self) -> CellType:
"""Get the cell type of the topology."""
return self._cpp_object.cell_type
[docs]
class Geometry(typing.Generic[Real]):
"""The geometry of a :class:`dolfinx.mesh.Mesh`."""
_cpp_object: _cpp.mesh.Geometry_float32 | _cpp.mesh.Geometry_float64
def __init__(self, geometry: _cpp.mesh.Geometry_float32 | _cpp.mesh.Geometry_float64):
"""Initialize a geometry from a C++ geometry.
Args:
geometry: The C++ geometry object.
Note:
Geometry objects should usually be constructed with the
:func:`create_geometry` and not using this class
initializer. This class is combined with different base
classes that depend on the scalar type used in the Geometry.
"""
self._cpp_object = geometry
@property
def cmaps(self) -> list[_CoordinateElement]:
"""The coordinate maps."""
return [_CoordinateElement(cm) for cm in self._cpp_object.cmaps]
@property
def cmap(self) -> _CoordinateElement:
"""The coordinate map (deprecated, use ``cmaps[0]``)."""
warnings.warn(
"cmap is deprecated. Use cmaps[0] instead.",
DeprecationWarning,
stacklevel=2,
)
return self.cmaps[0]
@property
def dim(self) -> int:
"""Dimension of the Euclidean coordinate system."""
return self._cpp_object.dim
@property
def dofmaps(self) -> list[npt.NDArray[np.int32]]:
"""The geometry dofmaps, one per cell type."""
return self._cpp_object.dofmaps
@property
def dofmap(self) -> npt.NDArray[np.int32]:
"""Dofmap for the geometry (deprecated, use ``dofmaps[0]``).
Shape is ``(num_cells, dofs_per_cell)``.
"""
warnings.warn(
"dofmap is deprecated. Use dofmaps[0] instead.",
DeprecationWarning,
stacklevel=2,
)
return self.dofmaps[0]
[docs]
def index_map(self) -> _IndexMap:
"""Index map for the geometry points (nodes) distribution."""
return self._cpp_object.index_map()
@property
def input_global_indices(self) -> npt.NDArray[np.int64]:
"""Global input indices of the geometry nodes."""
return self._cpp_object.input_global_indices
@property
def x(self) -> npt.NDArray[Real]:
"""Geometry coordinate points.
Shape is ``shape=(num_points, 3)``.
"""
return self._cpp_object.x # type: ignore[return-value]
[docs]
class Mesh(typing.Generic[Real]):
"""A mesh."""
_mesh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
_topology: Topology
_geometry: Geometry[Real]
_ufl_domain: ufl.Mesh | None
def __init__(
self,
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64,
domain: ufl.Mesh | None,
):
"""Initialize mesh from a C++ mesh.
Args:
msh: A C++ mesh object.
domain: A UFL domain.
Note:
Mesh objects should usually be constructed using
:func:`create_mesh` and not using this class initializer.
This class is combined with different base classes that
depend on the scalar type used in the Mesh.
"""
self._cpp_object = msh
self._topology = Topology(self._cpp_object.topology)
self._geometry = Geometry(self._cpp_object.geometry)
self._ufl_domain = domain
if self._ufl_domain is not None:
self._ufl_domain._ufl_cargo = self._cpp_object
@property
def comm(self) -> _MPI.Comm:
"""MPI communicator associated with the mesh."""
return self._cpp_object.comm
@property
def name(self) -> str:
"""Name of the mesh."""
return self._cpp_object.name
@name.setter
def name(self, value: str) -> None:
self._cpp_object.name = value
[docs]
def ufl_cell(self) -> ufl.Cell:
"""Return the UFL cell type.
Note:
This method is required for UFL compatibility.
"""
return ufl.Cell(self.topology.cell_name())
[docs]
def ufl_domain(self) -> ufl.Mesh | None:
"""Return the ufl domain corresponding to the mesh.
Returns:
The UFL domain. Is ``None`` if the domain has not been set.
Note:
This method is required for UFL compatibility.
"""
return self._ufl_domain
[docs]
def basix_cell(self) -> basix.CellType:
"""Return the Basix cell type."""
return typing.cast(basix.CellType, getattr(basix.CellType, self.topology.cell_name()))
[docs]
def h(self, dim: int, entities: npt.NDArray[np.int32]) -> npt.NDArray[Real]:
"""Geometric size measure of cell entities.
Args:
dim: Topological dimension of the entities to compute the
size measure of.
entities: Indices of entities of dimension ``dim`` to
compute size measure of.
Returns:
Size measure for each requested entity.
"""
return _cpp.mesh.h(self._cpp_object, dim, entities) # type: ignore[return-value]
@property
def topology(self) -> Topology:
"""Mesh topology."""
return self._topology
@property
def geometry(self) -> Geometry[Real]:
"""Mesh geometry."""
return self._geometry
[docs]
class EntityMap:
"""Bidirectional map between entities in two topologies."""
def __init__(self, entity_map: _cpp.mesh.EntityMap):
"""Initialise an entity map from a C++ `EntityMap` object.
Note:
`EntityMap` objects should not usually be created using this
initializer directly.
Args:
entity_map: A C++ `EntityMap` object.
"""
self._cpp_object = entity_map
self._topology = Topology(self._cpp_object.topology)
self._sub_topology = Topology(self._cpp_object.sub_topology)
[docs]
def sub_topology_to_topology(
self, entities: npt.NDArray[np.int32], inverse: bool
) -> npt.NDArray[np.int32]:
"""Map entities between the sub-topology and the parent topology.
If `inverse` is False, this function maps a list of
`self.dim()`-dimensional entities from `self.sub_topology()` to
the corresponding entities in `self.topology()`. If `inverse` is
True, it performs the inverse mapping from `self.topology()` to
`self.sub_topology()`. Entities that do not exist in the
sub-topology are marked as -1.
Note:
If `inverse` is `True`, this function recomputes the inverse
map on every call (it is not cached), which may be expensive
if called repeatedly.
Args:
entities:
A list of entity indices in the source topology.
inverse:
If False, maps from `self.sub_topology()` to
`self.topology()`. If True, maps from `self.topology()`
to `self.sub_topology()`.
Returns:
A list of mapped entity indices. Entities that don't exist
in the target topology are marked as -1.
"""
return self._cpp_object.sub_topology_to_topology(entities, inverse)
@property
def dim(self) -> int:
"""Topological dimension of the entities."""
return self._cpp_object.dim
@property
def topology(self) -> Topology:
"""The topology."""
return self._topology
@property
def sub_topology(self) -> Topology:
"""The sub-topology."""
return self._sub_topology
def entity_map(
topology: Topology,
sub_topology: Topology,
dim: int,
sub_topology_to_topology: npt.NDArray[np.int32],
) -> EntityMap:
"""Create a bidirectional map between (sub) topologies.
The map relates entities of dimension `dim` in `topology` and
`sub_topology`.
Args:
topology: A topology.
sub_topology: Topology of another mesh. This must be a
"sub-topology" of `topology` i.e. every entity in
`sub_topology` must also exist in `topology`.
dim: The dimension of the entities
sub_topology_to_topology: A list of entities in `topology` where
`sub_topology_to_topology[i]` is the index in `topology`
corresponding to entity `i` in `sub_topology`.
"""
return EntityMap(
_cpp.mesh.EntityMap(
topology._cpp_object, sub_topology._cpp_object, dim, sub_topology_to_topology
)
)
[docs]
def compute_incident_entities(
topology: Topology, entities: npt.NDArray[np.int32], d0: int, d1: int
) -> npt.NDArray[np.int32]:
"""Compute incident entities.
Computes entities of dimension ``d1`` connected to ``entities`` of
dimension ``d0``.
Args:
topology: The topology.
entities: List of entities of dimension ``d0``.
d0: Topological dimension.
d1: Topological dimension to map to.
Returns:
Incident entity indices.
"""
return _cpp.mesh.compute_incident_entities(topology._cpp_object, entities, d0, d1)
[docs]
def compute_midpoints(
msh: Mesh[Real], dim: int, entities: npt.NDArray[np.int32]
) -> npt.NDArray[Real]:
"""Compute the midpoints of a set of mesh entities.
Args:
msh: The mesh.
dim: Topological dimension of the mesh entities to consider.
entities: Indices of entities in ``mesh`` to consider.
Returns:
Midpoints of the entities, shape ``(num_entities, 3)``.
"""
return _cpp.mesh.compute_midpoints(msh._cpp_object, dim, entities) # type: ignore[return-value]
[docs]
def locate_entities(msh: Mesh, dim: int, marker: Callable) -> np.ndarray:
"""Compute mesh entities satisfying a geometric marking function.
Args:
msh: Mesh to locate entities on.
dim: Topological dimension of the mesh entities to consider.
marker: A function that takes an array of points ``x`` with
shape ``(gdim, num_points)`` and returns an array of
booleans of length ``num_points``, evaluating to `True` for
entities to be located.
Returns:
Indices (local to the process) of marked mesh entities.
"""
return _cpp.mesh.locate_entities(msh._cpp_object, dim, marker)
[docs]
def locate_entities_boundary(msh: Mesh, dim: int, marker: Callable) -> np.ndarray:
"""Mesh entities connected to a boundary facet a geometric condition.
Compute indices of all mesh entities that are attached to an owned
boundary facet and evaluate to true for the provided geometric
marking function. An entity is considered marked if the marker
function evaluates to true for all of its vertices.
For vertices and edges, in parallel this function will not
necessarily mark all entities that are on the exterior boundary. For
example, it is possible for a process to have a vertex that lies on
the boundary without any of the attached facets being a boundary
facet. When used to find degrees-of-freedom, e.g. using
:func:`dolfinx.fem.locate_dofs_topological`, the function that uses
the data returned by this function must typically perform some
parallel communication.
Args:
msh: Mesh to locate boundary entities on.
dim: Topological dimension of the mesh entities to consider
marker: Function that takes an array of points ``x`` with shape
``(gdim, num_points)`` and returns an array of booleans of
length ``num_points``, evaluating to ``True`` for entities
to be located.
Returns:
Indices (local to the process) of marked mesh entities.
"""
return _cpp.mesh.locate_entities_boundary(msh._cpp_object, dim, marker)
[docs]
def transfer_meshtag(
meshtag: MeshTags,
msh1: Mesh,
parent_cell: npt.NDArray[np.int32],
parent_facet: npt.NDArray[np.int8] | None = None,
) -> MeshTags:
"""Create cell mesh tags on a refined mesh from parent mesh tags.
Args:
meshtag: Mesh tags on the coarse, parent mesh.
msh1: The refined mesh.
parent_cell: Index of the parent cell for each cell in the
refined mesh.
parent_facet: Index of the local parent facet for each cell
in the refined mesh. Only required for transfer tags on
facets.
Returns:
Mesh tags on the refined mesh.
"""
if not isinstance(meshtag._cpp_object, _cpp.mesh.MeshTags_int32):
raise TypeError("transfer_meshtag only supports MeshTags with int32 values.")
if meshtag.dim == meshtag.topology.dim:
mt = _cpp.refinement.transfer_cell_meshtag(
meshtag._cpp_object, msh1.topology._cpp_object, parent_cell
)
return MeshTags(mt)
elif meshtag.dim == meshtag.topology.dim - 1:
if parent_facet is None:
raise ValueError("parent_facet is required for transferring facet tags.")
mt = _cpp.refinement.transfer_facet_meshtag(
meshtag._cpp_object, msh1.topology._cpp_object, parent_cell, parent_facet
)
return MeshTags(mt)
else:
raise RuntimeError("MeshTag transfer is supported on cells or facets.")
[docs]
def refine(
msh: Mesh,
edges: np.ndarray | None = None,
partitioner: Callable | IdentityPartitionerPlaceholder | None = (
IdentityPartitionerPlaceholder()
),
option: RefinementOption = RefinementOption.parent_cell,
ghost_mode: GhostMode = GhostMode.none,
) -> tuple[Mesh, npt.NDArray[np.int32] | None, npt.NDArray[np.int8] | None]:
"""Refine a mesh.
Note:
Using the `None` partitioner for the refined mesh, the refined mesh
will **not** include ghosts cells (cells connected by facet to an
owned cells) even if the parent mesh is ghosted. The refined cells
will be on the same process as the parent cell.
Args:
msh: Mesh from which to create the refined mesh.
edges: Indices of edges to split during refinement. If ``None``,
mesh refinement is uniform.
partitioner: Partitioner to distribute the refined mesh. If a
:py:class:`IdentityPartitionerPlaceholder` is passed (default)
no redistribution is performed, i.e. refined cells remain on
the same process as the parent cell, but the ghost layer is
updated. If a custom partitioner is passed, it will be used for
distributing the refined mesh. If ``None`` is passed no
redistribution will happen.
option: Controls whether parent cells and/or parent facets are
computed.
ghost_mode: Ghost mode used in the mesh partitioning, passed to
``partitioner`` at call time if it is a callable. Has no
effect if ``partitioner`` is an
:py:class:`IdentityPartitionerPlaceholder` (see the note
above).
Returns:
Refined mesh, (optional) parent cells, (optional) parent facets
"""
mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
msh._cpp_object, edges, partitioner, option, ghost_mode
)
if msh._ufl_domain is None:
raise ValueError("Cannot refine a mesh without a UFL domain.")
# Create new ufl domain as it will carry a reference to the C++ mesh
# in the ufl_cargo
ufl_domain = ufl.Mesh(msh._ufl_domain.ufl_coordinate_element())
return Mesh(mesh1, ufl_domain), parent_cell, parent_facet
_MeshPartitioner = Callable | tuple[Callable, npt.NDArray[np.int32] | None] | None
def _get_mesh_partitioner(
comm: _MPI.Comm, partitioner: _MeshPartitioner
) -> tuple[Callable | None, npt.NDArray[np.int32] | None]:
"""Get a cell partitioner and optional cell weights."""
partitioner_fn: Callable | None
if partitioner is None:
partitioner_fn = _cpp.graph.partitioner() if comm.size > 1 else None
return partitioner_fn, None
if isinstance(partitioner, tuple):
if len(partitioner) != 2:
raise TypeError("partitioner tuple must have two entries.")
partitioner_fn, cell_weights = partitioner
if partitioner_fn is None:
raise TypeError("The partitioner in the tuple must not be None.")
return partitioner_fn, cell_weights
return partitioner, None
def _create_mesh_coordinate_element(
e: ufl.Mesh | basix.finite_element.FiniteElement | basix.ufl._BasixElement | _CoordinateElement,
gdim: int,
) -> tuple[_CoordinateElement, ufl.Mesh | None]:
"""Create a coordinate element and UFL domain from a mesh element."""
if isinstance(e, ufl.domain.Mesh):
e_ufl = e.ufl_coordinate_element()
return _coordinate_element_from_basix(e_ufl.basix_element), e
if isinstance(e, basix.finite_element.FiniteElement):
cmap = _coordinate_element_from_basix(e)
e_ufl = basix.ufl.blocked_element(basix.ufl._BasixElement(e), shape=(gdim,))
domain = ufl.Mesh(e_ufl)
assert domain.geometric_dimension == gdim
return cmap, domain
if isinstance(e, ufl.finiteelement.AbstractFiniteElement):
cmap = _coordinate_element_from_basix(e.basix_element)
domain = ufl.Mesh(e)
assert domain.geometric_dimension == gdim
return cmap, domain
if isinstance(e, _CoordinateElement):
return e, None
raise ValueError(f"Unsupported element type {type(e)}.")
[docs]
def create_mesh(
comm: _MPI.Comm,
cells: npt.NDArray[np.int64],
e: ufl.Mesh | basix.finite_element.FiniteElement | basix.ufl._BasixElement | _CoordinateElement,
x: npt.NDArray[np.floating],
partitioner: _MeshPartitioner = None,
ghost_mode: GhostMode = GhostMode.none,
max_facet_to_cell_links: int = 2,
num_threads: int = 1,
reorder_fn: Callable | None = None,
) -> Mesh:
"""Create a mesh from topology and geometry arrays.
Args:
comm: MPI communicator to define the mesh on.
cells: Cells of the mesh. ``cells[i]`` are the 'nodes' of
cell ``i``.
e: UFL mesh. The mesh scalar type is determined by the scalar
type of ``e``.
x: Mesh geometry ('node' coordinates), with shape
``(num_nodes, gdim)``.
partitioner: ``None``, a partitioner function, or a pair
``(partitioner, cell_weights)``. Cell weights may be ``None``
or an array of weights associated with each cell in ``cells``.
ghost_mode: Ghost mode used in the mesh partitioning, passed to
``partitioner`` at call time.
max_facet_to_cell_links: Maximum number of cells a facet can
be connected to.
num_threads: Number of threads to use to build mesh. Must be
greater than 0.
reorder_fn: Function called with the local cell dual graph. It
must return an array that maps each cell index to its new
index. If ``None`` (default), reverse Cuthill-McKee ordering
is used.
Note:
If required, the coordinates ``x`` will be cast to the same
scalar type as the domain/element ``e``.
Returns:
A mesh.
"""
partitioner_fn, cell_weights = _get_mesh_partitioner(comm, partitioner)
x = np.asarray(x, order="C")
if x.ndim == 1:
gdim = 1
else:
gdim = x.shape[1]
cmap, domain = _create_mesh_coordinate_element(e, gdim)
x = np.asarray(x, dtype=cmap.dtype, order="C")
cells = np.asarray(cells, dtype=np.int64, order="C")
if cell_weights is not None:
cell_weights = np.asarray(cell_weights, dtype=np.int32, order="C")
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64 = _cpp.mesh.create_mesh(
comm=comm,
cells=cells,
element=cmap._cpp_object, # type: ignore[arg-type]
x=x,
partitioner=partitioner_fn,
ghost_mode=ghost_mode,
max_facet_to_cell_links=max_facet_to_cell_links,
num_threads=num_threads,
cell_weights=cell_weights,
reorder_fn=reorder_fn,
)
return Mesh(msh, domain)
[docs]
def create_submesh(
msh: Mesh, dim: int, entities: npt.NDArray[np.int32]
) -> tuple[Mesh, EntityMap, EntityMap, npt.NDArray[np.int32]]:
"""Create a mesh with specified entities from another mesh.
Args:
msh: Mesh to create the sub-mesh from.
dim: Topological dimension of the entities in ``msh`` to include
in the sub-mesh.
entities: Indices of entities in ``msh`` to include in the
sub-mesh.
Returns:
The (1) sub mesh, (2) entity map, (3) vertex map, and (4) node
map (geometry). Each of the maps maps a local index of the sub mesh
to a local index of ``msh``.
"""
submsh, entity_map, vertex_map, geom_map = _cpp.mesh.create_submesh(
msh._cpp_object, dim, entities
)
submsh_domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
to_string(submsh.topology.cell_type),
submsh.geometry.cmaps[0].degree,
basix.LagrangeVariant(submsh.geometry.cmaps[0].variant),
shape=(submsh.geometry.dim,),
dtype=submsh.geometry.x.dtype,
)
)
return (Mesh(submsh, submsh_domain), EntityMap(entity_map), EntityMap(vertex_map), geom_map)
[docs]
def create_interval(
comm: _MPI.Comm,
nx: int,
points: npt.ArrayLike,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
gdim: int = 1,
) -> Mesh:
"""Create an interval mesh.
Args:
comm: MPI communicator.
nx: Number of cells.
points: Coordinates of the end points.
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``).
ghost_mode: Ghost mode used in the mesh partitioning. Options
are ``GhostMode.none`` and ``GhostMode.shared_facet``.
partitioner: Partitioning function to use for determining the
parallel distribution of cells across MPI ranks.
gdim: Geometric dimension. The interval lies along the first
coordinate axis; remaining components are zero.
Returns:
An interval mesh.
"""
if partitioner is None and comm.size > 1:
partitioner = _cpp.graph.partitioner()
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
"interval",
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(gdim,),
dtype=dtype,
)
)
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner, gdim) # type: ignore[arg-type]
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner, gdim) # type: ignore[arg-type]
else:
raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
return Mesh(msh, domain)
[docs]
def create_unit_interval(
comm: _MPI.Comm,
nx: int,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
gdim: int = 1,
) -> Mesh:
"""Create a mesh on the unit interval.
Args:
comm: MPI communicator.
nx: Number of cells.
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``).
ghost_mode: Ghost mode used in the mesh partitioning. Options
are ``GhostMode.none`` and ``GhostMode.shared_facet``.
partitioner: Partitioning function to use for determining the
parallel distribution of cells across MPI ranks.
gdim: Geometric dimension. The interval lies along the first
coordinate axis; remaining components are zero.
Returns:
A unit interval mesh with end points at 0 and 1.
"""
return create_interval(
comm, nx, [0.0, 1.0], dtype=dtype, ghost_mode=ghost_mode, partitioner=partitioner, gdim=gdim
)
[docs]
def create_rectangle(
comm: _MPI.Comm,
points: npt.ArrayLike,
n: Sequence[int],
cell_type: CellType = CellType.triangle,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
diagonal: DiagonalType = DiagonalType.right,
gdim: int = 2,
) -> Mesh:
"""Create a rectangle mesh.
Args:
comm: MPI communicator.
points: Coordinates of the lower - left and upper - right
corners of the rectangle.
n: Number of cells in each direction.
cell_type: Mesh cell type.
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``)
ghost_mode: Ghost mode used in the mesh partitioning.
partitioner: Function that computes the parallel distribution of
cells across MPI ranks.
diagonal: Direction of diagonal of triangular meshes.
See :class:`DiagonalType` for available options.
gdim: Geometric dimension. The rectangle lies in the first 2
coordinate axes; remaining components are zero.
Returns:
A mesh of a rectangle.
"""
if partitioner is None and comm.size > 1:
partitioner = _cpp.graph.partitioner()
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
cell_type.name,
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(gdim,),
dtype=dtype,
)
)
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_rectangle_float32(
comm,
points, # type: ignore[arg-type]
n,
cell_type,
partitioner,
diagonal,
gdim,
ghost_mode,
)
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_rectangle_float64(
comm,
points, # type: ignore[arg-type]
n,
cell_type,
partitioner,
diagonal,
gdim,
ghost_mode,
)
else:
raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
return Mesh(msh, domain)
[docs]
def create_unit_square(
comm: _MPI.Comm,
nx: int,
ny: int,
cell_type: CellType = CellType.triangle,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
diagonal: DiagonalType = DiagonalType.right,
gdim: int = 2,
) -> Mesh:
"""Create a mesh of a unit square.
Args:
comm: MPI communicator.
nx: Number of cells in the "x" direction.
ny: Number of cells in the "y" direction.
cell_type: Mesh cell type.
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``).
ghost_mode: Ghost mode used in the mesh partitioning.
partitioner: Function that computes the parallel distribution of
cells across MPI ranks.
diagonal:
Direction of diagonal. See :class:`DiagonalType` for
available options.
gdim: Geometric dimension. The square lies in the first 2
coordinate axes; remaining components are zero.
Returns:
A mesh of a square with corners at ``(0, 0)`` and ``(1, 1)``.
"""
return create_rectangle(
comm,
[np.array([0.0, 0.0]), np.array([1.0, 1.0])],
(nx, ny),
cell_type,
dtype=dtype,
ghost_mode=ghost_mode,
partitioner=partitioner,
diagonal=diagonal,
gdim=gdim,
)
[docs]
def create_box(
comm: _MPI.Comm,
points: list[npt.ArrayLike],
n: Sequence[int],
cell_type: CellType = CellType.tetrahedron,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
) -> Mesh:
"""Create a box mesh.
Args:
comm: MPI communicator.
points: Coordinates of the 'lower-left' and 'upper-right'
corners of the box.
n: List of cells in each direction
cell_type: The cell type.
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``).
ghost_mode: The ghost mode used in the mesh partitioning.
partitioner: Function that computes the parallel distribution of
cells across MPI ranks.
Returns:
A mesh of a box domain.
"""
if partitioner is None and comm.size > 1:
partitioner = _cpp.graph.partitioner()
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
cell_type.name,
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(3,),
dtype=dtype,
)
)
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner, ghost_mode) # type: ignore[arg-type]
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner, ghost_mode) # type: ignore[arg-type]
else:
raise RuntimeError(f"Unsupported mesh geometry float type: {dtype}")
return Mesh(msh, domain)
[docs]
def create_unit_cube(
comm: _MPI.Comm,
nx: int,
ny: int,
nz: int,
cell_type: CellType = CellType.tetrahedron,
dtype: npt.DTypeLike = default_real_type,
ghost_mode: GhostMode = GhostMode.shared_facet,
partitioner: Callable | None = None,
) -> Mesh:
"""Create a mesh of a unit cube.
Args:
comm: MPI communicator.
nx: Number of cells in "x" direction.
ny: Number of cells in "y" direction.
nz: Number of cells in "z" direction.
cell_type: Mesh cell type
dtype: Float type for the mesh geometry(``numpy.float32``
or ``numpy.float64``).
ghost_mode: Ghost mode used in the mesh partitioning.
partitioner: Function that computes the parallel distribution of
cells across MPI ranks.
Returns:
A mesh of an axis-aligned unit cube with corners at ``(0, 0, 0)``
and ``(1, 1, 1)``.
"""
return create_box(
comm,
[np.array([0.0, 0.0, 0.0]), np.array([1.0, 1.0, 1.0])],
(nx, ny, nz),
cell_type,
dtype,
ghost_mode,
partitioner,
)
[docs]
def entities_to_geometry(
msh: Mesh, dim: int, entities: npt.NDArray[np.int32], permute: bool = False
) -> npt.NDArray[np.int32]:
"""Geometric DOFs associated with the closure of given mesh entities.
Args:
msh: The mesh.
dim: Topological dimension of the entities of interest.
entities: Entity indices (local to the process).
permute: Permute the DOFs such that they are consistent with the
orientation of `dim`-dimensional mesh entities. This
requires `create_entity_permutations` to be called first.
Returns:
The geometric DOFs associated with the closure of the entities
in `entities`.
"""
return _cpp.mesh.entities_to_geometry(msh._cpp_object, dim, entities, permute)
[docs]
def exterior_facet_indices(topology: Topology) -> npt.NDArray[np.int32]:
"""Compute the indices of exterior facets that are owned by the caller.
An exterior facet (co-dimension 1) is one that is connected globally
to only one cell (co-dimension 0).
Note:
This is a collective operation that should be called on all
processes.
Args:
topology: The topology
Returns:
Sorted list of owned facet indices that are exterior facets of
the mesh.
"""
return _cpp.mesh.exterior_facet_indices(topology._cpp_object)
[docs]
def create_geometry(
index_map: _IndexMap,
dofmap: npt.NDArray[np.int32],
element: _CoordinateElement,
x: np.ndarray,
input_global_indices: npt.NDArray[np.int64],
) -> Geometry:
"""Create a Geometry object.
Args:
index_map: Index map describing the layout of the geometry
points (nodes).
dofmap: The geometry (point) dofmap. For a cell, it gives the
row in the point coordinates ``x`` of each local geometry
node. ``shape=(num_cells, num_dofs_per_cell)``.
element: Element that describes the cell geometry map.
x: The point coordinates. The shape is
``(num_points, geometric_dimension).``
input_global_indices: The 'global' input index of each point,
commonly from a mesh input file.
"""
if x.dtype not in (np.float32, np.float64):
raise ValueError(f"Unknown floating type for geometry, got: {x.dtype}")
if (dtype := np.dtype(element.dtype)) != x.dtype:
raise ValueError(f"Mismatch in x dtype ({x.dtype}) and coordinate element ({dtype})")
cpp_element = element._cpp_object
if isinstance(cpp_element, _cpp.fem.CoordinateElement_float64):
return Geometry(
_cpp.mesh.Geometry_float64(index_map, dofmap, cpp_element, x, input_global_indices)
)
elif isinstance(cpp_element, _cpp.fem.CoordinateElement_float32):
return Geometry(
_cpp.mesh.Geometry_float32(index_map, dofmap, cpp_element, x, input_global_indices)
)
else:
raise ValueError(f"Unknown floating type for coordinate element, got: {dtype}")
[docs]
def create_point_mesh(comm: _MPI.Comm, points: npt.NDArray[np.float32 | np.float64]) -> Mesh:
"""Create a mesh consisting of points only.
Note:
No points are shared between processes.
Args:
comm: MPI communicator to create the mesh on.
points: Points local to the process in the mesh.
"""
# Create topology which only has a 0->0 connectivity and dim 0 entities
cells = np.arange(points.shape[0], dtype=np.int32).reshape(-1, 1)
num_nodes_local = cells.shape[0]
imap = _cpp.common.IndexMap(comm, num_nodes_local)
local_range = imap.local_range[0]
igi = np.arange(num_nodes_local, dtype=np.int64) + local_range
topology = _cpp.mesh.Topology(
cell_type=_cpp.mesh.CellType.point,
vertex_map=imap,
cell_map=imap,
cells=_cpp.graph.AdjacencyList_int32(cells),
original_index=igi,
)
e = basix.ufl.element("Lagrange", "point", 0, shape=(points.shape[1],), dtype=points.dtype)
c_el = _coordinate_element_from_basix(e.basix_element)
geometry = create_geometry(imap, cells, c_el, points, igi)
cpp_mesh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64
cpp_geometry = geometry._cpp_object
if isinstance(cpp_geometry, _cpp.mesh.Geometry_float64):
cpp_mesh = _cpp.mesh.Mesh_float64(comm, topology, cpp_geometry)
elif isinstance(cpp_geometry, _cpp.mesh.Geometry_float32):
cpp_mesh = _cpp.mesh.Mesh_float32(comm, topology, cpp_geometry)
else:
raise RuntimeError(f"Unsupported dtype for mesh {points.dtype}")
return Mesh(cpp_mesh, domain=ufl.Mesh(e))