rtx-cfd: the embedded-body pieces FSI1 stands on (belongs with b82f307)
EmbeddedMask::traction_at (the per-sample traction factored out of
surface_force, so a coupling loop can load a structure at its own
quadrature points), EmbeddedBody::polygon and the public
polygon_signed_distance (a deformable interface as a vertex list, usable
behind a lock through EmbeddedBody::from_sdf). Left unstaged by mistake
in b82f307 — that commit's FSI1 test needs these to compile.
Co-Authored-By: Claude Fable 5 <[email protected]>
This commit is contained in:
Omar Sobh
co-authored by
Claude Fable 5
Claude Fable 5
parent
b82f307cae
commit
c0f5a86f03
@@ -156,6 +156,61 @@ impl EmbeddedBody {
|
||||
body
|
||||
}
|
||||
|
||||
/// A closed polygon (vertices in order, either winding), at rest. The
|
||||
/// signed distance is exact (min distance to the edges, sign by even-odd
|
||||
/// ray crossing); the sampler walks the edges with outward normals. A
|
||||
/// coupling loop can rebuild the body each subiteration from a deformed
|
||||
/// structure boundary — or share the vertex list behind a lock and let
|
||||
/// the moving-body path pick the new shape up on its per-step rebuild.
|
||||
pub fn polygon(vertices: Vec<(f64, f64)>) -> Self {
|
||||
assert!(vertices.len() >= 3, "a polygon needs at least 3 vertices");
|
||||
// Signed area decides which perpendicular points outward.
|
||||
let signed_area: f64 = vertices
|
||||
.iter()
|
||||
.zip(vertices.iter().cycle().skip(1))
|
||||
.map(|(a, b)| a.0 * b.1 - b.0 * a.1)
|
||||
.take(vertices.len())
|
||||
.sum::<f64>()
|
||||
* 0.5;
|
||||
let ccw = signed_area > 0.0;
|
||||
let sdf_vertices = vertices.clone();
|
||||
let mut body =
|
||||
Self::from_sdf(move |x, y, _| polygon_signed_distance(&sdf_vertices, x, y));
|
||||
let sampler_vertices = vertices;
|
||||
body.sampler = Some(Box::new(move |ds| {
|
||||
let n = sampler_vertices.len();
|
||||
let mut out = Vec::new();
|
||||
for k in 0..n {
|
||||
let (ax, ay) = sampler_vertices[k];
|
||||
let (bx, by) = sampler_vertices[(k + 1) % n];
|
||||
let (ex, ey) = (bx - ax, by - ay);
|
||||
let len = (ex * ex + ey * ey).sqrt();
|
||||
if len == 0.0 {
|
||||
continue;
|
||||
}
|
||||
// Outward normal: right of the direction for CCW winding.
|
||||
let (mut nx, mut ny) = (ey / len, -ex / len);
|
||||
if !ccw {
|
||||
nx = -nx;
|
||||
ny = -ny;
|
||||
}
|
||||
let count = ((len / ds).ceil() as usize).max(1);
|
||||
for q in 0..count {
|
||||
let s = (q as f64 + 0.5) / count as f64;
|
||||
out.push(SurfaceSample {
|
||||
x: ax + s * ex,
|
||||
y: ay + s * ey,
|
||||
nx,
|
||||
ny,
|
||||
ds: len / count as f64,
|
||||
});
|
||||
}
|
||||
}
|
||||
out
|
||||
}));
|
||||
body
|
||||
}
|
||||
|
||||
/// Union of two bodies: the SDF is the minimum; the surface velocity and
|
||||
/// samples come from whichever body a point is closer to. Samples of one
|
||||
/// body lying inside the other are dropped.
|
||||
@@ -236,6 +291,39 @@ pub struct SurfaceForce {
|
||||
pub skipped: usize,
|
||||
}
|
||||
|
||||
|
||||
/// Signed distance to a closed polygon (negative inside, either winding):
|
||||
/// minimum distance over the edges, sign by the even-odd ray-crossing rule.
|
||||
/// Public so a coupling loop can build a time-dependent body from a shared,
|
||||
/// mutating vertex list via [`EmbeddedBody::from_sdf`].
|
||||
#[must_use]
|
||||
pub fn polygon_signed_distance(vertices: &[(f64, f64)], x: f64, y: f64) -> f64 {
|
||||
let n = vertices.len();
|
||||
let mut dist2 = f64::MAX;
|
||||
let mut inside = false;
|
||||
for k in 0..n {
|
||||
let (ax, ay) = vertices[k];
|
||||
let (bx, by) = vertices[(k + 1) % n];
|
||||
let (ex, ey) = (bx - ax, by - ay);
|
||||
let len2 = ex * ex + ey * ey;
|
||||
let s = if len2 > 0.0 {
|
||||
(((x - ax) * ex + (y - ay) * ey) / len2).clamp(0.0, 1.0)
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
let (qx, qy) = (ax + s * ex - x, ay + s * ey - y);
|
||||
dist2 = dist2.min(qx * qx + qy * qy);
|
||||
if (ay > y) != (by > y) {
|
||||
let x_cross = ax + (y - ay) / (by - ay) * ex;
|
||||
if x < x_cross {
|
||||
inside = !inside;
|
||||
}
|
||||
}
|
||||
}
|
||||
let dist = dist2.sqrt();
|
||||
if inside { -dist } else { dist }
|
||||
}
|
||||
|
||||
/// What a velocity face is.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub enum FaceKind {
|
||||
@@ -641,49 +729,17 @@ impl EmbeddedMask {
|
||||
t: f64,
|
||||
ds: f64,
|
||||
) -> SurfaceForce {
|
||||
let h = self.dx.min(self.dy);
|
||||
let samples = body.surface_samples(ds);
|
||||
let (mut fx, mut fy) = (0.0, 0.0);
|
||||
let mut skipped = 0;
|
||||
for s in &samples {
|
||||
let d1 = h;
|
||||
let d2 = 2.0 * h;
|
||||
let probes = (|| {
|
||||
let p1 = self.pressure_at(p, s.x + d1 * s.nx, s.y + d1 * s.ny)?;
|
||||
let p2 = self.pressure_at(p, s.x + d2 * s.nx, s.y + d2 * s.ny)?;
|
||||
let v1 = self.velocity_at(body, u, v, s.x + d1 * s.nx, s.y + d1 * s.ny, t)?;
|
||||
let v2 = self.velocity_at(body, u, v, s.x + d2 * s.nx, s.y + d2 * s.ny, t)?;
|
||||
Some((p1, p2, v1, v2))
|
||||
})();
|
||||
let Some((p1, p2, (u1, v1), (u2, v2))) = probes else {
|
||||
skipped += 1;
|
||||
continue;
|
||||
};
|
||||
let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1);
|
||||
|
||||
let (tx, ty) = (-s.ny, s.nx);
|
||||
let (u_s, v_s) = body.surface_velocity(s.x, s.y, t);
|
||||
let ut_wall = u_s * tx + v_s * ty;
|
||||
let un_wall = u_s * s.nx + v_s * s.ny;
|
||||
// Wall gradient of a quadratic `a s + b s^2` through the two
|
||||
// probes (values relative to the wall).
|
||||
let wall_gradient =
|
||||
|f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1));
|
||||
let dn_ut = wall_gradient(u1 * tx + v1 * ty - ut_wall, u2 * tx + v2 * ty - ut_wall);
|
||||
let dn_un = wall_gradient(
|
||||
u1 * s.nx + v1 * s.ny - un_wall,
|
||||
u2 * s.nx + v2 * s.ny - un_wall,
|
||||
);
|
||||
// Tangential derivative of (u . n) along the surface, n fixed.
|
||||
let eps = 1e-6 * h;
|
||||
let (up, vp) = body.surface_velocity(s.x + eps * tx, s.y + eps * ty, t);
|
||||
let (um, vm) = body.surface_velocity(s.x - eps * tx, s.y - eps * ty, t);
|
||||
let dt_un = ((up - um) * s.nx + (vp - vm) * s.ny) / (2.0 * eps);
|
||||
|
||||
let traction_n = -p_wall + 2.0 * mu * dn_un;
|
||||
let traction_t = mu * (dn_ut + dt_un);
|
||||
fx += (traction_n * s.nx + traction_t * tx) * s.ds;
|
||||
fy += (traction_n * s.ny + traction_t * ty) * s.ds;
|
||||
match self.traction_at(body, u, v, p, mu, t, s.x, s.y, s.nx, s.ny) {
|
||||
Some((tx, ty)) => {
|
||||
fx += tx * s.ds;
|
||||
fy += ty * s.ds;
|
||||
}
|
||||
None => skipped += 1,
|
||||
}
|
||||
}
|
||||
SurfaceForce {
|
||||
fx,
|
||||
@@ -693,7 +749,57 @@ impl EmbeddedMask {
|
||||
}
|
||||
}
|
||||
|
||||
/// Force on the body by a momentum balance over the rectangle of whole
|
||||
/// The reconstructed traction `sigma . n` (force per unit area) at one
|
||||
/// surface point with outward normal `(nx, ny)` — the per-sample core
|
||||
/// of [`Self::surface_force`], exposed so a coupling loop can hand the
|
||||
/// fluid load to a structure at its own quadrature points. `None` when
|
||||
/// a probe cannot be reconstructed (deep concave corner).
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
pub fn traction_at(
|
||||
&self,
|
||||
body: &EmbeddedBody,
|
||||
u: &DMatrix<f64>,
|
||||
v: &DMatrix<f64>,
|
||||
p: &DMatrix<f64>,
|
||||
mu: f64,
|
||||
t: f64,
|
||||
x: f64,
|
||||
y: f64,
|
||||
nx: f64,
|
||||
ny: f64,
|
||||
) -> Option<(f64, f64)> {
|
||||
let h = self.dx.min(self.dy);
|
||||
let d1 = h;
|
||||
let d2 = 2.0 * h;
|
||||
let p1 = self.pressure_at(p, x + d1 * nx, y + d1 * ny)?;
|
||||
let p2 = self.pressure_at(p, x + d2 * nx, y + d2 * ny)?;
|
||||
let (u1, v1) = self.velocity_at(body, u, v, x + d1 * nx, y + d1 * ny, t)?;
|
||||
let (u2, v2) = self.velocity_at(body, u, v, x + d2 * nx, y + d2 * ny, t)?;
|
||||
let p_wall = p1 + (p1 - p2) * d1 / (d2 - d1);
|
||||
|
||||
let (tx, ty) = (-ny, nx);
|
||||
let (u_s, v_s) = body.surface_velocity(x, y, t);
|
||||
let ut_wall = u_s * tx + v_s * ty;
|
||||
let un_wall = u_s * nx + v_s * ny;
|
||||
let wall_gradient =
|
||||
|f1: f64, f2: f64| (f1 * d2 * d2 - f2 * d1 * d1) / (d1 * d2 * (d2 - d1));
|
||||
let dn_ut = wall_gradient(u1 * tx + v1 * ty - ut_wall, u2 * tx + v2 * ty - ut_wall);
|
||||
let dn_un = wall_gradient(u1 * nx + v1 * ny - un_wall, u2 * nx + v2 * ny - un_wall);
|
||||
// Tangential derivative of (u . n) along the surface, n fixed.
|
||||
let eps = 1e-6 * h;
|
||||
let (up, vp) = body.surface_velocity(x + eps * tx, y + eps * ty, t);
|
||||
let (um, vm) = body.surface_velocity(x - eps * tx, y - eps * ty, t);
|
||||
let dt_un = ((up - um) * nx + (vp - vm) * ny) / (2.0 * eps);
|
||||
|
||||
let traction_n = -p_wall + 2.0 * mu * dn_un;
|
||||
let traction_t = mu * (dn_ut + dt_un);
|
||||
Some((
|
||||
traction_n * nx + traction_t * tx,
|
||||
traction_n * ny + traction_t * ty,
|
||||
))
|
||||
}
|
||||
|
||||
/// Force on the body by a momentum balance over the rectangle of whole /// Force on the body by a momentum balance over the rectangle of whole
|
||||
/// cells `[i0, i1) x [j0, j1)` (which must enclose the body and lie in
|
||||
/// the fluid on its boundary):
|
||||
/// `F = sum_outer (sigma.n - rho u (u.n)) A - d/dt int rho u dV + int f dV`,
|
||||
|
||||
@@ -40,7 +40,9 @@ pub use boundary_conditions::{
|
||||
BoundaryCondition, BoundaryConditions, BoundaryLocation, BoundaryType,
|
||||
};
|
||||
pub use embedded::{EmbeddedParameters, EmbeddedPisoSolver, EmbeddedResult};
|
||||
pub use embedded_body::{EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample};
|
||||
pub use embedded_body::{
|
||||
EmbeddedBody, EmbeddedMask, FaceKind, SurfaceForce, SurfaceSample, polygon_signed_distance,
|
||||
};
|
||||
pub use flow_field::FlowField;
|
||||
pub use piso::{PisoParameters, PisoResult, PisoSolver};
|
||||
#[cfg(feature = "cuda")]
|
||||
|
||||
Reference in New Issue
Block a user