From 5cde70315779fc33ea6b16bfbb983d8acc8ad061 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 12 Sep 2020 14:32:03 +0300
Subject: [PATCH 01/13] Start implementing ElementHex2

---
 skfem/element/element_hex/element_hex2.py | 131 ++++++++++++++++++++++
 1 file changed, 131 insertions(+)
 create mode 100644 skfem/element/element_hex/element_hex2.py

diff --git a/skfem/element/element_hex/element_hex2.py b/skfem/element/element_hex/element_hex2.py
new file mode 100644
index 000000000..09539355f
--- /dev/null
+++ b/skfem/element/element_hex/element_hex2.py
@@ -0,0 +1,131 @@
+import numpy as np
+from numpy.linalg import inv
+
+from ..element_h1 import ElementH1
+from ...mesh.mesh3d import MeshHex
+
+
+class ElementHexS2(ElementH1):
+    nodal_dofs = 1
+    facet_dofs = 1
+    edge_dofs = 1
+    interior_dofs = 1
+    dim = 3
+    maxdeg = 4
+    dofnames = ['u', 'u', 'u', 'u']
+    doflocs = np.array([[1., 1., 1.],
+                        [1., 1., 0.],
+                        [1., 0., 1.],
+                        [0., 1., 1.],
+                        [1., 0., 0.],
+                        [0., 1., 0.],
+                        [0., 0., 1.],
+                        [0., 0., 0.],
+                        [1., .5, .5],
+                        [.5, .5, 1.],
+                        [.5, 1., .5],
+                        [.5, .0, .5],
+                        [.5, .5, 0.],
+                        [0., .5, .5],
+                        [1., 1., .5],
+                        [1., .5, 1.],
+                        [.5, 1., 1.],
+                        [1., .5, 0.],
+                        [.5, 1., 0.],
+                        [1., 0., .5],
+                        [.5, 0., 1.],
+                        [0., 1., .5],
+                        [0., .5, 1.],
+                        [.5, 0., 0.],
+                        [0., .5, 0.],
+                        [0., 0., .5],
+                        [.5, .5, .5]])
+    mesh_type = MeshHex
+
+    def __init__(self):
+
+        X = self.doflocs.T
+        V = np.array([
+            1. + 0. * X[0],
+            X[0],
+            X[1],
+            X[2],
+            X[0] ** 2,
+            X[1] ** 2,
+            X[2] ** 2,
+            X[0] * X[1],
+            X[1] * X[2],
+            X[0] * X[2],
+            X[0] ** 2 * X[1],
+            X[1] ** 2 * X[2],
+            X[0] ** 2 * X[2],
+            X[0] * X[1] ** 2,
+            X[1] * X[2] ** 2,
+            X[0] * X[2] ** 2,
+            X[0] ** 2 * X[1] ** 2,
+            X[1] ** 2 * X[2] ** 2,
+            X[0] ** 2 * X[2] ** 2,
+            X[0] * X[1] * X[2],
+            X[0] ** 2 * X[1] * X[2],
+            X[0] * X[1] ** 2 * X[2],
+            X[0] * X[1] * X[2] ** 2,
+            X[0] ** 2 * X[1] ** 2 * X[2],
+            X[0] * X[1] ** 2 * X[2] ** 2,
+            X[0] ** 2 * X[1] * X[2] ** 2,
+            X[0] ** 2 * X[1] ** 2 * X[2] ** 2,
+        ])
+        self.invV = inv(V)
+
+    def lbasis(self, X, i):
+        x, y, z = 2 * X - 1
+
+        if i < 8:
+            s = [
+                (1, 1, 1),
+                (1, 1, -1),
+                (1, -1, 1),
+                (-1, 1, 1),
+                (1, -1, -1),
+                (-1, 1, -1),
+                (-1, -1, 1),
+                (-1, -1, -1),
+            ][i]
+            x *= s[0]
+            y *= s[1]
+            z *= s[2]
+            phi = (1 + x) * (1 + y) * (1 + z) * (x + y + z - 2) / 8
+            dphi = np.array([s[0] * (1 + y) * (1 + z) * (x + y + z - 2)
+                             + s[0] * (1 + x) * (1 + y) * (1 + z),
+                             s[1] * (1 + x) * (1 + z) * (x + y + z - 2)
+                             + s[1] * (1 + x) * (1 + y) * (1 + z),
+                             s[2] * (1 + x) * (1 + y) * (x + y + z - 2)
+                             + s[2] * (1 + x) * (1 + y) * (1 + z)]) / 8
+        elif i < 20:
+            s = [
+                (1, 1, -z),
+                (1, -y, 1),
+                (-x, 1, 1),
+                (1, -y, -1),
+                (-x, 1, -1),
+                (1, -1, -z),
+                (-x, -1, 1),
+                (-1, 1, -z),
+                (-1, -y, 1),
+                (-x, -1, -1),
+                (-1, -y, -1),
+                (-1, -1, -z),
+            ][i - 8]
+            x *= s[0]
+            y *= s[1]
+            z *= s[2]
+            phi = (1 + x) * (1 + y) * (1 + z) / 4
+            dx = 2.0 * s[0] if isinstance(s[0], np.ndarray) else s[0]
+            dy = 2.0 * s[1] if isinstance(s[1], np.ndarray) else s[1]
+            dz = 2.0 * s[2] if isinstance(s[2], np.ndarray) else s[2]
+            dphi = np.array([dx * (1 + y) * (1 + z),
+                             dy * (1 + x) * (1 + z),
+                             dz * (1 + x) * (1 + y)]) / 4
+        else:
+            self._index_error()
+
+        return phi, 2 * dphi

From f1881a91a40bf9da9fa24a1e3045c0b7a440e92f Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 00:28:22 +0300
Subject: [PATCH 02/13] Implement ElementHex2 local basis via Vandermonde
 matrix

---
 skfem/element/__init__.py                   |   3 +-
 skfem/element/element_hex/__init__.py       |   1 +
 skfem/element/element_hex/element_hex2.py   | 165 +++++++++++++-------
 skfem/element/element_hex/element_hex_s2.py |   6 +-
 tests/test_convergence.py                   |  28 +++-
 5 files changed, 142 insertions(+), 61 deletions(-)

diff --git a/skfem/element/__init__.py b/skfem/element/__init__.py
index 719861331..0f809b9b4 100644
--- a/skfem/element/__init__.py
+++ b/skfem/element/__init__.py
@@ -42,7 +42,7 @@
     ElementQuadBFS
 from .element_tet import ElementTetP0, ElementTetP1, ElementTetP2,\
     ElementTetRT0, ElementTetN0, ElementTetMini
-from .element_hex import ElementHex1, ElementHexS2  # noqa
+from .element_hex import ElementHex1, ElementHex2, ElementHexS2  # noqa
 from .element_line import ElementLineP1, ElementLineP2,\
     ElementLinePp, ElementLineHermite, ElementLineMini  # noqa
 from .element_composite import ElementComposite  # noqa
@@ -77,6 +77,7 @@
     "ElementTetN0",
     "ElementTetMini",
     "ElementHex1",
+    "ElementHex2",
     "ElementHexS2",
     "ElementLineP1",
     "ElementLineP2",
diff --git a/skfem/element/element_hex/__init__.py b/skfem/element/element_hex/__init__.py
index 33f38ec14..b263c59c9 100644
--- a/skfem/element/element_hex/__init__.py
+++ b/skfem/element/element_hex/__init__.py
@@ -1,2 +1,3 @@
 from .element_hex1 import ElementHex1  # noqa
 from .element_hex_s2 import ElementHexS2  # noqa
+from .element_hex2 import ElementHex2  # noqa
diff --git a/skfem/element/element_hex/element_hex2.py b/skfem/element/element_hex/element_hex2.py
index 09539355f..39edfe414 100644
--- a/skfem/element/element_hex/element_hex2.py
+++ b/skfem/element/element_hex/element_hex2.py
@@ -5,13 +5,14 @@
 from ...mesh.mesh3d import MeshHex
 
 
-class ElementHexS2(ElementH1):
+class ElementHex2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     edge_dofs = 1
     interior_dofs = 1
     dim = 3
