Cookbook

This page contains several examples of Capytaine’s features. The scripts can be downloaded individually as Python files from the Github folder with examples for version 2.2.

Base examples

Added mass of a rigid body

This example generates the mesh of an horizontal cylinder, solves radiation problems for the six rigid-body degrees of freedom and then plots the added mass.

#!/usr/bin/env python

import numpy as np
import capytaine as cpt

cpt.set_logging()

# Initialize floating body by generating a geometric mesh
cylinder = cpt.FloatingBody(
        mesh=cpt.mesh_horizontal_cylinder(
            length=10.0, radius=1.0,
            center=(0, 0, -2),
            resolution=(10, 20, 30)
            ))

# Automatically add the six degrees of freedom of a rigid body
cylinder.add_all_rigid_body_dofs()

# Define the range of frequencies as a Numpy array
omega_range = np.linspace(0.1, 3.0, 10)

# Set up the problems: we will solve a radiation problem for each
# degree of freedom of the body and for each frequency in the
# frequency range.
problems = [
    cpt.RadiationProblem(body=cylinder, radiating_dof=dof, omega=omega)
    for dof in cylinder.dofs
    for omega in omega_range
]
# Water density, gravity and water depth have not been specified.
# Default values are used.

# Solve all radiation problems
solver = cpt.BEMSolver()
results = solver.solve_all(problems)

# Gather the computed added mass into a labelled array.
data = cpt.assemble_dataset(results)

# Plot the added mass of each dofs as a function of the frequency
import matplotlib.pyplot as plt
plt.figure()
for dof in cylinder.dofs:
    plt.plot(
        omega_range,
        data['added_mass'].sel(radiating_dof=dof, influenced_dof=dof),
        label=dof,
        marker='o',
    )
plt.xlabel('omega')
plt.ylabel('added mass')
plt.legend()
plt.tight_layout()
plt.show()

Custom degree of freedom

This example defines arbitrary degrees of freedom for a sphere and solves a diffraction problem.

#!/usr/bin/env python

import numpy as np
import capytaine as cpt

# Initialize floating body
sphere = cpt.FloatingBody(
        mesh=cpt.mesh_sphere(
            radius=1.0,          # Dimension
            center=(0, 0, -2),   # Position
            resolution=(20, 20), # Fineness of the mesh
            ))

# DEFINE THE DOFS

# Manually defined heave,
# that is a vertical unit vector for all faces.
sphere.dofs["Heave"] = np.array(
    [(0, 0, 1) for face in sphere.mesh.faces]
)

# Bulging of the sphere,
# that is a deformation vector normal to the body at all faces.
# We can simply point to the normal vectors of the mesh for that.
sphere.dofs["Bulge"] = sphere.mesh.faces_normals

# Shearing of the sphere in the x direction.
# The deformation vector on a face is computed from the position of the face.
sphere.dofs["x-shear"] = np.array(
    [(np.cos(np.pi*z/2), 0, 0) for x, y, z in sphere.mesh.faces_centers]
)

# SOLVE DIFFRACTION PROBLEMS
solver = cpt.BEMSolver()

# Solve the problem for β=0 (incoming wave in the x direction).
problem_1 = cpt.DiffractionProblem(body=sphere, wave_direction=0, omega=1.0)
result_1 = solver.solve(problem_1)

# Solve the problem for β=π/2 (incoming wave in the y direction).
problem_2 = cpt.DiffractionProblem(body=sphere, wave_direction=np.pi/2, omega=1.0)
result_2 = solver.solve(problem_2)

# Print the generalized diffraction forces
# for the three dofs we defined
# for both values of the wave_direction β.
for result in [result_1, result_2]:
    print(f"Angle: {result.wave_direction:.2f}")
    for dof in sphere.dofs:
        force = result.forces[dof]
        print(f"{dof}: {np.abs(force):.2f}·exp({np.angle(force):.2f}i) N")
    print()

The diffraction force on the “Heave” and “Bulge” dofs should be the same for both incoming wave directions. The diffraction force on the “x-shear” dof is zero when the wave comes from the y direction.

