diff --git a/include/fesa/elements/mitc4_drilling_mode.h b/include/fesa/elements/mitc4_drilling_mode.h new file mode 100644 index 0000000..8df7201 --- /dev/null +++ b/include/fesa/elements/mitc4_drilling_mode.h @@ -0,0 +1,26 @@ +#ifndef FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ +#define FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ + +namespace fesa { + +/// @brief Identifies the immutable numerical MITC4 drilling formulation. +enum class Mitc4DrillingMode { kDirectSpring, kLegacyCorotatedTwist }; + +/// @brief Returns the stable storage name of a supported drilling formulation. +constexpr const char* Mitc4DrillingModeName(Mitc4DrillingMode mode) noexcept { + switch (mode) { + case Mitc4DrillingMode::kDirectSpring: + return "direct-initial-director-spring"; + case Mitc4DrillingMode::kLegacyCorotatedTwist: + return "legacy-corotated-twist"; + } + return ""; +} + +namespace mitc4_drilling_internal { +struct TestAccess; +} // namespace mitc4_drilling_internal + +} // namespace fesa + +#endif // FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ diff --git a/include/fesa/elements/mitc4_shell.h b/include/fesa/elements/mitc4_shell.h index 36e954d..657008e 100644 --- a/include/fesa/elements/mitc4_shell.h +++ b/include/fesa/elements/mitc4_shell.h @@ -8,6 +8,7 @@ #include "fesa/core/status.h" #include "fesa/elements/element.h" #include "fesa/elements/mitc4_shell_definition.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/materials/isotropic_linear_elastic_material.h" #include "fesa/math/matrix.h" #include "fesa/math/vector.h" @@ -92,6 +93,11 @@ class Mitc4Shell final : public Element { std::array, 4> initial_directors, const ShellSection& section, const LinearElasticMaterial& material); + /// @brief Returns this immutable kernel's production drilling formulation. + [[nodiscard]] Mitc4DrillingMode DrillingMode() const noexcept { + return drilling_mode_; + } + /// @brief Returns the factory-bound MITC4 DOF order. const ElementDofLayout& DofLayout() const noexcept override; @@ -163,6 +169,7 @@ class Mitc4Shell final : public Element { private: friend class ElementFactory; + friend struct mitc4_drilling_internal::TestAccess; /// @brief Binds Domain identity after the numerical candidate is valid. void BindRuntime(ElementDofLayout layout, EntityIndex element_index, @@ -203,6 +210,7 @@ class Mitc4Shell final : public Element { double thickness_; double youngs_modulus_; double poisson_ratio_; + Mitc4DrillingMode drilling_mode_; SourceLocation source_location_; std::string identity_; ElementDofLayout dof_layout_; diff --git a/include/fesa/elements/nonlinear_mitc4_shell.h b/include/fesa/elements/nonlinear_mitc4_shell.h index 4359400..e880fdc 100644 --- a/include/fesa/elements/nonlinear_mitc4_shell.h +++ b/include/fesa/elements/nonlinear_mitc4_shell.h @@ -54,6 +54,8 @@ struct NonlinearMitc4InternalEquilibrium { Matrix drilling_metric{0U, 0U}; double drilling_energy = 0.0; Vector force_directional_second{0U}; + /// @brief Stores D2(f_drill)[v,v] separately; direct springs give exact zero. + Vector drilling_force_directional_second{0U}; Vector drilling_constraint{0U}; Matrix drilling_constraint_gradient{0U, 0U}; std::array drilling_constraint_hessian{ @@ -104,6 +106,11 @@ class NonlinearMitc4Shell { std::array initial_frames, const ShellSection& section, const LinearElasticMaterial& material); + /// @brief Returns this immutable kernel's numerical drilling formulation. + [[nodiscard]] Mitc4DrillingMode DrillingMode() const noexcept { + return drilling_mode_; + } + /// @brief Evaluates trial physical kinematics and fixed recovery rows. /// @param trial Global nodal translation and fixed-global rotation vectors. /// @return Physical recovery or a fail-closed trial diagnostic. @@ -134,6 +141,7 @@ class NonlinearMitc4Shell { private: friend struct nonlinear_mitc4_internal::ConfigurationSnapshot; + friend struct mitc4_drilling_internal::TestAccess; /// @brief Retains validated fixed-reference data without analysis ownership. NonlinearMitc4Shell(std::array coordinates, @@ -155,6 +163,7 @@ class NonlinearMitc4Shell { double youngs_modulus_; double poisson_ratio_; double drilling_stiffness_; + Mitc4DrillingMode drilling_mode_; SourceLocation source_location_; std::string identity_; }; diff --git a/src/fesa/elements/mitc4_drilling_internal.h b/src/fesa/elements/mitc4_drilling_internal.h new file mode 100644 index 0000000..820e303 --- /dev/null +++ b/src/fesa/elements/mitc4_drilling_internal.h @@ -0,0 +1,28 @@ +#ifndef FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ +#define FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ + +#include "fesa/elements/nonlinear_mitc4_shell.h" + +namespace fesa::mitc4_drilling_internal { + +// The sole production selection. Change this constant and rebuild all targets. +inline constexpr Mitc4DrillingMode kProductionDrillingMode = + Mitc4DrillingMode::kDirectSpring; + +/// @brief Provides explicit alternate formulations only for internal tests. +/// @note Production factories use the single compile-time selection above. +struct TestAccess { + static Mitc4Shell WithMode(Mitc4Shell shell, Mitc4DrillingMode mode) { + shell.drilling_mode_ = mode; + return shell; + } + static NonlinearMitc4Shell WithMode(NonlinearMitc4Shell shell, + Mitc4DrillingMode mode) { + shell.drilling_mode_ = mode; + return shell; + } +}; + +} // namespace fesa::mitc4_drilling_internal + +#endif // FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp index 6aea36e..9c8c756 100644 --- a/src/fesa/elements/mitc4_shell.cpp +++ b/src/fesa/elements/mitc4_shell.cpp @@ -12,6 +12,7 @@ #include "math/finite_values.h" #include "math/vector3_operations.h" #include "elements/shell_kernel_values.h" +#include "elements/mitc4_drilling_internal.h" namespace fesa { @@ -1002,6 +1003,7 @@ Mitc4Shell::Mitc4Shell(std::array coordinates, thickness_{thickness}, youngs_modulus_{youngs_modulus}, poisson_ratio_{poisson_ratio}, + drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode}, source_location_{std::move(source_location)}, identity_{std::move(identity)} {} diff --git a/src/fesa/elements/nonlinear_mitc4_configuration_internal.h b/src/fesa/elements/nonlinear_mitc4_configuration_internal.h index 618eebc..304c50a 100644 --- a/src/fesa/elements/nonlinear_mitc4_configuration_internal.h +++ b/src/fesa/elements/nonlinear_mitc4_configuration_internal.h @@ -12,6 +12,7 @@ namespace fesa::nonlinear_mitc4_internal { /// key. struct ConfigurationSnapshot { std::array values{}; + Mitc4DrillingMode drilling_mode = Mitc4DrillingMode::kDirectSpring; SourceLocation source{}; std::string identity; diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index d6711cb..e8d41e9 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -12,6 +12,7 @@ #include "fesa/elements/finite_rotation.h" #include "elements/nonlinear_mitc4_configuration_internal.h" +#include "elements/mitc4_drilling_internal.h" #include "math/dense_blas_internal.h" #include "math/jet_arithmetic.h" #include "math/jet_linear_solve.h" @@ -39,6 +40,7 @@ nonlinear_mitc4_internal::ConfigurationSnapshot::Read( result.values[index++] = element.youngs_modulus_; result.values[index++] = element.poisson_ratio_; result.values[index] = element.drilling_stiffness_; + result.drilling_mode = element.drilling_mode_; result.source = element.source_location_; result.identity = element.identity_; return result; @@ -1227,8 +1229,91 @@ Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) { return Reciprocal(Sqrt(Dot(preimage, preimage))); } +/// @brief Uses precisely the current inertia kernel's cubic rotation maps. +/// @note This is deliberately distinct from the legacy analytic null. Keep the +/// cubic Exp/Jacobian order synchronized with EvaluateDynamicRotation; it is +/// not an approximation to be silently replaced by the exact exponential +/// inverse. +JetVector3 InertiaDrillNull(const JetVector3& phi, const Vector3& director) { + const JetMatrix3 skew = Skew(phi); + JetMatrix3 rotation = Identity(); + JetMatrix3 jacobian = Identity(); + JetMatrix3 term = Identity(); + double rotation_factorial = 1.0; + double jacobian_factorial = 1.0; + for (std::size_t order = 1U; order <= 3U; ++order) { + term = Multiply(term, skew); + rotation_factorial *= static_cast(order); + jacobian_factorial *= static_cast(order + 1U); + rotation = Add(rotation, Scale(1.0 / rotation_factorial, term)); + jacobian = Add(jacobian, Scale(1.0 / jacobian_factorial, term)); + } + const JetVector3 rhs = Multiply(rotation, ConstantVector(director)); + Matrix primal{3U, 3U}; + for (std::size_t row = 0U; row < 3U; ++row) + for (std::size_t column = 0U; column < 3U; ++column) + primal(row, column) = jacobian[row][column].value; + const auto factor = dense_blas_internal::DenseFactor::Lu(primal); + const std::array channels{1U, kGlobalDofCount + 1U, + kHessianSize + 1U + kGlobalDofCount, + kGlobalDofCount}; + JetVector3 preimage{}; + if (!factor.HasValue() || + !math_internal::SolveJetAction(jacobian, rhs, factor.Value(), channels, + Constant, ReadChannel, WriteChannel, + preimage)) { + preimage.fill(Constant(std::numeric_limits::quiet_NaN())); + return preimage; + } + return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage); +} + +/// @brief Evaluates the fixed initial-director spring without corotation. +Jet DirectDrillingEnergy( + const std::array& initial_directors, + const JetTrial& trial, double drilling_stiffness, Matrix* metric, + std::array* constraints) { + Jet energy{}; + double* tangent = energy.second.Data(); + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const auto d0 = ConstantVector(initial_directors[node]); + const Jet gamma = Dot(d0, trial.rotation[node]); + // Retain the validated spring scale through subnormal intermediates. + // Fill its linear derivatives explicitly, without differentiating a + // floating half-energy that can round kd before doubling it again. + energy.value += 0.5 * gamma.value * (drilling_stiffness * gamma.value); + energy.directional_first += + drilling_stiffness * gamma.value * gamma.directional_first; + energy.directional_second += + drilling_stiffness * gamma.directional_first * gamma.directional_first; + for (std::size_t i = 0U; i < 3U; ++i) { + const auto row = 6U * node + 3U + i; + energy.first[row] = + drilling_stiffness * gamma.value * initial_directors[node][i]; + energy.mixed_first[row] = drilling_stiffness * gamma.directional_first * + initial_directors[node][i]; + } + // D2 f_drill[v,v] is identically zero; mixed_second retains exact zeros. + if (constraints != nullptr) { + const auto null = + InertiaDrillNull(trial.rotation[node], initial_directors[node]); + (*constraints)[node] = drilling_stiffness * gamma * Dot(null, d0); + } + for (std::size_t i = 0U; i < 3U; ++i) + for (std::size_t j = 0U; j < 3U; ++j) { + const auto row = 6U * node + 3U + i; + const auto column = 6U * node + 3U + j; + const double value = drilling_stiffness * initial_directors[node][i] * + initial_directors[node][j]; + tangent[row * kGlobalDofCount + column] = value; + (*metric)(row, column) = value; + } + } + return energy; +} + /// @brief Evaluates the objective numerical drilling potential and variations. -Jet DrillingEnergy(const std::array& coordinates, +Jet LegacyDrillingEnergy(const std::array& coordinates, const std::array& initial_directors, const JetTrial& trial, const double thickness, const double drilling_stiffness, bool& valid, Matrix* metric, @@ -1323,6 +1408,7 @@ NonlinearMitc4Shell::NonlinearMitc4Shell( youngs_modulus_{youngs_modulus}, poisson_ratio_{poisson_ratio}, drilling_stiffness_{drilling_stiffness}, + drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode}, source_location_{std::move(source_location)}, identity_{std::move(identity)} {} @@ -1640,7 +1726,10 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint( jet_phi[component] = Variable(phi[component], component, phi_dot[component]); } - const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]); + const JetVector3 null_direction = + drilling_mode_ == Mitc4DrillingMode::kDirectSpring + ? InertiaDrillNull(jet_phi, initial_directors_[node]) + : DrillNull(jet_phi, initial_directors_[node]); const Jet scalar = -1.0 * Dot(null_direction, ConstantVector(force)); const Jet scalar_time = -1.0 * Dot(null_direction, ConstantVector(force_rate)); @@ -1741,10 +1830,17 @@ Result EvaluateInternalResidualAndTangent( element.poisson_ratio_, valid); Matrix drilling_metric{kGlobalDofCount, kGlobalDofCount}; std::array constraints{}; - const Jet drilling = DrillingEnergy( - element.coordinates_, element.initial_directors_, jet_trial, - element.thickness_, element.drilling_stiffness_, valid, &drilling_metric, - coordinate_direction != nullptr ? &constraints : nullptr); + const Jet drilling = + element.drilling_mode_ == Mitc4DrillingMode::kDirectSpring + ? DirectDrillingEnergy( + element.initial_directors_, jet_trial, + element.drilling_stiffness_, &drilling_metric, + coordinate_direction != nullptr ? &constraints : nullptr) + : LegacyDrillingEnergy( + element.coordinates_, element.initial_directors_, jet_trial, + element.thickness_, element.drilling_stiffness_, valid, + &drilling_metric, + coordinate_direction != nullptr ? &constraints : nullptr); const Jet internal_energy = physical + drilling; if (!valid || !IsFinite(internal_energy)) { return Failure( @@ -1815,10 +1911,13 @@ Result EvaluateInternalResidualAndTangent( equilibrium.drilling_metric = std::move(drilling_metric); equilibrium.drilling_energy = drilling.value; equilibrium.force_directional_second = Vector{kGlobalDofCount}; + equilibrium.drilling_force_directional_second = Vector{kGlobalDofCount}; for (std::size_t index = 0U; index < kGlobalDofCount; ++index) { equilibrium.internal_force[index] = equilibrium.internal_residual[index]; equilibrium.physical_force[index] = physical.first[index]; equilibrium.drilling_force[index] = drilling.first[index]; + equilibrium.drilling_force_directional_second[index] = + drilling.mixed_second[index]; equilibrium.force_directional_second[index] = internal_energy.mixed_second[index]; for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { diff --git a/tests/unit/elements/element_capabilities_test.cpp b/tests/unit/elements/element_capabilities_test.cpp index 585d749..7e6b556 100644 --- a/tests/unit/elements/element_capabilities_test.cpp +++ b/tests/unit/elements/element_capabilities_test.cpp @@ -133,6 +133,7 @@ TEST(ElementCapabilities, ShellTrialPreservesCompleteDirectionalResponseAndFailu Same(a.drilling_metric, b.drilling_metric); Same(a.drilling_energy, b.drilling_energy); Same(a.force_directional_second, b.force_directional_second); + Same(a.drilling_force_directional_second, b.drilling_force_directional_second); Same(a.drilling_constraint, b.drilling_constraint); Same(a.drilling_constraint_gradient, b.drilling_constraint_gradient); for (std::size_t n = 0; n < 4; ++n) { diff --git a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp index 5256a2b..20bcb00 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp +++ b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -13,6 +14,8 @@ #include "fesa/elements/nonlinear_mitc4_inertia.h" #include "fesa/elements/nonlinear_mitc4_shell.h" #include "fesa/model/model_types.h" +#include "elements/mitc4_drilling_internal.h" +#include "elements/nonlinear_mitc4_configuration_internal.h" namespace { @@ -44,7 +47,7 @@ std::array Frames() { {3U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}}}}; } -fesa::NonlinearMitc4Shell Shell(const std::array& nodes) { +fesa::NonlinearMitc4Shell ProductionShell(const std::array& nodes) { const fesa::ShellSection section{ "Section", 0.2, 0U, {"nonlinear-mitc4-tangent.inp", 20U}}; const fesa::LinearElasticMaterial material{ @@ -57,6 +60,11 @@ fesa::NonlinearMitc4Shell Shell(const std::array& nodes) { return shell.Value(); } +fesa::NonlinearMitc4Shell LegacyShell(const std::array& nodes) { + return fesa::mitc4_drilling_internal::TestAccess::WithMode( + ProductionShell(nodes), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + fesa::NonlinearMitc4TrialDofs Trial() { fesa::NonlinearMitc4TrialDofs trial{}; trial.displacement[0U] = {-0.02, 0.01, 0.03}; @@ -200,7 +208,7 @@ void ExpectRelative(const double actual, const double expected) { // NSNR-ELEM-TAN-001 TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto recovery = shell.RecoverState(trial); const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial); @@ -230,7 +238,7 @@ TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) { // NSNR-ELEM-TAN-002 and NSNR-ELEM-TAN-003 TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto direction = WorkConjugateDirection(nodes); const double length = ReferenceLength(nodes); @@ -295,7 +303,7 @@ TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) { TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { const auto nodes = Nodes(); const auto equilibrium = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), {}); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), {}); ASSERT_TRUE(equilibrium.HasValue()); const double length = ReferenceLength(nodes); @@ -326,9 +334,9 @@ TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { TEST(NonlinearMitc4Residual, DoesNotAcceptExternalLoadShares) { const auto nodes = Nodes(); const auto first = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial()); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial()); const auto second = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial()); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial()); ASSERT_TRUE(first.HasValue()); ASSERT_TRUE(second.HasValue()); for (std::size_t index = 0U; index < 24U; ++index) { @@ -357,7 +365,7 @@ TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) { trial.rotation_vector[node] = rotation; } const auto result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial); ASSERT_TRUE(result.HasValue()); EXPECT_LE(result.Value().internal_residual.Norm(), kNumericalTolerance); EXPECT_LE(result.Value().physical_strain_energy, kNumericalTolerance); @@ -366,7 +374,7 @@ TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) { // NSNR-ELEM-TAN-008 TEST(NonlinearMitc4Patch, RetainsPhysicalRecoveryIdentities) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); constexpr double kMembraneStrain = 1.0e-2; const double green_strain = kMembraneStrain + 0.5 * kMembraneStrain * kMembraneStrain; @@ -409,7 +417,7 @@ TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) { trial.rotation_vector[2U] = {std::numeric_limits::quiet_NaN(), 0.0, 0.0}; const auto result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial); ASSERT_FALSE(result.HasValue()); EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code, "invalid-nonlinear-shell-trial"); @@ -417,13 +425,13 @@ TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) { auto overflow = Trial(); overflow.displacement[1U] = {1.0e300, 0.0, 0.0}; const auto overflow_result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), overflow); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), overflow); EXPECT_FALSE(overflow_result.HasValue()); } void CheckThirdSweep(const fesa::NonlinearMitc4TrialDofs& trial) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); auto direction = Direction(); double norm = 0.0; const auto unscaled = ToVector(direction); @@ -488,7 +496,7 @@ TEST(NonlinearMitc4Residual, nodes[1].coordinates = {4.23076923, 2.30769231, 0.0}; nodes[2].coordinates = {4.23076923, 2.69230769, 0.0}; nodes[3].coordinates = {3.84615385, 2.69230769, 0.0}; - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto evaluated = fesa::EvaluateInternalResidualAndTangent(shell, {}); ASSERT_TRUE(evaluated.HasValue()); for (std::size_t i = 0; i < 24U; ++i) @@ -566,7 +574,7 @@ TEST(NonlinearMitc4Residual, for (bool warped : {false, true}) { auto nodes = Nodes(); if (warped) nodes[2U].coordinates[2U] = 0.1; - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); auto trial = Trial(); trial.rotation_vector[0U] = {.6, -.2, .1}; const auto direction = Direction(); @@ -636,7 +644,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExposesExactDrillingConstraintThirdDirectionalGradient) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto direction = Direction(); const auto velocity = ToVector(direction); @@ -679,7 +687,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExternalDrillingConstraintIncludesForceTimeAndCurvatureTerms) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const fesa::Vector3 phi{.6, -.2, .1}; const fesa::Vector3 phi_dot{.2, -.1, .3}; const fesa::Vector3 force{1.7, -.8, .4}; @@ -742,7 +750,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExternalDrillingConstraintRejectsInvalidRotationAndHasZeroForceValue) { - const auto shell = Shell(Nodes()); + const auto shell = LegacyShell(Nodes()); const auto zero = shell.EvaluateExternalDrillingConstraint( 0U, {.2, -.1, .3}, {.1, .2, -.3}, {}, {}, {}); ASSERT_TRUE(zero.HasValue()); @@ -838,4 +846,470 @@ TEST(NonlinearMitc4Residual, TinyRotationPreservesPhysicalForceAndRecoveryEnergy } } +// Independent fixed-initial-director spring oracle: a total-chart rigid spin +// carries numerical energy even though the physical shell is undeformed. +TEST(Mitc4DrillingSelection, ProductionUsesDirectInitialDirectorSpring) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + const fesa::ShellSection section{"Section", 0.2, 0U, {"spring.inp", 1U}}; + const fesa::LinearElasticMaterial material{ + "Material", 120.0, 0.25, {"spring.inp", 2U}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + const auto stiffness = linear.Value().Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double rotary_trace = 0.0; + for (std::size_t node = 0; node < 4; ++node) + for (std::size_t c = 3; c < 5; ++c) + rotary_trace += + stiffness.Value().physical_local20(5 * node + c, 5 * node + c); + const double kd = 1e-3 * rotary_trace / 8.; + fesa::NonlinearMitc4TrialDofs trial{}; + const fesa::Vector3 phi{0., 0., 0.4}; + const auto rotation = fesa::ExpSkew(phi); + for (std::size_t node = 0; node < 4; ++node) { + const fesa::Vector3 x{nodes[node].coordinates}; + trial.rotation_vector[node] = phi; + fesa::Vector3 rotated{ + rotation(0, 0) * x[0] + rotation(0, 1) * x[1] + rotation(0, 2) * x[2], + rotation(1, 0) * x[0] + rotation(1, 1) * x[1] + rotation(1, 2) * x[2], + rotation(2, 0) * x[0] + rotation(2, 1) * x[1] + rotation(2, 2) * x[2]}; + trial.displacement[node] = rotated - x; + } + const auto actual = fesa::EvaluateInternalResidualAndTangent(shell, trial); + ASSERT_TRUE(actual.HasValue()); + EXPECT_NEAR(actual.Value().physical_strain_energy, 0., 1e-25); + EXPECT_NEAR(actual.Value().drilling_energy, .5 * kd * 4 * .4 * .4, 1e-14); + for (std::size_t row = 0; row < 24; ++row) { + const double force = row % 6 == 5 ? kd * .4 : 0.; + EXPECT_NEAR(actual.Value().drilling_force[row], force, 1e-14); + for (std::size_t col = 0; col < 24; ++col) { + const double tangent = row == col && row % 6 == 5 ? kd : 0.; + EXPECT_NEAR(actual.Value().drilling_tangent(row, col), tangent, 1e-14); + } + } +} + +TEST(Mitc4DrillingSelection, + DirectSpringUsesEachInitialDirectorAndConstantBlocks) { + for (const bool warped : {false, true}) { + auto nodes = Nodes(); + if (warped) nodes[2].coordinates[2] = .1; + auto frames = Frames(); + const auto transform = [](const fesa::Matrix& r, const Coordinates& x) { + Coordinates y{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) y[i] += r(i, j) * x[j]; + return y; + }; + const auto global = fesa::ExpSkew({.4, -.3, .2}); + std::array directors{}; + for (std::size_t n = 0; n < 4; ++n) { + const auto r = + fesa::ExpSkew({0., warped ? .03 * static_cast(n) : 0., 0.}); + frames[n].director = transform(global, transform(r, frames[n].director)); + frames[n].tangent_a = + transform(global, transform(r, frames[n].tangent_a)); + frames[n].tangent_b = + transform(global, transform(r, frames[n].tangent_b)); + nodes[n].coordinates = transform(global, nodes[n].coordinates); + directors[n] = frames[n].director; + } + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + const fesa::LinearElasticMaterial material{ + "Material", 120., .25, {"spring.inp", 2}}; + const auto shell = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), frames, section, material); + const auto linear = fesa::Mitc4Shell::Create(NodePointers(nodes), directors, + section, material); + ASSERT_TRUE(shell.HasValue()); + ASSERT_TRUE(linear.HasValue()); + const auto k = linear.Value().Stiffness(); + ASSERT_TRUE(k.HasValue()); + double trace = 0.; + for (std::size_t n = 0; n < 4; ++n) + for (std::size_t c = 3; c < 5; ++c) + trace += k.Value().physical_local20(5 * n + c, 5 * n + c); + const double kd = 1e-3 * trace / 8.; + auto trial = Trial(); + trial.rotation_vector[0] = {.65, -.4, .3}; + const auto direction = ToVector(Direction()); + const auto response = fesa::EvaluateInternalResidualAndTangent( + shell.Value(), trial, &direction); + const auto shifted = fesa::EvaluateInternalResidualAndTangent( + shell.Value(), AddScaled(trial, Direction(), .2), &direction); + ASSERT_TRUE(response.HasValue()); + ASSERT_TRUE(shifted.HasValue()); + double energy = 0.; + for (std::size_t n = 0; n < 4; ++n) { + const double gamma = + fesa::Vector3{directors[n]}.Dot(trial.rotation_vector[n]); + energy += .5 * kd * gamma * gamma; + const fesa::Vector3 d{directors[n]}; + const auto rotation = fesa::EvaluateDynamicRotation( + trial.rotation_vector[n], d, Direction().rotation_vector[n], {}); + ASSERT_TRUE(rotation.HasValue()); + EXPECT_NEAR(response.Value().drilling_constraint[n], + kd * gamma * rotation.Value().drill_null.Dot(d), 1e-14); + for (std::size_t c = 0; c < 3; ++c) { + double ds = 0.; + for (std::size_t j = 0; j < 3; ++j) + ds += d[j] * rotation.Value().null_gradient(j, c); + EXPECT_NEAR( + response.Value().drilling_constraint_gradient(n, 6 * n + 3 + c), + kd * (d[c] * rotation.Value().drill_null.Dot(d) + gamma * ds), + 1e-14); + } + for (std::size_t c = 0; c < 3; ++c) + EXPECT_NEAR(response.Value().drilling_force[6 * n + 3 + c], + kd * gamma * directors[n][c], 1e-14); + } + EXPECT_NEAR(response.Value().drilling_energy, energy, 1e-14); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_DOUBLE_EQ(response.Value().drilling_force_directional_second[i], + 0.); + if (i % 6 < 3) EXPECT_DOUBLE_EQ(response.Value().drilling_force[i], 0.); + for (std::size_t j = 0; j < 24; ++j) { + const double expected = + i / 6 == j / 6 && i % 6 >= 3 && j % 6 >= 3 + ? kd * directors[i / 6][i % 6 - 3] * directors[i / 6][j % 6 - 3] + : 0.; + EXPECT_NEAR(response.Value().drilling_tangent(i, j), expected, 1e-14); + EXPECT_NEAR(response.Value().drilling_metric(i, j), expected, 1e-14); + EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j), + shifted.Value().drilling_tangent(i, j)); + if (expected == 0.) + EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j), 0.); + } + } + } +} + +TEST(Mitc4DrillingSelection, + DirectConstraintUsesPhysicalNullAndItsHighDerivatives) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + auto trial = Trial(); + trial.rotation_vector[0] = {.65, -.4, .3}; + const auto direction = Direction(); + const auto velocity = ToVector(direction); + const auto response = + fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity); + ASSERT_TRUE(response.HasValue()); + const auto zero = + fesa::EvaluateInternalResidualAndTangent(shell, {}, &velocity); + ASSERT_TRUE(zero.HasValue()); + const double kd = zero.Value().drilling_tangent(5, 5); + // Independent inertia API provides n, Dn and D2n, never a twist angle. + const fesa::Vector3 d{0., 0., 1.}; + const auto scalar_gradient = [&](const fesa::Vector3& phi, std::size_t c) { + const auto r = fesa::EvaluateDynamicRotation(phi, d, {}, {}); + EXPECT_TRUE(r.HasValue()); + double projection_derivative = 0.; + for (std::size_t k = 0; k < 3; ++k) + projection_derivative += r.Value().null_gradient(k, c) * d[k]; + return kd * (d[c] * r.Value().drill_null.Dot(d) + + d.Dot(phi) * projection_derivative); + }; + constexpr double h = .002; + for (std::size_t n = 0; n < 4; ++n) { + const auto phi = trial.rotation_vector[n]; + const auto v = direction.rotation_vector[n]; + const auto r = fesa::EvaluateDynamicRotation(phi, d, v, {}); + ASSERT_TRUE(r.HasValue()); + const double gamma = d.Dot(phi); + const double gamma_dot = d.Dot(v); + double projection_dot = 0.; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + projection_dot += d[i] * r.Value().null_gradient(i, j) * v[j]; + const double expected_second = + kd * (2. * gamma_dot * projection_dot + + gamma * r.Value().null_directional_second.Dot(d)); + double actual_second = 0.; + EXPECT_NEAR(response.Value().drilling_constraint[n], + kd * gamma * r.Value().drill_null.Dot(d), 1e-14); + for (std::size_t i = 0; i < 24; ++i) { + const double gradient = + i / 6 == n && i % 6 >= 3 ? scalar_gradient(phi, i % 6 - 3) : 0.; + EXPECT_NEAR(response.Value().drilling_constraint_gradient(n, i), gradient, + 1e-14); + for (std::size_t j = 0; j < 24; ++j) + actual_second += velocity[i] * + response.Value().drilling_constraint_hessian[n](i, j) * + velocity[j]; + if (i / 6 == n && i % 6 >= 3) { + const auto c = i % 6 - 3; + const double fd = + (scalar_gradient(phi + h * v, c) - 2. * scalar_gradient(phi, c) + + scalar_gradient(phi - h * v, c)) / + (h * h); + const double actual = + response.Value() + .drilling_constraint_third_directional_gradient[n][i]; + EXPECT_NEAR(actual, fd, + 2e-5 * (std::abs(actual) + std::abs(fd) + kd * 1e-4)); + } else { + EXPECT_DOUBLE_EQ( + response.Value() + .drilling_constraint_third_directional_gradient[n][i], + 0.); + } + } + EXPECT_NEAR(actual_second, expected_second, 1e-14); + } +} + +TEST(Mitc4DrillingSelection, DirectExternalConstraintUsesActualInertiaNull) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + const fesa::Vector3 phi{.65, -.4, .3}, velocity{.2, .1, -.15}; + const fesa::Vector3 force{.7, -.3, .2}, rate{-.1, .2, .4}, + acceleration{.5, -.2, .1}; + const auto actual = shell.EvaluateExternalDrillingConstraint( + 0, phi, velocity, force, rate, acceleration); + const auto rotation = + fesa::EvaluateDynamicRotation(phi, {0., 0., 1.}, velocity, {}); + ASSERT_TRUE(actual.HasValue()); + ASSERT_TRUE(rotation.HasValue()); + const double h = .002; + const auto plus = shell.EvaluateExternalDrillingConstraint( + 0, phi + h * velocity, velocity, + force + h * rate + .5 * h * h * acceleration, rate + h * acceleration, + acceleration); + const auto minus = shell.EvaluateExternalDrillingConstraint( + 0, phi - h * velocity, velocity, + force - h * rate + .5 * h * h * acceleration, rate - h * acceleration, + acceleration); + ASSERT_TRUE(plus.HasValue()); + ASSERT_TRUE(minus.HasValue()); + for (std::size_t c = 0; c < 3; ++c) { + EXPECT_NEAR( + actual.Value().gradient_directional_first[c], + (plus.Value().gradient[c] - minus.Value().gradient[c]) / (2. * h), + 1e-7); + EXPECT_NEAR(actual.Value().gradient_directional_second[c], + (plus.Value().gradient[c] - 2. * actual.Value().gradient[c] + + minus.Value().gradient[c]) / + (h * h), + 2e-7); + std::array basis{}; + basis[c] = 1.; + constexpr double coordinate_h = 1e-4; + const auto positive = shell.EvaluateExternalDrillingConstraint( + 0, phi + coordinate_h * fesa::Vector3{basis}, velocity, force, rate, + acceleration); + const auto negative = shell.EvaluateExternalDrillingConstraint( + 0, phi - coordinate_h * fesa::Vector3{basis}, velocity, force, rate, + acceleration); + ASSERT_TRUE(positive.HasValue()); + ASSERT_TRUE(negative.HasValue()); + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR( + actual.Value().hessian(j, c), + (positive.Value().gradient[j] - negative.Value().gradient[j]) / + (2. * coordinate_h), + 1e-7); + } + const auto& n = rotation.Value(); + EXPECT_NEAR(actual.Value().value, -n.drill_null.Dot(force), 1e-14); + fesa::Vector3 n_dot{}; + std::array values{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + values[i] += n.null_gradient(i, j) * velocity[j]; + n_dot = fesa::Vector3{values}; + EXPECT_NEAR(actual.Value().directional_first, + -n_dot.Dot(force) - n.drill_null.Dot(rate), 1e-14); + EXPECT_NEAR(actual.Value().directional_second, + -n.null_directional_second.Dot(force) - 2. * n_dot.Dot(rate) - + n.drill_null.Dot(acceleration), + 1e-14); + for (std::size_t c = 0; c < 3; ++c) { + double expected = 0.; + for (std::size_t i = 0; i < 3; ++i) + expected -= n.null_gradient(i, c) * force[i]; + EXPECT_NEAR(actual.Value().gradient[c], expected, 1e-14); + } +} + +TEST(Mitc4DrillingSelection, + ModesAreImmutableReferenceInputsAndPreservePhysicalResponse) { + const auto nodes = Nodes(); + const auto direct = ProductionShell(nodes); + const auto legacy = LegacyShell(nodes); + EXPECT_EQ(direct.DrillingMode(), fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(legacy.DrillingMode(), + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + EXPECT_STREQ(fesa::Mitc4DrillingModeName(direct.DrillingMode()), + "direct-initial-director-spring"); + EXPECT_STREQ(fesa::Mitc4DrillingModeName(legacy.DrillingMode()), + "legacy-corotated-twist"); + const auto a = + fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(direct); + const auto b = + fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(legacy); + EXPECT_EQ(a.values, b.values); + EXPECT_NE(a.drilling_mode, b.drilling_mode); + const auto direction = ToVector(Direction()); + const auto r = + fesa::EvaluateInternalResidualAndTangent(direct, Trial(), &direction); + const auto l = + fesa::EvaluateInternalResidualAndTangent(legacy, Trial(), &direction); + ASSERT_TRUE(r.HasValue()); + ASSERT_TRUE(l.HasValue()); + EXPECT_DOUBLE_EQ(r.Value().physical_strain_energy, + l.Value().physical_strain_energy); + EXPECT_NE(r.Value().drilling_energy, l.Value().drilling_energy); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_DOUBLE_EQ(r.Value().physical_force[i], l.Value().physical_force[i]); + for (std::size_t j = 0; j < 24; ++j) + EXPECT_NEAR( + r.Value().internal_tangent(i, j) - r.Value().drilling_tangent(i, j), + l.Value().internal_tangent(i, j) - l.Value().drilling_tangent(i, j), + 1e-14); + } + for (std::size_t gp = 0; gp < 4; ++gp) + EXPECT_EQ(r.Value().recovery.gauss_points[gp].section_resultant, + l.Value().recovery.gauss_points[gp].section_resultant); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + const fesa::LinearElasticMaterial material{ + "Material", 120., .25, {"spring.inp", 2}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + EXPECT_EQ(linear.Value().DrillingMode(), direct.DrillingMode()); + const auto alternate = fesa::mitc4_drilling_internal::TestAccess::WithMode( + linear.Value(), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + const auto k1 = linear.Value().Stiffness(); + const auto k2 = alternate.Stiffness(); + ASSERT_TRUE(k1.HasValue()); + ASSERT_TRUE(k2.HasValue()); + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + EXPECT_DOUBLE_EQ(k1.Value().stabilized_global24(i, j), + k2.Value().stabilized_global24(i, j)); +} + +TEST(Mitc4DrillingSelection, DirectRejectsInvalidScaleAndNonfiniteTrial) { + const auto nodes = Nodes(); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + for (double youngs_modulus : {1e-320, (std::numeric_limits::max)()}) { + const fesa::LinearElasticMaterial material{ + "Material", youngs_modulus, .25, {"spring.inp", 2}}; + const auto invalid = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), Frames(), section, material); + ASSERT_FALSE(invalid.HasValue()); + EXPECT_EQ(invalid.GetStatus().Diagnostics()[0].code, + youngs_modulus < 1.0 ? "invalid-nonlinear-shell-drilling" + : "invalid-shell-material"); + } + const auto shell = ProductionShell(nodes); + for (double value : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + auto trial = Trial(); + trial.rotation_vector[0] = {value, 0., 0.}; + EXPECT_FALSE( + fesa::EvaluateInternalResidualAndTangent(shell, trial).HasValue()); + } + const auto valid = fesa::EvaluateInternalResidualAndTangent(shell, {}); + ASSERT_TRUE(valid.HasValue()); + EXPECT_DOUBLE_EQ(valid.Value().drilling_energy, 0.); +} + +TEST(Mitc4DrillingSelection, DirectRetainsPositiveSubnormalSpringScale) { + const auto nodes = Nodes(); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + unsigned accepted = 0; + for (double youngs_modulus : + {1e-315, 1e-316, 1e-317, 1e-318, 2e-318, 3e-318, 1e-319}) { + const fesa::LinearElasticMaterial material{ + "Material", youngs_modulus, .25, {"spring.inp", 2}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + const auto stiffness = linear.Value().Stiffness(); + const auto shell = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), Frames(), section, material); + if (!stiffness.HasValue()) { + EXPECT_FALSE(shell.HasValue()); + continue; + } + ASSERT_TRUE(shell.HasValue()); + ++accepted; + fesa::NonlinearMitc4TrialDofs trial{}; + trial.rotation_vector[0] = {0., 0., 1.}; + const auto actual = + fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial); + ASSERT_TRUE(actual.HasValue()); + // Scale is already validated. No later half-scale rounding may change K. + const double kd = stiffness.Value().drilling_stiffness; + EXPECT_EQ(actual.Value().drilling_tangent(5, 5), kd); + EXPECT_EQ(actual.Value().drilling_force[5], kd); + trial.rotation_vector[0] = {0., 0., 2.}; + const auto double_angle = + fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial); + ASSERT_TRUE(double_angle.HasValue()); + // U(2) = 2*kd is representable even for the least positive spring scale. + EXPECT_EQ(double_angle.Value().drilling_energy, 2. * kd); + EXPECT_EQ(double_angle.Value().drilling_force[5], 2. * kd); + } + EXPECT_GT(accepted, 0U); +} + +TEST(Mitc4DrillingSelection, LegacyHeadFrozenEvidence) { + const auto nodes = Nodes(); + const auto shell = LegacyShell(nodes); + const auto direction = ToVector(Direction()); + const auto response = + fesa::EvaluateInternalResidualAndTangent(shell, Trial(), &direction); + ASSERT_TRUE(response.HasValue()); + const auto& r = response.Value(); + std::array hashes{}; + hashes.fill(14695981039346656037ULL); + const auto add = [&](std::size_t channel, double value) { + std::uint64_t bits{}; + std::memcpy(&bits, &value, sizeof(bits)); + for (unsigned b = 0; b < 8; ++b) { + hashes[channel] ^= (bits >> (8 * b)) & 255ULL; + hashes[channel] *= 1099511628211ULL; + } + }; + add(0, r.drilling_energy); + for (std::size_t i = 0; i < 24; ++i) { + add(1, r.drilling_force[i]); + add(4, r.force_directional_second[i]); + for (std::size_t j = 0; j < 24; ++j) { + add(2, r.drilling_tangent(i, j)); + add(3, r.drilling_metric(i, j)); + } + } + for (std::size_t n = 0; n < 4; ++n) { + add(5, r.drilling_constraint[n]); + for (std::size_t i = 0; i < 24; ++i) { + add(5, r.drilling_constraint_gradient(n, i)); + add(7, r.drilling_constraint_third_directional_gradient[n][i]); + for (std::size_t j = 0; j < 24; ++j) + add(6, r.drilling_constraint_hessian[n](i, j)); + } + } + // Captured before implementation from HEAD + // 52e321e15c70377d4318aace3f3bc64b3ea1a65c, native MSVC x64 Debug. Entire + // energy/force/tangent/metric/directional/psi inventories are pinned + // independently, including the historical polar1 constraint inconsistency. Do + // not regenerate these from the selected mode. + constexpr std::array frozen{ + 16395275265220000061ULL, 12405324468590059764ULL, 79034762627498281ULL, + 7305417493946306304ULL, 7191018337857616990ULL, 4845391414006668707ULL, + 11613163243225266145ULL, 3872161451216566265ULL}; + EXPECT_EQ(hashes, frozen); +} + } // namespace