Merge pull request #415 from adtzlr/add-quadratic-hex-region-boundary

Add boundary regions for quadratic hexahedrons

CHANGELOG.md

@@ -5,6 +5,7 @@ All notable changes to this project will be documented in this file. The format
 
 ### Added
 - Add boundary regions `RegionQuadraticQuadBoundary` and `RegionBiQuadraticQuadBoundary` for quadratic quads.
+- Add boundary regions `RegionQuadraticHexahedronBoundary` and `RegionTriQuadraticHexahedronBoundary` for quadratic hexahedrons.
 
 ### Changed
 - Change `einsumt` from an optional to a required dependency.

src/felupe/__init__.py

@@ -73,6 +73,7 @@
     RegionQuad,
     RegionQuadBoundary,
     RegionQuadraticHexahedron,
+    RegionQuadraticHexahedronBoundary,
     RegionQuadraticQuad,
     RegionQuadraticQuadBoundary,
     RegionQuadraticTetra,
@@ -82,6 +83,7 @@
     RegionTriangle,
     RegionTriangleMINI,
     RegionTriQuadraticHexahedron,
+    RegionTriQuadraticHexahedronBoundary,
 )
 
 try:

src/felupe/region/__init__.py

@@ -11,6 +11,7 @@
     RegionQuad,
     RegionQuadBoundary,
     RegionQuadraticHexahedron,
+    RegionQuadraticHexahedronBoundary,
     RegionQuadraticQuad,
     RegionQuadraticQuadBoundary,
     RegionQuadraticTetra,
@@ -20,4 +21,5 @@
     RegionTriangle,
     RegionTriangleMINI,
     RegionTriQuadraticHexahedron,
+    RegionTriQuadraticHexahedronBoundary,
 )

src/felupe/region/_boundary.py

@@ -119,10 +119,10 @@ def boundary_cells_hexahedron(mesh):
     "Convert the cells array of a hex mesh into a boundary cells array."
 
     # faces (boundary) of a hexahedron
