Merge pull request #702 from adtzlr/enhance-lagrange-element

Make `ArbitraryOrderLagrangeElement` VTK-compatible

CHANGELOG.md

@@ -20,6 +20,8 @@ All notable changes to this project will be documented in this file. The format
 - Add the `opacity=0.99` argument to `MeshContainer.plot()` and `MeshContainer.screenshot()`.
 - Pass the dpi-argument to the matplotlib figure in `imshow(dpi=None)` for solids, field- and mesh-containers.
 - Permute `GaussLegendre(order=2, dim=2)` according to the points of the `BiQuadraticQuad` element by default.
+- Permute the 2- and 3-dimensional `GaussLegendre` quadrature schemes for order > 2 according to the VTK-Lagrange element formulations. That means for linear and quadratic quads and hexahedrons, the points of `GaussLegendre` are sorted according to the default VTK elements and for all higher-order elements according to the Lagrange-elements.
+- Enable default point-permutations in `RegionLagrange(permute=True)` by default.
 
 ### Fixed
 - Fix mesh-expansion with one layer `mesh.expand(n=1)`. This expands the dimension of the points-array.

src/felupe/element/_lagrange.py

@@ -25,6 +25,131 @@
 from ._base import Element
 
 
+def lagrange_quad(order):
+    "Return the cell-connectivity for an arbitrary-order Lagrange quad."
+
+    # points on a unit rectangle
+    x = np.linspace(0, 1, order + 1)
+    points = np.vstack([p.ravel() for p in np.meshgrid(x, x, indexing="ij")][::-1]).T
+
+    # search vertices
+    xmin = min(points[:, 0])
+    ymin = min(points[:, 1])
+
+    xmax = max(points[:, 0])
+    ymax = max(points[:, 1])
+
+    def search_vertice(p, x, y):
+        return np.where(np.logical_and(p[:, 0] == x, p[:, 1] == y))[0][0]
+
+    def search_edge(p, a, x):
+        return np.where(p[:, a] == x)[0][1:-1]
+
+    v1 = search_vertice(points, xmin, ymin)
+    v2 = search_vertice(points, xmax, ymin)
+    v3 = search_vertice(points, xmax, ymax)
+    v4 = search_vertice(points, xmin, ymax)
+
+    vertices = [v1, v2, v3, v4]
+
+    mask1 = np.ones_like(points[:, 0], dtype=bool)
+    mask1[vertices] = 0
+
+    e1 = search_edge(points, 1, ymin)
+    e2 = search_edge(points, 0, xmax)
+    e3 = search_edge(points, 1, ymax)
+    e4 = search_edge(points, 0, xmin)
+
+    edges = np.hstack((e1, e2, e3, e4))
+
+    mask2 = np.ones_like(points[:, 0], dtype=bool)
+    mask2[np.hstack((vertices, edges))] = 0
+
+    face = np.arange(len(points))[mask2]
+    return np.hstack((vertices, edges, face))
+
+
+def lagrange_hexahedron(order):
+    "Return the cell-connectivity for an arbitrary-order Lagrange hexahedron."
+
+    x = np.linspace(0, 1, order + 1)
+    points = np.vstack([p.ravel() for p in np.meshgrid(x, x, x, indexing="ij")][::-1]).T
+
+    # search vertices
+    xmin = min(points[:, 0])
+    ymin = min(points[:, 1])
+    zmin = min(points[:, 2])
+
+    xmax = max(points[:, 0])
+    ymax = max(points[:, 1])
+    zmax = max(points[:, 2])
+
+    def search_vertice(p, x, y, z):
+        return np.where(
+            np.logical_and(np.logical_and(p[:, 0] == x, p[:, 1] == y), p[:, 2] == z)
+        )[0][0]
+
+    def search_edge(p, a, b, x, y):
+        return np.where(np.logical_and(p[:, a] == x, p[:, b] == y))[0][1:-1]
+
+    def search_face(p, a, x, vertices, edges):
+        face = np.where(points[:, a] == x)[0]
+        mask = np.zeros_like(p[:, 0], dtype=bool)
+        mask[face] = 1
+        mask[np.hstack((vertices, edges))] = 0
+        return np.arange(len(p[:, 0]))[mask]
+
+    v1 = search_vertice(points, xmin, ymin, zmin)
+    v2 = search_vertice(points, xmax, ymin, zmin)
+    v3 = search_vertice(points, xmax, ymax, zmin)
+    v4 = search_vertice(points, xmin, ymax, zmin)
+
+    v5 = search_vertice(points, xmin, ymin, zmax)
+    v6 = search_vertice(points, xmax, ymin, zmax)
+    v7 = search_vertice(points, xmax, ymax, zmax)
+    v8 = search_vertice(points, xmin, ymax, zmax)
+
+    vertices = [v1, v2, v3, v4, v5, v6, v7, v8]
+
+    mask1 = np.ones_like(points[:, 0], dtype=bool)
+    mask1[vertices] = 0
+
+    e1 = search_edge(points, 1, 2, ymin, zmin)
+    e2 = search_edge(points, 0, 2, xmax, zmin)
+    e3 = search_edge(points, 1, 2, ymax, zmin)
+    e4 = search_edge(points, 0, 2, xmin, zmin)
+
+    e5 = search_edge(points, 1, 2, ymin, zmax)
+    e6 = search_edge(points, 0, 2, xmax, zmax)
+    e7 = search_edge(points, 1, 2, ymax, zmax)
+    e8 = search_edge(points, 0, 2, xmin, zmax)
+
+    e9 = search_edge(points, 0, 1, xmin, ymin)
+    e10 = search_edge(points, 0, 1, xmax, ymin)
+    e11 = search_edge(points, 0, 1, xmax, ymax)
+    e12 = search_edge(points, 0, 1, xmin, ymax)
+
+    edges = np.hstack((e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12))
+
+    mask2 = np.ones_like(points[:, 0], dtype=bool)
+    mask2[np.hstack((vertices, edges))] = 0
+
+    f1 = search_face(points, 0, xmin, vertices, edges)
+    f2 = search_face(points, 0, xmax, vertices, edges)
+    f3 = search_face(points, 1, ymin, vertices, edges)
+    f4 = search_face(points, 1, ymax, vertices, edges)
+    f5 = search_face(points, 2, zmin, vertices, edges)
+    f6 = search_face(points, 2, zmax, vertices, edges)
+
+    faces = np.hstack((f1, f2, f3, f4, f5, f6))
+
+    mask3 = np.ones_like(points[:, 0], dtype=bool)
+    mask3[np.hstack((vertices, edges, faces))] = 0
+    volume = np.arange(len(points))[mask3]
+
+    return np.hstack((vertices, edges, faces, volume))
+
+
 class ArbitraryOrderLagrange(Element):
     r"""A n-dimensional Lagrange finite element (e.g. line, quad or hexahdron) of
     arbitrary order.
@@ -103,13 +228,17 @@ class ArbitraryOrderLagrange(Element):
             \otimes \boldsymbol{h}(t)
     """
 
