diff --git a/CHANGELOG.md b/CHANGELOG.md index 383c08fc..002a9fbe 100644 --- a/CHANGELOG.md +++ b/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)`. diff --git a/docs/felupe/quadrature.rst b/docs/felupe/quadrature.rst index d5132ab1..c781370e 100644 --- a/docs/felupe/quadrature.rst +++ b/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: diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index e1b623a6..58108f9c 100644 --- a/src/felupe/__init__.py +++ b/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", diff --git a/src/felupe/quadrature/__init__.py b/src/felupe/quadrature/__init__.py index 8ca019ab..158ff9b9 100644 --- a/src/felupe/quadrature/__init__.py +++ b/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", diff --git a/src/felupe/quadrature/_gausslegendre.py b/src/felupe/quadrature/_gauss_legendre.py similarity index 100% rename from src/felupe/quadrature/_gausslegendre.py rename to src/felupe/quadrature/_gauss_legendre.py diff --git a/src/felupe/quadrature/_gauss_lobatto.py b/src/felupe/quadrature/_gauss_lobatto.py new file mode 100644 index 00000000..c1223922 --- /dev/null +++ b/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 . +""" + +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).") diff --git a/tests/test_quadrature.py b/tests/test_quadrature.py index 721127f5..1e2a9e4f 100644 --- a/tests/test_quadrature.py +++ b/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()