Influence of the water depth

This example runs the same simulation for several water depth and plot the results.

#!/usr/bin/env python

import numpy as np
import capytaine as cpt

cpt.set_logging('INFO')

# Initialize floating body by generating a geometric mesh
mesh = cpt.mesh_horizontal_cylinder(
    length=10.0, radius=1.0,  # Dimensions
    center=(0, 0, -2),        # Position
    resolution=(5, 20, 40)    # Fineness of the mesh
    )
body = cpt.FloatingBody(mesh)

# Define a degree of freedom. The keyword "Heave"
# is recognized by the code and the vertical translation
# motion is automatically defined.
body.add_translation_dof(name="Heave")

# Define the range of water depth
depth_range = list(range(5, 25, 2)) + [np.inf]

# Set up the problems: we will solve a radiation problem for each
# water depth:
problems = [
    cpt.RadiationProblem(body=body, water_depth=depth, omega=2.0)
    for depth in depth_range
]
# Water density, gravity and radiating dof have not been specified.
# Default values are used. (For the radiating dof, the default value
# is usually the first one that has been defined. Here only one has
# been defined.)

# Solve all radiation problems
solver = cpt.BEMSolver(engine=cpt.HierarchicalToeplitzMatrixEngine())
results = [solver.solve(pb) for pb in sorted(problems)]

# Gather the computed added mass into a labelled array.
data = cpt.assemble_dataset(results)

# Plot the added mass of each dofs as a function of the water depth.
import matplotlib.pyplot as plt
plt.figure()
plt.plot(
    depth_range,
    data['added_mass'].sel(omega=2.0, radiating_dof="Heave", influenced_dof="Heave"),
    marker="s",
)
plt.xlabel('water depth')
plt.ylabel('added mass')
plt.show()

Irregular frequencies removal

import numpy as np
import capytaine as cpt

cpt.set_logging('INFO')

radius = 12.5
draught = 37.5
n_radius = 15
n_theta = 25
n_z = 30

# Define the range of frequencies to solve
omega_range = np.linspace(1.5, 2.25, 50)

solver = cpt.BEMSolver()

def compute_excitation_force(body, omega_range):
    problems = [
        cpt.DiffractionProblem(body=body, omega=omega)
        for omega in omega_range
    ]
    results = solver.solve_all(problems)
    data = cpt.assemble_dataset(results)
    return data['excitation_force']


# Initialize floating body by generating a geometric mesh
cylinder_mesh = cpt.mesh_vertical_cylinder(
    length=draught*2, radius=radius, center=(0, 0, 0), faces_max_radius=1.0,
    ).immersed_part()

### Body Without Lid (Body A)
body_1 = cpt.FloatingBody(mesh=cylinder_mesh)
body_1.add_translation_dof(name="Surge")

### Body With Lid at z = -3.699 (Body B)
lid_mesh = cylinder_mesh.generate_lid(z=-3.699)
body_2 = cpt.FloatingBody(mesh=cylinder_mesh, lid_mesh=lid_mesh)
body_2.add_translation_dof(name="Surge")

### Body With Lid at z = auto (Body C)
lid_mesh = cylinder_mesh.generate_lid(z=cylinder_mesh.lowest_lid_position(omega_max=omega_range.max()))
body_3 = cpt.FloatingBody(mesh=cylinder_mesh, lid_mesh=lid_mesh)
body_3.add_translation_dof(name="Surge")

force_1 = compute_excitation_force(body_1, omega_range)
force_2 = compute_excitation_force(body_2, omega_range)
force_3 = compute_excitation_force(body_3, omega_range)

# Plot the added mass of each dofs as a function of the frequency
import matplotlib.pyplot as plt
plt.figure()
plt.plot(
    omega_range,
    np.abs(force_1.sel(influenced_dof='Surge')),
    marker='+',
    label='no lid'
)
plt.plot(
    omega_range,
    np.abs(force_2.sel(influenced_dof='Surge')),
    marker='*',
    label=r'$z = -3.699$m'
)
plt.plot(
    omega_range,
    np.abs(force_3.sel(influenced_dof='Surge')),
    marker='x',
    label='auto'
)

