Meshes¶
Naming¶
Meshes in Capytaine have a name
attribute.
It is optional and is mostly used for clearer logging and outputs.
A name
optional argument can be provided to all methods below to initialize
a mesh or transform a mesh to set the name of the new mesh.
Initialization¶
Importing with included Meshmagick readers¶
To load an existing mesh file, use the following syntax:
import capytaine as cpt
mesh = cpt.load_mesh('path/to/mesh.dat', file_format='nemoh')
body = cpt.FloatingBody(mesh=mesh)
The above example uses Nemoh’s mesh format, as defined e.g. on page 19 of Nemoh v3.0.1 manual.
Thanks to inclusion of code from Meshmagick,
numerous other mesh format can be imported.
The file format can be given with the file_format
optional argument.
If no format is given, the code will try to infer it from the file extension:
mesh = cpt.load_mesh('path/to/mesh.msh') # gmsh file
The formats currently supported in reading are listed in the following table (adapted from the documentation of Meshmagick).
File extension |
Software |
Keywords |
Extra features |
---|---|---|---|
.mar |
NEMOH [1] |
nemoh, mar |
Symmetries |
.nem |
NEMOH [1] |
nemoh_mesh, nem |
|
.gdf |
WAMIT [2] |
wamit, gdf |
Symmetries |
.inp |
DIODORE [3] |
diodore-inp, inp |
|
.DAT |
DIODORE [3] |
diodore-dat |
|
.pnl |
HAMS |
pnl, hams |
Symmetries |
.hst |
HYDROSTAR [4] |
hydrostar, hst |
Symmetries |
.nat |
natural, nat |
||
.msh |
GMSH [5] |
gmsh, msh |
|
.rad |
RADIOSS |
rad, radioss |
|
.stl |
stl |
||
.vtu |
PARAVIEW [6] |
vtu |
|
.vtp |
PARAVIEW [6] |
vtp |
|
.vtk |
PARAVIEW [6] |
paraview-legacy, vtk |
|
.tec |
TECPLOT [7] |
tecplot, tec |
|
.med |
SALOME [8] |
med, salome |
Not all metadata is taken into account when reading the mesh file.
For instance, the body symmetry is taken into account only for the .mar
, .pnl
, .gdf
and .hst
file formats.
Feel free to open an issue on Github to suggest improvements.
Importing with Meshio¶
Mesh can also be imported using the meshio library. Unlike the Meshmagick mesh readers mentioned above, this library is not packaged with Capytaine and need to be installed independently:
pip install meshio
The meshio mesh object can converted to Capytaine’s mesh
format with the load_from_meshio()
function:
from capytaine.io.meshio import load_from_meshio
cpt_mesh = cpt.load_from_meshio(mesh, name="My mesh")
This features allows to use pygmsh to generate the mesh, since this library returns mesh in the same format as meshio. Below is an example of a mesh generation with pygmsh (which also needs to be installed independently):
import pygmsh
offset = 1e-2
T1 = 0.16
T2 = 0.37
r1 = 0.88
r2 = 0.35
with pygmsh.occ.Geometry() as geom:
cyl = geom.add_cylinder([0, 0, 0], [0, 0, -T1], r1)
cone = geom.add_cone([0, 0, -T1], [0, 0, -T2], r1, r2)
geom.translate(cyl, [0, 0, offset])
geom.translate(cone, [0, 0, offset])
geom.boolean_union([cyl, cone])
gmsh_mesh = geom.generate_mesh(dim=2)
mesh = cpt.load_mesh(gmsh_mesh, name="my_pygmsh_mesh")
Predefined simple shapes¶
Capytaine include mesh generators for a few simple shapes. They are mostly
meant for teaching (they are extensively used in the examples of this
documentation) as well as for testing.
The most useful ones are
mesh_sphere()
,
mesh_vertical_cylinder()
,
mesh_horizontal_cylinder()
,
mesh_parallelepiped()
.
Some applications may also make use of flat shapes
mesh_disk()
and
mesh_rectangle()
.
Refer to their documentation for details about the parameters they accepts.
Since version 2.1, their resolution can be set by the faces_max_radius
parameter which specifies the maximal size of a face in the mesh.
Note
There are several ways to measure the size of a face and the resolution of a mesh. In Capytaine, the size of faces is usually quatified with the radius of the face, that is the maximal distance between the center of the face and its vertices. The resolution of a mesh is estimated as the maximal radius among all the faces in the mesh, that is the radius of the biggest face.
Creating from scratch¶
Alternatively, a mesh can be defined by giving a list of vertices and faces:
mesh = cpt.Mesh(vertices=..., faces=..., name="my_mesh")
The vertices are expected to be provided as a Numpy array of floats with shape (nb_vertices, 3)
.
The faces are provided as a Numpy array of ints with shape (nb_faces, 4)
, such that the four integers on a line are the indices of the vertices composing that face:
v = np.array([[0.0, 0.0, -1.0],
[1.0, 0.0, -1.0],
[1.0, 1.0, -1.0],
[0.0, 1.0, -1.0]])
f = np.array([[0, 1, 2, 3]])
single_face_mesh = cpt.Mesh(vertices=v, faces=f)
The ordering of the vertices define the direction of the normal vector, using normal right rotation. In other words, the normal vector is towards you if you see the vertices as being in counterclockwise order. In the above example, the normal vector is going up.
Triangular faces are supported as quadrilateral faces with the same vertex repeated twice:
single_triangle_mesh = cpt.Mesh(vertices=v, faces=np.array([[0, 1, 2, 2]]))
Creating a symmetric mesh¶
Several mesh symmetries can be used by Capytaine to speed up the computation.
The most useful one is the vertical plane symmetry.
A mesh with such a symmetry is stored by Capytaine with the
ReflectionSymmetricMesh
class.
It is defined with an other mesh of the half and a plane (and optionally a name
like the usual meshes):
half_mesh = cpt.load_mesh(...)
mesh = cpt.ReflectionSymmetricMesh(half_mesh, cpt.xOz_Plane, name="my full mesh")
Two vertical plane symmetries can be nested to be used by Capytaine (assuming that the two planes are orthogonal):
quarter_mesh = cpt.load_mesh(...)
half_mesh = cpt.ReflectionSymmetricMesh(half_mesh, cpt.yOz_Plane)
mesh = cpt.ReflectionSymmetricMesh(half_mesh, cpt.xOz_Plane)
All the method defined afterwards in this documentation should be applicable
for ReflectionSymmetricMesh
as well as for standard Mesh
.
You can consider using the clipped
method discussed below to create a symmetric mesh:
half_mesh = original_mesh.clipped(plane=cpt.xOz_Plane)
mesh = cpt.ReflectionSymmetricMesh(half_mesh, cpt.xOz_Plane)
Display¶
Use the show
method to display the mesh in 3D using VTK (if installed)
with the show()
:
mesh.show()
or with Matplotlib (if installed) with
show_matplotlib()
:
mesh.show_matplotlib()
Geometric transformations¶
Several functions are available to transform existing meshes.
Below is a list of most of the available methods. All of them can be applied to both meshes or to floating bodies, in which case the degrees of freedom will also be transformed:
# TRANSLATIONS
mesh.translated_x(10.0)
mesh.translated_y(10.0)
mesh.translated_z(10.0)
mesh.translated([10.0, 5.0, 2.0])
# Translation such that point_a would become equal to point_b
mesh.translated_point_to_point(point_a=[5, 6, 7], point_b=[4, 3, 2])
# ROTATIONS
mesh.rotated_x(3.14/5) # Rotation of pi/5 around the Ox axis
mesh.rotated_y(3.14/5) # Rotation of pi/5 around the Oy axis
mesh.rotated_z(3.14/5) # Rotation of pi/5 around the Oz axis
# Rotation of pi/5 around an arbitrary axis.
from capytaine import Axis
my_axis = Axis(vector=[1, 1, 1], point=[3, 4, 5])
mesh.rotated(axis=my_axis, angle=3.14/5)
# Rotation around a point such that vec1 would become equal to vec2
mesh.rotated_around_center_to_align_vector(
center=(0, 0, 0),
vec1=(1, 4, 7),
vec2=(9, 2, 1)
)
# REFLECTIONS
from capytaine import Plane
mesh.mirrored(Plane(normal=[1, 2, 1], point=[0, 4, 5]))
All the above methods can also be applied to Plane
and Axis
objects.
Meshes can also be merged together with the +
operator:
larger_mesh = mesh_1 + mesh_2
Finally, meshes can be clipped with a Plane
.
The plane is defined by a point belonging to it and a normal vector:
xOy_Plane = Plane(point=(0, 0, 0), normal=(0, 0, 1))
clipped_mesh = mesh.clipped(xOy_Plane)
Beware that the orientation of the normal vector of the Plane
will
determine which part of the mesh will be returned:
higher_part = mesh.clipped(Plane(point=(0, 0, 0), normal=(0, 0, -1)))
lower_part = mesh.clipped(Plane(point=(0, 0, 0), normal=(0, 0, 1)))
# mesh = lower_part + higher_part
The method immersed_part()
will clip the body with respect to two
horizontal planes at \(z=0\) and \(z=-h\):
clipped_body = mesh.immersed_part(water_depth=10)
Note
Most transformation methods exist in two versions:
one, named as a infinitive verb (translate, rotate, clip, keep_immersed_part, …), is an in-place transformation;
the other, named as a past participle (translated, rotated, clipped, immersed_part, …), is the same transformation but returning a new object.
In most cases, performance is not significant and the method returning a new object should be preferred. In-place transformation are currently kept for backward compatibility, but they make the code significantly more complicated and their removal might be considered in the future.
Defining an integration quadrature¶
During the resolution of the BEM problem, the Green function has to be integrated on each panel of the mesh. Parts of the Green function (such as the \(1/r\) Rankine terms) are integrated using an exact analytical expression for the integral. Other parts of the Green function rely on numerical integration. By default, this numerical integration is done by taking the value at the center of the panel and multiplying by its area. For a more accurate intagration, an higher order quadrature can be defined.
To define a quadrature scheme for a mesh, run the following command:
mesh.compute_quadrature(method="Gauss-Legendre 2")
The quadrature data can then be accessed at:
mesh.quadrature_points
and will be used automatically when needed.
Warning
Transformations of the mesh (merging, clipping, …) may reset the quadrature. Compute it only on your final mesh.
Warning
Quadratures schemes have been designed with quadrilateral panels. They work on triangular panels, but might not be as optimal then.
Alternatively, the compute_quadrature()
also accepts methods from the Quadpy package:
import quadpy
mesh.compute_quadrature(method=quadpy.c2.get_good_scheme(8))
Extracting or generating a lid¶
If you loaded a mesh file already containing a lid on the \(z=0\) plane,
the hull and the lid can be split with the
extract_lid()
method:
full_mesh = cpt.load_mesh(...)
hull_mesh, lid_mesh = full_mesh.extract_lid()
If your mesh does not have a lid, and you’d like to have irregular frequencies
removal, you can generate a lid using
generate_lid()
as follows:
lid_mesh = hull_mesh.generate_lid(z=-0.1)
Here the mesh is generated below the free surface, as Capytaine currently more robustly supports lid meshes below the free surface although they might not totally cancel all irregular frequencies.
An estimate of the lowest position at which the lid can be put for a given frequency can be computed as follows:
lid_mesh = hull_mesh.generate_lid(z=hull_mesh.lowest_lid_position(omega_max=10.0))
The resolution of the lid mesh can be set with the faces_max_radius
argument. By default, the mean resolution of the hull mesh is used.
See Floating body for detail on how to assign a lid mesh when defining a floating body.