diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml
index 5b57220f0..509960977 100644
--- a/.github/workflows/pr.yml
+++ b/.github/workflows/pr.yml
@@ -101,4 +101,4 @@ jobs:
       run: |
         git fetch --force --prune --tags
         git fetch --force --prune origin '+refs/heads/*:refs/remotes/origin/*'
-        python3 docs/build_docs.py --build production
+        python3 docs/build_docs.py --build debug
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index bea33b1f1..da9e713e7 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -33,7 +33,6 @@
     "import numpy as np\n",
     "import jax.numpy as jnp\n",
     "import capytaine as cpy\n",
-    "from capytaine.io.meshio import load_from_meshio\n",
     "import matplotlib.pyplot as plt\n",
     "from scipy.optimize import brute\n",
     "\n",
@@ -155,10 +154,10 @@
     "mesh = wb.mesh(mesh_size_factor)\n",
     "\n",
     "# create mesh object for WaveBot and add internal lid\n",
-    "mesh_obj = load_from_meshio(mesh, 'WaveBot')\n",
+    "mesh_obj = cpy.load_mesh(mesh)\n",
     "lid_mesh = mesh_obj.generate_lid(-2e-2)\n",
     "\n",
-    "fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"WaveBot\")"
+    "fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"WaveBot\", center_of_mass=[0,0,0])"
    ]
   },
   {
diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
index e12193bdb..a8bd7e16f 100644
--- a/examples/tutorial_2_AquaHarmonics.ipynb
+++ b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -23,7 +23,6 @@
    "outputs": [],
    "source": [
     "import capytaine as cpy\n",
-    "from capytaine.io.meshio import load_from_meshio\n",
     "import numpy as np\n",
     "import jax.numpy as jnp\n",
     "import matplotlib.pyplot as plt\n",
@@ -80,10 +79,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mesh_obj = load_from_meshio(mesh, 'WaveBot')\n",
+    "mesh_obj = cpy.load_mesh(mesh)\n",
     "lid_mesh = mesh_obj.generate_lid(-5e-2)\n",
     "\n",
-    "fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"AquaHarmonics\")\n",
+    "fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"AquaHarmonics\", center_of_mass=[0,0,0])\n",
     "fb.add_translation_dof(name=\"Heave\")\n",
     "ndof = fb.nb_dofs\n",
     "fb.show_matplotlib()"
@@ -107,8 +106,6 @@
    "source": [
     "g = 9.81\n",
     "rho = 1025\n",
-    "fb.center_of_mass = [0, 0, 0]\n",
-    "fb.rotation_center = fb.center_of_mass\n",
     "displaced_mass = fb.compute_rigid_body_inertia(rho=rho).values # kg\n",
     "displacement = displaced_mass/rho # m^3\n",
     "\n",
diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
index 4cd34c1bc..689152c4f 100644
--- a/examples/tutorial_3_LUPA.ipynb
+++ b/examples/tutorial_3_LUPA.ipynb
@@ -39,7 +39,6 @@
     "\n",
     "import gmsh, pygmsh\n",
     "import capytaine as cpy\n",
-    "from capytaine.io.meshio import load_from_meshio\n",
     "import numpy as np\n",
     "import jax.numpy as jnp\n",
     "import matplotlib.pyplot as plt\n",
@@ -155,11 +154,14 @@
    "outputs": [],
    "source": [
     "# lid\n",
-    "mesh_obj = load_from_meshio(mesh_float, 'float')\n",
+    "mesh_obj = cpy.load_mesh(mesh_float)\n",
     "lid_mesh = mesh_obj.generate_lid(z=-1e-2)\n",
     "\n",
+    "# center of mass float\n",
+    "cm_float = np.array(float_mass_properties['CG'])\n",
+    "\n",
     "# floating body\n",
-    "float_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"float\")\n",
+    "float_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"float\", center_of_mass=cm_float)\n",
     "float_fb.add_translation_dof(name='Heave')"
    ]
   },
@@ -228,11 +230,14 @@
    "outputs": [],
    "source": [
     "# lid\n",
-    "mesh_obj = load_from_meshio(mesh_spar, 'float')\n",
+    "mesh_obj = cpy.load_mesh(mesh_spar)\n",
     "lid_mesh = mesh_obj.generate_lid(z=-1e-2)\n",
     "\n",
+    "# center of mass spar\n",
+    "cm_spar = np.array(spar_mass_properties['CG'])\n",
+    "\n",
     "# floating body\n",
-    "spar_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"spar\")\n",
+    "spar_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"spar\", center_of_mass=cm_spar)\n",
     "spar_fb.add_translation_dof(name='Heave')"
    ]
   },
@@ -242,15 +247,30 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "print(spar_fb.mesh)\n",
+    "print(dir(spar_fb.mesh))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def axis_aligned_bbox(mesh):\n",
+    "    x, y, z = mesh.vertices.T\n",
+    "    return (x.min(), x.max(), y.min(), y.max(), z.min(), z.max())\n",
+    "\n",
     "# show\n",
     "fig = plt.figure()\n",
-    "xmin, xmax, ymin, ymax, zmin, zmax = spar_fb.mesh.axis_aligned_bbox\n",
+    "xmin, xmax, ymin, ymax, zmin, zmax = axis_aligned_bbox(spar_fb.mesh)\n",
     "scalez = (zmax-zmin) / (xmax-xmin)\n",
     "ax = fig.add_subplot(111, projection='3d')\n",
     "ax.set_box_aspect(aspect=(1,1,scalez))\n",
     "spar_fb.show_matplotlib(ax=ax)\n",
     "ax.set_xlim(xmin, xmax)\n",
     "ax.set_ylim(ymin, ymax)\n",
+    "ax.set_zlim(zmin, zmax)\n",
     "\n",
     "# interactive, close pop-up image before being able to continue running the notebook\n",
     "# spar_fb.show()"
@@ -282,39 +302,32 @@
     "\n",
     "# mass properties float\n",
     "mass_float = float_mass_properties['mass']\n",
-    "cm_float = np.array(float_mass_properties['CG'])\n",
     "pitch_inertia_float = float_mass_properties['MOI'][1]\n",
-    "float_fb.center_of_mass = cm_float\n",
-    "float_fb.rotation_center = float_fb.center_of_mass\n",
     "\n",
     "# mass properties spar\n",
     "mass_spar = spar_mass_properties['mass']\n",
-    "cm_spar = np.array(spar_mass_properties['CG'])\n",
     "pitch_inertia_spar = spar_mass_properties['MOI'][1]\n",
-    "spar_fb.center_of_mass = cm_spar\n",
-    "spar_fb.rotation_center = spar_fb.center_of_mass\n",
     "\n",
     "# floating body\n",
-    "lupa_fb = float_fb + spar_fb\n",
-    "lupa_fb.name = 'LUPA'\n",
+    "lupa_fb = (float_fb + spar_fb).as_FloatingBody\n",
     "\n",
-    " # mass properties LUPA\n",
-    "lupa_fb.center_of_mass = ((mass_float*cm_float + mass_spar*cm_spar)\n",
+    "# center of mass LUPA\n",
+    "cm_lupa = ((mass_float*cm_float + mass_spar*cm_spar)\n",
     "                        / (mass_float + mass_spar))\n",
-    "lupa_fb.rotation_center = lupa_fb.center_of_mass\n",
+    "\n",
+    "# additional DOFs\n",
+    "lupa_fb.add_translation_dof(name='Surge')\n",
+    "lupa_fb.add_rotation_dof(name='Pitch', rotation_center = cm_lupa)\n",
     "\n",
     "# pitch moment of inertia of LUPA using the parallel axis theorem \n",
-    "d_float = cm_float[2] - lupa_fb.center_of_mass[2]\n",
-    "d_spar = cm_spar[2] - lupa_fb.center_of_mass[2]\n",
+    "d_float = cm_float[2] - cm_lupa[2]\n",
+    "d_spar = cm_spar[2] - cm_lupa[2]\n",
     "pitch_inertia = (\n",
     "    pitch_inertia_float + mass_float*d_float**2 + \n",
     "    pitch_inertia_spar + mass_spar*d_spar**2\n",
     ")\n",
-    "inertia = np.diag([mass_float, mass_spar, lupa_fb.disp_mass(), pitch_inertia])\n",
     "\n",
-    "# additional DOFs\n",
-    "lupa_fb.add_translation_dof(name='Surge')\n",
-    "lupa_fb.add_rotation_dof(name='Pitch')"
+    "inertia = np.diag([mass_float, mass_spar, lupa_fb.disp_mass(), pitch_inertia])"
    ]
   },
   {
diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
index 96a065370..3c11b90b5 100644
--- a/examples/tutorial_4_Pioneer.ipynb
+++ b/examples/tutorial_4_Pioneer.ipynb
@@ -35,7 +35,6 @@
    "outputs": [],
    "source": [
     "import capytaine as cpy\n",
-    "from capytaine.io.meshio import load_from_meshio\n",
     "import numpy as np\n",
     "import jax.numpy as jnp\n",
     "import matplotlib.pyplot as plt\n",
@@ -257,13 +256,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mesh_obj = load_from_meshio(mesh, 'Pioneer')\n",
+    "mesh_obj = cpy.load_mesh(mesh)\n",
     "lid_mesh = mesh_obj.generate_lid(-4e-2)\n",
     "\n",
-    "pnr_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"WaveBot\")\n",
-    "pnr_fb.add_rotation_dof(name='Pitch')\n",
-    "pnr_fb.center_of_mass = np.array([0., 0., buoy_props['CG']])\n",
-    "pnr_fb.rotation_center = pnr_fb.center_of_mass\n",
+    "pnr_fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name=\"WaveBot\", center_of_mass=np.array([0., 0., buoy_props['CG']]))\n",
+    "pnr_fb.add_rotation_dof(name='Pitch', rotation_center = pnr_fb.center_of_mass)\n",
     "ndof = pnr_fb.nb_dofs\n",
     "pnr_fb.show_matplotlib()\n",
     "\n",
diff --git a/pyproject.toml b/pyproject.toml
index be12f5832..856b33d4d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,7 +23,7 @@ dependencies = [
     "scipy",
     "xarray",
     "jax",
-    "capytaine==2.3.1",
+    "capytaine>=2.3.1",
     "joblib",
     "wavespectra>=4.0",
     "netcdf4",
@@ -56,4 +56,4 @@ geometry = [
 [tool.coverage.run]
 source = [
     "wecopttool",
-]
\ No newline at end of file
+]
diff --git a/tests/test_core.py b/tests/test_core.py
index c37e0eed6..1df7744ca 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -1265,6 +1265,8 @@ def test_it_runs(self,):
         )
         rect = cpy.FloatingBody(rect_mesh, name="rect")
         rect.add_translation_dof(name="Heave")
+        if cpy.__version__ < '3':
+            rect.rotation_center = (0, 0, 0)  # Workaround for a bug in Capytaine 2.x. Actually unused since the body has only translation dofs.
         bem_data = wot.run_bem(fb=rect, freq=[0.1, 0.2], wave_dirs=[0,])
         assert type(bem_data) == xr.Dataset
 
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 3c3bbabe0..247c92ecf 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -4,7 +4,6 @@
 from pytest import approx
 import wecopttool as wot
 import capytaine as cpy
-from capytaine.io.meshio import load_from_meshio
 import numpy as np
 import jax.numpy as jnp
 from scipy.optimize import Bounds
@@ -85,10 +84,12 @@ def fb():
     mesh_size_factor = 0.2
     wb = geom.WaveBot()
     mesh = wb.mesh(mesh_size_factor)
-    mesh_obj = load_from_meshio(mesh, 'WaveBot')
+    mesh_obj = cpy.load_mesh(mesh)
     lid_mesh = mesh_obj.generate_lid(-2e-2)
-    fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name="WaveBot")
+    fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, center_of_mass=(0, 0, 0), name="WaveBot")
     fb.add_translation_dof(name="Heave")