plt.xlabel(r'$\omega$')
plt.ylabel(r'$F_1 (abs)$')
plt.legend()
plt.tight_layout()
plt.show()

Simulation with several bodies

#!/usr/bin/env python
# coding: utf-8

import capytaine as cpt

# Define the first body
sphere = cpt.FloatingBody(
        cpt.mesh_sphere(radius=1.0, center=(0, 0, -2.0),
                        resolution=(20, 20)),
        name="sphere_1")
sphere.add_translation_dof(name="Surge")
sphere.add_translation_dof(name="Heave")

# Define the second body
other_sphere = cpt.FloatingBody(
        cpt.mesh_sphere(radius=0.5, center=(-2, -3, -1),
                        resolution=(20, 20)),
        name="sphere_2")
other_sphere.add_translation_dof(name="Surge")
other_sphere.add_translation_dof(name="Heave")

# Define the third body
cylinder = cpt.FloatingBody(
        cpt.mesh_horizontal_cylinder(
            length=5.0, radius=1.0,
            center=(1.5, 3.0, -3.0),
            resolution=(20, 20, 3)),
        name="cylinder")
cylinder.add_translation_dof(name="Surge")
cylinder.add_translation_dof(name="Heave")

# Combine the three individual bodies into a single body.
all_bodies = cylinder + sphere + other_sphere

print("Merged body name:", all_bodies.name)
print("Merged body dofs:", list(all_bodies.dofs.keys()))

# The merged body can be used to define the problems in the usual way
problems = [cpt.RadiationProblem(body=all_bodies, radiating_dof=dof, omega=1.0) for dof in all_bodies.dofs]
problems += [cpt.DiffractionProblem(body=all_bodies, wave_direction=0.0, omega=1.0)]

# Solves the problem
solver = cpt.BEMSolver()
results = solver.solve_all(problems)
data = cpt.assemble_dataset(results)

print(data)

Plot pressure on hull

import numpy as np
import capytaine as cpt
import matplotlib.pyplot as plt

mesh = cpt.mesh_sphere(faces_max_radius=0.1).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)))

pb = cpt.RadiationProblem(wavelength=1.0, body=body, radiating_dof="Surge", water_depth=10.0)
solver = cpt.BEMSolver()
res = solver.solve(pb, keep_details=True)

body.show_matplotlib(
        color_field=np.real(res.pressure),
        cmap=plt.get_cmap("viridis"),  # Colormap
        )

Free surface elevation

#!/usr/bin/env python

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
import matplotlib.pyplot as plt

mesh = cpt.mesh_sphere(resolution=(10, 10)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs())

pb = cpt.DiffractionProblem(body=body, wave_direction=np.pi/2, wavelength=3.0)
solver = cpt.BEMSolver()
res = solver.solve(pb)

grid = np.meshgrid(np.linspace(-10.0, 10.0, 100), np.linspace(-10.0, 10.0, 100))
fse = solver.compute_free_surface_elevation(grid, res)
incoming_fse = airy_waves_free_surface_elevation(grid, res)

plt.pcolormesh(grid[0], grid[1], np.real(fse + incoming_fse))
plt.xlabel("x")
plt.ylabel("y")
plt.show()

Comparison with hydrostatics from Meshmagick

Hydrostatic and inertia properties can be computed independently via Capytaine or Meshmagick. This script compares them both with analytical expression for a simple geometric object.

import numpy as np

import capytaine as cpt
import meshmagick.mesh
import meshmagick.hydrostatics

radius = 10
cog = np.array((0, 0, 0))
mesh = cpt.mesh_sphere(radius=radius, center=cog, resolution=(100, 100))
body = cpt.FloatingBody(mesh=mesh, center_of_mass=cog)

body.add_all_rigid_body_dofs()
body = body.immersed_part()

density = 1000
gravity = 9.81

