Source code for movement.tangling

import warnings

import firedrake
import ufl
from firedrake.__future__ import interpolate

__all__ = ["MeshTanglingChecker"]


[docs] class MeshTanglingChecker: """ A class for tracking whether a mesh has tangled, i.e. whether any of its elements have become inverted. """ def __init__(self, mesh, raise_error=True): """ :arg mesh: the mesh to track if tangled :type mesh: :class:`firedrake.mesh.MeshGeometry` :kwarg raise_error: if ``True``, an error is raised if any element is tangled, otherwise a warning is raised :type raise_error: :class:`bool` """ if mesh.topological_dimension() != 2: raise NotImplementedError("Tangling check only currently supported in 2D.") self.mesh = mesh self.raise_error = raise_error # Store initial signs of Jacobian determinant P0 = firedrake.FunctionSpace(mesh, "DG", 0) detJ = ufl.JacobianDeterminant(mesh) s = firedrake.assemble(interpolate(ufl.sign(detJ), P0)) # Get scaled Jacobian expression P0_ten = firedrake.TensorFunctionSpace(mesh, "DG", 0) J = firedrake.assemble(interpolate(ufl.Jacobian(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 norm1 = ufl.sqrt(ufl.dot(edge1, edge1)) norm2 = ufl.sqrt(ufl.dot(edge2, edge2)) norm3 = ufl.sqrt(ufl.dot(edge3, edge3)) prod1 = ufl.max_value(norm1 * norm2, norm1 * norm3) prod2 = ufl.max_value(norm2 * norm3, norm2 * norm1) prod3 = ufl.max_value(norm3 * norm1, norm3 * norm2) maxval = ufl.max_value(ufl.max_value(prod1, prod2), prod3) self._sj_expr = detJ / maxval * s self._sj = firedrake.Function(P0) @property def scaled_jacobian(self): self._sj.interpolate(self._sj_expr) return self._sj
[docs] def check(self): """ Check whether any element orientations have changed since the tangling checker was created. """ sj = self.scaled_jacobian.dat.data_with_halos num_tangled = len(sj[sj < 0]) if num_tangled > 0: plural = "s" if num_tangled > 1 else "" msg = f"Mesh has {num_tangled} tangled element{plural}." if self.raise_error: raise ValueError(msg) warnings.warn(msg, stacklevel=1) return num_tangled