+    if cpy.__version__ < '3':
+        fb.rotation_center = (0, 0, 0)  # Workaround for a bug in Capytaine 2.x. Actually unused since the body has only translation dofs.
     return fb
 
 
@@ -106,11 +107,13 @@ def fb_2d():
     mesh_size_factor = 0.2
     wb = geom.WaveBot()
     mesh = wb.mesh(mesh_size_factor)
-    mesh_obj = load_from_meshio(mesh, 'WaveBot')
+    mesh_obj = cpy.load_mesh(mesh)
     lid_mesh = mesh_obj.generate_lid(-2e-2)
-    fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name="WaveBot")
+    fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, center_of_mass=(0, 0, 0), name="WaveBot")
     fb.add_translation_dof(name="Heave")
     fb.add_translation_dof(name="Surge")
+    if cpy.__version__ < '3':
+        fb.rotation_center = (0, 0, 0)  # Workaround for a bug in Capytaine 2.x. Actually unused since the body has only translation dofs.
     return fb
 
 
@@ -198,8 +201,6 @@ def wec_from_impedance(hydro_data, pto, fb):
     w = np.expand_dims(omega, [1,2])
     A = bemc['added_mass'].values
     B = bemc['radiation_damping'].values
-    fb.center_of_mass = [0, 0, 0]
-    fb.rotation_center = fb.center_of_mass
     fb = fb.immersed_part()
     mass = bemc['inertia_matrix'].values
     hstiff = bemc['hydrostatic_stiffness'].values
