============
Hydrostatics
============

Capytaine can compute hydrostatic parameters of a given :code:`FloatingBody` or :code:`Multibody`.

Some of them may require the definition of the center of mass of the body.
It can be done by setting the ``center_of_mass`` at body initialization as in the example below::

    import capytaine as cpt
    import numpy as np

    rigid_sphere = cpt.FloatingBody(
            mesh=cpt.mesh_sphere(radius=1.0, center=(0, 0, 0)),
            dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, -0.2)),
            center_of_mass=(0, 0, -0.4),
            name="sphere"  # optional
            )

The hydrostatic parameters can be computed and stored in an ``xarray.Dataset``
with the :func:`~capytaine.compute_hydrostatics_dataset` function::

   hs_dataset = cpt.compute_hydrostatics_dataset(rigid_sphere, rho=1025)

   # <xarray.Dataset> Size: 956B
   # Dimensions:                (space_coordinate: 3, influenced_dof: 6,
   #                             radiating_dof: 6)
   # Coordinates:
   #   * space_coordinate       (space_coordinate) <U1 12B 'x' 'y' 'z'
   #     center_of_mass         (space_coordinate) float64 24B 0.0 0.0 -0.4
   #     rotation_center        (space_coordinate) float64 24B 0.0 0.0 -0.2
   #   * influenced_dof         (influenced_dof) <U5 120B 'Surge' 'Sway' ... 'Yaw'
   #   * radiating_dof          (radiating_dof) <U5 120B 'Surge' 'Sway' ... 'Yaw'
   #     g                      float64 8B 9.81
   #     rho                    int64 8B 1025
   #     body                   <U6 24B 'sphere'
   # Data variables:
   #     hydrostatic_stiffness  (influenced_dof, radiating_dof) float64 288B 0.0 ....
   #     inertia_matrix         (influenced_dof, radiating_dof) float64 288B 1.959...
   #     center_of_buoyancy     (space_coordinate) float64 24B -2.178e-17 ... -0.3688
   #     draught                float64 8B 1.0
   #     disp_mass              float64 8B 1.959e+03

This dataset is also by default computed and included in the hydrodynamical
dataset created by :meth:`~capytaine.bem.solver.BEMSolver.fill_dataset`.

These hydrostatic parameters and others can each be accessed by a separate method::

    print("Volume:", rigid_sphere.volume)
    print("Center of buoyancy:", rigid_sphere.center_of_buoyancy)
    print("Wet surface area:", rigid_sphere.wet_surface_area)
    print("Displaced mass:", rigid_sphere.disp_mass(rho=1025))
    print("Waterplane center:", rigid_sphere.waterplane_center)
    print("Waterplane area:", rigid_sphere.waterplane_area)
    print("Metacentric parameters:",
        rigid_sphere.transversal_metacentric_radius,
        rigid_sphere.longitudinal_metacentric_radius,
        rigid_sphere.transversal_metacentric_height,
        rigid_sphere.longitudinal_metacentric_height)


Hydrostatic stiffness
---------------------

The method :meth:`~capytaine.bodies.FloatingBody.compute_hydrostatic_stiffness`
computes the hydrostatic stiffness and returns a (DOF count x DOF count) 2D
matrix as an :code:`xarray.DataArray`. ::

    hs = rigid_sphere.compute_hydrostatic_stiffness(rho=1025)

(It is also called in the backend of :func:`~capytaine.compute_hydrostatics_dataset`.)

Capytaine has two built-in methods to compute the hydrostatic stiffness:

- for rigid-body degrees of freedom, exact analytical expressions are available.

- for generalized degrees of freedom (e.g. elastic bodies), the following expression is used

.. math::

    C_{ij} = \iint_S (\hat{n} \cdot V_j) (w_i + z D_i)  dS

where :math:`\hat{n}` is the surface normal,
:math:`V_i = u_i \hat{n}_x + v_i \hat{n}_y + w_i \hat{n}_z` is the dof vector and
:math:`D_i = \nabla \cdot V_i` is the divergence of the DOF.