-    def __init__(self, order, dim, interval=(-1, 1)):
+    def __init__(self, order, dim, interval=(-1, 1), permute=True):
         self._order = order
         self._nshape = order + 1
         self._npoints = self._nshape**dim
         self._nbasis = self._npoints
         self._interval = interval
 
+        self.permute = None
+        if permute:
+            self.permute = [None, None, lagrange_quad, lagrange_hexahedron][dim](order)
+
         super().__init__(shape=(self._npoints, dim))
 
         # init curve-parameter matrix
@@ -128,14 +257,21 @@ def __init__(self, order, dim, interval=(-1, 1)):
         grid = np.meshgrid(*np.tile(self._points(n), (dim, 1)), indexing="ij")[::-1]
         self.points = np.vstack([point.ravel() for point in grid]).T
 
+        if self.permute is not None:
+            self.points = self.points[self.permute]
+
     def function(self, r):
         "Return the shape functions at given coordinate vector r."
         n = self._nshape
 
         # 1d - basis function vectors per axis
-        h = [self._AT @ self._polynomial(ra, n) for ra in r]
+        fun = [self._AT @ self._polynomial(ra, n) for ra in r]
+        h = np.einsum(self._subscripts, *fun).ravel("F")
 
-        return np.einsum(self._subscripts, *h).ravel("F")
+        if self.permute is not None:
+            h = h[self.permute]
+
+        return h
 
     def gradient(self, r):
         "Return the gradient of shape functions at given coordinate vector r."