-    maxdeg = 4
+    maxdeg = 6
     dofnames = ['u', 'u', 'u', 'u']
     doflocs = np.array([[1., 1., 1.],
                         [1., 1., 0.],
@@ -43,9 +44,13 @@ class ElementHexS2(ElementH1):
     mesh_type = MeshHex
 
     def __init__(self):
-
         X = self.doflocs.T
-        V = np.array([
+        V = self._power_basis(X)
+        self.invV = inv(V).T
+        self.invV[np.abs(self.invV) < 1e-10] = 0.
+
+    def _power_basis(self, X):
+        return np.array([
             1. + 0. * X[0],
             X[0],
             X[1],
@@ -74,58 +79,106 @@ def __init__(self):
             X[0] ** 2 * X[1] * X[2] ** 2,
             X[0] ** 2 * X[1] ** 2 * X[2] ** 2,
         ])
-        self.invV = inv(V)
 
-    def lbasis(self, X, i):
-        x, y, z = 2 * X - 1
+    def _power_basis_dx(self, X):
+        return np.array([
+            0. * X[0],
+            1. + 0. * X[0],
+            0. * X[0],
+            0. * X[0],
+            2. * X[0],
+            0. * X[0],
+            0. * X[0],
+            X[1],
+            0. * X[0],
+            X[2],
+            2. * X[0] * X[1],
+            0. * X[0],
+            2. * X[0] * X[2],
+            X[1] ** 2,
+            0. * X[0],
+            X[2] ** 2,
+            2. * X[0] * X[1] ** 2,
+            0. * X[0],
+            2. * X[0] * X[2] ** 2,
+            X[1] * X[2],
+            2. * X[0] * X[1] * X[2],
+            X[1] ** 2 * X[2],
+            X[1] * X[2] ** 2,
+            2. * X[0] * X[1] ** 2 * X[2],
+            X[1] ** 2 * X[2] ** 2,
+            2. * X[0] * X[1] * X[2] ** 2,
+            2. * X[0] * X[1] ** 2 * X[2] ** 2,
+        ])
+
+    def _power_basis_dy(self, X):
+        return np.array([
+            0. * X[0],
+            0. * X[0],
+            1. + 0. * X[0],
+            0. * X[0],
+            0. * X[0],
+            2. * X[1],
+            0. * X[0],
+            X[0],
+            X[2],
+            0. * X[0],
+            X[0] ** 2,
+            2. * X[1] * X[2],
+            0. * X[0],
+            X[0] * 2. * X[1],
+            X[2] ** 2,
+            0. * X[0],
+            X[0] ** 2 * 2. * X[1],
+            2. * X[1] * X[2] ** 2,
+            0. * X[0],
+            X[0] * X[2],
+            X[0] ** 2 * X[2],
+            X[0] * 2. *  X[1] * X[2],
+            X[0] * X[2] ** 2,
+            X[0] ** 2 * 2. * X[1] * X[2],
+            X[0] * 2. * X[1] * X[2] ** 2,
+            X[0] ** 2 * X[2] ** 2,
+            X[0] ** 2 * 2. * X[1] * X[2] ** 2,
+        ])
 
-        if i < 8:
-            s = [
-                (1, 1, 1),
-                (1, 1, -1),
-                (1, -1, 1),
-                (-1, 1, 1),
-                (1, -1, -1),
-                (-1, 1, -1),
-                (-1, -1, 1),
-                (-1, -1, -1),
-            ][i]
-            x *= s[0]
-            y *= s[1]
-            z *= s[2]
-            phi = (1 + x) * (1 + y) * (1 + z) * (x + y + z - 2) / 8
-            dphi = np.array([s[0] * (1 + y) * (1 + z) * (x + y + z - 2)
-                             + s[0] * (1 + x) * (1 + y) * (1 + z),
-                             s[1] * (1 + x) * (1 + z) * (x + y + z - 2)
-                             + s[1] * (1 + x) * (1 + y) * (1 + z),
-                             s[2] * (1 + x) * (1 + y) * (x + y + z - 2)
-                             + s[2] * (1 + x) * (1 + y) * (1 + z)]) / 8
-        elif i < 20:
-            s = [
-                (1, 1, -z),
-                (1, -y, 1),
-                (-x, 1, 1),
-                (1, -y, -1),
-                (-x, 1, -1),
-                (1, -1, -z),
-                (-x, -1, 1),
-                (-1, 1, -z),
-                (-1, -y, 1),
-                (-x, -1, -1),
-                (-1, -y, -1),
-                (-1, -1, -z),
-            ][i - 8]
-            x *= s[0]
-            y *= s[1]
-            z *= s[2]
-            phi = (1 + x) * (1 + y) * (1 + z) / 4
-            dx = 2.0 * s[0] if isinstance(s[0], np.ndarray) else s[0]
-            dy = 2.0 * s[1] if isinstance(s[1], np.ndarray) else s[1]
-            dz = 2.0 * s[2] if isinstance(s[2], np.ndarray) else s[2]
-            dphi = np.array([dx * (1 + y) * (1 + z),
-                             dy * (1 + x) * (1 + z),
-                             dz * (1 + x) * (1 + y)]) / 4
-        else:
-            self._index_error()
+    def _power_basis_dz(self, X):
+        return np.array([
+            0. * X[0],
+            0. * X[0],
+            0. * X[0],
+            1. + 0. * X[0],
+            0. * X[0],
+            0. * X[0],
+            2. * X[2],
+            0. * X[0],
+            X[1],
+            X[0],
+            0. * X[0],
+            X[1] ** 2,
+            X[0] ** 2,
+            0. * X[0],
+            X[1] * 2. * X[2],
+            X[0] * 2. * X[2],
+            0. * X[0],
+            X[1] ** 2 * 2. * X[2],
+            X[0] ** 2 * 2. * X[2],
+            X[0] * X[1],
+            X[0] ** 2 * X[1],
+            X[0] * X[1] ** 2,
+            X[0] * X[1] * 2. * X[2],
+            X[0] ** 2 * X[1] ** 2,
+            X[0] * X[1] ** 2 * 2. * X[2],
+            X[0] ** 2 * X[1] * 2. * X[2],
+            X[0] ** 2 * X[1] ** 2 * 2. * X[2],
+        ])
 
-        return phi, 2 * dphi
+    def lbasis(self, X, i):
+        return (
+            np.sum(self.invV[i][:, None] * self._power_basis(X), axis=0),
+            np.array([
+                np.sum(self.invV[i][:, None] * self._power_basis_dx(X), axis=0),
+                np.sum(self.invV[i][:, None] * self._power_basis_dy(X), axis=0),
+                np.sum(self.invV[i][:, None] * self._power_basis_dz(X), axis=0),
+            ])
+        )
diff --git a/skfem/element/element_hex/element_hex_s2.py b/skfem/element/element_hex/element_hex_s2.py
index 06972a5bc..160daaad9 100644
--- a/skfem/element/element_hex/element_hex_s2.py
+++ b/skfem/element/element_hex/element_hex_s2.py
@@ -75,9 +75,9 @@ def lbasis(self, X, i):
             y *= s[1]
             z *= s[2]
             phi = (1 + x) * (1 + y) * (1 + z) / 4
-            dx = 2.0 * s[0] if isinstance(s[0], np.ndarray) else s[0]
-            dy = 2.0 * s[1] if isinstance(s[1], np.ndarray) else s[1]
-            dz = 2.0 * s[2] if isinstance(s[2], np.ndarray) else s[2]
+            dx = 2. * s[0] if isinstance(s[0], np.ndarray) else s[0]
+            dy = 2. * s[1] if isinstance(s[1], np.ndarray) else s[1]
+            dz = 2. * s[2] if isinstance(s[2], np.ndarray) else s[2]
             dphi = np.array([dx * (1 + y) * (1 + z),
                              dy * (1 + x) * (1 + z),
                              dz * (1 + x) * (1 + y)]) / 4
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index 1d15cbab9..a2693460f 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -4,7 +4,7 @@
 
 from skfem import BilinearForm, LinearForm, asm, solve, condense
 from skfem.models.poisson import laplace
-from skfem.element import (ElementHex1, ElementHexS2, ElementLineP1,
+from skfem.element import (ElementHex1, ElementHex2, ElementHexS2, ElementLineP1,
                            ElementLineP2, ElementLineMini, 
                            ElementQuad1, ElementQuad2,
                            ElementQuadS2, ElementTetMini, ElementTetP1,
@@ -15,6 +15,7 @@
 
 
 class ConvergenceQ1(unittest.TestCase):
+
     rateL2 = 2.0
     rateH1 = 1.0
     eps = 0.1
@@ -145,6 +146,7 @@ def setUp(self):
 
 
 class ConvergenceQ2(ConvergenceQ1):
+
     rateL2 = 3.0
     rateH1 = 2.0
 
@@ -176,6 +178,7 @@ def setUp(self):
 
 
 class ConvergenceTriP2(ConvergenceTriP1):
+
     rateL2 = 3.0
     rateH1 = 2.0
 
@@ -192,6 +195,7 @@ def create_basis(self, m):
 
 
 class ConvergenceHex1(ConvergenceQ1):
+
     rateL2 = 2.0
     rateH1 = 1.0
     eps = 0.11
@@ -206,6 +210,7 @@ def setUp(self):
 
 
 class ConvergenceHexS2(ConvergenceQ1):
+
     rateL2 = 3.05
     rateH1 = 2.21
     eps = 0.02
@@ -219,7 +224,23 @@ def setUp(self):
         self.mesh.refine(1)
 
 
+class ConvergenceHex2(ConvergenceQ1):
+
+    rateL2 = 1.67
+    rateH1 = 2.21
+    eps = 0.02
+
+    def create_basis(self, m):
+        e = ElementHex2()
+        return InteriorBasis(m, e)
+
+    def setUp(self):
+        self.mesh = MeshHex()
+        self.mesh.refine(1)
+
+
 class ConvergenceTetP1(ConvergenceQ1):
+
     rateL2 = 2.0
     rateH1 = 1.0
     eps = 0.13
@@ -234,6 +255,7 @@ def setUp(self):
 
 
 class ConvergenceTetP2(ConvergenceTetP1):
+
     rateL2 = 3.23
     rateH1 = 1.94
     eps = 0.01
@@ -266,6 +288,7 @@ def setUp(self):
 
 
 class ConvergenceLineP2(ConvergenceQ1):
+
     rateL2 = 3.0
     rateH1 = 2.0
 
@@ -390,18 +413,21 @@ def uz(x):
 
 
 class FacetConvergenceHex1(FacetConvergenceTetP2):
+
     case = (MeshHex, ElementHex1)
     limits = (0.9, 1.1)
     preref = 2
 
 
 class FacetConvergenceHexS2(FacetConvergenceTetP2):
+
     case = (MeshHex, ElementHexS2)
     limits = (1.9, 2.2)
     preref = 1
 
 
 class FacetConvergenceTetP1(FacetConvergenceTetP2):
+
     case = (MeshTet, ElementTetP1)
     limits = (0.9, 1.1)
     preref = 2

From bb29a1656a7e52d6d4f217fa2c623afe793ceb1e Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:36:01 +0300
Subject: [PATCH 03/13] More tests and fixes to ElementHex2

---
 skfem/element/element_hex/element_hex2.py | 12 ++++---
 tests/test_assembly.py                    | 43 ++++++++++++++++++++++-
 tests/test_basis.py                       |  7 +++-
 tests/test_convergence.py                 | 10 +-----
 tests/test_elements.py                    | 16 +++++----
 5 files changed, 66 insertions(+), 22 deletions(-)

diff --git a/skfem/element/element_hex/element_hex2.py b/skfem/element/element_hex/element_hex2.py
index 39edfe414..89b83e48b 100644
--- a/skfem/element/element_hex/element_hex2.py
+++ b/skfem/element/element_hex/element_hex2.py
@@ -46,7 +46,7 @@ class ElementHex2(ElementH1):
     def __init__(self):
         X = self.doflocs.T
         V = self._power_basis(X)
-        self.invV = inv(V).T
+        self.invV = inv(V)
         self.invV[np.abs(self.invV) < 1e-10] = 0.
 
     def _power_basis(self, X):
@@ -174,11 +174,13 @@ def _power_basis_dz(self, X):
         ])
 
     def lbasis(self, X, i):
+        if i >= 27 or i < 0:
+            self._index_error()
         return (
-            np.sum(self.invV[i][:, None] * self._power_basis(X), axis=0),
+            np.einsum('i...,i...', self.invV[i], self._power_basis(X)),
             np.array([
-                np.sum(self.invV[i][:, None] * self._power_basis_dx(X), axis=0),
-                np.sum(self.invV[i][:, None] * self._power_basis_dy(X), axis=0),
-                np.sum(self.invV[i][:, None] * self._power_basis_dz(X), axis=0),
+                np.einsum('i...,i...', self.invV[i], self._power_basis_dx(X)),
+                np.einsum('i...,i...', self.invV[i], self._power_basis_dy(X)),
+                np.einsum('i...,i...', self.invV[i], self._power_basis_dz(X)),
             ])
         )
diff --git a/tests/test_assembly.py b/tests/test_assembly.py
index 1c98528fe..aca53bdd5 100644
--- a/tests/test_assembly.py
+++ b/tests/test_assembly.py
@@ -6,13 +6,15 @@
 from skfem.element import (ElementQuad1, ElementQuadS2, ElementHex1,
                            ElementHexS2, ElementTetP0, ElementTetP1,
                            ElementTetP2, ElementTriP1, ElementQuad2,
-                           ElementTriMorley, ElementVectorH1, ElementQuadP)
+                           ElementTriMorley, ElementVectorH1, ElementQuadP,
+                           ElementHex2)
 from skfem.mesh import MeshQuad, MeshHex, MeshTet, MeshTri
 from skfem.assembly import FacetBasis, InteriorBasis
 from skfem.utils import project
 
 
 class IntegrateOneOverBoundaryQ1(unittest.TestCase):
+
     elem = ElementQuad1()
 
     def createBasis(self):
@@ -43,6 +45,7 @@ def gv(v, w):
 
 
 class IntegrateOneOverBoundaryS2(IntegrateOneOverBoundaryQ1):
+
     elem = ElementQuadS2()
 
 
@@ -64,6 +67,15 @@ def createBasis(self):
         self.boundary_area = 6.000
 
 
+class IntegrateOneOverBoundaryHex2(IntegrateOneOverBoundaryQ1):
+
+    def createBasis(self):
+        m = MeshHex()
+        m.refine(3)
+        self.fbasis = FacetBasis(m, ElementHex2())
+        self.boundary_area = 6.000
+
+
 class IntegrateFuncOverBoundary(unittest.TestCase):
 
     def runTest(self):
@@ -89,6 +101,7 @@ def uv(u, v, w):
 
 
 class IntegrateFuncOverBoundaryPart(unittest.TestCase):
+
     case = (MeshHex, ElementHex1)
 
     def runTest(self):
@@ -110,22 +123,32 @@ def uv(u, v, w):
 
 
 class IntegrateFuncOverBoundaryPartHexS2(IntegrateFuncOverBoundaryPart):
+
     case = (MeshHex, ElementHexS2)
 
 
+class IntegrateFuncOverBoundaryPartHex2(IntegrateFuncOverBoundaryPart):
+
+    case = (MeshHex, ElementHex2)
+
+
 class IntegrateFuncOverBoundaryPartTetP1(IntegrateFuncOverBoundaryPart):
+
     case = (MeshTet, ElementTetP1)
 
 
 class IntegrateFuncOverBoundaryPartTetP2(IntegrateFuncOverBoundaryPart):
+
     case = (MeshTet, ElementTetP2)
 
 
 class IntegrateFuncOverBoundaryPartTetP0(IntegrateFuncOverBoundaryPart):
+
     case = (MeshTet, ElementTetP0)
 
 
 class BasisInterpolator(unittest.TestCase):
+
     case = (MeshTri, ElementTriP1)
 
     def initOnes(self, basis):
@@ -146,18 +169,22 @@ def runTest(self):
 
 
 class BasisInterpolatorTriP2(BasisInterpolator):
+
     case = (MeshQuad, ElementQuad1)
 
 
 class BasisInterpolatorQuad1(BasisInterpolator):
+
     case = (MeshQuad, ElementQuad1)
 
 
 class BasisInterpolatorQuad2(BasisInterpolator):
+
     case = (MeshQuad, ElementQuad2)
 
 
 class BasisInterpolatorQuadS2(BasisInterpolator):
+
     case = (MeshQuad, ElementQuadS2)
 
 
@@ -182,6 +209,7 @@ def ones(v, w):
 
 
 class NormalVectorTestTri(unittest.TestCase):
+
     case = (MeshTri(), ElementTriP1())
     test_integrate_volume = True
     intorder = None
@@ -221,34 +249,47 @@ def linf(v, w):
 
 
 class NormalVectorTestTet(NormalVectorTestTri):
+
     case = (MeshTet(), ElementTetP1())
 
 
 class NormalVectorTestTetP2(NormalVectorTestTri):
+
     case = (MeshTet(), ElementTetP2())
     test_integrate_volume = False
 
 
 class NormalVectorTestQuad(NormalVectorTestTri):
+
     case = (MeshQuad(), ElementQuad1())
 
 
 class NormalVectorTestQuadP(NormalVectorTestTri):
+
     case = (MeshQuad(), ElementQuadP(3))
     test_integrate_volume = False
 
 
 class NormalVectorTestHex(NormalVectorTestTri):
+
     case = (MeshHex(), ElementHex1())
     intorder = 3
 
 
 class NormalVectorTestHexS2(NormalVectorTestTri):
+
     case = (MeshHex(), ElementHexS2())
     intorder = 3
     test_integrate_volume = False
 
 
+class NormalVectorTestHex2(NormalVectorTestTri):
+
+    case = (MeshHex(), ElementHex2())
+    intorder = 3
+    test_integrate_volume = False
+
+
 class EvaluateFunctional(unittest.TestCase):
 
     def runTest(self):
diff --git a/tests/test_basis.py b/tests/test_basis.py
index 1680cbdcc..31011519e 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -7,7 +7,7 @@
 from skfem.mesh import MeshTri, MeshTet, MeshHex
 from skfem.assembly import InteriorBasis, FacetBasis, Dofs
 from skfem.element import (ElementVectorH1, ElementTriP2, ElementTriP1,
-                           ElementTetP2, ElementHexS2)
+                           ElementTetP2, ElementHexS2, ElementHex2)
 
 
 class TestCompositeSplitting(TestCase):
@@ -117,3 +117,8 @@ class TestFacetExpansionHexS2(TestFacetExpansion):
 
     mesh_type = MeshHex
     elem_type = ElementHexS2
+
+
+class TestFacetExpansionHex2(TestFacetExpansionHexS2):
+
+    elem_type = ElementHex2
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index a2693460f..67d4d9cf4 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -224,20 +224,12 @@ def setUp(self):
         self.mesh.refine(1)
 
 
-class ConvergenceHex2(ConvergenceQ1):
-
-    rateL2 = 1.67
-    rateH1 = 2.21
-    eps = 0.02
+class ConvergenceHex2(ConvergenceHexS2):
 
     def create_basis(self, m):
         e = ElementHex2()
         return InteriorBasis(m, e)
 
-    def setUp(self):
-        self.mesh = MeshHex()
-        self.mesh.refine(1)
-
 
 class ConvergenceTetP1(ConvergenceQ1):
 
diff --git a/tests/test_elements.py b/tests/test_elements.py
index 2ea984716..6db14b845 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -1,7 +1,7 @@
 from unittest import TestCase, main
 
 import numpy as np
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_allclose, assert_array_equal
 
 from skfem.element import (ElementHex1, ElementHexS2, ElementLineP1, 
                            ElementLineP2, ElementLinePp, ElementLineMini,
@@ -9,7 +9,8 @@
                            ElementQuadP, ElementQuadS2, ElementTetMini, 
                            ElementTetP0, ElementTetP1, ElementTetP2,
                            ElementTriMini, ElementTriP0, ElementTriP1,
-                           ElementTriP2, ElementTriRT0, ElementVectorH1)
+                           ElementTriP2, ElementTriRT0, ElementVectorH1,
+                           ElementHex2)
 from skfem.mesh import MeshHex, MeshLine, MeshQuad, MeshTet, MeshTri
 from skfem.assembly.basis import InteriorBasis
 from skfem.mapping import MappingAffine
@@ -40,6 +41,7 @@ class TestNodality(TestCase):
         ElementTetMini(),
         ElementHex1(),
         ElementHexS2(),
+        ElementHex2(),
     ]
 
     def runTest(self):