capy_hsdb = body.compute_hydrostatics(rho=density, g=gravity)

stiff_compare_dofs = ["Heave", "Roll", "Pitch"]
capy_hsdb["stiffness_matrix"] = capy_hsdb["hydrostatic_stiffness"].sel(
    influenced_dof=stiff_compare_dofs, radiating_dof=stiff_compare_dofs
    ).values

mass_compare_dofs = ["Roll", "Pitch", "Yaw"]
capy_hsdb["inertia_matrix"] = capy_hsdb["inertia_matrix"].sel(
    influenced_dof=mass_compare_dofs, radiating_dof=mass_compare_dofs
    ).values


mm_mesh = meshmagick.mesh.Mesh(body.mesh.vertices, body.mesh.faces, name=body.mesh.name)

mm_hsdb = meshmagick.hydrostatics.compute_hydrostatics(mm_mesh, cog=cog, rho_water=density, grav=gravity)

mm_hsdb["inertia_matrix"] = mm_mesh.eval_plain_mesh_inertias(rho_medium=density).inertia_matrix
mm_hsdb["center_of_buoyancy"] = mm_hsdb["buoyancy_center"]


analytical = {}
analytical["waterplane_area"] = np.pi*radius**2
analytical["wet_surface_area"] = 2*np.pi*radius**2
analytical["disp_volume"] = (2/3)*np.pi*radius**3
analytical["center_of_buoyancy"] = np.array([0,0,-3*radius/8])
inertia_xx = np.pi*radius**4/4
inertia_yy = np.pi*radius**4/4
inertia_zz = np.pi*radius**4/2
analytical["transversal_metacentric_radius"] = inertia_xx / analytical["disp_volume"]
analytical["longitudinal_metacentric_radius"] = inertia_yy / analytical["disp_volume"]
analytical["transversal_metacentric_height"] = analytical["transversal_metacentric_radius"] + analytical["center_of_buoyancy"][2] - cog[2]
analytical["longitudinal_metacentric_height"] = analytical["longitudinal_metacentric_radius"] + analytical["center_of_buoyancy"][2] - cog[2]

analytical["stiffness_matrix"] = density * gravity * np.array([
    [analytical["waterplane_area"], 0, 0],
    [0, analytical["disp_volume"] * analytical["transversal_metacentric_height"], 0],
    [0, 0, analytical["disp_volume"] * analytical["transversal_metacentric_height"]],
    ])

vars_to_be_displayed = capy_hsdb.keys() & (mm_hsdb.keys() | analytical.keys())
for var in vars_to_be_displayed:
    print(f"{var}:")
    print(f"    Capytaine  :  {capy_hsdb[var]}")
    if var in mm_hsdb:
        print(f"    Meshmagick :  {mm_hsdb[var]}")
    if var in analytical:
        print(f"    Analytical :  {analytical[var]}")

Intermediate examples

Animated free surface elevation

This example solves a diffraction problem, it computes the free surface elevation and shows it as a 3D animation.

#!/usr/bin/env python

import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
from capytaine.ui.vtk.animation import Animation

cpt.set_logging('INFO')

# Generate the mesh of a sphere

full_mesh = cpt.mesh_sphere(radius=3, center=(0, 0, 0), resolution=(20, 20))
full_sphere = cpt.FloatingBody(mesh=full_mesh)
full_sphere.add_translation_dof(name="Heave")

# Keep only the immersed part of the mesh
sphere = full_sphere.immersed_part()

# Set up and solve problem
solver = cpt.BEMSolver()

diffraction_problem = cpt.DiffractionProblem(body=sphere, wave_direction=0.0, omega=2.0)
diffraction_result = solver.solve(diffraction_problem)

radiation_problem = cpt.RadiationProblem(body=sphere, radiating_dof="Heave", omega=2.0)
radiation_result = solver.solve(radiation_problem)

