rustytorch/crates/specialized/rtx-fsi/src/transfer.rs

//! Conservative load and motion transfer across a non-matching interface.
//!
//! # The two constraints that make it conservative
//!
//! Each fluid face's load is distributed onto nearby structure nodes with
//! weights `w_i`. Two conditions decide whether the transfer conserves
//! anything:
//!
//! - `sum(w_i) = 1` — **partition of unity**. Total force is preserved.
//! - `sum(w_i * x_i) = x_face` — **linear reproduction**. The load arrives
//!   where it left, in the weighted-average sense, so total **moment** is
//!   preserved too, about any point.
//!
//! The second is the one that gets skipped. Inverse-distance weighting
//! satisfies partition of unity and generally violates linear
//! reproduction, which conserves force while quietly corrupting moment —
//! a defect that shows up as a slow spurious rotation rather than as an
//! obvious error.
//!
//! Both constraints leave the weights underdetermined for more than four
//! nodes, so this takes the **minimum-norm** solution: `w = Aᵀ(AAᵀ)⁻¹ b`,
//! where `A` stacks the ones row and the coordinate rows. Minimum norm
//! keeps the load spread rather than concentrated on whichever node the
//! solver happened to favour.
//!
//! # Why motion transfer is the transpose
//!
//! Given the load operator `H` mapping face loads to nodal forces, using
//! `Hᵀ` to map nodal velocities to face velocities makes interface work
//! conserved *identically*:
//!
//! ```text
//! (H f)·v  =  f·(Hᵀ v)
//! ```
//!
//! which is just the definition of the transpose. Any other pairing leaks
//! energy across the interface every step, and the leak looks like physics
//! until it destabilises.

use nalgebra::{Matrix4, Vector3, Vector4};

use crate::error::FsiError;

/// Structure nodes recruited per fluid face.
///
/// Four is the minimum for linear reproduction in three dimensions. More
/// spreads the load and conditions the system better; too many turn a
/// local transfer into a global smear.
const NEIGHBOURS: usize = 8;

/// Singular values below this share of the largest are treated as zero
/// when forming the pseudo-inverse.
const SINGULAR_TOLERANCE: f64 = 1e-12;

/// How exactly the partition-of-unity and linear-reproduction constraints
/// must hold before a face is accepted. Tight, because these are the
/// properties the conservation guarantees rest on.
const CONSTRAINT_TOLERANCE: f64 = 1e-9;

/// One fluid-side boundary face on the wetted surface.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FluidFace {
    /// Face centroid.
    pub centroid: Vector3<f64>,
    /// Outward unit normal, pointing into the fluid.
    pub normal: Vector3<f64>,
    /// Face area, used to turn a traction into a force.
    pub area: f64,
}

/// The interface between the two meshes, with the transfer operator
/// already formed.
#[derive(Debug, Clone, PartialEq)]
pub struct WettedSurface {
    /// Per face: the recruited nodes and their weights.
    rows: Vec<Vec<(usize, f64)>>,
    node_count: usize,
}

