Repository: kinnala/scikit-fem

Issue:
FacetBasis and ElementTriRT0 fails
I think there is some issues with broadcasting because of the following failure:
```python
In [1]: from skfem import *
In [2]: FacetBasis(MeshTri(), ElementTriRT0())
...
ValueError: operand has more dimensions than subscripts given in einstein sum, but no '...' ellipsis provided to broadcast the extra dimensions.
```

Discussion:
kinnala: ![Kuvakaappaus - 2022-08-02 11-49-19](https://user-images.githubusercontent.com/973268/182335108-7f29b086-f9b2-42fb-ae9f-dd707d3808f8.png)

New pic for RT example.

Patch:
diff --git a/README.md b/README.md
index 3429e2d9..41adb7f1 100644
--- a/README.md
+++ b/README.md
@@ -212,18 +212,22 @@ with respect to documented and/or tested features.
 
 ### Unreleased
 
-- Changed: The optimization of using `DiscreteField.is_zero` in the helpers to
-  skip the evaluation of zero components in `ElementComposite` was removed to
+- Changed: Removed the optimization of using `DiscreteField.is_zero` in the
+  helpers to skip the evaluation of zero components in `ElementComposite` to
   improve type stability with respect to the size of the underlying numpy
-  arrays
-- Changed: Renamed `ElementTriRT0` to `ElementTriRT1`
-- Changed: Renamed `ElementTetRT0` to `ElementTetRT1`
-- Changed: Renamed `ElementQuadRT0` to `ElementQuadRT1`
-- Removed: The unused `DiscreteField.is_zero` method
+  arrays; this is technically a backwards incompatible change and might affect
+  self-created helper functions
 - Deprecated: `FacetBasis.trace` in favor of `Basis.interpolator` and `Basis.project`
 - Added: Output of `Basis.interpolator` supports trailing axes; can be now
   passed to `Basis.project` for (inexact) interpolation between meshes
+- Added: Renamed `ElementTriRT0` to `ElementTriRT1` and added alias for
+  backwards compatibility
+- Added: Renamed `ElementTetRT0` to `ElementTetRT1` and added alias for
+  backwards compatibility
+- Added: Renamed `ElementQuadRT0` to `ElementQuadRT1` and added alias for
+  backwards compatibility
 - Added: `ElementTriRT2`, the second order Raviart-Thomas element
+- Added: `ElementHexRT1`, the first order Raviart-Thomas element for hexahedral meshes
 - Added: `Basis.project` now better supports `ElementComposite`
 - Added: `solver_iter_cg`, a simple pure Python conjugate gradient solver for
   environments that do not have sparse solver libraries (e.g., Pyodide)
@@ -234,6 +238,7 @@ with respect to documented and/or tested features.
   (useful, e.g., for stabilized methods and the Stokes problem)
 - Added: `Basis.plot3`, a wrapper to `skfem.visuals.*.plot3`
 - Fixed: Calculation of size in `Basis.__repr__` was slow and incorrect
+- Fixed: Subclasses of `ElementHdiv` did not work together with `FacetBasis`
 
 ### [6.0.0] - 2022-03-15
 
@@ -270,6 +275,8 @@ with respect to documented and/or tested features.
 - Added: `FacetBasis` constructor now has the keyword argument `side` which
   allows changing the side of the facet used to calculate the basis function
   values and gradients
+- Added: `Basis.boundary` method to quickly initialize the corresponding
+  `FacetBasis`
 - Fixed: Improvements to backwards compatibility in `asm`/`assemble` keyword
   arguments
 - Fixed: Save format issue with meshio 5.3.0+
diff --git a/docs/examples/ex37.py b/docs/examples/ex37.py
index 17e80bd2..04d99d28 100644
--- a/docs/examples/ex37.py
+++ b/docs/examples/ex37.py
@@ -10,6 +10,7 @@ m = MeshTet.init_tensor(*(p,) * 3)
 
 e = ElementTetRT1() * ElementTetP0()
 basis = Basis(m, e)
+fbasis = basis.boundary()
 
 
 @BilinearForm
@@ -19,11 +20,14 @@ def bilinf(sigma, u, tau, v, w):
 
 @LinearForm
 def linf(tau, v, w):
-    return - 1. * v
+    return (dot(tau, w.n)
+            * np.sin(np.pi * w.x[1])
+            * np.sin(np.pi * w.x[2])
+            * np.isclose(w.x[0], 1))
 
 
 A = asm(bilinf, basis)
-b = asm(linf, basis)
+b = asm(linf, fbasis)
 x = solve(A, b)
 
 (sigma, rtbasis), (u, ubasis) = basis.split(x)
diff --git a/docs/listofexamples.rst b/docs/listofexamples.rst
index 0be2b726..a2a43578 100644
--- a/docs/listofexamples.rst
+++ b/docs/listofexamples.rst
@@ -142,7 +142,7 @@ Example 37: Mixed Poisson equation
 This example solves the mixed formulation of the Poisson equation
 using the lowest order Raviart-Thomas elements.
 
-.. figure:: https://user-images.githubusercontent.com/973268/93132097-c2862d00-f6dd-11ea-97ad-40aaf2732ad1.png
+.. figure:: https://user-images.githubusercontent.com/973268/182335108-7f29b086-f9b2-42fb-ae9f-dd707d3808f8.png
 
    The piecewise constant solution field.
    The figure was created using `ParaView <https://www.paraview.org/>`__.
diff --git a/skfem/assembly/basis/cell_basis.py b/skfem/assembly/basis/cell_basis.py
index 91c13750..b3920b04 100644
--- a/skfem/assembly/basis/cell_basis.py
+++ b/skfem/assembly/basis/cell_basis.py
@@ -47,8 +47,9 @@ class CellBasis(AbstractBasis):
                  elements: Optional[Any] = None,
                  quadrature: Optional[Tuple[ndarray, ndarray]] = None,
                  dofs: Optional[Dofs] = None):