-    i = [0, 1, 1, 2, 0, 4]
-    j = [3, 2, 0, 3, 1, 5]
-    k = [7, 6, 4, 7, 2, 6]
-    l = [4, 5, 5, 6, 3, 7]
+    i = [0, 2, 1, 3, 0, 5]
+    j = [3, 1, 0, 2, 1, 4]
+    k = [7, 5, 4, 6, 2, 7]
+    l = [4, 6, 5, 7, 3, 6]
 
     cells_faces = np.dstack(
         (
@@ -136,11 +136,10 @@ def boundary_cells_hexahedron(mesh):
     # complementary faces for the creation of "boundary" hexahedrons
     # (6 rotated hexahedrons with 1st face as n-th face of
     #  one original hexahedron)
-    t = [1, 0, 3, 2, 5, 4]
-    m = np.array(i)[t]
-    n = np.array(j)[t]
-    p = np.array(k)[t]
-    q = np.array(l)[t]
+    m = [1, 3, 2, 0, 4, 1]
+    n = [2, 0, 3, 1, 5, 0]
+    p = [6, 4, 7, 5, 6, 3]
+    q = [5, 7, 6, 4, 7, 2]
 
     cells = np.dstack(
         (
@@ -151,10 +150,103 @@ def boundary_cells_hexahedron(mesh):
             mesh.cells[:, q],
         )
     )
-    # ensure right-hand-side cell connectivity
-    for a in [1, 3, 5]:
-        cells[:, a, :4] = cells[:, a, :4].T[::-1].T
-        cells[:, a, 4:] = cells[:, a, 4:].T[::-1].T
+
+    return cells, cells_faces
+
+
+def boundary_cells_hexahedron20(mesh):
+    "Convert the cells array of a quadratic hex mesh into a boundary cells array."
+
+    cells_hexahedron, cells_faces_hexahedron = boundary_cells_hexahedron(mesh)
+
+    # midpoints on edges of faces (boundary) of a hexahedron
+    i = [11, 9, 8, 10, 8, 12]
+    j = [19, 17, 16, 18, 9, 15]
+    k = [15, 13, 12, 14, 10, 14]
+    l = [16, 18, 17, 19, 11, 13]
+
+    cells_faces = np.dstack(
+        (
+            cells_faces_hexahedron,
+            mesh.cells[:, i],
+            mesh.cells[:, j],
+            mesh.cells[:, k],
+            mesh.cells[:, l],
+        )
+    )
+
+    # complementary faces for the creation of "boundary" hexahedrons
+    # (6 rotated hexahedrons with 1st face as n-th face of
+    #  one original hexahedron)
+    m = [9, 11, 10, 8, 12, 8]
+    n = [18, 16, 19, 17, 13, 11]
+    p = [13, 15, 14, 12, 14, 10]
+    q = [17, 19, 18, 16, 15, 9]
+
+    r = [8, 10, 9, 11, 16, 17]
+    s = [10, 8, 11, 9, 17, 16]
+    t = [14, 12, 15, 13, 18, 19]
+    u = [12, 14, 13, 15, 19, 18]
+
+    cells = np.dstack(
+        (
+            cells_hexahedron,
+            mesh.cells[:, i],
+            mesh.cells[:, j],
+            mesh.cells[:, k],
+            mesh.cells[:, l],
+            mesh.cells[:, m],
+            mesh.cells[:, n],
+            mesh.cells[:, p],
+            mesh.cells[:, q],
+            mesh.cells[:, r],
+            mesh.cells[:, s],
+            mesh.cells[:, t],
+            mesh.cells[:, u],
+        )
+    )
+
+    return cells, cells_faces
+
+
+def boundary_cells_hexahedron27(mesh):
+    "Convert the cells array of a bi-quadratic hex mesh into a boundary cells array."
+
+    cells_hexahedron20, cells_faces_hexahedron20 = boundary_cells_hexahedron20(mesh)
+
+    # midpoints of faces (boundary) of a hexahedron
+    i = [20, 21, 22, 23, 24, 25]
+
+    cells_faces = np.dstack(
+        (
+            cells_faces_hexahedron20,
+            mesh.cells[:, i],
+        )
+    )
+
+    # complementary faces for the creation of "boundary" hexahedrons
+    # (6 rotated hexahedrons with 1st face as n-th face of
+    #  one original hexahedron)
+    j = [22, 23, 21, 20, 20, 21]
+    k = [23, 22, 20, 21, 21, 20]
+    l = [24, 25, 24, 25, 22, 23]
+    m = [25, 24, 25, 24, 23, 22]
+    n = [20, 21, 22, 23, 24, 25]
+    p = [21, 20, 23, 22, 25, 24]
+    q = [26, 26, 26, 26, 26, 26]
+
+    cells = np.dstack(
+        (
+            cells_hexahedron20,
+            mesh.cells[:, j],
+            mesh.cells[:, k],
+            mesh.cells[:, l],
+            mesh.cells[:, m],
+            mesh.cells[:, n],
+            mesh.cells[:, p],
+            mesh.cells[:, q],
+        )
+    )
 
     return cells, cells_faces
 
@@ -241,6 +333,10 @@ def __init__(
             cells, cells_faces = boundary_cells_quad9(mesh)
         elif mesh.cell_type == "hexahedron":
             cells, cells_faces = boundary_cells_hexahedron(mesh)
+        elif mesh.cell_type == "hexahedron20":
+            cells, cells_faces = boundary_cells_hexahedron20(mesh)
+        elif mesh.cell_type == "hexahedron27":
+            cells, cells_faces = boundary_cells_hexahedron27(mesh)
         else:
             raise NotImplementedError("Cell type not supported.")
 
@@ -296,7 +392,11 @@ def _init_faces(self):
             dA_1 = self.dXdr[:, 0][::-1]
             dA_1[0] = -dA_1[0]
 
-        elif self.mesh.cell_type == "hexahedron":
+        elif (
+            self.mesh.cell_type == "hexahedron"
+            or self.mesh.cell_type == "hexahedron20"
+            or self.mesh.cell_type == "hexahedron27"
+        ):
             dA_1 = cross(self.dXdr[:, 0], self.dXdr[:, 1])
 
         dA = -dA_1 * self.quadrature.weights.reshape(-1, 1)

src/felupe/region/_templates.py

@@ -259,6 +259,24 @@ def __init__(self, mesh, quadrature=GaussLegendre(order=2, dim=3), grad=True):
         super().__init__(mesh, element, quadrature, grad=grad)
 
 
+class RegionQuadraticHexahedronBoundary(RegionBoundary):
+    "A boundary region with a (serendipity) quadratic hexahedron element."
+
+    def __init__(
+        self,
+        mesh,
+        quadrature=GaussLegendreBoundary(order=2, dim=3),
+        grad=True,
+        only_surface=True,
+        mask=None,
+    ):
+
+        element = QuadraticHexahedron()
+        super().__init__(
+            mesh, element, quadrature, grad=grad, only_surface=only_surface, mask=mask
+        )
+
+
 class RegionTriQuadraticHexahedron(Region):
     "A region with a tri-quadratic (lagrange) hexahedron element."
 
@@ -268,6 +286,24 @@ def __init__(self, mesh, quadrature=GaussLegendre(order=2, dim=3), grad=True):
         super().__init__(mesh, element, quadrature, grad=grad)
 
 
+class RegionTriQuadraticHexahedronBoundary(RegionBoundary):
+    "A boundary region with a tri-quadratic (lagrange) hexahedron element."
+
+    def __init__(
+        self,
+        mesh,
+        quadrature=GaussLegendreBoundary(order=2, dim=3),
+        grad=True,
+        only_surface=True,
+        mask=None,
+    ):
+
+        element = TriQuadraticHexahedron()
+        super().__init__(
+            mesh, element, quadrature, grad=grad, only_surface=only_surface, mask=mask
+        )
+
+
 class RegionLagrange(Region):
     "A region with an arbitrary order lagrange element."
 

tests/test_region.py

@@ -72,10 +72,14 @@ def test_region():
     r = fe.RegionQuadraticHexahedron(mesh2)
     f = fe.FieldsMixed(r)
 
+    r = fe.RegionQuadraticHexahedronBoundary(mesh2)
+
     mesh3 = fe.mesh.convert(mesh, 2, True, True, True)
     r = fe.RegionTriQuadraticHexahedron(mesh3)
     f = fe.FieldsMixed(r)
 
+    r = fe.RegionTriQuadraticHexahedronBoundary(mesh3)
+
     triangle = fe.Triangle()
     points = triangle.points
     cells = np.arange(3).reshape(1, -1)