From ace89e29c2b78f34bf80d75d0e67b0c2d7fa2253 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 15 Mar 2024 14:05:41 +0100 Subject: [PATCH 1/7] Make `ArbitraryOrderLagrangeElement` VTK-compatible --- src/felupe/element/_lagrange.py | 145 +++++++++++++++++++++++++++++- src/felupe/mesh/_geometry.py | 155 ++++---------------------------- src/felupe/region/_templates.py | 4 +- 3 files changed, 161 insertions(+), 143 deletions(-) diff --git a/src/felupe/element/_lagrange.py b/src/felupe/element/_lagrange.py index 97a90516..47d75d6d 100644 --- a/src/felupe/element/_lagrange.py +++ b/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")]).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")]).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): diff --git a/src/felupe/mesh/_geometry.py b/src/felupe/mesh/_geometry.py index d79df8b1..44bd4de1 100644 --- a/src/felupe/mesh/_geometry.py +++ b/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 diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index 054a2828..5f6e70fa 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -388,11 +388,11 @@ def __init__( class RegionLagrange(Region): "A region with an arbitrary order lagrange element." - def __init__(self, mesh, order, dim, quadrature=None, grad=True): + def __init__(self, mesh, order, dim, quadrature=None, grad=True, permute=False): if quadrature is None: quadrature = GaussLegendre(order=order, dim=dim, permute=False) - element = ArbitraryOrderLagrange(order, dim) + element = ArbitraryOrderLagrange(order, dim, permute=permute) self.order = order super().__init__(mesh, element, quadrature, grad=grad) From 72b4738747d0951b0c7c7304861dd70f983fc65e Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Fri, 15 Mar 2024 14:07:35 +0100 Subject: [PATCH 2/7] Permute lagrange-based regions by default --- src/felupe/region/_templates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index 5f6e70fa..c25f970f 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -388,7 +388,7 @@ def __init__( class RegionLagrange(Region): "A region with an arbitrary order lagrange element." - def __init__(self, mesh, order, dim, quadrature=None, grad=True, permute=False): + def __init__(self, mesh, order, dim, quadrature=None, grad=True, permute=True): if quadrature is None: quadrature = GaussLegendre(order=order, dim=dim, permute=False) From 583a1a44b276433dc561dd5b2f6161a75415b406 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 16 Mar 2024 14:16:07 +0100 Subject: [PATCH 3/7] Fix lagrange-elements for VTK usage --- src/felupe/element/_lagrange.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/felupe/element/_lagrange.py b/src/felupe/element/_lagrange.py index 47d75d6d..d2a2d534 100644 --- a/src/felupe/element/_lagrange.py +++ b/src/felupe/element/_lagrange.py @@ -30,7 +30,7 @@ def lagrange_quad(order): # 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")]).T + points = np.vstack([p.ravel() for p in np.meshgrid(x, x, indexing="ij")][::-1]).T # search vertices xmin = min(points[:, 0]) @@ -73,7 +73,7 @@ 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")]).T + points = np.vstack([p.ravel() for p in np.meshgrid(x, x, x, indexing="ij")][::-1]).T # search vertices xmin = min(points[:, 0]) From 19149b1b7a64b36b4536437d5a955c1a9f46e123 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 16 Mar 2024 14:23:32 +0100 Subject: [PATCH 4/7] Permute `GaussLegendre` by default for all higher order lagrange elements too --- src/felupe/quadrature/_gausslegendre.py | 30 +++++++++++++++---------- 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/felupe/quadrature/_gausslegendre.py b/src/felupe/quadrature/_gausslegendre.py index 9c5eb861..c1befea4 100644 --- a/src/felupe/quadrature/_gausslegendre.py +++ b/src/felupe/quadrature/_gausslegendre.py @@ -20,6 +20,7 @@ import numpy as np +from ..element._lagrange import lagrange_hexahedron, lagrange_quad from ._scheme import Scheme @@ -72,29 +73,34 @@ def __init__(self, order: int, dim: int, permute: bool = True): idx = list(ascii_lowercase)[:dim] weights = np.einsum(", ".join(idx), *([w] * dim)).ravel() + sort = None if permute and order == 1 and dim == 2: - points = points[[0, 1, 3, 2]] - weights = weights[[0, 1, 3, 2]] + sort = [0, 1, 3, 2] - if permute and order == 1 and dim == 3: - points = points[[0, 1, 3, 2, 4, 5, 7, 6]] - weights = weights[[0, 1, 3, 2, 4, 5, 7, 6]] + elif permute and order == 1 and dim == 3: + sort = [0, 1, 3, 2, 4, 5, 7, 6] - if permute and order == 2 and dim == 2: - points = points[[0, 2, 8, 6, 1, 5, 7, 3, 4]] - weights = weights[[0, 2, 8, 6, 1, 5, 7, 3, 4]] + elif permute and order == 2 and dim == 2: + sort = [0, 2, 8, 6, 1, 5, 7, 3, 4] - if permute and order == 2 and dim == 3: + elif permute and order == 2 and dim == 3: vertices = np.array([0, 2, 8, 6, 18, 20, 26, 24]) edges = np.array([1, 5, 7, 3, 19, 23, 25, 21, 9, 11, 17, 15]) faces = np.array([12, 14, 10, 16, 4, 22]) volume = np.array([13]) - permute = np.concatenate((vertices, edges, faces, volume)) + sort = np.concatenate((vertices, edges, faces, volume)) - points = points[permute] - weights = weights[permute] + elif permute and order > 2 and dim == 2: + sort = lagrange_quad(order=order) + + elif permute and order > 2 and dim == 3: + sort = lagrange_hexahedron(order=order) + + if sort is not None: + points = points[sort] + weights = weights[sort] super().__init__(points, weights) From 79b8fe3c2c463ebb21e93b48146aa9bd042008d9 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 16 Mar 2024 14:26:19 +0100 Subject: [PATCH 5/7] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7398fd41..0e653b8c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ 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. ### Fixed - Fix mesh-expansion with one layer `mesh.expand(n=1)`. This expands the dimension of the points-array. From f6a7bdb080354e433dc6845dfc5657e31314b976 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 16 Mar 2024 14:34:32 +0100 Subject: [PATCH 6/7] Test new features (Lagrange element/region) --- src/felupe/region/_templates.py | 2 +- tests/test_region.py | 13 +++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/felupe/region/_templates.py b/src/felupe/region/_templates.py index c25f970f..b82eb802 100644 --- a/src/felupe/region/_templates.py +++ b/src/felupe/region/_templates.py @@ -390,7 +390,7 @@ class RegionLagrange(Region): def __init__(self, mesh, order, dim, quadrature=None, grad=True, permute=True): if quadrature is None: - quadrature = GaussLegendre(order=order, dim=dim, permute=False) + quadrature = GaussLegendre(order=order, dim=dim, permute=permute) element = ArbitraryOrderLagrange(order, dim, permute=permute) self.order = order diff --git a/tests/test_region.py b/tests/test_region.py index c69825bc..e082b06e 100644 --- a/tests/test_region.py +++ b/tests/test_region.py @@ -131,7 +131,20 @@ def test_region(): points = np.vstack((x.ravel(), y.ravel(), z.ravel())).T cells = np.arange(len(points)).reshape(1, -1) mesh = fem.Mesh(points, cells, "lagrange") + r = fem.RegionLagrange(mesh, order, dim, permute=False) + assert not np.any(r.dV <= 0) + + order = 5 + dim = 2 + mesh = fem.mesh.RectangleArbitraryOrderQuad(order=order) + r = fem.RegionLagrange(mesh, order, dim) + assert not np.any(r.dV <= 0) + + order = 5 + dim = 3 + mesh = fem.mesh.CubeArbitraryOrderHexahedron(order=order) r = fem.RegionLagrange(mesh, order, dim) + assert not np.any(r.dV <= 0) if __name__ == "__main__": From 5335ca36402e9916dba2e6ffe08e3b156740cc32 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 16 Mar 2024 14:36:30 +0100 Subject: [PATCH 7/7] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0e653b8c..4cdc25fa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ All notable changes to this project will be documented in this file. The format - 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.