-        """Combine :class:`~skfem.mesh.Mesh` and :class:`~skfem.element.Element`
-        into a set of precomputed global basis functions.
+        """Combine :class:`~skfem.mesh.Mesh` and
+        :class:`~skfem.element.Element` into a set of precomputed global basis
+        functions.
 
         Parameters
         ----------
@@ -230,6 +231,25 @@ class CellBasis(AbstractBasis):
             elements=self.tind,
         )
 
+    def boundary(self, facets: Optional[Any] = None):
+        """Return corresponding :class:`~skfem.assembly.basis.FacetBasis`.
+
+        Parameters
+        ----------
+        facets
+            Anything that can be passed to ``FacetBasis(..., facets=facets)``.
+
+        """
+        from skfem.assembly.basis.facet_basis import FacetBasis
+        if self.tind is not None:
+            raise NotImplementedError("Boundary of subdomain not supported.")
+        return FacetBasis(
+            self.mesh,
+            self.elem,
+            mapping=self.mapping,
+            facets=facets,
+        )
+
     def project(self, interp, elements=None):
         """Perform :math:`L^2` projection onto the basis.
 
diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py
index feaf76ae..c3231307 100644
--- a/skfem/element/__init__.py
+++ b/skfem/element/__init__.py
@@ -30,7 +30,7 @@ from .element_tet import (ElementTetP0, ElementTetP1, ElementTetP2,  # noqa
                           ElementTetRT1, ElementTetN0, ElementTetMini,
                           ElementTetCR, ElementTetCCR)
 from .element_hex import (ElementHex0, ElementHex1, ElementHex2,  # noqa
-                          ElementHexS2, ElementHex1DG)
+                          ElementHexS2, ElementHex1DG, ElementHexRT1)
 from .element_line import (ElementLineP0, ElementLineP1, ElementLineP2,  # noqa
                            ElementLinePp, ElementLineHermite,
                            ElementLineMini, ElementLineP1DG)
@@ -128,6 +128,7 @@ __all__ = [
     "ElementHexS2",
     "ElementHex1DG",
     "ElementHexDG",
+    "ElementHexRT1",
     "ElementLineP0",
     "ElementLineP1",
     "ElementLineP1DG",
diff --git a/skfem/element/element_h1.py b/skfem/element/element_h1.py
index 52fb043f..4892b21e 100644
--- a/skfem/element/element_h1.py
+++ b/skfem/element/element_h1.py
@@ -21,6 +21,7 @@ class ElementH1(Element):
                 value=np.broadcast_to(phi, (invDF.shape[2], invDF.shape[3])),
                 grad=np.einsum('ijkl,ikl->jkl', invDF, dphi)
             ),)
+        raise NotImplementedError
 
     def lbasis(self, X, i):
         raise NotImplementedError
diff --git a/skfem/element/element_hdiv.py b/skfem/element/element_hdiv.py
index 98e2570f..cc3776dd 100644
--- a/skfem/element/element_hdiv.py
+++ b/skfem/element/element_hdiv.py
@@ -8,16 +8,17 @@ class ElementHdiv(Element):
 
     def orient(self, mapping, i, tind=None):
         """Orientation based on facet-to-triangle indexing."""
-        if tind is not None:
-            # TODO fix
-            raise NotImplementedError("TODO: fix tind support in ElementHdiv")
+        if tind is None:
+            tind = slice(None)
         ix = int(i / self.facet_dofs)
         if ix >= self.refdom.nfacets:
             # interior dofs: no need for orientation
             # TODO support edge DOFs
-            return np.ones(mapping.mesh.t.shape[1], dtype=np.int64)
-        return -1 + 2 * (mapping.mesh.f2t[0, mapping.mesh.t2f[ix]]
-                         == np.arange(mapping.mesh.t.shape[1]))
+            return np.ones(len(np.arange(mapping.mesh.t.shape[1])[tind]),
+                           dtype=np.int64)
+        ori = -1 + 2 * (mapping.mesh.f2t[0, mapping.mesh.t2f[ix]]
+                        == np.arange(mapping.mesh.t.shape[1]))
+        return ori[tind]
 
     def gbasis(self, mapping, X, i, tind=None):
         """Covariant Piola transformation."""
@@ -25,11 +26,19 @@ class ElementHdiv(Element):
         DF = mapping.DF(X, tind)
         detDF = mapping.detDF(X, tind)
         orient = self.orient(mapping, i, tind)
-        return (DiscreteField(
-            value=np.einsum('ijkl,jl,kl->ikl', DF, phi,
-                            1. / np.abs(detDF) * orient[:, None]),
-            div=dphi / (np.abs(detDF) * orient[:, None])
-        ),)
+        if len(X.shape) == 2:
+            return (DiscreteField(
+                value=np.einsum('ijkl,jl,kl->ikl', DF, phi,
+                                1. / np.abs(detDF) * orient[:, None]),
+                div=dphi / (np.abs(detDF) * orient[:, None])
+            ),)
+        elif len(X.shape) == 3:
+            return (DiscreteField(
+                value=np.einsum('ijkl,jkl,kl->ikl', DF, phi,
+                                1. / np.abs(detDF) * orient[:, None]),
+                div=dphi / (np.abs(detDF) * orient[:, None])
+            ),)
+        raise NotImplementedError
 
     def lbasis(self, X, i):
         raise NotImplementedError
diff --git a/skfem/element/element_hex/__init__.py b/skfem/element/element_hex/__init__.py
index b0856099..d887abf1 100644
--- a/skfem/element/element_hex/__init__.py
+++ b/skfem/element/element_hex/__init__.py
@@ -2,3 +2,4 @@ from .element_hex0 import ElementHex0  # noqa
 from .element_hex1 import ElementHex1, ElementHex1DG  # noqa
 from .element_hex_s2 import ElementHexS2  # noqa
 from .element_hex2 import ElementHex2  # noqa
+from .element_hex_rt1 import ElementHexRT1  # noqa
diff --git a/skfem/element/element_hex/element_hex_rt1.py b/skfem/element/element_hex/element_hex_rt1.py
new file mode 100644
index 00000000..43e7c3cc
--- /dev/null
+++ b/skfem/element/element_hex/element_hex_rt1.py
@@ -0,0 +1,45 @@
+import numpy as np
+
+from ..element_hdiv import ElementHdiv
+from ...refdom import RefHex
+
+
+class ElementHexRT1(ElementHdiv):
+    """First-order Raviart-Thomas for cube."""
+
+    facet_dofs = 1
+    maxdeg = 1
+    dofnames = ['u_n']
+    doflocs = np.array([[1.0, 0.5, 0.5],
+                        [0.5, 0.5, 1.0],
+                        [0.5, 1.0, 0.5],
+                        [0.5, 0.0, 0.5],
+                        [0.5, 0.5, 0.0],
+                        [0.0, 0.5, 0.5]])
+    refdom = RefHex
+
+    def lbasis(self, X, i):
+        x, y, z = X
+
+        if i == 0:
+            phi = np.array([x, 0 * y, 0 * z])
+            dphi = 1. + 0. * x
+        elif i == 1:
+            phi = np.array([0 * x, 0 * y, z])
+            dphi = 1. + 0. * x
+        elif i == 2:
+            phi = np.array([0 * x, y, 0 * z])
+            dphi = 1. + 0. * x
+        elif i == 3:
+            phi = np.array([0 * x, y - 1., 0 * z])
+            dphi = 1. + 0. * x
+        elif i == 4:
+            phi = np.array([0 * x, 0 * y, z - 1.])
+            dphi = 1. + 0. * x
+        elif i == 5:
+            phi = np.array([x - 1., 0 * y, 0 * z])
+            dphi = 1. + 0. * x
+        else:
+            self._index_error()
+
+        return phi, dphi
diff --git a/skfem/mapping/mapping.py b/skfem/mapping/mapping.py
index 0629e4dc..f4d8ac6d 100644
--- a/skfem/mapping/mapping.py
+++ b/skfem/mapping/mapping.py
@@ -29,7 +29,8 @@ class Mapping():
     def invF(self,
              x: ndarray,
              tind: Optional[ndarray] = None) -> ndarray:
-        """Perform an inverse mapping from global elements to reference element.
+        """Perform an inverse mapping from global elements to reference
+        element.
 
         Parameters
         ----------
diff --git a/skfem/mapping/mapping_affine.py b/skfem/mapping/mapping_affine.py
index cf711b9c..26379ef1 100644
--- a/skfem/mapping/mapping_affine.py
+++ b/skfem/mapping/mapping_affine.py
@@ -128,8 +128,7 @@ class MappingAffine(Mapping):
             return (np.einsum('ijk,jl', A, X).T + b.T).T
         elif len(X.shape) == 3:
             return (np.einsum('ijk,jkl->ikl', A, X).T + b.T).T
-
-        raise Exception("Wrong dimension of input.")
+        raise NotImplementedError
 
     def invF(self, x, tind=None):
         if tind is None:
@@ -147,7 +146,7 @@ class MappingAffine(Mapping):
         else:
             detDF = self.detA[tind]
 
-        return np.tile(detDF, (X.shape[1], 1)).T
+        return np.tile(detDF, (X.shape[-1], 1)).T
 
     def DF(self, X, tind=None):
         if tind is None:
@@ -155,7 +154,11 @@ class MappingAffine(Mapping):
         else:
             DF = self.A[:, :, tind]
 
-        return np.einsum('ijk,l->ijkl', DF, 1 + 0 * X[0])
+        if len(X.shape) == 2:
+            return np.einsum('ijk,l->ijkl', DF, 1 + np.zeros_like(X[0]))
+        elif len(X.shape) == 3:
+            return np.einsum('ijk,kl->ijkl', DF, 1 + np.zeros_like(X[0]))
+        raise NotImplementedError
 
     def invDF(self, X, tind=None):
         if tind is None:
diff --git a/skfem/mapping/mapping_isoparametric.py b/skfem/mapping/mapping_isoparametric.py
index 1528a3c0..32edf095 100644
--- a/skfem/mapping/mapping_isoparametric.py
+++ b/skfem/mapping/mapping_isoparametric.py
@@ -160,8 +160,7 @@ class MappingIsoparametric(Mapping):
         if J is None:
             J = [[self.J(i, j, X, tind=tind) for j in range(self.dim)]
                  for i in range(self.dim)]
-
-        return np.array(J if tind is None else J[:, :, tind])
+        return np.array(J)
 
     def detDF(self, X, tind=None, J=None):
         if J is None:
diff --git a/skfem/mesh/mesh_quad_1.py b/skfem/mesh/mesh_quad_1.py
index 920e4b60..5f945250 100644
--- a/skfem/mesh/mesh_quad_1.py
+++ b/skfem/mesh/mesh_quad_1.py
@@ -184,8 +184,7 @@ class MeshQuad1(Mesh2D):
             for k in self.boundaries:
                 slots = enumerate(mesh.facets.T)
                 boundaries[k] = np.array([
-                    next(dropwhile(lambda slot: not(np.array_equal(f,
-                                                                   slot[1])),
+                    next(dropwhile(lambda s: (not np.array_equal(f, s[1])),
                                    slots))[0]
                     for f in self.facets.T[np.sort(self.boundaries[k])]])
 


Test patch:
diff --git a/tests/test_assembly.py b/tests/test_assembly.py
index 7f7a419d..6c31b795 100644
--- a/tests/test_assembly.py
+++ b/tests/test_assembly.py
@@ -13,7 +13,9 @@ from skfem.element import (ElementQuad1, ElementQuadS2, ElementHex1,
                            ElementTriMorley, ElementVectorH1, ElementQuadP,
                            ElementHex2, ElementTriArgyris, ElementTriP2,
                            ElementTriDG, ElementQuadDG, ElementHexDG,
-                           ElementTetDG, ElementTriHermite, ElementVector)
+                           ElementTetDG, ElementTriHermite, ElementVector,
+                           ElementTriRT1, ElementTriRT2, ElementTriBDM1,
+                           ElementQuadRT1, ElementTetRT1, ElementHexRT1)
 from skfem.mesh import (MeshQuad, MeshHex, MeshTet, MeshTri, MeshQuad2,
                         MeshTri2, MeshTet2, MeshHex2, MeshTri1DG, MeshQuad1DG,
                         MeshHex1DG)
@@ -588,3 +590,60 @@ def test_matrix_element_projection(m, e):
         return ddot(C(sym_grad(x)), v)
 
     y = projection(proj, basis1, basis0)
+
+
+@pytest.mark.parametrize(
+    "basis",
+    [
+        Basis(MeshTri().refined(6), ElementTriRT1()),
+        Basis(MeshTri().refined(4), ElementTriRT2()),
+        Basis(MeshTri().refined(5), ElementTriBDM1()),
+        Basis(MeshQuad().refined(4), ElementQuadRT1()),
+    ]
+)
+def test_hdiv_boundary_integration(basis):
+
+    y = basis.project(lambda x: np.array([x[0], 0 * x[1]]))
+    fbasis = basis.boundary('right')
+
+    @Functional
+    def test1(w):
+        return dot(w['y'], w.n)
+
+    @Functional
+    def test2(w):
+        return w['y'].div
+
+    assert_almost_equal(test1.assemble(fbasis, y=y), 1, decimal=2)
+
+    if not isinstance(basis.elem, ElementTriRT1):
+        assert_almost_equal(test2.assemble(fbasis, y=y), 1, decimal=5)
+
+
+@pytest.mark.parametrize(
+    "basis",
+    [
+        Basis(MeshTet().refined(4).with_boundaries({
+            'right': lambda x: x[0] == 1
+        }), ElementTetRT1()),
+        Basis(MeshHex().refined(4).with_boundaries({
+            'right': lambda x: x[0] == 1
+        }), ElementHexRT1()),
+    ]
+)
+def test_hdiv_boundary_integration_3d(basis):
+
+    y = basis.project(lambda x: np.array([x[0], 0 * x[1], 0 * x[2]]))
+    fbasis = basis.boundary('right')
+
+    @Functional
+    def test1(w):
+        return dot(w['y'], w.n)
+
+    @Functional
+    def test2(w):
+        return w['y'].div
+
+    assert_almost_equal(test1.assemble(fbasis, y=y), 1, decimal=2)
+    if not isinstance(basis.elem, ElementTetRT1):
+        assert_almost_equal(test2.assemble(fbasis, y=y), 1, decimal=5)
diff --git a/tests/test_convergence_hdiv.py b/tests/test_convergence_hdiv.py
index 3a27484a..1bc270c2 100644
--- a/tests/test_convergence_hdiv.py
+++ b/tests/test_convergence_hdiv.py
@@ -6,8 +6,8 @@ from skfem import BilinearForm, CellBasis, LinearForm, asm, solve
 from skfem.element import (ElementTetP0, ElementTetRT0, ElementTriP0,
                            ElementTriRT0, ElementTriBDM1, ElementDG,
                            ElementQuadRT0, ElementQuad0, ElementTriRT2,
-                           ElementTriP1)
-from skfem.mesh import MeshTet, MeshTri, MeshQuad
+                           ElementTriP1, ElementHexRT1, ElementHex0)
+from skfem.mesh import MeshTet, MeshTri, MeshQuad, MeshHex
 
 
 class ConvergenceRaviartThomas(unittest.TestCase):
@@ -198,3 +198,20 @@ class ConvergenceRaviartThomas3D(ConvergenceRaviartThomas):
 
     def setUp(self):
         self.mesh = MeshTet().refined(1)
+
+
+class ConvergenceRaviartThomas3DHex(ConvergenceRaviartThomas3D):
+
+    rateL2 = 1.0
+    L2bound = 0.1
+    rateHdiv = 1.0
+    Hdivbound = 0.31
+
+    def create_basis(self, m):
+        e = ElementHexRT1()
+        e0 = ElementHex0()
+        return (CellBasis(m, e, intorder=2),
+                CellBasis(m, e0, intorder=2))
+
+    def setUp(self):
+        self.mesh = MeshHex().refined(1)
diff --git a/tests/test_examples.py b/tests/test_examples.py
index a14797a9..b14df486 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -302,7 +302,7 @@ class TestEx37(TestCase):
 
     def runTest(self):
         from docs.examples.ex37 import u
-        self.assertAlmostEqual(np.max(u), 0.05594193697362236)
+        self.assertAlmostEqual(np.max(u), 0.8619044895314989)
 
 
 class TestEx38(TestCase):
diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index 74ad5571..e37afa37 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -4,7 +4,7 @@ import numpy as np
 from skfem.mesh import MeshHex, MeshQuad, MeshTri
 from skfem.element import ElementHex1, ElementQuad1, ElementHex2
 from skfem.assembly import FacetBasis
-from skfem.mapping.mapping_mortar import MappingMortar
+from skfem.mapping import MappingMortar
 
 
 class TestIsoparamNormals(unittest.TestCase):