From 1a6a9bfc715fc321907a6229a84bb4ef8747558a Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 12 Aug 2026 19:54:27 +0900 Subject: [PATCH] feat(linear-static-mitc4-shell): step 4 - mitc4-stiffness-drilling --- include/fesa/elements/mitc4_shell.hpp | 20 +- src/fesa/elements/mitc4_shell.cpp | 164 ++++++++- tests/unit/elements/mitc4_shell_test.cpp | 442 +++++++++++++++++++++++ 3 files changed, 620 insertions(+), 6 deletions(-) diff --git a/include/fesa/elements/mitc4_shell.hpp b/include/fesa/elements/mitc4_shell.hpp index c8593a4..8ba0db9 100644 --- a/include/fesa/elements/mitc4_shell.hpp +++ b/include/fesa/elements/mitc4_shell.hpp @@ -5,6 +5,7 @@ #include "fesa/model/model_types.hpp" #include +#include namespace fesa { @@ -30,8 +31,16 @@ struct Mitc4QuadraturePoint { double weight; }; -// Concrete small-rotation MITC4 kinematics and material value kernel. Global -// equation ownership, drilling stiffness, recovery, and assembly remain outside. +struct Mitc4Stiffness { + Matrix physicalLocal20; + Matrix physicalGlobal24; + Matrix drillingGlobal24; + Matrix stabilizedGlobal24; + double drillingStiffness; +}; + +// Concrete small-rotation MITC4 kinematics, constitutive, and stiffness kernel. +// Global equation ownership, physical recovery, and assembly remain outside. class Mitc4Shell { public: static Result create( @@ -63,6 +72,7 @@ public: [[nodiscard]] Matrix membraneSectionMatrix() const; [[nodiscard]] Matrix bendingSectionMatrix() const; [[nodiscard]] Matrix transverseShearSectionMatrix() const; + [[nodiscard]] Result stiffness() const; private: using Vector3 = std::array; @@ -82,7 +92,9 @@ private: Vector3 normalCandidate, double thickness, double youngsModulus, - double poissonRatio); + double poissonRatio, + SourceLocation sourceLocation, + std::string identity); bool evaluateGeometry( double xi, @@ -107,6 +119,8 @@ private: double thickness_; double youngsModulus_; double poissonRatio_; + SourceLocation sourceLocation_; + std::string identity_; }; } // namespace fesa diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp index 2f6c843..45e92f1 100644 --- a/src/fesa/elements/mitc4_shell.cpp +++ b/src/fesa/elements/mitc4_shell.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -213,6 +214,55 @@ Matrix scaledMatrix(const Matrix& source, double factor) { return result; } +Matrix congruence(const Matrix& local, const Matrix& transformation) { + if (local.rows() != local.columns() || + local.rows() != transformation.rows()) { + throw std::invalid_argument{"MITC4 congruence dimensions are incompatible."}; + } + Matrix result{transformation.columns(), transformation.columns()}; + for (std::size_t row = 0U; row < result.rows(); ++row) { + for (std::size_t column = row; column < result.columns(); ++column) { + double value = 0.0; + for (std::size_t localRow = 0U; localRow < local.rows(); ++localRow) { + for (std::size_t localColumn = 0U; + localColumn < local.columns(); ++localColumn) { + value += transformation(localRow, row) * + local(localRow, localColumn) * + transformation(localColumn, column); + } + } + result(row, column) = value; + result(column, row) = value; + } + } + return result; +} + +bool isFinite(const Matrix& matrix) { + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + for (std::size_t column = 0U; column < matrix.columns(); ++column) { + if (!std::isfinite(matrix(row, column))) { + return false; + } + } + } + return true; +} + +Result stiffnessFailure( + const SourceLocation& location, + const std::string& identity, + std::string message) { + return Result::failure(Status::failure( + FailureCategory::model, + {{Severity::error, + "invalid-shell-stiffness", + location, + "*ELEMENT", + identity, + std::move(message)}})); +} + } // namespace Result Mitc4Shell::create( @@ -302,7 +352,9 @@ Result Mitc4Shell::create( normalCandidate, section.thickness, material.youngsModulus, - material.poissonRatio}; + material.poissonRatio, + nodes[0]->location, + identity}; const double shearModulus = material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); @@ -519,6 +571,108 @@ Matrix Mitc4Shell::transverseShearSectionMatrix() const { return result; } +Result Mitc4Shell::stiffness() const { + const Matrix constitutive = materialConstitutive5(); + const Matrix tyingSamples = covariantTyingShearSamples20(); + Matrix physicalLocal{kPhysicalDofCount, kPhysicalDofCount}; + for (const auto& point : volumeQuadrature()) { + GeometryData geometry{}; + if (!evaluateGeometry( + point.naturalCoordinates[0], + point.naturalCoordinates[1], + point.naturalCoordinates[2], + geometry)) { + return stiffnessFailure( + sourceLocation_, identity_, + "Validated MITC4 quadrature geometry became invalid."); + } + const Matrix strain = strainDisplacement( + point.naturalCoordinates[0], + point.naturalCoordinates[1], + point.naturalCoordinates[2], + &tyingSamples); + for (std::size_t row = 0U; row < kPhysicalDofCount; ++row) { + for (std::size_t column = row; + column < kPhysicalDofCount; ++column) { + double integrand = 0.0; + for (std::size_t first = 0U; first < 5U; ++first) { + for (std::size_t second = 0U; second < 5U; ++second) { + integrand += strain(first, row) * + constitutive(first, second) * + strain(second, column); + } + } + const double contribution = + integrand * geometry.jacobian * point.weight; + physicalLocal(row, column) += contribution; + if (row != column) { + physicalLocal(column, row) += contribution; + } + } + } + } + + double drillingReference = (std::numeric_limits::max)(); + bool hasDrillingReference = false; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const std::size_t offset = node * kPhysicalDofsPerNode; + for (std::size_t rotation = 3U; rotation < 5U; ++rotation) { + const double diagonal = physicalLocal(offset + rotation, offset + rotation); + if (std::isfinite(diagonal) && diagonal > 0.0) { + drillingReference = (std::min)(drillingReference, diagonal); + hasDrillingReference = true; + } + } + } + if (!hasDrillingReference) { + return stiffnessFailure( + sourceLocation_, identity_, + "MITC4 drilling stabilization requires a finite positive physical " + "tangent-rotation diagonal."); + } + if (!isFinite(physicalLocal)) { + return stiffnessFailure( + sourceLocation_, identity_, + "MITC4 physical stiffness must contain only finite values."); + } + + const double drillingStiffness = 1.0e-3 * drillingReference; + if (!std::isfinite(drillingStiffness) || !(drillingStiffness > 0.0)) { + return stiffnessFailure( + sourceLocation_, identity_, + "MITC4 drilling stiffness must be finite and positive."); + } + + Matrix physicalGlobal = congruence( + physicalLocal, physicalTransformation20()); + Matrix drillingLocal{kNodeCount, kNodeCount}; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + drillingLocal(node, node) = drillingStiffness; + } + Matrix drillingGlobal = congruence( + drillingLocal, drillingTransformation4()); + Matrix stabilizedGlobal{kGlobalDofCount, kGlobalDofCount}; + for (std::size_t row = 0U; row < kGlobalDofCount; ++row) { + for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { + stabilizedGlobal(row, column) = + physicalGlobal(row, column) + drillingGlobal(row, column); + } + } + if (!isFinite(physicalGlobal) || !isFinite(drillingGlobal) || + !isFinite(stabilizedGlobal)) { + return stiffnessFailure( + sourceLocation_, identity_, + "MITC4 transformed stiffness must contain only finite values."); + } + + return Result::success(Mitc4Stiffness{ + std::move(physicalLocal), + std::move(physicalGlobal), + std::move(drillingGlobal), + std::move(stabilizedGlobal), + drillingStiffness}); +} + Mitc4Shell::Mitc4Shell( std::array coordinates, std::array directors, @@ -527,7 +681,9 @@ Mitc4Shell::Mitc4Shell( Vector3 normalCandidate, double thickness, double youngsModulus, - double poissonRatio) + double poissonRatio, + SourceLocation sourceLocation, + std::string identity) : coordinates_{std::move(coordinates)}, directors_{std::move(directors)}, tangentA_{std::move(tangentA)}, @@ -535,7 +691,9 @@ Mitc4Shell::Mitc4Shell( normalCandidate_{std::move(normalCandidate)}, thickness_{thickness}, youngsModulus_{youngsModulus}, - poissonRatio_{poissonRatio} {} + poissonRatio_{poissonRatio}, + sourceLocation_{std::move(sourceLocation)}, + identity_{std::move(identity)} {} bool Mitc4Shell::evaluateGeometry( double xi, diff --git a/tests/unit/elements/mitc4_shell_test.cpp b/tests/unit/elements/mitc4_shell_test.cpp index 2591b4d..e622023 100644 --- a/tests/unit/elements/mitc4_shell_test.cpp +++ b/tests/unit/elements/mitc4_shell_test.cpp @@ -2,12 +2,15 @@ #include +#include #include #include #include #include #include +#include #include +#include namespace { @@ -124,6 +127,194 @@ bool hasPositiveCholeskyPivots(const fesa::Matrix& matrix) { return true; } +double frobeniusNorm(const fesa::Matrix& matrix) { + double squaredNorm = 0.0; + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + for (std::size_t column = 0U; column < matrix.columns(); ++column) { + squaredNorm += matrix(row, column) * matrix(row, column); + } + } + return std::sqrt(squaredNorm); +} + +double scaledSymmetryError( + const fesa::Matrix& matrix, + std::size_t dofsPerNode, + double elementLength) { + fesa::Matrix difference{matrix.rows(), matrix.columns()}; + fesa::Matrix scaled{matrix.rows(), matrix.columns()}; + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + const double rowScale = row % dofsPerNode < 3U ? elementLength : 1.0; + for (std::size_t column = 0U; column < matrix.columns(); ++column) { + const double columnScale = + column % dofsPerNode < 3U ? elementLength : 1.0; + scaled(row, column) = + rowScale * matrix(row, column) * columnScale; + difference(row, column) = rowScale * + (matrix(row, column) - matrix(column, row)) * columnScale; + } + } + return frobeniusNorm(difference) / frobeniusNorm(scaled); +} + +fesa::Matrix scaledStiffness( + const fesa::Matrix& matrix, + std::size_t dofsPerNode, + double elementLength) { + fesa::Matrix scaled{matrix.rows(), matrix.columns()}; + for (std::size_t row = 0U; row < matrix.rows(); ++row) { + const double rowScale = row % dofsPerNode < 3U ? elementLength : 1.0; + for (std::size_t column = 0U; column < matrix.columns(); ++column) { + const double columnScale = + column % dofsPerNode < 3U ? elementLength : 1.0; + scaled(row, column) = + rowScale * matrix(row, column) * columnScale; + } + } + return scaled; +} + +std::vector symmetricEigenvalues(fesa::Matrix matrix) { + if (matrix.rows() != matrix.columns()) { + throw std::invalid_argument{"Symmetric eigensolve requires a square matrix."}; + } + const std::size_t size = matrix.rows(); + double matrixScale = 0.0; + for (std::size_t row = 0U; row < size; ++row) { + for (std::size_t column = 0U; column < size; ++column) { + matrixScale = (std::max)(matrixScale, std::abs(matrix(row, column))); + } + } + if (matrixScale != 0.0) { + const double convergenceTolerance = 1.0e-14 * matrixScale; + const std::size_t iterationLimit = 100U * size * size; + for (std::size_t iteration = 0U; iteration < iterationLimit; ++iteration) { + std::size_t pivotRow = 0U; + std::size_t pivotColumn = 0U; + double largestOffDiagonal = 0.0; + for (std::size_t row = 0U; row < size; ++row) { + for (std::size_t column = row + 1U; column < size; ++column) { + const double candidate = std::abs(matrix(row, column)); + if (candidate > largestOffDiagonal) { + largestOffDiagonal = candidate; + pivotRow = row; + pivotColumn = column; + } + } + } + if (largestOffDiagonal <= convergenceTolerance) { + break; + } + + const double pivot = matrix(pivotRow, pivotColumn); + const double tau = + (matrix(pivotColumn, pivotColumn) - matrix(pivotRow, pivotRow)) / + (2.0 * pivot); + const double tangent = tau >= 0.0 + ? 1.0 / (tau + std::sqrt(1.0 + tau * tau)) + : -1.0 / (-tau + std::sqrt(1.0 + tau * tau)); + const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent); + const double sine = tangent * cosine; + const double rowDiagonal = matrix(pivotRow, pivotRow); + const double columnDiagonal = matrix(pivotColumn, pivotColumn); + matrix(pivotRow, pivotRow) = rowDiagonal - tangent * pivot; + matrix(pivotColumn, pivotColumn) = columnDiagonal + tangent * pivot; + matrix(pivotRow, pivotColumn) = 0.0; + matrix(pivotColumn, pivotRow) = 0.0; + for (std::size_t index = 0U; index < size; ++index) { + if (index == pivotRow || index == pivotColumn) { + continue; + } + const double rowValue = matrix(index, pivotRow); + const double columnValue = matrix(index, pivotColumn); + const double rotatedRow = cosine * rowValue - sine * columnValue; + const double rotatedColumn = sine * rowValue + cosine * columnValue; + matrix(index, pivotRow) = rotatedRow; + matrix(pivotRow, index) = rotatedRow; + matrix(index, pivotColumn) = rotatedColumn; + matrix(pivotColumn, index) = rotatedColumn; + } + } + } + + std::vector eigenvalues(size); + for (std::size_t index = 0U; index < size; ++index) { + eigenvalues[index] = matrix(index, index); + } + return eigenvalues; +} + +std::size_t numericalRank(const fesa::Matrix& scaledMatrix) { + const auto eigenvalues = symmetricEigenvalues(scaledMatrix); + double spectralScale = 0.0; + for (double eigenvalue : eigenvalues) { + spectralScale = (std::max)(spectralScale, std::abs(eigenvalue)); + } + return static_cast(std::count_if( + eigenvalues.begin(), eigenvalues.end(), [spectralScale](double eigenvalue) { + return std::abs(eigenvalue) > 1.0e-9 * spectralScale; + })); +} + +double symmetricOperatorNorm(const fesa::Matrix& matrix) { + const auto eigenvalues = symmetricEigenvalues(matrix); + double result = 0.0; + for (double eigenvalue : eigenvalues) { + result = (std::max)(result, std::abs(eigenvalue)); + } + return result; +} + +double quadraticEnergy(const fesa::Matrix& stiffness, const fesa::Vector& vector) { + return 0.5 * vector.dot(stiffness.multiply(vector)); +} + +fesa::Vector physicalField(const std::array, 4>& values) { + fesa::Vector result{20U}; + for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) { + for (std::size_t component = 0U; component < 5U; ++component) { + result[5U * nodeIndex + component] = values[nodeIndex][component]; + } + } + return result; +} + +void expectStrain( + const fesa::Mitc4Shell& shell, + const fesa::Vector& field, + double xi, + double eta, + double zeta, + const std::array& expected) { + const auto actual = shell.strainDisplacement20(xi, eta, zeta).multiply(field); + for (std::size_t component = 0U; component < expected.size(); ++component) { + EXPECT_NEAR(actual[component], expected[component], 1.0e-12) + << "component " << component; + } +} + +fesa::Vector physicalRigidMode( + const std::array& nodes, + const Vector3& translation, + const Vector3& rotation) { + fesa::Vector mode{24U}; + for (std::size_t nodeIndex = 0U; nodeIndex < nodes.size(); ++nodeIndex) { + const Vector3 rotationalTranslation = cross( + rotation, nodes[nodeIndex].coordinates); + const std::size_t offset = 6U * nodeIndex; + for (std::size_t component = 0U; component < 3U; ++component) { + mode[offset + component] = + translation[component] + rotationalTranslation[component]; + } + // Remove the director-parallel component: it is numerical drilling, + // not part of the five-DOF physical rigid motion. + mode[offset + 3U] = rotation[0]; + mode[offset + 4U] = rotation[1]; + mode[offset + 5U] = 0.0; + } + return mode; +} + std::array planarNodes() { return { node(1, {-1.0, -1.0, 0.0}), @@ -381,3 +572,254 @@ TEST(Mitc4ShellKinematics, UsesOneFixedTwoByTwoByTwoQuadratureOrder) { EXPECT_DOUBLE_EQ(points[point].weight, 1.0); } } + +// MITC4-KERNEL-001 +TEST(Mitc4ShellKernel, FormsFiniteScaledSymmetricPhysicalAndStabilizedStiffness) { + const auto nodes = planarNodes(); + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + + const auto stiffnessCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value(); + EXPECT_EQ(stiffness.physicalLocal20.rows(), 20U); + EXPECT_EQ(stiffness.physicalLocal20.columns(), 20U); + EXPECT_EQ(stiffness.physicalGlobal24.rows(), 24U); + EXPECT_EQ(stiffness.drillingGlobal24.rows(), 24U); + EXPECT_EQ(stiffness.stabilizedGlobal24.rows(), 24U); + for (const fesa::Matrix* matrix : { + &stiffness.physicalLocal20, + &stiffness.physicalGlobal24, + &stiffness.drillingGlobal24, + &stiffness.stabilizedGlobal24}) { + for (std::size_t row = 0U; row < matrix->rows(); ++row) { + for (std::size_t column = 0U; column < matrix->columns(); ++column) { + EXPECT_TRUE(std::isfinite((*matrix)(row, column))); + } + } + } + EXPECT_LE(scaledSymmetryError(stiffness.physicalLocal20, 5U, 2.0), 1.0e-12); + EXPECT_LE(scaledSymmetryError(stiffness.physicalGlobal24, 6U, 2.0), 1.0e-12); + EXPECT_LE(scaledSymmetryError(stiffness.drillingGlobal24, 6U, 2.0), 1.0e-12); + EXPECT_LE(scaledSymmetryError(stiffness.stabilizedGlobal24, 6U, 2.0), 1.0e-12); + + const auto repeatedCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(repeatedCandidate.hasValue()); + const auto& repeated = repeatedCandidate.value(); + expectMatrixNear(repeated.physicalLocal20, stiffness.physicalLocal20, 0.0); + expectMatrixNear(repeated.physicalGlobal24, stiffness.physicalGlobal24, 0.0); + expectMatrixNear(repeated.drillingGlobal24, stiffness.drillingGlobal24, 0.0); + expectMatrixNear(repeated.stabilizedGlobal24, stiffness.stabilizedGlobal24, 0.0); + EXPECT_DOUBLE_EQ(repeated.drillingStiffness, stiffness.drillingStiffness); +} + +// MITC4-KERNEL-002 +TEST(Mitc4ShellKernel, PreservesPhysicalEnergyUnderTwentyToTwentyFourCongruence) { + const auto nodes = planarNodes(); + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + const auto stiffnessCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value(); + + fesa::Vector globalField{24U}; + for (std::size_t index = 0U; index < globalField.size(); ++index) { + globalField[index] = 0.125 * static_cast( + static_cast(index % 7U) - 3); + } + const auto physicalField20 = + shellCandidate.value().physicalTransformation20().multiply(globalField); + const double localEnergy = quadraticEnergy( + stiffness.physicalLocal20, physicalField20); + const double globalEnergy = quadraticEnergy( + stiffness.physicalGlobal24, globalField); + ASSERT_NE(localEnergy, 0.0); + ASSERT_NE(globalEnergy, 0.0); + EXPECT_LE( + std::abs(globalEnergy - localEnergy) / + (std::abs(globalEnergy) + std::abs(localEnergy)), + 1.0e-12); +} + +// MITC4-KERNEL-003 +TEST(Mitc4ShellKernel, RetainsSixRigidModesAndHasExpectedPhysicalAndStabilizedRank) { + const auto nodes = planarNodes(); + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + const auto stiffnessCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value(); + + const auto scaledPhysical20 = + scaledStiffness(stiffness.physicalLocal20, 5U, 2.0); + const auto scaledPhysical24 = + scaledStiffness(stiffness.physicalGlobal24, 6U, 2.0); + const auto scaledStabilized24 = + scaledStiffness(stiffness.stabilizedGlobal24, 6U, 2.0); + EXPECT_EQ(numericalRank(scaledPhysical20), 14U); + EXPECT_EQ(numericalRank(scaledStabilized24), 18U); + + const std::array rigidModes{ + physicalRigidMode(nodes, {1.0, 0.0, 0.0}, {}), + physicalRigidMode(nodes, {0.0, 1.0, 0.0}, {}), + physicalRigidMode(nodes, {0.0, 0.0, 1.0}, {}), + physicalRigidMode(nodes, {}, {1.0, 0.0, 0.0}), + physicalRigidMode(nodes, {}, {0.0, 1.0, 0.0}), + physicalRigidMode(nodes, {}, {0.0, 0.0, 1.0})}; + const double physicalNorm = symmetricOperatorNorm(scaledPhysical24); + const double stabilizedNorm = symmetricOperatorNorm(scaledStabilized24); + ASSERT_GT(physicalNorm, 0.0); + ASSERT_GT(stabilizedNorm, 0.0); + for (const auto& rigidMode : rigidModes) { + fesa::Vector scaledMode = rigidMode; + for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) { + for (std::size_t component = 0U; component < 3U; ++component) { + scaledMode[6U * nodeIndex + component] /= 2.0; + } + } + const double modeNorm = scaledMode.norm(); + ASSERT_GT(modeNorm, 0.0); + EXPECT_LE( + scaledPhysical24.multiply(scaledMode).norm() / + (physicalNorm * modeNorm), + 1.0e-10); + EXPECT_LE( + scaledStabilized24.multiply(scaledMode).norm() / + (stabilizedNorm * modeNorm), + 1.0e-10); + } +} + +// MITC4-KERNEL-004 +TEST(Mitc4ShellPatch, ReproducesIndependentMembraneBendingShearAndTwistFields) { + const auto nodes = planarNodes(); + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + const auto& shell = shellCandidate.value(); + const auto stiffnessCandidate = shell.stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value().physicalLocal20; + constexpr double magnitude = 0.2; + const double gauss = 1.0 / std::sqrt(3.0); + + std::array, 4> e11Values{}; + std::array, 4> e22Values{}; + std::array, 4> g12Values{}; + std::array, 4> k11Values{}; + std::array, 4> k22Values{}; + std::array, 4> g13Values{}; + std::array, 4> g23Values{}; + std::array, 4> k12Values{}; + for (std::size_t nodeIndex = 0U; nodeIndex < nodes.size(); ++nodeIndex) { + const double x = nodes[nodeIndex].coordinates[0]; + const double y = nodes[nodeIndex].coordinates[1]; + e11Values[nodeIndex][0] = magnitude * x; + e22Values[nodeIndex][1] = magnitude * y; + g12Values[nodeIndex][0] = 0.5 * magnitude * y; + g12Values[nodeIndex][1] = 0.5 * magnitude * x; + k11Values[nodeIndex][4] = magnitude * x; + k22Values[nodeIndex][3] = -magnitude * y; + g13Values[nodeIndex][2] = magnitude * x; + g23Values[nodeIndex][2] = magnitude * y; + k12Values[nodeIndex][2] = -0.5 * magnitude * x * y; + k12Values[nodeIndex][3] = -0.5 * magnitude * x; + k12Values[nodeIndex][4] = 0.5 * magnitude * y; + } + + const std::array fields{ + physicalField(e11Values), physicalField(e22Values), + physicalField(g12Values), physicalField(k11Values), + physicalField(k22Values), physicalField(g13Values), + physicalField(g23Values), physicalField(k12Values)}; + expectStrain(shell, fields[0], gauss, -gauss, gauss, {magnitude, 0.0, 0.0, 0.0, 0.0}); + expectStrain(shell, fields[1], gauss, -gauss, gauss, {0.0, magnitude, 0.0, 0.0, 0.0}); + expectStrain(shell, fields[2], gauss, -gauss, gauss, {0.0, 0.0, magnitude, 0.0, 0.0}); + expectStrain(shell, fields[3], gauss, -gauss, gauss, {gauss * magnitude, 0.0, 0.0, 0.0, 0.0}); + expectStrain(shell, fields[4], gauss, -gauss, gauss, {0.0, gauss * magnitude, 0.0, 0.0, 0.0}); + expectStrain(shell, fields[5], gauss, -gauss, gauss, {0.0, 0.0, 0.0, magnitude, 0.0}); + expectStrain(shell, fields[6], gauss, -gauss, gauss, {0.0, 0.0, 0.0, 0.0, magnitude}); + expectStrain(shell, fields[7], gauss, -gauss, gauss, {0.0, 0.0, gauss * magnitude, 0.0, 0.0}); + for (const auto& field : fields) { + EXPECT_GT(quadraticEnergy(stiffness, field), 0.0); + } +} + +// MITC4-KERNEL-005 +TEST(Mitc4ShellDrilling, UsesOnlyEightPositivePhysicalRotationDiagonalsAndFixedFactor) { + const std::array nodes{ + node(1, {-50.0, -50.0, 0.0}), node(2, {50.0, -50.0, 0.0}), + node(3, {50.0, 50.0, 0.0}), node(4, {-50.0, 50.0, 0.0})}; + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(0.1), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + const auto stiffnessCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value(); + + double expectedReference = (std::numeric_limits::max)(); + double allDiagonalMinimum = (std::numeric_limits::max)(); + for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) { + for (std::size_t rotation = 3U; rotation < 5U; ++rotation) { + const std::size_t index = 5U * nodeIndex + rotation; + const double diagonal = stiffness.physicalLocal20(index, index); + ASSERT_TRUE(std::isfinite(diagonal)); + ASSERT_GT(diagonal, 0.0); + expectedReference = (std::min)(expectedReference, diagonal); + } + for (std::size_t component = 0U; component < 5U; ++component) { + const std::size_t index = 5U * nodeIndex + component; + const double diagonal = stiffness.physicalLocal20(index, index); + if (std::isfinite(diagonal) && diagonal > 0.0) { + allDiagonalMinimum = (std::min)(allDiagonalMinimum, diagonal); + } + } + } + EXPECT_DOUBLE_EQ(stiffness.drillingStiffness, 1.0e-3 * expectedReference); + EXPECT_LT(allDiagonalMinimum, expectedReference); + EXPECT_NE(stiffness.drillingStiffness, 1.0e-3 * allDiagonalMinimum); +} + +// MITC4-KERNEL-006 +TEST(Mitc4ShellDrilling, FailsNonfiniteReferenceAndStabilizesEachPureDrillCoordinate) { + const auto nodes = planarNodes(); + const auto shellCandidate = fesa::Mitc4Shell::create( + nodePointers(nodes), directors(), section(), material()); + ASSERT_TRUE(shellCandidate.hasValue()); + const auto stiffnessCandidate = shellCandidate.value().stiffness(); + ASSERT_TRUE(stiffnessCandidate.hasValue()); + const auto& stiffness = stiffnessCandidate.value(); + + for (std::size_t nodeIndex = 0U; nodeIndex < 4U; ++nodeIndex) { + fesa::Vector pureDrill{24U}; + pureDrill[6U * nodeIndex + 5U] = 1.0; + EXPECT_DOUBLE_EQ( + stiffness.physicalGlobal24.multiply(pureDrill).norm(), 0.0); + const auto drillAction = stiffness.drillingGlobal24.multiply(pureDrill); + EXPECT_DOUBLE_EQ(drillAction[6U * nodeIndex + 5U], stiffness.drillingStiffness); + EXPECT_GT(quadraticEnergy(stiffness.drillingGlobal24, pureDrill), 0.0); + } + + const std::array extremeNodes{ + node(1, {-5.0e9, -5.0e9, 0.0}), node(2, {5.0e9, -5.0e9, 0.0}), + node(3, {5.0e9, 5.0e9, 0.0}), node(4, {-5.0e9, 5.0e9, 0.0})}; + const auto extremeShell = fesa::Mitc4Shell::create( + nodePointers(extremeNodes), directors(), section(1.0), material(1.0e300)); + ASSERT_TRUE(extremeShell.hasValue()); + const auto failure = extremeShell.value().stiffness(); + ASSERT_FALSE(failure.hasValue()); + ASSERT_EQ(failure.status().diagnostics().size(), 1U); + EXPECT_EQ(failure.status().diagnostics()[0].code, "invalid-shell-stiffness"); + const auto repeatedFailure = extremeShell.value().stiffness(); + ASSERT_FALSE(repeatedFailure.hasValue()); + ASSERT_EQ(repeatedFailure.status().diagnostics().size(), 1U); + EXPECT_EQ( + repeatedFailure.status().diagnostics()[0].code, + failure.status().diagnostics()[0].code); + EXPECT_EQ( + repeatedFailure.status().diagnostics()[0].message, + failure.status().diagnostics()[0].message); +}