FacetBasis.trace()

skfem/assembly/basis/facet_basis.py

@@ -108,10 +108,59 @@ def default_parameters(self):
                 'h': self.mesh_parameters(),
                 'n': self.normals}
 
-    def global_coordinates(self) -> ndarray:
+    def global_coordinates(self) -> DiscreteField:
         return DiscreteField(self.mapping.G(self.X, find=self.find))
 
-    def mesh_parameters(self) -> ndarray:
+    def mesh_parameters(self) -> DiscreteField:
         return DiscreteField((np.abs(self.mapping.detDG(self.X, self.find))
                               ** (1. / (self.mesh.dim() - 1.)))
                              if self.mesh.dim() != 1 else np.array([0.]))
+
+    def trace(self,
+              x: ndarray,
+              elem: Optional[Element] = None) -> Tuple[ndarray,
+                                                       ndarray,
+                                                       ndarray]:
+        """Restrict solution to a boundary/interface mesh.
+
+        Parameters
+        ----------
+        x
+            The solution vector.
+        elem
+            If given, perform L2 projection onto the new element before
+            restriction.  L2 projection is performed on the boundary.
+
+        Returns
+        -------
+        p
+            The nodes of the trace mesh.
+        t
+            The elements of the trace mesh.
+        y
+            The new (projected and) restricted solution vector.
+
+        """
+        from skfem.utils import project
+
+        cls = type(self)
+        fbasis = self if elem is None else cls(self.mesh,
+                                               elem,
+                                               facets=self.find,
+                                               quadrature=(self.X, self.W))
+        I = fbasis.get_dofs(self.find).all()
+        if len(I) == 0:  # special case: no facet DOFs
+            if fbasis.dofs.interior_dofs.shape[0] > 1:
+                # no one-to-one restriction: requires interpolation
+                raise NotImplementedError
+            # special case: piecewise constant elem
+            I = fbasis.dofs.interior_dofs[:, self.tind].flatten()
+        y = x[I] if elem is None else project(x, self, fbasis, I=I)
+
+        # build connectivity for lower dimensional mesh
+        ix = self.mesh.facets[:, self.find]
+        ixuniq = np.unique(ix)
+        b = np.zeros(np.max(ix) + 1, dtype=np.int64)
+        b[ixuniq] = np.arange(len(ixuniq), dtype=np.int64)
+
+        return self.mesh.p[:, ixuniq], b[ix], y

skfem/element/__init__.py

@@ -56,7 +56,8 @@
 from .element_tet import (ElementTetP0, ElementTetP1, ElementTetP2,
                           ElementTetRT0, ElementTetN0, ElementTetMini,
                           ElementTetCR)
-from .element_hex import ElementHex1, ElementHex2, ElementHexS2  # noqa
+from .element_hex import (ElementHex0, ElementHex1, ElementHex2,
+                          ElementHexS2)  # noqa
 from .element_line import (ElementLineP0, ElementLineP1, ElementLineP2,
                            ElementLinePp, ElementLineHermite,
                            ElementLineMini)  # noqa
@@ -95,6 +96,7 @@
     "ElementTetN0",
     "ElementTetMini",
     "ElementTetCR",
+    "ElementHex0",
     "ElementHex1",
     "ElementHex2",
     "ElementHexS2",

skfem/element/element_hex/__init__.py

@@ -1,4 +1,4 @@
-"""Yadda"""
+from .element_hex0 import ElementHex0  # noqa
 from .element_hex1 import ElementHex1  # noqa
 from .element_hex_s2 import ElementHexS2  # noqa
 from .element_hex2 import ElementHex2  # noqa

skfem/element/element_hex/element_hex0.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+from ..element_h1 import ElementH1
+from ...mesh.mesh3d import MeshHex
+
+
+class ElementHex0(ElementH1):
+
+    interior_dofs = 1
+    dim = 3
+    maxdeg = 0
+    dofnames = ['u']
+    doflocs = np.array([[.5, .5, .5]])
+    mesh_type = MeshHex
+
+    def lbasis(self, X, i):
+        if i == 0:
+            return np.ones(X.shape[1:]), np.zeros_like(X)
+        else:
+            self._index_error()

skfem/element/element_hex/element_hex1.py

@@ -5,6 +5,7 @@
 
 
 class ElementHex1(ElementH1):
+
     nodal_dofs = 1
     dim = 3
     maxdeg = 3

tests/test_basis.py

@@ -1,14 +1,18 @@
 from unittest import TestCase
 
+import pytest
 import numpy as np
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_almost_equal
 
 from skfem import BilinearForm, asm, solve, condense, project
 from skfem.mesh import MeshTri, MeshTet, MeshHex, MeshQuad, MeshLine
-from skfem.assembly import InteriorBasis, FacetBasis, Dofs
+from skfem.assembly import InteriorBasis, FacetBasis, Dofs, Functional
 from skfem.element import (ElementVectorH1, ElementTriP2, ElementTriP1,
                            ElementTetP2, ElementHexS2, ElementHex2,
-                           ElementQuad2, ElementLineP2)
+                           ElementQuad2, ElementLineP2, ElementTriP0,
+                           ElementLineP0, ElementQuad1, ElementQuad0,
+                           ElementTetP1, ElementTetP0, ElementHex1,
+                           ElementHex0)
 
 
 class TestCompositeSplitting(TestCase):
@@ -178,3 +182,33 @@ class TestInterpolatorLine2(TestInterpolatorTet):
     mesh_type = MeshLine
     element_type = ElementLineP2
     nrefs = 5
+
+
+@pytest.mark.parametrize(
+    "mtype,bndmtype,e1,e2,e3",
+    [
+        (MeshTri, MeshLine, ElementTriP1(), ElementTriP0(), ElementLineP0()),
+        (MeshQuad, MeshLine, ElementQuad1(), ElementQuad0(), ElementLineP0()),
+        (MeshTet, MeshTri, ElementTetP1(), ElementTetP0(), ElementTriP0()),
+        (MeshHex, MeshQuad, ElementHex1(), ElementHex0(), ElementQuad0()),
+    ]
+)
+def test_trace(mtype, bndmtype, e1, e2, e3):
+
+    m = mtype()
+    m.refine(3)
+
+    # use the boundary where last coordinate is zero
+    basis = FacetBasis(m, e1,
+                       facets=m.facets_satisfying(lambda x: x[x.shape[0] - 1] == 0.0))
+    p, t, y = basis.trace(m.p[0], e2)
+
+    # drop the last coordinate to create trace mesh
+    nbasis = InteriorBasis(bndmtype(p[0:(p.shape[0] - 1)], t), e3)
+
+    @Functional
+    def integ(w):
+        return w['funy']
+
+    # integrate f(x) = x_1 over trace mesh
+    assert_almost_equal(integ.assemble(nbasis, funy=nbasis.interpolate(y)), .5)