Add GaussLobatto quadrature scheme (#943)

* Add `GaussLobatto` quadrature

* Update test_quadrature.py

* Update test_quadrature.py

* Update _gauss_lobatto.py

CHANGELOG.md

@@ -17,6 +17,7 @@ All notable changes to this project will be documented in this file. The format
 - Add a vertex-region `RegionVertex`.
 - Add `MeshContainer.as_vertex_mesh()` to create a merged vertex mesh from the meshes of the mesh container.
 - Add `Field.from_mesh_container(container)` to create a top-level field on a vertex mesh.
+- Add `GaussLobatto` quadrature. The order-argument is equal to the 1d-sample points minus two. This ensures identical minimum order-arguments for Gauss-Legendre and Gauss-Lobatto schemes.
 
 ### Changed
 - The first Piola-Kirchhoff stress tensor is evaluated if `ViewSolid(stress_type=None)`.

docs/felupe/quadrature.rst

@@ -11,6 +11,8 @@ This module contains quadrature (numeric integration) schemes for different fini
 
    GaussLegendre
    GaussLegendreBoundary
+   GaussLobatto
+   GaussLobattoBoundary
 
 **Triangles and Tetrahedrons**
 
@@ -44,6 +46,16 @@ This module contains quadrature (numeric integration) schemes for different fini
    :undoc-members:
    :show-inheritance:
 
+.. autoclass:: felupe.GaussLobatto
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+.. autoclass:: felupe.GaussLobattoBoundary
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 .. autoclass:: felupe.quadrature.Triangle
    :members:
    :undoc-members:

src/felupe/__init__.py

@@ -86,7 +86,13 @@
     Step,
 )
 from .mesh import Circle, Cube, Grid, Mesh, MeshContainer, Point, Rectangle