Note that the latter method is an approximation and does not recover the former
expression for rigid-body dofs, even when passing a dof vector and a divergence
corresponding to the rigid-body dof.

.. warning::
   Currently, the detection of a rigid-body degree of freedom to use the exact
   analytical expression is not very robust. Changing the name of the dof or
   combining several rigid bodies usually makes Capytaine use the approximate
   method for generalize dofs instead.
   This is planned to be fixed in a later release of Capytaine.


If the divergence is not provided by the user, zero is used as an approximation.
Here is an example of setting the divergence of the dof::

    mesh = cpt.mesh_sphere(radius=1.0, center=(0, 0, 0)).immersed_part()
    dof = np.array([(0, 0, z) for (x, y, z) in elastic_sphere.mesh.faces_centers])
    elastic_sphere = cpt.FloatingBody(
            mesh=mesh,
            center_of_mass=(0, 0, -0.3),
            dofs={"elongate_in_z": dof},
            )

    dofs_divergence = {"elongate_in_z": np.ones(elastic_sphere.mesh.nb_faces)}

    density = 1000.0
    gravity = 9.81
    elongate_in_z_hs = elastic_sphere.compute_hydrostatic_stiffness(divergence=dofs_divergence, rho=density, g=gravity)

    analytical_hs = - density * gravity * (4 * elastic_sphere.volume * elastic_sphere.center_of_buoyancy[2])
    print(np.isclose(elongate_in_z_hs.values[0, 0], analytical_hs))
    # True


.. warning::
   In the example above, all computations are done on a mesh clipped with
   ``immersed_part()``.
   Defining a divergence for generalized dofs is one of the only cases where
   the user should make sure that the body contains only immersed panels, and
   that the divergence field is defined only on these panels.
   It is planned to remove this requirement in a future release of Capytaine.

If the mass is not specified (as in the examples above), the body is assumed to
be in buoyancy equilibrium. Its mass is the mass of the displaced volume of
water.

Non-neutrally buoyant bodies are partially implemented (only for a single rigid body).
In this case, the mass of the body can be given by setting the :code:`mass`
attribute of the :code:`FloatingBody`::

    mesh = cpt.mesh_sphere(radius=1.0, center=(0, 0, 0))
    rigid_sphere = cpt.FloatingBody(
            mesh=mesh,
            dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, -0.3)),
            center_of_mass=(0, 0, -0.3),
            mass=0.8*mesh.disp_mass(rho=1000),
            )
    hs = rigid_sphere.compute_hydrostatic_stiffness()


Inertia matrix
--------------

The method :meth:`~capytaine.bodies.FloatingBody.compute_rigid_body_inertia` is
able to computes the 6 x 6 inertia matrix of a body with 6 rigid dofs.
The inertia coefficient of other degrees of freedom are filled with :code:`NaN` by default.

::

    M = rigid_sphere.compute_rigid_body_inertia()

As for the hydrostatic stiffness, the mass is assumed to be the displaced mass
of water for a body at equilibrium, unless a :code:`mass` attribute has been specified when setting up the body.

When computing the inertia moments, the density of the body mass is assumed to be uniform within the body (which may or may not be consistent with the center of mass defined when initializing the floating body).
Note also that unlike most other hydrostatics properties, the rigid body inertia depends not only on the shape of immersed part of the body but also on the shape above the free surface.
Computing the inertia on only the immersed part leads to a different result since it neglects this contribution to the moment of inertia::

    M_ = rigid_sphere.immersed_part().compute_rigid_body_inertia()

The latter might still be a relevant approximation since most floating bodies have most of their mass below the free surface for equilibrium.

This feature of Capytaine is limited in scope, it is meant to provide examples of data for teaching and testing.
The best practice is usually to use another dedicated method or tool that is aware of the true distribution of mass within the body.
When providing a matrix computed externally, it is recommended to wrap it in a :code:`xarray.DataArray` with dof names as labels, for consistency with the data computed with Capytaine::

    elastic_sphere.inertia_matrix = elastic_sphere.add_dofs_labels_to_matrix(np.array([[1000.0]]))