@@ -62,8 +64,8 @@ def runTest(self):
             ixs = np.nonzero(~ix)[0]
             Ih = Ih[ixs].T[ixs].T
 
-            assert_array_equal(Ih, np.eye(N - Nnan),
-                               err_msg="{}".format(type(e)))
+            assert_allclose(Ih, np.eye(N - Nnan), atol=1e-13,
+                            err_msg="{}".format(type(e)))
 
 
 class TestNodalityTriRT0(TestCase):
@@ -83,8 +85,8 @@ def runTest(self):
             n = np.array([1., np.sqrt(2), 1.])
             Ih[itr] = A * n
 
-        assert_array_equal(Ih, np.eye(N),
-                           err_msg="{}".format(type(e)))
+        assert_allclose(Ih, np.eye(N),
+                        err_msg="{}".format(type(e)))
 
 
 class TestComposite(TestCase):
@@ -165,6 +167,7 @@ class TestDerivatives(TestCase):
         ElementTetMini(),
         ElementHex1(),
         ElementHexS2(),
+        ElementHex2(),
     ]
 
     def runTest(self):
@@ -218,6 +221,7 @@ class TestPartitionofUnity(TestCase):
         ElementTetP2(),
         ElementHex1(),
         ElementHexS2(),
+        ElementHex2(),
     ]
 
     def runTest(self):

