# 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.
"""Generate meshes of cylinders and disks"""
import logging
from itertools import product
import numpy as np
from numpy import pi, cos, sin
from capytaine.meshes.meshes import Mesh
from capytaine.meshes.symmetric_meshes import RotationSymmetricMesh, ReflectionSymmetricMesh
LOG = logging.getLogger(__name__)
[docs]
def mesh_disk(*, radius=1.0, center=(0, 0, 0), normal=(0, 0, 1),
resolution=(3, 6), faces_max_radius=None,
reflection_symmetry=False, axial_symmetry=False, name=None, _theta_max=2*pi):
"""(One-sided) disk.
Parameters
----------
radius : float, optional
radius of the disk
center : 3-ple or array of shape (3,), optional
position of the geometric center of the disk
normal: 3-ple of floats, optional
normal vector, default: along x axis
resolution : 2-ple of int, optional
number of panels along a radius and around the disk
faces_max_radius : float, optional
maximal radius of a panel. (Default: no maximal radius.)
If the provided resolution is too coarse, the number of panels is
changed to fit the constraint on the maximal radius.
axial_symmetry : bool, optional
if True, returns an AxialSymmetricMesh
reflection_symmetry : bool, optional
if True, returns a ReflectionSymmetricMesh
name : str, optional
a string naming the mesh
_theta_max: float, optional
internal parameter, to return an arc circle instead of a full circle
"""
assert radius > 0, "Radius of the disk mesh should be given as a positive value."
assert len(resolution) == 2, "Resolution of a disk should be given as a couple of values."
assert all([h > 0 for h in resolution]), "Resolution of the disk mesh should be given as positive values."
assert all([i == int(i) for i in resolution]), "Resolution of a disk should be given as integer values."
assert len(center) == 3, "Position of the center of a disk should be given a 3-ple of values."
nr, ntheta = resolution
if faces_max_radius is not None:
# The biggest cell is on the side of the disk.
# Assuming it is a rectangle of sides radius/nr and perimeter/ntheta
estimated_max_radius = np.hypot(radius/nr, 2*pi*radius/ntheta)/2
if estimated_max_radius > faces_max_radius:
nr = int(np.ceil(radius / (np.sqrt(2) * faces_max_radius)))
ntheta = int(np.ceil(2*pi*radius / (np.sqrt(2) * faces_max_radius)))
if reflection_symmetry and axial_symmetry:
raise NotImplementedError("Disks with both symmetries have not been implemented.")
LOG.debug(f"New disk of radius {radius} and resolution {resolution}, named {name}.")
if reflection_symmetry:
if ntheta % 2 == 1:
raise ValueError("To use the reflection symmetry of the mesh, "
"it should have an even number of panels in this direction.")
half_mesh = mesh_disk(radius=radius, _theta_max=_theta_max/2, resolution=(nr, ntheta//2),
center=(0, 0, 0), normal=(0, 0, 1),
reflection_symmetry=False, axial_symmetry=False, name=f"half_of_{name}")
mesh = ReflectionSymmetricMesh(half_mesh, plane='xOz', name=name)
elif axial_symmetry:
mesh_slice = mesh_disk(radius=radius, _theta_max=_theta_max/ntheta, resolution=(nr, 1),
center=(0, 0, 0), normal=(0, 0, 1),
reflection_symmetry=False, axial_symmetry=False, name=f"slice_of_{name}")
mesh = RotationSymmetricMesh(mesh_slice, axis='z+', n=ntheta, name=name)
else:
theta_range = np.linspace(0, _theta_max, ntheta+1)
r_range = np.linspace(0.0, radius, nr+1)
nodes = np.array([(0, r*sin(t), -r*cos(t)) for (r, t) in product(r_range, theta_range)])
panels = np.array([(j+i*(ntheta+1), j+1+i*(ntheta+1), j+1+(i+1)*(ntheta+1), j+(i+1)*(ntheta+1))
for (i, j) in product(range(0, nr), range(0, ntheta))])
mesh = Mesh(nodes, panels, name=name, auto_check=False)
mesh = mesh.rotated_such_that_vectors_are_aligned(mesh.faces_normals[0], normal)
mesh = mesh.translated(center, name=name)
return mesh
[docs]
def mesh_vertical_cylinder(*, length=10.0, radius=1.0, center=(0, 0, 0),
resolution=(2, 8, 10), faces_max_radius=None,
axial_symmetry=False, reflection_symmetry=False, name=None, _theta_max=2*pi):
"""Vertical cylinder.
Total number of panels = (2*resolution[0] + resolution[2])*resolution[1]
Parameters
----------
length : float, optional
length of the cylinder
radius : float, optional
radius of the cylinder
center : 3-ple or array of shape (3,), optional
position of the geometric center of the cylinder
resolution : 3-ple of int, optional
(number of panel along a radius at the end, number of panels around a slice, number of slices)
Mnemonic: same ordering as the cylindrical coordinates (nr, ntheta, nz)
faces_max_radius : float, optional
maximal radius of a panel. (Default: no maximal radius.)
If the provided resolution is too coarse, the number of panels is
changed to fit the constraint on the maximal radius.
axial_symmetry : bool, optional
if True, returns an AxialSymmetricMesh
reflection_symmetry : bool, optional
if True, returns a ReflectionSymmetricMesh
name : str, optional
a string naming the mesh
_theta_max: float, optional
internal parameter, to return an arc circle instead of a full circle
"""
assert length > 0, "Length of a cylinder should be given as a positive value."
assert radius > 0, "Radius of a cylinder should be given as a positive value."
assert len(resolution) == 3, "Resolution of a cylinder should be given as a 3-ple of values."
assert all([h >= 0 for h in resolution]), "Resolution of a cylinder should be given as positive values."
assert all([i == int(i) for i in resolution]), "Resolution of a cylinder should be given as integer values."
assert len(center) == 3, "Position of the center of a cylinder should be given a 3-ple of values."
nr, ntheta, nz = resolution
if faces_max_radius is not None:
dr, dtheta, dz = radius/nr, 2*pi*radius/ntheta, length/nz
estimated_max_radius = max(
np.hypot(dr, dtheta)/2, # Panel on disk side
np.hypot(dtheta, dz)/2, # Panel on cylinder itself
)
if estimated_max_radius > faces_max_radius:
nr = int(np.ceil(radius / (np.sqrt(2) * faces_max_radius)))
ntheta = 2*int(np.ceil(pi*radius / (np.sqrt(2) * faces_max_radius)))
nz = int(np.ceil(length / (np.sqrt(2) * faces_max_radius)))
LOG.debug(f"New vertical cylinder of length {length}, radius {radius} and resolution {resolution}, named {name}.")
if reflection_symmetry and axial_symmetry:
raise NotImplementedError("Vertical cylinders with both symmetries have not been implemented.")
if reflection_symmetry:
if ntheta % 2 == 1:
raise ValueError("To use the reflection symmetry of the mesh, "
"it should have an even number of panels in this direction.")
half_cylinder = mesh_vertical_cylinder(length=length, radius=radius, center=(0, 0, 0),
resolution=(nr, ntheta//2, nz), reflection_symmetry=False, axial_symmetry=False,
name=f"half_{name}", _theta_max=_theta_max/2)
mesh = ReflectionSymmetricMesh(half_cylinder, plane='yOz', name=name)
elif axial_symmetry:
mesh_slice = mesh_vertical_cylinder(length=length, radius=radius, resolution=(nr, 1, nz), center=(0, 0, 0),
reflection_symmetry=False, axial_symmetry=False, name=f"slice_of_{name}", _theta_max=_theta_max/ntheta)
mesh = RotationSymmetricMesh(mesh_slice, axis='z+', n=ntheta, name=name)
else:
theta_range = np.linspace(0, _theta_max, ntheta+1)
z_range = np.linspace(-length/2, length/2, nz+1)
if nr > 0:
r_range = np.linspace(0.0, radius, nr+1)
nodes = np.concatenate([
np.array([(r*sin(t), r*cos(t), -length/2) for (r, t) in product(r_range, theta_range)]),
np.array([(radius*sin(t), radius*cos(t), z) for (z, t) in product(z_range, theta_range)]),
np.array([(r*sin(t), r*cos(t), length/2) for (r, t) in product(r_range[::-1], theta_range)]),
])
else:
r_range = np.array([])
nodes = np.array([(radius*sin(t), radius*cos(t), z) for (z, t) in product(z_range, theta_range)])
panels = np.array([(j+i*(ntheta+1), j+(i+1)*(ntheta+1), j+1+(i+1)*(ntheta+1), j+1+i*(ntheta+1), )
for (i, j) in product(range(nz+2*(len(r_range))), range(ntheta))])
mesh = Mesh(nodes, panels, name=name, auto_check=False)
mesh = mesh.translated(center, name=name)
return mesh
[docs]
def mesh_horizontal_cylinder(*, length=10.0, radius=1.0, center=(0, 0, 0),
resolution=(2, 8, 10), faces_max_radius=None,
reflection_symmetry=False, name=None, _theta_max=2*pi):
"""Cylinder aligned along Ox axis.
Total number of panels = (2*resolution[0] + resolution[2])*resolution[1]
Parameters
----------
length : float, optional
length of the cylinder
radius : float, optional
radius of the cylinder
center : 3-ple or array of shape (3,), optional
position of the geometric center of the cylinder
resolution : 3-ple of int, optional
(number of panel along a radius at the end, number of panels around a slice, number of slices)
Mnemonic: same ordering as the cylindrical coordinates (nr, ntheta, nz)
faces_max_radius : float, optional
maximal radius of a panel. (Default: no maximal radius.)
If the provided resolution is too coarse, the number of panels is
changed to fit the constraint on the maximal radius.
reflection_symmetry : bool, optional
if True, returns a ReflectionSymmetricMesh
translation_symmetry : bool, optional
if True, uses a TranslationalSymmetricMesh internally for the main part of the cylinder
name : str, optional
a string naming the mesh
_theta_max: float, optional
internal parameter, to return an arc circle instead of a full circle
"""
assert length > 0, "Length of a cylinder should be given as a positive value."
assert radius > 0, "Radius of a cylinder should be given as a positive value."
assert len(resolution) == 3, "Resolution of a cylinder should be given as a 3-ple of values."
assert all([h >= 0 for h in resolution]), "Resolution of a cylinder should be given as positive values."
assert all([i == int(i) for i in resolution]), "Resolution of a cylinder should be given as integer values."
assert len(center) == 3, "Position of the center of a cylinder should be given a 3-ple of values."
LOG.debug(f"New horizontal cylinder of length {length}, radius {radius} and resolution {resolution}, named {name}.")
nr, ntheta, nx = resolution
if faces_max_radius is not None:
dr, dtheta, dx = radius/nr, 2*pi*radius/ntheta, length/nx
estimated_max_radius = max(
np.hypot(dr, dtheta)/2, # Panel on disk side
np.hypot(dtheta, dx)/2, # Panel on cylinder itself
)
if estimated_max_radius > faces_max_radius:
nr = int(np.ceil(radius / (np.sqrt(2) * faces_max_radius)))
ntheta = 2*int(np.ceil(pi*radius / (np.sqrt(2) * faces_max_radius)))
nx = int(np.ceil(length / (np.sqrt(2) * faces_max_radius)))
if reflection_symmetry:
if ntheta % 2 == 1:
raise ValueError("To use the reflection symmetry of the mesh, "
"it should have an even number of panels in this direction.")
half_cylinder = mesh_horizontal_cylinder(
length=length, radius=radius, center=(0, 0, 0),
resolution=(nr, ntheta//2, nx),
reflection_symmetry=False,
name=f"half_{name}", _theta_max=_theta_max/2,
)
mesh = ReflectionSymmetricMesh(half_cylinder, plane='xOz', name=name)
else: # General case
theta_range = np.linspace(0, _theta_max, ntheta+1)
x_range = np.linspace(-length/2, length/2, nx+1)
nodes = np.array([(x, radius*sin(t), -radius*cos(t)) for (x, t) in product(x_range, theta_range)])
panels = np.array([(i+j*(ntheta+1), i+1+j*(ntheta+1), i+1+(j+1)*(ntheta+1), i+(j+1)*(ntheta+1))
for (i, j) in product(range(ntheta), range(nx))])
mesh = Mesh(nodes, panels, name=f"open_{name}", auto_check=False)
if nr > 0:
side = mesh_disk(radius=radius, center=(-length/2, 0, 0), normal=(-1, 0, 0),
reflection_symmetry=False, resolution=(nr, ntheta), name=f"side_of_{name}",
_theta_max=_theta_max)
other_side = side.mirrored('yOz', name=f"other_side_of_{name}")
mesh = Mesh.join_meshes(mesh, side, other_side, name=name)
mesh = mesh.translated(center, name=name)
return mesh