@@ -223,8 +224,6 @@ def wec_from_impedance_2d(hydro_data_2d, pto_2d, fb_2d):
     w = np.expand_dims(omega, [1,2])
     A = bemc['added_mass'].values
     B = bemc['radiation_damping'].values
-    fb_2d.center_of_mass = [0, 0, 0]
-    fb_2d.rotation_center = fb_2d.center_of_mass
     fb_2d = fb_2d.immersed_part()
     mass = bemc['inertia_matrix'].values
     hstiff = bemc['hydrostatic_stiffness'].values
diff --git a/tests/test_utilities.py b/tests/test_utilities.py
index 4b7847e11..a65331eda 100644
--- a/tests/test_utilities.py
+++ b/tests/test_utilities.py
@@ -121,7 +121,7 @@ def fb():
     mesh_size_factor = 0.5
     wb = geom.WaveBot()
     mesh = wb.mesh(mesh_size_factor)
-    fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot")
+    fb = cpy.FloatingBody(mesh=cpy.load_mesh(mesh), name="WaveBot")
     fb.add_translation_dof(name="Heave")
     return fb
 
diff --git a/wecopttool/core.py b/wecopttool/core.py
index 81d27b732..56faf42ad 100644
--- a/wecopttool/core.py
+++ b/wecopttool/core.py
@@ -50,7 +50,6 @@
     "decompose_state",
     "frequency_parameters",
     "time_results",
