diff --git a/include/fesa/results/result_recovery.hpp b/include/fesa/results/result_recovery.hpp index 2344628..f85b7f6 100644 --- a/include/fesa/results/result_recovery.hpp +++ b/include/fesa/results/result_recovery.hpp @@ -17,7 +17,7 @@ struct NodeStationResultRow { std::array sectionResultant; }; -// Recovers full-space equilibrium and the V0 beam output rows without +// Recovers full-space equilibrium and the active concrete element rows without // exposing element or sparse-backend details to result consumers. class ResultRecovery { public: diff --git a/src/fesa/results/result_recovery.cpp b/src/fesa/results/result_recovery.cpp index 840c8f5..6623d8b 100644 --- a/src/fesa/results/result_recovery.cpp +++ b/src/fesa/results/result_recovery.cpp @@ -1,6 +1,7 @@ #include "fesa/results/result_recovery.hpp" #include "fesa/elements/euler_beam_3d.hpp" +#include "fesa/elements/mitc4_shell.hpp" #include #include @@ -21,7 +22,10 @@ namespace { constexpr std::size_t kDofsPerNode = 6U; constexpr std::size_t kElementDofCount = 12U; +constexpr std::size_t kShellElementDofCount = 24U; +constexpr std::size_t kShellLocationCount = 4U; constexpr double kFreeResidualTolerance = 1.0e-10; +constexpr double kGlobalEquilibriumTolerance = 1.0e-10; constexpr double kAxisTolerance = 1.0e-12; using AxisSet = std::array, 3>; @@ -56,13 +60,8 @@ bool sameSourceIdentity(const SourceEntityId& left, left.sourceLabelText == right.sourceLabelText; } -bool finite(const std::array& values) { - return std::all_of( - values.begin(), values.end(), - [](const double value) { return std::isfinite(value); }); -} - -bool finite(const std::array& values) { +template +bool finite(const std::array& values) { return std::all_of( values.begin(), values.end(), [](const double value) { return std::isfinite(value); }); @@ -252,6 +251,74 @@ Status validateRecoveryInputs(const AnalysisModel& model, "Every active element requires one twelve-DOF scatter map."); } } + + if (!model.activeElements().empty() && !domain.shellElements().empty()) { + return recoveryFailure( + "unsupported-mixed-element-model", + {domain.sourcePath(), 0U}, + "B33:FESA-MITC4", + "Result recovery does not support mixed beam and shell models."); + } + if (domain.shellElements().size() > + static_cast( + (std::numeric_limits::max)())) { + return recoveryFailure( + "invalid-recovery-dimensions", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "The shell element count cannot be represented by stable element identities."); + } + for (std::size_t elementOrder = 0U; + elementOrder < domain.shellElements().size(); + ++elementOrder) { + const auto& definition = domain.shellElements()[elementOrder]; + if (definition.materialIndex >= domain.materials().size() || + definition.sectionIndex >= domain.shellSections().size()) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Active shell material and section references must resolve before recovery."); + } + try { + const auto& scatter = dofs.shellElementScatter( + static_cast(elementOrder)); + for (std::size_t nodePosition = 0U; + nodePosition < definition.nodeIndices.size(); + ++nodePosition) { + const EntityIndex node = definition.nodeIndices[nodePosition]; + if (node >= domain.nodes().size()) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Active shell node references must resolve before recovery."); + } + for (std::size_t component = 0U; + component < kDofsPerNode; + ++component) { + const std::size_t local = + nodePosition * kDofsPerNode + component; + const std::size_t expected = + static_cast(node) * kDofsPerNode + + component; + if (scatter[local] != expected || expected >= fullCount) { + return recoveryFailure( + "invalid-recovery-order", + definition.location, + definition.sourceId.sourceLabelText, + "Shell scatter must preserve node/component full-DOF order."); + } + } + } + } catch (const std::out_of_range&) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Every active shell requires one twenty-four-DOF scatter map."); + } + } return Status::ok(); } @@ -346,6 +413,144 @@ double norm(const std::array& value) { return std::hypot(value[0U], value[1U], value[2U]); } +bool accumulateVectorAndScale( + std::array& total, + double& scale, + const std::array& contribution) { + const double magnitude = norm(contribution); + const double accumulatedScale = scale + magnitude; + if (!finite(contribution) || !std::isfinite(magnitude) || + !std::isfinite(accumulatedScale)) { + return false; + } + for (std::size_t component = 0U; + component < contribution.size(); + ++component) { + const double accumulated = total[component] + contribution[component]; + if (!std::isfinite(accumulated)) { + return false; + } + total[component] = accumulated; + } + scale = accumulatedScale; + return true; +} + +double normalizedBalance( + const std::array& balance, + const double scale) { + const double balanceNorm = norm(balance); + if (!std::isfinite(balanceNorm) || !std::isfinite(scale)) { + return (std::numeric_limits::infinity)(); + } + if (scale == 0.0) { + return balanceNorm == 0.0 + ? 0.0 + : (std::numeric_limits::infinity)(); + } + return balanceNorm / scale; +} + +Status populateShellGlobalEvidence( + const Domain& domain, + const DofManager& dofs, + const Vector& externalForce, + const Vector& residual, + const double normalizedResidual, + ShellStateCandidate& candidate) { + std::vector constrained(dofs.fullDofCount(), 0U); + for (const std::size_t fullDof : dofs.constrainedDofs()) { + constrained[fullDof] = 1U; + } + + std::array appliedForce{}; + std::array reactionForce{}; + std::array appliedMoment{}; + std::array reactionMoment{}; + double appliedForceScale = 0.0; + double reactionForceScale = 0.0; + double appliedMomentScale = 0.0; + double reactionMomentScale = 0.0; + for (std::size_t node = 0U; node < domain.nodes().size(); ++node) { + const std::size_t offset = node * kDofsPerNode; + std::array nodalAppliedForce{}; + std::array nodalReactionForce{}; + std::array nodalAppliedMoment{}; + std::array nodalReactionMoment{}; + for (std::size_t component = 0U; component < 3U; ++component) { + nodalAppliedForce[component] = externalForce[offset + component]; + nodalAppliedMoment[component] = + externalForce[offset + 3U + component]; + if (constrained[offset + component] != 0U) { + nodalReactionForce[component] = residual[offset + component]; + } + if (constrained[offset + 3U + component] != 0U) { + nodalReactionMoment[component] = + residual[offset + 3U + component]; + } + } + + const auto appliedForceMoment = + cross(domain.nodes()[node].coordinates, nodalAppliedForce); + const auto reactionForceMoment = + cross(domain.nodes()[node].coordinates, nodalReactionForce); + for (std::size_t component = 0U; component < 3U; ++component) { + nodalAppliedMoment[component] += appliedForceMoment[component]; + nodalReactionMoment[component] += reactionForceMoment[component]; + } + if (!accumulateVectorAndScale( + appliedForce, appliedForceScale, nodalAppliedForce) || + !accumulateVectorAndScale( + reactionForce, reactionForceScale, nodalReactionForce) || + !accumulateVectorAndScale( + appliedMoment, appliedMomentScale, nodalAppliedMoment) || + !accumulateVectorAndScale( + reactionMoment, reactionMomentScale, nodalReactionMoment)) { + return recoveryFailure( + "nonfinite-recovery-value", + domain.nodes()[node].location, + domain.nodes()[node].sourceId.sourceLabelText, + "Global force and moment evidence must remain finite in source-node order."); + } + } + + std::array forceBalance{}; + std::array momentBalance{}; + for (std::size_t component = 0U; component < 3U; ++component) { + forceBalance[component] = + appliedForce[component] + reactionForce[component]; + momentBalance[component] = + appliedMoment[component] + reactionMoment[component]; + candidate.equilibrium[component] = forceBalance[component]; + candidate.equilibrium[3U + component] = momentBalance[component]; + } + const double forceMetric = normalizedBalance( + forceBalance, + (std::max)(appliedForceScale, reactionForceScale)); + const double momentMetric = normalizedBalance( + momentBalance, + (std::max)(appliedMomentScale, reactionMomentScale)); + candidate.verificationMetrics = { + normalizedResidual, forceMetric, momentMetric}; + if (!finite(candidate.equilibrium) || + !finite(candidate.verificationMetrics)) { + return recoveryFailure( + "nonfinite-recovery-value", + {domain.sourcePath(), 0U}, + "global-equilibrium", + "Global equilibrium values and their physical normalization scales must be finite."); + } + if (forceMetric > kGlobalEquilibriumTolerance || + momentMetric > kGlobalEquilibriumTolerance) { + return recoveryFailure( + "global-equilibrium-tolerance-failure", + {domain.sourcePath(), 0U}, + "global-equilibrium", + "Normalized global force or moment balance exceeds 1e-10."); + } + return Status::ok(); +} + std::optional localAxes(const Domain& domain, const EulerBeam3DDefinition& element) { const auto& first = domain.nodes()[element.nodeIndices[0U]].coordinates; @@ -545,14 +750,156 @@ Status ResultRecovery::recover(const AnalysisModel& model, } } - // Commit only after all validation and element recovery succeeds so a - // failed recovery cannot leave a partially updated AnalysisState. - state.internalForce() = std::move(internalForce); - state.residual() = std::move(residual); - state.reaction() = std::move(reaction); - state.endpointResults() = std::move(endpointRows); - state.gaussResults() = std::move(gaussRows); - state.stressResults() = std::move(stressRows); + ShellStateCandidate shellCandidate{}; + std::vector expectedShellElements; + if (!domain.shellElements().empty()) { + if (domain.shellElements().size() > + (std::numeric_limits::max)() / + kShellLocationCount) { + return recoveryFailure( + "invalid-recovery-dimensions", + {domain.sourcePath(), 0U}, + domain.sourceContentIdentity(), + "The shell result-row inventory exceeds the addressable range."); + } + std::vector>> directorsByNode( + domain.nodes().size()); + for (const auto& frame : domain.shellNodeInitialFrames()) { + if (frame.nodeIndex >= directorsByNode.size() || + directorsByNode[frame.nodeIndex].has_value()) { + return recoveryFailure( + "invalid-recovery-entity", + {domain.sourcePath(), 0U}, + std::to_string(frame.nodeIndex), + "Shell initial directors must map uniquely to model nodes."); + } + directorsByNode[frame.nodeIndex] = frame.director; + } + + shellCandidate.rows.reserve( + domain.shellElements().size() * kShellLocationCount); + expectedShellElements.reserve(domain.shellElements().size()); + constexpr std::array + locations{ + ShellMidsurfaceLocation::gp1, + ShellMidsurfaceLocation::gp2, + ShellMidsurfaceLocation::gp3, + ShellMidsurfaceLocation::gp4}; + constexpr std::array positions{ + ShellSectionPosition::bottom, + ShellSectionPosition::middle, + ShellSectionPosition::top}; + constexpr std::array zeta{-1.0, 0.0, 1.0}; + for (std::size_t elementOrder = 0U; + elementOrder < domain.shellElements().size(); + ++elementOrder) { + const EntityIndex elementIndex = + static_cast(elementOrder); + const auto& definition = domain.shellElements()[elementOrder]; + std::array nodes{}; + std::array, 4> directors{}; + for (std::size_t nodePosition = 0U; + nodePosition < definition.nodeIndices.size(); + ++nodePosition) { + const EntityIndex node = definition.nodeIndices[nodePosition]; + if (!directorsByNode[node].has_value()) { + return recoveryFailure( + "invalid-recovery-entity", + definition.location, + definition.sourceId.sourceLabelText, + "Shell recovery requires one initial director per element node."); + } + nodes[nodePosition] = &domain.nodes()[node]; + directors[nodePosition] = *directorsByNode[node]; + } + + auto shell = Mitc4Shell::create( + nodes, + directors, + domain.shellSections()[definition.sectionIndex], + domain.materials()[definition.materialIndex]); + if (!shell.hasValue()) { + return shell.status(); + } + Vector elementDisplacement{kShellElementDofCount}; + const auto& scatter = dofs.shellElementScatter(elementIndex); + for (std::size_t localDof = 0U; + localDof < kShellElementDofCount; + ++localDof) { + elementDisplacement[localDof] = + state.displacement()[scatter[localDof]]; + } + auto recovered = shell.value().recoverPhysical( + elementDisplacement); + if (!recovered.hasValue()) { + return recovered.status(); + } + const double accumulatedEnergy = + shellCandidate.physicalStrainEnergy + + recovered.value().strainEnergy; + if (!std::isfinite(accumulatedEnergy)) { + return recoveryFailure( + "nonfinite-recovery-value", + definition.location, + definition.sourceId.sourceLabelText, + "Source-order physical shell energy reduction must remain finite."); + } + shellCandidate.physicalStrainEnergy = accumulatedEnergy; + expectedShellElements.push_back(elementIndex); + + for (std::size_t point = 0U; + point < recovered.value().points.size(); + ++point) { + const auto& physicalPoint = recovered.value().points[point]; + ShellResultRow row{}; + row.element = elementIndex; + row.location = locations[point]; + row.naturalCoordinates = physicalPoint.naturalCoordinates; + row.localFrame = { + physicalPoint.localFrame.e1, + physicalPoint.localFrame.e2, + physicalPoint.localFrame.e3}; + row.generalizedStrain = physicalPoint.generalizedStrain; + row.sectionResultant = physicalPoint.sectionResultant; + for (std::size_t position = 0U; + position < positions.size(); + ++position) { + row.stress[position] = { + positions[position], + zeta[position], + physicalPoint.inPlaneStress[position]}; + } + shellCandidate.rows.push_back(std::move(row)); + } + } + const Status evidenceStatus = populateShellGlobalEvidence( + domain, + dofs, + state.externalForce(), + residual, + normalizedResidual, + shellCandidate); + if (!evidenceStatus.isOk()) { + return evidenceStatus; + } + } + + // Build and validate a complete candidate state first. This preserves the + // prior full residual, beam rows, and shell rows if any later shell or + // candidate-inventory validation fails. + AnalysisState candidateState = state; + candidateState.internalForce() = std::move(internalForce); + candidateState.residual() = std::move(residual); + candidateState.reaction() = std::move(reaction); + candidateState.endpointResults() = std::move(endpointRows); + candidateState.gaussResults() = std::move(gaussRows); + candidateState.stressResults() = std::move(stressRows); + const Status shellCommitStatus = candidateState.commitShellResults( + expectedShellElements, std::move(shellCandidate)); + if (!shellCommitStatus.isOk()) { + return shellCommitStatus; + } + state = std::move(candidateState); return Status::ok(); } diff --git a/tests/unit/results/result_recovery_test.cpp b/tests/unit/results/result_recovery_test.cpp index c7ced70..5d7581f 100644 --- a/tests/unit/results/result_recovery_test.cpp +++ b/tests/unit/results/result_recovery_test.cpp @@ -2,6 +2,7 @@ #include "fesa/analysis/analysis_model.hpp" #include "fesa/analysis/analysis_state.hpp" +#include "fesa/assembly/load_assembler.hpp" #include "fesa/assembly/parallel_for.hpp" #include "fesa/assembly/sparse_assembler.hpp" #include "fesa/fem/dof_manager.hpp" @@ -34,6 +35,13 @@ struct RecoveryFixture { std::unique_ptr stiffness; }; +struct ShellRecoveryFixture { + std::unique_ptr domain; + std::unique_ptr model; + std::unique_ptr dofs; + std::unique_ptr stiffness; +}; + fesa::ModelDefinition makeDefinition( const bool twoElements = false, std::vector> sectionPoints = {}, @@ -144,6 +152,148 @@ RecoveryFixture makeFixture( std::move(stiffness)}; } +fesa::ModelDefinition makeShellDefinition( + const bool twoElements = false, + const bool constrainAll = false, + std::vector loads = {}) { + const std::filesystem::path source{"models/shell-result-recovery.inp"}; + fesa::ModelDefinition definition{}; + definition.sourcePath = source; + definition.sourceContentIdentity = "fnv1a64:fedcba9876543210"; + definition.nodes = { + {{"Shell-1", 1, "1"}, {-1.0, -1.0, 0.0}, {source, 10U}}, + {{"Shell-1", 2, "2"}, {1.0, -1.0, 0.0}, {source, 11U}}, + {{"Shell-1", 3, "3"}, {3.0, -1.0, 0.0}, {source, 12U}}, + {{"Shell-1", 4, "4"}, {-1.0, 1.0, 0.0}, {source, 13U}}, + {{"Shell-1", 5, "5"}, {1.0, 1.0, 0.0}, {source, 14U}}, + {{"Shell-1", 6, "6"}, {3.0, 1.0, 0.0}, {source, 15U}}}; + if (!twoElements) { + definition.nodes.erase( + definition.nodes.begin() + 2, + definition.nodes.begin() + 3); + definition.nodes.erase(definition.nodes.begin() + 4); + } + definition.materials = { + {"Material", 120.0, 0.25, {source, 20U}}}; + definition.shellSections = { + {"ShellSection", 2.0, 0U, {source, 30U}}}; + definition.shellElements = {{ + {"Shell-1", 10, "10"}, + fesa::ShellSourceElementType::s4, + {0U, 1U, twoElements ? 4U : 3U, twoElements ? 3U : 2U}, + 0U, + 0U, + {source, 40U}}}; + if (twoElements) { + definition.shellElements.push_back({ + {"Shell-1", 20, "20"}, + fesa::ShellSourceElementType::s4r, + {1U, 2U, 5U, 4U}, + 0U, + 0U, + {source, 41U}}); + } + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + definition.shellNodeInitialFrames.push_back({ + static_cast(node), + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}}); + } + if (constrainAll) { + std::vector allNodes; + allNodes.reserve(definition.nodes.size()); + for (std::size_t node = 0U; node < definition.nodes.size(); ++node) { + allNodes.push_back(static_cast(node)); + } + definition.nodeSets.push_back( + {"All", {}, std::move(allNodes), {source, 50U}}); + } + definition.steps = {{ + "Step-1", + constrainAll + ? std::vector{ + {"All", 1, 6, 0.0, {source, 60U}}} + : std::vector{}, + std::move(loads), + 0.1, + 1.0, + 0.01, + 1.0, + {source, 59U}}}; + return definition; +} + +ShellRecoveryFixture makeShellFixture(fesa::ModelDefinition definition) { + auto domainResult = fesa::Domain::create(std::move(definition)); + if (!domainResult.hasValue()) { + throw std::runtime_error{ + "Shell recovery fixture Domain construction failed."}; + } + auto domain = std::make_unique( + std::move(domainResult.value())); + + auto modelResult = fesa::AnalysisModel::create(*domain); + if (!modelResult.hasValue()) { + throw std::runtime_error{ + "Shell recovery fixture AnalysisModel construction failed."}; + } + auto model = std::make_unique( + std::move(modelResult.value())); + + auto dofsResult = fesa::DofManager::create(*model); + if (!dofsResult.hasValue()) { + throw std::runtime_error{ + "Shell recovery fixture DofManager construction failed."}; + } + auto dofs = std::make_unique( + std::move(dofsResult.value())); + + fesa::SerialParallelFor serial; + auto stiffnessResult = fesa::SparseAssembler::assembleStiffness( + *model, *dofs, serial); + if (!stiffnessResult.hasValue()) { + throw std::runtime_error{ + "Shell recovery fixture stiffness assembly failed."}; + } + auto stiffness = std::make_unique( + std::move(stiffnessResult.value())); + return { + std::move(domain), + std::move(model), + std::move(dofs), + std::move(stiffness)}; +} + +fesa::AnalysisState makeShellPhysicalState( + const ShellRecoveryFixture& fixture) { + constexpr std::array generalized{ + 0.1, -0.05, 0.2, 0.3, -0.15, 0.25, 0.4, -0.3}; + auto state = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; + node < fixture.domain->nodes().size(); + ++node) { + const double x = fixture.domain->nodes()[node].coordinates[0U]; + const double y = fixture.domain->nodes()[node].coordinates[1U]; + const std::size_t offset = 6U * node; + state.displacement()[offset] = + generalized[0U] * x + 0.5 * generalized[2U] * y; + state.displacement()[offset + 1U] = + generalized[1U] * y + 0.5 * generalized[2U] * x; + state.displacement()[offset + 2U] = + generalized[6U] * x + generalized[7U] * y - + 0.5 * generalized[5U] * x * y; + state.displacement()[offset + 3U] = + -generalized[4U] * y - 0.5 * generalized[5U] * x; + state.displacement()[offset + 4U] = + generalized[3U] * x + 0.5 * generalized[5U] * y; + } + state.externalForce() = + fixture.stiffness->multiply(state.displacement()); + return state; +} + fesa::AnalysisState makeAxialEquilibriumState(const RecoveryFixture& fixture) { auto state = fesa::AnalysisState::create( *fixture.dofs, {"Step-1", 0U}); @@ -217,6 +367,11 @@ std::vector makeStationRows( TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) { const auto fixture = makeFixture(); auto state = makeAxialEquilibriumState(fixture); + fesa::ShellStateCandidate staleShellEvidence{}; + staleShellEvidence.physicalStrainEnergy = 123.0; + staleShellEvidence.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + staleShellEvidence.verificationMetrics = {1.0e-11, 2.0e-11, 3.0e-11}; + ASSERT_TRUE(state.commitShellResults({}, staleShellEvidence).isOk()); const fesa::Status status = fesa::ResultRecovery::recover( *fixture.model, *fixture.dofs, *fixture.stiffness, state); @@ -230,6 +385,14 @@ TEST(ResultRecovery, ComputesResidualReactionForNonzeroPrescription) { for (const std::size_t fullDof : fixture.dofs->freeDofs()) { EXPECT_DOUBLE_EQ(state.reaction()[fullDof], 0.0); } + EXPECT_TRUE(state.shellResults().empty()); + EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), 0.0); + EXPECT_EQ( + state.equilibrium(), + (std::array{0.0, 0.0, 0.0, 0.0, 0.0, 0.0})); + EXPECT_EQ( + state.verificationMetrics(), + (std::array{0.0, 0.0, 0.0})); } TEST(ResultRecovery, EnforcesNormalizedFreeResidual) { @@ -469,3 +632,359 @@ TEST(ResultRecovery, RequiresInteriorEndpointConsistencyWithoutAveraging) { ASSERT_FALSE(jumped.hasValue()); expectStatusCode(jumped.status(), "ineligible-node-station"); } + +// MITC4-REC-001 +TEST(ResultRecovery, RecoversShellRowsInStableElementAndGpOrder) { + const auto fixture = makeShellFixture(makeShellDefinition(true)); + auto state = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; + node < fixture.domain->nodes().size(); + ++node) { + const double x = fixture.domain->nodes()[node].coordinates[0U]; + const double y = fixture.domain->nodes()[node].coordinates[1U]; + state.displacement()[node * 6U] = 0.1 * x + 0.1 * y; + state.displacement()[node * 6U + 1U] = -0.05 * y + 0.1 * x; + } + state.externalForce() = + fixture.stiffness->multiply(state.displacement()); + + const auto status = fesa::ResultRecovery::recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state); + + ASSERT_TRUE(status.isOk()); + ASSERT_EQ(state.shellResults().size(), 8U); + const double gauss = 1.0 / std::sqrt(3.0); + const std::array locations{ + fesa::ShellMidsurfaceLocation::gp1, + fesa::ShellMidsurfaceLocation::gp2, + fesa::ShellMidsurfaceLocation::gp3, + fesa::ShellMidsurfaceLocation::gp4}; + const std::array, 4> coordinates{ + std::array{-gauss, -gauss}, + std::array{gauss, -gauss}, + std::array{gauss, gauss}, + std::array{-gauss, gauss}}; + constexpr std::array expectedStrain{ + 0.1, -0.05, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0}; + constexpr std::array expectedResultant{ + 22.4, -6.4, 19.2, 0.0, 0.0, 0.0, 0.0, 0.0}; + const std::array, 3> expectedFrame{{ + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}}}; + for (std::size_t element = 0U; element < 2U; ++element) { + for (std::size_t point = 0U; point < locations.size(); ++point) { + const auto& row = state.shellResults()[element * 4U + point]; + EXPECT_EQ(row.element, element); + EXPECT_EQ(row.location, locations[point]); + EXPECT_EQ(row.naturalCoordinates, coordinates[point]); + EXPECT_EQ(row.localFrame, expectedFrame); + for (std::size_t component = 0U; + component < expectedStrain.size(); + ++component) { + EXPECT_NEAR( + row.generalizedStrain[component], + expectedStrain[component], + 1.0e-12); + EXPECT_NEAR( + row.sectionResultant[component], + expectedResultant[component], + 1.0e-12); + } + } + } +} + +// MITC4-REC-002 +TEST(ResultRecovery, RecoversDirectBottomMiddleTopShellStress) { + const auto fixture = makeShellFixture(makeShellDefinition()); + auto state = makeShellPhysicalState(fixture); + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + state) + .isOk()); + constexpr std::array positions{ + fesa::ShellSectionPosition::bottom, + fesa::ShellSectionPosition::middle, + fesa::ShellSectionPosition::top}; + constexpr std::array zeta{-1.0, 0.0, 1.0}; + constexpr std::array, 3> expectedStress{{ + {-22.4, 6.4, -2.4}, + {11.2, -3.2, 9.6}, + {44.8, -12.8, 21.6}}}; + ASSERT_EQ(state.shellResults().size(), 4U); + for (const auto& row : state.shellResults()) { + for (std::size_t position = 0U; + position < positions.size(); + ++position) { + EXPECT_EQ(row.stress[position].position, positions[position]); + EXPECT_DOUBLE_EQ(row.stress[position].zeta, zeta[position]); + for (std::size_t component = 0U; + component < expectedStress[position].size(); + ++component) { + EXPECT_NEAR( + row.stress[position].components[component], + expectedStress[position][component], + 1.0e-12); + } + } + } +} + +// MITC4-REC-003 +TEST(ResultRecovery, SumsOnlyPhysicalShellEnergyInSourceOrder) { + const auto fixture = makeShellFixture(makeShellDefinition()); + auto state = makeShellPhysicalState(fixture); + constexpr std::array drill{1.0, -1.0, 1.0, -1.0}; + for (std::size_t node = 0U; node < drill.size(); ++node) { + state.displacement()[node * 6U + 5U] = drill[node]; + } + state.externalForce() = + fixture.stiffness->multiply(state.displacement()); + const double stabilizedEnergy = 0.5 * state.displacement().dot( + fixture.stiffness->multiply(state.displacement())); + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + state) + .isOk()); + EXPECT_NEAR(state.physicalStrainEnergy(), 72.16, 1.0e-12); + EXPECT_GT(stabilizedEnergy, state.physicalStrainEnergy()); +} + +// MITC4-REC-004 +TEST(ResultRecovery, KeepsFullResidualAndComputesGlobalShellEquilibrium) { + const std::filesystem::path source{"models/shell-result-recovery.inp"}; + const std::vector loads{ + {"1", 1, 5.0, {source, 70U}}, + {"5", 1, -5.0, {source, 71U}}, + {"1", 4, 2.0, {source, 72U}}, + {"5", 4, -2.0, {source, 73U}}}; + const auto fixture = makeShellFixture( + makeShellDefinition(false, true, loads)); + auto state = fesa::AnalysisState::create( + *fixture.dofs, {"Step-1", 0U}); + auto fullLoad = fesa::LoadAssembler::assembleFullNodalLoad( + *fixture.model, *fixture.dofs); + ASSERT_TRUE(fullLoad.hasValue()); + state.externalForce() = std::move(fullLoad.value()); + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + state) + .isOk()); + ASSERT_EQ(state.shellResults().size(), 4U); + for (std::size_t fullDof = 0U; + fullDof < fixture.dofs->fullDofCount(); + ++fullDof) { + EXPECT_DOUBLE_EQ(state.internalForce()[fullDof], 0.0); + EXPECT_DOUBLE_EQ( + state.residual()[fullDof], -state.externalForce()[fullDof]); + EXPECT_DOUBLE_EQ(state.reaction()[fullDof], state.residual()[fullDof]); + } + EXPECT_EQ( + state.equilibrium(), + (std::array{0.0, 0.0, 0.0, 0.0, 0.0, 0.0})); + EXPECT_EQ( + state.verificationMetrics(), + (std::array{0.0, 0.0, 0.0})); + + const auto freeFixture = makeShellFixture(makeShellDefinition()); + auto perturbed = makeShellPhysicalState(freeFixture); + perturbed.externalForce()[0U] += 1.0e-9; + ASSERT_TRUE(fesa::ResultRecovery::recover( + *freeFixture.model, + *freeFixture.dofs, + *freeFixture.stiffness, + perturbed) + .isOk()); + for (const double metric : perturbed.verificationMetrics()) { + EXPECT_GT(metric, 0.0); + EXPECT_LE(metric, 1.0e-10); + } +} + +// MITC4-REC-004 +TEST(ResultRecovery, UsesGlobalOriginForShellMomentBalance) { + auto centeredDefinition = makeShellDefinition(); + auto translatedDefinition = centeredDefinition; + constexpr std::array translation{7.0, 11.0, 0.0}; + for (auto& node : translatedDefinition.nodes) { + for (std::size_t component = 0U; + component < translation.size(); + ++component) { + node.coordinates[component] += translation[component]; + } + } + const auto centeredFixture = makeShellFixture( + std::move(centeredDefinition)); + const auto translatedFixture = makeShellFixture( + std::move(translatedDefinition)); + auto centered = makeShellPhysicalState(centeredFixture); + auto translated = makeShellPhysicalState(translatedFixture); + centered.externalForce()[0U] += 1.0e-9; + translated.externalForce()[0U] += 1.0e-9; + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *centeredFixture.model, + *centeredFixture.dofs, + *centeredFixture.stiffness, + centered) + .isOk()); + ASSERT_TRUE(fesa::ResultRecovery::recover( + *translatedFixture.model, + *translatedFixture.dofs, + *translatedFixture.stiffness, + translated) + .isOk()); + std::array centeredForce{}; + for (std::size_t component = 0U; component < 3U; ++component) { + centeredForce[component] = centered.equilibrium()[component]; + EXPECT_NEAR( + translated.equilibrium()[component], + centeredForce[component], + 1.0e-12); + } + const std::array translatedMomentDelta{ + translation[1U] * centeredForce[2U] - + translation[2U] * centeredForce[1U], + translation[2U] * centeredForce[0U] - + translation[0U] * centeredForce[2U], + translation[0U] * centeredForce[1U] - + translation[1U] * centeredForce[0U]}; + for (std::size_t component = 0U; component < 3U; ++component) { + EXPECT_NEAR( + translated.equilibrium()[3U + component] - + centered.equilibrium()[3U + component], + translatedMomentDelta[component], + 5.0e-12); + } + EXPECT_GT(std::abs(translatedMomentDelta[2U]), 1.0e-9); +} + +// MITC4-REC-004 +TEST(ResultRecovery, UsesScaleAwareShellMetricsAndRejectsExcess) { + const auto fixture = makeShellFixture(makeShellDefinition()); + auto subunit = makeShellPhysicalState(fixture); + auto large = makeShellPhysicalState(fixture); + constexpr double subunitScale = 1.0e-6; + constexpr double largeScale = 1.0e6; + subunit.displacement().scale(subunitScale); + subunit.externalForce().scale(subunitScale); + subunit.externalForce()[0U] += 1.0e-9 * subunitScale; + large.displacement().scale(largeScale); + large.externalForce().scale(largeScale); + large.externalForce()[0U] += 1.0e-9 * largeScale; + + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + subunit) + .isOk()); + ASSERT_TRUE(fesa::ResultRecovery::recover( + *fixture.model, + *fixture.dofs, + *fixture.stiffness, + large) + .isOk()); + for (std::size_t metric = 0U; metric < 3U; ++metric) { + EXPECT_GT(subunit.verificationMetrics()[metric], 0.0); + EXPECT_GT(large.verificationMetrics()[metric], 0.0); + EXPECT_NEAR( + subunit.verificationMetrics()[metric], + large.verificationMetrics()[metric], + 1.0e-13); + } + + auto constrainedDefinition = makeShellDefinition(false, true); + constrainedDefinition.steps[0U].boundaries[0U].value = 1.0; + const auto constrainedFixture = makeShellFixture( + std::move(constrainedDefinition)); + auto rejected = fesa::AnalysisState::create( + *constrainedFixture.dofs, {"Step-1", 0U}); + for (std::size_t fullDof = 0U; + fullDof < constrainedFixture.dofs->fullDofCount(); + ++fullDof) { + rejected.displacement()[fullDof] = 1.0; + } + const std::vector unbalancedEntry{ + {0U, 0U, 1.0, 0U, 0U}}; + auto unbalancedStiffness = fesa::SparseMatrix::fromCoo( + constrainedFixture.dofs->fullDofCount(), + constrainedFixture.dofs->fullDofCount(), + unbalancedEntry, + constrainedFixture.dofs->sparsePattern()); + ASSERT_TRUE(unbalancedStiffness.hasValue()); + expectStatusCode( + fesa::ResultRecovery::recover( + *constrainedFixture.model, + *constrainedFixture.dofs, + unbalancedStiffness.value(), + rejected), + "global-equilibrium-tolerance-failure"); +} + +// MITC4-REC-005 +TEST(ResultRecovery, InvalidLaterShellLeavesEntirePriorStateUnchanged) { + const auto validFixture = makeShellFixture(makeShellDefinition(true)); + auto state = makeShellPhysicalState(validFixture); + ASSERT_TRUE(fesa::ResultRecovery::recover( + *validFixture.model, + *validFixture.dofs, + *validFixture.stiffness, + state) + .isOk()); + ASSERT_EQ(state.shellResults().size(), 8U); + const auto priorFirstRow = state.shellResults().front(); + const double priorEnergy = state.physicalStrainEnergy(); + const auto priorEquilibrium = state.equilibrium(); + const auto priorMetrics = state.verificationMetrics(); + state.internalForce()[0U] = 91.0; + state.residual()[0U] = 92.0; + state.reaction()[0U] = 93.0; + state.endpointResults().push_back({}); + state.displacement()[2U * 6U] = + (std::numeric_limits::max)(); + state.displacement()[5U * 6U] = + (std::numeric_limits::max)(); + state.externalForce() = fesa::Vector{validFixture.dofs->fullDofCount()}; + auto zeroStiffness = fesa::SparseMatrix::fromCoo( + validFixture.dofs->fullDofCount(), + validFixture.dofs->fullDofCount(), + {}, + validFixture.dofs->sparsePattern()); + ASSERT_TRUE(zeroStiffness.hasValue()); + + const auto status = fesa::ResultRecovery::recover( + *validFixture.model, + *validFixture.dofs, + zeroStiffness.value(), + state); + + expectStatusCode(status, "invalid-shell-recovery"); + EXPECT_DOUBLE_EQ(state.internalForce()[0U], 91.0); + EXPECT_DOUBLE_EQ(state.residual()[0U], 92.0); + EXPECT_DOUBLE_EQ(state.reaction()[0U], 93.0); + EXPECT_EQ(state.endpointResults().size(), 1U); + ASSERT_EQ(state.shellResults().size(), 8U); + EXPECT_EQ(state.shellResults().front().element, priorFirstRow.element); + EXPECT_EQ(state.shellResults().front().location, priorFirstRow.location); + EXPECT_EQ( + state.shellResults().front().generalizedStrain, + priorFirstRow.generalizedStrain); + EXPECT_EQ( + state.shellResults().front().sectionResultant, + priorFirstRow.sectionResultant); + EXPECT_DOUBLE_EQ(state.physicalStrainEnergy(), priorEnergy); + EXPECT_EQ(state.equilibrium(), priorEquilibrium); + EXPECT_EQ(state.verificationMetrics(), priorMetrics); +}