From b8ada4407d8cba62d403e772e67cf73878060e2f Mon Sep 17 00:00:00 2001 From: Matthieu Ancellin <31126826+mancellin@users.noreply.github.com> Date: Tue, 24 Oct 2023 15:09:20 +0200 Subject: [PATCH] Quadratures without dependencies (#416) --- capytaine/meshes/meshes.py | 64 +++-------------------- capytaine/meshes/quadratures.py | 80 +++++++++++++++++++++++++++++ docs/changelog.rst | 1 + docs/user_manual/mesh.rst | 35 +++++++------ pyproject.toml | 2 +- pytest/test_bem_with_quadratures.py | 76 ++++++++++++++------------- 6 files changed, 147 insertions(+), 111 deletions(-) create mode 100644 capytaine/meshes/quadratures.py diff --git a/capytaine/meshes/meshes.py b/capytaine/meshes/meshes.py index 7f1a0bfe..325613d6 100644 --- a/capytaine/meshes/meshes.py +++ b/capytaine/meshes/meshes.py @@ -8,7 +8,6 @@ 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 @@ -16,6 +15,7 @@ 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 # diff --git a/capytaine/meshes/quadratures.py b/capytaine/meshes/quadratures.py new file mode 100644 index 00000000..4498cf41 --- /dev/null +++ b/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 diff --git a/docs/changelog.rst b/docs/changelog.rst index 61662a13..f74a4516 100644 --- a/docs/changelog.rst +++ b/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 ~~~~~~~~~ diff --git a/docs/user_manual/mesh.rst b/docs/user_manual/mesh.rst index 6bffc821..d1dee5ec 100644 --- a/docs/user_manual/mesh.rst +++ b/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 -`_ 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)) diff --git a/pyproject.toml b/pyproject.toml index 7d85372d..ba21878f 100644 --- a/pyproject.toml +++ b/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"] diff --git a/pytest/test_bem_with_quadratures.py b/pytest/test_bem_with_quadratures.py index e83a5a4a..6daf5142 100644 --- a/pytest/test_bem_with_quadratures.py +++ b/pytest/test_bem_with_quadratures.py @@ -1,6 +1,3 @@ -#!/usr/bin/env python -# coding: utf-8 - import pytest import numpy as np import xarray as xr @@ -8,43 +5,42 @@ 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