Add shell-to-solid MPC validation and refactor coupling logic

engine/cpp/shell_core.cpp

@@ -9,7 +9,6 @@
 #include <array>
 #include <cmath>
 #include <cstdio>
-#include <limits>
 #include <map>
 #include <stdexcept>
 #include <string>
@@ -809,48 +808,72 @@ Rbe3Constraints build_rbe3_constraints(const CoupledInput& in, int nDof,
         for (size_t k = 0; k < 6; ++k) out.dep[R6 + k] = 1;
 
         if (cp.mpc) {
-            // Relaxed shell-to-solid MPC (Lu, Zhang & Yang 2023): rigid translation
-            // tie to the nearest solid node S — enforcing displacement CONTINUITY at
-            // the coincident junction — plus a ψ-scaled least-squares rotation of the
-            // coupled solid nodes about S (the unstructured generalisation of the
-            // paper's finite-difference rotation stencil). See Coupling::mpc.
+            // Relaxed shell-to-solid MPC (Lu, Zhang & Yang 2023), generalised to a
+            // NON-CONFORMING seam. The paper ties the shell node to a solid node
+            // coincident with it; an auto-detected seam has none, and tying to the
+            // nearest partner instead welds the mid-surface to one FACE of the very
+            // wall it idealises. The mid-surface sits t/2 from the footprint node on
+            // either face, so with both partners equidistant the pick is decided by
+            // floating-point noise and flips from node to node along the seam. Each
+            // flip is a spurious ±t/2 eccentricity, and an eccentric tie carries a
+            // moment the structure does not: on a 2 mm wall the shell facets touching
+            // a tied node ran 12× the von Mises of their neighbours, peaking exactly
+            // where the sign flipped (KOF-212).
+            //
+            // Anchor the tie on the patch's distance-weighted centroid instead of on
+            // one node of it. For the equidistant pair that is precisely the
+            // mid-surface, so the eccentricity cancels; a genuinely coincident node
+            // still dominates the kernel and recovers the paper's tie. The 1/(d²+ε²)
+            // weighting keeps the tie LOCAL, which is what separates it from the
+            // equal-weight distributing coupling — that one averages over the whole
+            // search ball and smears the junction stress away instead.
             const double psi = cp.relaxation;
-            int S = cp.solid_nodes[0];
-            double best = std::numeric_limits<double>::max();
-            for (int sn : cp.solid_nodes) {
-                const Vec3 d = sub(vtx(sn), pR);
-                const double dd = dot(d, d);
-                if (dd < best) { best = dd; S = sn; }
+            double rmax = 0.0;
+            for (int sn : cp.solid_nodes) rmax = std::max(rmax, norm(sub(vtx(sn), pR)));
+            if (rmax <= 0.0)
+                throw std::runtime_error(
+                    "coupled: every coupled node of the shell-to-solid MPC at node " +
+                    std::to_string(R) +
+                    " sits on the reference point, so the patch has no extent and no "
+                    "rotation can be measured from it.");
+            // ε softens the kernel at d → 0 and is a fraction of the patch's own
+            // radius, so the weighting is scale-free.
+            const double eps2 = (0.05 * rmax) * (0.05 * rmax);
+            double W = 0.0;
+            for (int i = 0; i < N; ++i) {
+                const Vec3 d = sub(vtx(cp.solid_nodes[i]), pR);
+                w[i] /= dot(d, d) + eps2;
+                W += w[i];
+            }
+            // Σ wᵢ rᵢ = 0 about the weighted centroid — that is what drops the u_A
+            // term out of the least-squares rotation below.
+            Vec3 pA = {0.0, 0.0, 0.0};
+            for (int i = 0; i < N; ++i) {
+                const Vec3 pi = vtx(cp.solid_nodes[i]);
+                for (int k = 0; k < 3; ++k) pA[k] += (w[i] / W) * pi[k];
             }
-            const Vec3 pS = vtx(S);
             std::vector<Vec3> r(N);
-            Vec3 Sp = {0.0, 0.0, 0.0};  // Σ w_i r_i, r_i measured from S
             std::array<std::array<double, 3>, 3> H{};
             for (int i = 0; i < N; ++i) {
-                const Vec3 pi = vtx(cp.solid_nodes[i]);
-                r[i] = {pi[0] - pS[0], pi[1] - pS[1], pi[2] - pS[2]};
-                for (int k = 0; k < 3; ++k) Sp[k] += w[i] * r[i][k];
+                r[i] = sub(vtx(cp.solid_nodes[i]), pA);
                 const double rr = dot(r[i], r[i]);
                 for (int a = 0; a < 3; ++a)
                     for (int b = 0; b < 3; ++b)
                         H[a][b] += w[i] * ((a == b ? rr : 0.0) - r[i][a] * r[i][b]);
             }
             const auto Hinv = mat3_inv(H);
-            // U_R = u_S  (translations follow the coincident solid node exactly).
-            for (int c = 0; c < 3; ++c) out.Cmap[R6 + c].emplace_back(6 * S + c, 1.0);
-            // Θ_R = ψ · Hinv · Σ w_i [r_i]× (u_i − u_S). Split u_i and u_S parts:
-            // the [S']× u_S term collects the −u_S contribution (skew is linear).
+            // U_R = Σ (wᵢ/W) uᵢ — the patch's weighted motion at the anchor.
+            for (int i = 0; i < N; ++i)
+                for (int c = 0; c < 3; ++c)
+                    out.Cmap[R6 + c].emplace_back(6 * cp.solid_nodes[i] + c, w[i] / W);
+            // Θ_R = ψ · Hinv · Σ wᵢ [rᵢ]× uᵢ.
             for (int i = 0; i < N; ++i) {
                 const int sn = cp.solid_nodes[i];
                 const auto M = mat3_mul(Hinv, skew(r[i]));
                 for (int a = 0; a < 3; ++a)
                     for (int c = 0; c < 3; ++c)
                         out.Cmap[R6 + 3 + a].emplace_back(6 * sn + c, psi * w[i] * M[a][c]);
             }
-            const auto Ms = mat3_mul(Hinv, skew(Sp));
-            for (int a = 0; a < 3; ++a)
-                for (int c = 0; c < 3; ++c)
-                    out.Cmap[R6 + 3 + a].emplace_back(6 * S + c, -psi * Ms[a][c]);
             continue;
         }
 

engine/cpp/shell_core.h

@@ -92,11 +92,23 @@ struct Coupling {
     std::vector<double> weights;         // per solid node (empty ⇒ equal weights)
     // Relaxed shell-to-solid MPC (Lu, Zhang & Yang, "A Relaxed MPC Method for
     // Non-rigid Shell to Solid Coupling", J. Phys.: Conf. Ser. 2528 012064, 2023).
-    // When true, the reference (shell) node's translations are rigidly tied to its
-    // nearest solid node — a coincident-node tie that enforces displacement
-    // CONTINUITY at the junction (unlike the distributing RBE3 average, which only
-    // matches the resultant and lets the shell separate) — and its rotations follow
-    // the relaxation-scaled least-squares rotation of the coupled solid nodes.
+    // When true, the reference (shell) node's translations are tied to the
+    // DISTANCE-WEIGHTED CENTROID of its solid patch, and its rotations follow the
+    // relaxation-scaled least-squares rotation of that patch. The 1/(d²+ε²) weight
+    // keeps the tie local, so it enforces displacement CONTINUITY at the junction —
+    // unlike the distributing RBE3 average, which spreads over the whole search
+    // ball, matches only the resultant, and lets the shell separate.
+    //
+    // The paper ties to a solid node COINCIDENT with the shell node. An
+    // auto-detected seam has none: the mid-surface sits t/2 from the footprint node
+    // on either wall face, equidistant to the last bit, so tying to "the nearest"
+    // picks a side by floating-point accident and welds the mid-surface to one face
+    // of the wall it idealises. That ±t/2 eccentricity carries a moment the
+    // structure does not, and it flips sign wherever the pick flips (KOF-212). The
+    // weighted centroid of the equidistant pair IS the mid-surface, so the
+    // eccentricity cancels; a genuinely coincident node still dominates the kernel
+    // and recovers the paper's tie.
+    //
     // Appropriate for a continuous-material seam (a thin wall idealised as shell,
     // tied back to the retained solid); the distributing coupling (mpc = false)
     // stays for genuinely gapped, non-conforming interfaces (a pin in a hole).

engine/tests/shell_validation.cpp

@@ -14,6 +14,7 @@
 
 #include "shell_core.h"
 
+#include <algorithm>
 #include <array>
 #include <cmath>
 #include <cstdio>
@@ -263,6 +264,110 @@ double coupled_mpc_clamped(double psi) {
     return w / ref;
 }
 
+// Shell-to-solid MPC on a NON-CONFORMING seam: the wall footprint has solid
+// nodes on both faces (y = W/2 ± t/2) and the shell mid-surface node exactly
+// between them, so the two candidate partners are equidistant to the last bit.
+// Which one an unstructured mesh lists first is an accident, so the solved
+// stress must not depend on it — anchoring the tie on one of them instead of on
+// the mid-surface carries a ±t/2 eccentricity, and a moment with it (KOF-212).
+//
+// `order`: 0 = the y−t/2 node is listed first for every seam node, 1 = y+t/2,
+// 2 = alternating along the seam (what a real mesh produces). Returns the peak
+// shell von Mises; the three must agree.
+double mpc_seam_peak(int order) {
+    const double L = 60, W = 40, H = 6, Hw = 30, t = 2.0;
+    const double E = 210e3, nu = 0.3, P = 5000.0;
+    const int nx = 12, nz = 2, nzw = 10;
+
+    std::vector<double> ys;
+    for (int j = 0; j <= 8; ++j) ys.push_back(W * j / 8);
+    ys.push_back(W / 2 - t / 2);
+    ys.push_back(W / 2 + t / 2);
+    ys.erase(std::remove_if(ys.begin(), ys.end(),
+                            [&](double y) { return std::fabs(y - W / 2) < 1e-12; }),
+             ys.end());
+    std::sort(ys.begin(), ys.end());
+    const int ny = (int)ys.size() - 1;
+
+    std::vector<double> V;
+    std::vector<int> tets;
+    auto sid = [&](int i, int j, int k) { return i * (ny + 1) * (nz + 1) + j * (nz + 1) + k; };
+    for (int i = 0; i <= nx; ++i)
+        for (int j = 0; j <= ny; ++j)
+            for (int k = 0; k <= nz; ++k) {
+                V.push_back(L * i / nx); V.push_back(ys[j]); V.push_back(H * k / nz);
+            }
+    for (int i = 0; i < nx; ++i)
+        for (int j = 0; j < ny; ++j)
+            for (int k = 0; k < nz; ++k) {
+                const int a=sid(i,j,k),b=sid(i+1,j,k),c=sid(i+1,j+1,k),d=sid(i,j+1,k),
+                          e=sid(i,j,k+1),f=sid(i+1,j,k+1),g=sid(i+1,j+1,k+1),h=sid(i,j+1,k+1);
+                const std::array<std::array<int,4>,6> q = {{
+                    {a,b,c,g},{a,c,d,g},{a,d,h,g},{a,h,e,g},{a,e,f,g},{a,f,b,g}}};
+                for (const auto& tt : q) for (int m = 0; m < 4; ++m) tets.push_back(tt[m]);
+            }
+    const int nSolid = (int)(V.size() / 3);
+
+    // Shell wall hanging below the block on the mid-surface y = W/2.
+    const int base = nSolid;
+    auto wid = [&](int i, int k) { return base + i * (nzw + 1) + k; };
+    for (int i = 0; i <= nx; ++i)
+        for (int k = 0; k <= nzw; ++k) {
+            V.push_back(L * i / nx); V.push_back(W / 2); V.push_back(-Hw * k / nzw);
+        }
+    std::vector<int> tris;
+    for (int i = 0; i < nx; ++i)
+        for (int k = 0; k < nzw; ++k) {
+            const int a=wid(i,k),b=wid(i+1,k),c=wid(i+1,k+1),d=wid(i,k+1);
+            tris.push_back(a); tris.push_back(b); tris.push_back(c);
+            tris.push_back(a); tris.push_back(c); tris.push_back(d);
+        }
+
+    CoupledInput in;
+    in.n_nodes = (int)(V.size() / 3);
+    in.vertices = V; in.triangles = tris;
+    in.shell_young = E; in.shell_poisson = nu; in.thickness = t;
+    in.solid_stiffness = tet_solid_stiffness(V, tets, E, nu);
+
+    const double radius = 1.6 * (L / nx);
+    for (int i = 0; i <= nx; ++i) {
+        const int rn = wid(i, 0);
+        std::vector<int> patch;
+        for (int sn = 0; sn < nSolid; ++sn) {
+            const size_t bs = 3 * static_cast<size_t>(sn), br = 3 * static_cast<size_t>(rn);
+            const double dx = V[bs] - V[br], dy = V[bs+1] - V[br+1], dz = V[bs+2] - V[br+2];
+            if (dx*dx + dy*dy + dz*dz <= radius*radius) patch.push_back(sn);
+        }
+        if (patch.size() < 3) continue;
+        const bool plusFirst = (order == 1) || (order == 2 && (i % 2 == 0));
+        std::stable_sort(patch.begin(), patch.end(), [&](int a, int b) {
+            auto key = [&](int n) {
+                const size_t b = 3 * static_cast<size_t>(n);
+                const double dy = V[b+1] - W / 2;
+                if (std::fabs(std::fabs(dy) - t/2) > 1e-9 || std::fabs(V[b+2]) > 1e-9) return 2;
+                return (dy > 0) == plusFirst ? 0 : 1;
+            };
+            return key(a) < key(b);
+        });
+        Coupling cp;
+        cp.ref_node = rn; cp.solid_nodes = patch;
+        cp.kind = CouplingKind::RelaxedMpc; cp.mpc = true; cp.relaxation = 1.0;
+        in.couplings.push_back(std::move(cp));
+    }
+    // Cantilever the block in its own plane so the seam sees a transverse
+    // gradient ∂u/∂y — what an eccentric tie converts into spurious bending.
+    for (int j = 0; j <= ny; ++j)
+        for (int k = 0; k <= nz; ++k)
+            for (int c = 0; c < 3; ++c) in.fixed_dofs.push_back(6 * sid(0, j, k) + c);
+    for (int j = 0; j <= ny; ++j)
+        for (int k = 0; k <= nz; ++k)
+            in.loads.emplace_back(6 * sid(nx, j, k) + 1, P / ((ny + 1.0) * (nz + 1.0)));
+
+    ShellResult r = solve_solid_shell_core(in);
+    const std::vector<double> vm = shell_von_mises(V, tris, t, {}, E, nu, r.dofs);
+    return *std::max_element(vm.begin(), vm.end());
+}
+
 // A fixed-DOF constraint that lands on an RBE3 coupling reference node (a
 // dependent, shell-side DOF) must be rejected loudly, not silently dropped
 // (issue #377). Reuses the cantilever-on-anchors coupling: every root node is a
@@ -325,6 +430,21 @@ int main() {
     check(failures, "coupled-mpc-clamped-rigid", coupled_mpc_clamped(1.0), 1.0, 6.0);
     check(failures, "coupled-mpc-clamped-relaxed", coupled_mpc_clamped(0.5), 1.0, 6.0);
 
+    printf("Non-conforming shell-to-solid seam (KOF-212):\n");
+    {
+        const double a0 = mpc_seam_peak(0), a1 = mpc_seam_peak(1), a2 = mpc_seam_peak(2);
+        const double lo = std::min({a0, a1, a2}), hi = std::max({a0, a1, a2});
+        const double spread = (hi - lo) / lo * 100.0;
+        // The tie anchors on the mid-surface, so which equidistant footprint node
+        // the patch happens to list first cannot change the answer at all. Before
+        // the fix these came out 24.6 / 26.3 / 63.7 — a 159 % spread, worst for
+        // the alternating pick a real mesh produces.
+        const bool ok = spread <= 0.1;
+        if (!ok) ++failures;
+        printf("  [%s] %-28s peak vM %.4e / %.4e / %.4e  spread %.2f%% (tol 0.1%%)\n",
+               ok ? "PASS" : "FAIL", "seam-tie-order-invariant", a0, a1, a2, spread);
+    }
+
     printf("Constraint-on-dependent-node rejection (issue #377):\n");
     {
         const bool threw = coupled_fixed_dependent_throws();