From 9cd69df6c666daf1b241124e3630c1b0e319401a Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:39:04 +0300
Subject: [PATCH 04/13] Make whitespace consistent

---
 skfem/element/element_line/element_line_hermite.py | 1 +
 skfem/element/element_line/element_line_mini.py    | 1 +
 skfem/element/element_line/element_line_p1.py      | 1 +
 skfem/element/element_line/element_line_p2.py      | 1 +
 skfem/element/element_line/element_line_pp.py      | 1 +
 skfem/element/element_quad/element_quad0.py        | 1 +
 skfem/element/element_quad/element_quad1.py        | 1 +
 skfem/element/element_quad/element_quad2.py        | 1 +
 skfem/element/element_quad/element_quad_bfs.py     | 1 +
 skfem/element/element_quad/element_quad_dg.py      | 1 +
 skfem/element/element_quad/element_quad_s2.py      | 1 +
 skfem/element/element_quad/element_quadp.py        | 1 +
 skfem/element/element_tet/element_tet_mini.py      | 1 +
 skfem/element/element_tet/element_tet_n0.py        | 1 +
 skfem/element/element_tet/element_tet_p0.py        | 1 +
 skfem/element/element_tet/element_tet_p1.py        | 1 +
 skfem/element/element_tet/element_tet_p2.py        | 1 +
 skfem/element/element_tet/element_tet_rt0.py       | 1 +
 skfem/element/element_tri/element_tri_argyris.py   | 1 +
 skfem/element/element_tri/element_tri_dg.py        | 1 +
 skfem/element/element_tri/element_tri_mini.py      | 1 +
 skfem/element/element_tri/element_tri_morley.py    | 1 +
 skfem/element/element_tri/element_tri_p0.py        | 1 +
 skfem/element/element_tri/element_tri_p1.py        | 1 +
 skfem/element/element_tri/element_tri_p2.py        | 1 +
 skfem/element/element_tri/element_tri_rt0.py       | 1 +
 26 files changed, 26 insertions(+)