# Define a mesh of the free surface and compute the free surface elevation
free_surface = cpt.FreeSurface(x_range=(-50, 50), y_range=(-50, 50), nx=150, ny=150)
diffraction_elevation_at_faces = solver.compute_free_surface_elevation(free_surface, diffraction_result)
radiation_elevation_at_faces = solver.compute_free_surface_elevation(free_surface, radiation_result)

# Add incoming waves
diffraction_elevation_at_faces = diffraction_elevation_at_faces + airy_waves_free_surface_elevation(free_surface, diffraction_problem)

# Run the animations
animation = Animation(loop_duration=diffraction_result.period)
animation.add_body(full_sphere, faces_motion=None)
animation.add_free_surface(free_surface, faces_elevation=0.5*diffraction_elevation_at_faces)
animation.run(camera_position=(-30, -30, 30))  # The camera is oriented towards (0, 0, 0) by default.
# animation.save("path/to/the/video/file.ogv", camera_position=(-30, -30, 30))

animation = Animation(loop_duration=radiation_result.period)
animation.add_body(full_sphere, faces_motion=full_sphere.dofs["Heave"])
animation.add_free_surface(free_surface, faces_elevation=3.0*radiation_elevation_at_faces)
animation.run(camera_position=(-30, -30, 30))
# animation.save("path/to/the/video/file.ogv", camera_position=(-30, -30, 30))

Animation of the RAO

This script generates the animation of the RAO motion for a wave incoming in front of a ship, such as the one used on the main page of this documentation. This script requires the mesh of the ship boat_200.mar. It can be downloaded from: https://raw.githubusercontent.com/capytaine/capytaine/master/docs/user_manual/examples/boat_200.mar

from numpy import pi

import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
from capytaine.ui.vtk import Animation

cpt.set_logging('INFO')

bem_solver = cpt.BEMSolver()


def generate_boat():
    boat_mesh = cpt.load_mesh("boat_200.mar", file_format="mar")
    boat = cpt.FloatingBody(
            mesh=boat_mesh,
            dofs=cpt.rigid_body_dofs(rotation_center=boat_mesh.center_of_buoyancy),
            center_of_mass = boat_mesh.center_of_buoyancy,
            name="pirate ship"
            )
    boat.inertia_matrix = boat.compute_rigid_body_inertia() / 10 # Artificially lower to have a more appealing animation
    boat.hydrostatic_stiffness = boat.immersed_part().compute_hydrostatic_stiffness()
    return boat


def setup_animation(body, fs, omega, wave_amplitude, wave_direction):
    # SOLVE BEM PROBLEMS
    radiation_problems = [cpt.RadiationProblem(omega=omega, body=body.immersed_part(), radiating_dof=dof) for dof in body.dofs]
    radiation_results = bem_solver.solve_all(radiation_problems)
    diffraction_problem = cpt.DiffractionProblem(omega=omega, body=body.immersed_part(), wave_direction=wave_direction)
    diffraction_result = bem_solver.solve(diffraction_problem)

    dataset = cpt.assemble_dataset(radiation_results + [diffraction_result])
    rao = cpt.post_pro.rao(dataset, wave_direction=wave_direction)

    # COMPUTE FREE SURFACE ELEVATION
    # Compute the diffracted wave pattern
    incoming_waves_elevation = airy_waves_free_surface_elevation(fs, diffraction_result)
    diffraction_elevation = bem_solver.compute_free_surface_elevation(fs, diffraction_result)

    # Compute the wave pattern radiated by the RAO
    radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs, res) for res in radiation_results}
    radiation_elevation = sum(rao.sel(omega=omega, radiating_dof=dof).data * radiation_elevations_per_dof[dof] for dof in body.dofs)

    # SET UP ANIMATION
    # Compute the motion of each face of the mesh for the animation
    rao_faces_motion = sum(rao.sel(omega=omega, radiating_dof=dof).data * body.dofs[dof] for dof in body.dofs)

    # Set up scene
    animation = Animation(loop_duration=2*pi/omega)
    animation.add_body(body, faces_motion=wave_amplitude*rao_faces_motion)
    animation.add_free_surface(fs, wave_amplitude * (incoming_waves_elevation + diffraction_elevation + radiation_elevation))
    return animation


