Add ElementHex2

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

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",

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

skfem/element/element_hex/element_hex2.py

@@ -0,0 +1,186 @@
+import numpy as np
+from numpy.linalg import inv
+
+from ..element_h1 import ElementH1
+from ...mesh.mesh3d import MeshHex
+
+
+class ElementHex2(ElementH1):
+
+    nodal_dofs = 1
+    facet_dofs = 1
+    edge_dofs = 1
+    interior_dofs = 1
+    dim = 3
+    maxdeg = 6
+    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., 1., .5],  # edges
+                        [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],
+                        [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
+
+    def __init__(self):
+        X = self.doflocs.T
+        V = self._power_basis(X)
+        self.invV = inv(V)
+        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],
+            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,
+        ])
+
+    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,
+        ])
+
+    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],
+        ])
+
+    def lbasis(self, X, i):
+        if i >= 27 or i < 0:
+            self._index_error()
+        return (
+            np.einsum('i...,i...', self.invV[i], self._power_basis(X)),
+            np.array([
+                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)),
+            ])
+        )

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

skfem/element/element_line/element_line_hermite.py

@@ -5,6 +5,7 @@
 
 
 class ElementLineHermite(ElementGlobal):
+
     nodal_dofs = 2
     dim = 1
     maxdeg = 3

skfem/element/element_line/element_line_mini.py

@@ -5,6 +5,7 @@
 
 
 class ElementLineMini(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 1

skfem/element/element_line/element_line_p1.py

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

skfem/element/element_line/element_line_p2.py

@@ -5,6 +5,7 @@
 
 
 class ElementLineP2(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 1

skfem/element/element_line/element_line_pp.py

@@ -8,6 +8,7 @@
 
 
 class ElementLinePp(ElementH1):
+
     nodal_dofs = 1
     dim = 1
     mesh_type = MeshLine

skfem/element/element_quad/element_quad0.py

@@ -5,6 +5,7 @@
 
 
 class ElementQuad0(ElementH1):
+
     interior_dofs = 1
     dim = 2
     maxdeg = 0

skfem/element/element_quad/element_quad1.py

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

skfem/element/element_quad/element_quad2.py

@@ -5,6 +5,7 @@
 
 
 class ElementQuad2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     interior_dofs = 1

skfem/element/element_quad/element_quad_bfs.py

@@ -5,6 +5,7 @@
 
 
 class ElementQuadBFS(ElementGlobal):
+
     nodal_dofs = 4
     dim = 2
     maxdeg = 6

skfem/element/element_quad/element_quad_dg.py

@@ -3,6 +3,7 @@
 
 
 class ElementQuadDG(ElementH1):
+
     dim = 2
     mesh_type = MeshQuad
 

skfem/element/element_quad/element_quad_s2.py

@@ -5,6 +5,7 @@
 
 
 class ElementQuadS2(ElementH1):
+
     nodal_dofs = 1
     facet_dofs = 1
     dim = 2

skfem/element/element_quad/element_quadp.py

@@ -5,6 +5,7 @@
 
 
 class ElementQuadP(ElementLinePp):
+
     nodal_dofs = 1
     dim = 2
     mesh_type = MeshQuad

skfem/element/element_tet/element_tet_mini.py

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

skfem/element/element_tet/element_tet_n0.py

@@ -4,6 +4,7 @@
 
 
 class ElementTetN0(ElementHcurl):
+
     edge_dofs = 1
     dim = 3
     maxdeg = 1

skfem/element/element_tet/element_tet_p0.py

@@ -4,6 +4,7 @@
 
 
 class ElementTetP0(ElementH1):
+
     interior_dofs = 1
     dim = 3
     maxdeg = 0

skfem/element/element_tet/element_tet_p1.py

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

skfem/element/element_tet/element_tet_p2.py

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

skfem/element/element_tet/element_tet_rt0.py

@@ -5,6 +5,7 @@
 
 
 class ElementTetRT0(ElementHdiv):
+
     facet_dofs = 1
     dim = 3
     maxdeg = 1

skfem/element/element_tri/element_tri_argyris.py

@@ -5,6 +5,7 @@
 
 
 class ElementTriArgyris(ElementGlobal):
+
     nodal_dofs = 6
     facet_dofs = 1
     dim = 2

skfem/element/element_tri/element_tri_dg.py

@@ -3,6 +3,7 @@
 
 
 class ElementTriDG(ElementH1):
+
     dim = 2
     mesh_type = MeshTri
 

skfem/element/element_tri/element_tri_mini.py

@@ -5,6 +5,7 @@
 
 
 class ElementTriMini(ElementH1):
+
     nodal_dofs = 1
     interior_dofs = 1
     dim = 2