# Copyright 2025 Mews Labs
#
# 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 __future__ import annotations
import logging
from abc import ABC, abstractmethod
from functools import cached_property, lru_cache
from typing import Literal, Tuple
import numpy as np
from capytaine.meshes.surface_integrals import SurfaceIntegralsMixin
from capytaine.tools.deprecation_handling import _get_water_depth
from capytaine.meshes.geometry import connected_components, connected_components_of_waterline
LOG = logging.getLogger(__name__)
[docs]
class AbstractMesh(SurfaceIntegralsMixin, ABC):
@property
@abstractmethod
def nb_vertices(self) -> int:
...
@property
@abstractmethod
def nb_faces(self) -> int:
...
@property
@abstractmethod
def faces_normals(self) -> np.ndarray:
...
@property
@abstractmethod
def faces_areas(self) -> np.ndarray:
...
@property
@abstractmethod
def faces_centers(self) -> np.ndarray:
...
@property
@abstractmethod
def faces_radiuses(self) -> np.ndarray:
...
@property
@abstractmethod
def faces(self) -> np.ndarray:
...
@property
@abstractmethod
def quadrature_points(self) -> np.ndarray:
...
@cached_property
def z_span(self) -> Tuple[float, float]:
return (self.vertices[:, 2].min(), self.vertices[:, 2].max())
@abstractmethod
def __str__(self) -> str:
...
@abstractmethod
def __short_str__(self) -> str:
...
[docs]
@abstractmethod
def with_quadrature(self, quadrature_method):
...
[docs]
@abstractmethod
def translated(self, shift, *, name=None) -> AbstractMesh:
...
[docs]
def translated_x(self, dx: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh translated in the x-direction along `dx`."""
return self.translated([dx, 0.0, 0.0], name=name)
[docs]
def translated_y(self, dy: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh translated in the y-direction along `dy`."""
return self.translated([0.0, dy, 0.0], name=name)
[docs]
def translated_z(self, dz: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh translated in the z-direction along `dz`."""
return self.translated([0.0, 0.0, dz], name=name)
[docs]
@abstractmethod
def rotated_with_matrix(self, R, *, name=None) -> AbstractMesh:
...
[docs]
def rotated_x(self, angle: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh rotated around the x-axis using the provided rotation angle in radians"""
c, s = np.cos(angle), np.sin(angle)
R = np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
return self.rotated_with_matrix(R, name=name)
[docs]
def rotated_y(self, angle: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh rotated around the y-axis using the provided rotation angle in radians"""
c, s = np.cos(angle), np.sin(angle)
R = np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]])
return self.rotated_with_matrix(R, name=name)
[docs]
def rotated_z(self, angle: float, *, name=None) -> AbstractMesh:
"""Return a new Mesh rotated around the z-axis using the provided rotation angle in radians"""
c, s = np.cos(angle), np.sin(angle)
R = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])
return self.rotated_with_matrix(R, name=name)
[docs]
def rotated_such_that_vectors_are_aligned(self, a, b, *, eps=1e-8, name=None) -> AbstractMesh:
a = np.asarray(a, dtype=float)
b = np.asarray(b, dtype=float)
# Normalize input vectors
a_norm = np.linalg.norm(a)
b_norm = np.linalg.norm(b)
if a_norm < eps or b_norm < eps:
raise ValueError("Input vectors must be non-zero")
a_hat = a / a_norm
b_hat = b / b_norm
# Cross and dot products
v = np.cross(a_hat, b_hat)
c = np.dot(a_hat, b_hat)
s = np.linalg.norm(v)
# Case 1: vectors are already aligned
if s < eps and c > 0:
return self.copy(name=name)
# Case 2: vectors are opposite
if s < eps and c < 0:
# Find an arbitrary orthogonal vector
# Prefer axis least aligned with a_hat
axis = np.array([1.0, 0.0, 0.0])
if abs(a_hat[0]) > abs(a_hat[1]):
axis = np.array([0.0, 1.0, 0.0])
axis = axis - a_hat * np.dot(a_hat, axis)
axis /= np.linalg.norm(axis)
# Rotation by pi around axis
K = np.array([[0, -axis[2], axis[1]],
[axis[2], 0, -axis[0]],
[-axis[1], axis[0], 0]])
return self.rotated_with_matrix(np.eye(3) + 2 * K @ K, name=name)
# General case: Rodrigues' rotation formula
K = np.array([[0, -v[2], v[1]],
[v[2], 0, -v[0]],
[-v[1], v[0], 0]])
R = np.eye(3) + K + K @ K * ((1 - c) / (s ** 2))
return self.rotated_with_matrix(R, name=name)
[docs]
def edges_faces_water_line(self):
"""Extract the water line of the mesh.
Returns
-------
(np.ndarray, np.ndarray)
A tuple (edges, faces) where the edges are the edges of the water line (pairs of vertex indices)
and the faces are the corresponding faces of the mesh that contain these edges (quadruplet of vertex indices).
"""
epsilon = 1e-6
edges_water_line = []
faces_water_line = []
indices_vertices = [v for v in range(self.nb_vertices) if np.abs(self.vertices[v,-1]) <= epsilon]
for k in range(self.nb_faces):
set_indices = {index for index in self.faces[k,:] if index in indices_vertices}
if len(set_indices) == 2:
edges_water_line.append(list(set_indices))
faces_water_line.append(self.faces[k,:])
return np.array(edges_water_line), np.array(faces_water_line)
@property
def nb_edges_waterline(self) -> int:
"""Number of edges in the water line."""
return np.shape(self.edges_waterline)[0]
@cached_property
def edges_waterline(self):
"""Return an array with the edges of the water line as pairs of vertex indices with shape (nb_edges_waterline,2)."""
return self.edges_faces_water_line()[0]
@cached_property
def faces_water_line(self):
"""Return an array with the faces of the water line as quadruplets of vertex indices with shape (nb_edges_waterline,4)."""
return self.edges_faces_water_line()[1]
@cached_property
def length_edges_water_line(self):
"""Return an array with the lengths of the edges of the water line with shape (nb_edges_waterline,)."""
edges = self.edges_waterline
vertices_left = self.vertices[edges[:,0],:]
vertices_right = self.vertices[edges[:,1],:]
length_water_line = np.linalg.norm(vertices_left - vertices_right, ord=2, axis=1)
return length_water_line
[docs]
def water_line_integral(self, data):
"""Returns integral of given data along the water line.
Parameters
----------
data : np.ndarray
Values of the function to integrate, with shape (nb_edges_waterline,).
Returns
-------
float
Value of the integral.
"""
return np.sum(self.length_edges_water_line*data)
[docs]
def mirrored(self, plane: Literal['xOz', 'yOz'], *, name=None) -> AbstractMesh:
...
[docs]
@abstractmethod
def join_meshes(*meshes, return_masks=False, name=None) -> AbstractMesh:
...
def _common_metadata_keys(*meshes):
metadata_keys = [set(m.faces_metadata.keys()) for m in meshes]
common_metadata_keys = set.intersection(*metadata_keys)
lost_metadata_keys = set.union(*metadata_keys) - common_metadata_keys
if len(lost_metadata_keys) > 0:
LOG.warning(f'The following metadata have been dropped when joining meshes: {lost_metadata_keys}')
return common_metadata_keys
def __add__(self, other: AbstractMesh) -> AbstractMesh:
"""Combine two meshes using the + operator.
Parameters
----------
other : Mesh
Another mesh to combine with this one.
Returns
-------
Mesh
New mesh containing vertices and faces from both meshes.
"""
if self.name is not None or other.name is not None:
name = f"{self.name}+{other.name}"
else:
name = None
return self.join_meshes(other, name=name)
[docs]
def lowest_lid_position(self, omega_max, *, g=9.81):
z_lid = 0.0
for comp in connected_components(self):
for ccomp in connected_components_of_waterline(comp):
x_span = ccomp.vertices[:, 0].max() - ccomp.vertices[:, 0].min()
y_span = ccomp.vertices[:, 1].max() - ccomp.vertices[:, 1].min()
p = np.hypot(1/x_span, 1/y_span)
z_lid_comp = -np.arctanh(np.pi*g*p/omega_max**2) / (np.pi * p)
z_lid = min(z_lid, z_lid_comp)
return 0.9*z_lid # Add a small safety margin
[docs]
@abstractmethod
def generate_lid(self, z=0.0, faces_max_radius=None, name=None):
...
[docs]
@abstractmethod
def with_normal_vector_going_down(self, **kwargs) -> AbstractMesh:
...
[docs]
@abstractmethod
def copy(self) -> AbstractMesh:
...
[docs]
@abstractmethod
def merged(self) -> AbstractMesh:
...
[docs]
@abstractmethod
def clipped(self, *, origin, normal, name=None) -> AbstractMesh:
...
[docs]
@lru_cache
def immersed_part(self, free_surface=0.0, *, sea_bottom=None, water_depth=None) -> AbstractMesh:
"""
Clip the mesh to keep only the part below the free surface.
Parameters
----------
free_surface: float
The :math:`z` coordinate of the free surface (default: 0.0)
water_depth: Optional[float]
The water depth, as a positive value (default: infinity)
Returns
-------
Mesh
A new Mesh instance that has been clipped.
"""
water_depth = _get_water_depth(free_surface, water_depth, sea_bottom,
default_water_depth=np.inf)
if (free_surface - water_depth <= self.z_span[0]
and self.z_span[1] <= free_surface): # Already clipped
return self # Shortcut for performance
clipped = self.clipped(origin=(0, 0, 0), normal=(0, 0, 1))
if water_depth < np.inf:
clipped = clipped.clipped(origin=(0, 0, free_surface-water_depth), normal=(0, 0, -1))
return clipped
[docs]
@abstractmethod
def show(self, *, backend=None, **kwargs):
...
[docs]
def show_pyvista(self, **kwargs):
"""
Equivalent to show(backend="pyvista").
See also :func:`~capytaine.meshes.visualization.show_pyvista`
"""
return self.show(backend="pyvista", **kwargs)
[docs]
def show_matplotlib(self, **kwargs):
"""
Equivalent to show(backend="matplotlib").
See also :func:`~capytaine.meshes.visualization.show_matplotlib`
"""
return self.show(backend="matplotlib", **kwargs)
[docs]
@abstractmethod
def export(self, format, **kwargs):
...
[docs]
def export_to_pyvista(self, **kwargs):
return self.export(format="pyvista", **kwargs)
[docs]
def export_to_xarray(self, **kwargs):
return self.export(format="xarray", **kwargs)
[docs]
def export_to_meshio(self, **kwargs):
return self.export(format="meshio", **kwargs)
[docs]
def export_to_trimesh(self, **kwargs):
return self.export(format="trimesh", **kwargs)