if __name__ == '__main__':
    body = generate_boat()
    fs = cpt.FreeSurface(x_range=(-100, 75), y_range=(-100, 75), nx=100, ny=100)

    anim = setup_animation(body, fs, omega=1.5, wave_amplitude=0.5, wave_direction=pi)
    anim.run(camera_position=(70, 70, 100), resolution=(800, 600))
    anim.save("animated_boat.ogv", camera_position=(70, 70, 100), resolution=(800, 600))

Kochin function

This example computes the Kochin function for a surging buoy and plot the results.

import numpy as np
from numpy import pi
import xarray as xr
import matplotlib.pyplot as plt

import capytaine as cpt

cpt.set_logging('INFO')

# SET UP THE MESH
buoy = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            radius=1.0, length=3.0, center=(0, 0, -1.0),
            resolution=(15, 15, 3)
            ))
buoy.add_translation_dof(name="Surge")

# SOLVE THE BEM PROBLEMS AND COMPUTE THE KOCHIN FUNCTIONS
theta_range = np.linspace(0, 2*pi, 40)
omega_range = np.linspace(1.0, 5.0, 3)

test_matrix = xr.Dataset(coords={
    'omega': omega_range, 'theta': theta_range, 'radiating_dof': ["Surge"],
})
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, buoy.immersed_part(), wavenumber=True)

# PLOT THE KOCHIN FUNCTION
plt.figure()
dataset = dataset.swap_dims({'omega': 'wavenumber'})  # Use wavenumber to index data
for k in dataset.coords['wavenumber']:
    (dataset['kochin']
     .sel(radiating_dof="Surge", wavenumber=k)
     .real
     .plot(x='theta', label=f"k={float(k):.2f} m¯¹", ax=plt.gca()))
plt.legend()
plt.title("Real part of the Kochin function as a function of theta")
plt.tight_layout()
plt.show()

Haskind’s relation

This example computes the excitation force from the radiation potential using Haskind’s relation. The result is compared with the one obtained by direct integration of the potentials from incident waves and from the diffraction problem.

#! /usr/bin/env python3

import numpy as np
import capytaine as cpt


r = 1.
draft = 0.5
depth = 10.
omega = 1.
rho = 1000

# Define geometry and heave degree of freedom
body = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            length=2*draft, radius=r, center=(0.,0.,0.),
            resolution=(10, 50, 10)
            ))
body.add_translation_dof(name='Heave')
body = body.immersed_part()

solver = cpt.BEMSolver()

# Define and solve the diffraction and radiation problems
diff_problem = cpt.DiffractionProblem(body=body, water_depth=depth,
                                      omega=omega, wave_direction=0.)
rad_problem = cpt.RadiationProblem(body=body, water_depth=depth,
                                   omega=omega, radiating_dof='Heave')

diff_solution = solver.solve(diff_problem)
rad_solution = solver.solve(rad_problem)


# Read mesh properties
faces_centers = body.mesh.faces_centers
faces_normals = body.mesh.faces_normals
faces_areas = body.mesh.faces_areas

from capytaine.bem.airy_waves import airy_waves_potential, airy_waves_velocity, froude_krylov_force

# Computation from the diffraction solution (Capytaine)
FK = froude_krylov_force(diff_problem)['Heave']
diff_force = diff_solution.forces['Heave']

# Get potentials
phi_inc = airy_waves_potential(faces_centers, diff_problem)
v_inc = airy_waves_velocity(faces_centers, diff_problem)
phi_rad = rad_solution.potential

# Indirect computation from the radiation solution, via the Haskind relation
integrand = - (phi_inc * faces_normals[:,2]
              - 1j/omega * phi_rad * np.diag(v_inc@faces_normals.T))
F_hask = 1j * omega * rho * np.sum(integrand*faces_areas)

# Direct integration of the potentials
diff_force_recomputed = - 1j * omega * rho * np.sum(
                        diff_solution.potential*faces_areas*faces_normals[:,2])
