feat(linear-static-3d-euler-beam): step 16 - euler-beam-element-review-fix

2026-08-09 18:37:41 +09:00
parent c4ffe13477
commit cfdac70756
3 changed files with 431 additions and 30 deletions
@@ -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<double>(elementCount);
+    const std::size_t systemSize = 2U * (elementCount + 1U);
+    Matrix assembledStiffness{systemSize, systemSize};
+    Vector assembledLoad{systemSize};
+    const std::array<std::size_t, 4> 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<std::size_t, 4> 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<double>(elementCount);
+    const std::array<double, 5> gaussPoints = {
+        -0.9061798459386640,
+        -0.5384693101056831,
+        0.0,
+        0.5384693101056831,
+        0.9061798459386640};
+    const std::array<double, 5> 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<double>(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<std::array<double, 3>, 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<std::size_t, 3> elementCounts = {1U, 2U, 4U};
+    std::array<double, 3> relativeErrors{};
+    for (std::size_t mesh = 0; mesh < elementCounts.size(); ++mesh) {
+        const double elementLength =
+            length / static_cast<double>(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<EulerBeam3D>& 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<double>::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<double>::denorm_min();
+    underflowSection.i11 = std::numeric_limits<double>::denorm_min();
+    underflowSection.i22 = std::numeric_limits<double>::denorm_min();
+    underflowSection.torsionalConstant =
+        std::numeric_limits<double>::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(