impl WettedSurface {
    /// Form the transfer operator between a fluid boundary and a
    /// structural interface mesh.
    ///
    /// # Errors
    /// - [`FsiError::EmptyInterface`] if either side is empty.
    /// - [`FsiError::InsufficientNodes`] with fewer than four nodes.
    /// - [`FsiError::DegenerateNeighbourhood`] if a face's nearest nodes
    ///   are collinear or coplanar in a way that makes the constraint
    ///   system singular.
    pub fn build(
        fluid_faces: &[FluidFace],
        structure_nodes: &[Vector3<f64>],
    ) -> Result<Self, FsiError> {
        if fluid_faces.is_empty() {
            return Err(FsiError::EmptyInterface { side: "fluid" });
        }
        if structure_nodes.is_empty() {
            return Err(FsiError::EmptyInterface { side: "structure" });
        }
        if structure_nodes.len() < 4 {
            return Err(FsiError::InsufficientNodes {
                got: structure_nodes.len(),
                needed: 4,
            });
        }

        let mut rows = Vec::with_capacity(fluid_faces.len());
        for (face_index, face) in fluid_faces.iter().enumerate() {
            // Adaptive recruitment: a NEARLY collinear neighbourhood (the
            // nearest nodes of a face on a smoothly DEFORMED edge — y is
            // almost a linear function of x, off by the curvature sagitta)
            // cannot satisfy exact centroid reproduction with bounded
            // weights: the offending singular value is too large to
            // truncate and too small to invert accurately, and any
            // truncated solve misses reproduction by the sagitta. No
            // weight choice fixes that; a two-dimensional neighbourhood
            // does. So on failure the recruitment widens (8 → 16 → 32 →
            // everything) until the verified constraints hold — for a thin
            // structure that pulls in the opposite face's nodes, which is
            // exactly the transverse spread the constraint system needs.
            // Found by the Turek–Hron FSI1 flag's deformed bottom edge.
            let mut solved = None;
            let mut count = NEIGHBOURS;
            loop {
                let recruited = nearest(structure_nodes, face.centroid, count);
                if let Some(weights) = solve_weights(structure_nodes, &recruited, face.centroid) {
                    solved = Some(recruited.into_iter().zip(weights).collect());
                    break;
                }
                if count >= structure_nodes.len() {
                    break;
                }
                count = (count * 2).min(structure_nodes.len());
            }
            let row = solved.ok_or(FsiError::DegenerateNeighbourhood { face: face_index })?;
            rows.push(row);
        }

        Ok(Self {
            rows,
            node_count: structure_nodes.len(),
        })
    }

    /// Number of fluid faces on the interface.
    #[must_use]
    pub fn face_count(&self) -> usize {
        self.rows.len()
    }

    /// Number of structure nodes on the interface.
    #[must_use]
    pub const fn node_count(&self) -> usize {
        self.node_count
    }

    /// The `(node, weight)` pairs a face distributes onto.
    #[must_use]
    pub fn weights_for(&self, face: usize) -> &[(usize, f64)] {
        self.rows.get(face).map_or(&[], Vec::as_slice)
    }

    /// Fluid tractions to structural nodal forces.
    ///
    /// Conserves total force and total moment exactly.
    ///
    /// # Errors
    /// [`FsiError::CountMismatch`] if the field lengths disagree with the
    /// interface, or [`FsiError::NonFinite`] for a NaN or infinite value.
    pub fn transfer_load(
        &self,
        fluid_faces: &[FluidFace],
        tractions: &[Vector3<f64>],
    ) -> Result<Vec<Vector3<f64>>, FsiError> {
        if fluid_faces.len() != self.rows.len() {
            return Err(FsiError::CountMismatch {
                field: "fluid_faces",
                got: fluid_faces.len(),
                expected: self.rows.len(),
            });
        }
        if tractions.len() != self.rows.len() {
            return Err(FsiError::CountMismatch {
                field: "tractions",
                got: tractions.len(),
                expected: self.rows.len(),
            });
        }
        check_finite("tractions", tractions)?;

        let mut nodal = vec![Vector3::zeros(); self.node_count];
        for ((face, traction), row) in fluid_faces.iter().zip(tractions).zip(&self.rows) {
            let force = traction * face.area;
            for (node, weight) in row {
                nodal[*node] += force * *weight;
            }
        }
        Ok(nodal)
    }

    /// Structural nodal velocities to fluid face velocities.
    ///
    /// Uses the transpose of the load operator, which is what makes
    /// interface work conserved identically.
    ///
    /// # Errors
    /// [`FsiError::CountMismatch`] or [`FsiError::NonFinite`].
    pub fn transfer_motion(
        &self,
        node_velocities: &[Vector3<f64>],
    ) -> Result<Vec<Vector3<f64>>, FsiError> {
        if node_velocities.len() != self.node_count {
            return Err(FsiError::CountMismatch {
                field: "node_velocities",
                got: node_velocities.len(),
                expected: self.node_count,
            });
        }
        check_finite("node_velocities", node_velocities)?;

        Ok(self
            .rows
            .iter()
            .map(|row| {
                row.iter()
                    .map(|(node, weight)| node_velocities[*node] * *weight)
                    .sum()
            })
            .collect())
    }
}

/// Indices of the `count` nodes nearest a point.
fn nearest(nodes: &[Vector3<f64>], point: Vector3<f64>, count: usize) -> Vec<usize> {
    let mut ordered: Vec<(f64, usize)> = nodes
        .iter()
        .enumerate()
        .map(|(index, node)| ((node - point).norm_squared(), index))
        .collect();
    ordered.sort_by(|a, b| a.0.total_cmp(&b.0));
    ordered
        .into_iter()
        .take(count.min(nodes.len()))
        .map(|(_, index)| index)
        .collect()
}

