# 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 functools import singledispatch
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 import coordinate_element as _coordinate_element
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_incident_entities",
"compute_midpoints",
"create_box",
"create_cell_partitioner",
"create_geometry",
"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",
]
[docs]
@singledispatch
def create_cell_partitioner(
part: Callable, mode: GhostMode, max_facet_to_cell_links: int
) -> Callable:
"""Create a function to partition a mesh.
Args:
part: Partition function.
mode: Ghosting mode to use.
max_facet_to_cell_links: Maximum number of cells connected to a
facet. Equal to 2 for non-branching manifold meshes.
``None`` corresponds to no upper bound on the number of
possible connections.
Return:
Partitioning function.
"""
return _cpp.mesh.create_cell_partitioner(part, mode, max_facet_to_cell_links)
@create_cell_partitioner.register(GhostMode)
def _(mode: GhostMode, max_facet_to_cell_links: int) -> Callable:
"""Create a function to partition a mesh.
Args:
mode: Ghosting mode to use
max_facet_to_cell_links: Maximum number of cells connected to a
facet. Equal to 2 for non-branching manifold meshes.
``None`` corresponds to no upper bound on the number of
possible connections.
Return:
Partitioning function.
"""
return _cpp.mesh.create_cell_partitioner(mode, max_facet_to_cell_links)
[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):
"""Return the MPI communicator associated with the topology."""
return self._cpp_object.comm
[docs]
def create_connectivity(self, d0: int, d1: int):
"""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 = 0) -> bool:
"""Create entities of given topological dimension.
Args:
dim: Topological dimension of entities to create.
num_threads: Number of CPU threads to use when creating. If
0, threads are not spawned.
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):
"""Compute entity permutations and reflections."""
self._cpp_object.create_entity_permutations()
@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):
"""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._cpp_object.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
[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 # type: ignore
@property
def comm(self):
"""MPI communicator associated with the mesh."""
return self._cpp_object.comm
@property
def name(self):
"""Name of the mesh."""
return self._cpp_object.name
@name.setter
def name(self, value):
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 getattr(basix.CellType, self.topology.cell_name())
[docs]
def h(self, dim: int, entities: npt.NDArray[np.int32]) -> npt.NDArray[np.float64]:
"""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)
@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):
"""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, inverse):
"""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):
"""Topological dimension of the entities."""
return self._cpp_object.dim
@property
def topology(self):
"""The topology."""
return self._topology
@property
def sub_topology(self):
"""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, dim: int, entities: npt.NDArray[np.int32]):
"""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)
[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 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:
assert parent_facet is not None
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 on cells or facets.")
[docs]
def refine(
msh: Mesh,
edges: np.ndarray | None = None,
partitioner: Callable | IdentityPartitionerPlaceholder = IdentityPartitionerPlaceholder(),
option: RefinementOption = RefinementOption.parent_cell,
) -> tuple[Mesh, npt.NDArray[np.int32], npt.NDArray[np.int8]]:
"""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.
Returns:
Refined mesh, (optional) parent cells, (optional) parent facets
"""
mesh1, parent_cell, parent_facet = _cpp.refinement.refine(
msh._cpp_object, edges, partitioner, option
)
# 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()) # type: ignore
return Mesh(mesh1, ufl_domain), parent_cell, parent_facet
[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: Callable | None = None,
max_facet_to_cell_links: int = 2,
) -> 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: Function that determines the parallel distribution of
cells across MPI ranks.
max_facet_to_cell_links: Maximum number of cells a facet can
be connected to.
Note:
If required, the coordinates ``x`` will be cast to the same
scalar type as the domain/element ``e``.
Returns:
A mesh.
"""
if partitioner is None and comm.size > 1:
partitioner = create_cell_partitioner(GhostMode.none, 2) # type: ignore
x = np.asarray(x, order="C")
if x.ndim == 1:
gdim = 1
else:
gdim = x.shape[1]
dtype = None
if isinstance(e, ufl.domain.Mesh):
# e is a UFL domain
e_ufl = e.ufl_coordinate_element() # type: ignore
cmap = _coordinate_element(e_ufl.basix_element) # type: ignore
domain = e
dtype = cmap.dtype
# TODO: Resolve UFL vs Basix geometric dimension issue
# assert domain.geometric_dimension() == gdim
elif isinstance(e, basix.finite_element.FiniteElement):
# e is a Basix element
# TODO: Resolve geometric dimension vs shape for manifolds
cmap = _coordinate_element(e) # type: ignore
e_ufl = basix.ufl._BasixElement(e) # type: ignore
e_ufl = basix.ufl.blocked_element(e_ufl, shape=(gdim,))
domain = ufl.Mesh(e_ufl)
dtype = cmap.dtype
assert domain.geometric_dimension == gdim
elif isinstance(e, ufl.finiteelement.AbstractFiniteElement):
# e is a Basix 'UFL' element
cmap = _coordinate_element(e.basix_element) # type: ignore
domain = ufl.Mesh(e)
dtype = cmap.dtype
assert domain.geometric_dimension == gdim
elif isinstance(e, _CoordinateElement):
# e is a CoordinateElement
cmap = e
domain = None
dtype = cmap.dtype # type: ignore
else:
raise ValueError(f"Unsupported element type {type(e)}.")
x = np.asarray(x, dtype=dtype, order="C")
cells = np.asarray(cells, dtype=np.int64, order="C")
msh: _cpp.mesh.Mesh_float32 | _cpp.mesh.Mesh_float64 = _cpp.mesh.create_mesh(
comm, cells, cmap._cpp_object, x, partitioner, max_facet_to_cell_links
)
return Mesh(msh, domain) # type: ignore
[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.shared_facet,
partitioner=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.mesh.create_cell_partitioner(ghost_mode, 2)
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
"interval",
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(gdim,),
dtype=dtype,
)
) # type: ignore
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_interval_float32(comm, nx, points, ghost_mode, partitioner, gdim)
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_interval_float64(comm, nx, points, ghost_mode, partitioner, gdim)
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.shared_facet,
partitioner=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.triangle,
dtype: npt.DTypeLike = default_real_type,
ghost_mode=GhostMode.shared_facet,
partitioner=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.mesh.create_cell_partitioner(ghost_mode, 2)
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
cell_type.name,
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(gdim,),
dtype=dtype,
)
) # type: ignore
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_rectangle_float32(
comm, points, n, cell_type, partitioner, diagonal, gdim
)
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_rectangle_float64(
comm, points, n, cell_type, partitioner, diagonal, gdim
)
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.triangle,
dtype: npt.DTypeLike = default_real_type,
ghost_mode=GhostMode.shared_facet,
partitioner=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.tetrahedron,
dtype: npt.DTypeLike = default_real_type,
ghost_mode=GhostMode.shared_facet,
partitioner=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.mesh.create_cell_partitioner(ghost_mode, 2)
domain = ufl.Mesh(
basix.ufl.element(
"Lagrange",
cell_type.name,
1,
lagrange_variant=basix.LagrangeVariant.unset,
shape=(3,),
dtype=dtype,
)
) # type: ignore
if np.issubdtype(dtype, np.float32):
msh = _cpp.mesh.create_box_float32(comm, points, n, cell_type, partitioner)
elif np.issubdtype(dtype, np.float64):
msh = _cpp.mesh.create_box_float64(comm, points, n, cell_type, partitioner)
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.tetrahedron,
dtype: npt.DTypeLike = default_real_type,
ghost_mode=GhostMode.shared_facet,
partitioner=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=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 of 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 == np.float64:
ftype = _cpp.mesh.Geometry_float64
elif x.dtype == np.float32:
ftype = _cpp.mesh.Geometry_float32
else:
raise ValueError("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})")
return Geometry(ftype(index_map, dofmap, element._cpp_object, x, input_global_indices))
[docs]
def create_point_mesh(comm: _MPI.Intracomm, 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(e.basix_element) # type: ignore[call-arg]
geometry = create_geometry(imap, cells, c_el, points, igi)
if points.dtype == np.float64:
cpp_mesh = _cpp.mesh.Mesh_float64(comm, topology, geometry._cpp_object)
elif points.dtype == np.float32:
cpp_mesh = _cpp.mesh.Mesh_float32(comm, topology, geometry._cpp_object)
else:
raise RuntimeError(f"Unsupported dtype for mesh {points.dtype}")
return Mesh(cpp_mesh, domain=ufl.Mesh(e))