"""
Functions for computing mesh quality measures.
"""
import os
import firedrake
import ufl
from firedrake.__future__ import interpolate
from firedrake.petsc import PETSc
from pyop2 import op2
from pyop2.utils import get_petsc_dir
PETSC_DIR, PETSC_ARCH = get_petsc_dir()
include_dir = ["%s/include/eigen3" % PETSC_ARCH]
__all__ = ["QualityMeasure"]
[docs]
class QualityMeasure:
"""
Class for computing quality measures associated with a given mesh.
Choices of quality measure:
* ``min_angle``: the minimum angle of each cell
* ``area``: the area of each cell in a 2D triangular mesh
* ``volume``: the volume of each cell in a 3D tetrahedral mesh
* ``facet_area``: the area of each *facet*.
* ``aspect_ratio``: the aspect ratio of each cell
* ``eskew``: the equiangle skew of each cell
* ``skewness``: the skewness of each cell in a 2D triangular mesh
* ``scaled_jacobian``: the scaled Jacobian of each cell
* ``metric``: given a Riemannian metric, this function outputs the
value of the quality measure :eq:`Q_M` based on the transformation
encoded by the metric.
"""
_measures = (
"min_angle",
"area",
"volume",
"facet_area",
"aspect_ratio",
"eskew",
"skewness",
"scaled_jacobian",
"metric",
)
@PETSc.Log.EventDecorator()
def __init__(self, mesh, metric=None, python=False):
"""
:arg mesh: the input mesh to do computations on
:arg metric: the tensor field representing the metric space transformation
:kwarg python: compute the measure using Python?
"""
self.mesh = mesh
self.metric = metric
self.python = python
self.dim = mesh.topological_dimension()
self.coords = mesh.coordinates
self.P0 = firedrake.FunctionSpace(mesh, "DG", 0)
src_dir = os.path.join(os.path.dirname(__file__), "cxx")
self.fname = os.path.join(src_dir, f"quality{self.dim}d.cxx")
def _get_dats(self, func):
dats = (
func.dat(op2.WRITE, func.cell_node_map()),
self.coords.dat(op2.READ, self.coords.cell_node_map()),
)
if self.metric is not None:
dats += (self.metric.dat(op2.READ, self.metric.cell_node_map()),)
return dats
@PETSc.Log.EventDecorator()
def __call__(self, name):
if name not in QualityMeasure._measures:
raise ValueError(f"Quality measure '{name}' not recognised.")
msg = (
f"Quality measure '{name}' not implemented in the {self.dim}D case in C++."
)
if self.python:
return self._call_python(name)
elif name == "facet_area":
raise NotImplementedError(msg)
elif name == "skewness" and self.dim == 3:
raise NotImplementedError(msg)
with open(self.fname, "r") as f:
code = f.read()
func = firedrake.Function(self.P0, name=name)
kwargs = {"cpp": True, "include_dirs": include_dir}
kernel = op2.Kernel(code, f"get_{name}", **kwargs)
op2.par_loop(kernel, self.mesh.cell_set, *self._get_dats(func))
return func
@PETSc.Log.EventDecorator()
def _call_python(self, name):
if name in ("area", "volume"):
return firedrake.assemble(interpolate(ufl.CellVolume(self.mesh), self.P0))
elif name == "facet_area":
HDivTrace = firedrake.FunctionSpace(self.mesh, "HDiv Trace", 0)
v = firedrake.TestFunction(HDivTrace)
u = firedrake.TrialFunction(HDivTrace)
facet_area = firedrake.Function(HDivTrace, name="Facet areas")
mass_term = v("+") * u("+") * ufl.dS + v * u * ufl.ds
rhs = (
v("+") * ufl.FacetArea(self.mesh) * ufl.dS
+ v * ufl.FacetArea(self.mesh) * ufl.ds
)
sp = {
"snes_type": "ksponly",
"ksp_type": "preonly",
"pc_type": "jacobi",
}
firedrake.solve(mass_term == rhs, facet_area, solver_parameters=sp)
return facet_area
elif name == "aspect_ratio" and self.dim == 2:
P0_ten = firedrake.TensorFunctionSpace(self.mesh, "DG", 0)
J = firedrake.assemble(interpolate(ufl.Jacobian(self.mesh), P0_ten))
edge1 = ufl.as_vector([J[0, 0], J[1, 0]])
edge2 = ufl.as_vector([J[0, 1], J[1, 1]])
edge3 = edge1 - edge2
a = ufl.sqrt(ufl.dot(edge1, edge1))
b = ufl.sqrt(ufl.dot(edge2, edge2))
c = ufl.sqrt(ufl.dot(edge3, edge3))
ar = firedrake.Function(self.P0)
ar.interpolate(a * b * c / ((a + b - c) * (b + c - a) * (c + a - b)))
return ar
elif name == "scaled_jacobian" and self.dim == 2:
P0_ten = firedrake.TensorFunctionSpace(self.mesh, "DG", 0)
J = firedrake.assemble(interpolate(ufl.Jacobian(self.mesh), P0_ten))
edge1 = ufl.as_vector([J[0, 0], J[1, 0]])
edge2 = ufl.as_vector([J[0, 1], J[1, 1]])
edge3 = edge1 - edge2
a = ufl.sqrt(ufl.dot(edge1, edge1))
b = ufl.sqrt(ufl.dot(edge2, edge2))
c = ufl.sqrt(ufl.dot(edge3, edge3))
detJ = ufl.JacobianDeterminant(self.mesh)
jacobian_sign = ufl.sign(detJ)
max_product = ufl.max_value(
ufl.max_value(ufl.max_value(a * b, a * c), ufl.max_value(b * c, b * a)),
ufl.max_value(c * a, c * b),
)
return firedrake.assemble(
interpolate(detJ / max_product * jacobian_sign, self.P0)
)
else:
raise NotImplementedError(
f"Quality measure '{name}' not implemented in the {self.dim}D case in Python."
)