@@ -156,6 +292,9 @@ def gradient(self, r):
             g[i] = k[i]
             dhdr[:, i] = np.einsum(self._subscripts, *g).ravel("F")
 
+        if self.permute is not None:
+            dhdr = dhdr[self.permute]
+
         return dhdr
 
     def _points(self, n):

src/felupe/mesh/_geometry.py

@@ -18,6 +18,7 @@
 
 import numpy as np
 
+from ..element._lagrange import lagrange_hexahedron, lagrange_quad
 from ._line_rectangle_cube import cube_hexa, line_line, rectangle_quad
 from ._mesh import Mesh
 from ._tools import concatenate
@@ -281,153 +282,31 @@ def __init__(self, *xi, indexing="ij", **kwargs):
         )
 
 
-class RectangleArbitraryOrderQuad(Mesh):
-    """A rectangular 2d-mesh with arbitrary order quads between ``a`` and ``b`` with
-    ``n`` points per axis."""
+class RectangleArbitraryOrderQuad(Rectangle):
+    """A rectangular 2d-mesh with an arbitrarr-order Lagrange quad between ``a`` and
+    ``b``.
+    """
 
     def __init__(self, a=(0.0, 0.0), b=(1.0, 1.0), order=2):
-        yv, xv = np.meshgrid(
-            np.linspace(a[1], b[1], order + 1),
-            np.linspace(a[0], b[0], order + 1),
-            indexing="ij",
+        super().__init__(a=a, b=b, n=order + 1)
+        self.update(
+            cells=lagrange_quad(order=order).reshape(1, -1),
+            cell_type="VTK_LAGRANGE_QUADRILATERAL",
         )
 
-        points = np.vstack((xv.flatten(), yv.flatten())).T
-
-        # search vertices
-        xmin = min(points[:, 0])
-        ymin = min(points[:, 1])
-
-        xmax = max(points[:, 0])
-        ymax = max(points[:, 1])
-
-        def search_vertice(p, x, y):
-            return np.where(np.logical_and(p[:, 0] == x, p[:, 1] == y))[0][0]
-
-        def search_edge(p, a, x):
-            return np.where(p[:, a] == x)[0][1:-1]
-
-        v1 = search_vertice(points, xmin, ymin)
-        v2 = search_vertice(points, xmax, ymin)
-        v3 = search_vertice(points, xmax, ymax)
-        v4 = search_vertice(points, xmin, ymax)
-
-        vertices = [v1, v2, v3, v4]
-
-        mask1 = np.ones_like(points[:, 0], dtype=bool)
-        mask1[vertices] = 0
-        # points_no_verts = points[mask1]
-
-        e1 = search_edge(points, 1, ymin)
-        e2 = search_edge(points, 0, xmax)
-        e3 = search_edge(points, 1, ymax)
-        e4 = search_edge(points, 0, xmin)
 
-        edges = np.hstack((e1, e2, e3, e4))
-
-        mask2 = np.ones_like(points[:, 0], dtype=bool)
-        mask2[np.hstack((vertices, edges))] = 0
-        # points_no_verts_edges = points[mask2]
-
-        face = np.arange(len(points))[mask2]
-
-        cells = np.hstack((vertices, edges, face)).reshape(1, -1)
-
-        super().__init__(points, cells, cell_type="VTK_LAGRANGE_QUADRILATERAL")
-
-
-class CubeArbitraryOrderHexahedron(Mesh):
-    """A cube shaped 3d-mesh with arbitrary order hexahedrons between ``a`` and ``b``
-    with ``n`` points per axis."""
+class CubeArbitraryOrderHexahedron(Cube):
+    """A rectangular 2d-mesh with an arbitrarr-order Lagrange hexahedron between ``a``
+    and ``b``.
+    """
 
     def __init__(self, a=(0.0, 0.0, 0.0), b=(1.0, 1.0, 1.0), order=2):
-        zv, yv, xv = np.meshgrid(
-            np.linspace(a[2], b[2], order + 1),
-            np.linspace(a[1], b[1], order + 1),
-            np.linspace(a[0], b[0], order + 1),
-            indexing="ij",
+        super().__init__(a=a, b=b, n=order + 1)
+        self.update(
+            cells=lagrange_hexahedron(order=order).reshape(1, -1),
+            cell_type="VTK_LAGRANGE_HEXAHEDRON",
         )
 
-        points = np.vstack((xv.flatten(), yv.flatten(), zv.flatten())).T
-
-        # search vertices
-        xmin = min(points[:, 0])
-        ymin = min(points[:, 1])
-        zmin = min(points[:, 2])
-
-        xmax = max(points[:, 0])
-        ymax = max(points[:, 1])
-        zmax = max(points[:, 2])
-
-        def search_vertice(p, x, y, z):
-            return np.where(
-                np.logical_and(np.logical_and(p[:, 0] == x, p[:, 1] == y), p[:, 2] == z)
-            )[0][0]
-
-        def search_edge(p, a, b, x, y):
-            return np.where(np.logical_and(p[:, a] == x, p[:, b] == y))[0][1:-1]
-
-        def search_face(p, a, x, vertices, edges):
-            face = np.where(points[:, a] == x)[0]
-            mask = np.zeros_like(p[:, 0], dtype=bool)
-            mask[face] = 1
-            mask[np.hstack((vertices, edges))] = 0
-            return np.arange(len(p[:, 0]))[mask]
-
-        v1 = search_vertice(points, xmin, ymin, zmin)
-        v2 = search_vertice(points, xmax, ymin, zmin)
-        v3 = search_vertice(points, xmax, ymax, zmin)
-        v4 = search_vertice(points, xmin, ymax, zmin)
-
-        v5 = search_vertice(points, xmin, ymin, zmax)
-        v6 = search_vertice(points, xmax, ymin, zmax)
-        v7 = search_vertice(points, xmax, ymax, zmax)
-        v8 = search_vertice(points, xmin, ymax, zmax)
-
-        vertices = [v1, v2, v3, v4, v5, v6, v7, v8]
-
-        mask1 = np.ones_like(points[:, 0], dtype=bool)
-        mask1[vertices] = 0
-        # points_no_verts = points[mask1]
-
-        e1 = search_edge(points, 1, 2, ymin, zmin)
-        e2 = search_edge(points, 0, 2, xmax, zmin)
-        e3 = search_edge(points, 1, 2, ymax, zmin)
-        e4 = search_edge(points, 0, 2, xmin, zmin)
-
-        e5 = search_edge(points, 1, 2, ymin, zmax)
-        e6 = search_edge(points, 0, 2, xmax, zmax)
-        e7 = search_edge(points, 1, 2, ymax, zmax)
-        e8 = search_edge(points, 0, 2, xmin, zmax)
-
-        e9 = search_edge(points, 0, 1, xmin, ymin)
-        e10 = search_edge(points, 0, 1, xmax, ymin)
-        e11 = search_edge(points, 0, 1, xmax, ymax)
-        e12 = search_edge(points, 0, 1, xmin, ymax)
-
-        edges = np.hstack((e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12))
-
-        mask2 = np.ones_like(points[:, 0], dtype=bool)
-        mask2[np.hstack((vertices, edges))] = 0
-        # points_no_verts_edges = points[mask2]
-
-        f1 = search_face(points, 0, xmin, vertices, edges)
-        f2 = search_face(points, 0, xmax, vertices, edges)
-        f3 = search_face(points, 1, ymin, vertices, edges)
-        f4 = search_face(points, 1, ymax, vertices, edges)
-        f5 = search_face(points, 2, zmin, vertices, edges)
-        f6 = search_face(points, 2, zmax, vertices, edges)
-
-        faces = np.hstack((f1, f2, f3, f4, f5, f6))
-
-        mask3 = np.ones_like(points[:, 0], dtype=bool)
-        mask3[np.hstack((vertices, edges, faces))] = 0
-        volume = np.arange(len(points))[mask3]
-
-        cells = np.hstack((vertices, edges, faces, volume)).reshape(1, -1)
-
-        super().__init__(points, cells, cell_type="VTK_LAGRANGE_HEXAHEDRON")
-
 
 class Circle(Mesh):
     """A circular shaped 2d-mesh with quads and ``n`` points on the circumferential