From c6188b3ad71e30e4483455d90ce1711a44332a20 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Tue, 20 Oct 2026 09:21:39 +0900 Subject: [PATCH 1/2] Connect MITC4 drilling evidence across analysis procedures --- include/fesa/analysis/analysis_state.h | 4 + .../fesa/analysis/dynamic_analysis_state.h | 2 + include/fesa/analysis/eigen_analysis_state.h | 2 + .../fesa/analysis/nonlinear_analysis_state.h | 3 + include/fesa/elements/element.h | 5 + include/fesa/elements/element_capabilities.h | 3 + include/fesa/results/result_records.h | 12 + src/fesa/analysis/analysis_state.cpp | 16 ++ src/fesa/analysis/dynamic_analysis_state.cpp | 8 + src/fesa/analysis/hht_dynamic_algebra.cpp | 15 +- src/fesa/analysis/hht_dynamic_assembly.cpp | 44 +++- src/fesa/analysis/hht_dynamic_internal.h | 3 + src/fesa/analysis/hht_reduced_operators.cpp | 3 +- .../analysis/linear_frequency_procedure.cpp | 7 + .../analysis/nonlinear_analysis_state.cpp | 10 +- .../analysis/nonlinear_static_procedure.cpp | 43 +++- src/fesa/elements/element_capabilities.cpp | 4 + src/fesa/elements/mitc4_shell.cpp | 16 ++ src/fesa/results/analysis_state_commit.cpp | 3 +- .../results/global_equilibrium_recovery.cpp | 16 +- src/fesa/results/mitc4_drilling_evidence.h | 54 +++++ src/fesa/results/shell_result_recovery.cpp | 23 ++ .../linear_frequency_procedure_test.cpp | 40 +++- .../nonlinear_static_procedure_test.cpp | 61 +++++ .../analysis/dynamic_analysis_state_test.cpp | 24 ++ .../analysis/hht_dynamic_procedure_test.cpp | 10 + .../unit/analysis/hht_reduced_shell_test.cpp | 222 +++++++++++++++++- .../unit/analysis/hht_shell_physics_test.cpp | 14 ++ .../nonlinear_analysis_state_test.cpp | 25 ++ .../result_recovery_components_test.cpp | 99 ++++++++ 30 files changed, 763 insertions(+), 28 deletions(-) create mode 100644 src/fesa/results/mitc4_drilling_evidence.h diff --git a/include/fesa/analysis/analysis_state.h b/include/fesa/analysis/analysis_state.h index 0222ae2..e1baff6 100644 --- a/include/fesa/analysis/analysis_state.h +++ b/include/fesa/analysis/analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/status.h" @@ -72,6 +73,8 @@ class AnalysisState { double PhysicalStrainEnergy() const noexcept; /// @brief Returns global force and moment equilibrium components. const std::array& Equilibrium() const noexcept; + /// @brief Returns actual MITC4 mode and artificial spatial moment evidence. + const std::optional& Mitc4Drilling() const noexcept; /// @brief Returns normalized shell verification metrics. const std::array& VerificationMetrics() const noexcept; @@ -96,6 +99,7 @@ class AnalysisState { double physical_strain_energy_{0.0}; std::array equilibrium_{}; std::array verification_metrics_{}; + std::optional mitc4_drilling_; }; } // namespace fesa diff --git a/include/fesa/analysis/dynamic_analysis_state.h b/include/fesa/analysis/dynamic_analysis_state.h index 96bca40..d96f83c 100644 --- a/include/fesa/analysis/dynamic_analysis_state.h +++ b/include/fesa/analysis/dynamic_analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/source_identity.h" @@ -54,6 +55,7 @@ struct DynamicAcceptedFrameState { double external_work = 0.0; double energy_defect = 0.0; std::array global_equilibrium{}; + std::optional mitc4_drilling; std::array linear_momentum{}; std::array angular_momentum{}; std::array linear_momentum_rate{}; diff --git a/include/fesa/analysis/eigen_analysis_state.h b/include/fesa/analysis/eigen_analysis_state.h index d5aecd0..cbe376c 100644 --- a/include/fesa/analysis/eigen_analysis_state.h +++ b/include/fesa/analysis/eigen_analysis_state.h @@ -8,6 +8,7 @@ #include #include "fesa/math/matrix.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/solvers/eigen/eigen_solver.h" namespace fesa { @@ -53,6 +54,7 @@ struct EigenAnalysisState { explicit EigenAnalysisState(FeastEigenResult backend_result) : backend{std::move(backend_result)} {} + std::optional drilling_mode; std::size_t full_dof_count{0U}; std::size_t free_dof_count{0U}; std::size_t mass_null_dimension{0U}; diff --git a/include/fesa/analysis/nonlinear_analysis_state.h b/include/fesa/analysis/nonlinear_analysis_state.h index 7efee97..c0e3127 100644 --- a/include/fesa/analysis/nonlinear_analysis_state.h +++ b/include/fesa/analysis/nonlinear_analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/source_identity.h" @@ -46,6 +47,7 @@ struct NonlinearTrialState { std::vector shell_results; double physical_strain_energy; std::array global_equilibrium; + std::optional mitc4_drilling; }; /// @brief Stores one authoritative nonlinear increment frame. @@ -69,6 +71,7 @@ struct AcceptedFrameState { std::vector shell_results; double physical_strain_energy; std::array global_equilibrium; + std::optional mitc4_drilling; }; /// @brief Owns committed nonlinear state and validated accepted frames. diff --git a/include/fesa/elements/element.h b/include/fesa/elements/element.h index 5430172..6a1ec5d 100644 --- a/include/fesa/elements/element.h +++ b/include/fesa/elements/element.h @@ -1,11 +1,14 @@ #ifndef FESA_ELEMENTS_ELEMENT_H_ #define FESA_ELEMENTS_ELEMENT_H_ +#include +#include #include #include #include "fesa/core/status.h" #include "fesa/elements/element_dof_layout.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/math/matrix.h" #include "fesa/math/vector.h" #include "fesa/results/result_records.h" @@ -30,6 +33,8 @@ struct BeamElementResultRows { struct ShellElementResultRows { std::vector rows; double physical_strain_energy{0.0}; + std::optional drilling_mode; + std::array drilling_moment{}; }; /// @brief Selects the physical recovery shape of one runtime element. diff --git a/include/fesa/elements/element_capabilities.h b/include/fesa/elements/element_capabilities.h index d6c610d..c68ed6a 100644 --- a/include/fesa/elements/element_capabilities.h +++ b/include/fesa/elements/element_capabilities.h @@ -178,6 +178,8 @@ class ShellReferenceElement : public Element { virtual Result Stiffness() const = 0; /// @copydoc Mitc4Shell::Mass virtual Result Mass(double density) const = 0; + /// @brief Returns the actual immutable reference-kernel drilling formulation. + virtual Mitc4DrillingMode DrillingMode() const noexcept = 0; }; /// @brief Borrows the same linear MITC4 kernel for stiffness and recovery. @@ -190,6 +192,7 @@ class Mitc4ReferenceAdapter final : public ShellReferenceElement { Result Recover(const Vector& u) const override; Result Stiffness() const override; Result Mass(double density) const override; + Mitc4DrillingMode DrillingMode() const noexcept override; private: const Mitc4Shell& kernel_; diff --git a/include/fesa/results/result_records.h b/include/fesa/results/result_records.h index c1a5848..b867cab 100644 --- a/include/fesa/results/result_records.h +++ b/include/fesa/results/result_records.h @@ -4,9 +4,11 @@ #include #include #include +#include #include #include "fesa/core/source_identity.h" +#include "fesa/elements/mitc4_drilling_mode.h" namespace fesa { @@ -70,6 +72,15 @@ struct ShellResultRow { std::array stress; }; +/// @brief Separates artificial internal spatial moment from raw equilibrium. +/// @note The actual immutable kernel mode accompanies every production frame. +struct Mitc4DrillingEvidence { + Mitc4DrillingMode mode; + std::array drilling_moment{}; + // Raw force balance and raw moment balance minus positive internal drilling. + std::array equilibrium_with_drilling{}; +}; + /// @brief Carries a complete shell result candidate for atomic validation. struct ShellStateCandidate { std::vector rows; @@ -79,6 +90,7 @@ struct ShellStateCandidate { // [FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED, // MOMENT_BALANCE_NORMALIZED]. std::array verification_metrics{}; + std::optional mitc4_drilling; }; } // namespace fesa diff --git a/src/fesa/analysis/analysis_state.cpp b/src/fesa/analysis/analysis_state.cpp index 786be82..dabb8a5 100644 --- a/src/fesa/analysis/analysis_state.cpp +++ b/src/fesa/analysis/analysis_state.cpp @@ -1,5 +1,7 @@ #include "fesa/analysis/analysis_state.h" +#include "results/mitc4_drilling_evidence.h" + #include #include #include @@ -167,9 +169,18 @@ Status AnalysisState::CommitShellResults( "Shell energy, equilibrium, and normalized metrics must be finite."); } + if (!results_internal::ValidDrillingEvidence(candidate.mitc4_drilling, + candidate.equilibrium) || + (mitc4_drilling_ && !results_internal::SameDrillingInventory( + mitc4_drilling_, candidate.mitc4_drilling))) + return ShellCandidateFailure("invalid-shell-drilling-evidence", + identity_.step_name, + "Drilling evidence must be finite, consistent, and retain its mode."); + shell_results_ = std::move(candidate.rows); physical_strain_energy_ = candidate.physical_strain_energy; equilibrium_ = candidate.equilibrium; + mitc4_drilling_ = candidate.mitc4_drilling; verification_metrics_ = candidate.verification_metrics; return Status::Ok(); } @@ -187,6 +198,11 @@ const std::array& AnalysisState::Equilibrium() const noexcept { return equilibrium_; } +const std::optional& AnalysisState::Mitc4Drilling() + const noexcept { + return mitc4_drilling_; +} + const std::array& AnalysisState::VerificationMetrics() const noexcept { return verification_metrics_; diff --git a/src/fesa/analysis/dynamic_analysis_state.cpp b/src/fesa/analysis/dynamic_analysis_state.cpp index 3e374c0..06bb1f7 100644 --- a/src/fesa/analysis/dynamic_analysis_state.cpp +++ b/src/fesa/analysis/dynamic_analysis_state.cpp @@ -1,5 +1,7 @@ #include "fesa/analysis/dynamic_analysis_state.h" +#include "results/mitc4_drilling_evidence.h" + #include #include #include @@ -273,6 +275,12 @@ Status DynamicAnalysisState::CommitAcceptedFrame( return Failure(FailureCategory::kSolver, "Frame physical evidence is nonfinite or has invalid size."); } + if (!results_internal::ValidDrillingEvidence(frame.mitc4_drilling, + frame.global_equilibrium) || + (!accepted_frames_.empty() && !results_internal::SameDrillingInventory( + frame.mitc4_drilling, accepted_frames_.back().mitc4_drilling))) + return Failure(FailureCategory::kSolver, + "Accepted drilling evidence must be finite, consistent, and retain its mode."); for (const std::size_t dof : constrained_dofs_) { if (frame.generalized_coordinates[dof] != 0.0 || frame.velocity[dof] != 0.0 || frame.acceleration[dof] != 0.0 || diff --git a/src/fesa/analysis/hht_dynamic_algebra.cpp b/src/fesa/analysis/hht_dynamic_algebra.cpp index 09d6803..7b4c8cd 100644 --- a/src/fesa/analysis/hht_dynamic_algebra.cpp +++ b/src/fesa/analysis/hht_dynamic_algebra.cpp @@ -64,6 +64,16 @@ bool Within(const double value, const double scale, const double tolerance) { Status Validate(const Domain& domain, const DofManager& dofs, const DynamicAssembly& assembly, const double length) { const std::size_t count = dofs.FullDofCount(); + if (assembly.frame.mitc4_drilling && + assembly.frame.mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) { + if (assembly.hidden_constraint.Size() != assembly.rotations.size()) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must cover every node."); + for (std::size_t i = 0U; i < assembly.hidden_constraint.Size(); ++i) + if (!std::isfinite(assembly.hidden_constraint[i])) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must remain finite."); + } if (!std::isfinite(length) || !(length > 0.0) || count % 6U != 0U || assembly.rotations.size() != count / 6U) { return Failure(domain, "dynamic-invalid-algebra-state", @@ -547,6 +557,8 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const SparseMatrix& gradient, DynamicAcceptedFrameState* frame, const bool require_acceleration) { + const bool direct = frame->mitc4_drilling && + frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring; Vector psi{nodes.size()}; Vector psi_scale{nodes.size()}; Vector velocity_defect{nodes.size()}; @@ -575,9 +587,10 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const double g = frame->internal_force[full] - frame->external_force[full]; const double n = rotation.drill_null[c]; - psi[row] += n * g; + if (!direct) psi[row] += n * g; psi_scale[row] += std::abs(n) * frame->residual_scale[full]; } + if (direct) psi[row] = assembly.hidden_constraint[node]; velocity_defect[row] += velocity_correction; AddCompensated(assembly.hidden_second[node], &acceleration_defect[row], &acceleration_correction); diff --git a/src/fesa/analysis/hht_dynamic_assembly.cpp b/src/fesa/analysis/hht_dynamic_assembly.cpp index f4a5e1a..8faf890 100644 --- a/src/fesa/analysis/hht_dynamic_assembly.cpp +++ b/src/fesa/analysis/hht_dynamic_assembly.cpp @@ -18,6 +18,7 @@ #include "fesa/elements/finite_rotation.h" #include "fesa/loads/load.h" #include "fesa/model/source_target_resolver.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa::hht_internal { @@ -498,6 +499,10 @@ Result EvaluateAssembly( internal[e] = std::move(f.Value()); inertia[e] = std::move(m.Value()); }); + if (shells.empty()) + return Result::Failure(Failure("HHT requires a shell kernel.")); + frame.mitc4_drilling = + Mitc4DrillingEvidence{shells.front().internal_kernel->DrillingMode()}; frame.inertia_force = Vector(count); frame.internal_force = Vector(count); frame.residual = Vector(count); @@ -515,6 +520,7 @@ Result EvaluateAssembly( std::vector> constraint_gradients(domain.Nodes().size()); std::vector external_done(domain.Nodes().size(), false); + Vector hidden_constraint(directional_third ? domain.Nodes().size() : 0U); Vector third(count); Vector third_scale(count); std::array, 6> entries; @@ -532,9 +538,31 @@ Result EvaluateAssembly( return Result::Failure(*failures[e]); const auto& f = *internal[e]; const auto& m = *inertia[e]; + if (shells[e].internal_kernel->DrillingMode() != frame.mitc4_drilling->mode || + f.drilling_force.Size() != 24U) + return Result::Failure(Failure("Invalid actual drilling mode or force.")); + for (std::size_t node = 0U; node < 4U; ++node) { + const auto local = 6U * node; + const auto global = shells[e].scatter[local]; + const Vector3 force = Triple(f.drilling_force, local); + const auto spatial = SpatialMomentSO3( + Triple(frame.generalized_coordinates, global + 3U), + Triple(f.drilling_force, local + 3U)); + if (!spatial.HasValue()) + return Result::Failure(spatial.GetStatus()); + const Vector3 position = Vector3(domain.Nodes()[global / 6U].coordinates) + + Triple(frame.generalized_coordinates, global); + const Vector3 moment = position.Cross(force) + spatial.Value(); + for (std::size_t c = 0U; c < 3U; ++c) + frame.mitc4_drilling->drilling_moment[c] += moment[c]; + } + for (const double value : frame.mitc4_drilling->drilling_moment) + if (!std::isfinite(value)) + return Result::Failure(Failure("Drilling spatial moment overflowed.")); if (directional_third) { for (std::size_t local = 0U; local < 4U; ++local) { const auto node = shells[e].scatter[6U * local] / 6U; + hidden_constraint[node] += f.drilling_constraint[local]; auto& gradient = constraint_gradients[node]; auto& hessian = constraint_hessians[node]; for (std::size_t i = 0U; i < 24U; ++i) @@ -555,6 +583,7 @@ Result EvaluateAssembly( Triple(frame.external_force, first), {}, {}); if (!external.HasValue()) return Result::Failure(external.GetStatus()); + hidden_constraint[node] += external.Value().value; for (std::size_t i = 0U; i < 3U; ++i) { gradient[first + i] += external.Value().gradient[i]; for (std::size_t j = 0U; j < 3U; ++j) @@ -740,7 +769,8 @@ Result EvaluateAssembly( std::move(third), std::move(third_scale), std::move(rotations), momentum_scale, angular_scale, std::move(hidden_second), std::move(hidden_second_scale), std::move(hidden_gradient), - std::move(trial_constraints), std::move(field_flux)}); + std::move(trial_constraints), std::move(field_flux), + std::move(hidden_constraint)}); } Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs, @@ -776,10 +806,20 @@ Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs, frame->global_equilibrium[c] = force_error[c]; frame->global_equilibrium[3U + c] = moment_error[c]; } + Vector3 acceptance_moment = moment_error; + if (frame->mitc4_drilling) { + results_internal::CompleteDrillingEvidence(frame->global_equilibrium, + &*frame->mitc4_drilling); + if (!results_internal::ValidDrillingEvidence(frame->mitc4_drilling, + frame->global_equilibrium)) + return Failure("Invalid dynamic drilling equilibrium evidence."); + if (frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) + acceptance_moment = moment_error - Vector3(frame->mitc4_drilling->drilling_moment); + } if (!force_error.IsFinite() || !moment_error.IsFinite() || !std::isfinite(force_scale) || !std::isfinite(moment_scale) || Ratio(force_error.Norm(), force_scale) > 1e-8 || - Ratio(moment_error.Norm(), moment_scale) > 1e-8) + Ratio(acceptance_moment.Norm(), moment_scale) > 1e-8) return Failure("Dynamic physical force or moment balance exceeds 1e-8."); return Status::Ok(); } diff --git a/src/fesa/analysis/hht_dynamic_internal.h b/src/fesa/analysis/hht_dynamic_internal.h index e18dd62..5c5642c 100644 --- a/src/fesa/analysis/hht_dynamic_internal.h +++ b/src/fesa/analysis/hht_dynamic_internal.h @@ -97,6 +97,9 @@ struct DynamicAssembly { std::optional hidden_gradient; std::optional trial_constraints; std::array, 2> field_flux; + // Source-ordered scalar primitive, including the external term exactly once. + // Empty when directional primitives were not evaluated; never a zero fallback. + Vector hidden_constraint{0U}; }; /// @brief Builds checked shell evaluators in stable source order. diff --git a/src/fesa/analysis/hht_reduced_operators.cpp b/src/fesa/analysis/hht_reduced_operators.cpp index be4781f..e94d382 100644 --- a/src/fesa/analysis/hht_reduced_operators.cpp +++ b/src/fesa/analysis/hht_reduced_operators.cpp @@ -82,7 +82,8 @@ const hht_internal::DirectionalShellConstraintData* MatchingConstraints( nonlinear_mitc4_internal::ConfigurationSnapshot::Read(*shell.internal_kernel); if (record.configuration.source.file != configuration.source.file || record.configuration.source.line != configuration.source.line || - record.configuration.identity != configuration.identity) + record.configuration.identity != configuration.identity || + record.configuration.drilling_mode != configuration.drilling_mode) return nullptr; for (std::size_t i = 0U; i < configuration.values.size(); ++i) if (!std::isfinite(record.configuration.values[i]) || diff --git a/src/fesa/analysis/linear_frequency_procedure.cpp b/src/fesa/analysis/linear_frequency_procedure.cpp index d351098..3c9a0c9 100644 --- a/src/fesa/analysis/linear_frequency_procedure.cpp +++ b/src/fesa/analysis/linear_frequency_procedure.cpp @@ -217,6 +217,13 @@ Result ExecuteCandidate(const Domain& domain, if (!backend.HasValue()) return Result::Failure(backend.GetStatus()); EigenAnalysisState state(std::move(backend.Value())); + for (const auto& reference : references) { + if (state.drilling_mode && *state.drilling_mode != reference->DrillingMode()) + return Result::Failure(Failure( + "invalid-shell-drilling-evidence", "mode", + "All frequency kernels must retain one actual drilling mode.")); + state.drilling_mode = reference->DrillingMode(); + } state.full_dof_count = dofs.FullDofCount(); state.free_dof_count = reduction.Value().free_count; state.mass_null_dimension = reduction.Value().null_count; diff --git a/src/fesa/analysis/nonlinear_analysis_state.cpp b/src/fesa/analysis/nonlinear_analysis_state.cpp index 23c3ee4..d96f3f0 100644 --- a/src/fesa/analysis/nonlinear_analysis_state.cpp +++ b/src/fesa/analysis/nonlinear_analysis_state.cpp @@ -12,6 +12,7 @@ #include "fesa/elements/element_definition.h" #include "math/finite_values.h" #include "results/shell_result_recovery.h" +#include "results/mitc4_drilling_evidence.h" #include "analysis/shell_state_inventory.h" namespace fesa { @@ -110,6 +111,13 @@ Status ValidateAcceptedFrame(const AcceptedFrameState& frame, "Accepted frame physical result inventory is invalid."); } + if (!results_internal::ValidDrillingEvidence(frame.mitc4_drilling, + frame.global_equilibrium) || + (accepted_frame_count != 0U && !results_internal::SameDrillingInventory( + frame.mitc4_drilling, committed.mitc4_drilling))) + return StateFailure(FailureCategory::kSolver, "invalid-shell-drilling-evidence", {}, + "Accepted drilling evidence must be finite and retain its mode."); + const double gauss = 1.0 / std::sqrt(3.0); const std::array expected_locations{ShellMidsurfaceLocation::kGp1, @@ -243,7 +251,7 @@ Status NonlinearAnalysisState::CommitAcceptedFrame( frame.reaction, frame.shell_results, frame.physical_strain_energy, - frame.global_equilibrium}; + frame.global_equilibrium, frame.mitc4_drilling}; std::vector accepted_candidate = accepted_frames_; accepted_candidate.push_back(frame); diff --git a/src/fesa/analysis/nonlinear_static_procedure.cpp b/src/fesa/analysis/nonlinear_static_procedure.cpp index 4412f22..a47760e 100644 --- a/src/fesa/analysis/nonlinear_static_procedure.cpp +++ b/src/fesa/analysis/nonlinear_static_procedure.cpp @@ -33,6 +33,7 @@ #include "assembly/ordered_contribution_internal.h" #include "math/finite_values.h" #include "elements/shell_kernel_values.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa { namespace { @@ -57,6 +58,7 @@ struct InternalAssembly { SparseMatrix tangent; std::vector shell_rows; double physical_strain_energy; + Mitc4DrillingEvidence mitc4_drilling; }; Status ProcedureFailure(const FailureCategory category, const std::string& code, @@ -175,6 +177,10 @@ Result BuildBaseline(const Domain& domain) { Result> BuildRuntimeShells(const Domain& domain, const DofManager& dofs) { std::vector shells; + if (domain.Elements().Empty()) + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "unsupported-nonlinear-element", {}, "ELEMENT", + {}, "The nonlinear static procedure requires a MITC4 shell element.")); shells.reserve(domain.Elements().Size()); for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) { const auto* definition = @@ -320,7 +326,7 @@ NonlinearMitc4TrialDofs ExtractElementTrial( Result AssembleInternal( const std::vector& shells, const DofManager& dofs, - const Vector& coordinates, const ParallelFor& parallel_for) { + const Domain& domain, const Vector& coordinates, const ParallelFor& parallel_for) { std::vector> local( shells.size()); std::vector> failures(shells.size()); @@ -341,6 +347,7 @@ Result AssembleInternal( std::vector rows; rows.reserve(shells.size() * 4U); double energy = 0.0; + Mitc4DrillingEvidence drilling{shells.front().kernel->DrillingMode()}; for (std::size_t index = 0U; index < shells.size(); ++index) { if (failures[index].has_value()) { return Result::Failure(*failures[index]); @@ -354,6 +361,28 @@ Result AssembleInternal( "NONLINEAR_STATIC_PROCEDURE", std::to_string(index), "A nonlinear shell internal residual or tangent is nonfinite.")); } + if (shells[index].kernel->DrillingMode() != drilling.mode || + local[index]->drilling_force.Size() != kShellDofCount || + !IsFinite(local[index]->drilling_force)) + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "invalid-shell-drilling-evidence", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, "Drilling mode or force is invalid.")); + const auto element_trial = ExtractElementTrial(shells[index], coordinates); + for (std::size_t node = 0U; node < 4U; ++node) { + const auto first = 6U * node; + const auto global_node = shells[index].scatter[first] / 6U; + const auto& values = local[index]->drilling_force; + const Vector3 drilling_force{values[first], values[first + 1U], values[first + 2U]}; + const auto spatial = SpatialMomentSO3(element_trial.rotation_vector[node], + {values[first + 3U], values[first + 4U], values[first + 5U]}); + if (!spatial.HasValue()) + return Result::Failure(spatial.GetStatus()); + const Vector3 position = Vector3(domain.Nodes()[global_node].coordinates) + + element_trial.displacement[node]; + const Vector3 moment = position.Cross(drilling_force) + spatial.Value(); + for (std::size_t c = 0U; c < 3U; ++c) + drilling.drilling_moment[c] += moment[c]; + } for (std::size_t row = 0U; row < kShellDofCount; ++row) { force[shells[index].scatter[row]] += local[index]->internal_residual[row]; assembly_internal::AppendLocalMatrixRow( @@ -373,7 +402,8 @@ Result AssembleInternal( } energy += local[index]->physical_strain_energy; } - if (!IsFinite(force) || !std::isfinite(energy)) { + if (!IsFinite(force) || !std::isfinite(energy) || + !IsFinite(drilling.drilling_moment)) { return Result::Failure(ProcedureFailure( FailureCategory::kSolver, "nonfinite-newton-data", {}, "NONLINEAR_STATIC_PROCEDURE", {}, @@ -387,7 +417,7 @@ Result AssembleInternal( return Result::Failure(tangent.GetStatus()); } return Result::Success( - {std::move(force), std::move(tangent.Value()), std::move(rows), energy}); + {std::move(force), std::move(tangent.Value()), std::move(rows), energy, drilling}); } Result Difference(const Vector& left, const Vector& right) { @@ -699,7 +729,7 @@ AcceptedFrameState AcceptedFrame(const std::uint32_t frame_index, trial.reaction, trial.shell_results, trial.physical_strain_energy, - trial.global_equilibrium}; + trial.global_equilibrium, trial.mitc4_drilling}; } } // namespace @@ -807,7 +837,7 @@ Result ExecuteNonlinearStaticCore( NewtonIterationStage::kKinematicsEvaluated; Result internal = AssembleInternal( - shells.Value(), dofs.Value(), + shells.Value(), dofs.Value(), domain, trial.Value().generalized_coordinates, parallel_for); if (!internal.HasValue()) { return Result>::Failure( @@ -879,6 +909,9 @@ Result ExecuteNonlinearStaticCore( NewtonIterationStage::kExternalAssembled)); } trial.Value().global_equilibrium = global_equilibrium.Value(); + trial.Value().mitc4_drilling = internal.Value().mitc4_drilling; + results_internal::CompleteDrillingEvidence( + trial.Value().global_equilibrium, &*trial.Value().mitc4_drilling); Result partition = constraints.Partition(tangent.Value(), dofs.Value()); diff --git a/src/fesa/elements/element_capabilities.cpp b/src/fesa/elements/element_capabilities.cpp index 772fc0c..f780864 100644 --- a/src/fesa/elements/element_capabilities.cpp +++ b/src/fesa/elements/element_capabilities.cpp @@ -95,6 +95,10 @@ Result Ds4ThermalAdapter::Recover( Mitc4ReferenceAdapter::Mitc4ReferenceAdapter(const Mitc4Shell& kernel) : kernel_{kernel} {} +Mitc4DrillingMode Mitc4ReferenceAdapter::DrillingMode() const noexcept { + return kernel_.DrillingMode(); +} + const ElementDofLayout& Mitc4ReferenceAdapter::DofLayout() const noexcept { return kernel_.DofLayout(); } diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp index 9c8c756..ef7f002 100644 --- a/src/fesa/elements/mitc4_shell.cpp +++ b/src/fesa/elements/mitc4_shell.cpp @@ -958,6 +958,22 @@ Result Mitc4Shell::Recover( ShellElementResultRows rows{}; rows.rows.reserve(kNodeCount); rows.physical_strain_energy = recovered.Value().strain_energy; + const auto stiffness = Stiffness(); + if (!stiffness.HasValue()) + return Result::Failure(stiffness.GetStatus()); + const Vector drilling_force = + stiffness.Value().drilling_global24.Multiply(element_displacement); + rows.drilling_mode = DrillingMode(); + Vector3 drilling_moment; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const std::size_t first = 6U * node; + const Vector3 force{drilling_force[first], drilling_force[first + 1U], + drilling_force[first + 2U]}; + const Vector3 moment{drilling_force[first + 3U], drilling_force[first + 4U], + drilling_force[first + 5U]}; + drilling_moment = drilling_moment + coordinates_[node].Cross(force) + moment; + } + rows.drilling_moment = drilling_moment.Components(); for (std::size_t point = 0U; point < recovered.Value().points.size(); ++point) { const auto& physical_point = recovered.Value().points[point]; diff --git a/src/fesa/results/analysis_state_commit.cpp b/src/fesa/results/analysis_state_commit.cpp index df6ed96..6058e83 100644 --- a/src/fesa/results/analysis_state_commit.cpp +++ b/src/fesa/results/analysis_state_commit.cpp @@ -48,7 +48,8 @@ Status CommitRecoveryCandidate(RecoveryCandidate candidate, !candidate.shell.rows.empty() || candidate.shell.physical_strain_energy != 0.0 || candidate.shell.equilibrium != zero_equilibrium || - candidate.shell.verification_metrics != zero_metrics; + candidate.shell.verification_metrics != zero_metrics || + candidate.shell.mitc4_drilling.has_value(); const bool has_complete_shell_inventory = !candidate.expected_shell_elements.empty() && !candidate.shell.rows.empty(); diff --git a/src/fesa/results/global_equilibrium_recovery.cpp b/src/fesa/results/global_equilibrium_recovery.cpp index 3f36bca..cd44b7b 100644 --- a/src/fesa/results/global_equilibrium_recovery.cpp +++ b/src/fesa/results/global_equilibrium_recovery.cpp @@ -10,6 +10,7 @@ #include #include "fesa/math/vector3.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa::results_internal { namespace { @@ -243,11 +244,24 @@ Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs, candidate.shell.equilibrium[component] = force_balance[component]; candidate.shell.equilibrium[3U + component] = moment_balance[component]; } + Vector3 acceptance_moment = moment_balance; + if (candidate.shell.mitc4_drilling) { + CompleteDrillingEvidence(candidate.shell.equilibrium, + &*candidate.shell.mitc4_drilling); + if (!ValidDrillingEvidence(candidate.shell.mitc4_drilling, + candidate.shell.equilibrium)) + return RecoveryFailure("invalid-shell-drilling-evidence", + {domain.SourcePath(), 0U}, "global-equilibrium", + "Shell drilling diagnostics must be finite and consistent."); + if (candidate.shell.mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) + acceptance_moment = moment_balance - + Vector3(candidate.shell.mitc4_drilling->drilling_moment); + } const double force_metric = NormalizedBalance(force_balance.Components(), (std::max)(applied_force_scale, reaction_force_scale)); const double moment_metric = NormalizedBalance( - moment_balance.Components(), + acceptance_moment.Components(), (std::max)(applied_moment_scale, reaction_moment_scale)); candidate.shell.verification_metrics = {global.normalized_free_residual, force_metric, moment_metric}; diff --git a/src/fesa/results/mitc4_drilling_evidence.h b/src/fesa/results/mitc4_drilling_evidence.h new file mode 100644 index 0000000..29ed8b3 --- /dev/null +++ b/src/fesa/results/mitc4_drilling_evidence.h @@ -0,0 +1,54 @@ +#ifndef FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ +#define FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ + +#include +#include +#include +#include + +#include "fesa/results/result_records.h" + +namespace fesa::results_internal { + +inline bool ValidDrillingMode(const Mitc4DrillingMode mode) { + return mode == Mitc4DrillingMode::kDirectSpring || + mode == Mitc4DrillingMode::kLegacyCorotatedTwist; +} + +/// @brief Validates exact raw-minus-internal diagnostics without a tolerance. +inline bool ValidDrillingEvidence( + const std::optional& evidence, + const std::array& raw) { + if (!evidence) return true; // Older manually constructed fixtures. + if (!ValidDrillingMode(evidence->mode)) return false; + for (std::size_t i = 0U; i < 6U; ++i) { + const double expected = i < 3U ? raw[i] + : raw[i] - evidence->drilling_moment[i - 3U]; + if (!std::isfinite(raw[i]) || !std::isfinite(expected) || + !std::isfinite(evidence->equilibrium_with_drilling[i]) || + evidence->equilibrium_with_drilling[i] != expected) + return false; + } + for (const double value : evidence->drilling_moment) + if (!std::isfinite(value)) return false; + return true; +} + +inline bool SameDrillingInventory( + const std::optional& left, + const std::optional& right) { + return left.has_value() == right.has_value() && + (!left || left->mode == right->mode); +} + +/// @brief Keeps force balance raw and subtracts positive artificial moment once. +inline void CompleteDrillingEvidence( + const std::array& raw, Mitc4DrillingEvidence* evidence) { + evidence->equilibrium_with_drilling = raw; + for (std::size_t i = 0U; i < 3U; ++i) + evidence->equilibrium_with_drilling[3U + i] -= evidence->drilling_moment[i]; +} + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ diff --git a/src/fesa/results/shell_result_recovery.cpp b/src/fesa/results/shell_result_recovery.cpp index b7c00da..b443b0f 100644 --- a/src/fesa/results/shell_result_recovery.cpp +++ b/src/fesa/results/shell_result_recovery.cpp @@ -4,6 +4,8 @@ #include #include +#include "results/mitc4_drilling_evidence.h" + namespace fesa::results_internal { bool IsFiniteShellRow(const ShellResultRow& row) { @@ -29,6 +31,27 @@ Status ShellResultRecovery(const SourceLocation& location, const EntityIndex expected_element, const ShellElementResultRows& rows, RecoveryCandidate& candidate) { + if ((rows.drilling_mode && !ValidDrillingMode(*rows.drilling_mode)) || + !IsFiniteArray(rows.drilling_moment) || + (candidate.has_shell_results && + (candidate.shell.mitc4_drilling.has_value() != rows.drilling_mode.has_value() || + (rows.drilling_mode && + candidate.shell.mitc4_drilling->mode != *rows.drilling_mode)))) { + return RecoveryFailure("invalid-shell-drilling-evidence", location, + source_id.source_label_text, + "Shell recovery requires one finite actual drilling mode."); + } + if (rows.drilling_mode) { + if (!candidate.shell.mitc4_drilling) + candidate.shell.mitc4_drilling = Mitc4DrillingEvidence{*rows.drilling_mode}; + for (std::size_t i = 0U; i < 3U; ++i) { + candidate.shell.mitc4_drilling->drilling_moment[i] += rows.drilling_moment[i]; + if (!std::isfinite(candidate.shell.mitc4_drilling->drilling_moment[i])) + return RecoveryFailure("nonfinite-recovery-value", location, + source_id.source_label_text, + "Source-ordered drilling moment overflowed."); + } + } const double accumulated_energy = candidate.shell.physical_strain_energy + rows.physical_strain_energy; if (!std::isfinite(rows.physical_strain_energy) || diff --git a/tests/integration/analysis/linear_frequency_procedure_test.cpp b/tests/integration/analysis/linear_frequency_procedure_test.cpp index cd15396..fab87a8 100644 --- a/tests/integration/analysis/linear_frequency_procedure_test.cpp +++ b/tests/integration/analysis/linear_frequency_procedure_test.cpp @@ -18,6 +18,9 @@ #include "fesa/analysis/prepared_analysis.h" #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" +#include "fesa/elements/element_factory.h" +#include "fesa/results/result_recovery.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/fem/dof_manager.h" #include "fesa/model/domain.h" #include "fesa/model/model_types.h" @@ -132,6 +135,8 @@ TEST(LinearFrequencyProcedureExecutor, fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value()); ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); const auto& state = result.Value(); + ASSERT_TRUE(state.drilling_mode); + EXPECT_EQ(*state.drilling_mode, fesa::Mitc4DrillingMode::kDirectSpring); ASSERT_EQ(state.modes.size(), 1U); EXPECT_EQ(state.mass_null_dimension, 0U); EXPECT_EQ(state.finite_dimension, 1U); @@ -189,10 +194,10 @@ TEST(LinearFrequencyProcedureExecutor, } } -// E2-C09 is a pencil oracle. The public static contract conflict is recorded -// separately; direct covariance does not certify spatial moment equilibrium. +// E2-C09 keeps raw spatial imbalance visible while the direct-spring gate +// subtracts only the separately recovered artificial internal moment. TEST(LinearFrequencyProcedureExecutor, - PartialBoundaryStaticRejectsUnbalancedDrillingMoment) { + PartialBoundaryStaticSeparatesArtificialMomentAndPreservesLegacyGate) { const auto domain = PanelDomain(false, false, true); ASSERT_TRUE(domain.HasValue()); const fesa::SerialParallelFor serial; @@ -200,9 +205,32 @@ TEST(LinearFrequencyProcedureExecutor, const auto result = fesa::LinearStaticProcedureExecutor{serial, solver}.Execute( domain.Value()); - ASSERT_FALSE(result.HasValue()); - EXPECT_NE(StatusText(result.GetStatus()) - .find("global-equilibrium-tolerance-failure"), + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + const auto& state = result.Value().state; + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_NE(state.Equilibrium()[5U], 0.0); + for (std::size_t c = 0U; c < 3U; ++c) { + EXPECT_NEAR(state.Equilibrium()[3U + c], state.Mitc4Drilling()->drilling_moment[c], 1e-14); + EXPECT_NEAR(state.Mitc4Drilling()->equilibrium_with_drilling[3U + c], 0.0, 1e-14); + } + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + auto element = fesa::ElementFactory{}.Create(domain.Value().Elements()[0U], domain.Value()); + ASSERT_TRUE(element.HasValue()); + auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + *shell = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell, fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + fesa::ElementView view{std::cref(*shell)}; + auto stiffness = fesa::SparseAssembler::Assemble(view, dofs.Value(), serial); + ASSERT_TRUE(stiffness.HasValue()); + auto legacy_state = fesa::AnalysisState::Create(dofs.Value(), {"Step-1", 0U}); + const auto legacy = fesa::ResultRecovery::Recover(model.Value(), view, dofs.Value(), + stiffness.Value(), state.Displacement(), state.ExternalForce(), legacy_state); + ASSERT_FALSE(legacy.IsOk()); + EXPECT_NE(StatusText(legacy).find("global-equilibrium-tolerance-failure"), std::string::npos); } diff --git a/tests/integration/analysis/nonlinear_static_procedure_test.cpp b/tests/integration/analysis/nonlinear_static_procedure_test.cpp index 8f41c34..cf4bcd6 100644 --- a/tests/integration/analysis/nonlinear_static_procedure_test.cpp +++ b/tests/integration/analysis/nonlinear_static_procedure_test.cpp @@ -16,6 +16,13 @@ #include #include "fesa/elements/finite_rotation.h" +#include "fesa/analysis/prepared_analysis.h" +#include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/mitc4_shell.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" #include "fesa/math/vector3.h" namespace { @@ -343,3 +350,57 @@ TEST(NonlinearStaticProcedure, } } // namespace + +TEST(Mitc4DrillingEvidence, NonlinearAndPressureFramesRetainSignedTorqueOracle) { + for (const bool pressure : {false, true}) { + SCOPED_TRACE(pressure); + TempDirectory directory{"drilling"}; + auto deck = AllConstrainedNonlinearShellDeck(); + const auto boundary = deck.find("All, 1, 6"); + ASSERT_NE(boundary, std::string::npos); + deck.replace(boundary, std::string{"All, 1, 6"}.size(), "All, 1, 5\nAll, 6, 6, 0.125"); + if (pressure) { + const auto assembly_end = deck.find("*End Assembly"); + deck.insert(assembly_end, "*Surface, TYPE=ELEMENT, NAME=Pressure\nShell-1.1, SPOS\n"); + deck.insert(deck.find("*End Step"), "*Dsload\nPressure, P, 4.0\n"); + } + const auto path = directory.Path() / "drilling.inp"; + WriteText(path, deck); + auto prepared = fesa::PrepareAnalysisInput(path); + ASSERT_TRUE(prepared.HasValue()); + const auto& domain = prepared.Value().domain; + auto element = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain); + ASSERT_TRUE(element.HasValue()); + const auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + const auto stiffness = shell->Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double trace = 0.0; + for (std::size_t node = 0U; node < 4U; ++node) + for (std::size_t c = 3U; c < 5U; ++c) + trace += stiffness.Value().physical_local20(5U * node + c, 5U * node + c); + const double kd = 1e-3 * trace / 8.0; + fesa::SerialParallelFor serial; + fesa::MklPardisoTangentSolver tangent; + fesa::NewtonIterationLogWriter log; + ASSERT_TRUE(log.Open(directory.Path() / "drilling.log").IsOk()); + auto result = fesa::NonlinearStaticProcedureExecutor{serial, tangent}.Execute(domain, log); + ASSERT_TRUE(log.Close().IsOk()); + ASSERT_TRUE(result.HasValue()) << (result.GetStatus().Diagnostics().empty() ? "" : result.GetStatus().Diagnostics()[0U].message); + ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 10U); + for (const auto& frame : result.Value().state.AcceptedFrames()) { + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_EQ(frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + const double torque = 4.0 * kd * 0.125 * frame.load_factor; + EXPECT_NEAR(frame.global_equilibrium[5U], torque, 1e-13); + EXPECT_NEAR(frame.mitc4_drilling->drilling_moment[2U], torque, 1e-13); + for (const double value : frame.mitc4_drilling->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-13); + EXPECT_NEAR(frame.physical_strain_energy, 0.0, 1e-20); + if (pressure) EXPECT_NE(frame.external_force[2U], 0.0); + } + ASSERT_TRUE(result.Value().state.CommittedState().mitc4_drilling); + EXPECT_EQ(result.Value().state.CommittedState().mitc4_drilling->drilling_moment, + result.Value().state.AcceptedFrames().back().mitc4_drilling->drilling_moment); + } +} diff --git a/tests/unit/analysis/dynamic_analysis_state_test.cpp b/tests/unit/analysis/dynamic_analysis_state_test.cpp index 9a7faf9..efea398 100644 --- a/tests/unit/analysis/dynamic_analysis_state_test.cpp +++ b/tests/unit/analysis/dynamic_analysis_state_test.cpp @@ -320,3 +320,27 @@ TEST(DynamicAnalysisState, CreateRejectsInvalidScheduleAndNonzeroBoundary) { } } // namespace + +TEST(DynamicAnalysisState, DrillingEvidenceRejectsPartialCommitAndModeChanges) { + auto fixture = MakeFixture(); + auto initial = MakeFrame(*fixture.domain, 0U); + initial.mitc4_drilling = fesa::Mitc4DrillingEvidence{fesa::Mitc4DrillingMode::kDirectSpring}; + ASSERT_TRUE(fixture.state.CommitAcceptedFrame(initial).IsOk()); + for (const int kind : {0, 1, 2, 3, 4}) { + auto next = MakeFrame(*fixture.domain, 1U); + next.mitc4_drilling = initial.mitc4_drilling; + if (kind == 0) next.mitc4_drilling.reset(); + if (kind == 1) next.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) next.mitc4_drilling->mode = static_cast(99); + if (kind == 3) next.mitc4_drilling->equilibrium_with_drilling[3U] = 1.0; + if (kind == 4) next.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::quiet_NaN(); + EXPECT_FALSE(fixture.state.CommitAcceptedFrame(next).IsOk()) << kind; + ASSERT_EQ(fixture.state.AcceptedFrames().size(), 1U); + EXPECT_EQ(fixture.state.CommittedState().frame_index, 0U); + EXPECT_EQ(fixture.state.CommittedState().mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + } + auto valid = MakeFrame(*fixture.domain, 1U); + valid.mitc4_drilling = initial.mitc4_drilling; + EXPECT_TRUE(fixture.state.CommitAcceptedFrame(valid).IsOk()); +} diff --git a/tests/unit/analysis/hht_dynamic_procedure_test.cpp b/tests/unit/analysis/hht_dynamic_procedure_test.cpp index 75594e9..4b9ca29 100644 --- a/tests/unit/analysis/hht_dynamic_procedure_test.cpp +++ b/tests/unit/analysis/hht_dynamic_procedure_test.cpp @@ -216,6 +216,12 @@ void ExpectSameDynamicHistory(const RigidTranslationRun& left, EXPECT_DOUBLE_EQ(lhs.drilling_energy, rhs.drilling_energy); EXPECT_DOUBLE_EQ(lhs.external_work, rhs.external_work); EXPECT_DOUBLE_EQ(lhs.energy_defect, rhs.energy_defect); + ASSERT_TRUE(lhs.mitc4_drilling); + ASSERT_TRUE(rhs.mitc4_drilling); + EXPECT_EQ(lhs.mitc4_drilling->mode, rhs.mitc4_drilling->mode); + EXPECT_EQ(lhs.mitc4_drilling->drilling_moment, rhs.mitc4_drilling->drilling_moment); + EXPECT_EQ(lhs.mitc4_drilling->equilibrium_with_drilling, + rhs.mitc4_drilling->equilibrium_with_drilling); } } @@ -357,6 +363,10 @@ TEST(HhtDynamicProcedure, ClampedStepLoadHasPhysicalReactionAndZeroMotion) { EXPECT_EQ(frame.time_integration_formulation, fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht); EXPECT_EQ(frame.iteration_count, 0U); + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_EQ(frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(frame.mitc4_drilling->drilling_moment, (std::array{})); + EXPECT_EQ(frame.mitc4_drilling->equilibrium_with_drilling, frame.global_equilibrium); EXPECT_DOUBLE_EQ(frame.reaction[2U], -4.0); EXPECT_DOUBLE_EQ(frame.external_force[2U], 4.0); EXPECT_DOUBLE_EQ(frame.kinetic_energy, 0.0); diff --git a/tests/unit/analysis/hht_reduced_shell_test.cpp b/tests/unit/analysis/hht_reduced_shell_test.cpp index 752198a..4e0d477 100644 --- a/tests/unit/analysis/hht_reduced_shell_test.cpp +++ b/tests/unit/analysis/hht_reduced_shell_test.cpp @@ -22,10 +22,13 @@ #include "analysis/hht_reduced_internal.h" #include "math/dense_blas_internal.h" +#include "fesa/elements/mitc4_shell.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/analysis/analysis_model.h" #include "fesa/analysis/hht_time_integration.h" #include "fesa/assembly/parallel_for.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" namespace { @@ -625,6 +628,10 @@ TEST(HhtReducedOperators, [](Assembly& a) { a.trial_constraints->shells[0U].configuration.identity += "x"; }, + [](Assembly& a) { + a.trial_constraints->shells[0U].configuration.drilling_mode = + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + }, [](Assembly& a) { a.trial_constraints->shells[0U].q[0U] = -0.0; }, [](Assembly& a) { a.trial_constraints->shells[0U].v[0U] = -0.0; }, [](Assembly& a) { @@ -935,11 +942,15 @@ TEST(HhtReducedOperators, NonfiniteMatchingTimeKeyMissPreservesFreshSuccess) { } } -TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { +void CheckShellConstraintScatter(const fesa::Mitc4DrillingMode mode) { const auto fixture = MakeFixture(false); const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); ASSERT_TRUE(shells.HasValue()); + SCOPED_TRACE(static_cast(mode)); + for (const auto& shell : shells.Value()) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, mode); fesa::DynamicAcceptedFrameState frame{}; for (auto* values : {&frame.generalized_coordinates, &frame.velocity, &frame.acceleration, &frame.algorithmic_acceleration, @@ -980,7 +991,10 @@ TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { assembly.Value().frame.residual_scale[6U * node + 3U + component]; EXPECT_DOUBLE_EQ(row.value_scale, expected_scale); EXPECT_FALSE(row.gradient.empty()); - EXPECT_FALSE(row.hessian.empty()); + if (mode == fesa::Mitc4DrillingMode::kLegacyCorotatedTwist) + EXPECT_FALSE(row.hessian.empty()); + else + EXPECT_TRUE(row.hessian.empty()); EXPECT_LT(row.gradient.size(), fixture.dofs.FullDofCount()); for (const auto& entry : row.hessian) { EXPECT_LT(entry.row, fixture.dofs.FullDofCount()); @@ -989,6 +1003,14 @@ TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { } } +TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { + CheckShellConstraintScatter(fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + +TEST(HhtReducedOperators, DirectSpringScattersZeroInitialHessianByConstraintRow) { + CheckShellConstraintScatter(fesa::Mitc4DrillingMode::kDirectSpring); +} + TEST(HhtReducedEmbedding, RejectsNonfiniteGraphResidualBeforeAcceptance) { const auto fixture = MakeFixture(false); const auto chart = fesa::hht_reduced_internal::BuildFixedChart( @@ -1184,11 +1206,15 @@ TEST(HhtReducedEmbedding, EXPECT_DOUBLE_EQ(rates.Value().z_acceleration[0U], 1.0); } -TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { +void CheckActualShellRatesWithMixedTimeLoad(const fesa::Mitc4DrillingMode mode) { const auto fixture = MakeFixture(false, true); const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); ASSERT_TRUE(shells.HasValue()); + SCOPED_TRACE(static_cast(mode)); + for (const auto& shell : shells.Value()) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, mode); const fesa::Vector zero(fixture.dofs.FullDofCount()); const auto initial = Assemble(fixture, shells.Value(), zero, zero, zero, zero, true); @@ -1470,13 +1496,26 @@ TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { ++row) { auto& omitted = omitted_external_gradient.rows[row].gradient; const auto& no_force = no_force_operators.Value().rows[row].gradient; - ASSERT_EQ(omitted.size(), no_force.size()); - for (std::size_t entry = 0U; entry < omitted.size(); ++entry) { - ASSERT_EQ(omitted[entry].dof, no_force[entry].dof); - saw_external_position_gradient = - saw_external_position_gradient || - omitted[entry].value != no_force[entry].value; - omitted[entry].value = no_force[entry].value; + if (mode == fesa::Mitc4DrillingMode::kLegacyCorotatedTwist) { + ASSERT_EQ(omitted.size(), no_force.size()); + for (std::size_t entry = 0U; entry < omitted.size(); ++entry) { + ASSERT_EQ(omitted[entry].dof, no_force[entry].dof); + saw_external_position_gradient = + saw_external_position_gradient || + omitted[entry].value != no_force[entry].value; + omitted[entry].value = no_force[entry].value; + } + } else { + // Direct internal/external derivatives can cancel exact sparse entries. + // Compare the full union before replacing only the external contribution. + fesa::Vector with_force(fixture.dofs.FullDofCount()); + fesa::Vector without_force(fixture.dofs.FullDofCount()); + for (const auto& entry : omitted) with_force[entry.dof] = entry.value; + for (const auto& entry : no_force) without_force[entry.dof] = entry.value; + for (std::size_t dof = 0U; dof < with_force.Size(); ++dof) + saw_external_position_gradient = saw_external_position_gradient || + with_force[dof] != without_force[dof]; + omitted = no_force; } } ASSERT_TRUE(saw_external_position_gradient); @@ -1560,6 +1599,14 @@ TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { EXPECT_GT(omitted_second_time_hidden, 1.0e-8); } +TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { + CheckActualShellRatesWithMixedTimeLoad(fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + +TEST(HhtReducedEmbedding, DirectSpringReconstructsActualShellRatesWithMixedTimeLoad) { + CheckActualShellRatesWithMixedTimeLoad(fesa::Mitc4DrillingMode::kDirectSpring); +} + TEST(HhtReducedChannelNewton, RequiresEveryTranslationGateAndThaws) { fesa::hht_reduced_internal::FixedChart chart{}; chart.full_dof_count = 6U; @@ -1737,3 +1784,158 @@ TEST(HhtReducedOperators, NullMetricKeepsFiniteLeftAssociatedProducts) { } } // namespace + +TEST(Mitc4DrillingEvidence, HhtActualMomentMatchesClosedFormAndSourceOrder) { + const auto fixture = MakeFixture(false, true); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + for (std::size_t node = 0U; node < fixture.domain->Nodes().size(); ++node) { + q[6U * node + 0U] = 0.03 * static_cast(node); + q[6U * node + 3U] = 0.13 + 0.01 * static_cast(node); + q[6U * node + 4U] = -0.09; + q[6U * node + 5U] = 0.21 - 0.01 * static_cast(node); + } + const auto serial = Assemble(fixture, shells.Value(), q, zero, zero, zero, true); + ASSERT_TRUE(serial.HasValue()); + ASSERT_TRUE(serial.Value().frame.mitc4_drilling); + fesa::Vector3 expected; + for (std::size_t element = 0U; element < fixture.domain->ShellElements().Size(); ++element) { + const auto& definition = fixture.domain->ShellElements()[element]; + std::array nodes{}; + std::array, 4> directors{}; + for (std::size_t i = 0U; i < 4U; ++i) { + nodes[i] = &fixture.domain->Nodes()[definition.node_indices[i]]; + directors[i] = fixture.domain->ShellNodeInitialFrames()[definition.node_indices[i]].director; + } + const auto* section = dynamic_cast( + &fixture.domain->Properties()[definition.section_index]); + const auto* material = dynamic_cast( + &fixture.domain->Materials()[definition.material_index]); + ASSERT_NE(section, nullptr); + ASSERT_NE(material, nullptr); + const auto kernel = fesa::Mitc4Shell::Create(nodes, directors, *section, *material); + ASSERT_TRUE(kernel.HasValue()); + const auto stiffness = kernel.Value().Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double trace = 0.0; + for (std::size_t i = 0U; i < 4U; ++i) + for (std::size_t c = 3U; c < 5U; ++c) + trace += stiffness.Value().physical_local20(5U * i + c, 5U * i + c); + const double kd = 1e-3 * trace / 8.0; + for (std::size_t i = 0U; i < 4U; ++i) { + const auto node = definition.node_indices[i]; + const fesa::Vector3 phi{q[6U * node + 3U], q[6U * node + 4U], q[6U * node + 5U]}; + const fesa::Vector3 director{directors[i]}; + const double angle = phi.Norm(); + const double coefficient = (1.0 - 0.5 * angle / std::tan(0.5 * angle)) / (angle * angle); + // Independent closed inverse-transpose Rodrigues identity; no dense solve. + expected = expected + kd * phi.Dot(director) * + (director + 0.5 * phi.Cross(director) + coefficient * phi.Cross(phi.Cross(director))); + } + } + for (std::size_t c = 0U; c < 3U; ++c) + EXPECT_NEAR(serial.Value().frame.mitc4_drilling->drilling_moment[c], expected[c], 2e-14); + EXPECT_EQ(serial.Value().frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + auto frame = serial.Value().frame; + const auto parameters = fesa::MakeHhtParameters(-0.05); + ASSERT_TRUE(parameters.HasValue()); + fesa::TbbParallelFor parallel; + const auto tbb = fesa::hht_internal::EvaluateAssembly(*fixture.domain, fixture.dofs, + shells.Value(), frame, zero, 0.1, parameters.Value(), parallel, true); + ASSERT_TRUE(tbb.HasValue()); + ASSERT_TRUE(tbb.Value().frame.mitc4_drilling); + for (std::size_t c = 0U; c < 3U; ++c) + ExpectSameBits(tbb.Value().frame.mitc4_drilling->drilling_moment[c], + serial.Value().frame.mitc4_drilling->drilling_moment[c]); +} + +TEST(Mitc4DrillingEvidence, HhtCorrectedGateRejectsOppositeSignAndRealImbalance) { + const auto fixture = MakeFixture(false); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + for (std::size_t node = 0U; node < 4U; ++node) q[6U * node + 5U] = 0.125; + const auto assembly = Assemble(fixture, shells.Value(), q, zero, zero, zero, false); + ASSERT_TRUE(assembly.HasValue()); + EXPECT_EQ(assembly.Value().hidden_constraint.Size(), 0U); + const auto directional = Assemble(fixture, shells.Value(), q, zero, zero, zero, true); + ASSERT_TRUE(directional.HasValue()); + ASSERT_EQ(directional.Value().hidden_constraint.Size(), 4U); + for (std::size_t node = 0U; node < 4U; ++node) + EXPECT_GT(directional.Value().hidden_constraint[node], 0.0); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + auto unavailable_frame = assembly.Value().frame; + const auto unavailable = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, + fesa::hht_internal::NullNodes(fixture.dofs, assembly.Value().rotations), + positive, tangent, &unavailable_frame); + EXPECT_FALSE(unavailable.IsOk()); + EXPECT_EQ(unavailable_frame.hidden_constraint.Size(), 0U); + auto frame = assembly.Value().frame; + // A support supplying the complete internal force gives exact static balance. + frame.reaction = frame.internal_force; + ASSERT_TRUE(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs, + assembly.Value(), &frame).IsOk()); + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_GT(frame.global_equilibrium[5U], 0.0); + EXPECT_NEAR(frame.global_equilibrium[5U], frame.mitc4_drilling->drilling_moment[2U], 1e-14); + for (const double value : frame.mitc4_drilling->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-14); + for (const int kind : {0, 1, 2, 3}) { + auto invalid = frame; + if (kind == 0) invalid.mitc4_drilling->drilling_moment[2U] *= -1.0; + if (kind == 1) invalid.reaction[0U] += 1.0; + if (kind == 2) invalid.reaction[4U] += 1.0; + if (kind == 3) invalid.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + EXPECT_FALSE(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs, + assembly.Value(), &invalid).IsOk()) << kind; + } +} + +TEST(Mitc4DrillingEvidence, DirectHiddenCertificateUsesCanonicalPrimitive) { + const auto fixture = MakeFixture(false); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + auto external = zero; + fesa::NonlinearMitc4TrialDofs trial{}; + for (std::size_t node = 0U; node < 4U; ++node) { + q[6U * node] = 0.03 * static_cast(node); + trial.displacement[node] = {q[6U * node], 0.0, 0.0}; + trial.rotation_vector[node] = {0.00005, 0.000035, 0.00002}; + for (std::size_t c = 0U; c < 3U; ++c) + q[6U * node + 3U + c] = trial.rotation_vector[node][c]; + } + const auto local = shells.Value()[0U].internal->EvaluateTrial(trial, &zero); + ASSERT_TRUE(local.HasValue()); + external = local.Value().drilling_force; + const auto assembly = Assemble(fixture, shells.Value(), q, zero, zero, external, true); + ASSERT_TRUE(assembly.HasValue()); + const auto length = fesa::hht_internal::ReferenceLength(*fixture.domain); + ASSERT_TRUE(length.HasValue()); + const auto null_nodes = fesa::hht_internal::NullNodes(fixture.dofs, assembly.Value().rotations); + ASSERT_EQ(null_nodes.size(), 4U); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + auto frame = assembly.Value().frame; + const auto status = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), length.Value(), null_nodes, + positive, tangent, &frame); + ASSERT_TRUE(status.IsOk()) << (status.Diagnostics().empty() ? "" : status.Diagnostics()[0U].message); + const auto chart = fesa::hht_reduced_internal::BuildFixedChart(fixture.dofs, assembly.Value().rotations); + ASSERT_TRUE(chart.HasValue()); + const auto canonical = fesa::hht_reduced_internal::EvaluateConstraintOperators( + *fixture.domain, fixture.dofs, shells.Value(), assembly.Value(), chart.Value(), zero, zero); + ASSERT_TRUE(canonical.HasValue()); + for (std::size_t row = 0U; row < 4U; ++row) { + EXPECT_DOUBLE_EQ(frame.hidden_constraint[row], canonical.Value().rows[row].value); + EXPECT_NEAR(frame.hidden_constraint[row], 0.0, 1e-17); + EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[row], 0.0); + EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[row], 0.0); + } +} diff --git a/tests/unit/analysis/hht_shell_physics_test.cpp b/tests/unit/analysis/hht_shell_physics_test.cpp index abac0f0..733ec91 100644 --- a/tests/unit/analysis/hht_shell_physics_test.cpp +++ b/tests/unit/analysis/hht_shell_physics_test.cpp @@ -28,6 +28,7 @@ #include "fesa/analysis/newton_iteration_log.h" #include "fesa/assembly/parallel_for.h" #include "fesa/elements/finite_rotation.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" @@ -120,6 +121,14 @@ struct Fixture { double length; }; +// The finite rigid-orbit and free-spin C08 fixtures prescribe the historical +// objective potential and its hidden velocity equation explicitly. +void UseLegacyDrilling(const Fixture& fixture) { + for (const auto& shell : fixture.shells) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + Fixture MakeFixture(const std::size_t divisions, const bool lifecycle_translation = false) { fesa::ModelDefinition definition{}; @@ -643,6 +652,7 @@ TEST(HhtShellPhysics, ManufacturedHiddenCountsExternalHessianOnce) { TEST(HhtShellPhysics, C08aRigidOrbitMatchesIndependentVolumeForcePowerAndMomentum) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); const auto parameters = Require(fesa::MakeHhtParameters(-0.05)); const auto initial = Evaluate(fixture, EmptyFrame(fixture), 1.0 / 16.0, parameters, false); @@ -3078,6 +3088,7 @@ std::array AccelerationSplit(const Frame& basis, TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); const auto grid = RunPlate(fixture, 1.0 / 64.0, true); ASSERT_EQ(grid.history.size(), 65U); const auto exact = RigidOrbit(fixture, 0.0); @@ -3115,6 +3126,7 @@ TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) { TEST(HhtShellPhysics, C08aRigidManufacturedOrbitReducedCoordinateTimeConvergence) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); std::vector grids; std::vector> errors; for (const unsigned denominator : {16U, 32U, 64U, 128U, 256U, 512U, 1024U}) { @@ -3297,6 +3309,7 @@ TEST(HhtShellPhysics, C08bInitialPhysicalSpectrumProbe) { TEST(HhtShellPhysics, PredictorStartCompletesPreviouslyFailingLumpedFreePlateGrid) { const auto fixture = MakeFixture(2U); + UseLegacyDrilling(fixture); const auto evidence = RunPlate(fixture, 1.0 / 256.0, false); ASSERT_EQ(evidence.history.size(), 65U); EXPECT_DOUBLE_EQ(evidence.history.back().time, kEndTime); @@ -3307,6 +3320,7 @@ TEST(HhtShellPhysics, TEST(HhtShellPhysics, C08bFreeFiniteRotationPlateSelfConvergenceAndPhysicalInvariants) { const auto fixture = MakeFixture(2U); + UseLegacyDrilling(fixture); std::vector grids; for (const unsigned denominator : {64U, 128U, 256U, 512U, 1024U, 32768U, 65536U, 131072U}) { diff --git a/tests/unit/analysis/nonlinear_analysis_state_test.cpp b/tests/unit/analysis/nonlinear_analysis_state_test.cpp index 4fc9ad3..7c2aa06 100644 --- a/tests/unit/analysis/nonlinear_analysis_state_test.cpp +++ b/tests/unit/analysis/nonlinear_analysis_state_test.cpp @@ -454,3 +454,28 @@ TEST(NonlinearAnalysisState, StoresNoUnapprovedVerificationMetrics) { EXPECT_FALSE((HasVerificationMetrics::value)); EXPECT_FALSE((HasVerificationMetrics::value)); } + +TEST(NonlinearAnalysisState, DrillingEvidenceRetainsModeAndRejectsPartialCommit) { + const auto fixture = MakeFixture(); + auto result = fesa::NonlinearAnalysisState::Create(*fixture.domain, *fixture.dofs, + MakeBaseline(*fixture.domain)); + ASSERT_TRUE(result.HasValue()); + auto& state = result.Value(); + auto first = MakeFrame(state.CommittedState(), 0U, 1U, 0.5, 1U); + first.mitc4_drilling = fesa::Mitc4DrillingEvidence{fesa::Mitc4DrillingMode::kDirectSpring}; + ASSERT_TRUE(state.CommitAcceptedFrame(first).IsOk()); + for (const int kind : {0, 1, 2, 3, 4}) { + auto next = MakeFrame(state.CommittedState(), 1U, 2U, 1.0, 1U); + next.mitc4_drilling = first.mitc4_drilling; + if (kind == 0) next.mitc4_drilling.reset(); + if (kind == 1) next.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) next.mitc4_drilling->mode = static_cast(99); + if (kind == 3) next.mitc4_drilling->equilibrium_with_drilling[3U] = 1.0; + if (kind == 4) next.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::infinity(); + EXPECT_FALSE(state.CommitAcceptedFrame(next).IsOk()) << kind; + ASSERT_EQ(state.AcceptedFrames().size(), 1U); + EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.5); + EXPECT_EQ(state.CommittedState().mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + } +} diff --git a/tests/unit/results/result_recovery_components_test.cpp b/tests/unit/results/result_recovery_components_test.cpp index 1f9a057..8c2785d 100644 --- a/tests/unit/results/result_recovery_components_test.cpp +++ b/tests/unit/results/result_recovery_components_test.cpp @@ -17,9 +17,12 @@ #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" #include "fesa/elements/element.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/mitc4_shell.h" #include "fesa/fem/dof_manager.h" #include "fesa/math/vector3.h" #include "fesa/model/domain.h" +#include "fesa/results/result_recovery.h" #include "results/analysis_state_commit.h" #include "results/beam_result_recovery.h" #include "results/global_equilibrium_recovery.h" @@ -627,3 +630,99 @@ TEST(ResultRecoveryComponents, ExpectStatusCode(status, "nonfinite-recovery-value"); ExpectAnalysisStateEqual(state, prior); } + +TEST(Mitc4DrillingEvidence, LinearDirectSpringAcceptsKnownArtificialTorque) { + auto definition = MakeShellDefinition(); + for (std::size_t node = 0U; node < 4U; ++node) { + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 1, 5, 0.0, {definition.source_path, 51U}}); + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 6, 6, 0.125, {definition.source_path, 52U}}); + } + const auto fixture = MakeShellFixture(std::move(definition)); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; node < 4U; ++node) + state.Displacement()[6U * node + 5U] = 0.125; + const auto status = fesa::ResultRecovery::Recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state); + ASSERT_TRUE(status.IsOk()) << (status.Diagnostics().empty() ? "" : status.Diagnostics()[0U].message); +} + +TEST(Mitc4DrillingEvidence, LinearTorqueOracleAndSignedGateKeepRealImbalance) { + auto definition = MakeShellDefinition(); + for (std::size_t node = 0U; node < 4U; ++node) { + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 1, 5, 0.0, {definition.source_path, 51U}}); + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 6, 6, 0.125, {definition.source_path, 52U}}); + } + const auto fixture = MakeShellFixture(std::move(definition)); + auto element = fesa::ElementFactory{}.Create(fixture.domain->Elements()[0U], + *fixture.domain); + ASSERT_TRUE(element.HasValue()); + const auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + const auto stiffness = shell->Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double rotary_trace = 0.0; + for (std::size_t node = 0U; node < 4U; ++node) + for (std::size_t c = 3U; c < 5U; ++c) + rotary_trace += stiffness.Value().physical_local20(5U * node + c, 5U * node + c); + const double expected_torque = 4.0 * 0.125 * (1e-3 * rotary_trace / 8.0); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; node < 4U; ++node) + state.Displacement()[6U * node + 5U] = 0.125; + ASSERT_TRUE(fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs, + *fixture.stiffness, state).IsOk()); + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_EQ(state.Mitc4Drilling()->mode, fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_NEAR(state.Equilibrium()[5U], expected_torque, 1e-14); + EXPECT_NEAR(state.Mitc4Drilling()->drilling_moment[2U], expected_torque, 1e-14); + EXPECT_DOUBLE_EQ(state.Mitc4Drilling()->drilling_moment[0U], 0.0); + EXPECT_DOUBLE_EQ(state.Mitc4Drilling()->drilling_moment[1U], 0.0); + for (double value : state.Mitc4Drilling()->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-14); + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), 0.0); + + auto global = fesa::results_internal::GlobalEquilibriumRecovery( + *fixture.domain, *fixture.dofs, *fixture.stiffness, state.Displacement(), + state.ExternalForce()); + ASSERT_TRUE(global.HasValue()); + for (const int kind : {0, 1, 2, 3}) { + auto changed = global.Value(); + fesa::results_internal::RecoveryCandidate candidate{24U}; + candidate.shell.mitc4_drilling = *state.Mitc4Drilling(); + if (kind == 0) candidate.shell.mitc4_drilling->drilling_moment[2U] *= -1.0; + if (kind == 1) changed.residual[0U] += 1.0; + if (kind == 2) changed.residual[4U] += 1.0; + if (kind == 3) candidate.shell.mitc4_drilling->mode = + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + const auto rejected = fesa::results_internal::RecoverShellGlobalEvidence( + *fixture.domain, *fixture.dofs, state.ExternalForce(), changed, candidate); + EXPECT_FALSE(rejected.IsOk()) << kind; + } +} + +TEST(Mitc4DrillingEvidence, StaticCommitRejectsInvalidEvidenceAtomically) { + const auto fixture = MakeShellFixture(MakeShellDefinition()); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + fesa::ShellStateCandidate candidate{}; + candidate.rows = MakeShellRows(0U, 1.0).rows; + candidate.mitc4_drilling = fesa::Mitc4DrillingEvidence{ + fesa::Mitc4DrillingMode::kDirectSpring, {0.0, 0.0, 2.0}, {0, 0, 0, 0, 0, -2}}; + ASSERT_TRUE(state.CommitShellResults({0U}, candidate).IsOk()); + for (const int kind : {0, 1, 2, 3, 4}) { + auto invalid = candidate; + invalid.physical_strain_energy = 77.0; + if (kind == 0) invalid.mitc4_drilling->mode = static_cast(99); + if (kind == 1) invalid.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) invalid.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::infinity(); + if (kind == 3) invalid.mitc4_drilling->equilibrium_with_drilling[5U] = 2.0; + if (kind == 4) invalid.mitc4_drilling.reset(); + EXPECT_FALSE(state.CommitShellResults({0U}, invalid).IsOk()) << kind; + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), 0.0); + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_EQ(state.Mitc4Drilling()->mode, candidate.mitc4_drilling->mode); + EXPECT_EQ(state.Mitc4Drilling()->drilling_moment, candidate.mitc4_drilling->drilling_moment); + } +} From df1bf94369e19a3bb1db936bfbca55046afe3e0a Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Tue, 20 Oct 2026 09:27:58 +0900 Subject: [PATCH 2/2] Preserve empty MITC4 hidden-constraint boundary --- src/fesa/analysis/hht_dynamic_algebra.cpp | 20 +++--- .../analysis/hht_dynamic_algebra_test.cpp | 62 ++++++++++++++++++- .../unit/analysis/hht_reduced_shell_test.cpp | 20 ++++++ 3 files changed, 90 insertions(+), 12 deletions(-) diff --git a/src/fesa/analysis/hht_dynamic_algebra.cpp b/src/fesa/analysis/hht_dynamic_algebra.cpp index 7b4c8cd..9311989 100644 --- a/src/fesa/analysis/hht_dynamic_algebra.cpp +++ b/src/fesa/analysis/hht_dynamic_algebra.cpp @@ -64,16 +64,6 @@ bool Within(const double value, const double scale, const double tolerance) { Status Validate(const Domain& domain, const DofManager& dofs, const DynamicAssembly& assembly, const double length) { const std::size_t count = dofs.FullDofCount(); - if (assembly.frame.mitc4_drilling && - assembly.frame.mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) { - if (assembly.hidden_constraint.Size() != assembly.rotations.size()) - return Failure(domain, "dynamic-invalid-algebra-state", - "Canonical hidden values must cover every node."); - for (std::size_t i = 0U; i < assembly.hidden_constraint.Size(); ++i) - if (!std::isfinite(assembly.hidden_constraint[i])) - return Failure(domain, "dynamic-invalid-algebra-state", - "Canonical hidden values must remain finite."); - } if (!std::isfinite(length) || !(length > 0.0) || count % 6U != 0U || assembly.rotations.size() != count / 6U) { return Failure(domain, "dynamic-invalid-algebra-state", @@ -559,6 +549,16 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const bool require_acceleration) { const bool direct = frame->mitc4_drilling && frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring; + if (direct) { + for (const std::size_t node : nodes) { + if (node >= assembly.hidden_constraint.Size()) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must cover every admissible null node."); + if (!std::isfinite(assembly.hidden_constraint[node])) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must remain finite."); + } + } Vector psi{nodes.size()}; Vector psi_scale{nodes.size()}; Vector velocity_defect{nodes.size()}; diff --git a/tests/unit/analysis/hht_dynamic_algebra_test.cpp b/tests/unit/analysis/hht_dynamic_algebra_test.cpp index a6da082..06e1c87 100644 --- a/tests/unit/analysis/hht_dynamic_algebra_test.cpp +++ b/tests/unit/analysis/hht_dynamic_algebra_test.cpp @@ -28,7 +28,8 @@ struct Fixture { Fixture MakeFixture(const std::size_t free_nodes = 1U, const bool partial_rotation = false, const bool translations_free = true, - const double youngs_modulus = 1000.0) { + const double youngs_modulus = 1000.0, + const bool all_rotations_constrained = false) { fesa::ModelDefinition definition{}; definition.source_path = "hht-algebra.inp"; definition.nodes = { @@ -66,7 +67,10 @@ Fixture MakeFixture(const std::size_t free_nodes = 1U, step.boundaries.push_back( {name, 1, 6, 0.0, {definition.source_path, 10U}}); } else { - if (partial_rotation) + if (all_rotations_constrained) + step.boundaries.push_back( + {name, 4, 6, 0.0, {definition.source_path, 10U}}); + else if (partial_rotation) step.boundaries.push_back( {name, 4, 5, 0.0, {definition.source_path, 10U}}); if (!translations_free) @@ -705,3 +709,57 @@ TEST(HhtEndpointFields, MeasuresActualFullChartAcrossNonlinearRootAndPartialBoun } } // namespace + +TEST(HhtDynamicAlgebra, NonDirectionalAssemblyRetainsValidEmptyNullBoundary) { + for (const std::size_t free_nodes : {4U, 0U}) { + SCOPED_TRACE(free_nodes); + const auto fixture = MakeFixture(free_nodes, false, true, 1000.0, true); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, + fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + auto trial = EmptyTrial(fixture); + trial.external_force[0U] = 0.4; + fesa::Vector load_scale(24U); + load_scale[0U] = 0.4; + fesa::SerialParallelFor serial; + const auto parameters = fesa::MakeHhtParameters(-0.05); + ASSERT_TRUE(parameters.HasValue()); + const auto assembly = fesa::hht_internal::EvaluateAssembly( + *fixture.domain, fixture.dofs, shells.Value(), trial, load_scale, + 0.1, parameters.Value(), serial, false); + ASSERT_TRUE(assembly.HasValue()); + ASSERT_TRUE(assembly.Value().frame.mitc4_drilling); + EXPECT_EQ(assembly.Value().frame.mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(assembly.Value().hidden_constraint.Size(), 0U); + EXPECT_FALSE(assembly.Value().hidden_gradient.has_value()); + const auto nodes = fesa::hht_internal::NullNodes(fixture.dofs, + assembly.Value().rotations); + ASSERT_TRUE(nodes.empty()); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + const auto acceleration = fesa::hht_internal::InitializeAcceleration( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, positive, tangent); + EXPECT_TRUE(acceleration.HasValue()) + << (acceleration.GetStatus().Diagnostics().empty() ? "" : + acceleration.GetStatus().Diagnostics()[0U].message); + auto frame = assembly.Value().frame; + if (acceleration.HasValue()) { + frame.acceleration = acceleration.Value(); + // rho*t*area/4 = .2 per node; the .4 force gives acceleration 2. + EXPECT_NEAR(acceleration.Value()[0U], free_nodes == 0U ? 0.0 : 2.0, 1e-12); + for (const auto full : fixture.dofs.ConstrainedDofs()) + EXPECT_DOUBLE_EQ(acceleration.Value()[full], 0.0); + } + const auto complete = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, nodes, + positive, tangent, &frame); + EXPECT_TRUE(complete.IsOk()) + << (complete.Diagnostics().empty() ? "" : complete.Diagnostics()[0U].message); + if (complete.IsOk()) { + EXPECT_EQ(frame.hidden_constraint.Size(), 0U); + EXPECT_EQ(frame.hidden_velocity_defect.Size(), 0U); + EXPECT_EQ(frame.hidden_acceleration_defect.Size(), 0U); + } + } +} diff --git a/tests/unit/analysis/hht_reduced_shell_test.cpp b/tests/unit/analysis/hht_reduced_shell_test.cpp index 4e0d477..10577c7 100644 --- a/tests/unit/analysis/hht_reduced_shell_test.cpp +++ b/tests/unit/analysis/hht_reduced_shell_test.cpp @@ -1938,4 +1938,24 @@ TEST(Mitc4DrillingEvidence, DirectHiddenCertificateUsesCanonicalPrimitive) { EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[row], 0.0); EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[row], 0.0); } + const auto accepted = frame; + for (const bool missing : {true, false}) { + auto invalid = assembly.Value(); + if (missing) + invalid.hidden_constraint = fesa::Vector(0U); + else + invalid.hidden_constraint[0U] = std::numeric_limits::infinity(); + const auto rejected = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, invalid, length.Value(), null_nodes, + positive, tangent, &frame); + EXPECT_FALSE(rejected.IsOk()); + for (const auto& pair : { + std::make_pair(&frame.hidden_constraint, &accepted.hidden_constraint), + std::make_pair(&frame.hidden_velocity_defect, &accepted.hidden_velocity_defect), + std::make_pair(&frame.hidden_acceleration_defect, &accepted.hidden_acceleration_defect)}) { + ASSERT_EQ(pair.first->Size(), pair.second->Size()); + for (std::size_t row = 0U; row < pair.first->Size(); ++row) + ExpectSameBits((*pair.first)[row], (*pair.second)[row]); + } + } }