-from .quadrature import BazantOh, GaussLegendre, GaussLegendreBoundary
+from .quadrature import (
+    BazantOh,
+    GaussLegendre,
+    GaussLegendreBoundary,
+    GaussLobatto,
+    GaussLobattoBoundary,
+)
 from .quadrature import Tetrahedron as TetrahedronQuadrature
 from .quadrature import Triangle as TriangleQuadrature
 from .region import (
@@ -211,6 +217,8 @@
     "MultiPointContact",
     "GaussLegendre",
     "GaussLegendreBoundary",
+    "GaussLobatto",
+    "GaussLobattoBoundary",
     "TetrahedronQuadrature",
     "TriangleQuadrature",
     "BazantOh",

src/felupe/quadrature/__init__.py

@@ -1,4 +1,5 @@
-from ._gausslegendre import GaussLegendre, GaussLegendreBoundary
+from ._gauss_legendre import GaussLegendre, GaussLegendreBoundary
+from ._gauss_lobatto import GaussLobatto, GaussLobattoBoundary
 from ._scheme import Scheme
 from ._sphere import BazantOh
 from ._tetra import Tetrahedron
@@ -8,6 +9,8 @@
     "Scheme",
     "GaussLegendre",
     "GaussLegendreBoundary",
+    "GaussLobatto",
+    "GaussLobattoBoundary",
     "Tetrahedron",
     "Triangle",
     "BazantOh",

src/felupe/quadrature/_gauss_lobatto.py

@@ -0,0 +1,171 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from string import ascii_lowercase
+
+import numpy as np
+
+from ._scheme import Scheme
+
+
+def gauss_lobatto(deg):
+    r"""Gauss-Lobatto quadrature.
+
+    Computes the sample points and weights for Gauss-Lobatto quadrature. These sample
+    points and weights will correctly integrate polynomials of degree
+    :math:`2 \cdot deg - 3` or less over the interval :math:`[-1, 1]` with the weight
+    function :math:`f(x) = 1`.
+
+    Parameters
+    ----------
+    deg : int
+        Number of sample points and weights. It must be >= 2.
+
+    Returns
+    -------
+    x : ndarray
+        1-D ndarray containing the sample points.
+    y : ndarray
+        1-D ndarray containing the weights.
+    """
+    if deg == 2:
+        x = np.array([-1.0, 1.0])
+        y = np.array([1.0, 1.0])
+
+    elif deg == 3:
+        x = np.array([-1.0, 0.0, 1.0])
+        y = np.array([1.0, 4.0, 1.0]) / 3
+
+    elif deg == 4:
+        a = np.sqrt(0.2)
+        x = np.array([-1.0, -a, a, 1.0])
+        y = np.array([1.0, 5.0, 5.0, 1.0]) / 6
+
+    elif deg == 5:
+        a = np.sqrt(3 / 7)
+        x = np.array([-1.0, -a, 0.0, a, 1.0])
+        y = np.array([0.1, 49 / 90, 32 / 45, 49 / 90, 0.1])
+
+    elif deg == 6:
+        a = np.sqrt(1 / 3 - 2 * np.sqrt(7) / 21)
+        b = np.sqrt(1 / 3 + 2 * np.sqrt(7) / 21)
+        c = (14 + np.sqrt(7)) / 30
+        d = (14 - np.sqrt(7)) / 30
+        x = np.array([-1.0, -b, -a, a, b, 1.0])
+        y = np.array([1 / 15, d, c, c, d, 1 / 15])
+
+    elif deg == 7:
+        a = np.sqrt(5 / 11 - 2 / 11 * np.sqrt(5 / 3))
+        b = np.sqrt(5 / 11 + 2 / 11 * np.sqrt(5 / 3))
+        c = (124 + 7 * np.sqrt(15)) / 350
+        d = (124 - 7 * np.sqrt(15)) / 350
+        x = np.array([-1.0, -b, -a, 0, a, b, 1.0])
+        y = np.array([1 / 21, d, c, 256 / 525, c, d, 1 / 21])
+
+    else:
+        raise ValueError("deg must be a positive integer (2 <= deg <= 7)")
+
+    return x, y
+
+
+class GaussLobatto(Scheme):
+    r"""An arbitrary-`order` Gauss-Lobatto quadrature rule of dimension 1, 2 or 3 on
+    the interval :math:`[-1, 1]`.
+
+    Parameters
+    ----------
+    order : int
+        The number of sample points :math:`n` minus two. The quadrature rule integrates
+        degree :math:`2n-3` polynomials exactly.
+    dim : int
+        The dimension of the quadrature region.
+    permute : bool, optional
+        Permute the quadrature points according to the cell point orderings (default is
+        True). This is supported for two and three dimensions as well as first and
+        second order schemes. Otherwise this flag is silently ignored.
+
+    Notes
+    -----
+
+    The approximation is given by
+
+    ..  math::
+
+        \int_{-1}^1 f(x) dx \approx \sum f(x_q) w_q
+
+    with quadrature points :math:`x_q` and corresponding weights :math:`w_q`.
+
+    """
+
+    def __init__(self, order: int, dim: int):
+        if dim not in [1, 2, 3]:
+            raise ValueError("Wrong dimension.")
+
+        x, w = gauss_lobatto(2 + order)
+
+        points = (
+            np.stack(np.meshgrid(*([x] * dim), indexing="ij"))[::-1].reshape(dim, -1).T
+        )
+
+        idx = list(ascii_lowercase)[:dim]
+        weights = np.einsum(", ".join(idx), *([w] * dim)).ravel()
+
+        super().__init__(points, weights)
+
+
+class GaussLobattoBoundary(GaussLobatto):
+    r"""An arbitrary-`order` Gauss-Lobatto quadrature rule of `dim` 1, 2 or 3 on the
+    interval ``[-1, 1]``.
+
+    Parameters
+    ----------
+    order : int
+        The number of sample points :math:`n` minus two. The quadrature rule integrates
+        degree :math:`2n-3` polynomials exactly.
+    dim : int
+        The dimension of the quadrature region.
+    permute : bool, optional
+        Permute the quadrature points according to the cell point orderings (default is
+        True).
+
+    Notes
+    -----
+
+    The approximation is given by
+
+    ..  math::
+
+        \int_{-1}^1 f(x) dx \approx \sum f(x_q) w_q
+
+    with quadrature points :math:`x_q` and corresponding weights :math:`w_q`.
+
+    """
+
+    def __init__(self, order: int, dim: int):
+        super().__init__(order=order, dim=dim - 1)
+
+        # reset dimension
+        self.dim = dim
+
+        if self.dim == 2 or self.dim == 3:
+            # quadrature points projected onto first edge of a quad
+            #                          or onto first face of a hexahedron
+            self.points = np.hstack((self.points, -np.ones((len(self.points), 1))))
+
+        else:
+            raise ValueError("Given dimension not implemented (must be 2 or 3).")

tests/test_quadrature.py

@@ -1,11 +1,11 @@
 # -*- coding: utf-8 -*-
 """
- _______  _______  ___      __   __  _______  _______ 
+ _______  _______  ___      __   __  _______  _______
 |       ||       ||   |    |  | |  ||       ||       |
 |    ___||    ___||   |    |  | |  ||    _  ||    ___|
-|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___ 
+|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___
 |    ___||    ___||   |___ |       ||    ___||    ___|
-|   |    |   |___ |       ||       ||   |    |   |___ 
+|   |    |   |___ |       ||       ||   |    |   |___
 |___|    |_______||_______||_______||___|    |_______|
 
 This file is part of felupe.
@@ -86,6 +86,30 @@ def test_gausslegendre_boundary():
         fem.GaussLegendreBoundary(order=1, dim=4)
 
 
+def test_gausslobatto():
+    for order in [0, 1, 2, 3, 4, 5]:
+        for dim in [1, 2, 3]:
+            fem.GaussLobatto(order=order, dim=dim)
+
+    with pytest.raises(ValueError):
+        fem.GaussLobatto(order=1, dim=4)
+
+    with pytest.raises(ValueError):
+        fem.GaussLobatto(order=9, dim=4)
+
+
+def test_gausslobatto_boundary():
+    for order in [0, 1, 2, 3, 4, 5]:
+        for dim in [2, 3]:
+            fem.GaussLobattoBoundary(order=order, dim=dim)
+
+    with pytest.raises(ValueError):
+        fem.GaussLobattoBoundary(order=1, dim=4)
+
+    with pytest.raises(ValueError):
+        fem.GaussLobattoBoundary(order=9, dim=4)
+
+
 def test_triangle():
     q = fem.TriangleQuadrature(order=1)
     assert q.points.shape == (1, 2)
@@ -129,6 +153,8 @@ def test_sphere():
 if __name__ == "__main__":
     test_gausslegendre()
     test_gausslegendre_boundary()
+    test_gausslobatto_boundary()
+    test_gausslobatto()
     test_triangle()
     test_tetra()
     test_sphere()