import numpy as np
import warnings
from fractions import Fraction
from collections import defaultdict
from collections.abc import Sequence
import pyop2
from pyop2.datatypes import IntType
import petsctools
import firedrake
from functools import cached_property
from firedrake import utils
from firedrake.cython import mgimpl as impl
import firedrake.cython.dmcommon as dmcommon
from .utils import set_level
__all__ = ("HierarchyBase", "MeshHierarchy", "ExtrudedMeshHierarchy", "NonNestedHierarchy",
"SemiCoarsenedExtrudedHierarchy", "SubmeshHierarchy")
def make_unoverlapped_dm(dm):
"""Effectively invert dm.distributeOverlap().
The resulting plex has the identical data structure as the one before
distributeOverlap(). This is algorithmically guaranteed.
"""
tdim = dm.getDimension()
dm = dmcommon.submesh_create(dm, tdim, "depth", tdim, ignore_label_halo=True)
dm.removeLabel("pyop2_core")
dm.removeLabel("pyop2_owned")
dm.removeLabel("pyop2_ghost")
return dm
[docs]
class HierarchyBase(object):
"""Create an encapsulation of an hierarchy of meshes.
Parameters
----------
meshes :
List of meshes (coarse to fine).
coarse_to_fine_cells :
Optional dictionary of numpy arrays keyed by level pair, mapping each
coarse cell to its fine children. Every row is as wide as the busiest
coarse cell's count. After adaptive refinement, a coarse cell that
was not refined has fewer children, so its row is right-padded with
-1. Here -1 is only padding. Omit this and ``fine_to_coarse_cells``
together to derive both maps from ``fine_to_coarse_points``.
fine_to_coarse_cells :
Optional dictionary of numpy arrays keyed by level, mapping each fine
cell to its coarse parent cell. Here -1 marks a fine cell that has
no parent in the coarse mesh. This happens in a `SubmeshHierarchy`
that contains interior facets. A fine facet inside a coarse cell
comes from that volume cell, which is not a cell of the coarse
submesh. Omit this and ``coarse_to_fine_cells`` together to derive
both maps from ``fine_to_coarse_points``.
refinements_per_level :
Number of mesh refinements each multigrid level should "see".
nested :
Is this mesh hierarchy nested?
fine_to_coarse_points :
Dict of numpy arrays for each level that was refined from the level
below it. Each array maps fine DMPlex points to their coarse DMPlex
source points, or to -1 when no corresponding coarse point is present.
Cell maps derived from this point map follow Firedrake cell numbering.
Notes
-----
Most of the time, you do not need to create this object yourself, instead
using `MeshHierarchy`, `ExtrudedMeshHierarchy`, or `NonNestedHierarchy`.
"""
def __init__(self, meshes, coarse_to_fine_cells=None, fine_to_coarse_cells=None,
refinements_per_level=1, nested=False,
fine_to_coarse_points=None):
petsctools.cite("Mitchell2016")
self._meshes = list(meshes)
self.meshes = self._meshes[::refinements_per_level]
self.refinements_per_level = refinements_per_level
self.nested = nested
self.fine_to_coarse_points = dict(fine_to_coarse_points or {})
if (coarse_to_fine_cells is None) != (fine_to_coarse_cells is None):
raise ValueError("coarse_to_fine_cells and fine_to_coarse_cells must be provided together")
if coarse_to_fine_cells is None:
coarse_to_fine_cells = {}
fine_to_coarse_cells = {Fraction(0, 1): None}
for i, (coarse, fine) in enumerate(zip(self._meshes[:-1], self._meshes[1:])):
c2f, f2c = impl.coarse_to_fine_cells(
coarse, fine, self.fine_to_coarse_points[Fraction(i+1, refinements_per_level)])
coarse_to_fine_cells[Fraction(i, refinements_per_level)] = c2f
fine_to_coarse_cells[Fraction(i+1, refinements_per_level)] = f2c
self.coarse_to_fine_cells = coarse_to_fine_cells
self.fine_to_coarse_cells = fine_to_coarse_cells
for level, m in enumerate(meshes):
set_level(m, self, Fraction(level, refinements_per_level))
for level, m in enumerate(self):
set_level(m, self, level)
self._shared_data_cache = defaultdict(dict)
[docs]
@cached_property
def comm(self):
comm = self[0].comm
if not all(m.comm == comm for m in self):
raise NotImplementedError("All meshes in hierarchy must be on same communicator")
return comm
def __iter__(self):
"""Iterate over the hierarchy of meshes from coarsest to finest"""
for m in self.meshes:
yield m
def __len__(self):
"""Return the size of hierarchy"""
return len(self.meshes)
def __getitem__(self, idx):
"""Return a mesh in the hierarchy
:arg idx: The :func:`~.Mesh` to return"""
return self.meshes[idx]
[docs]
def add_mesh(self, mesh, coarse_to_fine_cells=None, fine_to_coarse_cells=None):
"""Add a mesh on top of the finest level of the hierarchy.
Only supported for hierarchies with ``refinements_per_level == 1``.
Parameters
----------
mesh :
The mesh to add, usually obtained by calling
:meth:`~firedrake.mesh.MeshGeometry.refine_marked_elements` on the
current finest mesh.
coarse_to_fine_cells :
Map from the cells of the current finest mesh to the cells of
``mesh``. Defaults to the map ``mesh`` recorded when it was
adaptively refined.
fine_to_coarse_cells :
Map from the cells of ``mesh`` to the cells of the current finest
mesh. Defaults the same way as ``coarse_to_fine_cells``.
Returns
-------
MeshGeometry
The mesh that was added.
"""
if self.refinements_per_level != 1:
raise NotImplementedError("Cannot add a mesh to a hierarchy with "
"refinements_per_level > 1")
fine_to_coarse_points = None
if mesh._adaptive_parent is self[-1]:
fine_to_coarse_points = mesh._adaptive_fine_to_coarse_points
if coarse_to_fine_cells is None or fine_to_coarse_cells is None:
if fine_to_coarse_points is not None:
coarse_to_fine_cells, fine_to_coarse_cells = impl.coarse_to_fine_cells(
self[-1], mesh, fine_to_coarse_points)
elif self.nested:
raise ValueError("Expecting a mesh adaptively refined from the finest "
"level of this hierarchy, or explicit cell maps")
level = len(self.meshes)
self._meshes.append(mesh)
self.meshes.append(mesh)
set_level(mesh, self, level)
mesh.topology_dm.setRefineLevel(level)
self.coarse_to_fine_cells[Fraction(level - 1, 1)] = coarse_to_fine_cells
self.fine_to_coarse_cells[Fraction(level, 1)] = fine_to_coarse_cells
self.fine_to_coarse_points[Fraction(level, 1)] = fine_to_coarse_points
return mesh
[docs]
def adapt(self, eta, theta: float):
"""Add a new mesh to the hierarchy by locally refining the finest mesh
with a simplified variant of Dorfler marking.
Parameters
----------
eta :
A DG0 :class:`~firedrake.function.Function` with the local error estimator.
theta :
The threshold for marking as a fraction of the maximum error.
Returns
-------
MeshGeometry
The mesh that was added.
Note
----
Dorfler marking involves sorting all of the elements by decreasing
error estimator and taking the minimal set that exceeds some fixed
fraction of the total error. What this code implements is the simpler
variant that doesn't have a proof of convergence (as far as I know)
but works as well in practice.
"""
if not isinstance(eta, (firedrake.Function, firedrake.Cofunction)):
raise TypeError(f"eta must be a Function or Cofunction, not a {type(eta).__name__}")
M = eta.function_space()
if M.finat_element.space_dimension() != 1:
raise ValueError("eta must be a Function or Cofunction in DG0")
mesh = self[-1]
if M.mesh() is not mesh:
raise ValueError("eta must be defined on the finest mesh of the hierarchy")
# Take the maximum over all processes
with eta.dat.vec_ro as evec:
_, eta_max = evec.max()
markers = firedrake.Function(M)
markers.dat.data_wo[eta.dat.data_ro > theta * eta_max] = 1
return self.add_mesh(mesh.refine_marked_elements(markers))
[docs]
def MeshHierarchy(mesh, refinement_levels=0,
refinements_per_level=1,
netgen_flags=None,
reorder=None,
distribution_parameters=None, callbacks=None,
mesh_builder=firedrake.Mesh, nested=True):
"""Build a hierarchy of meshes by uniformly refining a coarse mesh.
Parameters
----------
mesh : MeshGeometry
the coarse mesh to refine
refinement_levels : int
the number of levels of uniform refinement. This may be dynamically
increased by :meth:`HierarchyBase.adapt` or
:meth:`HierarchyBase.add_mesh`.
refinements_per_level : int
the number of refinements for each level in the hierarchy.
Adaptive refinement only supports one refinement per level.
netgen_flags : dict or None
Options for a mesh generated by Netgen. If ``mesh`` was generated by
Netgen and this value is ``None``, the hierarchy reuses its Netgen
flags. The vertices of each refined mesh are snapped onto the Netgen
geometry. The ``"degree"`` option sets the degree of the curved
coordinates, either as an integer or as a sequence with one entry per
level in the hierarchy. The ``"cg"`` option sets their continuity.
distribution_parameters : dict
options controlling mesh distribution, see :py:func:`.Mesh`
for details. If ``None``, use the same distribution
parameters as were used to distribute the coarse mesh,
otherwise, these options override the default.
reorder : bool
optional flag indicating whether to reorder the
refined meshes.
callbacks : tuple
A 2-tuple of callbacks to call before and
after refinement of the DM. The before callback receives
the DM to be refined (and the current level), the after
callback receives the refined DM (and the current level).
mesh_builder
Function to turn a DM into a ``Mesh``. Used by pyadjoint.
nested : bool
Are the meshes added to this hierarchy required to be nested? If
`False`, :meth:`HierarchyBase.add_mesh` accepts a mesh that was not
adaptively refined from the finest level.
Returns
-------
HierarchyBase
The mesh hierarchy.
"""
nlevels = refinement_levels + 1
is_netgen = hasattr(mesh, "netgen_mesh")
if netgen_flags is not None and not isinstance(netgen_flags, dict):
raise TypeError("netgen_flags must be a dictionary or None")
if netgen_flags is not None and not is_netgen:
raise RuntimeError("Cannot pass netgen_flags to a mesh that was not generated by Netgen.")
if is_netgen:
utils.check_netgen_installed()
from firedrake.netgen import _snap_to_netgen, _curve_netgen_mesh
if netgen_flags is None:
netgen_flags = mesh.netgen_flags
degree = netgen_flags.get("degree", 1)
if isinstance(degree, int):
degree = (degree,) * nlevels
if len(degree) != nlevels:
raise ValueError(f"Expecting one coordinate degree per level in the hierarchy ({nlevels}), "
f"got {len(degree)}.")
cg_field = netgen_flags.get("cg")
ngmeshes = [mesh.netgen_mesh]
mesh = _curve_netgen_mesh(mesh, degree[0], cg_field=cg_field)
if callbacks is not None:
before, after = callbacks
else:
before = after = lambda dm, i: None
# Refine an unoverlapped plex at each level. Keeping every dm here
# unoverlapped means overlap only ever needs to be added once, by
# mesh_builder below.
cdm = mesh.topology_dm
if refinement_levels > 0:
cdm = make_unoverlapped_dm(cdm)
cdm.setRefinementUniform(True)
dms = [cdm]
for i in range(refinement_levels*refinements_per_level):
if i % refinements_per_level == 0:
before(cdm, i)
cdm.setSaveTransform()
rdm = cdm.refine()
if i % refinements_per_level == 0:
after(rdm, i)
if is_netgen:
ngmeshes.append(_snap_to_netgen(rdm, mesh.netgen_mesh))
# Fix up coords if refining embedded circle or sphere
if hasattr(mesh, '_radius'):
# FIXME, really we need some CAD-like representation
# of the boundary we're trying to conform to. This
# doesn't DTRT really for cubed sphere meshes (the
# refined meshes are no longer gnonomic).
coords = rdm.getCoordinatesLocal().array.reshape(-1, mesh.geometric_dimension)
scale = mesh._radius / np.linalg.norm(coords, axis=1).reshape(-1, 1)
coords *= scale
dms.append(rdm)
cdm = rdm
# mesh_builder adds overlap to each DMPlex in place, so first capture the
# lgmaps and read the fine-to-coarse point map from each unoverlapped DM.
lgmaps = [impl.create_lgmap(dm) for dm in dms]
points = [impl.transform_source_points(dm) for dm in dms[1:]]
# Build a mesh for each level, adding overlap here.
parameters = {}
if distribution_parameters is not None:
parameters.update(distribution_parameters)
else:
parameters.update(mesh._distribution_parameters)
parameters["partition"] = False
meshes = [mesh]
for rdm in dms[1:]:
fmesh = mesh_builder(
rdm,
dim=mesh.geometric_dimension,
distribution_parameters=parameters,
reorder=reorder,
comm=mesh.comm,
)
meshes.append(fmesh)
with pyop2.mpi.temp_internal_comm(meshes[0].comm) as icomm:
num_halo_cells = icomm.allreduce(
sum(m.cell_set.total_size - m.cell_set.size for m in meshes)
)
if num_halo_cells == 0:
lgmaps = [None] * len(dms)
if is_netgen:
for i in range(1, len(meshes)):
meshes[i].netgen_mesh = ngmeshes[i]
meshes[i].netgen_flags = netgen_flags
level, remainder = divmod(i, refinements_per_level)
if remainder == 0:
meshes[i] = _curve_netgen_mesh(meshes[i], degree[level], cg_field=cg_field)
for i, m in enumerate(meshes):
# Firedrake counts multigrid levels, PETSc counts refinements
m.topology_dm.setRefineLevel(i)
fine_to_coarse_points = {
Fraction(i+1, refinements_per_level): impl.overlapped_fine_to_coarse_points(
meshes[i], meshes[i+1], points[i], lgmaps[i], lgmaps[i+1])
for i in range(len(points))
}
return HierarchyBase(meshes, refinements_per_level=refinements_per_level,
nested=nested, fine_to_coarse_points=fine_to_coarse_points)
[docs]
def ExtrudedMeshHierarchy(base_hierarchy: HierarchyBase,
height: float,
base_layer: int = -1,
refinement_ratio: int = 2,
layers: Sequence[int] | None = None,
kernel=None,
extrusion_type: str = 'uniform',
periodic: bool = False,
gdim: int | None = None,
mesh_builder=firedrake.ExtrudedMesh
) -> HierarchyBase:
"""Build a hierarchy of extruded meshes by extruding a hierarchy of meshes.
Parameters
----------
base_hierarchy : HierarchyBase
The unextruded base mesh hierarchy to extrude.
height : float
The height of the domain to extrude to. This is in contrast to the
extrusion routines, which take ``layer_height``, the height of an
individual layer. The height of an individual layer varies when the
hierarchy is refined in the extruded dimension.
base_layer : int
The number of layers to use for the coarsest grid.
refinement_ratio : int
The ratio by which to increase ``base_layer`` at every refinement. A
value of 2 gives standard uniform refinement, and a value of 1 does
not refine in the extruded dimension, which gives a semicoarsened
hierarchy.
layers : Sequence[int] or None
An alternative to specifying ``base_layer`` and ``refinement_ratio``.
This sequence gives the number of layers at each level in the
extruded hierarchy. Do not combine this option with ``base_layer`` and
``refinement_ratio``. The ratio of successive entries must be an
integer for the multigrid transfer operators to work.
kernel : pyop2.op2.Kernel or None
An optional kernel that computes the coordinates of the extruded mesh.
extrusion_type : str
The algorithm that computes the extruded coordinates. Supported
values include ``"uniform"``, ``"radial"``,
``"radial_hedgehog"``, and ``"custom"``.
periodic : bool
Whether to identify the top and bottom boundaries of each extruded
mesh. Periodic extrusion requires a constant number of layers.
gdim : int or None
The number of spatial dimensions in the resulting mesh. This value is
used only when ``kernel`` is provided.
mesh_builder : collections.abc.Callable
The function that turns a ``Mesh`` into an extruded mesh. This is used
by pyadjoint.
Returns
-------
HierarchyBase
The hierarchy of extruded meshes.
"""
if not isinstance(base_hierarchy, HierarchyBase):
raise ValueError("Expecting a HierarchyBase, not a %r" % type(base_hierarchy))
if any(m.cell_set._extruded for m in base_hierarchy):
raise ValueError("Meshes in base hierarchy must not be extruded")
if layers is None:
if base_layer == -1:
raise ValueError("Must specify number of layers for coarsest grid with base_layer=N")
layers = [base_layer * refinement_ratio**idx for idx in range(len(base_hierarchy._meshes))]
else:
if base_layer != -1:
raise ValueError("Can't specify both layers and base_layer")
meshes = [mesh_builder(m, layer, kernel=kernel,
layer_height=height/layer,
extrusion_type=extrusion_type,
periodic=periodic,
gdim=gdim)
for (m, layer) in zip(base_hierarchy._meshes, layers)]
return HierarchyBase(meshes,
base_hierarchy.coarse_to_fine_cells,
base_hierarchy.fine_to_coarse_cells,
refinements_per_level=base_hierarchy.refinements_per_level,
nested=base_hierarchy.nested)
[docs]
def SemiCoarsenedExtrudedHierarchy(base_mesh, height, nref=1, base_layer=-1, refinement_ratio=2, layers=None,
kernel=None, extrusion_type='uniform', gdim=None,
mesh_builder=firedrake.ExtrudedMesh):
"""Build a hierarchy of extruded meshes with refinement only in the extruded dimension.
:arg base_mesh: the unextruded base mesh to extrude.
:arg nref: Number of refinements.
:arg height: the height of the domain to extrude to. This is in contrast
to the extrusion routines, which take in layer_height, the height of
an individual layer. This is because when refining in the extruded
dimension, the height of an individual layer will vary.
:arg base_layer: the number of layers to use the extrusion of the coarsest
grid.
:arg refinement_ratio: the ratio by which base_layer should be increased
on every refinement. refinement_ratio = 2 means standard uniform
refinement. refinement_ratio = 1 means to not refine in the extruded
dimension, i.e. the multigrid hierarchy will use semicoarsening.
:arg layers: as an alternative to specifying base_layer and refinement_ratio,
one may specify directly the number of layers to be used by each level
in the extruded hierarchy. This option cannot be combined with base_layer
and refinement_ratio. Note that the ratio of successive entries in this
iterable must be an integer for the multigrid transfer operators to work.
:arg mesh_builder: function used to turn a ``Mesh`` into an
extruded mesh. Used by pyadjoint.
See :func:`~.ExtrudedMesh` for the meaning of the remaining parameters.
See also :func:`~.ExtrudedMeshHierarchy` if you want to extruded a
hierarchy of unstructured meshes.
"""
if not isinstance(base_mesh, firedrake.mesh.MeshGeometry):
raise ValueError(f"Can only extruded a mesh, not a {type(base_mesh)}")
if base_mesh.cell_set._extruded:
raise ValueError("Base mesh must not be extruded")
if layers is None:
if base_layer == -1:
raise ValueError("Must specify number of layers for coarsest grid with base_layer=N")
layers = [base_layer * refinement_ratio**idx for idx in range(nref+1)]
else:
if base_layer != -1:
raise ValueError("Can't specify both layers and base_layer")
if len(layers) == nref+1:
raise ValueError("Need to provide a number of layers for every refined mesh. "
f"Got {len(layers)}, needed {nref+1}")
meshes = [mesh_builder(base_mesh, layer, kernel=kernel,
layer_height=height/layer,
extrusion_type=extrusion_type,
gdim=gdim)
for layer in layers]
refinements_per_level = 1
identity = np.arange(base_mesh.cell_set.size, dtype=IntType).reshape(-1, 1)
coarse_to_fine_cells = dict((Fraction(i, refinements_per_level), identity)
for i in range(nref))
fine_to_coarse_cells = dict((Fraction(i+1, refinements_per_level), identity)
for i in range(nref))
return HierarchyBase(meshes, coarse_to_fine_cells, fine_to_coarse_cells,
refinements_per_level=refinements_per_level,
nested=True)
def AdaptiveMeshHierarchy(*args, **kwargs):
"""Deprecated alias for `MeshHierarchy`."""
warnings.warn(
"The ``AdaptiveMeshHierarchy`` constructor is deprecated and will be removed in a future "
"release. Please use the ``MeshHierarchy`` constructor instead.", FutureWarning
)
return MeshHierarchy(*args, **kwargs)
[docs]
def NonNestedHierarchy(*meshes):
return HierarchyBase(meshes, [None for _ in meshes], [None for _ in meshes],
nested=False)
[docs]
def SubmeshHierarchy(parent_hierarchy: HierarchyBase,
subdim: int | None = None,
subdomain_id: int | Sequence | None = None,
label_name: str | None = None,
name: str | None = None,
ignore_halo: bool = False,
reorder: bool | None = None,
comm=None):
"""Build a hierarchy of submeshes from an existing mesh hierarchy.
Parameters
----------
parent_hierarchy :
Parent mesh hierarchy.
subdim :
Topological dimension of the submesh.
Defaults to ``mesh.topological_dimension``.
subdomain_id :
Subdomain ID representing the submesh.
If `None` the submesh will cover the entire domain.
This is useful to obtain a codim-1 submesh over all facets or
a submesh over a different communicator.
label_name :
Name of the label to search ``subdomain_id`` in.
Defaults to ``'Cell Sets'`` or ``'Face Sets'`` depending on ``subdim``.
name : str | None
Name of the submesh.
Defaults to ``mesh.name + "_submesh"``ยท
ignore_halo :
Whether to exclude the halo from the submesh.
reorder :
Whether to reorder the mesh entities. By default,
the submesh will be reordered if the parent mesh was reordered.
comm : PETSc.Comm | None
An optional sub-communicator to define the submesh.
By default, the submesh is defined on `mesh.comm`.
Returns
-------
HierarchyBase
The submesh hierarchy.
Notes
-----
A facet submesh that contains interior facets, such as the skeleton
obtained with ``subdim=tdim-1`` and ``subdomain_id=None``, is only
partially nested. A fine facet inside a coarse cell has no parent in the
coarse submesh.
Injection is still possible, because the fine children of each coarse
facet cover it. `prolong` and `restrict` raise `NotImplementedError`
when a fine node lies only on facets without a parent.
"""
meshes = [firedrake.Submesh(mesh,
subdim=subdim, subdomain_id=subdomain_id,
label_name=label_name, name=name,
ignore_halo=ignore_halo, reorder=reorder,
comm=comm)
for mesh in parent_hierarchy._meshes]
refinements_per_level = parent_hierarchy.refinements_per_level
fine_to_coarse_points = {}
for i, (coarse, fine) in enumerate(zip(meshes[:-1], meshes[1:])):
level = Fraction(i+1, refinements_per_level)
parent_points = parent_hierarchy.fine_to_coarse_points.get(level)
if parent_points is None:
raise ValueError("SubmeshHierarchy requires a hierarchy whose DMPlex "
f"points on level {level} are refined from the level below")
# Map fine submesh points through the parent hierarchy, then translate
# them into coarse-submesh numbering. A source outside the coarse
# submesh remains -1; for example, a new interior facet can have a
# coarse volume cell as its source.
coarse_subpoints = coarse.topology_dm.getSubpointIS().indices
parent_to_coarse = np.full(parent_hierarchy._meshes[i].topology_dm.getChart()[1], -1, dtype=IntType)
parent_to_coarse[coarse_subpoints] = np.arange(len(coarse_subpoints), dtype=IntType)
points = impl.compose_points(parent_points, fine.topology_dm.getSubpointIS().indices)
fine_to_coarse_points[level] = impl.compose_points(parent_to_coarse, points)
for i, m in enumerate(meshes):
m.topology_dm.setRefineLevel(i)
return HierarchyBase(meshes, refinements_per_level=refinements_per_level,
nested=parent_hierarchy.nested,
fine_to_coarse_points=fine_to_coarse_points)
