From 84bb97949b64fbc532430ea09efadc4bc6ed9bd4 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 4 Oct 2024 10:45:03 +0300 Subject: [PATCH 1/3] Add affine mapping for edges and EdgeBasis --- skfem/assembly/basis/edge_basis.py | 105 +++++++++++++++++++++++++++++ skfem/mapping/mapping_affine.py | 62 +++++++++++++++++ 2 files changed, 167 insertions(+) create mode 100644 skfem/assembly/basis/edge_basis.py diff --git a/skfem/assembly/basis/edge_basis.py b/skfem/assembly/basis/edge_basis.py new file mode 100644 index 00000000..6123c809 --- /dev/null +++ b/skfem/assembly/basis/edge_basis.py @@ -0,0 +1,105 @@ +import logging +from typing import Callable, Optional, Tuple, Any + +import numpy as np +from numpy import ndarray +from skfem.element import (BOUNDARY_ELEMENT_MAP, DiscreteField, Element, + ElementHex0, ElementQuad0, ElementTetP0, + ElementTriP0) +from skfem.mapping import Mapping +from skfem.mesh import Mesh, MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri +from skfem.generic_utils import OrientedBoundary, deprecated +from skfem.refdom import RefLine + +from .abstract_basis import AbstractBasis +from .cell_basis import CellBasis +from ..dofs import Dofs + + +logger = logging.getLogger(__name__) + + +class EdgeBasis(AbstractBasis): + """For integrating over edges of the mesh.""" + + def __init__(self, + mesh: Mesh, + elem: Element, + mapping: Optional[Mapping] = None, + intorder: Optional[int] = None, + quadrature: Optional[Tuple[ndarray, ndarray]] = None, + edges: Optional[Any] = None, + dofs: Optional[Dofs] = None, + side: int = 0, + disable_doflocs: bool = False): + """Precomputed global basis on edges.""" + typestr = ("{}({}, {})".format(type(self).__name__, + type(mesh).__name__, + type(elem).__name__)) + logger.info("Initializing {}".format(typestr)) + super(EdgeBasis, self).__init__( + mesh, + elem, + mapping, + intorder, + quadrature, + RefLine, + dofs, + disable_doflocs, + ) + + # by default use boundary edges (for testing only) + if edges is None: + self.eind = self.mesh.boundary_edges() + else: + self.eind = edges + + # TODO fix the orientation + # element corresponding to the edge + # is the first matching element from e2t + tmp = mesh.e2t[:, self.eind] + self.tind = tmp.indices[tmp.indptr[:-1]] + + if len(self.eind) == 0: + logger.warning("Initializing {} with no edges.".format(typestr)) + + # edge refdom to global edge + x = self.mapping.H(self.X, eind=self.eind) + # global edge to refdom edge + Y = self.mapping.invF(x, tind=self.tind) + + # TODO calculate tangents + self.tangents = ... + + self.nelems = len(self.eind) + + self.basis = [self.elem.gbasis(self.mapping, Y, j, tind=self.tind) + for j in range(self.Nbfun)] + + self.dx = (np.abs(self.mapping.detDH(self.X, eind=self.eind)) + * np.broadcast_to(self.W, (self.nelems, self.W.shape[-1]))) + logger.info("Initializing finished.") + + def default_parameters(self): + """Return default parameters for `~skfem.assembly.asm`.""" + return { + 'x': self.global_coordinates(), + 'h': self.mesh_parameters(), + 't': self.tangents, + } + + def global_coordinates(self) -> DiscreteField: + return DiscreteField(self.mapping.H(self.X, eind=self.eind)) + + def mesh_parameters(self) -> DiscreteField: + return DiscreteField(np.abs(self.mapping.detDH(self.X, self.eind))) + + def with_element(self, elem: Element) -> 'EdgeBasis': + """Return a similar basis using a different element.""" + return type(self)( + self.mesh, + elem, + mapping=self.mapping, + quadrature=self.quadrature, + edges=self.eind, + ) diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py index 150b3b1a..03ef42c5 100644 --- a/skfem/mapping/mapping_affine.py +++ b/skfem/mapping/mapping_affine.py @@ -151,6 +151,27 @@ def detB(self): self._init_boundary_mapping() return self._detB + @property + def C(self): + """Edge affine mapping matrix.""" + if not hasattr(self, '_C'): + self._init_edge_mapping() + return self._C + + @property + def d(self): + """Edge affine mapping constant.""" + if not hasattr(self, '_d'): + self._init_edge_mapping() + return self._d + + @property + def detC(self): + """Edge affine mapping determinant.""" + if not hasattr(self, '_detC'): + self._init_edge_mapping() + return self._detC + def _init_boundary_mapping(self): """For lazy evaluation of boundary mapping.""" dim = self.dim @@ -180,6 +201,26 @@ def _init_boundary_mapping(self): else: raise Exception("Not implemented for the given dimension.") + def _init_edge_mapping(self): + """For lazy evaluation of edge mapping.""" + dim = self.dim + assert dim == 3 + + ne = self.mesh.edges.shape[1] + # initialize the boundary mapping + self._C = np.empty((dim, 1, ne)) + self._d = np.empty((dim, ne)) + + for i in range(dim): + self._d[i] = self.mesh.p[i, self.mesh.edges[0]] + self._C[i, 0] = (self.mesh.p[i, self.mesh.edges[1]] - + self.mesh.p[i, self.mesh.edges[0]]) + + # length scaling + self._detC = np.sqrt(self._C[0, 0] ** 2 + + self._C[1, 0] ** 2 + + self._C[2, 0] ** 2) + def F(self, X, tind=None): if tind is None or self.tind is not None: A, b = self.A, self.b @@ -252,6 +293,27 @@ def detDG(self, X, find=None): return np.tile(detDG, (X.shape[-1], 1)).T + def H(self, X, eind=None): + if eind is None: + C, d = self.C, self.d + else: + C, d = self.C[:, :, eind], self.d[:, eind] + + if len(X.shape) == 2: + return (np.einsum('ijk,jl', C, X).T + d.T).T + elif len(X.shape) == 3: + return (np.einsum('ijk,jkl->ikl', C, X).T + d.T).T + + raise Exception("Wrong dimension of input.") + + def detDH(self, X, eind=None): + if eind is None: + detDH = self.detC + else: + detDH = self.detC[eind] + + return np.tile(detDH, (X.shape[-1], 1)).T + def normals(self, X, tind, find, t2f): if self.dim == 1: Nref = np.array([[-1.], From 243a3cfe9d4d1cafccacf08a6996fd167a4a4ff5 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 4 Oct 2024 10:48:32 +0300 Subject: [PATCH 2/3] fix flake8 --- skfem/assembly/basis/edge_basis.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/skfem/assembly/basis/edge_basis.py b/skfem/assembly/basis/edge_basis.py index 6123c809..2cae84e2 100644 --- a/skfem/assembly/basis/edge_basis.py +++ b/skfem/assembly/basis/edge_basis.py @@ -1,18 +1,14 @@ import logging -from typing import Callable, Optional, Tuple, Any +from typing import Optional, Tuple, Any import numpy as np from numpy import ndarray -from skfem.element import (BOUNDARY_ELEMENT_MAP, DiscreteField, Element, - ElementHex0, ElementQuad0, ElementTetP0, - ElementTriP0) +from skfem.element import DiscreteField, Element from skfem.mapping import Mapping -from skfem.mesh import Mesh, MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri -from skfem.generic_utils import OrientedBoundary, deprecated +from skfem.mesh import Mesh from skfem.refdom import RefLine from .abstract_basis import AbstractBasis -from .cell_basis import CellBasis from ..dofs import Dofs From a1435dacbe661463205f29aa077585b699ea1b95 Mon Sep 17 00:00:00 2001 From: Tom Gustafsson Date: Fri, 4 Oct 2024 11:01:51 +0300 Subject: [PATCH 3/3] add prototypes to satisfy mypy --- skfem/mapping/mapping.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/skfem/mapping/mapping.py b/skfem/mapping/mapping.py index f4d8ac6d..c08369aa 100644 --- a/skfem/mapping/mapping.py +++ b/skfem/mapping/mapping.py @@ -87,6 +87,18 @@ def detDG(self, """ raise NotImplementedError + def H(self, + X: ndarray, + eind: Optional[ndarray] = None) -> ndarray: + """Perform a mapping from the reference to global edge.""" + raise NotImplementedError + + def detDH(self, + X: ndarray, + eind: Optional[ndarray] = None) -> ndarray: + """The jacobian determinant of H.""" + raise NotImplementedError + def normals(self, X: ndarray, tind: ndarray,