feat(linear-static-mitc4-shell): step 12 - shell-linear-static-flow

src/fesa/results/result_recovery.cpp

@@ -85,6 +85,39 @@ double indexedNorm(const Vector& values,
     return result;
 }
 
+std::array<double, 2> freeEquationInternalTermNorms(
+    const SparseMatrix& stiffness,
+    const Vector& displacement,
+    const DofManager& dofs) {
+    std::vector<unsigned char> freeColumn(dofs.fullDofCount(), 0U);
+    for (const std::size_t fullDof : dofs.freeDofs()) {
+        freeColumn[fullDof] = 1U;
+    }
+
+    double freeTermNorm = 0.0;
+    double constrainedTermNorm = 0.0;
+    for (const std::size_t row : dofs.freeDofs()) {
+        double freeTerm = 0.0;
+        double constrainedTerm = 0.0;
+        for (std::size_t position = stiffness.rowOffsets()[row];
+             position < stiffness.rowOffsets()[row + 1U];
+             ++position) {
+            const std::size_t column = stiffness.columnIndices()[position];
+            const double contribution =
+                stiffness.values()[position] * displacement[column];
+            if (freeColumn[column] != 0U) {
+                freeTerm += contribution;
+            } else {
+                constrainedTerm += contribution;
+            }
+        }
+        freeTermNorm = std::hypot(freeTermNorm, freeTerm);
+        constrainedTermNorm =
+            std::hypot(constrainedTermNorm, constrainedTerm);
+    }
+    return {freeTermNorm, constrainedTermNorm};
+}
+
 bool strictlyIncreasing(const std::vector<std::size_t>& values) {
     return std::adjacent_find(
                values.begin(), values.end(),
@@ -639,10 +672,15 @@ Status ResultRecovery::recover(const AnalysisModel& model,
     }
 
     const double residualNorm = indexedNorm(residual, dofs.freeDofs());
-    const double internalNorm = indexedNorm(internalForce, dofs.freeDofs());
+    const auto internalTermNorms = freeEquationInternalTermNorms(
+        fullStiffness, state.displacement(), dofs);
     const double externalNorm =
         indexedNorm(state.externalForce(), dofs.freeDofs());
-    const double denominator = (std::max)(internalNorm, externalNorm);
+    // Normalize against the three terms of
+    // Kff*df + Kfc*dc - Ff. Using only the already-cancelled K*d term would
+    // classify prescribed-only equilibrium roundoff as a unit residual.
+    const double denominator = (std::max)({
+        internalTermNorms[0U], internalTermNorms[1U], externalNorm});
     if (!std::isfinite(residualNorm) || !std::isfinite(denominator)) {
         return recoveryFailure(
             "nonfinite-recovery-value",