feat(linear-static-mitc4-shell): step 8 - shell-load-validation

src/fesa/assembly/load_assembler.cpp

@@ -17,6 +17,7 @@ namespace fesa {
 namespace {
 
 constexpr std::size_t dofsPerNode = 6U;
+constexpr double shellMomentProjectionTolerance = 1.0e-12;
 
 Status loadFailure(
     const std::string& code,
@@ -215,6 +216,70 @@ Status validateFiniteVector(
     return Status::ok();
 }
 
+Status validateShellMoments(
+    const Domain& domain,
+    const Vector& fullLoad) {
+    if (domain.shellElements().empty()) {
+        return Status::ok();
+    }
+
+    std::vector<const ShellNodeInitialFrame*> frameByNode(
+        domain.nodes().size(), nullptr);
+    for (const auto& frame : domain.shellNodeInitialFrames()) {
+        if (frame.nodeIndex >= frameByNode.size() ||
+            frameByNode[frame.nodeIndex] != nullptr) {
+            return loadFailure(
+                "invalid-shell-director",
+                {domain.sourcePath(), 0U},
+                "NODE",
+                std::to_string(frame.nodeIndex),
+                "Shell nodal directors must have unique in-range node identities.");
+        }
+        frameByNode[frame.nodeIndex] = &frame;
+    }
+
+    for (std::size_t node = 0U; node < domain.nodes().size(); ++node) {
+        const double momentX = fullLoad[node * dofsPerNode + 3U];
+        const double momentY = fullLoad[node * dofsPerNode + 4U];
+        const double momentZ = fullLoad[node * dofsPerNode + 5U];
+        if (momentX == 0.0 && momentY == 0.0 && momentZ == 0.0) {
+            continue;
+        }
+
+        const auto* const frame = frameByNode[node];
+        if (frame == nullptr) {
+            return loadFailure(
+                "invalid-shell-director",
+                domain.nodes()[node].location,
+                "NODE",
+                domain.nodes()[node].sourceId.sourceLabelText,
+                "A loaded shell node must have an approved initial director.");
+        }
+
+        const double momentScale = std::max(
+            std::abs(momentX),
+            std::max(std::abs(momentY), std::abs(momentZ)));
+        const double scaledX = momentX / momentScale;
+        const double scaledY = momentY / momentScale;
+        const double scaledZ = momentZ / momentScale;
+        const double scaledNorm = std::hypot(scaledX, scaledY, scaledZ);
+        const double scaledDot =
+            frame->director[0U] * scaledX +
+            frame->director[1U] * scaledY +
+            frame->director[2U] * scaledZ;
+        const double projectionRatio = std::abs(scaledDot) / scaledNorm;
+        if (!(projectionRatio <= shellMomentProjectionTolerance)) {
+            return loadFailure(
+                "unsupported-drilling-load",
+                domain.nodes()[node].location,
+                "CLOAD",
+                domain.nodes()[node].sourceId.sourceLabelText,
+                "The aggregate nodal moment has an unsupported director-parallel component.");
+        }
+    }
+    return Status::ok();
+}
+
 } // namespace
 
 Result<Vector> LoadAssembler::assembleFullNodalLoad(
@@ -328,6 +393,10 @@ Result<Vector> LoadAssembler::assembleFullNodalLoad(
             fullLoad[fullDof] = accumulated;
         }
     }
+    const Status shellMomentStatus = validateShellMoments(domain, fullLoad);
+    if (!shellMomentStatus.isOk()) {
+        return Result<Vector>::failure(shellMomentStatus);
+    }
     return Result<Vector>::success(std::move(fullLoad));
 }