# Copyright 2026 Capytaine developers
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from abc import ABC, abstractmethod
from functools import lru_cache
import logging
import numpy as np
import xarray as xr
from capytaine.meshes.abstract_meshes import AbstractMesh
LOG = logging.getLogger(__name__)
[docs]
class AbstractDof(ABC):
def __str__(self):
return repr(self)
[docs]
@abstractmethod
def evaluate_motion_at_points(self, points: np.ndarray) -> np.ndarray:
# points is an array of shape (..., 3)
# output is of shape (..., 3)
...
[docs]
@lru_cache
def evaluate_motion(self, mesh: AbstractMesh) -> np.ndarray:
if mesh.nb_faces == 0:
return np.empty((mesh.nb_faces, 3))
else:
return self.evaluate_motion_at_points(mesh.faces_centers)
[docs]
@lru_cache
def evaluate_motion_at_quad_points(self, mesh: AbstractMesh) -> np.ndarray:
if mesh.nb_faces == 0:
return np.empty(mesh.quadrature_points[0].shape)
else:
return self.evaluate_motion_at_points(mesh.quadrature_points[0])
[docs]
@abstractmethod
def evaluate_gradient_of_motion_at_points(self, points: np.ndarray) -> np.ndarray:
# points is an array of shape (..., 3)
# output is of shape (..., 3, 3)
# output is a Jacobian matrix, such that output[..., i_dir, i_deriv_dir]
# is the derivative with respect to `i_deriv_dir` of the `i_dir` component of the motion.
# In other words, output[..., 0, :] is the gradient of the x-component of the motion on each face and
# output[..., :, 0] is the derivative with respect to x of the motion vector on each face.
...
[docs]
@lru_cache
def evaluate_gradient_of_motion(self, mesh: AbstractMesh) -> np.ndarray:
# output is of shape (nb_faces, 3, 3)
if mesh.nb_faces == 0:
return np.empty((mesh.nb_faces, 3, 3))
else:
return self.evaluate_gradient_of_motion_at_points(mesh.faces_centers)
[docs]
@lru_cache
def evaluate_gradient_of_motion_at_quad_points(self, mesh: AbstractMesh) -> np.ndarray:
# output is of shape (nb_faces, nb_quad_points, 3, 3)
if mesh.nb_faces == 0:
return np.empty((*mesh.quadrature_points[1].shape, 3, 3))
else:
return self.evaluate_gradient_of_motion_at_points(mesh.quadrature_points[0])
[docs]
class TranslationDof(AbstractDof):
def __init__(self, direction):
self.direction = np.asarray(direction)
assert self.direction.shape == (3,)
def __repr__(self):
return f"TranslationDof(direction={self.direction})"
[docs]
def evaluate_motion_at_points(self, points: np.ndarray) -> np.ndarray:
return np.tile(self.direction, (*points.shape[:-1], 1))
[docs]
def evaluate_gradient_of_motion_at_points(self, points: np.ndarray) -> np.ndarray:
return np.zeros((*points.shape[:-1], 3, 3))
[docs]
class RotationDof(AbstractDof):
def __init__(self, rotation_center, direction):
self.direction = np.asarray(direction)
assert self.direction.shape == (3,)
if rotation_center is None:
self.rotation_center = np.array([0.0, 0.0, 0.0])
LOG.warning("Rigid body rotation dof has been initialized "
"around the origin of the domain (0, 0, 0).")
else:
self.rotation_center = np.asarray(rotation_center, dtype=float)
assert self.rotation_center.shape == (3,)
def __repr__(self):
return f"RotationDof(rotation_center={self.rotation_center}, direction={self.direction})"
[docs]
def evaluate_motion_at_points(self, points: np.ndarray) -> np.ndarray:
return np.cross(self.direction, points - self.rotation_center)
[docs]
def evaluate_gradient_of_motion_at_points(self, points: np.ndarray) -> np.ndarray:
grad = np.cross(self.direction, np.eye(3)).T
# Transposing because np.cross compute the cross product row-wise,
# but we want it column-wise for the conventions of the Jacobian matrix.
return np.tile(grad, (*points.shape[:-1], 1, 1))
[docs]
class DofOnSubmesh(AbstractDof):
"""Defines a dof that is zeros everywhere except on a given set of faces
on which another dof object is used for the definition.
Parameters
----------
dof: AbstractDof
Some other dof
faces: boolean array (or slice or range)
The indices of the faces on which the `dof` should be evaluated.
Can be provided as a list of indices or a slice,
or as a boolean array of size `nb_faces` the total size of the mesh on which this dof is defined.
"""
def __init__(self, dof: AbstractDof, faces):
self.dof = dof
self.faces = faces
def __repr__(self):
return f"DofOnSubmesh(dof={self.dof}, faces={self.faces})"
[docs]
def evaluate_motion_at_points(self, points: np.ndarray) -> np.ndarray:
raise NotImplementedError(
"DofOnSubmesh.evaluate_motion_at_points() is not defined for arbitrary points: "
"membership in the submesh is only meaningful relative to the faces of the mesh "
"this dof was built from. Use evaluate_motion(mesh) instead."
)
[docs]
def evaluate_gradient_of_motion_at_points(self, points: np.ndarray) -> np.ndarray:
raise NotImplementedError(
"DofOnSubmesh.evaluate_gradient_of_motion_at_points() is not defined for arbitrary points: "
"membership in the submesh is only meaningful relative to the faces of the mesh "
"this dof was built from. Use evaluate_gradient_of_motion(mesh) instead."
)
[docs]
@lru_cache
def evaluate_motion(self, mesh: AbstractMesh) -> np.ndarray:
motion = np.zeros((mesh.nb_faces, 3))
motion[self.faces, :] = self.dof.evaluate_motion_at_points(mesh.faces_centers[self.faces, :])
return motion
[docs]
@lru_cache
def evaluate_gradient_of_motion(self, mesh: AbstractMesh) -> np.ndarray:
grad = np.zeros((mesh.nb_faces, 3, 3))
grad[self.faces, :, :] = self.dof.evaluate_gradient_of_motion_at_points(mesh.faces_centers[self.faces, :])
return grad
[docs]
@lru_cache
def evaluate_motion_at_quad_points(self, mesh: AbstractMesh) -> np.ndarray:
motion = np.zeros((*mesh.quadrature_points[1].shape, 3))
motion[self.faces, :, :] = self.dof.evaluate_motion_at_points(mesh.quadrature_points[0][self.faces, :, :])
return motion
[docs]
@lru_cache
def evaluate_gradient_of_motion_at_quad_points(self, mesh: AbstractMesh) -> np.ndarray:
grad = np.zeros((*mesh.quadrature_points[1].shape, 3, 3))
grad[self.faces, :, :, :] = self.dof.evaluate_gradient_of_motion_at_points(mesh.quadrature_points[0][self.faces, :, :])
return grad
[docs]
class CustomDof(AbstractDof):
"""Defines a fully custom dof with the same interface as the other AbstractDof.
To be used for elastic dofs.
Parameters
----------
motion: Callable
A function returning the motion at a given location. Required.
The function is expected to take a 3-element vector as input (location in space)
and return a 3-element vector (displacement at this location).
gradient_of_motion: Callable, optional
A function returning the gradient (actually jacobian matrix) of the motion
The function is expected to take a 3-element vector as input (location in space)
and return a 3×3-matrix (jacobian of the displacement at this location).
If none is provided, zero value is used, which makes no difference for first order hydrodynamics,
but loses some accuracy for hydrostatics and second-order forces.
"""
def __init__(self, motion, gradient_of_motion=None):
self.motion = motion
self.gradient_of_motion = gradient_of_motion
def __repr__(self):
grad = f", gradient_of_motion={self.gradient_of_motion}" if self.gradient_of_motion is not None else ""
return f"CustomDof(motion={self.motion}{grad})"
[docs]
def evaluate_motion_at_points(self, points: np.ndarray) -> np.ndarray:
flat_points = points.reshape(-1, 3)
motion = np.zeros((flat_points.shape[0], 3))
for i, p in enumerate(flat_points):
motion[i, :] = self.motion(p)
return motion.reshape(*points.shape[:-1], 3)
[docs]
def evaluate_gradient_of_motion_at_points(self, points: np.ndarray) -> np.ndarray:
flat_points = points.reshape(-1, 3)
grad = np.zeros((flat_points.shape[0], 3, 3))
if self.gradient_of_motion is not None:
for i, p in enumerate(flat_points):
grad[i, :, :] = self.gradient_of_motion(p)
return grad.reshape(*points.shape[:-1], 3, 3)
[docs]
def is_rigid_body_dof(dof):
# is_single_rigid_body_dof for now...
return (
isinstance(dof, TranslationDof)
or isinstance(dof, RotationDof)
# or (isinstance(dof, DofOnSubmesh) and is_rigid_body_dof(dof.dof))
)
[docs]
def rigid_body_dofs(only=None, rotation_center=None):
"""Pass this to FloatingBody initializer to give it rigid body dofs.
Parameters
----------
only: sequence of str, optional
list of the name of the rigid body dofs to be included.
By default: all six of them
rotation_center: np.array, optional
the center for the definition of the rotations
"""
if rotation_center is None:
rotation_center = np.array([0, 0, 0])
if not (only is not None
and set(only).issubset({"Surge", "Sway", "Heave"})
): # Skip the warning if only translations are required.
LOG.warning("Rigid body rotation dofs have been initialized "
"around the origin of the domain (0, 0, 0).")
# This warning is redundant with the one in RotationDof.__init__,
# but it is done here to have a single warning displayed on screen
# when a rigid body is initialized.
dofs = {
"Surge": TranslationDof(direction=(1, 0, 0)),
"Sway": TranslationDof(direction=(0, 1, 0)),
"Heave": TranslationDof(direction=(0, 0, 1)),
"Roll": RotationDof(rotation_center=rotation_center, direction=(1, 0, 0)),
"Pitch": RotationDof(rotation_center=rotation_center, direction=(0, 1, 0)),
"Yaw": RotationDof(rotation_center=rotation_center, direction=(0, 0, 1)),
}
if only is not None:
dofs = {k: v for k, v in dofs.items() if k in only}
return dofs
[docs]
def normalize_name(name):
return name[0].upper() + name[1:].lower()
[docs]
def add_dofs_labels_to_vector(dof_names, vector):
"""Helper function turning a bare vector into a vector labelled by the name of the dofs,
to be used for instance for the computation of RAO."""
return xr.DataArray(data=np.asarray(vector), dims=['influenced_dof'],
coords={'influenced_dof': list(dof_names)},
)
[docs]
def add_dofs_labels_to_matrix(dof_names, matrix):
"""Helper function turning a bare matrix into a matrix labelled by the name of the dofs,
to be used for instance for the computation of RAO."""
return xr.DataArray(data=np.asarray(matrix), dims=['influenced_dof', 'radiating_dof'],
coords={'influenced_dof': list(dof_names), 'radiating_dof': list(dof_names)},
)