R8-c: the 3D flag's interface — deformed-plate body (host + device φ/ub) and the conservative load transfer onto a Hex20 plate
- embedded3/plate.rs: PlateSurface (mid-surface on span stations) and the host evaluation of DeviceSdf (phi_host / velocity_host / is_flag_host), expression for expression e3_geom.cu; a span-uniform plate is the polyline capsule to the bit; first-order spanwise-slope correction. - e3_geom.cu / device/geom.rs: plate branch of geom_phi_at and body_velocity (nst = 0 keeps the polyline path unchanged). - cutwall.rs / exchange.rs: the load loops observed through a sink (sums unchanged); interface.rs: Mask::cut_wall_loads (every summand of cut_wall_force with its foot), HexPlate (R8-b's Hex20 lattice numbering), consistent point-force transfer conserving force and moment to round-off, locate() for the transpose, mid_surface() for the fluid body. - flag test: RTX_E3_FLAG_BODY=plate, _STATIONS, _TWIST, _TRANSFER(_EVERY, _CSV, _NODAL); all default off. Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
This commit is contained in:
Omar Sobh
co-authored by
Claude Opus 5.5
Claude Opus 5.5
parent
d63806c0e6
commit
9c2ae32061
@@ -24,6 +24,7 @@ struct GeomSdf {
|
|||||||
double fillet; /* root fillet radius (0: min) */
|
double fillet; /* root fillet radius (0: min) */
|
||||||
int cyl_cut, flag_cut; /* cut to the span */
|
int cyl_cut, flag_cut; /* cut to the span */
|
||||||
int npts; /* polyline points (x, y interleaved) */
|
int npts; /* polyline points (x, y interleaved) */
|
||||||
|
int nst, ns; /* R8-c plate: stations and points per station (nst 0: the polyline) */
|
||||||
};
|
};
|
||||||
|
|
||||||
struct GeomGrid {
|
struct GeomGrid {
|
||||||
@@ -45,15 +46,107 @@ __device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s
|
|||||||
return outside + rs_min(rs_max(q1, q2), 0.0) - r;
|
return outside + rs_min(rs_max(q1, q2), 0.0) - r;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/*
|
||||||
|
* R8-c: the deformed plate (`plate.rs`, the host twin expression for
|
||||||
|
* expression): `poly` holds the stations' rows (x, y interleaved, row-major
|
||||||
|
* by station) followed by the stations' z. The stations bracketing z are
|
||||||
|
* interpolated linearly into one polyline (outside their range: the end
|
||||||
|
* station as it is), the in-plane closest point taken as the polyline's,
|
||||||
|
* the distance corrected for the spanwise slope d / sqrt(1 + (e·c)²);
|
||||||
|
* with `vel`, the velocity at the closest point the same way.
|
||||||
|
*/
|
||||||
|
__device__ double plate_dist(double x, double y, double z, const GeomSdf& s,
|
||||||
|
const double* __restrict__ P, const double* __restrict__ V,
|
||||||
|
double* vx, double* vy)
|
||||||
|
{
|
||||||
|
int nst = s.nst, ns = s.ns;
|
||||||
|
const double* Z = P + 2 * (long long) nst * ns;
|
||||||
|
int k = 0, interp = 0;
|
||||||
|
double w = 0.0;
|
||||||
|
if (nst > 1) {
|
||||||
|
if (z <= Z[0]) {
|
||||||
|
k = 0;
|
||||||
|
} else if (z >= Z[nst - 1]) {
|
||||||
|
k = nst - 2;
|
||||||
|
} else {
|
||||||
|
while (k + 1 < nst && Z[k + 1] <= z) ++k;
|
||||||
|
if (k > nst - 2) k = nst - 2;
|
||||||
|
}
|
||||||
|
interp = 1;
|
||||||
|
w = (z - Z[k]) / (Z[k + 1] - Z[k]);
|
||||||
|
/* beyond the end stations: at most half an interval extrapolated */
|
||||||
|
if (w < -0.5) w = -0.5;
|
||||||
|
if (w > 1.5) w = 1.5;
|
||||||
|
}
|
||||||
|
const double* R0 = P + 2 * (long long) k * ns;
|
||||||
|
const double* R1 = interp ? R0 + 2 * ns : R0;
|
||||||
|
double best = 1.0 / 0.0, ub = 0.0, qbx = 0.0, qby = 0.0;
|
||||||
|
int mb = 0;
|
||||||
|
for (int m = 0; m + 1 < ns; ++m) {
|
||||||
|
double ax, ay, bx, by;
|
||||||
|
if (interp) {
|
||||||
|
ax = R0[2 * m] + w * (R1[2 * m] - R0[2 * m]);
|
||||||
|
ay = R0[2 * m + 1] + w * (R1[2 * m + 1] - R0[2 * m + 1]);
|
||||||
|
bx = R0[2 * m + 2] + w * (R1[2 * m + 2] - R0[2 * m + 2]);
|
||||||
|
by = R0[2 * m + 3] + w * (R1[2 * m + 3] - R0[2 * m + 3]);
|
||||||
|
} else {
|
||||||
|
ax = R0[2 * m]; ay = R0[2 * m + 1];
|
||||||
|
bx = R0[2 * m + 2]; by = R0[2 * m + 3];
|
||||||
|
}
|
||||||
|
double ex = bx - ax, ey = by - ay;
|
||||||
|
double l2 = ex * ex + ey * ey;
|
||||||
|
double u = ((x - ax) * ex + (y - ay) * ey) / l2;
|
||||||
|
if (u < 0.0) u = 0.0;
|
||||||
|
if (u > 1.0) u = 1.0;
|
||||||
|
double px = ax + u * ex, py = ay + u * ey;
|
||||||
|
double qx = x - px, qy = y - py;
|
||||||
|
double d = sqrt(qx * qx + qy * qy);
|
||||||
|
if (d < best) {
|
||||||
|
best = d;
|
||||||
|
mb = m;
|
||||||
|
ub = u;
|
||||||
|
qbx = qx;
|
||||||
|
qby = qy;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (interp && best > 0.0) {
|
||||||
|
double dz = Z[k + 1] - Z[k];
|
||||||
|
double cax = R1[2 * mb] - R0[2 * mb], cay = R1[2 * mb + 1] - R0[2 * mb + 1];
|
||||||
|
double cbx = R1[2 * mb + 2] - R0[2 * mb + 2], cby = R1[2 * mb + 3] - R0[2 * mb + 3];
|
||||||
|
double cx = (cax + ub * (cbx - cax)) / dz;
|
||||||
|
double cy = (cay + ub * (cby - cay)) / dz;
|
||||||
|
double qn = (qbx * cx + qby * cy) / best;
|
||||||
|
best = best / sqrt(1.0 + qn * qn);
|
||||||
|
}
|
||||||
|
if (V) {
|
||||||
|
const double* V0 = V + 2 * (long long) k * ns;
|
||||||
|
const double* V1 = interp ? V0 + 2 * ns : V0;
|
||||||
|
double avx, avy, bvx, bvy;
|
||||||
|
if (interp) {
|
||||||
|
avx = V0[2 * mb] + w * (V1[2 * mb] - V0[2 * mb]);
|
||||||
|
avy = V0[2 * mb + 1] + w * (V1[2 * mb + 1] - V0[2 * mb + 1]);
|
||||||
|
bvx = V0[2 * mb + 2] + w * (V1[2 * mb + 2] - V0[2 * mb + 2]);
|
||||||
|
bvy = V0[2 * mb + 3] + w * (V1[2 * mb + 3] - V0[2 * mb + 3]);
|
||||||
|
} else {
|
||||||
|
avx = V0[2 * mb]; avy = V0[2 * mb + 1];
|
||||||
|
bvx = V0[2 * mb + 2]; bvy = V0[2 * mb + 3];
|
||||||
|
}
|
||||||
|
*vx = avx + ub * (bvx - avx);
|
||||||
|
*vy = avy + ub * (bvy - avy);
|
||||||
|
}
|
||||||
|
return best;
|
||||||
|
}
|
||||||
|
|
||||||
__device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, const double* __restrict__ poly)
|
__device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, const double* __restrict__ poly)
|
||||||
{
|
{
|
||||||
/* the circle */
|
/* the circle */
|
||||||
double ex0 = x - s.cx, ey0 = y - s.cy;
|
double ex0 = x - s.cx, ey0 = y - s.cy;
|
||||||
double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
|
double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
|
||||||
if (s.cyl_cut) dc = span_cut(dc, z, s);
|
if (s.cyl_cut) dc = span_cut(dc, z, s);
|
||||||
/* the capsule: distance to the polyline */
|
/* the capsule: distance to the polyline (or the plate, R8-c) */
|
||||||
double best = 1.0 / 0.0;
|
double best = 1.0 / 0.0;
|
||||||
for (int m = 0; m + 1 < s.npts; ++m) {
|
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr);
|
||||||
|
else for (int m = 0; m + 1 < s.npts; ++m) {
|
||||||
double ax = poly[2 * m], ay = poly[2 * m + 1];
|
double ax = poly[2 * m], ay = poly[2 * m + 1];
|
||||||
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
|
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
|
||||||
double ex = bx - ax, ey = by - ay;
|
double ex = bx - ax, ey = by - ay;
|
||||||
@@ -299,7 +392,8 @@ __device__ double body_velocity(double x, double y, double z, int c, const GeomS
|
|||||||
const double* __restrict__ poly, const double* __restrict__ vel)
|
const double* __restrict__ poly, const double* __restrict__ vel)
|
||||||
{
|
{
|
||||||
double best = 1.0 / 0.0, vx = 0.0, vy = 0.0;
|
double best = 1.0 / 0.0, vx = 0.0, vy = 0.0;
|
||||||
for (int m = 0; m + 1 < s.npts; ++m) {
|
if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy);
|
||||||
|
else for (int m = 0; m + 1 < s.npts; ++m) {
|
||||||
double ax = poly[2 * m], ay = poly[2 * m + 1];
|
double ax = poly[2 * m], ay = poly[2 * m + 1];
|
||||||
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
|
double bx = poly[2 * m + 2], by = poly[2 * m + 3];
|
||||||
double ex = bx - ax, ey = by - ay;
|
double ex = bx - ax, ey = by - ay;
|
||||||
|
|||||||
@@ -40,6 +40,10 @@ pub struct DeviceSdf {
|
|||||||
/// no z component — the flag test's host closure. Empty: the surface
|
/// no z component — the flag test's host closure. Empty: the surface
|
||||||
/// velocity stays on the host.
|
/// velocity stays on the host.
|
||||||
pub vel: Vec<[f64; 2]>,
|
pub vel: Vec<[f64; 2]>,
|
||||||
|
/// R8-c: the flag as a deformed plate (a mid-surface on span stations,
|
||||||
|
/// `plate.rs`) in place of the polyline: when `Some`, `poly` and `vel`
|
||||||
|
/// are ignored. A span-uniform plate is the polyline capsule to the bit.
|
||||||
|
pub plate: Option<super::plate::PlateSurface>,
|
||||||
}
|
}
|
||||||
|
|
||||||
/// One surface quadrature point: position, unit normal out of the solid,
|
/// One surface quadrature point: position, unit normal out of the solid,
|
||||||
|
|||||||
@@ -908,6 +908,26 @@ impl Mask {
|
|||||||
t: f64,
|
t: f64,
|
||||||
planes: Option<(usize, usize)>,
|
planes: Option<(usize, usize)>,
|
||||||
) -> Option<([f64; 3], [f64; 3])> {
|
) -> Option<([f64; 3], [f64; 3])> {
|
||||||
|
self.cut_wall_force_parts_sink(body, f, mu, t, planes, &mut |_, _, _| {})
|
||||||
|
}
|
||||||
|
|
||||||
|
/// `cut_wall_force_parts_in` with every summand also handed to `sink`
|
||||||
|
/// (R8-c: the per-cell pressure `p_c W_c` at the cell centre, the
|
||||||
|
/// per-face implicit shear at the face position; the sums are computed
|
||||||
|
/// exactly as before — the sink only observes them).
|
||||||
|
pub(super) fn cut_wall_force_parts_sink<S>(
|
||||||
|
&self,
|
||||||
|
body: &Body,
|
||||||
|
f: &Field,
|
||||||
|
mu: f64,
|
||||||
|
t: f64,
|
||||||
|
planes: Option<(usize, usize)>,
|
||||||
|
sink: &mut S,
|
||||||
|
) -> Option<([f64; 3], [f64; 3])>
|
||||||
|
where
|
||||||
|
S: FnMut(super::interface::LoadKind, [f64; 3], [f64; 3]),
|
||||||
|
{
|
||||||
|
use super::interface::LoadKind;
|
||||||
let cut = self.cut.as_ref()?;
|
let cut = self.cut.as_ref()?;
|
||||||
let g = self.grid;
|
let g = self.grid;
|
||||||
let (nx, ny, nz) = (g.nx, g.ny, g.nz);
|
let (nx, ny, nz) = (g.nx, g.ny, g.nz);
|
||||||
@@ -915,12 +935,20 @@ impl Mask {
|
|||||||
let mut pressure = [0.0; 3];
|
let mut pressure = [0.0; 3];
|
||||||
let mut force = [0.0; 3];
|
let mut force = [0.0; 3];
|
||||||
for (idx, w) in cut.wall.iter().enumerate() {
|
for (idx, w) in cut.wall.iter().enumerate() {
|
||||||
let (k, _, i) = g.kji(idx);
|
let (k, j, i) = g.kji(idx);
|
||||||
if self.cell_fluid[idx] && k >= k0 && k < k1 && in_load_window((i as f64 + 0.5) * g.dx)
|
if self.cell_fluid[idx] && k >= k0 && k < k1 && in_load_window((i as f64 + 0.5) * g.dx)
|
||||||
{
|
{
|
||||||
|
let mut v = [0.0; 3];
|
||||||
for c in 0..3 {
|
for c in 0..3 {
|
||||||
pressure[c] += f.p[idx] * w[c];
|
v[c] = f.p[idx] * w[c];
|
||||||
|
pressure[c] += v[c];
|
||||||
}
|
}
|
||||||
|
let x = [
|
||||||
|
(i as f64 + 0.5) * g.dx,
|
||||||
|
(j as f64 + 0.5) * g.dy,
|
||||||
|
(k as f64 + 0.5) * g.dz,
|
||||||
|
];
|
||||||
|
sink(LoadKind::Pressure, x, v);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
let gw = self.gradient_weight_force(&f.p, Some((k0, k1)));
|
let gw = self.gradient_weight_force(&f.p, Some((k0, k1)));
|
||||||
@@ -967,7 +995,11 @@ impl Mask {
|
|||||||
};
|
};
|
||||||
let (c1, c2, nb) = self.wall_gradient(c, p, &cv, xi);
|
let (c1, c2, nb) = self.wall_gradient(c, p, &cv, xi);
|
||||||
let un = nb.map_or(ub, |f| values[c][f]);
|
let un = nb.map_or(ub, |f| values[c][f]);
|
||||||
force[c] += mu * a_w * (c1 * (values[c][idx] - ub) + c2 * (un - ub));
|
let v = mu * a_w * (c1 * (values[c][idx] - ub) + c2 * (un - ub));
|
||||||
|
force[c] += v;
|
||||||
|
let mut fv = [0.0; 3];
|
||||||
|
fv[c] = v;
|
||||||
|
sink(LoadKind::Shear, x, fv);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -92,6 +92,33 @@ impl Mask {
|
|||||||
t: f64,
|
t: f64,
|
||||||
planes: Option<(usize, usize)>,
|
planes: Option<(usize, usize)>,
|
||||||
) -> Option<([f64; 3], [f64; 3])> {
|
) -> Option<([f64; 3], [f64; 3])> {
|
||||||
|
self.cut_wall_exchange_parts_sink(body, f, mu, rho, t, planes, &mut |_, _, _| {})
|
||||||
|
}
|
||||||
|
|
||||||
|
/// `cut_wall_exchange_parts` with every summand also handed to `sink`
|
||||||
|
/// (R8-c: each prescribed neighbour's diffusive and convective exchange
|
||||||
|
/// at the fluid face's position, as a force on the body; the sums are
|
||||||
|
/// computed exactly as before — the sink only observes them).
|
||||||
|
#[allow(clippy::too_many_arguments)]
|
||||||
|
pub(super) fn cut_wall_exchange_parts_sink<S>(
|
||||||
|
&self,
|
||||||
|
body: &Body,
|
||||||
|
f: &Field,
|
||||||
|
mu: f64,
|
||||||
|
rho: f64,
|
||||||
|
t: f64,
|
||||||
|
planes: Option<(usize, usize)>,
|
||||||
|
sink: &mut S,
|
||||||
|
) -> Option<([f64; 3], [f64; 3])>
|
||||||
|
where
|
||||||
|
S: FnMut(super::interface::LoadKind, [f64; 3], [f64; 3]),
|
||||||
|
{
|
||||||
|
use super::interface::LoadKind;
|
||||||
|
let comp = |c: usize, v: f64| {
|
||||||
|
let mut out = [0.0; 3];
|
||||||
|
out[c] = v;
|
||||||
|
out
|
||||||
|
};
|
||||||
let _ = body;
|
let _ = body;
|
||||||
let _ = t;
|
let _ = t;
|
||||||
self.cut.as_ref()?;
|
self.cut.as_ref()?;
|
||||||
@@ -139,6 +166,7 @@ impl Mask {
|
|||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
let cv = self.cv_geometry(c, p);
|
let cv = self.cv_geometry(c, p);
|
||||||
|
let xf = lat.face_position(c, p);
|
||||||
let u0 = vals[c][idx];
|
let u0 = vals[c][idx];
|
||||||
let shift0 = self.face_shift(c, p);
|
let shift0 = self.face_shift(c, p);
|
||||||
let cell_minus = add(p, ec, -1);
|
let cell_minus = add(p, ec, -1);
|
||||||
@@ -192,10 +220,13 @@ impl Mask {
|
|||||||
};
|
};
|
||||||
let u_face = upwind(m_plus, u0, un) + delta;
|
let u_face = upwind(m_plus, u0, un) + delta;
|
||||||
if !self.exchange_convection_off {
|
if !self.exchange_convection_off {
|
||||||
convective[c] -= -rho * m_plus * (u_face - u0);
|
let v = -rho * m_plus * (u_face - u0);
|
||||||
|
convective[c] -= v;
|
||||||
|
sink(LoadKind::ExchangeConvective, xf, comp(c, -v));
|
||||||
}
|
}
|
||||||
force[c] -=
|
let v = mu * cv.ap[d][1] * a_d * (un - u0) / solid_spacing(1.0);
|
||||||
mu * cv.ap[d][1] * a_d * (un - u0) / solid_spacing(1.0);
|
force[c] -= v;
|
||||||
|
sink(LoadKind::ExchangeDiffusive, xf, comp(c, -v));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
// Minus side.
|
// Minus side.
|
||||||
@@ -211,10 +242,14 @@ impl Mask {
|
|||||||
};
|
};
|
||||||
let u_face = upwind(m_minus, ud, u0) + delta;
|
let u_face = upwind(m_minus, ud, u0) + delta;
|
||||||
if !self.exchange_convection_off {
|
if !self.exchange_convection_off {
|
||||||
convective[c] -= rho * m_minus * (u_face - u0);
|
let v = rho * m_minus * (u_face - u0);
|
||||||
|
convective[c] -= v;
|
||||||
|
sink(LoadKind::ExchangeConvective, xf, comp(c, -v));
|
||||||
}
|
}
|
||||||
force[c] -=
|
let v =
|
||||||
mu * cv.ap[d][0] * a_d * (ud - u0) / solid_spacing(-1.0);
|
mu * cv.ap[d][0] * a_d * (ud - u0) / solid_spacing(-1.0);
|
||||||
|
force[c] -= v;
|
||||||
|
sink(LoadKind::ExchangeDiffusive, xf, comp(c, -v));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@@ -0,0 +1,718 @@
|
|||||||
|
//! R8-c: the fluid–structure interface's LOAD side for a 3D flag — the
|
||||||
|
//! cut-cell wall's force as the list of its summands (each located), and
|
||||||
|
//! their consistent, conservative distribution onto a structured Hex20
|
||||||
|
//! plate (the R8-b structure's node layout).
|
||||||
|
//!
|
||||||
|
//! # The loads
|
||||||
|
//!
|
||||||
|
//! [`Mask::cut_wall_loads`] returns every summand of the operator load
|
||||||
|
//! route `cut_wall_force` — the per-cell pressure `p_c W_c` (at the cell
|
||||||
|
//! centre), the per-face implicit wall shear and the per-face wall exchange
|
||||||
|
//! (diffusive and convective; at the fluid face's position) — each with
|
||||||
|
//! its FOOT on the body's surface (`x − φ n`, two projection steps on the
|
||||||
|
//! body's φ). The sums are the route's own: the same loops, observed.
|
||||||
|
//! The route's total is returned with them; the summands' sum differs from
|
||||||
|
//! it by the summation order only (round-off).
|
||||||
|
//!
|
||||||
|
//! # The transfer
|
||||||
|
//!
|
||||||
|
//! [`HexPlate::transfer`] hands each load `(point, F)` to the Hex20 element
|
||||||
|
//! of the deformed plate that contains `point` (the isoparametric map
|
||||||
|
//! inverted by Newton to round-off; the nearest element's extrapolation
|
||||||
|
//! when no element contains it) and distributes it with the element's
|
||||||
|
//! shape functions, `f_a = N_a(ξ) F`: the consistent nodal load of a point
|
||||||
|
//! force (its virtual work). The serendipity functions are a partition of
|
||||||
|
//! unity and reproduce the coordinates (`Σ N_a x_a = x(ξ) = point` once
|
||||||
|
//! Newton has converged), so the total force AND the total moment about
|
||||||
|
//! any point are conserved to round-off, whatever the element. A foot on
|
||||||
|
//! the plate's faces loads that face's nodes only; a foot inside the solid
|
||||||
|
//! (the capsule's rounded span edges and tip lie inside the Hex box) also
|
||||||
|
//! loads the element's other nodes — reported as the interior share.
|
||||||
|
//! The motion side is the transpose: [`HexPlate::locate`] gives the same
|
||||||
|
//! weights for the velocity `Σ N_a v_a` at a fluid point (work-conjugate).
|
||||||
|
//!
|
||||||
|
//! Host only (phase 1): the loads are read from the host mirror of the
|
||||||
|
//! device state at the sample steps, as the load routes are today.
|
||||||
|
|
||||||
|
use super::body::Body;
|
||||||
|
use super::field::Field;
|
||||||
|
use super::plate::PlateSurface;
|
||||||
|
use super::wall::Mask;
|
||||||
|
|
||||||
|
/// Which summand of the operator load route a [`WallLoad`] is.
|
||||||
|
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||||
|
pub enum LoadKind {
|
||||||
|
/// `p_c W_c` of a cut cell.
|
||||||
|
Pressure,
|
||||||
|
/// The implicit wall shear of a fluid face.
|
||||||
|
Shear,
|
||||||
|
/// The diffusive exchange with a prescribed neighbour face.
|
||||||
|
ExchangeDiffusive,
|
||||||
|
/// The convective exchange with a prescribed neighbour face.
|
||||||
|
ExchangeConvective,
|
||||||
|
}
|
||||||
|
|
||||||
|
/// One summand of the cut-cell wall force (a force ON the body).
|
||||||
|
#[derive(Debug, Clone, Copy, PartialEq)]
|
||||||
|
pub struct WallLoad {
|
||||||
|
pub kind: LoadKind,
|
||||||
|
/// Where the operator evaluates it: the cell centre or the face position.
|
||||||
|
pub x: [f64; 3],
|
||||||
|
/// Its foot on the body's surface (the application point).
|
||||||
|
pub foot: [f64; 3],
|
||||||
|
/// The force on the body.
|
||||||
|
pub f: [f64; 3],
|
||||||
|
}
|
||||||
|
|
||||||
|
impl Mask {
|
||||||
|
/// The operator load route's summands (see the module doc) and the
|
||||||
|
/// route's total as `cut_wall_force` computes it. `None` without a cut
|
||||||
|
/// geometry, or with the S2-5 gradient weights on (their force is not
|
||||||
|
/// decomposed).
|
||||||
|
pub fn cut_wall_loads(
|
||||||
|
&self,
|
||||||
|
body: &Body,
|
||||||
|
f: &Field,
|
||||||
|
mu: f64,
|
||||||
|
t: f64,
|
||||||
|
) -> Option<(Vec<WallLoad>, [f64; 3])> {
|
||||||
|
if self.grad_weights.is_some() {
|
||||||
|
return None;
|
||||||
|
}
|
||||||
|
let mut raw: Vec<(LoadKind, [f64; 3], [f64; 3])> = Vec::new();
|
||||||
|
let (p, s) = self
|
||||||
|
.cut_wall_force_parts_sink(body, f, mu, t, None, &mut |k, x, v| raw.push((k, x, v)))?;
|
||||||
|
let (d, c) = self.cut_wall_exchange_parts_sink(
|
||||||
|
body,
|
||||||
|
f,
|
||||||
|
mu,
|
||||||
|
self.density,
|
||||||
|
t,
|
||||||
|
None,
|
||||||
|
&mut |k, x, v| raw.push((k, x, v)),
|
||||||
|
)?;
|
||||||
|
let xch = [d[0] + c[0], d[1] + c[1], d[2] + c[2]];
|
||||||
|
let total = [
|
||||||
|
p[0] + s[0] + xch[0],
|
||||||
|
p[1] + s[1] + xch[1],
|
||||||
|
p[2] + s[2] + xch[2],
|
||||||
|
];
|
||||||
|
let g = self.grid;
|
||||||
|
let eps = 1e-6 * g.dx.min(g.dy).min(g.dz);
|
||||||
|
let foot = |x: [f64; 3]| {
|
||||||
|
let mut q = x;
|
||||||
|
for _ in 0..2 {
|
||||||
|
let s = body.phi(q[0], q[1], q[2], t);
|
||||||
|
let n = body.normal(q[0], q[1], q[2], t, eps);
|
||||||
|
q = [q[0] - s * n.0, q[1] - s * n.1, q[2] - s * n.2];
|
||||||
|
}
|
||||||
|
q
|
||||||
|
};
|
||||||
|
use rayon::prelude::*;
|
||||||
|
let loads = raw
|
||||||
|
.par_iter()
|
||||||
|
.map(|&(kind, x, f)| WallLoad {
|
||||||
|
kind,
|
||||||
|
x,
|
||||||
|
foot: foot(x),
|
||||||
|
f,
|
||||||
|
})
|
||||||
|
.collect();
|
||||||
|
Some((loads, total))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The Hex20 node order of R8-b's `Flag3d` (rtx-fea `analysis/flag3d.rs`):
|
||||||
|
/// the natural coordinates `(ξ, η, ζ)` ↔ lattice `(i, j, k)` = (length,
|
||||||
|
/// thickness, span); corners of `ζ = −1` then `ζ = +1` counter-clockwise
|
||||||
|
/// from `(−1, −1)`, the mid-edges of `ζ = −1`, of `ζ = +1`, then the four
|
||||||
|
/// span edges.
|
||||||
|
const HEX20: [[i8; 3]; 20] = [
|
||||||
|
[-1, -1, -1],
|
||||||
|
[1, -1, -1],
|
||||||
|
[1, 1, -1],
|
||||||
|
[-1, 1, -1],
|
||||||
|
[-1, -1, 1],
|
||||||
|
[1, -1, 1],
|
||||||
|
[1, 1, 1],
|
||||||
|
[-1, 1, 1],
|
||||||
|
[0, -1, -1],
|
||||||
|
[1, 0, -1],
|
||||||
|
[0, 1, -1],
|
||||||
|
[-1, 0, -1],
|
||||||
|
[0, -1, 1],
|
||||||
|
[1, 0, 1],
|
||||||
|
[0, 1, 1],
|
||||||
|
[-1, 0, 1],
|
||||||
|
[-1, -1, 0],
|
||||||
|
[1, -1, 0],
|
||||||
|
[1, 1, 0],
|
||||||
|
[-1, 1, 0],
|
||||||
|
];
|
||||||
|
|
||||||
|
/// The Hex20 serendipity shape functions and their natural derivatives.
|
||||||
|
#[must_use]
|
||||||
|
pub fn hex20(xi: [f64; 3]) -> ([f64; 20], [[f64; 3]; 20]) {
|
||||||
|
let mut n = [0.0; 20];
|
||||||
|
let mut d = [[0.0; 3]; 20];
|
||||||
|
for (a, c) in HEX20.iter().enumerate() {
|
||||||
|
let ca = [f64::from(c[0]), f64::from(c[1]), f64::from(c[2])];
|
||||||
|
let lin = |m: usize| 1.0 + xi[m] * ca[m];
|
||||||
|
if c.iter().all(|&v| v != 0) {
|
||||||
|
let (p, q, r) = (lin(0), lin(1), lin(2));
|
||||||
|
let s = xi[0] * ca[0] + xi[1] * ca[1] + xi[2] * ca[2];
|
||||||
|
n[a] = 0.125 * p * q * r * (s - 2.0);
|
||||||
|
d[a][0] = 0.125 * q * r * ca[0] * (s - 2.0 + p);
|
||||||
|
d[a][1] = 0.125 * p * r * ca[1] * (s - 2.0 + q);
|
||||||
|
d[a][2] = 0.125 * p * q * ca[2] * (s - 2.0 + r);
|
||||||
|
} else {
|
||||||
|
// The zero direction `m0`; the other two linear.
|
||||||
|
let m0 = c.iter().position(|&v| v == 0).expect("mid-edge");
|
||||||
|
let (m1, m2) = ((m0 + 1) % 3, (m0 + 2) % 3);
|
||||||
|
let bub = 1.0 - xi[m0] * xi[m0];
|
||||||
|
n[a] = 0.25 * bub * lin(m1) * lin(m2);
|
||||||
|
d[a][m0] = 0.25 * (-2.0 * xi[m0]) * lin(m1) * lin(m2);
|
||||||
|
d[a][m1] = 0.25 * bub * ca[m1] * lin(m2);
|
||||||
|
d[a][m2] = 0.25 * bub * lin(m1) * ca[m2];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
(n, d)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// A structured Hex20 plate on the `(2nx+1) × (2ny+1) × (2nz+1)`
|
||||||
|
/// serendipity lattice — `nx` elements along the length, `ny` through the
|
||||||
|
/// thickness, `nz` along the span — numbered exactly as R8-b's `Flag3d`:
|
||||||
|
/// lattice points with at most one odd index, scanned `i` (length)
|
||||||
|
/// outermost, then `j` (thickness), then `k` (span); node `n` is the
|
||||||
|
/// `n`-th such point. The elements scan `(ex, ey, ez)` the same way.
|
||||||
|
#[derive(Debug, Clone)]
|
||||||
|
pub struct HexPlate {
|
||||||
|
pub nx: usize,
|
||||||
|
pub ny: usize,
|
||||||
|
pub nz: usize,
|
||||||
|
dims: [usize; 3],
|
||||||
|
lattice: Vec<Option<usize>>,
|
||||||
|
points: Vec<[usize; 3]>,
|
||||||
|
elements: Vec<[usize; 20]>,
|
||||||
|
}
|
||||||
|
|
||||||
|
/// A point located in the plate: its element, natural coordinates, the
|
||||||
|
/// nodes and weights `N_a(ξ)`, the Newton residual `|x(ξ) − p|` and how
|
||||||
|
/// far outside the element it is (`max |ξ| − 1`, ≤ 0 inside).
|
||||||
|
#[derive(Debug, Clone)]
|
||||||
|
pub struct Location {
|
||||||
|
pub element: usize,
|
||||||
|
pub xi: [f64; 3],
|
||||||
|
pub nodes: [usize; 20],
|
||||||
|
pub weights: [f64; 20],
|
||||||
|
pub residual: f64,
|
||||||
|
pub outside: f64,
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The result of one load transfer (see [`HexPlate::transfer`]).
|
||||||
|
#[derive(Debug, Clone)]
|
||||||
|
pub struct PlateTransfer {
|
||||||
|
/// The nodal forces, in the node numbering.
|
||||||
|
pub nodal: Vec<[f64; 3]>,
|
||||||
|
/// Σ F of the loads and Σ f_a of the nodes.
|
||||||
|
pub force_in: [f64; 3],
|
||||||
|
pub force_out: [f64; 3],
|
||||||
|
/// The moments about `origin`: Σ (p − o) × F and Σ (x_a − o) × f_a.
|
||||||
|
pub moment_in: [f64; 3],
|
||||||
|
pub moment_out: [f64; 3],
|
||||||
|
/// The largest Newton residual (m) and the largest outside-ness.
|
||||||
|
pub max_residual: f64,
|
||||||
|
pub max_outside: f64,
|
||||||
|
/// Loads located outside every element (extrapolated), and their Σ|F|.
|
||||||
|
pub extrapolated: usize,
|
||||||
|
pub extrapolated_load: f64,
|
||||||
|
/// The share of Σ|f_a| on nodes off the wetted surface (the interior
|
||||||
|
/// layers and the clamped root face).
|
||||||
|
pub interior_share: f64,
|
||||||
|
}
|
||||||
|
|
||||||
|
fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
|
||||||
|
[
|
||||||
|
a[1] * b[2] - a[2] * b[1],
|
||||||
|
a[2] * b[0] - a[0] * b[2],
|
||||||
|
a[0] * b[1] - a[1] * b[0],
|
||||||
|
]
|
||||||
|
}
|
||||||
|
|
||||||
|
fn norm(a: [f64; 3]) -> f64 {
|
||||||
|
(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt()
|
||||||
|
}
|
||||||
|
|
||||||
|
impl HexPlate {
|
||||||
|
/// The lattice and the elements of an `nx × ny × nz` plate.
|
||||||
|
#[must_use]
|
||||||
|
pub fn new(nx: usize, ny: usize, nz: usize) -> Self {
|
||||||
|
assert!(
|
||||||
|
nx > 0 && ny > 0 && nz > 0,
|
||||||
|
"element counts must be positive"
|
||||||
|
);
|
||||||
|
let dims = [2 * nx + 1, 2 * ny + 1, 2 * nz + 1];
|
||||||
|
let mut lattice = vec![None; dims[0] * dims[1] * dims[2]];
|
||||||
|
let mut points = Vec::new();
|
||||||
|
for i in 0..dims[0] {
|
||||||
|
for j in 0..dims[1] {
|
||||||
|
for k in 0..dims[2] {
|
||||||
|
if (i % 2) + (j % 2) + (k % 2) > 1 {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
lattice[(i * dims[1] + j) * dims[2] + k] = Some(points.len());
|
||||||
|
points.push([i, j, k]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let mut plate = Self {
|
||||||
|
nx,
|
||||||
|
ny,
|
||||||
|
nz,
|
||||||
|
dims,
|
||||||
|
lattice,
|
||||||
|
points,
|
||||||
|
elements: Vec::new(),
|
||||||
|
};
|
||||||
|
for ex in 0..nx {
|
||||||
|
for ey in 0..ny {
|
||||||
|
for ez in 0..nz {
|
||||||
|
let (a, b, c) = (2 * ex, 2 * ey, 2 * ez);
|
||||||
|
let e = HEX20.map(|o| {
|
||||||
|
let at = |base: usize, v: i8| (base as i64 + 1 + i64::from(v)) as usize;
|
||||||
|
plate
|
||||||
|
.lattice_node(at(a, o[0]), at(b, o[1]), at(c, o[2]))
|
||||||
|
.expect("serendipity node")
|
||||||
|
});
|
||||||
|
plate.elements.push(e);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
plate
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The node at lattice `(i, j, k)`, if the point carries one.
|
||||||
|
#[must_use]
|
||||||
|
pub fn lattice_node(&self, i: usize, j: usize, k: usize) -> Option<usize> {
|
||||||
|
if i >= self.dims[0] || j >= self.dims[1] || k >= self.dims[2] {
|
||||||
|
return None;
|
||||||
|
}
|
||||||
|
self.lattice[(i * self.dims[1] + j) * self.dims[2] + k]
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The lattice point of node `n`.
|
||||||
|
#[must_use]
|
||||||
|
pub fn lattice_of(&self, n: usize) -> [usize; 3] {
|
||||||
|
self.points[n]
|
||||||
|
}
|
||||||
|
|
||||||
|
#[must_use]
|
||||||
|
pub fn node_count(&self) -> usize {
|
||||||
|
self.points.len()
|
||||||
|
}
|
||||||
|
|
||||||
|
#[must_use]
|
||||||
|
pub fn elements(&self) -> &[[usize; 20]] {
|
||||||
|
&self.elements
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Whether node `n` lies on the wetted surface (the two faces, the tip,
|
||||||
|
/// the span edges; not the clamped root face `i = 0` unless also on one
|
||||||
|
/// of those).
|
||||||
|
#[must_use]
|
||||||
|
pub fn is_wetted(&self, n: usize) -> bool {
|
||||||
|
let [i, j, k] = self.points[n];
|
||||||
|
j == 0 || j == 2 * self.ny || i == 2 * self.nx || k == 0 || k == 2 * self.nz
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Node positions from a placement `(s, η, ζ) → x` of the lattice's
|
||||||
|
/// fractions: `s = i/(2nx) ∈ [0, 1]` along the length, `η = j/ny − 1 ∈
|
||||||
|
/// [−1, 1]` through the thickness, `ζ = k/(2nz) ∈ [0, 1]` along the span.
|
||||||
|
#[must_use]
|
||||||
|
pub fn place<F: Fn(f64, f64, f64) -> [f64; 3]>(&self, f: F) -> Vec<[f64; 3]> {
|
||||||
|
self.points
|
||||||
|
.iter()
|
||||||
|
.map(|&[i, j, k]| {
|
||||||
|
f(
|
||||||
|
i as f64 / (2 * self.nx) as f64,
|
||||||
|
j as f64 / self.ny as f64 - 1.0,
|
||||||
|
k as f64 / (2 * self.nz) as f64,
|
||||||
|
)
|
||||||
|
})
|
||||||
|
.collect()
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Newton on the isoparametric map of element `e` for the point `p`.
|
||||||
|
fn invert(&self, pos: &[[f64; 3]], e: usize, p: [f64; 3]) -> Location {
|
||||||
|
let nodes = self.elements[e];
|
||||||
|
let x: Vec<[f64; 3]> = nodes.iter().map(|&n| pos[n]).collect();
|
||||||
|
let mut xi = [0.0f64; 3];
|
||||||
|
let mut residual = f64::INFINITY;
|
||||||
|
let mut weights = [0.0; 20];
|
||||||
|
for it in 0..60 {
|
||||||
|
let (n, d) = hex20(xi);
|
||||||
|
let mut r = [-p[0], -p[1], -p[2]];
|
||||||
|
let mut jac = [[0.0f64; 3]; 3];
|
||||||
|
for a in 0..20 {
|
||||||
|
for c in 0..3 {
|
||||||
|
r[c] += n[a] * x[a][c];
|
||||||
|
for m in 0..3 {
|
||||||
|
jac[c][m] += x[a][c] * d[a][m];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
weights = n;
|
||||||
|
let rn = norm(r);
|
||||||
|
// Converged: two more iterations past the first residual at
|
||||||
|
// round-off level make the last one a no-op.
|
||||||
|
if rn <= residual && rn < 1e-15 && it > 2 {
|
||||||
|
residual = rn;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
residual = rn;
|
||||||
|
// δ = −J⁻¹ r by the adjugate.
|
||||||
|
let det = jac[0][0] * (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1])
|
||||||
|
- jac[0][1] * (jac[1][0] * jac[2][2] - jac[1][2] * jac[2][0])
|
||||||
|
+ jac[0][2] * (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]);
|
||||||
|
if det == 0.0 || !det.is_finite() {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
let inv = [
|
||||||
|
[
|
||||||
|
jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1],
|
||||||
|
jac[0][2] * jac[2][1] - jac[0][1] * jac[2][2],
|
||||||
|
jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1],
|
||||||
|
],
|
||||||
|
[
|
||||||
|
jac[1][2] * jac[2][0] - jac[1][0] * jac[2][2],
|
||||||
|
jac[0][0] * jac[2][2] - jac[0][2] * jac[2][0],
|
||||||
|
jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2],
|
||||||
|
],
|
||||||
|
[
|
||||||
|
jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0],
|
||||||
|
jac[0][1] * jac[2][0] - jac[0][0] * jac[2][1],
|
||||||
|
jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0],
|
||||||
|
],
|
||||||
|
];
|
||||||
|
for m in 0..3 {
|
||||||
|
let dm = (inv[m][0] * r[0] + inv[m][1] * r[1] + inv[m][2] * r[2]) / det;
|
||||||
|
// Keep the iterate bounded (a far point in a distorted element).
|
||||||
|
xi[m] = (xi[m] - dm).clamp(-4.0, 4.0);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let outside = xi.iter().fold(f64::NEG_INFINITY, |m, v| m.max(v.abs())) - 1.0;
|
||||||
|
Location {
|
||||||
|
element: e,
|
||||||
|
xi,
|
||||||
|
nodes,
|
||||||
|
weights,
|
||||||
|
residual,
|
||||||
|
outside,
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Locate `p` in the deformed plate `pos`: the containing element
|
||||||
|
/// (Newton on the elements nearest by centroid), else the nearest
|
||||||
|
/// element's extrapolation (the smallest outside-ness).
|
||||||
|
#[must_use]
|
||||||
|
pub fn locate(&self, pos: &[[f64; 3]], centroids: &[[f64; 3]], p: [f64; 3]) -> Location {
|
||||||
|
let mut order: Vec<(f64, usize)> = centroids
|
||||||
|
.iter()
|
||||||
|
.enumerate()
|
||||||
|
.map(|(e, c)| {
|
||||||
|
let d = [c[0] - p[0], c[1] - p[1], c[2] - p[2]];
|
||||||
|
(d[0] * d[0] + d[1] * d[1] + d[2] * d[2], e)
|
||||||
|
})
|
||||||
|
.collect();
|
||||||
|
let take = 8.min(order.len());
|
||||||
|
order.select_nth_unstable_by(take - 1, |a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
|
||||||
|
order[..take].sort_by(|a, b| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1)));
|
||||||
|
let mut best: Option<Location> = None;
|
||||||
|
for &(_, e) in &order[..take] {
|
||||||
|
let loc = self.invert(pos, e, p);
|
||||||
|
if loc.outside <= 1e-9 && loc.residual < 1e-12 {
|
||||||
|
return loc;
|
||||||
|
}
|
||||||
|
if best
|
||||||
|
.as_ref()
|
||||||
|
.is_none_or(|b| (loc.outside, loc.residual) < (b.outside, b.residual))
|
||||||
|
{
|
||||||
|
best = Some(loc);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
best.expect("an element")
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The elements' corner centroids in the deformed plate.
|
||||||
|
#[must_use]
|
||||||
|
pub fn centroids(&self, pos: &[[f64; 3]]) -> Vec<[f64; 3]> {
|
||||||
|
self.elements
|
||||||
|
.iter()
|
||||||
|
.map(|e| {
|
||||||
|
let mut c = [0.0; 3];
|
||||||
|
for &n in &e[..8] {
|
||||||
|
for m in 0..3 {
|
||||||
|
c[m] += 0.125 * pos[n][m];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
c
|
||||||
|
})
|
||||||
|
.collect()
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The consistent nodal load of point forces `(p, F)` on the deformed
|
||||||
|
/// plate `pos` (see the module doc), with the conservation budget about
|
||||||
|
/// `origin`.
|
||||||
|
#[must_use]
|
||||||
|
pub fn transfer(
|
||||||
|
&self,
|
||||||
|
pos: &[[f64; 3]],
|
||||||
|
loads: &[([f64; 3], [f64; 3])],
|
||||||
|
origin: [f64; 3],
|
||||||
|
) -> PlateTransfer {
|
||||||
|
use rayon::prelude::*;
|
||||||
|
assert_eq!(pos.len(), self.node_count(), "one position per node");
|
||||||
|
let centroids = self.centroids(pos);
|
||||||
|
let located: Vec<Location> = loads
|
||||||
|
.par_iter()
|
||||||
|
.map(|&(p, _)| self.locate(pos, ¢roids, p))
|
||||||
|
.collect();
|
||||||
|
let mut nodal = vec![[0.0f64; 3]; self.node_count()];
|
||||||
|
let (mut force_in, mut moment_in) = ([0.0f64; 3], [0.0f64; 3]);
|
||||||
|
let (mut max_residual, mut max_outside) = (0.0f64, f64::NEG_INFINITY);
|
||||||
|
let (mut extrapolated, mut extrapolated_load) = (0usize, 0.0f64);
|
||||||
|
for (&(p, f), loc) in loads.iter().zip(&located) {
|
||||||
|
for (&n, &w) in loc.nodes.iter().zip(&loc.weights) {
|
||||||
|
for c in 0..3 {
|
||||||
|
nodal[n][c] += w * f[c];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let m = cross([p[0] - origin[0], p[1] - origin[1], p[2] - origin[2]], f);
|
||||||
|
for c in 0..3 {
|
||||||
|
force_in[c] += f[c];
|
||||||
|
moment_in[c] += m[c];
|
||||||
|
}
|
||||||
|
max_residual = max_residual.max(loc.residual);
|
||||||
|
max_outside = max_outside.max(loc.outside);
|
||||||
|
if loc.outside > 1e-9 {
|
||||||
|
extrapolated += 1;
|
||||||
|
extrapolated_load += norm(f);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let (mut force_out, mut moment_out) = ([0.0f64; 3], [0.0f64; 3]);
|
||||||
|
let (mut total_abs, mut interior_abs) = (0.0f64, 0.0f64);
|
||||||
|
for (n, f) in nodal.iter().enumerate() {
|
||||||
|
let x = pos[n];
|
||||||
|
let m = cross([x[0] - origin[0], x[1] - origin[1], x[2] - origin[2]], *f);
|
||||||
|
for c in 0..3 {
|
||||||
|
force_out[c] += f[c];
|
||||||
|
moment_out[c] += m[c];
|
||||||
|
}
|
||||||
|
let a = norm(*f);
|
||||||
|
total_abs += a;
|
||||||
|
if !self.is_wetted(n) || self.points[n][0] == 0 {
|
||||||
|
interior_abs += a;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
PlateTransfer {
|
||||||
|
nodal,
|
||||||
|
force_in,
|
||||||
|
force_out,
|
||||||
|
moment_in,
|
||||||
|
moment_out,
|
||||||
|
max_residual,
|
||||||
|
max_outside,
|
||||||
|
extrapolated,
|
||||||
|
extrapolated_load,
|
||||||
|
interior_share: if total_abs > 0.0 {
|
||||||
|
interior_abs / total_abs
|
||||||
|
} else {
|
||||||
|
0.0
|
||||||
|
},
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The mid-surface of the deformed plate for the fluid's body
|
||||||
|
/// ([`PlateSurface`]): the lattice's middle layer `j = ny` (`ny` even:
|
||||||
|
/// a node layer), the stations at the element boundaries `k` even (their
|
||||||
|
/// z the nodes' at the root, `i = 0`), the points along the length every
|
||||||
|
/// lattice step; the last point of each station pulled back along its
|
||||||
|
/// last segment by `tip_inset` (the capsule's apex on the structure's
|
||||||
|
/// tip: the flag test's `RTX_E3_FLAG_TIP_INSET`, one half-thickness).
|
||||||
|
/// With nodal velocities, the matching station velocities.
|
||||||
|
#[must_use]
|
||||||
|
pub fn mid_surface(
|
||||||
|
&self,
|
||||||
|
pos: &[[f64; 3]],
|
||||||
|
vel: Option<&[[f64; 3]]>,
|
||||||
|
tip_inset: f64,
|
||||||
|
) -> PlateSurface {
|
||||||
|
assert!(
|
||||||
|
self.ny % 2 == 0,
|
||||||
|
"the mid-surface is a node layer for an even thickness count"
|
||||||
|
);
|
||||||
|
let j = self.ny;
|
||||||
|
let ns = self.dims[0];
|
||||||
|
let mut z = Vec::new();
|
||||||
|
let mut xy = Vec::new();
|
||||||
|
let mut vv = Vec::new();
|
||||||
|
for k in (0..self.dims[2]).step_by(2) {
|
||||||
|
z.push(pos[self.lattice_node(0, j, k).expect("root node")][2]);
|
||||||
|
let row0 = xy.len();
|
||||||
|
for i in 0..ns {
|
||||||
|
let n = self.lattice_node(i, j, k).expect("mid-layer node");
|
||||||
|
xy.push([pos[n][0], pos[n][1]]);
|
||||||
|
if let Some(v) = vel {
|
||||||
|
vv.push([v[n][0], v[n][1]]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if tip_inset > 0.0 {
|
||||||
|
let [ax, ay] = xy[row0 + ns - 2];
|
||||||
|
let [bx, by] = xy[row0 + ns - 1];
|
||||||
|
let len = ((bx - ax).powi(2) + (by - ay).powi(2)).sqrt();
|
||||||
|
let f = (1.0 - tip_inset / len).max(0.0);
|
||||||
|
xy[row0 + ns - 1] = [ax + f * (bx - ax), ay + f * (by - ay)];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
PlateSurface { z, ns, xy, vel: vv }
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#[cfg(test)]
|
||||||
|
mod tests {
|
||||||
|
use super::*;
|
||||||
|
|
||||||
|
#[test]
|
||||||
|
fn hex20_is_a_partition_of_unity_with_nodal_interpolation() {
|
||||||
|
for (a, c) in HEX20.iter().enumerate() {
|
||||||
|
let (n, _) = hex20([f64::from(c[0]), f64::from(c[1]), f64::from(c[2])]);
|
||||||
|
for (b, v) in n.iter().enumerate() {
|
||||||
|
assert!((v - if a == b { 1.0 } else { 0.0 }).abs() < 1e-15);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let (n, d) = hex20([0.3, -0.7, 0.45]);
|
||||||
|
assert!((n.iter().sum::<f64>() - 1.0).abs() < 1e-15);
|
||||||
|
for m in 0..3 {
|
||||||
|
assert!(d.iter().map(|v| v[m]).sum::<f64>().abs() < 1e-14);
|
||||||
|
}
|
||||||
|
// Derivatives against central differences.
|
||||||
|
let h = 1e-6;
|
||||||
|
for m in 0..3 {
|
||||||
|
let mut p = [0.3, -0.7, 0.45];
|
||||||
|
let mut q = p;
|
||||||
|
p[m] += h;
|
||||||
|
q[m] -= h;
|
||||||
|
let (np, _) = hex20(p);
|
||||||
|
let (nq, _) = hex20(q);
|
||||||
|
for a in 0..20 {
|
||||||
|
assert!(((np[a] - nq[a]) / (2.0 * h) - d[a][m]).abs() < 1e-8);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#[test]
|
||||||
|
fn layout_matches_the_structure() {
|
||||||
|
let p = HexPlate::new(35, 2, 10);
|
||||||
|
// Serendipity count: (2nx+1)(2ny+1)(2nz+1) minus the points with ≥ 2 odd.
|
||||||
|
let full = 71 * 5 * 21;
|
||||||
|
let two_odd = 35 * 2 * 21 + 35 * 5 * 10 + 71 * 2 * 10 - 2 * 35 * 2 * 10;
|
||||||
|
assert_eq!(p.node_count(), full - two_odd);
|
||||||
|
assert_eq!(p.elements().len(), 35 * 2 * 10);
|
||||||
|
// Node 0 is (0,0,0), the next along the span.
|
||||||
|
assert_eq!(p.lattice_of(0), [0, 0, 0]);
|
||||||
|
assert_eq!(p.lattice_of(1), [0, 0, 1]);
|
||||||
|
assert_eq!(p.elements()[0][1], p.lattice_node(2, 0, 0).unwrap());
|
||||||
|
}
|
||||||
|
|
||||||
|
/// A bent, twisted plate: arbitrary point loads inside and just outside
|
||||||
|
/// are distributed with the total force and moment conserved to
|
||||||
|
/// round-off.
|
||||||
|
#[test]
|
||||||
|
fn transfer_conserves_force_and_moment() {
|
||||||
|
let plate = HexPlate::new(35, 2, 10);
|
||||||
|
let bend = |s: f64, zeta: f64| 0.06 * s * s * (1.0 + 0.5 * (2.0 * zeta - 1.0));
|
||||||
|
let pos = plate.place(|s, eta, zeta| {
|
||||||
|
let x = 0.25 + 0.35 * s;
|
||||||
|
let y = 0.2 + bend(s, zeta) + 0.01 * eta;
|
||||||
|
[x, y, 0.105 + 0.2 * zeta]
|
||||||
|
});
|
||||||
|
let mut loads = Vec::new();
|
||||||
|
let mut state = 12345u64;
|
||||||
|
let mut rnd = || {
|
||||||
|
state = state
|
||||||
|
.wrapping_mul(6364136223846793005)
|
||||||
|
.wrapping_add(1442695040888963407);
|
||||||
|
(state >> 11) as f64 / (1u64 << 53) as f64
|
||||||
|
};
|
||||||
|
for _ in 0..3000 {
|
||||||
|
let (s, zeta, eta) = (rnd(), rnd(), 2.4 * rnd() - 1.2);
|
||||||
|
let p = [
|
||||||
|
0.25 + 0.35 * s,
|
||||||
|
0.2 + bend(s, zeta) + 0.01 * eta,
|
||||||
|
0.105 + 0.2 * zeta,
|
||||||
|
];
|
||||||
|
loads.push((p, [rnd() - 0.5, 10.0 * (rnd() - 0.5), 0.1 * (rnd() - 0.5)]));
|
||||||
|
}
|
||||||
|
let tr = plate.transfer(&pos, &loads, [0.25, 0.2, 0.205]);
|
||||||
|
let scale_f = loads.iter().map(|l| norm(l.1)).sum::<f64>();
|
||||||
|
for c in 0..3 {
|
||||||
|
assert!((tr.force_in[c] - tr.force_out[c]).abs() < 1e-13 * scale_f);
|
||||||
|
assert!((tr.moment_in[c] - tr.moment_out[c]).abs() < 1e-13 * scale_f);
|
||||||
|
}
|
||||||
|
assert!(tr.max_residual < 1e-14);
|
||||||
|
assert!(tr.extrapolated > 0, "the ±1.2 band reaches outside");
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The mid-surface of the deformed Hex plate, thickened by the half
|
||||||
|
/// thickness, passes through the plate's face nodes (φ ≈ 0 there): the
|
||||||
|
/// motion side's geometry and the structure agree.
|
||||||
|
#[test]
|
||||||
|
fn mid_surface_passes_the_face_nodes() {
|
||||||
|
use super::super::body::DeviceSdf;
|
||||||
|
let plate = HexPlate::new(35, 2, 10);
|
||||||
|
let half = 0.01;
|
||||||
|
let bend = |s: f64, zeta: f64| 0.08 * s * s * (1.0 + 0.3 * (2.0 * zeta - 1.0));
|
||||||
|
let pos = plate.place(|s, eta, zeta| {
|
||||||
|
// Offsets along the 3D normal of the mid-surface y = w(x, z).
|
||||||
|
let ds = 1e-7;
|
||||||
|
let (x, y) = (0.25 + 0.35 * s, 0.2 + bend(s, zeta));
|
||||||
|
let wx = (bend(s + ds, zeta) - bend(s - ds, zeta)) / (2.0 * ds) / 0.35;
|
||||||
|
let wz = (bend(s, zeta + ds) - bend(s, zeta - ds)) / (2.0 * ds) / 0.2;
|
||||||
|
let r = (1.0 + wx * wx + wz * wz).sqrt();
|
||||||
|
[
|
||||||
|
x - half * eta * wx / r,
|
||||||
|
y + half * eta / r,
|
||||||
|
0.105 + 0.2 * zeta - half * eta * wz / r,
|
||||||
|
]
|
||||||
|
});
|
||||||
|
let surf = plate.mid_surface(&pos, None, 0.0);
|
||||||
|
surf.validate().unwrap();
|
||||||
|
let sdf = DeviceSdf {
|
||||||
|
cyl: [-10.0, -10.0, 0.05],
|
||||||
|
cyl_cut: false,
|
||||||
|
flag_cut: false,
|
||||||
|
zc: 0.205,
|
||||||
|
span: 0.2,
|
||||||
|
r_edge: 0.0,
|
||||||
|
half,
|
||||||
|
fillet: 0.0,
|
||||||
|
poly: Vec::new(),
|
||||||
|
vel: Vec::new(),
|
||||||
|
plate: Some(surf),
|
||||||
|
};
|
||||||
|
let mut worst = 0.0f64;
|
||||||
|
for n in 0..plate.node_count() {
|
||||||
|
let [i, j, _] = plate.lattice_of(n);
|
||||||
|
if (j == 0 || j == 4) && i > 1 && i < 69 {
|
||||||
|
let p = pos[n];
|
||||||
|
worst = worst.max(sdf.phi_host(p[0], p[1], p[2]).abs());
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// Measured 4.5e-7 m (the mid-surface's 5 mm chords, the first-order
|
||||||
|
// slope correction) against the 10 mm half-thickness.
|
||||||
|
assert!(
|
||||||
|
worst < 2e-6,
|
||||||
|
"face nodes off the level set by {worst:.3e} m"
|
||||||
|
);
|
||||||
|
}
|
||||||
|
}
|
||||||
@@ -15,8 +15,10 @@ pub mod export_vtk;
|
|||||||
pub mod field;
|
pub mod field;
|
||||||
pub mod grid;
|
pub mod grid;
|
||||||
pub mod impose;
|
pub mod impose;
|
||||||
|
pub mod interface;
|
||||||
pub mod loads;
|
pub mod loads;
|
||||||
pub mod maskupdate;
|
pub mod maskupdate;
|
||||||
|
pub mod plate;
|
||||||
pub mod poisson;
|
pub mod poisson;
|
||||||
pub mod reconstruct;
|
pub mod reconstruct;
|
||||||
pub mod step;
|
pub mod step;
|
||||||
@@ -27,6 +29,8 @@ pub use cut::CutGeometry;
|
|||||||
pub use export_vtk::write_vtk;
|
pub use export_vtk::write_vtk;
|
||||||
pub use field::Field;
|
pub use field::Field;
|
||||||
pub use grid::Grid;
|
pub use grid::Grid;
|
||||||
|
pub use interface::{HexPlate, LoadKind, PlateTransfer, WallLoad};
|
||||||
pub use loads::SurfaceForce;
|
pub use loads::SurfaceForce;
|
||||||
|
pub use plate::PlateSurface;
|
||||||
pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult};
|
pub use step::{Boundaries, Fluid, Parameters, Side, Solver, StepResult};
|
||||||
pub use wall::{FaceKind, Mask, WallScheme};
|
pub use wall::{FaceKind, Mask, WallScheme};
|
||||||
|
|||||||
@@ -0,0 +1,414 @@
|
|||||||
|
//! R8-c: the flag as a DEFORMED PLATE — a mid-surface given on a structured
|
||||||
|
//! grid of points `(s_m, z_k) → (x, y)` at span stations `z_k`, thickened by
|
||||||
|
//! the capsule's half-thickness — and the host evaluation of a
|
||||||
|
//! [`DeviceSdf`], expression for expression the device kernel's
|
||||||
|
//! (`e3_geom.cu` `geom_phi_at` / `body_velocity`), for both the polyline
|
||||||
|
//! capsule (R6-1) and the plate.
|
||||||
|
//!
|
||||||
|
//! The plate's φ at a point `(x, y, z)`: the stations bracketing `z` are
|
||||||
|
//! interpolated linearly in `z` into one centreline polyline, the capsule
|
||||||
|
//! distance to that polyline is taken in the plane `z = const` exactly as
|
||||||
|
//! the polyline capsule does, and that in-plane distance `d` is corrected
|
||||||
|
//! to first order for the mid-surface's spanwise slope: with `c = ∂C/∂z`
|
||||||
|
//! the centreline's rate along the span at the closest point and `e` the
|
||||||
|
//! unit in-plane offset, the distance to the ruled surface's tangent plane
|
||||||
|
//! is `d / sqrt(1 + (e·c)²)`. Beyond the end stations the end interval is
|
||||||
|
//! extrapolated by at most half an interval. When the stations are identical (a
|
||||||
|
//! span-uniform mid-surface) every interpolation is `a + w·0 = a` and the
|
||||||
|
//! slope `c` is exactly zero, so φ is the polyline capsule's to the bit —
|
||||||
|
//! the G1 identity. The span cut (the flag's span edges rounded to
|
||||||
|
//! `r_edge`) is the polyline capsule's: the stations carry no spanwise
|
||||||
|
//! displacement (phase 1; a structure's `u_z` is not represented).
|
||||||
|
//!
|
||||||
|
//! The surface velocity: the stations' velocities interpolated the same
|
||||||
|
//! way at the closest point (no z component, as the polyline capsule's).
|
||||||
|
|
||||||
|
use super::body::DeviceSdf;
|
||||||
|
|
||||||
|
/// The deformed plate's mid-surface: `z.len()` span stations (ascending),
|
||||||
|
/// each a centreline polyline of `ns` points `(x, y)`; the velocities
|
||||||
|
/// `(vx, vy)` per point (empty: the surface velocity stays on the host).
|
||||||
|
#[derive(Debug, Clone, PartialEq)]
|
||||||
|
pub struct PlateSurface {
|
||||||
|
/// The stations' span coordinates, ascending (at least one).
|
||||||
|
pub z: Vec<f64>,
|
||||||
|
/// Points per station (at least two).
|
||||||
|
pub ns: usize,
|
||||||
|
/// Row-major by station: point `m` of station `k` is `xy[k·ns + m]`.
|
||||||
|
pub xy: Vec<[f64; 2]>,
|
||||||
|
/// The points' velocities, the same layout (or empty).
|
||||||
|
pub vel: Vec<[f64; 2]>,
|
||||||
|
}
|
||||||
|
|
||||||
|
impl PlateSurface {
|
||||||
|
/// A span-uniform plate: the polyline `row` (and its velocities) at
|
||||||
|
/// every station `z`.
|
||||||
|
#[must_use]
|
||||||
|
pub fn uniform(z: Vec<f64>, row: &[[f64; 2]], vel: &[[f64; 2]]) -> Self {
|
||||||
|
let n = z.len();
|
||||||
|
Self {
|
||||||
|
z,
|
||||||
|
ns: row.len(),
|
||||||
|
xy: row.iter().copied().cycle().take(n * row.len()).collect(),
|
||||||
|
vel: vel.iter().copied().cycle().take(n * vel.len()).collect(),
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Stations × points.
|
||||||
|
#[must_use]
|
||||||
|
pub fn stations(&self) -> usize {
|
||||||
|
self.z.len()
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The station bracket of `z`: `(k, Some(w))` interpolates stations `k`
|
||||||
|
/// and `k + 1` with weight `w` — beyond the end stations the end
|
||||||
|
/// interval extrapolated by at most half an interval (`w ∈ [−½, 3/2]`:
|
||||||
|
/// a face point of a twisted plate lies a little beyond its station's
|
||||||
|
/// z; the span cut takes over further out); `(0, None)` for a single
|
||||||
|
/// station.
|
||||||
|
#[must_use]
|
||||||
|
pub fn bracket(&self, z: f64) -> (usize, Option<f64>) {
|
||||||
|
let n = self.z.len();
|
||||||
|
if n == 1 {
|
||||||
|
return (0, None);
|
||||||
|
}
|
||||||
|
let k = if z <= self.z[0] {
|
||||||
|
0
|
||||||
|
} else if z >= self.z[n - 1] {
|
||||||
|
n - 2
|
||||||
|
} else {
|
||||||
|
// The last station at or below z.
|
||||||
|
(self.z.partition_point(|&zk| zk <= z) - 1).min(n - 2)
|
||||||
|
};
|
||||||
|
let w = (z - self.z[k]) / (self.z[k + 1] - self.z[k]);
|
||||||
|
(k, Some(w.clamp(-0.5, 1.5)))
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Check the layout (the caller's contract); `Err` names the defect.
|
||||||
|
pub fn validate(&self) -> Result<(), String> {
|
||||||
|
if self.z.is_empty() || self.ns < 2 {
|
||||||
|
return Err("a plate needs at least one station of at least two points".into());
|
||||||
|
}
|
||||||
|
if self.xy.len() != self.z.len() * self.ns {
|
||||||
|
return Err(format!(
|
||||||
|
"plate: {} points for {} stations × {}",
|
||||||
|
self.xy.len(),
|
||||||
|
self.z.len(),
|
||||||
|
self.ns
|
||||||
|
));
|
||||||
|
}
|
||||||
|
if !self.vel.is_empty() && self.vel.len() != self.xy.len() {
|
||||||
|
return Err("plate: the velocities do not match the points".into());
|
||||||
|
}
|
||||||
|
if self.z.windows(2).any(|w| w[1] <= w[0]) {
|
||||||
|
return Err("plate: the stations must ascend strictly".into());
|
||||||
|
}
|
||||||
|
Ok(())
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// `rs_min` of the kernel: `b < a ? b : a`.
|
||||||
|
#[inline]
|
||||||
|
fn rs_min(a: f64, b: f64) -> f64 {
|
||||||
|
if b < a { b } else { a }
|
||||||
|
}
|
||||||
|
|
||||||
|
#[inline]
|
||||||
|
fn rs_max(a: f64, b: f64) -> f64 {
|
||||||
|
if b > a { b } else { a }
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The span cut with rounded edges (the kernel's `span_cut`).
|
||||||
|
#[inline]
|
||||||
|
fn span_cut(d2: f64, z: f64, s: &DeviceSdf) -> f64 {
|
||||||
|
let r = s.r_edge;
|
||||||
|
let q1 = d2 + r;
|
||||||
|
let q2 = (z - s.zc).abs() - 0.5 * s.span + r;
|
||||||
|
let m1 = rs_max(q1, 0.0);
|
||||||
|
let m2 = rs_max(q2, 0.0);
|
||||||
|
let outside = (m1 * m1 + m2 * m2).sqrt();
|
||||||
|
outside + rs_min(rs_max(q1, q2), 0.0) - r
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The closest segment of a polyline given by `pt(m)`, `m < n`: the
|
||||||
|
/// in-plane distance, the segment, its parameter and the offset `(qx, qy)`.
|
||||||
|
#[inline]
|
||||||
|
fn closest<P: Fn(usize) -> [f64; 2]>(
|
||||||
|
x: f64,
|
||||||
|
y: f64,
|
||||||
|
n: usize,
|
||||||
|
pt: P,
|
||||||
|
) -> (f64, usize, f64, [f64; 2]) {
|
||||||
|
let mut best = f64::INFINITY;
|
||||||
|
let (mut mb, mut ub, mut qb) = (0usize, 0.0, [0.0; 2]);
|
||||||
|
for m in 0..n.saturating_sub(1) {
|
||||||
|
let [ax, ay] = pt(m);
|
||||||
|
let [bx, by] = pt(m + 1);
|
||||||
|
let ex = bx - ax;
|
||||||
|
let ey = by - ay;
|
||||||
|
let l2 = ex * ex + ey * ey;
|
||||||
|
let mut u = ((x - ax) * ex + (y - ay) * ey) / l2;
|
||||||
|
if u < 0.0 {
|
||||||
|
u = 0.0;
|
||||||
|
}
|
||||||
|
if u > 1.0 {
|
||||||
|
u = 1.0;
|
||||||
|
}
|
||||||
|
let px = ax + u * ex;
|
||||||
|
let py = ay + u * ey;
|
||||||
|
let qx = x - px;
|
||||||
|
let qy = y - py;
|
||||||
|
let d = (qx * qx + qy * qy).sqrt();
|
||||||
|
if d < best {
|
||||||
|
best = d;
|
||||||
|
mb = m;
|
||||||
|
ub = u;
|
||||||
|
qb = [qx, qy];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
(best, mb, ub, qb)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The plate's in-plane closest point at `(x, y, z)`, the slope-corrected
|
||||||
|
/// distance to the mid-surface, and the interpolated velocity there
|
||||||
|
/// (zero without velocities).
|
||||||
|
fn plate_closest(p: &PlateSurface, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) {
|
||||||
|
let ns = p.ns;
|
||||||
|
let (k, w) = p.bracket(z);
|
||||||
|
let row = |k: usize, m: usize| p.xy[k * ns + m];
|
||||||
|
let lerp2 =
|
||||||
|
|a: [f64; 2], b: [f64; 2], w: f64| [a[0] + w * (b[0] - a[0]), a[1] + w * (b[1] - a[1])];
|
||||||
|
let (best, m, u, q) = match w {
|
||||||
|
None => closest(x, y, ns, |m| row(k, m)),
|
||||||
|
Some(w) => closest(x, y, ns, |m| lerp2(row(k, m), row(k + 1, m), w)),
|
||||||
|
};
|
||||||
|
let mut d = best;
|
||||||
|
if let Some(_w) = w {
|
||||||
|
if d > 0.0 {
|
||||||
|
let dz = p.z[k + 1] - p.z[k];
|
||||||
|
let ca = [
|
||||||
|
row(k + 1, m)[0] - row(k, m)[0],
|
||||||
|
row(k + 1, m)[1] - row(k, m)[1],
|
||||||
|
];
|
||||||
|
let cb = [
|
||||||
|
row(k + 1, m + 1)[0] - row(k, m + 1)[0],
|
||||||
|
row(k + 1, m + 1)[1] - row(k, m + 1)[1],
|
||||||
|
];
|
||||||
|
let cx = (ca[0] + u * (cb[0] - ca[0])) / dz;
|
||||||
|
let cy = (ca[1] + u * (cb[1] - ca[1])) / dz;
|
||||||
|
let qn = (q[0] * cx + q[1] * cy) / d;
|
||||||
|
d /= (1.0 + qn * qn).sqrt();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let v = if p.vel.is_empty() {
|
||||||
|
[0.0, 0.0]
|
||||||
|
} else {
|
||||||
|
let vrow = |k: usize, m: usize| p.vel[k * ns + m];
|
||||||
|
let at = |m: usize| match w {
|
||||||
|
None => vrow(k, m),
|
||||||
|
Some(w) => lerp2(vrow(k, m), vrow(k + 1, m), w),
|
||||||
|
};
|
||||||
|
let (a, b) = (at(m), at(m + 1));
|
||||||
|
[a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])]
|
||||||
|
};
|
||||||
|
(d, v)
|
||||||
|
}
|
||||||
|
|
||||||
|
impl DeviceSdf {
|
||||||
|
/// The circle's distance (cut to the span when `cyl_cut`).
|
||||||
|
#[must_use]
|
||||||
|
pub fn circle_distance_host(&self, x: f64, y: f64, z: f64) -> f64 {
|
||||||
|
let ex0 = x - self.cyl[0];
|
||||||
|
let ey0 = y - self.cyl[1];
|
||||||
|
let mut dc = (ex0 * ex0 + ey0 * ey0).sqrt() - self.cyl[2];
|
||||||
|
if self.cyl_cut {
|
||||||
|
dc = span_cut(dc, z, self);
|
||||||
|
}
|
||||||
|
dc
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The flag's distance (the capsule around the polyline, or the plate)
|
||||||
|
/// and its surface velocity at the closest point.
|
||||||
|
#[must_use]
|
||||||
|
pub fn flag_distance_host(&self, x: f64, y: f64, z: f64) -> (f64, [f64; 2]) {
|
||||||
|
let (best, v) = match self.plate.as_ref() {
|
||||||
|
Some(p) => plate_closest(p, x, y, z),
|
||||||
|
None => {
|
||||||
|
let (best, m, u, _) = closest(x, y, self.poly.len(), |m| self.poly[m]);
|
||||||
|
let v = if self.vel.len() == self.poly.len() && !self.vel.is_empty() {
|
||||||
|
let (a, b) = (self.vel[m], self.vel[m + 1]);
|
||||||
|
[a[0] + u * (b[0] - a[0]), a[1] + u * (b[1] - a[1])]
|
||||||
|
} else {
|
||||||
|
[0.0, 0.0]
|
||||||
|
};
|
||||||
|
(best, v)
|
||||||
|
}
|
||||||
|
};
|
||||||
|
let mut df = best - self.half;
|
||||||
|
if self.flag_cut {
|
||||||
|
df = span_cut(df, z, self);
|
||||||
|
}
|
||||||
|
(df, v)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// φ on the host, the kernel's arithmetic (`geom_phi_at`).
|
||||||
|
#[must_use]
|
||||||
|
pub fn phi_host(&self, x: f64, y: f64, z: f64) -> f64 {
|
||||||
|
let dc = self.circle_distance_host(x, y, z);
|
||||||
|
let (df, _) = self.flag_distance_host(x, y, z);
|
||||||
|
let r = self.fillet;
|
||||||
|
if r > 0.0 && dc < r && df < r {
|
||||||
|
let a = r - dc;
|
||||||
|
let b = r - df;
|
||||||
|
return r - (a * a + b * b).sqrt();
|
||||||
|
}
|
||||||
|
rs_min(dc, df)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Whether `(x, y, z)` belongs to the flag rather than the circle (the
|
||||||
|
/// surface velocity's rule: the flag where it is not farther).
|
||||||
|
#[must_use]
|
||||||
|
pub fn is_flag_host(&self, x: f64, y: f64, z: f64) -> bool {
|
||||||
|
self.flag_distance_host(x, y, z).0 <= self.circle_distance_host(x, y, z)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The surface velocity on the host (the kernel's `body_velocity`):
|
||||||
|
/// the flag's at the closest point where the flag is not farther than
|
||||||
|
/// the circle, zero on the circle, no z component.
|
||||||
|
#[must_use]
|
||||||
|
pub fn velocity_host(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
|
||||||
|
let (df, v) = self.flag_distance_host(x, y, z);
|
||||||
|
if df <= self.circle_distance_host(x, y, z) {
|
||||||
|
(v[0], v[1], 0.0)
|
||||||
|
} else {
|
||||||
|
(0.0, 0.0, 0.0)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The device buffers: the points `(x, y)` interleaved — for a plate
|
||||||
|
/// the stations' rows followed by the stations' z — and the
|
||||||
|
/// velocities interleaved (empty when the body has none).
|
||||||
|
#[must_use]
|
||||||
|
pub fn device_buffers(&self) -> (Vec<f64>, Vec<f64>) {
|
||||||
|
match self.plate.as_ref() {
|
||||||
|
Some(p) => {
|
||||||
|
let mut pts: Vec<f64> = p.xy.iter().flat_map(|q| [q[0], q[1]]).collect();
|
||||||
|
pts.extend_from_slice(&p.z);
|
||||||
|
(pts, p.vel.iter().flat_map(|q| [q[0], q[1]]).collect())
|
||||||
|
}
|
||||||
|
None => (
|
||||||
|
self.poly.iter().flat_map(|q| [q[0], q[1]]).collect(),
|
||||||
|
self.vel.iter().flat_map(|q| [q[0], q[1]]).collect(),
|
||||||
|
),
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Whether the body carries a device velocity matching its points.
|
||||||
|
#[must_use]
|
||||||
|
pub fn has_device_velocity(&self) -> bool {
|
||||||
|
match self.plate.as_ref() {
|
||||||
|
Some(p) => !p.vel.is_empty() && p.vel.len() == p.xy.len(),
|
||||||
|
None => !self.vel.is_empty() && self.vel.len() == self.poly.len(),
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
#[cfg(test)]
|
||||||
|
mod tests {
|
||||||
|
use super::*;
|
||||||
|
|
||||||
|
fn polyline(t: f64) -> (Vec<[f64; 2]>, Vec<[f64; 2]>) {
|
||||||
|
let n = 41;
|
||||||
|
let pts = (0..n)
|
||||||
|
.map(|m| {
|
||||||
|
let s = m as f64 / (n - 1) as f64;
|
||||||
|
[0.25 + 0.35 * s, 0.2 + 0.08 * s * s * (3.0 * t).sin()]
|
||||||
|
})
|
||||||
|
.collect();
|
||||||
|
let vel = (0..n)
|
||||||
|
.map(|m| {
|
||||||
|
let s = m as f64 / (n - 1) as f64;
|
||||||
|
[0.0, 0.24 * s * s * (3.0 * t).cos()]
|
||||||
|
})
|
||||||
|
.collect();
|
||||||
|
(pts, vel)
|
||||||
|
}
|
||||||
|
|
||||||
|
fn sdf(poly: Vec<[f64; 2]>, vel: Vec<[f64; 2]>, plate: Option<PlateSurface>) -> DeviceSdf {
|
||||||
|
DeviceSdf {
|
||||||
|
cyl: [0.2, 0.2, 0.05],
|
||||||
|
cyl_cut: false,
|
||||||
|
flag_cut: true,
|
||||||
|
zc: 0.205,
|
||||||
|
span: 0.2,
|
||||||
|
r_edge: 0.0066,
|
||||||
|
half: 0.01,
|
||||||
|
fillet: 0.0,
|
||||||
|
poly,
|
||||||
|
vel,
|
||||||
|
plate,
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/// G1 (host): a span-uniform plate is the polyline capsule to the bit —
|
||||||
|
/// φ, the part and the surface velocity — on a lattice through the
|
||||||
|
/// flag, its span edges, the stations and beyond them.
|
||||||
|
#[test]
|
||||||
|
fn uniform_plate_is_the_polyline_capsule_bit_for_bit() {
|
||||||
|
let (poly, vel) = polyline(0.37);
|
||||||
|
let z: Vec<f64> = (0..11).map(|k| 0.105 + 0.02 * k as f64).collect();
|
||||||
|
let plate = PlateSurface::uniform(z, &poly, &vel);
|
||||||
|
plate.validate().unwrap();
|
||||||
|
let a = sdf(poly.clone(), vel.clone(), None);
|
||||||
|
let b = sdf(Vec::new(), Vec::new(), Some(plate));
|
||||||
|
let mut n = 0;
|
||||||
|
for i in 0..90 {
|
||||||
|
for j in 0..40 {
|
||||||
|
for k in 0..45 {
|
||||||
|
let (x, y, zz) = (
|
||||||
|
0.18 + 0.0051 * i as f64,
|
||||||
|
0.13 + 0.0037 * j as f64,
|
||||||
|
0.07 + 0.0063 * k as f64,
|
||||||
|
);
|
||||||
|
assert_eq!(
|
||||||
|
a.phi_host(x, y, zz).to_bits(),
|
||||||
|
b.phi_host(x, y, zz).to_bits()
|
||||||
|
);
|
||||||
|
let (va, vb) = (a.velocity_host(x, y, zz), b.velocity_host(x, y, zz));
|
||||||
|
assert_eq!(va.0.to_bits(), vb.0.to_bits());
|
||||||
|
assert_eq!(va.1.to_bits(), vb.1.to_bits());
|
||||||
|
n += 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
assert_eq!(n, 90 * 40 * 45);
|
||||||
|
}
|
||||||
|
|
||||||
|
/// A plate tilted rigidly along the span (y = y0 + a (z − zc)) has
|
||||||
|
/// the exact distance `|y − y(z)| / sqrt(1 + a²)` from its mid-plane:
|
||||||
|
/// the slope correction recovers it where the in-plane distance is
|
||||||
|
/// `sqrt(1 + a²)` too large.
|
||||||
|
#[test]
|
||||||
|
fn spanwise_slope_correction_is_exact_on_a_tilted_plane() {
|
||||||
|
let a = 0.4;
|
||||||
|
let z: Vec<f64> = vec![0.0, 0.1, 0.2];
|
||||||
|
let ns = 3;
|
||||||
|
let mut xy = Vec::new();
|
||||||
|
for &zk in &z {
|
||||||
|
for m in 0..ns {
|
||||||
|
xy.push([0.1 * m as f64, 0.5 + a * (zk - 0.1)]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let p = PlateSurface {
|
||||||
|
z,
|
||||||
|
ns,
|
||||||
|
xy,
|
||||||
|
vel: Vec::new(),
|
||||||
|
};
|
||||||
|
for &(y, zz) in &[(0.53, 0.05), (0.47, 0.15), (0.6, 0.12)] {
|
||||||
|
let (d, _) = plate_closest(&p, 0.1, y, zz);
|
||||||
|
let exact = (y - (0.5 + a * (zz - 0.1))).abs() / (1.0 + a * a).sqrt();
|
||||||
|
assert!((d - exact).abs() < 1e-15, "{d} vs {exact}");
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
+14
-12
@@ -73,6 +73,8 @@ struct GeomSdf {
|
|||||||
cyl_cut: i32,
|
cyl_cut: i32,
|
||||||
flag_cut: i32,
|
flag_cut: i32,
|
||||||
npts: i32,
|
npts: i32,
|
||||||
|
nst: i32,
|
||||||
|
ns: i32,
|
||||||
}
|
}
|
||||||
unsafe impl DeviceRepr for GeomSdf {}
|
unsafe impl DeviceRepr for GeomSdf {}
|
||||||
unsafe impl ValidAsZeroBits for GeomSdf {}
|
unsafe impl ValidAsZeroBits for GeomSdf {}
|
||||||
@@ -147,19 +149,18 @@ fn geom_sdf(sdf: &DeviceSdf) -> GeomSdf {
|
|||||||
cyl_cut: i32::from(sdf.cyl_cut),
|
cyl_cut: i32::from(sdf.cyl_cut),
|
||||||
flag_cut: i32::from(sdf.flag_cut),
|
flag_cut: i32::from(sdf.flag_cut),
|
||||||
npts: sdf.poly.len() as i32,
|
npts: sdf.poly.len() as i32,
|
||||||
|
nst: sdf.plate.as_ref().map_or(0, |p| p.z.len() as i32),
|
||||||
|
ns: sdf.plate.as_ref().map_or(0, |p| p.ns as i32),
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Interleaved (x, y) pairs on the device (one dummy entry when empty).
|
/// A flat f64 buffer on the device (one dummy entry when empty): the
|
||||||
fn upload_pairs(v: &[[f64; 2]]) -> CudaSlice<f64> {
|
/// body's points / velocities as `DeviceSdf::device_buffers` lays them out
|
||||||
let flat: Vec<f64> = v.iter().flat_map(|p| [p[0], p[1]]).collect();
|
/// (the polyline's `(x, y)` pairs interleaved, as before R8-c).
|
||||||
|
fn upload_flat(flat: &[f64]) -> CudaSlice<f64> {
|
||||||
runtime()
|
runtime()
|
||||||
.stream
|
.stream
|
||||||
.memcpy_stod(if flat.is_empty() {
|
.memcpy_stod(if flat.is_empty() { &[0.0f64][..] } else { flat })
|
||||||
&[0.0f64][..]
|
|
||||||
} else {
|
|
||||||
&flat
|
|
||||||
})
|
|
||||||
.expect("upload pairs")
|
.expect("upload pairs")
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -188,13 +189,14 @@ impl UbBody {
|
|||||||
return None;
|
return None;
|
||||||
}
|
}
|
||||||
let sdf = solver.body()?.device_sdf(t)?;
|
let sdf = solver.body()?.device_sdf(t)?;
|
||||||
if sdf.vel.is_empty() || sdf.vel.len() != sdf.poly.len() {
|
if !sdf.has_device_velocity() {
|
||||||
return None;
|
return None;
|
||||||
}
|
}
|
||||||
|
let (pts, vel) = sdf.device_buffers();
|
||||||
Some(Self {
|
Some(Self {
|
||||||
gs: geom_sdf(&sdf),
|
gs: geom_sdf(&sdf),
|
||||||
poly: upload_pairs(&sdf.poly),
|
poly: upload_flat(&pts),
|
||||||
vel: upload_pairs(&sdf.vel),
|
vel: upload_flat(&vel),
|
||||||
})
|
})
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@@ -402,7 +404,7 @@ impl DeviceGeom {
|
|||||||
}
|
}
|
||||||
let gg = geom_grid(g);
|
let gg = geom_grid(g);
|
||||||
let gs = geom_sdf(&sdf);
|
let gs = geom_sdf(&sdf);
|
||||||
let d_poly = upload_pairs(&sdf.poly);
|
let d_poly = upload_flat(&sdf.device_buffers().0);
|
||||||
let (has_prev, band, motion) = match prev {
|
let (has_prev, band, motion) = match prev {
|
||||||
Some((_, band, motion)) => (1i32, band, motion),
|
Some((_, band, motion)) => (1i32, band, motion),
|
||||||
None => (0i32, 0.0, 0.0),
|
None => (0i32, 0.0, 0.0),
|
||||||
|
|||||||
@@ -26,8 +26,8 @@ use embedded3_flag_kinematics::{Recorded, recorded};
|
|||||||
use rtx_cfd::solvers::incompressible::ConvectionScheme;
|
use rtx_cfd::solvers::incompressible::ConvectionScheme;
|
||||||
use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep;
|
use rtx_cfd::solvers::incompressible::embedded3::step::device::DeviceStep;
|
||||||
use rtx_cfd::solvers::incompressible::embedded3::{
|
use rtx_cfd::solvers::incompressible::embedded3::{
|
||||||
Body, Boundaries, DeviceSdf, Field, Fluid, Grid, Parameters, Side, Solver, WallScheme,
|
Body, Boundaries, DeviceSdf, Field, Fluid, Grid, HexPlate, Parameters, PlateSurface, Side,
|
||||||
write_vtk,
|
Solver, WallScheme, write_vtk,
|
||||||
};
|
};
|
||||||
use std::io::Write as _;
|
use std::io::Write as _;
|
||||||
|
|
||||||
@@ -287,6 +287,99 @@ fn cylinder_3d(d2: f64, z: f64, r: f64) -> f64 {
|
|||||||
outside + q1.max(q2).min(0.0) - r
|
outside + q1.max(q2).min(0.0) - r
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/// R8-c: the flag as a deformed plate (`RTX_E3_FLAG_BODY=plate`): the
|
||||||
|
/// span stations (`RTX_E3_FLAG_STATIONS`, default 21, spread over the
|
||||||
|
/// flag's span) each carry the centreline polyline. `RTX_E3_FLAG_TWIST=κ`
|
||||||
|
/// (default 0; analytic mode only) scales each station's deflection and
|
||||||
|
/// velocity by `1 + κ ζ`, `ζ = (z − z_c)/(span/2)` — the first bending mode
|
||||||
|
/// times a span-linear twist. At κ = 0 every station is the polyline as it
|
||||||
|
/// is (the G1 identity with the polyline capsule).
|
||||||
|
fn plate_body() -> bool {
|
||||||
|
std::env::var("RTX_E3_FLAG_BODY").is_ok_and(|v| v == "plate")
|
||||||
|
}
|
||||||
|
|
||||||
|
fn twist() -> f64 {
|
||||||
|
env_f("RTX_E3_FLAG_TWIST", 0.0)
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The span stations' z (ascending) and their ζ.
|
||||||
|
fn stations() -> Vec<f64> {
|
||||||
|
let n = (env_f("RTX_E3_FLAG_STATIONS", 21.0) as usize).max(2);
|
||||||
|
let (zc, span) = (0.5 * duct_depth(), flag_span());
|
||||||
|
(0..n)
|
||||||
|
.map(|k| zc - 0.5 * span + span * k as f64 / (n - 1) as f64)
|
||||||
|
.collect()
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The span factor `1 + κ ζ` of the deflection at `z` (clamped to the span).
|
||||||
|
fn span_factor(z: f64) -> f64 {
|
||||||
|
let (zc, span) = (0.5 * duct_depth(), flag_span());
|
||||||
|
let zeta = ((z - zc) / (0.5 * span)).clamp(-1.0, 1.0);
|
||||||
|
1.0 + twist() * zeta
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The plate at `t`: per station the centreline (and its velocity).
|
||||||
|
fn plate_at(t: f64) -> PlateSurface {
|
||||||
|
let z = stations();
|
||||||
|
let (row, vel): (Vec<[f64; 2]>, Vec<[f64; 2]>) = match recorded() {
|
||||||
|
Some(rec) => recorded_polyline(rec, t)
|
||||||
|
.iter()
|
||||||
|
.map(|p| ([p.0, p.1], [p.2, p.3]))
|
||||||
|
.unzip(),
|
||||||
|
None => analytic_polyline(t)
|
||||||
|
.iter()
|
||||||
|
.map(|p| ([p.0, p.1], [0.0, p.2]))
|
||||||
|
.unzip(),
|
||||||
|
};
|
||||||
|
if twist() == 0.0 {
|
||||||
|
return PlateSurface::uniform(z, &row, &vel);
|
||||||
|
}
|
||||||
|
assert!(
|
||||||
|
recorded().is_none(),
|
||||||
|
"RTX_E3_FLAG_TWIST: the analytic mode only"
|
||||||
|
);
|
||||||
|
let ns = N + 1;
|
||||||
|
let (mut xy, mut vv) = (Vec::new(), Vec::new());
|
||||||
|
for &zk in &z {
|
||||||
|
let f = span_factor(zk);
|
||||||
|
let mut pts: Vec<(f64, f64, f64, f64)> = (0..ns)
|
||||||
|
.map(|m| {
|
||||||
|
let s = m as f64 / N as f64;
|
||||||
|
let (d, v) = deflection(s, t);
|
||||||
|
(FLAG_X0 + s * FLAG_LEN, body_cy() + d * f, 0.0, v * f)
|
||||||
|
})
|
||||||
|
.collect();
|
||||||
|
inset_last(&mut pts, tip_inset());
|
||||||
|
xy.extend(pts.iter().map(|p| [p.0, p.1]));
|
||||||
|
vv.extend(pts.iter().map(|p| [p.2, p.3]));
|
||||||
|
}
|
||||||
|
PlateSurface { z, ns, xy, vel: vv }
|
||||||
|
}
|
||||||
|
|
||||||
|
/// The deformed flag's mid-surface point `y = w(x, z)` and its 3D unit
|
||||||
|
/// normal at arc fraction `s` and span `z` (the analytic kinematics with the
|
||||||
|
/// span factor; the structure's placement for the load transfer: the
|
||||||
|
/// thickness along the mid-surface's normal, as a solid plate carries it).
|
||||||
|
fn mid_point(s: f64, z: f64, t: f64) -> ([f64; 3], [f64; 3]) {
|
||||||
|
let f = span_factor(z);
|
||||||
|
let (d, _) = deflection(s, t);
|
||||||
|
let ds = 1e-6;
|
||||||
|
let (s0, s1) = ((s - ds).max(0.0), (s + ds).min(1.0));
|
||||||
|
let wx = (deflection(s1, t).0 - deflection(s0, t).0) / (s1 - s0) / FLAG_LEN * f;
|
||||||
|
// The span factor's rate: κ / (span/2) inside the span.
|
||||||
|
let (zc, span) = (0.5 * duct_depth(), flag_span());
|
||||||
|
let wz = if ((z - zc) / (0.5 * span)).abs() < 1.0 {
|
||||||
|
d * twist() / (0.5 * span)
|
||||||
|
} else {
|
||||||
|
0.0
|
||||||
|
};
|
||||||
|
let r = (1.0 + wx * wx + wz * wz).sqrt();
|
||||||
|
(
|
||||||
|
[FLAG_X0 + s * FLAG_LEN, body_cy() + d * f, z],
|
||||||
|
[-wx / r, 1.0 / r, -wz / r],
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
fn inflow(y: f64, z: f64) -> f64 {
|
fn inflow(y: f64, z: f64) -> f64 {
|
||||||
let (hd, d) = (duct_height(), duct_depth());
|
let (hd, d) = (duct_height(), duct_depth());
|
||||||
16.0 * U_M * y * z * (hd - y) * (d - z) / (hd * hd * d * d)
|
16.0 * U_M * y * z * (hd - y) * (d - z) / (hd * hd * d * d)
|
||||||
@@ -385,20 +478,10 @@ fn flag_wake_on_the_device() {
|
|||||||
let cy = body_cy();
|
let cy = body_cy();
|
||||||
let cyl = move |x: f64, y: f64| ((x - CX).powi(2) + (y - cy).powi(2)).sqrt() - R_CYL;
|
let cyl = move |x: f64, y: f64| ((x - CX).powi(2) + (y - cy).powi(2)).sqrt() - R_CYL;
|
||||||
let r_fillet = root_fillet();
|
let r_fillet = root_fillet();
|
||||||
let body = Body::from_sdf(move |x, y, z, t| {
|
// The device form of φ: the polyline capsule (R6-1), or the plate (R8-c).
|
||||||
fillet_union(cylinder_3d(cyl(x, y), z, r_edge), flag_3d(x, y, z, t, r_edge).0, r_fillet)
|
let device_sdf = move |t: f64| {
|
||||||
})
|
let plate = plate_body().then(|| plate_at(t));
|
||||||
.with_surface_velocity(move |x, y, z, t| {
|
DeviceSdf {
|
||||||
let (df, (vx, vy)) = flag_3d(x, y, z, t, r_edge);
|
|
||||||
if df <= cylinder_3d(cyl(x, y), z, r_edge) {
|
|
||||||
(vx, vy, 0.0)
|
|
||||||
} else {
|
|
||||||
(0.0, 0.0, 0.0)
|
|
||||||
}
|
|
||||||
});
|
|
||||||
// R6-1: the same φ in the device's form (the device geometry, default ON): the
|
|
||||||
// circle, the capsule around the step's centreline, the span cuts.
|
|
||||||
let body = body.with_device_sdf(move |t| DeviceSdf {
|
|
||||||
cyl: [CX, cy, R_CYL],
|
cyl: [CX, cy, R_CYL],
|
||||||
cyl_cut: !(flag_span() >= duct_depth()
|
cyl_cut: !(flag_span() >= duct_depth()
|
||||||
|| !std::env::var("RTX_E3_FLAG_CYL_SPAN").is_ok_and(|v| v == "flag")),
|
|| !std::env::var("RTX_E3_FLAG_CYL_SPAN").is_ok_and(|v| v == "flag")),
|
||||||
@@ -408,22 +491,71 @@ fn flag_wake_on_the_device() {
|
|||||||
r_edge,
|
r_edge,
|
||||||
half: FLAG_HALF,
|
half: FLAG_HALF,
|
||||||
fillet: r_fillet,
|
fillet: r_fillet,
|
||||||
poly: match recorded() {
|
poly: match (&plate, recorded()) {
|
||||||
Some(rec) => recorded_polyline(rec, t)
|
(Some(_), _) => Vec::new(),
|
||||||
|
(None, Some(rec)) => recorded_polyline(rec, t)
|
||||||
.iter()
|
.iter()
|
||||||
.map(|p| [p.0, p.1])
|
.map(|p| [p.0, p.1])
|
||||||
.collect(),
|
.collect(),
|
||||||
None => analytic_polyline(t).iter().map(|p| [p.0, p.1]).collect(),
|
(None, None) => analytic_polyline(t).iter().map(|p| [p.0, p.1]).collect(),
|
||||||
},
|
},
|
||||||
// R6-2 step 2: the centreline's velocity per point (the analytic mode is transverse).
|
// R6-2 step 2: the centreline's velocity per point (the analytic mode is transverse).
|
||||||
vel: match recorded() {
|
vel: match (&plate, recorded()) {
|
||||||
Some(rec) => recorded_polyline(rec, t)
|
(Some(_), _) => Vec::new(),
|
||||||
|
(None, Some(rec)) => recorded_polyline(rec, t)
|
||||||
.iter()
|
.iter()
|
||||||
.map(|p| [p.2, p.3])
|
.map(|p| [p.2, p.3])
|
||||||
.collect(),
|
.collect(),
|
||||||
None => analytic_polyline(t).iter().map(|p| [0.0, p.2]).collect(),
|
(None, None) => analytic_polyline(t).iter().map(|p| [0.0, p.2]).collect(),
|
||||||
},
|
},
|
||||||
});
|
plate,
|
||||||
|
}
|
||||||
|
};
|
||||||
|
// The device form at `t`, once per thread and time (the host closures
|
||||||
|
// of the plate body evaluate it ~10⁶ times per step).
|
||||||
|
let sdf_at = move |t: f64| -> std::sync::Arc<DeviceSdf> {
|
||||||
|
thread_local! {
|
||||||
|
static SDF: std::cell::RefCell<(u64, Option<std::sync::Arc<DeviceSdf>>)> =
|
||||||
|
const { std::cell::RefCell::new((u64::MAX, None)) };
|
||||||
|
}
|
||||||
|
SDF.with(|cell| {
|
||||||
|
let mut c = cell.borrow_mut();
|
||||||
|
if c.0 != t.to_bits() || c.1.is_none() {
|
||||||
|
c.1 = Some(std::sync::Arc::new(device_sdf(t)));
|
||||||
|
c.0 = t.to_bits();
|
||||||
|
}
|
||||||
|
c.1.clone().expect("sdf")
|
||||||
|
})
|
||||||
|
};
|
||||||
|
let body = if plate_body() {
|
||||||
|
// R8-c: the host φ and surface velocity ARE the device form's (the
|
||||||
|
// kernel's arithmetic on the host).
|
||||||
|
Body::from_sdf(move |x, y, z, t| sdf_at(t).phi_host(x, y, z))
|
||||||
|
.with_surface_velocity(move |x, y, z, t| sdf_at(t).velocity_host(x, y, z))
|
||||||
|
} else {
|
||||||
|
Body::from_sdf(move |x, y, z, t| {
|
||||||
|
fillet_union(
|
||||||
|
cylinder_3d(cyl(x, y), z, r_edge),
|
||||||
|
flag_3d(x, y, z, t, r_edge).0,
|
||||||
|
r_fillet,
|
||||||
|
)
|
||||||
|
})
|
||||||
|
.with_surface_velocity(move |x, y, z, t| {
|
||||||
|
let (df, (vx, vy)) = flag_3d(x, y, z, t, r_edge);
|
||||||
|
if df <= cylinder_3d(cyl(x, y), z, r_edge) {
|
||||||
|
(vx, vy, 0.0)
|
||||||
|
} else {
|
||||||
|
(0.0, 0.0, 0.0)
|
||||||
|
}
|
||||||
|
})
|
||||||
|
};
|
||||||
|
assert!(
|
||||||
|
twist() == 0.0 || plate_body(),
|
||||||
|
"RTX_E3_FLAG_TWIST needs RTX_E3_FLAG_BODY=plate"
|
||||||
|
);
|
||||||
|
// R6-1: the same φ in the device's form (the device geometry, default ON): the
|
||||||
|
// circle, the capsule around the step's centreline (or the plate), the span cuts.
|
||||||
|
let body = body.with_device_sdf(device_sdf);
|
||||||
solver.set_moving_body(body);
|
solver.set_moving_body(body);
|
||||||
let g = Grid::cubic(nx, ny_grid, nz, h);
|
let g = Grid::cubic(nx, ny_grid, nz, h);
|
||||||
let mut field = Field::new(g);
|
let mut field = Field::new(g);
|
||||||
@@ -484,6 +616,32 @@ fn flag_wake_on_the_device() {
|
|||||||
(mid - 2, mid + 2)
|
(mid - 2, mid + 2)
|
||||||
};
|
};
|
||||||
let width = nz as f64 * h;
|
let width = nz as f64 * h;
|
||||||
|
// R8-c: the load transfer onto the structure's Hex20 plate (35 × 2 × n
|
||||||
|
// with n = `RTX_E3_FLAG_TRANSFER`; off when unset) every
|
||||||
|
// `RTX_E3_FLAG_TRANSFER_EVERY`-th sample (default 1), the budget to
|
||||||
|
// `RTX_E3_FLAG_TRANSFER_CSV`, the last transfer's nodal forces to
|
||||||
|
// `RTX_E3_FLAG_TRANSFER_NODAL`. The plate is placed on the prescribed
|
||||||
|
// kinematics (analytic mode): the centreline with the span factor, the
|
||||||
|
// thickness along its in-plane normal, the flag's span.
|
||||||
|
let transfer_nz = env_f("RTX_E3_FLAG_TRANSFER", 0.0) as usize;
|
||||||
|
let transfer_every = (env_f("RTX_E3_FLAG_TRANSFER_EVERY", 1.0) as usize).max(1);
|
||||||
|
let hex = (transfer_nz > 0).then(|| {
|
||||||
|
assert!(
|
||||||
|
recorded().is_none(),
|
||||||
|
"RTX_E3_FLAG_TRANSFER: the analytic mode only"
|
||||||
|
);
|
||||||
|
HexPlate::new(35, 2, transfer_nz)
|
||||||
|
});
|
||||||
|
let mut transfer_csv = std::env::var("RTX_E3_FLAG_TRANSFER_CSV").ok().map(|p| {
|
||||||
|
let mut f = std::fs::File::create(p).expect("transfer csv");
|
||||||
|
writeln!(
|
||||||
|
f,
|
||||||
|
"t,loads,flag_loads,route_x,route_y,route_z,sum_dx,sum_dy,sum_dz,fin_x,fin_y,fin_z,dfx,dfy,dfz,min_x,min_y,min_z,dmx,dmy,dmz,rel_force,rel_moment,max_newton,max_outside,extrapolated,extrap_share,interior_share,lever_x,lever_y,lever_z,ms"
|
||||||
|
)
|
||||||
|
.unwrap();
|
||||||
|
f
|
||||||
|
});
|
||||||
|
let mut samples_seen = 0usize;
|
||||||
let start = std::time::Instant::now();
|
let start = std::time::Instant::now();
|
||||||
let mut drag_rec_sum = 0.0;
|
let mut drag_rec_sum = 0.0;
|
||||||
// The routes' PARTS over the whole body (x, per unit width): operator
|
// The routes' PARTS over the whole body (x, per unit width): operator
|
||||||
@@ -520,6 +678,128 @@ fn flag_wake_on_the_device() {
|
|||||||
let ft = mask
|
let ft = mask
|
||||||
.cut_wall_force(body, &field, RHO * NU, t)
|
.cut_wall_force(body, &field, RHO * NU, t)
|
||||||
.expect("wall");
|
.expect("wall");
|
||||||
|
if sample {
|
||||||
|
samples_seen += 1;
|
||||||
|
}
|
||||||
|
if let (Some(hex), true) = (hex.as_ref(), sample && samples_seen % transfer_every == 0)
|
||||||
|
{
|
||||||
|
let lap = std::time::Instant::now();
|
||||||
|
let (loads, route) = mask
|
||||||
|
.cut_wall_loads(body, &field, RHO * NU, t)
|
||||||
|
.expect("loads");
|
||||||
|
// The route's total is cut_wall_force's, to the bit (same loops).
|
||||||
|
assert_eq!(route.map(f64::to_bits), ft.map(f64::to_bits), "route total");
|
||||||
|
let mut sum = [0.0f64; 3];
|
||||||
|
for l in &loads {
|
||||||
|
for c in 0..3 {
|
||||||
|
sum[c] += l.f[c];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let sdf = body.device_sdf(t).expect("device form");
|
||||||
|
let flag: Vec<_> = loads
|
||||||
|
.iter()
|
||||||
|
.filter(|l| sdf.is_flag_host(l.foot[0], l.foot[1], l.foot[2]))
|
||||||
|
.collect();
|
||||||
|
let (zc, span) = (0.5 * duct_depth(), flag_span());
|
||||||
|
let pos = hex.place(|s, eta, zeta| {
|
||||||
|
let z = zc - 0.5 * span + span * zeta;
|
||||||
|
let (c, n) = mid_point(s, z, t);
|
||||||
|
[
|
||||||
|
c[0] + FLAG_HALF * eta * n[0],
|
||||||
|
c[1] + FLAG_HALF * eta * n[1],
|
||||||
|
c[2] + FLAG_HALF * eta * n[2],
|
||||||
|
]
|
||||||
|
});
|
||||||
|
let pairs: Vec<([f64; 3], [f64; 3])> = flag.iter().map(|l| (l.foot, l.f)).collect();
|
||||||
|
let origin = [FLAG_X0, body_cy(), zc];
|
||||||
|
let tr = hex.transfer(&pos, &pairs, origin);
|
||||||
|
// The lever the foot adds over the operator point: Σ (foot − x) × F.
|
||||||
|
let mut lever = [0.0f64; 3];
|
||||||
|
for l in &flag {
|
||||||
|
let d = [l.foot[0] - l.x[0], l.foot[1] - l.x[1], l.foot[2] - l.x[2]];
|
||||||
|
let m = [
|
||||||
|
d[1] * l.f[2] - d[2] * l.f[1],
|
||||||
|
d[2] * l.f[0] - d[0] * l.f[2],
|
||||||
|
d[0] * l.f[1] - d[1] * l.f[0],
|
||||||
|
];
|
||||||
|
for c in 0..3 {
|
||||||
|
lever[c] += m[c];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
let nrm = |a: [f64; 3]| (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt();
|
||||||
|
let df = [0, 1, 2].map(|c| tr.force_out[c] - tr.force_in[c]);
|
||||||
|
let dm = [0, 1, 2].map(|c| tr.moment_out[c] - tr.moment_in[c]);
|
||||||
|
let abs_f: f64 = flag.iter().map(|l| nrm(l.f)).sum();
|
||||||
|
let rel_f = nrm(df) / abs_f.max(1e-300);
|
||||||
|
let rel_m = nrm(dm) / (abs_f * FLAG_LEN).max(1e-300);
|
||||||
|
let extrap_share = tr.extrapolated_load / abs_f.max(1e-300);
|
||||||
|
let ms = lap.elapsed().as_secs_f64() * 1e3;
|
||||||
|
println!(
|
||||||
|
" transfer t {t:.4}: {} loads ({} flag); route − Σ {:.1e} N; flag F {:+.4} {:+.4} {:+.4} N, M {:+.5} {:+.5} {:+.5} N m; |ΔF|/Σ|F| {rel_f:.1e}, |ΔM|/(Σ|F| L) {rel_m:.1e}; Newton ≤ {:.1e} m, outside ≤ {:.2e} ({} extrapolated, {:.1e} of Σ|F|), interior share {:.3}; {ms:.0} ms",
|
||||||
|
loads.len(),
|
||||||
|
flag.len(),
|
||||||
|
nrm([route[0] - sum[0], route[1] - sum[1], route[2] - sum[2]]),
|
||||||
|
tr.force_in[0],
|
||||||
|
tr.force_in[1],
|
||||||
|
tr.force_in[2],
|
||||||
|
tr.moment_in[0],
|
||||||
|
tr.moment_in[1],
|
||||||
|
tr.moment_in[2],
|
||||||
|
tr.max_residual,
|
||||||
|
tr.max_outside,
|
||||||
|
tr.extrapolated,
|
||||||
|
extrap_share,
|
||||||
|
tr.interior_share
|
||||||
|
);
|
||||||
|
if let Some(f) = transfer_csv.as_mut() {
|
||||||
|
writeln!(
|
||||||
|
f,
|
||||||
|
"{t:.6},{},{},{:.9e},{:.9e},{:.9e},{:.3e},{:.3e},{:.3e},{:.9e},{:.9e},{:.9e},{:.3e},{:.3e},{:.3e},{:.9e},{:.9e},{:.9e},{:.3e},{:.3e},{:.3e},{rel_f:.3e},{rel_m:.3e},{:.3e},{:.3e},{},{extrap_share:.3e},{:.4e},{:.4e},{:.4e},{:.4e},{ms:.1}",
|
||||||
|
loads.len(),
|
||||||
|
flag.len(),
|
||||||
|
route[0],
|
||||||
|
route[1],
|
||||||
|
route[2],
|
||||||
|
route[0] - sum[0],
|
||||||
|
route[1] - sum[1],
|
||||||
|
route[2] - sum[2],
|
||||||
|
tr.force_in[0],
|
||||||
|
tr.force_in[1],
|
||||||
|
tr.force_in[2],
|
||||||
|
df[0],
|
||||||
|
df[1],
|
||||||
|
df[2],
|
||||||
|
tr.moment_in[0],
|
||||||
|
tr.moment_in[1],
|
||||||
|
tr.moment_in[2],
|
||||||
|
dm[0],
|
||||||
|
dm[1],
|
||||||
|
dm[2],
|
||||||
|
tr.max_residual,
|
||||||
|
tr.max_outside,
|
||||||
|
tr.extrapolated,
|
||||||
|
tr.interior_share,
|
||||||
|
lever[0],
|
||||||
|
lever[1],
|
||||||
|
lever[2]
|
||||||
|
)
|
||||||
|
.unwrap();
|
||||||
|
}
|
||||||
|
if let Ok(path) = std::env::var("RTX_E3_FLAG_TRANSFER_NODAL") {
|
||||||
|
let mut f = std::fs::File::create(path).expect("nodal csv");
|
||||||
|
writeln!(f, "node,i,j,k,x,y,z,fx,fy,fz").unwrap();
|
||||||
|
for (n, fv) in tr.nodal.iter().enumerate() {
|
||||||
|
let [i, j, k] = hex.lattice_of(n);
|
||||||
|
let x = pos[n];
|
||||||
|
writeln!(
|
||||||
|
f,
|
||||||
|
"{n},{i},{j},{k},{:.9e},{:.9e},{:.9e},{:.9e},{:.9e},{:.9e}",
|
||||||
|
x[0], x[1], x[2], fv[0], fv[1], fv[2]
|
||||||
|
)
|
||||||
|
.unwrap();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
// The tip's transverse deflection: the record's last station in
|
// The tip's transverse deflection: the record's last station in
|
||||||
// recorded mode (until 2026-09-21 this column held the analytic
|
// recorded mode (until 2026-09-21 this column held the analytic
|
||||||
// first mode even then — R2's fits use the record directly).
|
// first mode even then — R2's fits use the record directly).
|
||||||
|
|||||||
Reference in New Issue
Block a user