"""
This module provides various predefined geometric shapes in 2D and 3D,
including lines, arcs, rectangles, circles, blocks, cylinders, tori, and spheres.
"""
from warnings import warn
from .utils import *
from .ops import *
from typing import Sequence
[docs]class Line(Geometry):
"""
1D line aligned to a principal axis, then reoriented.
Shortest import: `from geoparticle import Line`
"""
[docs] def __init__(self, length: float, direction: str, dl: float,
anchor: Sequence[float] = (0, 0, 0),
name=None):
"""
Args:
length (float): Length of the line.
direction (str): Principal axis direction ('x', 'y', or 'z').
dl (float): Spacing between points along the line.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the line. The anchor is the line end with
a smaller coordinate along the principal axis. Defaults to (0,0,0).
name (str, optional): Name of the line. Defaults to None.
"""
super().__init__(name=name or f'Line {self.get_counter()}', dimension=2)
zs = _arange0_quantized(length, dl)
self.length = float(zs[-1])
_check_size_change(length, self.length, self.name, 'length')
xs = np.zeros_like(zs)
ys = np.zeros_like(zs)
self.set_coord(*_transform_coordinate(xs, ys, zs, axis=direction))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class SymmLines(Geometry):
"""
Two symmetrical lines centered at origin and aligned along `direction`.
Shortest import: `from geoparticle import SymmLines`
"""
[docs] def __init__(self, length: float, direction: str, symm_plane: str,
dist_half: float,
dl: float,
anchor: Sequence[float] = (0, 0, 0),
name=None):
"""
Args:
length (float): Length of each symmetrical line.
direction (str): Principal axis direction ('x', 'y', or 'z').
symm_plane (str): Symmetrical plane ('XOY', 'YOZ', or 'XOZ').
dist_half (float): Half the distance between the two lines.
dl (float): Spacing between points along the lines.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the symmetrical lines. The anchor is one of the
symmetrical line ends with a smaller coordinate along the principal axis. Defaults to (0,0,0).
name (str, optional): Name of the symmetrical lines. Defaults to None.
"""
super().__init__(name=name or f'SymmLines {self.get_counter()}', dimension=2)
upper = Line(length, direction, dl)
normal = _resolve_axis_or_plane(plane=symm_plane)
eval(f'upper.shift({normal}=dist_half, inplace=True)')
lower = upper.mirror(symm_plane, 0)
self.load_from(Union((upper, lower)))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class Arc(Geometry):
"""
2D circular arc on a plane.
Shortest import: `from geoparticle import Arc`
"""
[docs] def __init__(self, r: float, phi_range: str, dl: float, plane: str = 'XOY', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r (float): Radius of the arc.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full circle.
dl (float): Spacing between points along the arc.
plane (str): Plane in which the arc lies ('XOY', 'YOZ', or 'XOZ').
name (str, optional): Name of the arc. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the arc. The anchor is the arc center.
"""
super().__init__(name=name or f'Arc {self.get_counter()}', dimension=2)
a_deg, b_deg, incl_min, incl_max, total_deg = _parse_interval_deg(phi_range)
self.phi_tot_deg = total_deg
phis = _discretize_arc_by_dl(r, dl, a_deg, b_deg, incl_min, incl_max)
xs = r * np.cos(phis)
ys = r * np.sin(phis)
zs = np.zeros_like(xs)
self.set_coord(*_transform_coordinate(xs, ys, zs, plane=plane))
# apply anchor translation if provided
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class ConcentricArc(Geometry):
"""
Two concentric arcs. Can be used to create a 2D torus section.
Shortest import: `from geoparticle import ConcentricArc`
"""
[docs] def __init__(self, r_out: float, r_in: float, dl: float, plane: str = 'XOY',
phi_range='[0,360)', name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_out (float): Outer radius.
r_in (float): Inner radius.
dl (float): Spacing between points along the arcs.
plane (str, optional): Plane in which the torus lies. Defaults to 'XOY'.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full circle.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the arc. The anchor is the arc center.
name (str, optional): Name of the torus. Defaults to None.
"""
super().__init__(name=name or f'ConcentricArc {self.get_counter()}', dimension=2)
if r_in >= r_out:
raise ValueError('r_in must be smaller than r_out')
outer = Arc(r_out, phi_range, dl, plane)
inner = Arc(r_in, phi_range, dl, plane)
me = Union((inner, outer))
self.set_coord(*_transform_coordinate(me.xs, me.ys, me.zs, plane=plane))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class Circle(Arc):
"""
Full circle as an arc over [0,360).
Shortest import: `from geoparticle import Circle`
"""
[docs] def __init__(self, r: float, dl: float, plane: str='XOY',name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a Circle object.
Args:
r (float): Radius of the circle.
plane (str): Plane in which the circle lies ('XOY', 'YOZ', or 'XOZ').
dl (float): Spacing between points along the circle.
name (str, optional): Name of the circle. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the arc. The anchor is the circle center.
"""
super().__init__(r, '[0,360)', dl, plane,
name=name or f'Circle {self.get_counter()}',
anchor=anchor)
class ThickArc(Geometry):
"""
2D thick arc.
Shortest import: `from geoparticle import ThickArc`
"""
def __init__(self, r_out: float, r_in: float, dl: float, phi_range: str, plane: str = 'XOY',
name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_out (float): Outer radius.
r_in (float): Inner radius.
dl (float): Spacing between points along the arcs.
plane (str, optional): Plane in which the torus lies. Defaults to 'XOY'.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full circle.
name (str, optional): Name of the torus. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the arc. The anchor is the arc center.
"""
super().__init__(name=name or f'ThickArc {self.get_counter()}', dimension=2)
if r_in >= r_out:
raise ValueError('r_in must be smaller than r_out')
n_layers = int((r_out - r_in) / dl) + 1
self.r_in = r_in
self.r_out = r_in + (n_layers - 1) * dl
_check_size_change(r_out, self.r_out, self.name, 'r_out')
layers = []
for i in range(n_layers):
r_layer = r_in + i * dl
layer = Arc(r_layer, phi_range, dl, 'XOY')
layers.append(layer)
me = Union(layers)
self.set_coord(*_transform_coordinate(me.xs, me.ys, me.zs, plane=plane))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
[docs]class ThickRing(ThickArc):
"""
2D thick ring (full circle).
Shortest import: `from geoparticle import ThickRing`
"""
[docs] def __init__(self, r_out: float, r_in: float, dl: float, plane: str = 'XOY', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a ThickRing object.
Args:
r_out (float): Outer radius of the ring.
r_in (float): Inner radius of the ring.
dl (float): Spacing between points along the ring.
plane (str, optional): Plane in which the ring lies. Defaults to 'XOY'.
name (str, optional): Name of the thick ring. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the ring. The anchor is the ring center.
"""
super().__init__(r_out, r_in, dl, phi_range='[0,360)', plane=plane,
name=name or f'ThickRing {self.get_counter()}',
anchor=anchor)
[docs]class FilledCircle(ThickRing):
"""
2D filled circle (ring with r_in=0).
Shortest import: `from geoparticle import FilledCircle`
"""
[docs] def __init__(self, r: float, dl: float, plane: str = 'XOY', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a FilledCircle object.
Args:
r (float): Radius of the circle.
dl (float): Spacing between points within the circle.
plane (str, optional): Plane in which the circle lies. Defaults to 'XOY'.
name (str, optional): Name of the filled circle. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the circle. The anchor is the circle center.
"""
super().__init__(r, 0, dl, plane=plane, anchor=anchor,
name=name or f'FilledCircle {self.get_counter()}')
[docs]class Rectangle(Geometry):
"""
2D rectangle boundary (inner dimensions).
Shortest import: `from geoparticle import Rectangle`
"""
[docs] def __init__(self, length: float, width: float, dl: float, plane: str = 'XOY',
name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a Rectangle object.
Args:
length (float): Length of the rectangle.
width (float): Width of the rectangle.
dl (float): Spacing between points along the rectangle boundary.
plane (str): Plane in which the rectangle lies ('XOY', 'YOZ', or 'XOZ'). Defaults to 'XOY'.
name (str, optional): Name of the rectangle. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the rectangle. The anchor is the vertex
with the smallest (x,y,z) coordinates.
"""
super().__init__(name=name or f'Rectangle {self.get_counter()}', dimension=2)
x_bot = _arange0_quantized(length, dl)
y_left = _arange0_quantized(width - dl, dl) + dl
self.length, self.width = float(x_bot[-1]), float(y_left[-1])
_check_size_change(length, self.length, self.name, 'length')
_check_size_change(width, self.width, self.name, 'width')
# Build boundary without duplicating corners
y_left = y_left[:-1]
x_left = np.full_like(y_left, 0)
y_right = np.copy(y_left)
x_right = np.full_like(y_right, self.length)
y_bot = np.full_like(x_bot, 0)
x_top = np.copy(x_bot)
y_top = np.full_like(x_top, self.width)
xs = np.r_[x_bot, x_left, x_top, x_right]
ys = np.r_[y_bot, y_left, y_top, y_right]
zs = np.zeros_like(xs)
self.set_coord(*_transform_coordinate(xs, ys, zs, plane=plane))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class ThickRectangle(Geometry):
"""
Rectangle wall with thickness inwards/outwards by dl-layers. Inner dims are (length,width).
Shortest import: `from geoparticle import ThickRectangle`
"""
[docs] def __init__(self, length: float, width: float, n_thick: int,
dl: float, plane: str = 'XOY',
name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a ThickRectangle object.
Args:
length (float): Inner length of the rectangle.
width (float): Inner width of the rectangle.
n_thick (int): Number of thickness layers.
dl (float): Spacing between points along the rectangle boundary.
plane (str): Plane in which the rectangle lies ('XOY', 'YOZ', or 'XOZ'). Defaults to 'XOY'.
name (str, optional): Name of the thick rectangle. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the thick rectangle. The anchor is the vertex
with the smallest (x,y,z) coordinates of the innermost layer.
"""
super().__init__(name=name or f'ThickRectangle {self.get_counter()}', dimension=2)
layers = []
for i in range(n_thick):
L = length + 2 * i * dl
W = width + 2 * i * dl
layer = Rectangle(L, W, dl, plane).shift(x=-i * dl, y=-i * dl)
layers.append(layer)
self.load_from(Union(layers))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class FilledRectangle(Geometry):
"""
2D filled rectangle.
Shortest import: `from geoparticle import FilledRectangle`
"""
[docs] def __init__(self, length: float, width: float, dl: float, plane: str = 'XOY',
name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a FilledRectangle object.
Args:
length (float): Length of the rectangle.
width (float): Width of the rectangle.
dl (float): Spacing between points within the rectangle.
plane (str): Plane in which the rectangle lies ('XOY', 'YOZ', or 'XOZ'). Defaults to 'XOY'.
name (str, optional): Name of the filled rectangle. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the rectangle. The anchor is the vertex
with the smallest (x,y,z) coordinates.
"""
super().__init__(name=name or f'FilledRectangle {self.get_counter()}', dimension=2)
x = _arange0_quantized(length, dl)
y = _arange0_quantized(width, dl)
self.length, self.width = float(x[-1]), float(y[-1])
_check_size_change(length, self.length, self.name, 'length')
_check_size_change(width, self.width, self.name, 'width')
X, Y = np.meshgrid(x, y, indexing='xy')
xs, ys = X.ravel(), Y.ravel()
zs = np.zeros_like(xs)
self.set_coord(*_transform_coordinate(xs, ys, zs, plane=plane))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class Block(Geometry):
"""
3D block by stacking a filled rectangle along the z-axis.
Shortest import: `from geoparticle import Block`
"""
[docs] def __init__(self, length: float, width: float, height: float, dl: float, name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Initialize a Block object.
Args:
length (float): Length of the block along the x-axis.
width (float): Width of the block along the y-axis.
height (float): Height of the block along the z-axis.
dl (float): Spacing between points in the grid.
name (str, optional): Name of the block. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the block. The anchor is the vertex
with the smallest (x,y,z) coordinates.
"""
super().__init__(name=name or f'Block {self.get_counter()}', dimension=3)
layer = FilledRectangle(length, width, dl, 'XOY')
n_height = int(height / dl) + 1
self.height = (n_height - 1) * dl
_check_size_change(height, self.height, self.name, 'height')
self.load_from(Stack(layer, 'z', n_height, dl, dimension=3))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class ThickBlockWall(Geometry):
"""
3D thick walls of a rectangular box, including side walls and top/bottom lids.
Shortest import: `from geoparticle import ThickBlockWall`
"""
[docs] def __init__(self, length: float, width: float, height: float, n_thick: int, dl: float,
name=None, anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
length (float): Inner length of the box.
width (float): Inner width of the box.
height (float): Height of the box.
n_thick (int): Number of thickness layers for the walls.
dl (float): Spacing between points in the grid.
name (str, optional): Name of the thick block wall. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the thick block wall. The anchor is the vertex
with the smallest (x,y,z) coordinates of the innermost layer.
"""
super().__init__(name=name or f'ThickBlockWall {self.get_counter()}', dimension=3)
side_layer = ThickRectangle(length, width, n_thick, dl, 'XOY')
n_height = int(height / dl) + 1
side = Stack(side_layer, 'z', n_height + n_thick, dl, dimension=3).shift(z=-(n_thick - 1) * dl)
lid_layer = FilledRectangle(length - 2 * dl, width - 2 * dl, dl, 'XOY').shift(x=dl, y=dl)
lid_lower = Stack(lid_layer, 'z', -n_thick, dl, dimension=3)
z_mid = (n_height - 1) * dl / 2
lid_upper = lid_lower.mirror('XOY', z_mid)
self.load_from(Union((side, lid_lower, lid_upper)))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class CylinderSide(Geometry):
"""
Side surface of a cylinder.
Shortest import: `from geoparticle import CylinderSide`
"""
[docs] def __init__(self, r: float, l_axis: float, dl: float, axis: str = 'y', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r (float): Radius of the cylinder.
l_axis (float): Length of the cylinder along the specified axis.
dl (float): Spacing between points in the grid.
axis (str, optional): Axis along which the cylinder is oriented. Defaults to 'y'.
name (str, optional): Name of the cylinder side. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the cylinder side. The anchor is the
bottom/left/back center of the cylinder.
"""
super().__init__(name=name or f'CylinderSide {self.get_counter()}', dimension=3)
self.l_axis = float(l_axis)
self.r = float(r)
self.radius = float(r)
self.n_axis = int(self.l_axis / dl) + 1
z = np.arange(0, self.n_axis) * dl
self.l_axis = float(z[-1])
_check_size_change(l_axis, self.l_axis, self.name, 'l_axis')
self.n_ring = int(n_per_ring(self.radius, dl))
x_ring, y_ring = _ring_xy(self.n_ring, self.radius)
xs = np.tile(x_ring, self.n_axis)
ys = np.tile(y_ring, self.n_axis)
zs = np.repeat(z, self.n_ring)
self.set_coord(*_transform_coordinate(xs, ys, zs, axis=axis))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
@property
def dl_in_ring(self) -> float:
"""
Get the actual spacing between points on a ring.
Returns:
float: Spacing between points on the ring.
"""
return float(spacing_ring(self.r, self.n_ring))
[docs]class ThickCylinderSide(Geometry):
"""
Thick wall of a cylinder.
Shortest import: `from geoparticle import ThickCylinderSide`
"""
[docs] def __init__(self, r_out, r_in, l_axis: float, dl: float,
axis: str = 'z', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_out (float): Outer radius of the cylinder.
r_in (float): Inner radius of the cylinder.
l_axis (float): Length of the cylinder along the specified axis.
dl (float): Spacing between points in the grid.
axis (str, optional): Axis along which the cylinder is oriented. Defaults to 'y'.
name (str, optional): Name of the thick cylinder wall. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the cylinder side. The anchor is the
bottom/left/back center of the cylinder.
"""
super().__init__(name=name or f'ThickCylinderWall {self.get_counter()}', dimension=3)
n_axis = int(l_axis / dl) + 1
self.l_axis = (n_axis - 1) * dl
_check_size_change(l_axis, self.l_axis, self.name, 'l_axis')
layer = ThickRing(r_out, r_in, dl)
me = Stack(layer, 'z', n_axis, dl, dimension=3)
self.radius = np.sqrt(me.xs ** 2 + me.ys ** 2)
self.set_coord(*_transform_coordinate(me.xs, me.ys, me.zs, axis=axis))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class FilledCylinder(Geometry):
"""
Filled cylinder volume.
Shortest import: `from geoparticle import FilledCylinder`
"""
[docs] def __init__(self, r: float, l_axis: float, dl: float, axis: str = 'y', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r (float): Radius of the cylinder.
l_axis (float): Length of the cylinder along the specified axis.
dl (float): Spacing between points in the grid.
axis (str, optional): Axis along which the cylinder is oriented. Defaults to 'y'.
name (str, optional): Name of the filled cylinder. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the cylinder side. The anchor is the
bottom/left/back center of the cylinder.
"""
super().__init__(name=name or f'FilledCylinder {self.get_counter()}', dimension=3)
radial_layer = FilledCircle(r, dl)
n_axis = int(l_axis / dl) + 1
self.l_axis = (n_axis - 1) * dl
_check_size_change(l_axis, self.l_axis, self.name, 'l_axis')
me = Stack(radial_layer, 'z', n_axis, dl, dimension=3)
self.radius = np.sqrt(me.xs ** 2 + me.ys ** 2)
self.set_coord(*_transform_coordinate(me.xs, me.ys, me.zs, axis=axis))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class TorusSurface(Geometry):
"""
3D torus surface (tube radius r_minor around circle radius r_major).
Shortest import: `from geoparticle import TorusSurface`
"""
[docs] def __init__(self, r_minor: float, r_major: float, dl: float, n_ring: int = None,
plane='XOY', phi_range='[180,270)', regular_id=False, name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_minor (float): Radius of the torus tube.
r_major (float): Radius of the torus centerline.
dl (float): Spacing between points in the grid.
n_ring (int): Number of points along the tube's cross-section.
plane (str, optional): Plane in which the torus lies. Defaults to 'XOY'.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full torus.
regular_id (bool, optional): Whether to use regular indexing for points. Defaults to False.
name (str, optional): Name of the torus. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the center of the major circle.
"""
super().__init__(name=name or f'TorusSurface {self.get_counter()}', dimension=3)
if not (r_minor < r_major):
raise ValueError('r_minor must be smaller than r_major')
if r_minor == 0:
self.load_from(Arc(r_major, phi_range, dl, plane))
return
if n_ring is None:
n_ring = n_per_ring(r_minor, dl)
thetas = np.linspace(0, 2 * np.pi, int(n_ring), endpoint=False) # section angle
a_deg, b_deg, incl_min, incl_max, total_deg = _parse_interval_deg(phi_range)
a = np.deg2rad(a_deg)
b = np.deg2rad(b_deg)
phi_tot = b - a
if regular_id:
n_large = int(n_per_ring(r_major, dl, phi_tot))
d_phi = phi_tot / max(n_large, 1)
phis = d_phi * np.arange(int(not incl_min), max(n_large, 1) + int(incl_max)) + a
all_phi = np.repeat(phis, int(n_ring))
all_theta = np.tile(thetas, phis.size)
else:
r_P = r_major - r_minor * np.cos(thetas) # path radius of section point
nLs = np.maximum(n_per_ring(r_P, dl, phi_tot), 1)
all_phi, all_theta = [], []
for theta, nL in zip(thetas, nLs):
d_phi = phi_tot / nL
cur_phi = d_phi * np.arange(int(not incl_min), nL + int(incl_max)) + a
all_phi.append(cur_phi)
all_theta.append(np.full_like(cur_phi, theta))
all_phi = np.hstack(all_phi) if all_phi else np.array([], dtype=float)
all_theta = np.hstack(all_theta) if all_theta else np.array([], dtype=float)
zs = r_minor * np.sin(all_theta)
ys = (r_major - r_minor * np.cos(all_theta)) * np.sin(all_phi)
xs = (r_major - r_minor * np.cos(all_theta)) * np.cos(all_phi)
self.set_coord(*_transform_coordinate(xs, ys, zs, plane=plane))
self.phi_tot = phi_tot
self.theta = all_theta
self.n_theta = int(n_ring)
self.phi = all_phi
self.n_phi = None if not regular_id else int(n_per_ring(r_major, dl, phi_tot))
self.phi_tot_deg = total_deg
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class ThickTorusWall(Geometry):
"""
3D thick wall of a torus.
Shortest import: `from geoparticle import ThickTorusWall`
"""
[docs] def __init__(self, r_in: float, r_major: float, n_thick: int, dl: float,
plane='XOY', phi_range='[180,270)', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_in (float): Inner radius of the torus tube.
n_thick (int): Number of thickness layers.
r_major (float): Radius of the torus centerline.
dl (float): Spacing between points in the grid.
plane (str, optional): Plane in which the torus lies. Defaults to 'XOY'.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full torus.
name (str, optional): Name of the thick torus wall. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the center of the major circle.
"""
super().__init__(name=name or f'ThickTorusWall {self.get_counter()}', dimension=3)
layers = []
radii = []
for i in range(n_thick):
r_minor = r_in + i * dl
layer = TorusSurface(r_minor, r_major, dl, phi_range=phi_range)
layers.append(layer)
radii.append(np.full_like(layer.xs, r_minor))
self.radius = np.hstack(radii)
me = Union(layers)
self.set_coord(*_transform_coordinate(me.xs, me.ys, me.zs, plane=plane))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class FilledTorus(ThickTorusWall):
"""
3D filled torus volume.
Shortest import: `from geoparticle import FilledTorus`
"""
[docs] def __init__(self, r_minor: float, r_major: float, dl: float,
plane='XOY', phi_range='[180,270)', name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_minor (float): Outer radius of the torus tube.
r_major (float): Radius of the torus centerline.
dl (float): Spacing between points in the grid.
plane (str, optional): Plane in which the torus lies. Defaults to 'XOY'.
phi_range (str, optional): Angular range of the torus in degrees. Defaults to '[180,270)'. Interval
notation should be used, where `[` and `]` denote inclusion, whereas `(` and `)` denote exclusion,
e.g., `[180,270)`, `(0, 90)`, etc. Use `[0,360)` for a full torus.
name (str, optional): Name of the filled torus. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the center of the major circle.
"""
n_thick = int(r_minor / dl) + 1
self.r_minor = (n_thick - 1) * dl
super().__init__(0, r_major, n_thick, dl, plane=plane, phi_range=phi_range,
name=name or f'FilledTorus {self.get_counter()}',
anchor=anchor)
[docs]class SphereSurface(Geometry):
"""
3D sphere surface.
Shortest import: `from geoparticle import SphereSurface`
"""
[docs] def __init__(self, r: float, dl: float, name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r (float): Sphere radius
dl (float): Particle spacing
name (str, optional): Geometry name
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the sphere center.
"""
super().__init__(name=name or f'SphereSurface {self.get_counter()}', dimension=3)
self.r = float(r)
self.dl = float(dl)
if r < 0:
raise ValueError('radius must not be negative')
if dl <= 0:
raise ValueError('particle spacing must be positive')
if r < dl:
self.set_coord(np.array([0]), 0, 0)
return
# Compute number of latitude layers
n_lat = max(1, int(np.pi * r / dl) + 1)
lat_angles = np.linspace(0, np.pi, n_lat)
xs, ys, zs = [], [], []
for i, theta in enumerate(lat_angles):
# Radius of the current latitude ring
r_lat = r * np.sin(theta)
z = r * np.cos(theta) # z coordinate corresponding to this latitude
# Number of points on the current latitude ring
if r_lat < 1e-10: # near the poles
if i == 0 or i == n_lat - 1: # north or south pole
xs.append(0.0)
ys.append(0.0)
zs.append(z if i == 0 else -r) # north or south pole
continue
n_lon = n_per_ring(r_lat, dl)
# Generate angles in the longitude direction
lon_angles = np.linspace(0, 2 * np.pi, n_lon, endpoint=False)
# Compute coordinates on the current latitude ring
x_ring = r_lat * np.cos(lon_angles)
y_ring = r_lat * np.sin(lon_angles)
z_ring = np.full_like(lon_angles, z)
xs.extend(x_ring)
ys.extend(y_ring)
zs.extend(z_ring)
# Set coordinates
self.set_coord(xs, ys, zs)
# Check and adjust pole positions
self._adjust_poles()
# Validate spacing
self._validate_spacing()
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
@property
def radius(self):
return np.sqrt((self.matrix_coords ** 2).sum(axis=1))
def _adjust_poles(self):
"""Adjust pole positions to ensure spacing to neighboring points is close to dl."""
if self.size == 0:
return
coords = self.matrix_coords
# Find potential pole candidates (points with z near ±r)
north_pole_candidates = np.where(coords[:, 2] > self.r - 1e-6)[0]
# Adjust north pole
if len(north_pole_candidates) == 0:
# Add a north pole point
new_coords = np.vstack([coords, [0, 0, self.r]])
self.set_coord(new_coords[:, 0], new_coords[:, 1], new_coords[:, 2])
elif len(north_pole_candidates) > 1:
# Keep the point closest to the north pole, remove others
distances = np.linalg.norm(coords[north_pole_candidates] - [0, 0, self.r], axis=1)
keep_idx = north_pole_candidates[np.argmin(distances)]
mask = np.ones(self.size, dtype=bool)
mask[north_pole_candidates] = False
mask[keep_idx] = True
self.set_coord(self.xs[mask], self.ys[mask], self.zs[mask])
# Adjust south pole (similar handling)
coords = self.matrix_coords # re-fetch coordinates
south_pole_candidates = np.where(coords[:, 2] < -self.r + 1e-6)[0]
if len(south_pole_candidates) == 0:
new_coords = np.vstack([coords, [0, 0, -self.r]])
self.set_coord(new_coords[:, 0], new_coords[:, 1], new_coords[:, 2])
elif len(south_pole_candidates) > 1:
distances = np.linalg.norm(coords[south_pole_candidates] - [0, 0, -self.r], axis=1)
keep_idx = south_pole_candidates[np.argmin(distances)]
mask = np.ones(self.size, dtype=bool)
mask[south_pole_candidates] = False
mask[keep_idx] = True
self.set_coord(self.xs[mask], self.ys[mask], self.zs[mask])
def _validate_spacing(self):
"""Validate spacing of points on the sphere surface."""
from scipy.spatial import KDTree
if self.size < 2:
return
coords = self.matrix_coords
tree = KDTree(coords)
# Check distance to the nearest neighbor for each point
distances, indices = tree.query(coords, k=2)
nearest_distances = distances[:, 1] # exclude self
avg_distance = np.mean(nearest_distances)
min_distance = np.min(nearest_distances)
max_distance = np.max(nearest_distances)
tolerance = 0.2 * self.dl # 20% tolerance
if abs(avg_distance - self.dl) > tolerance:
warn(f'{self.name}: average spacing {avg_distance:.4f} '
f'differs significantly from target spacing {self.dl:.4f}')
if min_distance < self.dl * 0.8:
warn(f'{self.name}: found minimum spacing {min_distance:.4f} '
f'much smaller than target spacing {self.dl:.4f}')
if max_distance > self.dl * 1.2:
warn(f'{self.name}: found maximum spacing {max_distance:.4f} '
f'much larger than target spacing {self.dl:.4f}')
[docs]class ThickSphere(Geometry):
"""
Thick spherical shell between inner and outer spheres.
Shortest import: `from geoparticle import ThickSphere`
"""
[docs] def __init__(self, r_out: float, r_in: float, dl: float, name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r_out (float): Outer radius of the spherical shell.
r_in (float): Inner radius of the spherical shell.
dl (float): Spacing between points in the grid.
name (str, optional): Name of the thick sphere. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the sphere center.
"""
super().__init__(name=name or f'ThickSphere {self.get_counter()}', dimension=3)
if r_out < r_in:
raise ValueError('r_out must be >= r_in')
layers = []
for r in np.arange(r_in, r_out + dl * 0.1, dl):
layer = SphereSurface(r, dl)
layers.append(layer)
self.load_from(Union(layers))
self.rs = np.sqrt((self.matrix_coords ** 2).sum(axis=1))
self.shift(x=anchor[0], y=anchor[1], z=anchor[2], inplace=True)
self.check_overlap()
[docs]class FilledSphere(ThickSphere):
"""
Filled sphere volume.
Shortest import: `from geoparticle import FilledSphere`
"""
[docs] def __init__(self, r: float, dl: float, name=None,
anchor: Sequence[float] = (0, 0, 0)):
"""
Args:
r (float): Radius of the filled sphere.
dl (float): Spacing between points in the grid.
name (str, optional): Name of the filled sphere. Defaults to None.
anchor (Sequence[float]): The anchor coordinate
to determine the absolute position of the torus. The anchor is the sphere center.
"""
super().__init__(
r_out=r, r_in=0, dl=dl,
name=name or f'FilledSphere {self.get_counter()}',
anchor=anchor
)