Particle-based immersed boundaries for pymrm.
Instead of sampling one global signed-distance field, an assembly of
Particle objects generates the immersed-boundary data directly:
each particle classifies the cells it covers (bounding-box window only),
provides the exact wall position on every cut face (analytic where the
shape allows it, root-finding otherwise), and evaluates the exact outward
normal at every wall crossing. This avoids the two error sources of the
union-SDF route near particle contacts:
the union
min_i(level_i)has a gradient kink on the contact medial axis, degrading SDF-gradient normals to O(1) locally (observed to reduce the flux/conjugate IBM from 2nd to ~1st order for touching particles);the fractional wall position
thetainterpolated from the union is polluted when the ghost cell’s union value comes from a different particle.
Particle protocol¶
A particle answers three geometric questions, all in world coordinates:
Particle.level— signed level function (< 0 inside, > 0 outside; approximately a distance near the surface),Particle.intersect— surface crossing on a straight segment whose endpoints straddle the surface,Particle.normal— outward (solid → fluid) unit normal at surface points.
The base class supplies world↔body transforms (position plus an
orientation: an angle in 2-D, a scipy.spatial.transform.Rotation
in 3-D) and numeric defaults for intersect (vectorised bisection) and
normal (finite differences), so a new shape only has to implement the
body-frame level function and bounding box. Sphere (alias
Circle), Box, AnalyticParticle and
GridParticle (B-spline interpolated local samples) are provided.
Contact policy¶
When two particles are closer than one grid cell, the cell-centre
classification alone cannot see the gap: two adjacent cells are solid but
belong to different particles. construct_ibm_particles detects these
contact faces (the segment between the two cell centres leaves one particle
before entering the other) and applies a policy:
"no_flux"(default): a contact crossing is created with each side’s Lagrange reconstruction anchored at its own particle surface; impose a Neumann–Neumann interface condition on these crossings (seecontact_conditions) so there is no direct transport between the particles."merge": no crossing is created; the particles are numerically connected (the union-SDF behaviour)."error": raise, for setups where contact indicates a bug.
Faces where the particles genuinely overlap (the segment does not leave one particle before entering the other) are always merged — interpenetrating particles form one body.
Usage¶
::
particles = [Sphere(c, r) for c, r in zip(centers, radii)]
ibm, info = construct_ibm_particles(particles, x_c)
recon = construct_ibm_normal_derivative(ibm, info.pseudo_sdf, x_c,
normals=info.normals)
ic = contact_conditions(base_ic, ibm, info)
A, g = apply_ibm_interface(L, ibm, recon, ic)Public API¶
| Symbol | Type | Summary |
|---|---|---|
AnalyticParticle | class | Particle from a user-supplied body-frame level function. |
Box | class | Axis-aligned (in body frame) box; rotate via orientation. |
GridParticle | class | Particle from level-function samples on its own body-frame grid. |
Particle | class | Abstract particle: a shape at a position with an orientation. |
ParticleIBMInfo | class | Side-car information produced by construct_ibm_particles. |
Sphere | class | Sphere (any dimension; in 2-D this is a disk — see Circle). |
construct_ibm_particles | function | Build immersed-boundary data directly from a particle assembly. |
contact_conditions | function | Per-crossing ic: base_ic everywhere, a contact condition on contacts. |
AnalyticParticle(level_func, bounding_box, position, orientation = None, normal_func = None)¶
Particle from a user-supplied body-frame level function.
Parameters¶
level_func(callable)level_func(coords) -> valueswithcoordsshaped (..., ndim); negative inside, positive outside, approximately a distance near the surface.bounding_box(tuple of (lo, hi)) Body-frame box containing the particle surface.position, orientation(seeParticle.)normal_func(callable, optional) Body-frame outward normal direction (need not be normalised); default: finite differences of level_func.
Members¶
__init__(level_func, bounding_box, position, orientation = None, normal_func = None)¶
bounding_box(pad = 0.0)¶
World axis-aligned bounding box ((lo, hi), ...) per axis.
bounding_box_body()¶
Body-frame bounding box ((lo, hi), ...) per axis.
intersect(p0, p1)¶
Surface crossing fraction t on the segments p0 → p1.
p0/p1 are (n, ndim) batches whose endpoints straddle the
surface (level(p0) and level(p1) of opposite sign); returns
t in (0, 1) with level(p0 + t (p1 - p0)) == 0. Default:
vectorised bisection on level; shapes with closed-form
intersections override this.
level(coords)¶
Signed level function at world coords shaped (..., ndim).
level_body(coords)¶
Signed level function at body-frame coords shaped (..., ndim).
normal(coords)¶
Outward (solid→fluid) unit normal at world surface points.
normal_body(coords)¶
Gradient direction of level_body (finite differences).
to_body(coords)¶
vec_to_world(vecs)¶
__slots__¶
Box(position, half_extents, orientation = None)¶
Axis-aligned (in body frame) box; rotate via orientation.
half_extents are the half side lengths per axis. Exact SDF and
face normals; segment intersections by the default bisection.
Members¶
__init__(position, half_extents, orientation = None)¶
bounding_box(pad = 0.0)¶
World axis-aligned bounding box ((lo, hi), ...) per axis.
bounding_box_body()¶
Body-frame bounding box ((lo, hi), ...) per axis.
intersect(p0, p1)¶
Surface crossing fraction t on the segments p0 → p1.
p0/p1 are (n, ndim) batches whose endpoints straddle the
surface (level(p0) and level(p1) of opposite sign); returns
t in (0, 1) with level(p0 + t (p1 - p0)) == 0. Default:
vectorised bisection on level; shapes with closed-form
intersections override this.
level(coords)¶
Signed level function at world coords shaped (..., ndim).
level_body(coords)¶
Signed level function at body-frame coords shaped (..., ndim).
normal(coords)¶
Outward (solid→fluid) unit normal at world surface points.
normal_body(coords)¶
Gradient direction of level_body (finite differences).
to_body(coords)¶
vec_to_world(vecs)¶
__slots__¶
GridParticle(values, x_local, position, orientation = None, method = 'cubic')¶
Particle from level-function samples on its own body-frame grid.
The samples are interpolated with a cubic B-spline
(scipy.interpolate.RegularGridInterpolator), so the particle can
be translated and rotated for free. The local grid must extend beyond
the particle surface (positive samples all around); queries outside the
local grid return the clamped boundary value plus the clamping distance,
keeping the sign correct far away.
Parameters¶
values(ndarray) Level-function samples, negative inside.x_local(list of 1-D arrays) Body-frame cell coordinates of the sample grid, one per axis.position, orientation(seeParticle.)method(str, optional) Interpolation method (default"cubic").
Members¶
__init__(values, x_local, position, orientation = None, method = 'cubic')¶
bounding_box(pad = 0.0)¶
World axis-aligned bounding box ((lo, hi), ...) per axis.
bounding_box_body()¶
Body-frame bounding box ((lo, hi), ...) per axis.
intersect(p0, p1)¶
Surface crossing fraction t on the segments p0 → p1.
p0/p1 are (n, ndim) batches whose endpoints straddle the
surface (level(p0) and level(p1) of opposite sign); returns
t in (0, 1) with level(p0 + t (p1 - p0)) == 0. Default:
vectorised bisection on level; shapes with closed-form
intersections override this.
level(coords)¶
Signed level function at world coords shaped (..., ndim).
level_body(coords)¶
Signed level function at body-frame coords shaped (..., ndim).
normal(coords)¶
Outward (solid→fluid) unit normal at world surface points.
normal_body(coords)¶
Gradient direction of level_body (finite differences).
to_body(coords)¶
vec_to_world(vecs)¶
__slots__¶
Particle(position, orientation = None)¶
Abstract particle: a shape at a position with an orientation.
Subclasses implement the body-frame interface
(level_body, bounding_box_body, optionally
normal_body); the world-frame API used by the IBM assembly
(level, normal, intersect,
bounding_box) is provided here, including the numeric fallbacks.
Parameters¶
position(array_like, shape (ndim,)) World position of the particle (body-frame origin).orientation(optional)None(default), an angle in radians (2-D), ascipy.spatial.transform.Rotation(3-D), or an explicit rotation matrix.
Members¶
__init__(position, orientation = None)¶
bounding_box(pad = 0.0)¶
World axis-aligned bounding box ((lo, hi), ...) per axis.
bounding_box_body()¶
Body-frame bounding box ((lo, hi), ...) per axis.
intersect(p0, p1)¶
Surface crossing fraction t on the segments p0 → p1.
p0/p1 are (n, ndim) batches whose endpoints straddle the
surface (level(p0) and level(p1) of opposite sign); returns
t in (0, 1) with level(p0 + t (p1 - p0)) == 0. Default:
vectorised bisection on level; shapes with closed-form
intersections override this.
level(coords)¶
Signed level function at world coords shaped (..., ndim).
level_body(coords)¶
Signed level function at body-frame coords shaped (..., ndim).
normal(coords)¶
Outward (solid→fluid) unit normal at world surface points.
normal_body(coords)¶
Gradient direction of level_body (finite differences).
to_body(coords)¶
vec_to_world(vecs)¶
__slots__¶
ParticleIBMInfo()¶
Side-car information produced by construct_ibm_particles.
Attributes¶
owner(ndarray of int) Per spatial cell, index of the owning particle (deepest level function) or-1for fluid. Shaped like the spatial grid.crossing_particle(ndarray of int, shape (n_crossings,)) Owning particle of the solid (in) side of each crossing.contact(ndarray of bool, shape (n_crossings,)) True for contact crossings (solid–solid faces between two particles under the"no_flux"policy). For these,crossing_particleis thein-side particle andcontact_partnertheout-side one.contact_partner(ndarray of int, shape (n_crossings,)) Theout-side particle of a contact crossing,-1elsewhere.normals(ndarray, shape (n_crossings, ndim)) Exact outward (solid→fluid) unit normals from the owning particle — pass toconstruct_ibm_normal_derivative(..., normals=...). For contact crossings: thein-side particle’s outward normal.pseudo_sdf(ndarray) Sign-correct union level field on the spatial grid (positive filler far from every particle) for region classification in the reconstruction.segmentation(Segmentation) Per-particle labels (owner + 1) — valid even at exact contact, wherepymrm.segment_domainwould merge the bodies.
Members¶
contact¶
contact_partner¶
crossing_particle¶
normals¶
owner¶
pseudo_sdf¶
segmentation¶
Sphere(center, radius)¶
Sphere (any dimension; in 2-D this is a disk — see Circle).
Fully analytic: exact level function, normals, and segment intersections.
Members¶
__init__(center, radius)¶
bounding_box(pad = 0.0)¶
World axis-aligned bounding box ((lo, hi), ...) per axis.
bounding_box_body()¶
Body-frame bounding box ((lo, hi), ...) per axis.
intersect(p0, p1)¶
Surface crossing fraction t on the segments p0 → p1.
p0/p1 are (n, ndim) batches whose endpoints straddle the
surface (level(p0) and level(p1) of opposite sign); returns
t in (0, 1) with level(p0 + t (p1 - p0)) == 0. Default:
vectorised bisection on level; shapes with closed-form
intersections override this.
level(coords)¶
Signed level function at world coords shaped (..., ndim).
level_body(coords)¶
Signed level function at body-frame coords shaped (..., ndim).
normal(coords)¶
Outward (solid→fluid) unit normal at world surface points.
normal_body(coords)¶
Gradient direction of level_body (finite differences).
to_body(coords)¶
vec_to_world(vecs)¶
__slots__¶
construct_ibm_particles(particles, x_c, *, axes = None, shape = None, rescale = True, contact = 'no_flux', halo = 2, fill_value = None)¶
Build immersed-boundary data directly from a particle assembly.
Parameters¶
particles(sequence of Particle)x_c(array_like or list of array_like) Cell-centre coordinates, one 1-D array per spatial axis.axes, shape, rescale(seepymrm.construct_ibm.)contact({‘no_flux’, ‘merge’, ‘error’}, optional) Policy for solid–solid faces between two different particles whose surfaces are separated by a sub-cell gap (see the module docstring). Genuinely overlapping particles are always merged.halo(int, optional) Extra cells around each particle’s bounding box when classifying.fill_value(float, optional)pseudo_sdfvalue for cells not covered by any particle window (default: 4× the largest cell spacing — any positive value works, only the sign is used downstream).
Returns¶
ibm(IBM) Standard immersed-boundary container; all ofpymrm.apply_ibm,pymrm.apply_ibm_interface, etc. apply unchanged.info(ParticleIBMInfo) Ownership, exact normals, contact bookkeeping,pseudo_sdfand a per-particle~pymrm.segmentation.Segmentation.
contact_conditions(base_ic, ibm, info, *, contact_ic = None)¶
Per-crossing ic: base_ic everywhere, a contact condition on contacts.
Parameters¶
base_ic(tuple of two dicts) Interface condition for the regular (fluid–solid) crossings, in the format ofpymrm.apply_ibm_interface. Coefficients may be scalars or ns-broadcastable arrays (no per-crossing arrays).ibm(IBM)info(ParticleIBMInfo)contact_ic(tuple of two dicts, optional) Condition imposed on the contact crossings. Default: independent homogeneous Neumann on both sides (q_out = 0andq_in = 0) — no transport between the particles.
Returns¶
tuple of two dicts with per-crossing coefficient arrays.