Add affine mapping for edges and EdgeBasis

skfem/assembly/basis/edge_basis.py

@@ -0,0 +1,101 @@
+import logging
+from typing import Optional, Tuple, Any
+
+import numpy as np
+from numpy import ndarray
+from skfem.element import DiscreteField, Element
+from skfem.mapping import Mapping
+from skfem.mesh import Mesh
+from skfem.refdom import RefLine
+
+from .abstract_basis import AbstractBasis
+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,
+        )

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,

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.],