diff --git a/capytaine/meshes/meshes.py b/capytaine/meshes/meshes.py
index 325613d6..d1341b92 100644
--- a/capytaine/meshes/meshes.py
+++ b/capytaine/meshes/meshes.py
@@ -34,11 +34,14 @@ class Mesh(ClippableMixin, SurfaceIntegralsMixin, Abstract3DObject):
         description.
     name : str, optional
         The name of the mesh. If None, the mesh is given an automatic name based on its internal ID.
+    quadrature_method: None or str or Quadpy quadrature, optional
+        The method used to compute quadrature points in each cells.
+        By default: None, that is a one-point first order scheme is used.
     """
 
     _ids = count(0)  # A counter for automatic naming of new meshes.
 
-    def __init__(self, vertices=None, faces=None, name=None):
+    def __init__(self, vertices=None, faces=None, name=None, *, quadrature_method=None):
 
         if vertices is None or len(vertices) == 0:
             vertices = np.zeros((0, 3))
@@ -57,6 +60,8 @@ def __init__(self, vertices=None, faces=None, name=None):
 
         LOG.debug(f"New mesh: {repr(self)}")
 
+        self.quadrature_method = quadrature_method
+
     def __short_str__(self):
         return (f"{self.__class__.__name__}(..., name=\"{self.name}\")")
 
@@ -314,28 +319,20 @@ def faces_radiuses(self) -> np.ndarray:
 
     @property
     def quadrature_points(self):
-        if 'quadrature' in self.__internals__:
-            return self.__internals__['quadrature']
-        else:
-            # Default: first order quadrature
-            return (
-                self.faces_centers.reshape((self.nb_faces, 1, 3)),  # Points
-                self.faces_areas.reshape((self.nb_faces, 1))        # Weights
-            )
-
-    @property
-    def quadrature_method(self):
-        if 'quadrature_method' in self.__internals__:
-            return self.__internals__['quadrature_method']
-        else:
-            return None
+        if 'quadrature' not in self.__internals__:
+            self.compute_quadrature(self.quadrature_method)
+        return self.__internals__['quadrature']
 
     def compute_quadrature(self, method):
         self.heal_triangles()
         all_faces = self.vertices[self.faces[:, :], :]
-        points, weights = compute_quadrature_on_faces(all_faces, method)
+        if method is None:
+            points = self.faces_centers.reshape((self.nb_faces, 1, 3))
+            weights = self.faces_areas.reshape((self.nb_faces, 1))
+        else:
+            points, weights = compute_quadrature_on_faces(all_faces, method)
         self.__internals__['quadrature'] = (points, weights)
-        self.__internals__['quadrature_method'] = method
+        self.quadrature_method = method
         return points, weights
 
 
diff --git a/docs/changelog.rst b/docs/changelog.rst
index aafcdde7..1cec4bf9 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -52,6 +52,8 @@ Bug fixes
 
 * Fix error message when computing hydrostatic stiffness of non-neutrally-buoyant body that is not a single rigid body. (:issue:`413` and :pull:`414`)
 
+* Fix bug causing the quadrature method of a mesh to be forgotten when the mesh was put in a body. ``quadrature_method`` can now be passed as argument when initializing a new mesh. (:pull:`417`)
+
 Internals
 ~~~~~~~~~
 
diff --git a/pytest/test_bem_engines.py b/pytest/test_bem_engines.py
index f4483c8d..c4680ae4 100644
--- a/pytest/test_bem_engines.py
+++ b/pytest/test_bem_engines.py
@@ -87,7 +87,7 @@ def test_gradient_G_with_collocation_points():
     mesh = cpt.mesh_sphere(radius=1, center=(0, 0, 0), resolution=(10, 10)).immersed_part()
     _, gradG_1 = cpt.Delhommeau().evaluate(mesh.faces_centers, mesh, 0.0, np.infty, 1.0, early_dot_product=False)
     _, gradG_2 = cpt.Delhommeau().evaluate(mesh.copy(), mesh, 0.0, np.infty, 1.0, early_dot_product=False)
-    np.testing.assert_allclose(gradG_1, gradG_2)
+    np.testing.assert_allclose(gradG_1, gradG_2, atol=1e-10, rtol=1e-10)
 
 
 def test_gradient_G_diagonal_term():
@@ -101,7 +101,7 @@ def test_gradient_G_diagonal_term():
     for i in range(mesh.nb_faces):
         diag_normal[i, i, :] = mesh.faces_normals[i, :]
 
-    np.testing.assert_allclose(gradG_1, gradG_2 - 0.5*diag_normal)
+    np.testing.assert_allclose(gradG_1, gradG_2 - 0.5*diag_normal, atol=1e-10, rtol=1e-10)
 
 
 def test_a_posteriori_scalar_product_direct_method():
diff --git a/pytest/test_bem_potential_velocity_and_free_surface_elevation.py b/pytest/test_bem_potential_velocity_and_free_surface_elevation.py
index 20e3968d..61d936d7 100644
--- a/pytest/test_bem_potential_velocity_and_free_surface_elevation.py
+++ b/pytest/test_bem_potential_velocity_and_free_surface_elevation.py
@@ -221,7 +221,7 @@ def test_compute_free_surface_elevation_on_free_surface(solver, result):
 
 def test_pressure_integration(solver, result):
     f = result.body.integrate_pressure(solver.compute_pressure(result.body.mesh, result))
-    assert f == result.forces
+    assert f == pytest.approx(result.forces)
 
 def test_reconstruction_of_given_boundary_condition(solver, result):
     velocities = solver.compute_velocity(result.body.mesh, result)
diff --git a/pytest/test_bem_with_quadratures.py b/pytest/test_bem_with_quadratures.py
index 6daf5142..5277038d 100644
--- a/pytest/test_bem_with_quadratures.py
+++ b/pytest/test_bem_with_quadratures.py
@@ -32,6 +32,16 @@ def test_area(method):
         assert np.isclose(np.sum(mesh.quadrature_points[1][i_face, :]), mesh.faces_areas[i_face], rtol=1e-2)
 
 
+def test_mesh_in_body():
+    mesh = cpt.mesh_sphere()
+    mesh.compute_quadrature("Gauss-Legendre 2")
+    body = cpt.FloatingBody(mesh, cpt.rigid_body_dofs())
+    assert body.mesh.quadrature_method == "Gauss-Legendre 2"
+    points, weights = body.mesh.quadrature_points
+    assert points.shape == (mesh.nb_faces, 4, 3)
+    assert weights.shape == (mesh.nb_faces, 4)
+
+
 @pytest.mark.parametrize('method', list(builtin_methods) + quadpy_methods)
 def test_resolution(method):
     mesh = cpt.mesh_horizontal_cylinder(