-    "set_fb_centers",
     "block_diag_jax",
 ]
 
@@ -1882,7 +1881,7 @@ def write_netcdf(fpath: Union[str, Path], data: Dataset) -> None:
     --------
     read_netcdf,
     """
-    cpy.io.xarray.separate_complex_values(data).to_netcdf(fpath)
+    cpy.export_dataset(fpath, data, format="netcdf")
 
 
 def check_radiation_damping(
@@ -2200,7 +2199,9 @@ def run_bem(
                       'wavenumber': False,
                      }
     wec_im = fb.copy()
-    wec_im = set_fb_centers(wec_im, rho=rho)
+    if getattr(wec_im, "center_of_mass", None) is None:
+        wec_im.center_of_mass = wec_im.center_of_buoyancy
+        _log.warning("Using the geometric centroid as the center of gravity (COG).")
     if not hasattr(wec_im, 'inertia_matrix'):
         if wec_im.mass is None:
             wec_im.mass = rho*wec_im.immersed_part().volume
@@ -2600,41 +2601,6 @@ def time_results(fd: DataArray, time: DataArray) -> ndarray:
     return out
 
 
-def set_fb_centers(
-    fb: FloatingBody,
-    rho: float = _default_parameters["rho"],
-) -> FloatingBody:
-    """Sets default properties if not provided by the user:
-        - `center_of_mass` is set to the geometric centroid
-        - `rotation_center` is set to the center of mass
-    """
-    valid_properties = ['center_of_mass', 'rotation_center']
-
-    for property in valid_properties:
-        if not hasattr(fb, property):
-            setattr(fb, property, None)
-        if getattr(fb, property) is None:
-            if property == 'center_of_mass':
-                def_val = fb.center_of_buoyancy
-                log_str = (
-                    "Using the geometric centroid as the center of gravity (COG).")
-            elif property == 'rotation_center':
-                def_val = fb.center_of_mass
-                log_str = (
-                    "Using the center of gravity (COG) as the rotation center " +
-                    "for hydrostatics. Note that the hydrostatics do not use the " +
-                    "axes defined by the FloatingBody degrees of freedom, and the " +
-                    "rotation center should be set manually when using Capytaine to " +
-                    "calculate hydrostatics about an axis other than the COG.")
-            setattr(fb, property, def_val)
-            _log.warning(log_str)
-        elif getattr(fb, property) is not None:
-            _log.warning(
-                f'{property} already defined as {getattr(fb, property)}.')
-
-    return fb
-
-
 def block_diag_jax(*arrays: ArrayLike) -> ndarray:
     """Creates a block diagonal matrix from provided arrays.
 
@@ -2681,4 +2647,4 @@ def block_diag_jax(*arrays: ArrayLike) -> ndarray:
         r += rr
         c += cc
 
-    return block_diag_mat
\ No newline at end of file
+    return block_diag_mat