FK_recomputed = - 1j * omega * rho * np.sum(
                        phi_inc*faces_areas*faces_normals[:,2])

print('Result from Capytaine = ', FK + diff_force)
print('Result from recomputed direct integration',
        FK_recomputed + diff_force_recomputed)
print('Result from Haskind relation = ', F_hask)

Axisymmetric body

This example generates an axisymmetric mesh from a profile function and solves radiation problems for this floating body.

#!/usr/bin/env python

import numpy as np
import capytaine as cpt
import capytaine.io.xarray

cpt.set_logging('INFO')

# Profile of the axisymmetric body
def shape(z):
    return 0.1*(-(z+1)**2 + 16)

# Generate the mesh and display it with VTK.
buoy = cpt.FloatingBody(
    mesh=cpt.AxialSymmetricMesh.from_profile(shape, z_range=np.linspace(-5, 0, 30), nphi=40)
)
buoy.add_translation_dof(name="Heave")
buoy.show()

# Set up problems
omega_range = np.linspace(0.1, 5.0, 60)
problems = [cpt.RadiationProblem(body=buoy, radiating_dof='Heave', omega=omega)
            for omega in omega_range]

# Solve the problems using the axial symmetry
solver = cpt.BEMSolver(engine=cpt.HierarchicalToeplitzMatrixEngine())
results = [solver.solve(pb) for pb in sorted(problems)]
dataset = capytaine.io.xarray.assemble_dataset(results)

# Plot results
import matplotlib.pyplot as plt
plt.figure()
plt.plot(
    omega_range,
    dataset['added_mass'].sel(radiating_dof='Heave',
                              influenced_dof='Heave'),
    label="Added mass",
)
plt.plot(
    omega_range,
    dataset['radiation_damping'].sel(radiating_dof='Heave',
                                     influenced_dof='Heave'),
    label="Radiation damping",
)
plt.xlabel('omega')
plt.grid()
plt.legend()
plt.show()

Advanced examples

Convergence study

This example runs a mesh convergence study for a submerged cylinder.

import numpy as np
from numpy import pi
import xarray as xr
import matplotlib.pyplot as plt
import capytaine as cpt

cpt.set_logging('INFO')

def make_cylinder(resolution):
    """Make cylinder with a mesh of a given resolution in panels/meter."""
    radius = 1.0
    length = 5.0
    mesh = cpt.mesh_horizontal_cylinder(
        length=length, radius=radius,
        center=(0, 0, -1.5*radius),
        resolution=(int(resolution*radius),
                    int(2*pi*length*resolution),
                    int(length*resolution))
    )
    body = cpt.FloatingBody(mesh, name=f"cylinder_{mesh.nb_faces:04d}")
    # Store the number of panels in the name of the body
    body.add_translation_dof(name="Heave")
    return body

test_matrix = xr.Dataset(coords={
    'omega': [1.0],
    'radiating_dof': ['Heave'],
})

bodies = [make_cylinder(n) for n in np.linspace(1, 5, 10)]
ds1 = cpt.BEMSolver().fill_dataset(test_matrix, bodies)

def read_nb_faces_in_mesh_name(ds):
    """Read the name of the body to guess the resolution of the mesh."""
    ds.coords['nb_faces'] = xr.DataArray([int(name[9:]) for name in ds['body_name'].values], coords=[ds['body_name']])
    ds = ds.swap_dims({'body_name': 'nb_faces'})
    return ds
ds1 = read_nb_faces_in_mesh_name(ds1)

ds1['added_mass'].plot(x='nb_faces')

plt.show()

Plot the influence matrix

This example plots the influence matrix for an horizontal cylinder.

#!/usr/bin/env python

import numpy as np
import capytaine as cpt

# Generate the mesh of a cylinder
cylinder = cpt.FloatingBody(
        mesh=cpt.mesh_horizontal_cylinder(
            length=10.0, radius=1.0,
            center=(0, 0, -2),
            resolution=(1, 6, 8)
            ))

