Quadratures without dependencies (capytaine#416)

capytaine/meshes/meshes.py

@@ -8,14 +8,14 @@
 from itertools import count
 
 import numpy as np
-from numpy.linalg import norm
 
 from capytaine.meshes.geometry import Abstract3DObject, ClippableMixin, Plane, inplace_transformation
 from capytaine.meshes.properties import compute_faces_properties
 from capytaine.meshes.surface_integrals import SurfaceIntegralsMixin
 from capytaine.meshes.quality import (merge_duplicates, heal_normals, remove_unused_vertices,
                                       heal_triangles, remove_degenerated_faces)
 from capytaine.tools.optional_imports import import_optional_dependency
+from capytaine.meshes.quadratures import compute_quadrature_on_faces
 
 LOG = logging.getLogger(__name__)
 
@@ -331,63 +331,13 @@ def quadrature_method(self):
             return None
 
     def compute_quadrature(self, method):
-        quadpy = import_optional_dependency("quadpy")
-        transform = quadpy.c2.transform
-        get_detJ = quadpy.cn._helpers.get_detJ
-
-        if method is None:
-            # No quadrature (i.e. default first order quadrature)
-            if 'quadrature' in self.__internals__:
-                del self.__internals__['quadrature']
-                del self.__internals__['quadrature_method']
-            else:
-                pass
-
-        elif isinstance(method, quadpy.c2._helpers.C2Scheme):
-            assert method.points.shape[0] == method.dim == 2
-            nb_points = method.points.shape[1]
-            points = np.empty((self.nb_faces, nb_points, 3))
-            weights = np.empty((self.nb_faces, nb_points))
-
-            self.heal_triangles()
-
-            for i_face in range(self.nb_faces):
-                # Define a local frame (Oxyz) such that
-                # * the corner A of the quadrilateral panel is the origin of the local frame
-                # * the edge AB of the quadrilateral panel is along the local x-axis,
-                # * the quadrilateral panel is within the local xy-plane (that is, its normal is along the local z-axis).
-                # Hence, the corners of the panels all have 0 as z-coordinate in the local frame.
-
-                # Coordinates in global frame
-                global_A, global_B, global_C, global_D = self.vertices[self.faces[i_face, :], :]
-                n = self.faces_normals[i_face, :]
-
-                ex = (global_B-global_A)/norm(global_B-global_A)  # unit vector of the local x-axis
-                ez = n/norm(n)                                    # unit vector of the local z-axis
-                ey = np.cross(ex, ez)                             # unit vector of the local y-axis, such that the basis is orthonormal
-
-                R = np.array([ex, ey, ez])
-                local_A = np.zeros((3,))             # coordinates of A in local frame, should be zero by construction
-                local_B = R @ (global_B - global_A)  # coordinates of B in local frame
-                local_C = R @ (global_C - global_A)  # coordinates of C in local frame
-                local_D = R @ (global_D - global_A)  # coordinates of D in local frame
-
-                local_quadrilateral = np.array([[local_A, local_D], [local_B, local_C]])[:, :, :-1]
-                # Removing last index in last dimension because not interested in z-coordinate which is 0.
-
-                local_quadpoints = transform(method.points, local_quadrilateral)
-
-                local_quadpoints_in_3d = np.concatenate([local_quadpoints, np.zeros((nb_points, 1))], axis=1)
-                global_quadpoints = np.array([R.T @ p for p in local_quadpoints_in_3d]) + global_A
-                points[i_face, :, :] = global_quadpoints
-
-                weights[i_face, :] = method.weights * 4 * np.abs(get_detJ(method.points, local_quadrilateral))
-
-            self.__internals__['quadrature'] = (points, weights)
-            self.__internals__['quadrature_method'] = method
+        self.heal_triangles()
+        all_faces = self.vertices[self.faces[:, :], :]
+        points, weights = compute_quadrature_on_faces(all_faces, method)
+        self.__internals__['quadrature'] = (points, weights)
+        self.__internals__['quadrature_method'] = method
+        return points, weights
 
-        else:
-            raise NotImplementedError
 
     ###############################
     #  Triangles and quadrangles  #

capytaine/meshes/quadratures.py

@@ -0,0 +1,80 @@
+import logging
+
+import numpy as np
+
+from capytaine.tools.optional_imports import silently_import_optional_dependency
+
+LOG = logging.getLogger(__name__)
+
+# The builtin methods are stored as a list of 2D-points in [-1, 1]² and a list
+# of corresponding weights. The 2D points will be remapped to the actual shape
+# of the faces. They are only defined for quadrilaterals. They also work for
+# triangles (although they might be subobtimal).
+
+builtin_methods = {
+        "First order": (np.array([(0.0, 0.0)]), np.array([1.0])),
+        "Gauss-Legendre 2": (
+            np.array([(+1/np.sqrt(3), +1/np.sqrt(3)),
+             (+1/np.sqrt(3), -1/np.sqrt(3)),
+             (-1/np.sqrt(3), +1/np.sqrt(3)),
+             (-1/np.sqrt(3), -1/np.sqrt(3))]),
+            np.array([1/4, 1/4, 1/4, 1/4])
+            )
+        }
+
+
+def compute_quadrature_on_faces(faces, method):
+    """
+    Compute the quadrature points and weight for numerical integration over the faces.
+
+    Parameters
+    ----------
+    faces: array of shape (nb_faces, 4, 3)
+        The 3D-coordinates of each of the 4 corners of each face.
+    method: string or quadpy object
+        The method used to compute the quadrature scheme
+
+    Returns
+    -------
+    points: array of shape (nb_faces, nb_quad_points, 3)
+        The 3D-coordinates of each of the quadrature points (their number depends on the method) of each face.
+    weights: array of shape (nb_faces, nb_quad_points)
+        Weights associated to each of the quadrature points.
+    """
+
+    if method in builtin_methods:
+        LOG.debug("Quadrature method found in builtin methods.")
+        local_points, local_weights = builtin_methods[method]
+
+    elif ((quadpy := silently_import_optional_dependency("quadpy")) is not None
+                and isinstance(method, quadpy.c2._helpers.C2Scheme)):
+        LOG.debug("Quadrature method is a Quadpy scheme: %s", method.name)
+        local_points = method.points.T
+        local_weights = method.weights
+
+    else:
+        raise ValueError(f"Unrecognized quadrature scheme: {method}.\n"
+                         f"Consider using one of the following: {set(builtin_methods.keys())}")
+
+    nb_faces = faces.shape[0]
+    nb_quad_points = len(local_weights)
+    points = np.empty((nb_faces, nb_quad_points, 3))
+    weights = np.empty((nb_faces, nb_quad_points))
+
+    for i_face in range(nb_faces):
+        for k_quad in range(nb_quad_points):
+            xk, yk = local_points[k_quad, :]
+            points[i_face, k_quad, :] = (
+                      (1+xk)*(1+yk) * faces[i_face, 0, :]
+                    + (1+xk)*(1-yk) * faces[i_face, 1, :]
+                    + (1-xk)*(1-yk) * faces[i_face, 2, :]
+                    + (1-xk)*(1+yk) * faces[i_face, 3, :]
+                    )/4
+            dxidx = ((1+yk)*faces[i_face, 0, :] + (1-yk)*faces[i_face, 1, :]
+                     - (1-yk)*faces[i_face, 2, :] - (1+yk)*faces[i_face, 3, :])/4
+            dxidy = ((1+xk)*faces[i_face, 0, :] - (1+xk)*faces[i_face, 1, :]
+                     - (1-xk)*faces[i_face, 2, :] + (1-xk)*faces[i_face, 3, :])/4
+            detJ = np.linalg.norm(np.cross(dxidx, dxidy))
+            weights[i_face, k_quad] = local_weights[k_quad] * 4 * detJ
+
+    return points, weights

docs/changelog.rst

@@ -31,6 +31,7 @@ Minor changes
 
 * Computing the RAO with :func:`cpt.post_pro.rao.rao` is not restricted to a single wave direction (or a single value of any other extra parameter) at the time anymore. (:issue:`405` and :pull:`406`)
 
+* New computation of quadrature schemes without relying on Quadpy. (:pull:`416`)
 
 Bug fixes
 ~~~~~~~~~

docs/user_manual/mesh.rst

@@ -266,28 +266,31 @@ If none of them are define, the rotation is defined around the origin of the dom
 Defining an integration quadrature
 ----------------------------------
 
-.. warning:: This feature is experimental.
-             Only quadrilaterals panels are supported at the moment.
-
 During the resolution of the BEM problem, the Green function has to be
-integrated on the mesh. By default, the integration is approximated by taking
-the value at the center of the panel and multiplying by its area. For a more
-accurate intagration, an higher order quadrature can be defined.
-
-This feature relies on the external package `quadpy` to compute the quadrature.
-You can install it with::
+integrated on each panel of the mesh. Parts of the Green function (such as the
+:math:`1/r` Rankine terms) are integrated using an exact analytical expression
+for the integral. Other parts of the Green function rely on numerical
+integration. By default, this numerical integration is done by taking the value
+at the center of the panel and multiplying by its area. For a more accurate
+intagration, an higher order quadrature can be defined.
 
-    pip install quadpy
+To define a quadrature scheme for a mesh, run the following command::
 
-Then chose one of the `available quadratures
-<https://github.com/nschloe/quadpy#quadrilateral>`_ and give it to the
-:code:`compute_quadrature` method::
+    body.mesh.compute_quadrature(method="Gauss-Legendre 2")
 
-    from quadpy.quadrilateral import stroud_c2_7_2
+The quadrature data can then be accessed at::
 
-    body.mesh.compute_quadrature(method=stroud_c2_7_2())
+    body.mesh.quadrature_points
 
-It will then be used automatically when needed.
+and will be used automatically when needed.
 
 .. warning:: Transformations of the mesh (merging, clipping, ...) may reset the quadrature.
              Compute it only on your final mesh.
+
+.. warning:: Quadratures schemes have been designed with quadrilateral panels.
+             They work on triangular panels, but might not be as optimal then.
+
+Alternatively, the :meth:`~capytaine.meshes.quadratures.compute_quadrature` method also accepts methods from the `Quadpy` package::
+
+    import quadpy
+    body.mesh.compute_quadrature(method=quadpy.c2.get_good_scheme(8))

pyproject.toml

@@ -19,7 +19,7 @@ scripts = {capytaine = "capytaine.ui.cli:main"}
 dynamic = ['version']
 
 [project.optional-dependencies]
-optional = ["matplotlib", "joblib>=1.3", "meshio", "bemio @ https://github.com/mancellin/bemio/archive/fa3a6d3d4b0862491c3723919ec8ee45cc7d055a.zip"]
+optional = ["matplotlib", "joblib>=1.3", "meshio", "quadpy-gpl", "bemio @ https://github.com/mancellin/bemio/archive/fa3a6d3d4b0862491c3723919ec8ee45cc7d055a.zip"]
 build = ["numpy", "ninja", "meson-python", "charset-normalizer"]
 test = ["pytest"]
 docs = ["sphinx", "sphinx-toolbox", "sphinxcontrib-proof", "sphinxcontrib-mermaid"]

pytest/test_bem_with_quadratures.py

@@ -1,50 +1,46 @@
-#!/usr/bin/env python
-# coding: utf-8
-
 import pytest
 import numpy as np
 import xarray as xr
 import capytaine as cpt
 
 try:
     import quadpy
+    quadpy_methods = [quadpy.c2.schemes['sommariva_05']()]
 except ImportError:
-    quadpy = None
+    quadpy_methods = []
 
+from capytaine.meshes.quadratures import builtin_methods
 
-@pytest.mark.skipif(quadpy is None, reason="quadpy is not installed")
-def test_quadrature_points_are_coplanar():
+@pytest.mark.parametrize('method', list(builtin_methods) + quadpy_methods)
+def test_quadrature_points_are_coplanar(method):
     # Check that all quadrature points are within the plane containing the panel
-    mesh = cpt.Sphere().mesh
-    for method in [quadpy.c2.schemes['sommariva_05'](), quadpy.c2.schemes["stroud_c2_7_4"]()]:
-        mesh.compute_quadrature(method)
-        for i_face in range(mesh.nb_faces):
-            A, B, C, D = mesh.vertices[mesh.faces[i_face, :], :]
-            for j_quad_point in range(mesh.quadrature_points[0].shape[1]):
-                X = mesh.quadrature_points[0][i_face, j_quad_point, :]
-                assert np.isclose(np.dot(np.cross(B-A, C-A), X-A), 0.0, atol=1e-8)
+    mesh = cpt.mesh_sphere()
+    mesh.compute_quadrature(method)
+    for i_face in range(mesh.nb_faces):
+        A, B, C, D = mesh.vertices[mesh.faces[i_face, :], :]
+        for j_quad_point in range(mesh.quadrature_points[0].shape[1]):
+            X = mesh.quadrature_points[0][i_face, j_quad_point, :]
+            assert np.isclose(np.dot(np.cross(B-A, C-A), X-A), 0.0, atol=1e-8)
 
 
-@pytest.mark.skipif(quadpy is None, reason="quadpy is not installed")
-def test_area():
+@pytest.mark.parametrize('method', list(builtin_methods) + quadpy_methods)
+def test_area(method):
     # Check that quadrature weights sum to the area of the panel
-    mesh = cpt.Sphere().mesh
-    for method in [quadpy.c2.schemes['sommariva_05'](), quadpy.c2.schemes["stroud_c2_7_4"]()]:
-        mesh.compute_quadrature(method)
-        for i_face in range(mesh.nb_faces):
-            assert np.isclose(np.sum(mesh.quadrature_points[1][i_face, :]), mesh.faces_areas[i_face], rtol=1e-2)
+    mesh = cpt.mesh_sphere()
+    mesh.compute_quadrature(method)
+    for i_face in range(mesh.nb_faces):
+        assert np.isclose(np.sum(mesh.quadrature_points[1][i_face, :]), mesh.faces_areas[i_face], rtol=1e-2)
 
 
-@pytest.mark.skipif(quadpy is None, reason="quadpy is not installed")
-def test_resolution():
-    cylinder = cpt.HorizontalCylinder(
+@pytest.mark.parametrize('method', list(builtin_methods) + quadpy_methods)
+def test_resolution(method):
+    mesh = cpt.mesh_horizontal_cylinder(
         length=5.0, radius=1.0,
-        center=(0, 0, -2),
-        reflection_symmetry=False,
-        nr=2, nx=10, ntheta=5,
-    )
-    # cylinder.show()
-    cylinder.add_translation_dof(name="Heave")
+        center=(0, 0, 0),
+        resolution=(2, 5, 10),
+    ).immersed_part()
+    body = cpt.FloatingBody(mesh=mesh)
+    body.add_translation_dof(name="Heave")
 
     test_matrix = xr.Dataset(coords={
         "omega": np.linspace(0.5, 3.0, 2),
@@ -53,12 +49,18 @@ def test_resolution():
 
     solver = cpt.BEMSolver()
 
-    cylinder.mesh.compute_quadrature(quadpy.c2.schemes['sommariva_01']())
-    data_1 = solver.fill_dataset(test_matrix, [cylinder], mesh=True)
+    data_0 = solver.fill_dataset(test_matrix, [body], mesh=True)
+
+    body.mesh.compute_quadrature(method)
+    data_1 = solver.fill_dataset(test_matrix, [body], mesh=True)
+
+    assert np.allclose(data_0["added_mass"].data, data_1["added_mass"].data, rtol=1e-1)
 
-    cylinder.mesh.compute_quadrature(quadpy.c2.schemes['sommariva_03']())
-    data_3 = solver.fill_dataset(test_matrix, [cylinder], mesh=True)
 
-    assert data_1['quadrature_method'] == "Sommariva 1"
-    assert data_3['quadrature_method'] == "Sommariva 3"
-    assert np.allclose(data_1["added_mass"].data, data_3["added_mass"].data, rtol=1e-2)
+# def test_triangular_mesh()
+#     mesh = cpt.Mesh(
+#         vertices=np.array([[0.0, 0.0, -1.0], [1.0, 0.0, -1.0], [0.0, 1.0, -1.0], [0.0, 0.0, -1.0]]),
+#         faces=np.array([[0, 1, 2, 3]])
+#     )
+#     mesh.compute_quadrature("GL2")
+#     points, weights = mesh.quadrature_points