/// Minimum-norm weights satisfying partition of unity and linear
/// reproduction of `centroid`.
///
/// Solves `w = Aᵀ(AAᵀ)⁺ b` with `A` the 4×k constraint matrix, using the
/// **pseudo-inverse** rather than an inverse.
///
/// That is not defensive programming, it is the common case. A wetted
/// surface is a surface, so its nodes are usually planar — and for a
/// planar patch the `z` constraint row is an affine multiple of the ones
/// row, leaving `AAᵀ` genuinely rank-deficient. The constraint is not
/// unsatisfiable there, it is *redundant*: every partition of unity
/// reproduces a coordinate that is the same at every node. An ordinary
/// inverse would reject the most ordinary interface there is.
///
/// The constraints are then checked against the weights actually
/// obtained, because a pseudo-inverse returns a least-squares answer
/// whether or not the system was consistent. If the target genuinely
/// cannot be reproduced — a face outside the span of its recruited
/// nodes — the residual reveals it and the face is refused.
fn solve_weights(
    nodes: &[Vector3<f64>],
    recruited: &[usize],
    centroid: Vector3<f64>,
) -> Option<Vec<f64>> {
    // Centre on the face and scale by the neighbourhood radius before
    // forming the constraint Gram. The constraints are unchanged —
    // `sum w = 1` and `sum w (node - centroid)/scale = 0` is exactly
    // partition of unity plus reproduction of the centroid — but the
    // conditioning is O(1) instead of growing with (position / spacing)^2:
    // with raw coordinates, nodes near x ~ 0.26 spaced 0.005 apart put the
    // Gram's condition number past 1e4, the 4x4 SVD pseudo-inverse lost
    // enough accuracy that the verification below rejected a perfectly
    // healthy neighbourhood, and the operator's behaviour depended on WHERE
    // the interface sat — found by the Turek-Hron FSI1 flag, whose bottom
    // edge is exactly such a neighbourhood. A transfer operator must be
    // translation-invariant; centring makes it so.
    let scale = recruited
        .iter()
        .map(|index| (nodes[*index] - centroid).norm())
        .fold(0.0_f64, f64::max)
        .max(1e-300);
    let mut gram = Matrix4::zeros();
    for index in recruited {
        let local = (nodes[*index] - centroid) / scale;
        let row = Vector4::new(1.0, local.x, local.y, local.z);
        gram += row * row.transpose();
    }

    let target = Vector4::new(1.0, 0.0, 0.0, 0.0);
    let lambda = gram.pseudo_inverse(SINGULAR_TOLERANCE).ok()? * target;

    let weights: Vec<f64> = recruited
        .iter()
        .map(|index| {
            let local = (nodes[*index] - centroid) / scale;
            Vector4::new(1.0, local.x, local.y, local.z).dot(&lambda)
        })
        .collect();

    // Verify what was asked for, rather than trusting the solve.
    let unity: f64 = weights.iter().sum();
    let reproduced: Vector3<f64> = recruited
        .iter()
        .zip(&weights)
        .map(|(index, weight)| nodes[*index] * *weight)
        .sum();
    if (unity - 1.0).abs() > CONSTRAINT_TOLERANCE
        || (reproduced - centroid).norm() > CONSTRAINT_TOLERANCE
    {
        return None;
    }

    Some(weights)
}