engine = cpt.BasicMatrixEngine()
green_function = cpt.Delhommeau()

S, K = engine.build_matrices(
    cylinder.mesh, cylinder.mesh,
    free_surface=0.0, water_depth=np.inf,
    wavenumber=1.0,
    green_function=green_function,
)

# Plot the absolute value of the matrix S
import matplotlib.pyplot as plt
plt.imshow(abs(np.array(S)))
plt.colorbar()
plt.title("$|S|$")
plt.tight_layout()
plt.show()

Compare two implementations of the Green function

This is an example of comparison of two implementations of the Green function.

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import capytaine as cpt

cpt.set_logging('INFO')

# Generate body
mesh = cpt.mesh_horizontal_cylinder(
    length=3.0, radius=1.0, center=(0, 0, -1.01), resolution=(5, 30, 15)
)
body = cpt.FloatingBody(mesh)
body.add_translation_dof(name="Heave")

test_matrix = xr.Dataset(coords={
    'omega': np.linspace(0.5, 4, 40),
    'radiating_dof': list(body.dofs.keys()),
})

ds2 = cpt.BEMSolver(green_function=cpt.Delhommeau(gf_singularities="high_freq")).fill_dataset(test_matrix, body)
ds1 = cpt.BEMSolver(green_function=cpt.Delhommeau(gf_singularities="low_freq")).fill_dataset(test_matrix, body)

plt.figure()
ds1['added_mass'].plot(x='omega', label='High freq singularities')
ds2['added_mass'].plot(x='omega', label='Low freq singularities')
plt.legend()

plt.figure()
ds1['radiation_damping'].plot(x='omega', label='High freq singularities')
ds2['radiation_damping'].plot(x='omega', label='Low freq singularities')
plt.legend()

plt.show()

Use a custom Green function

This is an example of how to implement a custom Green function.

#!/usr/bin/env python

import numpy as np
from numpy import dot, pi
from numpy.linalg import norm

import capytaine as cpt
from capytaine.green_functions.abstract_green_function import AbstractGreenFunction

class MyGreenFunction(AbstractGreenFunction):
    """An example of a custom routine to evaluate the Green function."""

    def evaluate(self, mesh1, mesh2, free_surface, water_depth, wavenumber, **kwargs):
        """The main method that needs to be implemented in the class."""

        if free_surface == np.inf and water_depth == np.inf:

            # Initialize the matrices
            S = np.zeros((mesh1.nb_faces, mesh2.nb_faces))
            K = np.zeros((mesh1.nb_faces, mesh2.nb_faces))

            for i in range(mesh1.nb_faces):
                for j in range(mesh2.nb_faces):
                    area = mesh1.faces_areas[i]
                    if i != j:
                        # The integral on the face is computed using the value at the center of the face.
                        r = mesh1.faces_centers[i, :] - mesh2.faces_centers[j, :]

                        S[j, i] = -area/(4*pi*norm(r))
                        K[j, i] = -area*dot(r, mesh1.faces_normals[j, :])/(4*pi*norm(r)**2)

                    else:
                        # The kernel is singular.
                        # When i == j, we can't use the approximation above for the integral.
                        # Below, we take a rough approximation.
                        S[j, i] = -area/(4*pi)*3
                        K[j, i] = -area/(4*pi)*3

            if mesh1 is mesh2:
                for i in range(mesh1.nb_faces):
                    K[i, i] += 1/2

            return S, K

        else:
            raise NotImplementedError

# Define a solver using the Green function above.
solver = cpt.BEMSolver(green_function=MyGreenFunction())

# Example of a problem
sphere = cpt.FloatingBody(
        mesh=cpt.mesh_sphere(
            radius=1.0,          # Dimension
            center=(0, 0, -2),   # Position
            resolution=(4, 4),   # Fineness of the mesh
            ))
sphere.add_translation_dof(name="Heave")
problem = cpt.RadiationProblem(body=sphere, free_surface=np.inf, water_depth=np.inf, radiating_dof='Heave')

result = solver.solve(problem)
print(result.added_masses)