diff --git a/skfem/element/element_line/element_line_hermite.py b/skfem/element/element_line/element_line_hermite.py
index 98e156d8f..4f7ddc103 100644
--- a/skfem/element/element_line/element_line_hermite.py
+++ b/skfem/element/element_line/element_line_hermite.py
@@ -5,6 +5,7 @@
 
 
 class ElementLineHermite(ElementGlobal):
+
     nodal_dofs = 2
     dim = 1
     maxdeg = 3
diff --git a/skfem/element/element_line/element_line_mini.py b/skfem/element/element_line/element_line_mini.py
index eba83d806..686813599 100644
--- a/skfem/element/element_line/element_line_mini.py
+++ b/skfem/element/element_line/element_line_mini.py
@@ -5,6 +5,7 @@
 
 
 class ElementLineMini(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 1
diff --git a/skfem/element/element_line/element_line_p1.py b/skfem/element/element_line/element_line_p1.py
index 6f0182bf5..e9001021e 100644
--- a/skfem/element/element_line/element_line_p1.py
+++ b/skfem/element/element_line/element_line_p1.py
@@ -5,6 +5,7 @@
 
 
 class ElementLineP1(ElementH1):
+
     nodal_dofs = 1
     dim = 1
     maxdeg = 1
diff --git a/skfem/element/element_line/element_line_p2.py b/skfem/element/element_line/element_line_p2.py
index ece0e68e1..4f02e0f76 100644
--- a/skfem/element/element_line/element_line_p2.py
+++ b/skfem/element/element_line/element_line_p2.py
@@ -5,6 +5,7 @@
 
 
 class ElementLineP2(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 1
diff --git a/skfem/element/element_line/element_line_pp.py b/skfem/element/element_line/element_line_pp.py
index 682d33e0a..09453ae96 100644
--- a/skfem/element/element_line/element_line_pp.py
+++ b/skfem/element/element_line/element_line_pp.py
@@ -8,6 +8,7 @@
 
 
 class ElementLinePp(ElementH1):
+
     nodal_dofs = 1
     dim = 1
     mesh_type = MeshLine
diff --git a/skfem/element/element_quad/element_quad0.py b/skfem/element/element_quad/element_quad0.py
index ef814e25c..540daf7d3 100644
--- a/skfem/element/element_quad/element_quad0.py
+++ b/skfem/element/element_quad/element_quad0.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuad0(ElementH1):
+
     interior_dofs = 1
     dim = 2
     maxdeg = 0
diff --git a/skfem/element/element_quad/element_quad1.py b/skfem/element/element_quad/element_quad1.py
index ba30a243d..ec08c455d 100644
--- a/skfem/element/element_quad/element_quad1.py
+++ b/skfem/element/element_quad/element_quad1.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuad1(ElementH1):
+
     nodal_dofs = 1
     dim = 2
     maxdeg = 2
diff --git a/skfem/element/element_quad/element_quad2.py b/skfem/element/element_quad/element_quad2.py
index 3da131564..ec9f038d8 100644
--- a/skfem/element/element_quad/element_quad2.py
+++ b/skfem/element/element_quad/element_quad2.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuad2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     interior_dofs = 1
diff --git a/skfem/element/element_quad/element_quad_bfs.py b/skfem/element/element_quad/element_quad_bfs.py
index 49b5d2525..41e1499ca 100644
--- a/skfem/element/element_quad/element_quad_bfs.py
+++ b/skfem/element/element_quad/element_quad_bfs.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuadBFS(ElementGlobal):
+
     nodal_dofs = 4
     dim = 2
     maxdeg = 6
diff --git a/skfem/element/element_quad/element_quad_dg.py b/skfem/element/element_quad/element_quad_dg.py
index 3826a771e..8eb2c2488 100644
--- a/skfem/element/element_quad/element_quad_dg.py
+++ b/skfem/element/element_quad/element_quad_dg.py
@@ -3,6 +3,7 @@
 
 
 class ElementQuadDG(ElementH1):
+
     dim = 2
     mesh_type = MeshQuad
 
diff --git a/skfem/element/element_quad/element_quad_s2.py b/skfem/element/element_quad/element_quad_s2.py
index aa41338aa..a75ba1d31 100644
--- a/skfem/element/element_quad/element_quad_s2.py
+++ b/skfem/element/element_quad/element_quad_s2.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuadS2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     dim = 2
diff --git a/skfem/element/element_quad/element_quadp.py b/skfem/element/element_quad/element_quadp.py
index c451b8dbf..93b55b607 100644
--- a/skfem/element/element_quad/element_quadp.py
+++ b/skfem/element/element_quad/element_quadp.py
@@ -5,6 +5,7 @@
 
 
 class ElementQuadP(ElementLinePp):
+
     nodal_dofs = 1
     dim = 2
     mesh_type = MeshQuad
diff --git a/skfem/element/element_tet/element_tet_mini.py b/skfem/element/element_tet/element_tet_mini.py
index 4887adbce..ab278fbc6 100644
--- a/skfem/element/element_tet/element_tet_mini.py
+++ b/skfem/element/element_tet/element_tet_mini.py
@@ -5,6 +5,7 @@
 
 
 class ElementTetMini(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 3
diff --git a/skfem/element/element_tet/element_tet_n0.py b/skfem/element/element_tet/element_tet_n0.py
index 903d9f291..141f6b933 100644
--- a/skfem/element/element_tet/element_tet_n0.py
+++ b/skfem/element/element_tet/element_tet_n0.py
@@ -4,6 +4,7 @@
 
 
 class ElementTetN0(ElementHcurl):
+
     edge_dofs = 1
     dim = 3
     maxdeg = 1
diff --git a/skfem/element/element_tet/element_tet_p0.py b/skfem/element/element_tet/element_tet_p0.py
index 355b2f186..21a117c97 100644
--- a/skfem/element/element_tet/element_tet_p0.py
+++ b/skfem/element/element_tet/element_tet_p0.py
@@ -4,6 +4,7 @@
 
 
 class ElementTetP0(ElementH1):
+
     interior_dofs = 1
     dim = 3
     maxdeg = 0
diff --git a/skfem/element/element_tet/element_tet_p1.py b/skfem/element/element_tet/element_tet_p1.py
index 191cbd7ce..3e107976a 100644
--- a/skfem/element/element_tet/element_tet_p1.py
+++ b/skfem/element/element_tet/element_tet_p1.py
@@ -5,6 +5,7 @@
 
 
 class ElementTetP1(ElementH1):
+
     nodal_dofs = 1
     dim = 3
     maxdeg = 1
diff --git a/skfem/element/element_tet/element_tet_p2.py b/skfem/element/element_tet/element_tet_p2.py
index 8477f3321..a0c1c57a1 100644
--- a/skfem/element/element_tet/element_tet_p2.py
+++ b/skfem/element/element_tet/element_tet_p2.py
@@ -5,6 +5,7 @@
 
 
 class ElementTetP2(ElementH1):
+
     nodal_dofs = 1
     edge_dofs = 1
     dim = 3
diff --git a/skfem/element/element_tet/element_tet_rt0.py b/skfem/element/element_tet/element_tet_rt0.py
index 19a6fc119..8ee00f0fe 100644
--- a/skfem/element/element_tet/element_tet_rt0.py
+++ b/skfem/element/element_tet/element_tet_rt0.py
@@ -5,6 +5,7 @@
 
 
 class ElementTetRT0(ElementHdiv):
+
     facet_dofs = 1
     dim = 3
     maxdeg = 1
diff --git a/skfem/element/element_tri/element_tri_argyris.py b/skfem/element/element_tri/element_tri_argyris.py
index 5c97f1f62..73da0e2e0 100644
--- a/skfem/element/element_tri/element_tri_argyris.py
+++ b/skfem/element/element_tri/element_tri_argyris.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriArgyris(ElementGlobal):
+
     nodal_dofs = 6
     facet_dofs = 1
     dim = 2
diff --git a/skfem/element/element_tri/element_tri_dg.py b/skfem/element/element_tri/element_tri_dg.py
index dedcd243e..85b9de840 100644
--- a/skfem/element/element_tri/element_tri_dg.py
+++ b/skfem/element/element_tri/element_tri_dg.py
@@ -3,6 +3,7 @@
 
 
 class ElementTriDG(ElementH1):
+
     dim = 2
     mesh_type = MeshTri
 
diff --git a/skfem/element/element_tri/element_tri_mini.py b/skfem/element/element_tri/element_tri_mini.py
index 3c8db7626..b20b77a2f 100644
--- a/skfem/element/element_tri/element_tri_mini.py
+++ b/skfem/element/element_tri/element_tri_mini.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriMini(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 2
diff --git a/skfem/element/element_tri/element_tri_morley.py b/skfem/element/element_tri/element_tri_morley.py
index 15187661c..88d506684 100644
--- a/skfem/element/element_tri/element_tri_morley.py
+++ b/skfem/element/element_tri/element_tri_morley.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriMorley(ElementGlobal):
+
     nodal_dofs = 1
     facet_dofs = 1
     dim = 2
diff --git a/skfem/element/element_tri/element_tri_p0.py b/skfem/element/element_tri/element_tri_p0.py
index 6a63f0a96..c754424bf 100644
--- a/skfem/element/element_tri/element_tri_p0.py
+++ b/skfem/element/element_tri/element_tri_p0.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriP0(ElementH1):
+
     interior_dofs = 1
     dim = 2
     maxdeg = 0
diff --git a/skfem/element/element_tri/element_tri_p1.py b/skfem/element/element_tri/element_tri_p1.py
index c548844f5..8c101d539 100644
--- a/skfem/element/element_tri/element_tri_p1.py
+++ b/skfem/element/element_tri/element_tri_p1.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriP1(ElementH1):
+
     nodal_dofs = 1
     dim = 2
     maxdeg = 1
diff --git a/skfem/element/element_tri/element_tri_p2.py b/skfem/element/element_tri/element_tri_p2.py
index 898fd60f2..ccf243828 100644
--- a/skfem/element/element_tri/element_tri_p2.py
+++ b/skfem/element/element_tri/element_tri_p2.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriP2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     dim = 2
diff --git a/skfem/element/element_tri/element_tri_rt0.py b/skfem/element/element_tri/element_tri_rt0.py
index dc7d81d37..1446df175 100644
--- a/skfem/element/element_tri/element_tri_rt0.py
+++ b/skfem/element/element_tri/element_tri_rt0.py
@@ -5,6 +5,7 @@
 
 
 class ElementTriRT0(ElementHdiv):
+
     facet_dofs = 1
     dim = 2
     maxdeg = 1

From f7db38749c5e0bc33f5b1f4ba79ac95304c4c757 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:41:04 +0300
Subject: [PATCH 05/13] Make whitespace consistent

---
 tests/test_examples.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/tests/test_examples.py b/tests/test_examples.py
index c376dcf0e..47dbb384c 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -4,24 +4,28 @@
 
 
 class TestEx01(TestCase):
+
     def runTest(self):
         import docs.examples.ex01 as ex01
         self.assertAlmostEqual(np.max(ex01.x), 0.07344576657)
 
 
 class TestEx02(TestCase):
+
     def runTest(self):
         import docs.examples.ex02 as ex02
         self.assertAlmostEqual(np.max(ex02.x), 0.001217973811129439)
 
 
 class TestEx03(TestCase):
+
     def runTest(self):
         import docs.examples.ex03 as ex03
         self.assertAlmostEqual(ex03.L[0], 0.00418289)
 
 
 class TestEx04(TestCase):
+
     def runTest(self):
         import docs.examples.ex04 as ex04
         self.assertAlmostEqual(np.max(ex04.vonmises1), 64.57919816083539)
@@ -29,42 +33,49 @@ def runTest(self):
 
 
 class TestEx05(TestCase):
+
     def runTest(self):
         import docs.examples.ex05 as ex05
         self.assertAlmostEqual(np.max(ex05.x), 0.93570751751091152)
 
 
 class TestEx06(TestCase):
+
     def runTest(self):
         import docs.examples.ex06 as ex06
         self.assertAlmostEqual(np.max(ex06.x), 0.073651530833125131)
 
 
 class TestEx07(TestCase):
+
     def runTest(self):
         import docs.examples.ex07 as ex07
         self.assertAlmostEqual(np.max(ex07.x), 0.07869083767545548)
 
 
 class TestEx08(TestCase):
+
     def runTest(self):
         import docs.examples.ex08 as ex08  # noqa
         # only run the initialization, nothing to test
 
 
 class TestEx09(TestCase):
+
     def runTest(self):
         import docs.examples.ex09 as ex09
         self.assertAlmostEqual(np.max(ex09.x), 0.05528520791811886)
 
 
 class TestEx10(TestCase):
+
     def runTest(self):
         import docs.examples.ex10 as ex10
         self.assertAlmostEqual(np.mean(ex10.x), 0.277931521728906)
 
 
 class TestEx11(TestCase):
+
     def runTest(self):
         import docs.examples.ex11 as ex11
         u = ex11.u
@@ -75,6 +86,7 @@ def runTest(self):
 
 
 class TestEx12(TestCase):
+
     def runTest(self):
         import docs.examples.ex12 as ex
         self.assertAlmostEqual(ex.area, np.pi, delta=1e-2)
@@ -83,6 +95,7 @@ def runTest(self):
 
 
 class TestEx13(TestCase):
+
     def runTest(self):
         import docs.examples.ex13 as ex
         u = ex.u
@@ -91,6 +104,7 @@ def runTest(self):
 
 
 class TestEx14(TestCase):
+
     def runTest(self):
         import docs.examples.ex14
         u = docs.examples.ex14.u
@@ -99,12 +113,14 @@ def runTest(self):
 
 
 class TestEx15(TestCase):
+
     def runTest(self):
         import docs.examples.ex15 as ex15
         self.assertTrue(np.max(ex15.x) - 0.1234567 < 1e-5)
 
 
 class TestEx16(TestCase):
+
     def runTest(self):
         import docs.examples.ex16 as ex16
         self.assertTrue(np.linalg.norm(np.array([0, 2, 6, 12, 20, 30])
@@ -113,24 +129,28 @@ def runTest(self):
 
 
 class TestEx17(TestCase):
+
     def runTest(self):
         import docs.examples.ex17 as ex  # noqa
         # TODO improve
 
 
 class TestEx18(TestCase):
+
     def runTest(self):
         import docs.examples.ex18 as ex  # noqa
         # TODO improve
 
 
 class TestEx19(TestCase):
+
     def runTest(self):
         import docs.examples.ex19 as ex  # noqa
         # TODO improve
 
 
 class TestEx20(TestCase):
+
     def runTest(self):
         import docs.examples.ex20 as ex
         psi0 = ex.psi0
@@ -138,6 +158,7 @@ def runTest(self):
 
 
 class TestEx21(TestCase):
+
     def runTest(self):
         import docs.examples.ex21 as ex
         x = ex.x
@@ -148,6 +169,7 @@ def runTest(self):
 
 
 class TestEx22(TestCase):
+
     def runTest(self):
         import docs.examples.ex22 as ex
         u = ex.u
@@ -156,6 +178,7 @@ def runTest(self):
 
 
 class TestEx23(TestCase):
+
     def runTest(self):
         import docs.examples.ex23 as ex
         self.assertAlmostEqual(max(ex.lmbda_list), ex.turning_point,
@@ -163,11 +186,13 @@ def runTest(self):
 
 
 class TestEx24(TestCase):
+
     def runTest(self):
         import docs.examples.ex24 as ex24  # noqa
 
 
 class TestEx25(TestCase):
+
     def runTest(self):
         import docs.examples.ex25 as ex25
         mu = np.mean(ex25.t)
@@ -176,11 +201,13 @@ def runTest(self):
 
 
 class TestEx26(TestCase):
+
     def runTest(self):
         import docs.examples.ex26 as ex26  # noqa
 
 
 class TestEx27(TestCase):
+
     def runTest(self):
         import docs.examples.ex27 as ex
         _, psi = ex.psi.popitem()
@@ -189,12 +216,14 @@ def runTest(self):
 
 
 class TestEx28(TestCase):
+
     def runTest(self):
         from docs.examples.ex28 import exit_interface_temperature as t
         self.assertAlmostEqual(*t.values(), delta=2e-4)
 
 
 class TestEx29(TestCase):
+
     def runTest(self):
         from docs.examples.ex29 import c
         wavespeed = tuple(
@@ -205,6 +234,7 @@ def runTest(self):
 
 
 class TestEx30(TestCase):
+
     def runTest(self):
         from docs.examples.ex30 import psi0
         self.assertAlmostEqual(psi0, 0.162/128, delta=1e-6)
@@ -239,6 +269,7 @@ def runTest(self):
 
 
 class TestEx35(TestCase):
+
     def runTest(self):
         from docs.examples.ex35 import Z
         # exact value depends also on mesh generation,
@@ -247,6 +278,7 @@ def runTest(self):
         self.assertAlmostEqual(Z, 52.563390368494424, delta=1e-1)
 
 class TestEx36(TestCase):
+
     def runTest(self):
         from docs.examples.ex36 import dp
         # make test more robust or this is fine?

From c72f888f9d3ef6a5fa6a4cc765bca4350463d473 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:44:37 +0300
Subject: [PATCH 06/13] Add exactness test for ElementHex2

---
 tests/test_manufactured.py | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/tests/test_manufactured.py b/tests/test_manufactured.py
index b6f79631d..886cfc004 100644
--- a/tests/test_manufactured.py
+++ b/tests/test_manufactured.py
@@ -9,7 +9,7 @@
 from skfem.element import (ElementHex1, ElementHexS2,
                            ElementLineP1, ElementLineP2, ElementLineMini, 
                            ElementQuad1, ElementQuad2, ElementTetP1,
-                           ElementTriP2)
+                           ElementTriP2, ElementHex2)
 from skfem.assembly import FacetBasis, InteriorBasis
 from skfem import asm, condense, solve, LinearForm
 
@@ -133,6 +133,7 @@ class LineNeumann1DMini(LineNeumann1D):
     
 
 class TestExactHexElement(unittest.TestCase):
+
     mesh = MeshHex
     elem = ElementHex1
     funs = [
@@ -146,7 +147,7 @@ def set_bc(self, fun, basis):
     def runTest(self):
 
         m = self.mesh()
-        m.refine(4)
+        m.refine(3)
 
         ib = InteriorBasis(m, self.elem())
 
@@ -163,10 +164,22 @@ def runTest(self):
 
 
 class TestExactHexS2(TestExactHexElement):
+
     elem = ElementHexS2
+    funs = [
+        lambda x: 1 + 0 * x[0],
+    ]
+
+    def set_bc(self, fun, basis):
+        return fun(basis.doflocs)
 
+
+class TestExactHex2(TestExactHexElement):
+
+    elem = ElementHex2
     funs = [
         lambda x: 1 + 0 * x[0],
+        lambda x: 1 + x[0] * x[1] * x[2],
     ]
 
     def set_bc(self, fun, basis):
@@ -174,6 +187,7 @@ def set_bc(self, fun, basis):
 
 
 class TestExactQuadElement(TestExactHexElement):
+
     mesh = MeshQuad
     elem = ElementQuad1
     funs = [
@@ -183,6 +197,7 @@ class TestExactQuadElement(TestExactHexElement):
 
 
 class TestExactTetElement(TestExactHexElement):
+
     mesh = MeshTet
     elem = ElementTetP1
     funs = [
@@ -192,6 +207,7 @@ class TestExactTetElement(TestExactHexElement):
 
 
 class TestExactTriElementP2(TestExactHexElement):
+
     mesh = MeshTri
     elem = ElementTriP2
     funs = [
@@ -204,6 +220,7 @@ def set_bc(self, fun, basis):
 
 
 class TestExactQuadElement2(TestExactTriElementP2):
+
     mesh = MeshQuad
     elem = ElementQuad2
 

From cc1aa4279c13ba899e4c5781bb0de18b3f16f263 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:45:58 +0300
Subject: [PATCH 07/13] Improve whitespace

---
 tests/test_mapping.py | 6 ++++++
 tests/test_mesh.py    | 2 ++
 2 files changed, 8 insertions(+)

diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index 2e33b47fa..2b1865dca 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -9,6 +9,7 @@
 
 class TestIsoparamNormals(unittest.TestCase):
     """Test that normals on x[i] == 0 are correct."""
+
     mesh = MeshHex
     elem = ElementHex1
 
@@ -30,6 +31,7 @@ def runTest(self):
 
 
 class TestIsoparamNormalsQuad(TestIsoparamNormals):
+
     mesh = MeshQuad
     elem = ElementQuad1
 
@@ -106,22 +108,26 @@ def runTest(self):
 
 
 class TestMortarPairTriQuad(TestMortarPair):
+
     mesh1_type = MeshTri
     mesh2_type = MeshQuad
 
 
 class TestMortarPairQuadQuad(TestMortarPair):
+
     mesh1_type = MeshQuad
     mesh2_type = MeshQuad
 
 
 class TestMortarPairNoMatch1(TestMortarPair):
+
     mesh1_type = MeshQuad
     mesh2_type = MeshTri
     translate_y = 0.1
 
 
 class TestMortarPairNoMatch2(TestMortarPair):
+
     mesh1_type = MeshQuad
     mesh2_type = MeshTri
     translate_y = -np.pi / 10.
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 96d591824..0857e3014 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -118,11 +118,13 @@ def runTest(self):
 
 
 class SaveLoadCycleHex(SaveLoadCycle):
+
     cls = MeshHex
 
 
 class SerializeUnserializeCycle(unittest.TestCase):
     """Check to_dict/initialize cycles."""
+
     clss = [MeshTet,
             MeshTri,
             MeshHex,

From b444e025aa8a6de1eb1cca20cf8a3b1fea8f8aaf Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 13:50:29 +0300
Subject: [PATCH 08/13] Move ElementLinePp initialization test to test_elements

---
 tests/test_element_line_pp.py | 12 ------------
 tests/test_elements.py        |  8 ++++++++
 2 files changed, 8 insertions(+), 12 deletions(-)
 delete mode 100644 tests/test_element_line_pp.py

diff --git a/tests/test_element_line_pp.py b/tests/test_element_line_pp.py
deleted file mode 100644
index e576697b8..000000000
--- a/tests/test_element_line_pp.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from unittest import TestCase
-from skfem.element.element_line.element_line_pp import ElementLinePp
-
-
-class TestElementLinePp(TestCase):
-
-    def test_p_less_than_1_error(self):
-        """ Tests that exception is thrown when trying to initialize element
-        with p < 1. """
-        with self.assertRaises(ValueError):
-            ElementLinePp(0)
-
diff --git a/tests/test_elements.py b/tests/test_elements.py
index 6db14b845..ca76d4426 100644
--- a/tests/test_elements.py
+++ b/tests/test_elements.py
@@ -241,5 +241,13 @@ def runTest(self):
             self.assertAlmostEqual(out, 1, msg='failed for {}'.format(elem))
 
 
+class TestElementLinePp(TestCase):
+
+    def test_p_less_than_1_error(self):
+        """Tests that exception is thrown when initializing with p < 1."""
+        with self.assertRaises(ValueError):
+            ElementLinePp(0)
+
+
 if __name__ == '__main__':
     main()

From a5d2dcf06a1905aeae7f7f1d987fd8a923876db6 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 14:29:37 +0300
Subject: [PATCH 09/13] Clip inverse mapping Newton iteration to a range for
 robustness

---
 skfem/mapping/mapping_isoparametric.py |  4 ++--
 tests/test_mapping.py                  | 16 ++++++++++++++--
 2 files changed, 16 insertions(+), 4 deletions(-)

diff --git a/skfem/mapping/mapping_isoparametric.py b/skfem/mapping/mapping_isoparametric.py
index 86149421c..78410423e 100644
--- a/skfem/mapping/mapping_isoparametric.py
+++ b/skfem/mapping/mapping_isoparametric.py
@@ -129,14 +129,14 @@ def detDG(self, X: ndarray, find: Optional[ndarray] = None):
         else:
             raise NotImplementedError
 
-    def invF(self, x, tind=None, newton_max_iters=50, newton_tol=1e-8):
+    def invF(self, x, tind=None, newton_max_iters=50, newton_tol=1e-12):
         """Newton iteration for evaluating inverse isoparametric mapping."""
         X = np.zeros(x.shape) + .5
         for _ in range(newton_max_iters):
             F = self.F(X, tind)
             invDF = self.invDF(X, tind)
             dX = np.einsum('ijkl,jkl->ikl', invDF, x - F)
-            X = X + dX
+            X = np.clip(X + dX, 0., 1.)
             if (np.linalg.norm(dX, 1, (0, 2)) < newton_tol).all():
                 return X
         raise Exception(("Newton iteration didn't converge "
diff --git a/tests/test_mapping.py b/tests/test_mapping.py
index 2b1865dca..497a8664c 100644
--- a/tests/test_mapping.py
+++ b/tests/test_mapping.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from skfem.mesh import MeshHex, MeshQuad, MeshTri
-from skfem.element import ElementHex1, ElementQuad1
+from skfem.element import ElementHex1, ElementQuad1, ElementHex2
 from skfem.assembly import FacetBasis
 from skfem.mapping.mapping_mortar import MappingMortar
 
@@ -36,6 +36,11 @@ class TestIsoparamNormalsQuad(TestIsoparamNormals):
     elem = ElementQuad1
 
 
+class TestIsoparamNormalsHex2(TestIsoparamNormals):
+
+    elem = ElementHex2
+
+
 class TestInverseMapping(unittest.TestCase):
     """Test that inverse mapping works for non-rectangular elements."""
 
@@ -49,7 +54,8 @@ def initialize_meshes(self):
         return m
 
     def within_refelem(self, y):
-        return (np.abs(y) < 1.0 + 1e-12).all()
+        return ((np.abs(y) < 1. + 1e-12).all()
+                and (np.abs(y) > 0. - 1e-12).all())
 
     def runTest(self):
         m = self.initialize_meshes()
@@ -79,6 +85,12 @@ def initialize_meshes(self):
         return m
 
 
+class TestInverseMappingHex2(TestInverseMappingHex):
+    """This should be equivalent to TestInverseMappingHex."""
+
+    element = ElementHex2
+
+
 class TestMortarPair(unittest.TestCase):
     """Check that mapped points match."""
 

From cdf7384db2174df09102a2952ecb1033b618c862 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 15:47:53 +0300
Subject: [PATCH 10/13] Fix order of doflocs

---
 skfem/element/element_hex/element_hex2.py | 14 +++++++-------
 tests/test_convergence.py                 |  3 +++
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/skfem/element/element_hex/element_hex2.py b/skfem/element/element_hex/element_hex2.py
index 89b83e48b..ccd9a7203 100644
--- a/skfem/element/element_hex/element_hex2.py
+++ b/skfem/element/element_hex/element_hex2.py
@@ -22,13 +22,7 @@ class ElementHex2(ElementH1):
                         [0., 1., 0.],
                         [0., 0., 1.],
                         [0., 0., 0.],
-                        [1., .5, .5],
-                        [.5, .5, 1.],
-                        [.5, 1., .5],
-                        [.5, .0, .5],
-                        [.5, .5, 0.],
-                        [0., .5, .5],
-                        [1., 1., .5],
+                        [1., 1., .5], # edges
                         [1., .5, 1.],
                         [.5, 1., 1.],
                         [1., .5, 0.],
@@ -40,6 +34,12 @@ class ElementHex2(ElementH1):
                         [.5, 0., 0.],
                         [0., .5, 0.],
                         [0., 0., .5],
+                        [1., .5, .5], # facets
+                        [.5, .5, 1.],
+                        [.5, 1., .5],
+                        [.5, .0, .5],
+                        [.5, .5, 0.],
+                        [0., .5, .5],
                         [.5, .5, .5]])
     mesh_type = MeshHex
 
diff --git a/tests/test_convergence.py b/tests/test_convergence.py
index 67d4d9cf4..8dc6f540c 100644
--- a/tests/test_convergence.py
+++ b/tests/test_convergence.py
@@ -226,6 +226,9 @@ def setUp(self):
 
 class ConvergenceHex2(ConvergenceHexS2):
 
+    rateL2 = 2.92
+    rateH1 = 2.01
+
     def create_basis(self, m):
         e = ElementHex2()
         return InteriorBasis(m, e)

From 0381cf72e54f69eb790e7761e38be2025456f224 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 15:48:52 +0300
Subject: [PATCH 11/13] Add changelog entry

---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index 4429be368..514ffb54c 100644
--- a/README.md
+++ b/README.md
@@ -168,6 +168,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Unreleased
 
+#### Added
+- `ElementHex2`, a triquadratic hexahedral element
+
 ### [2.0.0] - 2020-08-21
 
 #### Added

From dc59bd288ccafb7b274b5273b80e2ab4c2f2c659 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 16:09:31 +0300
Subject: [PATCH 12/13] Fix flake8 issues

---
 skfem/element/element_hex/element_hex2.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skfem/element/element_hex/element_hex2.py b/skfem/element/element_hex/element_hex2.py
index ccd9a7203..5777cf7cc 100644
--- a/skfem/element/element_hex/element_hex2.py
+++ b/skfem/element/element_hex/element_hex2.py
@@ -22,7 +22,7 @@ class ElementHex2(ElementH1):
                         [0., 1., 0.],
                         [0., 0., 1.],
                         [0., 0., 0.],
-                        [1., 1., .5], # edges
+                        [1., 1., .5],  # edges
                         [1., .5, 1.],
                         [.5, 1., 1.],
                         [1., .5, 0.],
@@ -34,7 +34,7 @@ class ElementHex2(ElementH1):
                         [.5, 0., 0.],
                         [0., .5, 0.],
                         [0., 0., .5],
-                        [1., .5, .5], # facets
+                        [1., .5, .5],  # facets
                         [.5, .5, 1.],
                         [.5, 1., .5],
                         [.5, .0, .5],
@@ -134,7 +134,7 @@ def _power_basis_dy(self, X):
             0. * X[0],
             X[0] * X[2],
             X[0] ** 2 * X[2],
-            X[0] * 2. *  X[1] * X[2],
+            X[0] * 2. * X[1] * X[2],
             X[0] * X[2] ** 2,
             X[0] ** 2 * 2. * X[1] * X[2],
             X[0] * 2. * X[1] * X[2] ** 2,

From 996b0b2bc0b863ce7a04bd356e7d26500fc343c4 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sun, 13 Sep 2020 20:42:57 +0300
Subject: [PATCH 13/13] Raise error if meshio mesh type not supported

---
 skfem/io/meshio.py | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/skfem/io/meshio.py b/skfem/io/meshio.py
index f12782dae..b361be3ce 100644
--- a/skfem/io/meshio.py
+++ b/skfem/io/meshio.py
@@ -26,22 +26,32 @@ def from_meshio(m, force_mesh_type=None):
     Parameters
     ----------
     m
-        The mesh object from meshio.
+        The mesh from meshio.
     force_mesh_type
         An optional string forcing the mesh type if automatic detection
         fails. See skfem.importers.meshio.MESH_TYPE_MAPPING for possible
         values.
 
+    Returns
+    -------
+    A :class:`~skfem.mesh.Mesh` object.
+
     """
 
     cells = m.cells_dict
 
     if force_mesh_type is None:
+        meshio_type = None
+
         for k, v in MESH_TYPE_MAPPING.items():
             # find first if match
             if k in cells:
                 meshio_type, mesh_type = k, v
                 break
+
+        if meshio_type is None:
+            raise NotImplementedError("Mesh type(s) not supported "
+                                      "in import: {}.".format(cells.keys()))
     else:
         meshio_type, mesh_type = (force_mesh_type,
                                   MESH_TYPE_MAPPING[force_mesh_type])