fn check_finite(field: &'static str, values: &[Vector3<f64>]) -> Result<(), FsiError> {
    for (index, value) in values.iter().enumerate() {
        if !value.iter().all(|component| component.is_finite()) {
            return Err(FsiError::NonFinite { field, index });
        }
    }
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    use nalgebra::Vector3;

    /// A structure-side patch that deliberately does not match the fluid
    /// side: 4x4 nodes at 0.25 spacing on the z = 0 plane. Matching meshes
    /// hide every interesting failure, and a planar patch is the ordinary
    /// case a wetted surface presents.
    fn structure_nodes() -> Vec<Vector3<f64>> {
        let mut nodes = Vec::new();
        for i in 0..4 {
            for j in 0..4 {
                nodes.push(Vector3::new(f64::from(i) * 0.25, f64::from(j) * 0.25, 0.0));
            }
        }
        nodes
    }

    /// Fluid faces at 0.35 spacing, offset so no face sits on a node.
    fn fluid_faces() -> Vec<FluidFace> {
        let mut faces = Vec::new();
        for i in 0..3 {
            for j in 0..3 {
                faces.push(FluidFace {
                    centroid: Vector3::new(0.1 + f64::from(i) * 0.3, 0.1 + f64::from(j) * 0.3, 0.0),
                    normal: Vector3::new(0.0, 0.0, 1.0),
                    area: 0.09,
                });
            }
        }
        faces
    }

    // ---- the conservation properties ----

    #[test]
    fn the_weights_form_a_partition_of_unity() {
        // Sum to one is what conserves total force. Without it the
        // coupling quietly gains or loses load every step.
        let surface = WettedSurface::build(&fluid_faces(), &structure_nodes()).expect("buildable");
        for face in 0..surface.face_count() {
            let total: f64 = surface.weights_for(face).iter().map(|(_, w)| w).sum();
            assert!(
                (total - 1.0).abs() < 1e-12,
                "face {face} weights sum to {total}"
            );
        }
    }

    #[test]
    fn the_weights_reproduce_the_face_centroid() {
        // Linear reproduction. This is what conserves *moment*: the load
        // must arrive at the same place it left, in the weighted-average
        // sense. Partition of unity alone conserves force and silently
        // corrupts the moment.
        let nodes = structure_nodes();
        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
        for (face, expected) in fluid_faces().iter().enumerate() {
            let reproduced: Vector3<f64> = surface
                .weights_for(face)
                .iter()
                .map(|(node, w)| nodes[*node] * *w)
                .sum();
            let error = (reproduced - expected.centroid).norm();
            assert!(error < 1e-10, "face {face} centroid off by {error}");
        }
    }

    #[test]
    fn total_force_is_preserved_across_the_interface() {
        let nodes = structure_nodes();
        let faces = fluid_faces();
        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");

        // A non-uniform traction, so a bug cannot hide behind symmetry.
        let tractions: Vec<Vector3<f64>> = (0..faces.len())
            .map(|i| Vector3::new(1.0 + f64::from(i as i32), -2.0, 0.5))
            .collect();

        let fluid_total: Vector3<f64> = faces
            .iter()
            .zip(&tractions)
            .map(|(face, traction)| traction * face.area)
            .sum();
        let nodal = surface
            .transfer_load(&faces, &tractions)
            .expect("transferable");
        let structure_total: Vector3<f64> = nodal.iter().sum();

        let error = (structure_total - fluid_total).norm();
        assert!(error < 1e-10, "force lost across interface: {error}");
    }

    #[test]
    fn total_moment_is_preserved_across_the_interface() {
        let nodes = structure_nodes();
        let faces = fluid_faces();
        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
        let tractions: Vec<Vector3<f64>> = (0..faces.len())
            .map(|i| Vector3::new(1.0 + f64::from(i as i32), -2.0, 0.5))
            .collect();

        // About an arbitrary point, deliberately not the origin: a scheme
        // that only conserves moment about one special point is not
        // conserving moment.
        let about = Vector3::new(-0.7, 1.3, 2.1);
        let fluid_moment: Vector3<f64> = faces
            .iter()
            .zip(&tractions)
            .map(|(face, traction)| (face.centroid - about).cross(&(traction * face.area)))
            .sum();

        let nodal = surface
            .transfer_load(&faces, &tractions)
            .expect("transferable");
        let structure_moment: Vector3<f64> = nodes
            .iter()
            .zip(&nodal)
            .map(|(position, force)| (position - about).cross(force))
            .sum();

        let error = (structure_moment - fluid_moment).norm();
        assert!(error < 1e-10, "moment lost across interface: {error}");
    }

    #[test]
    fn interface_work_is_conserved() {
        // The decisive one. Motion transfer uses the transpose of the load
        // transfer, so the work the fluid does equals the work the
        // structure receives, exactly. A coupling that fails this injects
        // or drains energy every step and will eventually destabilise for
        // reasons that look like physics.
        let nodes = structure_nodes();
        let faces = fluid_faces();
        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");

        let tractions: Vec<Vector3<f64>> = (0..faces.len())
            .map(|i| Vector3::new(0.3 * f64::from(i as i32), 1.7, -0.4))
            .collect();
        let node_velocities: Vec<Vector3<f64>> = (0..nodes.len())
            .map(|i| Vector3::new(0.1, 0.2 * f64::from(i as i32), -0.05))
            .collect();

        let nodal_forces = surface
            .transfer_load(&faces, &tractions)
            .expect("transferable");
        let face_velocities = surface
            .transfer_motion(&node_velocities)
            .expect("transferable");

        let fluid_work: f64 = faces
            .iter()
            .zip(&tractions)
            .zip(&face_velocities)
            .map(|((face, traction), velocity)| (traction * face.area).dot(velocity))
            .sum();
        let structure_work: f64 = nodal_forces
            .iter()
            .zip(&node_velocities)
            .map(|(force, velocity)| force.dot(velocity))
            .sum();

        assert!(
            (fluid_work - structure_work).abs() < 1e-10,
            "interface work not conserved: fluid {fluid_work}, structure {structure_work}"
        );
    }

    // ---- limits ----

    #[test]
    fn zero_traction_moves_nothing() {
        let nodes = structure_nodes();
        let faces = fluid_faces();
        let surface = WettedSurface::build(&faces, &nodes).expect("buildable");
        let tractions = vec![Vector3::zeros(); faces.len()];
        let nodal = surface
            .transfer_load(&faces, &tractions)
            .expect("transferable");
        assert!(nodal.iter().all(|force| force.norm() < 1e-15));
    }

    #[test]
    fn a_motionless_structure_leaves_the_fluid_boundary_still() {
        let nodes = structure_nodes();
        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
        let velocities = vec![Vector3::zeros(); nodes.len()];
        let face_velocities = surface.transfer_motion(&velocities).expect("transferable");
        assert!(face_velocities.iter().all(|v| v.norm() < 1e-15));
    }

    #[test]
    fn a_rigid_translation_transfers_unchanged() {
        // Every structure node moving together must give every fluid face
        // the same velocity. This follows from partition of unity, and it
        // is the check that catches a normalisation slip.
        let nodes = structure_nodes();
        let surface = WettedSurface::build(&fluid_faces(), &nodes).expect("buildable");
        let rigid = Vector3::new(0.4, -1.1, 0.9);
        let velocities = vec![rigid; nodes.len()];

        for velocity in surface.transfer_motion(&velocities).expect("transferable") {
            assert!(
                (velocity - rigid).norm() < 1e-12,
                "rigid translation distorted: {velocity:?}"
            );
        }
    }

    // ---- refusals ----

    #[test]
    fn an_interface_with_too_few_nodes_is_refused() {
        // Linear reproduction in three dimensions needs four independent
        // nodes. Fewer cannot satisfy the constraint, and pretending
        // otherwise would silently break moment conservation.
        let nodes: Vec<Vector3<f64>> = (0..3)
            .map(|i| Vector3::new(f64::from(i), f64::from(i) * 0.5, 0.0))
            .collect();
        assert!(matches!(
            WettedSurface::build(&fluid_faces(), &nodes),
            Err(FsiError::InsufficientNodes { .. })
        ));
    }

    #[test]
    fn an_empty_interface_is_refused() {
        assert!(WettedSurface::build(&[], &structure_nodes()).is_err());
        assert!(WettedSurface::build(&fluid_faces(), &[]).is_err());
    }

    #[test]
    fn a_traction_count_mismatch_is_refused() {
        let surface = WettedSurface::build(&fluid_faces(), &structure_nodes()).expect("buildable");
        let wrong = vec![Vector3::new(1.0, 0.0, 0.0); 2];
        assert!(matches!(
            surface.transfer_load(&fluid_faces(), &wrong),
            Err(FsiError::CountMismatch { .. })
        ));
    }

    #[test]
    fn a_non_finite_traction_is_refused() {
        let faces = fluid_faces();
        let surface = WettedSurface::build(&faces, &structure_nodes()).expect("buildable");
        let mut tractions = vec![Vector3::new(1.0, 0.0, 0.0); faces.len()];
        tractions[1].y = f64::NAN;
        assert!(matches!(
            surface.transfer_load(&faces, &tractions),
            Err(FsiError::NonFinite { .. })
        ));
    }
}