From cfdac70756948481eaa07a6b0385843f0da541b8 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Sun, 9 Aug 2026 18:37:41 +0900 Subject: [PATCH] feat(linear-static-3d-euler-beam): step 16 - euler-beam-element-review-fix --- ...tic-3d-euler-beam-implementation-report.md | 50 +++ src/fesa/elements/euler_beam_3d.cpp | 65 +++- tests/unit/elements/euler_beam_3d_test.cpp | 346 ++++++++++++++++-- 3 files changed, 431 insertions(+), 30 deletions(-) diff --git a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md index 89cbc52..d094810 100644 --- a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md +++ b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md @@ -580,3 +580,53 @@ local stiffness, constant local equivalent load, and distinct recovery data to Step 18 sparse assembly and Step 22 result recovery without owning element identity, equation numbering, parser coupling, or result persistence. + +### Root numerical review — Fix Round 1 + +- scope: all three `Important` findings addressed in + `src/fesa/elements/euler_beam_3d.cpp` and the existing exact T16 tests; the + two `Minor` findings were intentionally not changed in this round. +- defect_RED: after adding finite-positive derived-overflow, + derived-underflow, and length-scaled-underflow fixtures to existing + `T16-BEAM-007`, + `ctest --test-dir .harness/build -C Debug -R EulerBeam3D.RejectsInvalidGeometryAndProperties --output-on-failure` + exited 8. All three invalid fixtures were accepted; the diagnostic path + observed respectively a nonfinite stiffness with maximum `inf`, an all-zero + stiffness, and a stiffness with a required bending coefficient underflowed + to zero. +- defect_GREEN: `cmake --build .harness/build --config Debug --target fesa_tests` + exited 0 and the same focused CTest exited 0 with 1/1 passed. `create()` now + requires derived `EA`, `GJ`, `EIy`, `EIz` and all ten distinct positive + length-scaled axial/torsion/bending block magnitudes to be finite and + strictly positive. This rejects arithmetic overflow/underflow as + `invalid-beam-property` without a conditioning threshold. +- invariant_hardening: production `normalizedMatrixError` now returns infinite + error for any nonfinite operand, difference, scale, or quotient, so the + two-point/closed-form `1e-12` invariant fails closed. +- analytical_hardening: every nonzero reference in existing `T16-BEAM-006` + now checks `abs(actual-reference)/abs(reference) <= 1e-9`; only the explicitly + zero free-end resultants retain an absolute check. This evidence-only change + was expected to remain GREEN and its focused test passed without a production + correction. +- convergence_hardening: existing `T16-BEAM-005` directly assembles aligned + 1/2/4-element `[v,rz]` cantilever systems from public `localStiffness()` and + `localEquivalentLoad()` only. It checks each exact signed constant transverse + element load, solves a test-only dense reduced system, and uses independent + five-point integration of interior cubic-Hermite displacement error against + the quartic uniform-load solution. Both refinements are monotone and each + observed order is 4 within `1e-8`. No `DLOAD`, parser, Domain, or production + assembly support was introduced. This evidence-only change was expected to + remain GREEN and the T16-005/T16-006 focused run passed 2/2. + +| stage | exact command | exit_code | observed_result | +| --- | --- | ---: | --- | +| FIX1-RED | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D.RejectsInvalidGeometryAndProperties --output-on-failure` | 8 | 0/1 passed; unexpected acceptance exposed nonfinite and zero stiffness outcomes | +| FIX1-GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | corrected source and unchanged ten test names compiled under `/W4 /WX` without a warning | +| FIX1-GREEN-defect | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D.RejectsInvalidGeometryAndProperties --output-on-failure` | 0 | 1/1 passed | +| FIX1-GREEN-evidence | `ctest --test-dir .harness/build -C Debug -R 'EulerBeam3D.(ConstantLineLoadMatchesAllSignedComponents|AnalyticalAxialTorsionAndTwoPlaneBendingRecover)' --output-on-failure` | 0 | 2/2 analytical/convergence hardening tests passed without another production change | +| FIX1-VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | exact MSVC x64 Debug tree configured and generated | +| FIX1-VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | full Debug build completed without a new warning | +| FIX1-VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D --output-on-failure` | 0 | exact 10/10 passed | +| FIX1-VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | 39 total tests discovered, including the same ten EulerBeam3D names | +| FIX1-VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | 39/39 passed | +| FIX1-VERIFY-scans | fail-on-match numerical invariant, dependency, scope, exact-test-count, reference, phase-index, and diff/whitespace scans | 0 | derived rigidity/coefficient and fail-closed checks present; tests 10; leaks 0; reference and phase index unchanged; diff clean | diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp index 6d2a1ac..be9d216 100644 --- a/src/fesa/elements/euler_beam_3d.cpp +++ b/src/fesa/elements/euler_beam_3d.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -149,14 +150,30 @@ double normalizedMatrixError(const Matrix& lhs, const Matrix& rhs) { double scale = 1.0; for (std::size_t row = 0; row < lhs.rows(); ++row) { for (std::size_t column = 0; column < lhs.columns(); ++column) { + const double lhsValue = lhs(row, column); + const double rhsValue = rhs(row, column); + if (!std::isfinite(lhsValue) || !std::isfinite(rhsValue)) { + return std::numeric_limits::infinity(); + } + const double difference = std::abs(lhsValue - rhsValue); + if (!std::isfinite(difference)) { + return std::numeric_limits::infinity(); + } maximumDifference = (std::max)( maximumDifference, - std::abs(lhs(row, column) - rhs(row, column))); - scale = (std::max)(scale, std::abs(lhs(row, column))); - scale = (std::max)(scale, std::abs(rhs(row, column))); + difference); + scale = (std::max)(scale, std::abs(lhsValue)); + scale = (std::max)(scale, std::abs(rhsValue)); } } - return maximumDifference / scale; + if (!std::isfinite(maximumDifference) || !std::isfinite(scale) || + !(scale > 0.0)) { + return std::numeric_limits::infinity(); + } + const double normalizedError = maximumDifference / scale; + return std::isfinite(normalizedError) + ? normalizedError + : std::numeric_limits::infinity(); } std::array generalizedStrain( @@ -288,6 +305,46 @@ Result EulerBeam3D::create( "E, G, A, Iy, Iz, and J must be finite and positive."); } + const auto derivedRigidity = constitutiveDiagonal( + material.youngsModulus, + shearModulus, + section.area, + section.i11, + section.i22, + section.torsionalConstant); + const double lengthSquared = length * length; + const double lengthCubed = lengthSquared * length; + // These are every distinct positive magnitude used by the exact axial, + // torsion, and two bending closed-form blocks. Reject arithmetic + // overflow/underflow without introducing a conditioning threshold. + const std::array requiredStiffnessMagnitudes = { + derivedRigidity[0U], + derivedRigidity[1U], + derivedRigidity[2U], + derivedRigidity[3U], + derivedRigidity[0U] / length, + derivedRigidity[1U] / length, + 12.0 * derivedRigidity[2U] / lengthCubed, + 6.0 * derivedRigidity[2U] / lengthSquared, + 4.0 * derivedRigidity[2U] / length, + 2.0 * derivedRigidity[2U] / length, + 12.0 * derivedRigidity[3U] / lengthCubed, + 6.0 * derivedRigidity[3U] / lengthSquared, + 4.0 * derivedRigidity[3U] / length, + 2.0 * derivedRigidity[3U] / length}; + if (std::any_of( + requiredStiffnessMagnitudes.begin(), + requiredStiffnessMagnitudes.end(), + [](double magnitude) { + return !std::isfinite(magnitude) || !(magnitude > 0.0); + })) { + return modelFailure( + "invalid-beam-property", + section.location, + identity, + "Derived beam stiffness coefficients must be finite and positive."); + } + const Vector3 ey = { eyTrial[0] / eyTrialNorm, eyTrial[1] / eyTrialNorm, diff --git a/tests/unit/elements/euler_beam_3d_test.cpp b/tests/unit/elements/euler_beam_3d_test.cpp index 7c9addf..3f9e4d4 100644 --- a/tests/unit/elements/euler_beam_3d_test.cpp +++ b/tests/unit/elements/euler_beam_3d_test.cpp @@ -76,6 +76,17 @@ double maximumAbsoluteEntry(const Matrix& matrix) { return maximum; } +bool matrixIsFinite(const Matrix& matrix) { + for (std::size_t row = 0; row < matrix.rows(); ++row) { + for (std::size_t column = 0; column < matrix.columns(); ++column) { + if (!std::isfinite(matrix(row, column))) { + return false; + } + } + } + return true; +} + double normalizedMatrixError(const Matrix& actual, const Matrix& expected) { if (actual.rows() != expected.rows() || actual.columns() != expected.columns()) { throw std::invalid_argument{"Matrix comparison requires equal shapes."}; @@ -118,6 +129,11 @@ void expectScaledNear(double actual, double expected, double relativeTolerance) EXPECT_LE(std::abs(actual - expected), relativeTolerance * scale); } +void expectRelativeNear(double actual, double expected, double relativeTolerance) { + ASSERT_NE(expected, 0.0); + EXPECT_LE(std::abs(actual - expected) / std::abs(expected), relativeTolerance); +} + Matrix expectedClosedStiffness(double length, const GeneralBeamSection& section, const LinearElasticMaterial& material) { @@ -311,6 +327,153 @@ Vector solveFixedFirstNode(const Matrix& stiffness, return displacement; } +Vector solveDenseSystem(Matrix matrix, Vector rightHandSide) { + if (matrix.rows() != matrix.columns() || + matrix.rows() != rightHandSide.size()) { + throw std::invalid_argument{"Dense test solve requires a square system."}; + } + + for (std::size_t pivot = 0; pivot < matrix.rows(); ++pivot) { + std::size_t pivotRow = pivot; + for (std::size_t row = pivot + 1U; row < matrix.rows(); ++row) { + if (std::abs(matrix(row, pivot)) > + std::abs(matrix(pivotRow, pivot))) { + pivotRow = row; + } + } + if (!(std::abs(matrix(pivotRow, pivot)) > 0.0) || + !std::isfinite(matrix(pivotRow, pivot))) { + throw std::runtime_error{"Uniform-load test fixture is singular."}; + } + for (std::size_t column = pivot; column < matrix.columns(); ++column) { + std::swap(matrix(pivot, column), matrix(pivotRow, column)); + } + std::swap(rightHandSide[pivot], rightHandSide[pivotRow]); + + const double pivotValue = matrix(pivot, pivot); + for (std::size_t column = pivot; column < matrix.columns(); ++column) { + matrix(pivot, column) /= pivotValue; + } + rightHandSide[pivot] /= pivotValue; + for (std::size_t row = 0; row < matrix.rows(); ++row) { + if (row == pivot) { + continue; + } + const double factor = matrix(row, pivot); + for (std::size_t column = pivot; column < matrix.columns(); ++column) { + matrix(row, column) -= factor * matrix(pivot, column); + } + rightHandSide[row] -= factor * rightHandSide[pivot]; + } + } + return rightHandSide; +} + +Vector solveUniformTransverseCantilever(std::size_t elementCount, + double length, + double lineLoad, + const GeneralBeamSection& section, + const LinearElasticMaterial& material) { + const double elementLength = length / static_cast(elementCount); + const std::size_t systemSize = 2U * (elementCount + 1U); + Matrix assembledStiffness{systemSize, systemSize}; + Vector assembledLoad{systemSize}; + const std::array bendingDofs = {1U, 5U, 7U, 11U}; + + // Test-only direct assembly keeps this evidence at the formulation boundary: + // two [v,rz] DOFs per node, with no Domain, parser, or DLOAD path. + for (std::size_t element = 0; element < elementCount; ++element) { + const EulerBeam3D beam = alignedBeam(elementLength, section, material); + const Matrix elementStiffness = beam.localStiffness(); + const Vector elementLoad = beam.localEquivalentLoad( + {0.0, lineLoad, 0.0, 0.0}); + const std::array assembledDofs = { + 2U * element, + 2U * element + 1U, + 2U * (element + 1U), + 2U * (element + 1U) + 1U}; + for (std::size_t row = 0; row < bendingDofs.size(); ++row) { + assembledLoad[assembledDofs[row]] += elementLoad[bendingDofs[row]]; + for (std::size_t column = 0; column < bendingDofs.size(); ++column) { + assembledStiffness(assembledDofs[row], assembledDofs[column]) += + elementStiffness(bendingDofs[row], bendingDofs[column]); + } + } + } + + const std::size_t freeSize = systemSize - 2U; + Matrix freeStiffness{freeSize, freeSize}; + Vector freeLoad{freeSize}; + for (std::size_t row = 0; row < freeSize; ++row) { + freeLoad[row] = assembledLoad[row + 2U]; + for (std::size_t column = 0; column < freeSize; ++column) { + freeStiffness(row, column) = + assembledStiffness(row + 2U, column + 2U); + } + } + + const Vector freeDisplacement = + solveDenseSystem(std::move(freeStiffness), std::move(freeLoad)); + Vector nodalDisplacement{systemSize}; + for (std::size_t dof = 0; dof < freeSize; ++dof) { + nodalDisplacement[dof + 2U] = freeDisplacement[dof]; + } + return nodalDisplacement; +} + +double uniformLoadInteriorDisplacementError( + const Vector& nodalDisplacement, + std::size_t elementCount, + double length, + double lineLoad, + double flexuralRigidity) { + const double elementLength = length / static_cast(elementCount); + const std::array gaussPoints = { + -0.9061798459386640, + -0.5384693101056831, + 0.0, + 0.5384693101056831, + 0.9061798459386640}; + const std::array gaussWeights = { + 0.2369268850561891, + 0.4786286704993665, + 0.5688888888888889, + 0.4786286704993665, + 0.2369268850561891}; + double squaredError = 0.0; + double squaredReference = 0.0; + + // Five-point integration is independent of production and exactly integrates + // the squared error between cubic Hermite interpolation and the quartic beam solution. + for (std::size_t element = 0; element < elementCount; ++element) { + for (std::size_t point = 0; point < gaussPoints.size(); ++point) { + const double r = 0.5 * (1.0 + gaussPoints[point]); + const double rSquared = r * r; + const double rCubed = rSquared * r; + const double h1 = 1.0 - 3.0 * rSquared + 2.0 * rCubed; + const double h2 = elementLength * (r - 2.0 * rSquared + rCubed); + const double h3 = 3.0 * rSquared - 2.0 * rCubed; + const double h4 = elementLength * (-rSquared + rCubed); + const double interpolated = + h1 * nodalDisplacement[2U * element] + + h2 * nodalDisplacement[2U * element + 1U] + + h3 * nodalDisplacement[2U * (element + 1U)] + + h4 * nodalDisplacement[2U * (element + 1U) + 1U]; + const double x = elementLength * + (static_cast(element) + r); + const double analytical = + lineLoad * x * x * + (6.0 * length * length - 4.0 * length * x + x * x) / + (24.0 * flexuralRigidity); + const double weight = 0.5 * elementLength * gaussWeights[point]; + const double difference = interpolated - analytical; + squaredError += weight * difference * difference; + squaredReference += weight * analytical * analytical; + } + } + return std::sqrt(squaredError / squaredReference); +} + Matrix transformationFromKnownRows( const std::array, 3>& rotation) { Matrix transformation{kElementDofCount, kElementDofCount}; @@ -587,6 +750,61 @@ TEST(EulerBeam3D, ConstantLineLoadMatchesAllSignedComponents) { for (std::size_t index = 0; index < expected.size(); ++index) { expectScaledNear(equivalent[index], expected[index], kMatrixTolerance); } + + auto convergenceSection = makeSection(); + convergenceSection.area = 1.0; + convergenceSection.i11 = 1.0; + convergenceSection.i22 = 1.0; + convergenceSection.torsionalConstant = 1.0; + const auto convergenceMaterial = makeMaterial(5.0, 0.25); + const double transverseLoad = -3.0; + const std::array elementCounts = {1U, 2U, 4U}; + std::array relativeErrors{}; + for (std::size_t mesh = 0; mesh < elementCounts.size(); ++mesh) { + const double elementLength = + length / static_cast(elementCounts[mesh]); + const Vector elementLoad = alignedBeam( + elementLength, + convergenceSection, + convergenceMaterial) + .localEquivalentLoad( + {0.0, transverseLoad, 0.0, 0.0}); + expectRelativeNear( + elementLoad[1U], transverseLoad * elementLength / 2.0, kMatrixTolerance); + expectRelativeNear( + elementLoad[5U], + transverseLoad * elementLength * elementLength / 12.0, + kMatrixTolerance); + expectRelativeNear( + elementLoad[7U], transverseLoad * elementLength / 2.0, kMatrixTolerance); + expectRelativeNear( + elementLoad[11U], + -transverseLoad * elementLength * elementLength / 12.0, + kMatrixTolerance); + + const Vector nodalDisplacement = solveUniformTransverseCantilever( + elementCounts[mesh], + length, + transverseLoad, + convergenceSection, + convergenceMaterial); + relativeErrors[mesh] = uniformLoadInteriorDisplacementError( + nodalDisplacement, + elementCounts[mesh], + length, + transverseLoad, + convergenceMaterial.youngsModulus * convergenceSection.i22); + } + EXPECT_GT(relativeErrors[0U], relativeErrors[1U]); + EXPECT_GT(relativeErrors[1U], relativeErrors[2U]); + EXPECT_NEAR( + std::log(relativeErrors[0U] / relativeErrors[1U]) / std::log(2.0), + 4.0, + 1.0e-8); + EXPECT_NEAR( + std::log(relativeErrors[1U] / relativeErrors[2U]) / std::log(2.0), + 4.0, + 1.0e-8); } TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) { @@ -600,75 +818,108 @@ TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) { const double axialForce = 1250.0; const Vector axial = solveFixedFirstNode(stiffness, {axialForce, 0.0, 0.0, 0.0, 0.0, 0.0}); - expectScaledNear( + expectRelativeNear( axial[6U], axialForce * length / (material.youngsModulus * section.area), kAnalyticalTolerance); const BeamRecovery axialRecovery = beam.recover(axial); - expectScaledNear(axialRecovery.equilibriumEndActions[0U][0U], -axialForce, kMatrixTolerance); - expectScaledNear(axialRecovery.equilibriumEndActions[1U][0U], axialForce, kMatrixTolerance); - expectScaledNear(axialRecovery.endpointSectionResultants[0U][0U], axialForce, kMatrixTolerance); - expectScaledNear(axialRecovery.endpointSectionResultants[1U][0U], axialForce, kMatrixTolerance); + expectRelativeNear( + axialRecovery.equilibriumEndActions[0U][0U], + -axialForce, + kAnalyticalTolerance); + expectRelativeNear( + axialRecovery.equilibriumEndActions[1U][0U], + axialForce, + kAnalyticalTolerance); + expectRelativeNear( + axialRecovery.endpointSectionResultants[0U][0U], + axialForce, + kAnalyticalTolerance); + expectRelativeNear( + axialRecovery.endpointSectionResultants[1U][0U], + axialForce, + kAnalyticalTolerance); const double torque = -870.0; const Vector torsion = solveFixedFirstNode(stiffness, {0.0, 0.0, 0.0, torque, 0.0, 0.0}); - expectScaledNear( + expectRelativeNear( torsion[9U], torque * length / (shearModulus * section.torsionalConstant), kAnalyticalTolerance); const BeamRecovery torsionRecovery = beam.recover(torsion); - expectScaledNear(torsionRecovery.equilibriumEndActions[0U][3U], -torque, kMatrixTolerance); - expectScaledNear(torsionRecovery.equilibriumEndActions[1U][3U], torque, kMatrixTolerance); - expectScaledNear(torsionRecovery.endpointSectionResultants[0U][1U], torque, kMatrixTolerance); + expectRelativeNear( + torsionRecovery.equilibriumEndActions[0U][3U], + -torque, + kAnalyticalTolerance); + expectRelativeNear( + torsionRecovery.equilibriumEndActions[1U][3U], + torque, + kAnalyticalTolerance); + expectRelativeNear( + torsionRecovery.endpointSectionResultants[0U][1U], + torque, + kAnalyticalTolerance); const double localYForce = 640.0; const Vector localY = solveFixedFirstNode(stiffness, {0.0, localYForce, 0.0, 0.0, 0.0, 0.0}); - expectScaledNear( + expectRelativeNear( localY[7U], localYForce * length * length * length / (3.0 * material.youngsModulus * section.i22), kAnalyticalTolerance); - expectScaledNear( + expectRelativeNear( localY[11U], localYForce * length * length / (2.0 * material.youngsModulus * section.i22), kAnalyticalTolerance); const BeamRecovery localYRecovery = beam.recover(localY); - expectScaledNear(localYRecovery.equilibriumEndActions[0U][1U], -localYForce, kMatrixTolerance); - expectScaledNear(localYRecovery.equilibriumEndActions[1U][1U], localYForce, kMatrixTolerance); - expectScaledNear( + expectRelativeNear( + localYRecovery.equilibriumEndActions[0U][1U], + -localYForce, + kAnalyticalTolerance); + expectRelativeNear( + localYRecovery.equilibriumEndActions[1U][1U], + localYForce, + kAnalyticalTolerance); + expectRelativeNear( localYRecovery.equilibriumEndActions[0U][5U], -localYForce * length, - kMatrixTolerance); - expectScaledNear( + kAnalyticalTolerance); + expectRelativeNear( localYRecovery.endpointSectionResultants[0U][3U], localYForce * length, - kMatrixTolerance); + kAnalyticalTolerance); EXPECT_NEAR(localYRecovery.endpointSectionResultants[1U][3U], 0.0, 1.0e-8); const double localZForce = -510.0; const Vector localZ = solveFixedFirstNode(stiffness, {0.0, 0.0, localZForce, 0.0, 0.0, 0.0}); - expectScaledNear( + expectRelativeNear( localZ[8U], localZForce * length * length * length / (3.0 * material.youngsModulus * section.i11), kAnalyticalTolerance); - expectScaledNear( + expectRelativeNear( localZ[10U], -localZForce * length * length / (2.0 * material.youngsModulus * section.i11), kAnalyticalTolerance); const BeamRecovery localZRecovery = beam.recover(localZ); - expectScaledNear(localZRecovery.equilibriumEndActions[0U][2U], -localZForce, kMatrixTolerance); - expectScaledNear(localZRecovery.equilibriumEndActions[1U][2U], localZForce, kMatrixTolerance); - expectScaledNear( + expectRelativeNear( + localZRecovery.equilibriumEndActions[0U][2U], + -localZForce, + kAnalyticalTolerance); + expectRelativeNear( + localZRecovery.equilibriumEndActions[1U][2U], + localZForce, + kAnalyticalTolerance); + expectRelativeNear( localZRecovery.equilibriumEndActions[0U][4U], localZForce * length, - kMatrixTolerance); - expectScaledNear( + kAnalyticalTolerance); + expectRelativeNear( localZRecovery.endpointSectionResultants[0U][2U], -localZForce * length, - kMatrixTolerance); + kAnalyticalTolerance); EXPECT_NEAR(localZRecovery.endpointSectionResultants[1U][2U], 0.0, 1.0e-8); } @@ -680,7 +931,14 @@ TEST(EulerBeam3D, RejectsInvalidGeometryAndProperties) { const auto expectFailure = [](const Result& result, const std::string& code) { - ASSERT_FALSE(result.hasValue()); + if (result.hasValue()) { + const Matrix stiffness = result.value().localStiffness(); + ADD_FAILURE() + << "Invalid fixture was accepted; local stiffness finite=" + << matrixIsFinite(stiffness) + << ", maximum absolute entry=" << maximumAbsoluteEntry(stiffness); + return; + } EXPECT_EQ(result.status().failureCategory(), FailureCategory::model); ASSERT_EQ(result.status().diagnostics().size(), 1U); EXPECT_EQ(result.status().diagnostics()[0U].code, code); @@ -746,6 +1004,42 @@ TEST(EulerBeam3D, RejectsInvalidGeometryAndProperties) { "invalid-beam-property"); } + auto overflowMaterial = makeMaterial( + std::numeric_limits::max() / 4.0, + 0.25); + auto overflowSection = validSection; + overflowSection.area = 8.0; + expectFailure( + EulerBeam3D::create(origin, unitX, overflowSection, overflowMaterial), + "invalid-beam-property"); + + auto underflowSection = validSection; + underflowSection.area = std::numeric_limits::denorm_min(); + underflowSection.i11 = std::numeric_limits::denorm_min(); + underflowSection.i22 = std::numeric_limits::denorm_min(); + underflowSection.torsionalConstant = + std::numeric_limits::denorm_min(); + expectFailure( + EulerBeam3D::create( + origin, + unitX, + underflowSection, + makeMaterial(0.5, 0.25)), + "invalid-beam-property"); + + auto lengthScaledSection = validSection; + lengthScaledSection.area = 1.0; + lengthScaledSection.i11 = 1.0; + lengthScaledSection.i22 = 1.0; + lengthScaledSection.torsionalConstant = 1.0; + expectFailure( + EulerBeam3D::create( + origin, + makeNode({1.0e103, 0.0, 0.0}, 2U), + lengthScaledSection, + makeMaterial(1.0, 0.25)), + "invalid-beam-property"); + auto coupledSection = validSection; coupledSection.i12 = 1.0e-9; expectFailure(