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 be0b0ec..89cbc52 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 @@ -513,3 +513,70 @@ - handoff: the exact backend-free result records and `AnalysisState` ledger API provide full-space mutable solution vectors, step/frame identity, and stable recovery-row collections to Step 16 and later result recovery/output tasks. + +## Step 16 — euler-beam-element + +- task_id: `TASK-16` +- status: `completed` +- changed_files: `include/fesa/elements/euler_beam_3d.hpp`, + `src/fesa/elements/euler_beam_3d.cpp`, + `tests/unit/elements/euler_beam_3d_test.cpp`, `src/fesa/CMakeLists.txt`, + `tests/CMakeLists.txt`, + `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`, + `phases/linear-static-3d-euler-beam/index.json`, + `.superpowers/sdd/linear-static-3d-euler-beam/task-16-report.md` +- requirement_ids: `FESA-REQ-LS3DEB-002`, `FESA-REQ-LS3DEB-004`, + `FESA-REQ-LS3DEB-008`, `FESA-REQ-LS3DEB-012`, + `FESA-REQ-LS3DEB-013`, `FESA-REQ-LS3DEB-014`, + `FESA-REQ-LS3DEB-031`, `FESA-REQ-LS3DEB-032`, + `FESA-REQ-LS3DEB-035`, `FESA-REQ-LS3DEB-044` +- test_ids: `T16-BEAM-001`, `T16-BEAM-002`, `T16-BEAM-003`, + `T16-BEAM-004`, `T16-BEAM-005`, `T16-BEAM-006`, `T16-BEAM-007`, + `T16-BEAM-008`, `T16-BEAM-009`, `T16-BEAM-010` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 1 | Exactly ten planned tests were registered before production and the EulerBeam3D public API was absent | MSVC C1083 reported missing `fesa/elements/euler_beam_3d.hpp` from `euler_beam_3d_test.cpp` after successful CMake regeneration | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | The minimum beam kernel, ten tests, solver library, and unit executable compile and link | `euler_beam_3d.cpp`, `euler_beam_3d_test.cpp`, `fesa_solver.lib`, and `fesa_unit_tests.exe` built without a FESA warning under `/W4 /WX` | +| GREEN-test | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D --output-on-failure` | 0 | Reviewed signs, Gauss integration, rank, transform, loads, analytical modes, validation, recovery, patches, and negative control pass | 10/10 exact `EulerBeam3D` tests passed | +| 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 | Approved explicit-dependency MSVC x64 build tree generates | Windows SDK, oneMKL 2026.1, oneTBB, and HDF5 resolved; configure and generate completed | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R EulerBeam3D --output-on-failure` | 0 | Focused Step 16 suite remains green | 10/10 exact `EulerBeam3D` tests passed | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated suite and all ten exact EulerBeam3D names | 39 tests discovered, including 10 `EulerBeam3D` tests, with feature/unit labels | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 39/39 tests passed | +| VERIFY-contract-scans | Backend public-header, upward-dependency, out-of-scope formulation, production one-point path, two-point Gauss/invariant, exact-test-count, and CMake-registration scans using fail-on-match `rg` wrappers | 0 | The kernel implements only the approved exact ledger and numerical integration contract | backend leaks 0; upward dependencies 0; out-of-scope couplings 0; production one-point paths 0; two-point Gauss consumers 3; normalized closed-form invariant `1e-12`; tests 10; registrations 1/1 | +| VERIFY-diff | `git diff --check` plus trailing-whitespace scan over the three new files | 0 | Tracked and untracked Step 16 files have no whitespace errors | Diff check exit 0; new-file trailing whitespace matches 0 | +| VERIFY-reference | `git diff --exit-code -- reference/`; `git status --short -- reference/` | 0 | Approved legacy reference artifacts remain read-only and unchanged | Reference diff exit 0 and reference status empty | + +- green_triage: the first implementation build passed and the first focused + run was 9/10. The single failure was isolated to a large-coordinate test + fixture whose requested equality rounded to a represented length above the + threshold. The test fixture was replaced by an exactly representable unit + threshold equality plus below/above large-coordinate cases; the production + strict `>` comparison and formulation were unchanged before the 10/10 run. +- contract_checks: local DOF order is exactly + `[u1,v1,w1,rx1,ry1,rz1,u2,v2,w2,rx2,ry2,rz2]`. Rotation rows are + `(ex,ey,ez)`, `dl=T*dg`, and `Kg=T^T*Kl*T`. The B-matrix follows the reviewed + `theta_y=-w'`, `theta_z=v'` signs; two-point Gauss `B^T D B` is the sole + production stiffness path and is checked against the independent closed + matrix with normalized tolerance `1e-12`. The test-only one-point rule has + rank four while the production rule has rank six and six rigid modes. +- contract_checks: length and guide-vector tests use the exact scale-aware + strict thresholds with no fallback axis. `E`, derived `G`, `A`, `Iy`, `Iz`, + and `J` must be finite and positive and `I12` must be exactly zero. Constant + local `px,py,pz,mx` receives the consistent signed load only; no parser, + `*DLOAD`, B31, Timoshenko, shear correction, or transverse/torsional stress + support was added. +- contract_checks: recovery keeps local outward equilibrium end actions, + endpoint section resultants, and two-Gauss generalized strain/resultants + distinct. Axial stress is + `S11=E(epsilon0+x2*kappa_y-x1*kappa_z)` in input section-point order; absent + input points produce section point 0 at the centroid with source + `fesa-default`. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`, + `.harness/build/tests/Debug/fesa_unit_tests.exe` +- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0 +- handoff: the exact backend-free `EulerBeam3D` ledger API supplies global and + 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. diff --git a/include/fesa/elements/euler_beam_3d.hpp b/include/fesa/elements/euler_beam_3d.hpp new file mode 100644 index 0000000..4867c55 --- /dev/null +++ b/include/fesa/elements/euler_beam_3d.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include "fesa/core/status.hpp" +#include "fesa/math/matrix.hpp" +#include "fesa/math/vector.hpp" +#include "fesa/model/model_types.hpp" + +#include +#include +#include +#include + +namespace fesa { + +struct ConstantLocalLineLoad { + double px; + double py; + double pz; + double mx; +}; + +struct BeamStressPoint { + int gaussPoint; + std::size_t sectionPoint; + double x1; + double x2; + double s11; + std::string source; +}; + +struct BeamRecovery { + std::array, 2> equilibriumEndActions; + std::array, 2> endpointSectionResultants; + std::array, 2> gaussGeneralizedStrains; + std::array, 2> gaussGeneralizedResultants; + std::vector stressPoints; +}; + +// Implements the approved two-node straight prismatic B33 Euler-Bernoulli +// kernel. Equation numbering and element identity remain outside this type. +class EulerBeam3D { +public: + static Result create(const Node& firstNode, + const Node& secondNode, + const GeneralBeamSection& section, + const LinearElasticMaterial& material); + Matrix localStiffness() const; + Matrix globalStiffness() const; + Vector localEquivalentLoad(const ConstantLocalLineLoad& load) const; + BeamRecovery recover(const Vector& globalElementDisplacement) const; + +private: + EulerBeam3D(double length, + double youngsModulus, + double shearModulus, + double area, + double iy, + double iz, + double torsionalConstant, + std::array rotation, + std::vector> sectionPoints); + + double length_; + double youngsModulus_; + double shearModulus_; + double area_; + double iy_; + double iz_; + double torsionalConstant_; + std::array rotation_; + std::vector> sectionPoints_; +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index be98364..2971e04 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -6,6 +6,7 @@ add_library( build_info.cpp core/diagnostic.cpp core/status.cpp + elements/euler_beam_3d.cpp fem/dof_manager.cpp io/abaqus/domain_mapper.cpp io/abaqus/input_reader.cpp diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp new file mode 100644 index 0000000..6d2a1ac --- /dev/null +++ b/src/fesa/elements/euler_beam_3d.cpp @@ -0,0 +1,458 @@ +#include "fesa/elements/euler_beam_3d.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr std::size_t kElementDofCount = 12U; +constexpr std::size_t kGeneralizedComponentCount = 4U; +constexpr double kGeometryTolerance = 1.0e-12; +constexpr double kStiffnessInvariantTolerance = 1.0e-12; + +using Vector3 = std::array; + +double norm(const Vector3& value) { + return std::hypot(value[0], value[1], value[2]); +} + +double dot(const Vector3& lhs, const Vector3& rhs) { + return lhs[0] * rhs[0] + lhs[1] * rhs[1] + lhs[2] * rhs[2]; +} + +Vector3 cross(const Vector3& lhs, const Vector3& rhs) { + return { + lhs[1] * rhs[2] - lhs[2] * rhs[1], + lhs[2] * rhs[0] - lhs[0] * rhs[2], + lhs[0] * rhs[1] - lhs[1] * rhs[0]}; +} + +bool isFinite(const Vector3& value) { + return std::isfinite(value[0]) && std::isfinite(value[1]) && + std::isfinite(value[2]); +} + +std::string elementIdentity(const Node& firstNode, const Node& secondNode) { + return firstNode.sourceId.instanceName + ":" + + firstNode.sourceId.sourceLabelText + "-" + + secondNode.sourceId.sourceLabelText; +} + +Result modelFailure(const std::string& code, + const SourceLocation& location, + const std::string& identity, + const std::string& message) { + return Result::failure(Status::failure( + FailureCategory::model, + {{Severity::error, code, location, "*ELEMENT", identity, message}})); +} + +Matrix transformation(const std::array& rotation) { + Matrix result{kElementDofCount, kElementDofCount}; + // Blocks preserve [translation, rotation] at node 1 then node 2. + for (std::size_t block = 0; block < 4U; ++block) { + for (std::size_t row = 0; row < 3U; ++row) { + for (std::size_t column = 0; column < 3U; ++column) { + result(block * 3U + row, block * 3U + column) = + rotation[row * 3U + column]; + } + } + } + return result; +} + +Matrix strainDisplacement(double xi, double length) { + Matrix b{kGeneralizedComponentCount, kElementDofCount}; + const double r = 0.5 * (1.0 + xi); + const double inverseLength = 1.0 / length; + const double inverseLengthSquared = inverseLength * inverseLength; + + b(0U, 0U) = -inverseLength; + b(0U, 6U) = inverseLength; + b(1U, 3U) = -inverseLength; + b(1U, 9U) = inverseLength; + + // theta_y=-w' makes kappa_y=-w''; theta_z=v' makes kappa_z=v''. + b(2U, 2U) = (6.0 - 12.0 * r) * inverseLengthSquared; + b(2U, 4U) = (-4.0 + 6.0 * r) * inverseLength; + b(2U, 8U) = (-6.0 + 12.0 * r) * inverseLengthSquared; + b(2U, 10U) = (-2.0 + 6.0 * r) * inverseLength; + + b(3U, 1U) = (-6.0 + 12.0 * r) * inverseLengthSquared; + b(3U, 5U) = (-4.0 + 6.0 * r) * inverseLength; + b(3U, 7U) = (6.0 - 12.0 * r) * inverseLengthSquared; + b(3U, 11U) = (-2.0 + 6.0 * r) * inverseLength; + return b; +} + +std::array constitutiveDiagonal( + double youngsModulus, + double shearModulus, + double area, + double iy, + double iz, + double torsionalConstant) { + return { + youngsModulus * area, + shearModulus * torsionalConstant, + youngsModulus * iy, + youngsModulus * iz}; +} + +Matrix closedStiffness(double length, + const std::array& diagonal) { + Matrix closed{kElementDofCount, kElementDofCount}; + const auto addBlock = [&closed](const std::array& indices, + double coefficient) { + closed(indices[0], indices[0]) = coefficient; + closed(indices[0], indices[1]) = -coefficient; + closed(indices[1], indices[0]) = -coefficient; + closed(indices[1], indices[1]) = coefficient; + }; + addBlock({0U, 6U}, diagonal[0U] / length); + addBlock({3U, 9U}, diagonal[1U] / length); + + const auto addBendingBlock = [&closed, length]( + const std::array& indices, + double flexuralRigidity, + double rotationSign) { + const double value = 12.0 * flexuralRigidity / + (length * length * length); + const double coupling = rotationSign * 6.0 * flexuralRigidity / + (length * length); + const double diagonalRotation = 4.0 * flexuralRigidity / length; + const double offDiagonalRotation = 2.0 * flexuralRigidity / length; + const std::array block = { + value, coupling, -value, coupling, + coupling, diagonalRotation, -coupling, offDiagonalRotation, + -value, -coupling, value, -coupling, + coupling, offDiagonalRotation, -coupling, diagonalRotation}; + for (std::size_t row = 0; row < indices.size(); ++row) { + for (std::size_t column = 0; column < indices.size(); ++column) { + closed(indices[row], indices[column]) = block[row * indices.size() + column]; + } + } + }; + addBendingBlock({1U, 5U, 7U, 11U}, diagonal[3U], 1.0); + addBendingBlock({2U, 4U, 8U, 10U}, diagonal[2U], -1.0); + return closed; +} + +double normalizedMatrixError(const Matrix& lhs, const Matrix& rhs) { + double maximumDifference = 0.0; + double scale = 1.0; + for (std::size_t row = 0; row < lhs.rows(); ++row) { + for (std::size_t column = 0; column < lhs.columns(); ++column) { + 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))); + } + } + return maximumDifference / scale; +} + +std::array generalizedStrain( + const Matrix& b, + const Vector& localDisplacement) { + std::array strain{}; + for (std::size_t component = 0; component < strain.size(); ++component) { + for (std::size_t dof = 0; dof < localDisplacement.size(); ++dof) { + strain[component] += b(component, dof) * localDisplacement[dof]; + } + } + return strain; +} + +std::array generalizedResultant( + const std::array& strain, + const std::array& diagonal) { + std::array resultant{}; + for (std::size_t component = 0; component < resultant.size(); ++component) { + resultant[component] = diagonal[component] * strain[component]; + } + return resultant; +} + +Matrix kinematicInterpolation(double xi, double length) { + Matrix interpolation{4U, kElementDofCount}; + const double r = 0.5 * (1.0 + xi); + const double rSquared = r * r; + const double rCubed = rSquared * r; + const double n1 = 1.0 - r; + const double n2 = r; + const double h1 = 1.0 - 3.0 * rSquared + 2.0 * rCubed; + const double h2 = length * (r - 2.0 * rSquared + rCubed); + const double h3 = 3.0 * rSquared - 2.0 * rCubed; + const double h4 = length * (-rSquared + rCubed); + + interpolation(0U, 0U) = n1; + interpolation(0U, 6U) = n2; + interpolation(1U, 1U) = h1; + interpolation(1U, 5U) = h2; + interpolation(1U, 7U) = h3; + interpolation(1U, 11U) = h4; + interpolation(2U, 2U) = h1; + interpolation(2U, 4U) = -h2; + interpolation(2U, 8U) = h3; + interpolation(2U, 10U) = -h4; + interpolation(3U, 3U) = n1; + interpolation(3U, 9U) = n2; + return interpolation; +} + +} // namespace + +Result EulerBeam3D::create( + const Node& firstNode, + const Node& secondNode, + const GeneralBeamSection& section, + const LinearElasticMaterial& material) { + const std::string identity = elementIdentity(firstNode, secondNode); + const Vector3& first = firstNode.coordinates; + const Vector3& second = secondNode.coordinates; + const Vector3 delta = { + second[0] - first[0], second[1] - first[1], second[2] - first[2]}; + const double length = norm(delta); + const double coordinateScale = + (std::max)({1.0, norm(first), norm(second)}); + if (!isFinite(first) || !isFinite(second) || !isFinite(delta) || + !std::isfinite(length) || !std::isfinite(coordinateScale) || + !(length > kGeometryTolerance * coordinateScale)) { + return modelFailure( + "invalid-beam-length", + firstNode.location, + identity, + "Beam length must exceed the scale-aware geometry threshold."); + } + + const Vector3 ex = {delta[0] / length, delta[1] / length, delta[2] / length}; + const Vector3& guide = section.firstAxis; + const double guideNorm = norm(guide); + const double guideProjection = dot(guide, ex); + const Vector3 eyTrial = { + guide[0] - guideProjection * ex[0], + guide[1] - guideProjection * ex[1], + guide[2] - guideProjection * ex[2]}; + const double eyTrialNorm = norm(eyTrial); + if (!isFinite(guide) || !std::isfinite(guideNorm) || !isFinite(eyTrial) || + !std::isfinite(eyTrialNorm) || + !(eyTrialNorm > kGeometryTolerance * (std::max)(1.0, guideNorm))) { + return modelFailure( + "invalid-beam-guide-vector", + section.location, + identity, + "Beam guide vector must define a scale-aware transverse direction."); + } + + if (!std::isfinite(section.i12)) { + return modelFailure( + "invalid-beam-property", + section.location, + identity, + "Beam section properties must be finite and positive."); + } + if (section.i12 != 0.0) { + return modelFailure( + "unsupported-coupled-section", + section.location, + identity, + "The Euler beam kernel requires exact I12=0."); + } + + const double shearModulus = + material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); + const std::array positiveProperties = { + material.youngsModulus, + shearModulus, + section.area, + section.i11, + section.i22, + section.torsionalConstant}; + if (!std::isfinite(material.poissonRatio) || + std::any_of( + positiveProperties.begin(), + positiveProperties.end(), + [](double property) { return !std::isfinite(property) || !(property > 0.0); })) { + return modelFailure( + "invalid-beam-property", + section.location, + identity, + "E, G, A, Iy, Iz, and J must be finite and positive."); + } + + const Vector3 ey = { + eyTrial[0] / eyTrialNorm, + eyTrial[1] / eyTrialNorm, + eyTrial[2] / eyTrialNorm}; + const Vector3 ez = cross(ex, ey); + // Rows map global vectors to the approved right-handed local (ex,ey,ez) basis. + const std::array rotation = { + ex[0], ex[1], ex[2], + ey[0], ey[1], ey[2], + ez[0], ez[1], ez[2]}; + + return Result::success(EulerBeam3D{ + length, + material.youngsModulus, + shearModulus, + section.area, + section.i11, + section.i22, + section.torsionalConstant, + rotation, + section.sectionPoints}); +} + +Matrix EulerBeam3D::localStiffness() const { + const auto diagonal = constitutiveDiagonal( + youngsModulus_, shearModulus_, area_, iy_, iz_, torsionalConstant_); + Matrix stiffness{kElementDofCount, kElementDofCount}; + const double inverseSqrtThree = 1.0 / std::sqrt(3.0); + const std::array gaussPoints = {-inverseSqrtThree, inverseSqrtThree}; + const double jacobian = 0.5 * length_; + + // Both Gauss points are required: a one-point bending rule loses two ranks. + for (const double xi : gaussPoints) { + const Matrix b = strainDisplacement(xi, length_); + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = 0; column < kElementDofCount; ++column) { + for (std::size_t component = 0; component < diagonal.size(); ++component) { + stiffness(row, column) += + b(component, row) * diagonal[component] * + b(component, column) * jacobian; + } + } + } + } + + const Matrix closed = closedStiffness(length_, diagonal); + if (normalizedMatrixError(stiffness, closed) > kStiffnessInvariantTolerance) { + throw std::logic_error{ + "Two-point Euler beam stiffness violated the closed-form invariant."}; + } + return stiffness; +} + +Matrix EulerBeam3D::globalStiffness() const { + const Matrix local = localStiffness(); + const Matrix transform = transformation(rotation_); + const Matrix localTimesTransform = local.multiply(transform); + Matrix global{kElementDofCount, kElementDofCount}; + // Kg=T^T*Kl*T while dl=T*dg. + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = 0; column < kElementDofCount; ++column) { + for (std::size_t inner = 0; inner < kElementDofCount; ++inner) { + global(row, column) += + transform(inner, row) * localTimesTransform(inner, column); + } + } + } + return global; +} + +Vector EulerBeam3D::localEquivalentLoad(const ConstantLocalLineLoad& load) const { + const std::array components = {load.px, load.py, load.pz, load.mx}; + Vector equivalent{kElementDofCount}; + const double inverseSqrtThree = 1.0 / std::sqrt(3.0); + const std::array gaussPoints = {-inverseSqrtThree, inverseSqrtThree}; + const double jacobian = 0.5 * length_; + for (const double xi : gaussPoints) { + const Matrix interpolation = kinematicInterpolation(xi, length_); + for (std::size_t dof = 0; dof < equivalent.size(); ++dof) { + for (std::size_t component = 0; component < components.size(); ++component) { + equivalent[dof] += + interpolation(component, dof) * components[component] * jacobian; + } + } + } + return equivalent; +} + +BeamRecovery EulerBeam3D::recover(const Vector& globalElementDisplacement) const { + const Matrix transform = transformation(rotation_); + const Vector localDisplacement = transform.multiply(globalElementDisplacement); + const auto diagonal = constitutiveDiagonal( + youngsModulus_, shearModulus_, area_, iy_, iz_, torsionalConstant_); + BeamRecovery recovery{}; + + // With parser/CLI distributed loading excluded, Kl*dl is the local outward end action. + const Vector endAction = localStiffness().multiply(localDisplacement); + for (std::size_t endpoint = 0; endpoint < 2U; ++endpoint) { + for (std::size_t component = 0; component < 6U; ++component) { + recovery.equilibriumEndActions[endpoint][component] = + endAction[endpoint * 6U + component]; + } + + const double xi = endpoint == 0U ? -1.0 : 1.0; + recovery.endpointSectionResultants[endpoint] = generalizedResultant( + generalizedStrain(strainDisplacement(xi, length_), localDisplacement), + diagonal); + } + + const double inverseSqrtThree = 1.0 / std::sqrt(3.0); + const std::array gaussPoints = {-inverseSqrtThree, inverseSqrtThree}; + for (std::size_t point = 0; point < gaussPoints.size(); ++point) { + recovery.gaussGeneralizedStrains[point] = generalizedStrain( + strainDisplacement(gaussPoints[point], length_), localDisplacement); + recovery.gaussGeneralizedResultants[point] = generalizedResultant( + recovery.gaussGeneralizedStrains[point], diagonal); + + if (sectionPoints_.empty()) { + recovery.stressPoints.push_back({ + static_cast(point + 1U), + 0U, + 0.0, + 0.0, + youngsModulus_ * recovery.gaussGeneralizedStrains[point][0U], + "fesa-default"}); + continue; + } + + for (std::size_t sectionPoint = 0; sectionPoint < sectionPoints_.size(); + ++sectionPoint) { + const double x1 = sectionPoints_[sectionPoint][0U]; + const double x2 = sectionPoints_[sectionPoint][1U]; + const auto& strain = recovery.gaussGeneralizedStrains[point]; + // x1=y and x2=z: S11=E(epsilon0+x2*kappa_y-x1*kappa_z). + recovery.stressPoints.push_back({ + static_cast(point + 1U), + sectionPoint + 1U, + x1, + x2, + youngsModulus_ * + (strain[0U] + x2 * strain[2U] - x1 * strain[3U]), + "input"}); + } + } + return recovery; +} + +EulerBeam3D::EulerBeam3D( + double length, + double youngsModulus, + double shearModulus, + double area, + double iy, + double iz, + double torsionalConstant, + std::array rotation, + std::vector> sectionPoints) + : length_{length}, + youngsModulus_{youngsModulus}, + shearModulus_{shearModulus}, + area_{area}, + iy_{iy}, + iz_{iz}, + torsionalConstant_{torsionalConstant}, + rotation_{rotation}, + sectionPoints_{std::move(sectionPoints)} {} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c5a6c21..87b9737 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,6 +8,7 @@ add_executable( unit/core/diagnostic_test.cpp unit/core/source_identity_test.cpp unit/core/status_test.cpp + unit/elements/euler_beam_3d_test.cpp unit/fem/dof_manager_test.cpp unit/math/matrix_test.cpp unit/math/vector_test.cpp diff --git a/tests/unit/elements/euler_beam_3d_test.cpp b/tests/unit/elements/euler_beam_3d_test.cpp new file mode 100644 index 0000000..7c9addf --- /dev/null +++ b/tests/unit/elements/euler_beam_3d_test.cpp @@ -0,0 +1,885 @@ +#include "fesa/elements/euler_beam_3d.hpp" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fesa { +namespace { + +constexpr std::size_t kElementDofCount = 12U; +constexpr double kMatrixTolerance = 1.0e-12; +constexpr double kRigidTolerance = 1.0e-10; +constexpr double kAnalyticalTolerance = 1.0e-9; + +Node makeNode(std::array coordinates, std::size_t line) { + return {{"Beam-1", static_cast(line), std::to_string(line)}, + coordinates, + {"beam-test.inp", line}}; +} + +LinearElasticMaterial makeMaterial(double youngsModulus = 210.0e9, + double poissonRatio = 0.3) { + return {"Steel", youngsModulus, poissonRatio, {"beam-test.inp", 20U}}; +} + +GeneralBeamSection makeSection( + std::array firstAxis = {0.0, 1.0, 0.0}, + std::vector> sectionPoints = {}) { + return {"Section-1", + 0.012, + 2.5e-5, + 0.0, + 4.0e-5, + 1.5e-5, + firstAxis, + std::move(sectionPoints), + {"beam-test.inp", 30U}}; +} + +EulerBeam3D requireBeam(const Node& firstNode, + const Node& secondNode, + const GeneralBeamSection& section, + const LinearElasticMaterial& material) { + auto result = EulerBeam3D::create(firstNode, secondNode, section, material); + if (!result.hasValue()) { + throw std::runtime_error{"Expected a valid EulerBeam3D fixture."}; + } + return std::move(result.value()); +} + +EulerBeam3D alignedBeam(double length, + const GeneralBeamSection& section = makeSection(), + const LinearElasticMaterial& material = makeMaterial()) { + return requireBeam( + makeNode({0.0, 0.0, 0.0}, 1U), + makeNode({length, 0.0, 0.0}, 2U), + section, + material); +} + +double maximumAbsoluteEntry(const Matrix& matrix) { + double maximum = 0.0; + for (std::size_t row = 0; row < matrix.rows(); ++row) { + for (std::size_t column = 0; column < matrix.columns(); ++column) { + maximum = (std::max)(maximum, std::abs(matrix(row, column))); + } + } + return maximum; +} + +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."}; + } + + double maximumDifference = 0.0; + for (std::size_t row = 0; row < actual.rows(); ++row) { + for (std::size_t column = 0; column < actual.columns(); ++column) { + maximumDifference = (std::max)( + maximumDifference, + std::abs(actual(row, column) - expected(row, column))); + } + } + + const double scale = (std::max)( + 1.0, + (std::max)(maximumAbsoluteEntry(actual), maximumAbsoluteEntry(expected))); + return maximumDifference / scale; +} + +double vectorNorm(const Vector& vector) { + double sum = 0.0; + for (std::size_t index = 0; index < vector.size(); ++index) { + sum += vector[index] * vector[index]; + } + return std::sqrt(sum); +} + +double quadraticEnergy(const Matrix& matrix, const Vector& vector) { + const Vector product = matrix.multiply(vector); + double value = 0.0; + for (std::size_t index = 0; index < vector.size(); ++index) { + value += vector[index] * product[index]; + } + return value; +} + +void expectScaledNear(double actual, double expected, double relativeTolerance) { + const double scale = (std::max)(1.0, std::abs(expected)); + EXPECT_LE(std::abs(actual - expected), relativeTolerance * scale); +} + +Matrix expectedClosedStiffness(double length, + const GeneralBeamSection& section, + const LinearElasticMaterial& material) { + Matrix expected{kElementDofCount, kElementDofCount}; + const double shearModulus = + material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); + + const auto addBlock = [&expected](const std::vector& indices, + const std::vector& values) { + const std::size_t width = indices.size(); + for (std::size_t row = 0; row < width; ++row) { + for (std::size_t column = 0; column < width; ++column) { + expected(indices[row], indices[column]) = values[row * width + column]; + } + } + }; + + const double axial = material.youngsModulus * section.area / length; + addBlock({0U, 6U}, {axial, -axial, -axial, axial}); + + const double torsion = shearModulus * section.torsionalConstant / length; + addBlock({3U, 9U}, {torsion, -torsion, -torsion, torsion}); + + const auto bendingBlock = [length](double flexuralRigidity, double rotationSign) { + const double v = 12.0 * flexuralRigidity / (length * length * length); + const double c = rotationSign * 6.0 * flexuralRigidity / (length * length); + const double d = 4.0 * flexuralRigidity / length; + const double e = 2.0 * flexuralRigidity / length; + return std::vector{ + v, c, -v, c, + c, d, -c, e, + -v, -c, v, -c, + c, e, -c, d}; + }; + + addBlock( + {1U, 5U, 7U, 11U}, + bendingBlock(material.youngsModulus * section.i22, 1.0)); + addBlock( + {2U, 4U, 8U, 10U}, + bendingBlock(material.youngsModulus * section.i11, -1.0)); + return expected; +} + +std::array symmetricEigenvalues(Matrix matrix) { + for (std::size_t iteration = 0; iteration < 100U * kElementDofCount; ++iteration) { + std::size_t p = 0U; + std::size_t q = 1U; + double maximumOffDiagonal = 0.0; + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = row + 1U; column < kElementDofCount; ++column) { + const double candidate = std::abs(matrix(row, column)); + if (candidate > maximumOffDiagonal) { + maximumOffDiagonal = candidate; + p = row; + q = column; + } + } + } + if (maximumOffDiagonal <= + 1.0e-14 * (std::max)(1.0, maximumAbsoluteEntry(matrix))) { + break; + } + + const double app = matrix(p, p); + const double aqq = matrix(q, q); + const double apq = matrix(p, q); + const double angle = 0.5 * std::atan2(2.0 * apq, aqq - app); + const double cosine = std::cos(angle); + const double sine = std::sin(angle); + + for (std::size_t index = 0; index < kElementDofCount; ++index) { + if (index == p || index == q) { + continue; + } + const double aip = matrix(index, p); + const double aiq = matrix(index, q); + matrix(index, p) = cosine * aip - sine * aiq; + matrix(p, index) = matrix(index, p); + matrix(index, q) = sine * aip + cosine * aiq; + matrix(q, index) = matrix(index, q); + } + + matrix(p, p) = cosine * cosine * app - 2.0 * sine * cosine * apq + + sine * sine * aqq; + matrix(q, q) = sine * sine * app + 2.0 * sine * cosine * apq + + cosine * cosine * aqq; + matrix(p, q) = 0.0; + matrix(q, p) = 0.0; + } + + std::array eigenvalues{}; + for (std::size_t index = 0; index < kElementDofCount; ++index) { + eigenvalues[index] = matrix(index, index); + } + return eigenvalues; +} + +std::size_t symmetricRank(const Matrix& matrix, double relativeTolerance) { + const auto eigenvalues = symmetricEigenvalues(matrix); + double maximum = 0.0; + for (const double value : eigenvalues) { + maximum = (std::max)(maximum, std::abs(value)); + } + return static_cast(std::count_if( + eigenvalues.begin(), + eigenvalues.end(), + [maximum, relativeTolerance](double value) { + return std::abs(value) > relativeTolerance * maximum; + })); +} + +Matrix testOnlyOnePointStiffness(double length, + const GeneralBeamSection& section, + const LinearElasticMaterial& material) { + // At xi=0 the two bending curvature rows retain only the nodal rotations. + // This deliberately under-integrated negative control is independent of production. + Matrix b{4U, kElementDofCount}; + b(0U, 0U) = -1.0 / length; + b(0U, 6U) = 1.0 / length; + b(1U, 3U) = -1.0 / length; + b(1U, 9U) = 1.0 / length; + b(2U, 4U) = -1.0 / length; + b(2U, 10U) = 1.0 / length; + b(3U, 5U) = -1.0 / length; + b(3U, 11U) = 1.0 / length; + + const double shearModulus = + material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); + const std::array diagonal = { + material.youngsModulus * section.area, + shearModulus * section.torsionalConstant, + material.youngsModulus * section.i11, + material.youngsModulus * section.i22}; + + Matrix stiffness{kElementDofCount, kElementDofCount}; + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = 0; column < kElementDofCount; ++column) { + for (std::size_t component = 0; component < diagonal.size(); ++component) { + stiffness(row, column) += + b(component, row) * diagonal[component] * b(component, column) * length; + } + } + } + return stiffness; +} + +Vector solveFixedFirstNode(const Matrix& stiffness, + const std::array& freeEndLoad) { + std::array, 6> augmented{}; + for (std::size_t row = 0; row < 6U; ++row) { + for (std::size_t column = 0; column < 6U; ++column) { + augmented[row][column] = stiffness(row + 6U, column + 6U); + } + augmented[row][6U] = freeEndLoad[row]; + } + + for (std::size_t pivot = 0; pivot < 6U; ++pivot) { + std::size_t pivotRow = pivot; + for (std::size_t row = pivot + 1U; row < 6U; ++row) { + if (std::abs(augmented[row][pivot]) > + std::abs(augmented[pivotRow][pivot])) { + pivotRow = row; + } + } + if (std::abs(augmented[pivotRow][pivot]) <= + std::numeric_limits::min()) { + throw std::runtime_error{"Cantilever fixture is singular."}; + } + std::swap(augmented[pivot], augmented[pivotRow]); + + const double pivotValue = augmented[pivot][pivot]; + for (std::size_t column = pivot; column < 7U; ++column) { + augmented[pivot][column] /= pivotValue; + } + for (std::size_t row = 0; row < 6U; ++row) { + if (row == pivot) { + continue; + } + const double factor = augmented[row][pivot]; + for (std::size_t column = pivot; column < 7U; ++column) { + augmented[row][column] -= factor * augmented[pivot][column]; + } + } + } + + Vector displacement{kElementDofCount}; + for (std::size_t component = 0; component < 6U; ++component) { + displacement[component + 6U] = augmented[component][6U]; + } + return displacement; +} + +Matrix transformationFromKnownRows( + const std::array, 3>& rotation) { + Matrix transformation{kElementDofCount, kElementDofCount}; + for (std::size_t block = 0; block < 4U; ++block) { + for (std::size_t row = 0; row < 3U; ++row) { + for (std::size_t column = 0; column < 3U; ++column) { + transformation(block * 3U + row, block * 3U + column) = + rotation[row][column]; + } + } + } + return transformation; +} + +Vector transposeMultiply(const Matrix& matrix, const Vector& vector) { + if (matrix.rows() != vector.size()) { + throw std::invalid_argument{"Transpose multiply dimension mismatch."}; + } + Vector result{matrix.columns()}; + for (std::size_t column = 0; column < matrix.columns(); ++column) { + for (std::size_t row = 0; row < matrix.rows(); ++row) { + result[column] += matrix(row, column) * vector[row]; + } + } + return result; +} + +double determinant(const std::array, 3>& matrix) { + return matrix[0][0] * + (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) - + matrix[0][1] * + (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) + + matrix[0][2] * + (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]); +} + +TEST(EulerBeam3D, HermiteAndBMatrixMatchReviewedSigns) { + const double length = 2.5; + const auto section = makeSection(); + const auto material = makeMaterial(); + const auto beam = alignedBeam(length, section, material); + + const double axialStrain = 0.012; + const double twist = -0.021; + const std::array v = {1.2, -0.4, 0.3, -0.07}; + const std::array w = {-0.8, 0.6, -0.2, 0.05}; + const auto value = [](const std::array& coefficients, double x) { + return coefficients[0] + coefficients[1] * x + coefficients[2] * x * x + + coefficients[3] * x * x * x; + }; + const auto slope = [](const std::array& coefficients, double x) { + return coefficients[1] + 2.0 * coefficients[2] * x + + 3.0 * coefficients[3] * x * x; + }; + const auto curvature = [](const std::array& coefficients, double x) { + return 2.0 * coefficients[2] + 6.0 * coefficients[3] * x; + }; + + Vector displacement{kElementDofCount}; + displacement[0U] = 0.2; + displacement[1U] = value(v, 0.0); + displacement[2U] = value(w, 0.0); + displacement[3U] = -0.1; + displacement[4U] = -slope(w, 0.0); + displacement[5U] = slope(v, 0.0); + displacement[6U] = displacement[0U] + axialStrain * length; + displacement[7U] = value(v, length); + displacement[8U] = value(w, length); + displacement[9U] = displacement[3U] + twist * length; + displacement[10U] = -slope(w, length); + displacement[11U] = slope(v, length); + + const BeamRecovery recovery = beam.recover(displacement); + const double inverseSqrtThree = 1.0 / std::sqrt(3.0); + const std::array gaussXi = {-inverseSqrtThree, inverseSqrtThree}; + for (std::size_t point = 0; point < gaussXi.size(); ++point) { + const double x = 0.5 * length * (1.0 + gaussXi[point]); + EXPECT_NEAR(recovery.gaussGeneralizedStrains[point][0U], axialStrain, 1.0e-14); + EXPECT_NEAR(recovery.gaussGeneralizedStrains[point][1U], twist, 1.0e-14); + EXPECT_NEAR( + recovery.gaussGeneralizedStrains[point][2U], + -curvature(w, x), + 1.0e-13); + EXPECT_NEAR( + recovery.gaussGeneralizedStrains[point][3U], + curvature(v, x), + 1.0e-13); + } + + EXPECT_NEAR( + recovery.endpointSectionResultants[0U][2U], + material.youngsModulus * section.i11 * -curvature(w, 0.0), + 1.0e-5); + EXPECT_NEAR( + recovery.endpointSectionResultants[1U][2U], + material.youngsModulus * section.i11 * -curvature(w, length), + 1.0e-5); + EXPECT_NEAR( + recovery.endpointSectionResultants[0U][3U], + material.youngsModulus * section.i22 * curvature(v, 0.0), + 1.0e-5); + EXPECT_NEAR( + recovery.endpointSectionResultants[1U][3U], + material.youngsModulus * section.i22 * curvature(v, length), + 1.0e-5); +} + +TEST(EulerBeam3D, TwoPointGaussMatchesClosedStiffness) { + const double length = 3.7; + auto section = makeSection(); + section.area = 0.019; + section.i11 = 3.1e-5; + section.i22 = 7.4e-5; + section.torsionalConstant = 2.2e-5; + const auto material = makeMaterial(73.0e9, 0.27); + const auto beam = alignedBeam(length, section, material); + + const Matrix actual = beam.localStiffness(); + const Matrix closed = expectedClosedStiffness(length, section, material); + EXPECT_LE(normalizedMatrixError(actual, closed), kMatrixTolerance); + + Matrix transpose{actual.rows(), actual.columns()}; + for (std::size_t row = 0; row < actual.rows(); ++row) { + for (std::size_t column = 0; column < actual.columns(); ++column) { + transpose(row, column) = actual(column, row); + } + } + EXPECT_LE(normalizedMatrixError(actual, transpose), kMatrixTolerance); +} + +TEST(EulerBeam3D, HasSixRigidModesRankSixAndPositiveDeformationEnergy) { + const double length = 2.0; + auto section = makeSection(); + section.area = 1.4; + section.i11 = 0.8; + section.i22 = 1.1; + section.torsionalConstant = 0.6; + const auto material = makeMaterial(5.0, 0.25); + const Matrix stiffness = alignedBeam(length, section, material).localStiffness(); + + std::array rigidModes = { + Vector{kElementDofCount}, Vector{kElementDofCount}, Vector{kElementDofCount}, + Vector{kElementDofCount}, Vector{kElementDofCount}, Vector{kElementDofCount}}; + rigidModes[0U][0U] = rigidModes[0U][6U] = 1.0; + rigidModes[1U][1U] = rigidModes[1U][7U] = 1.0; + rigidModes[2U][2U] = rigidModes[2U][8U] = 1.0; + rigidModes[3U][3U] = rigidModes[3U][9U] = 1.0; + rigidModes[4U][4U] = rigidModes[4U][10U] = 1.0; + rigidModes[4U][8U] = -length; + rigidModes[5U][5U] = rigidModes[5U][11U] = 1.0; + rigidModes[5U][7U] = length; + + const double stiffnessScale = (std::max)(1.0, maximumAbsoluteEntry(stiffness)); + for (const Vector& mode : rigidModes) { + const double normalizedResidual = + vectorNorm(stiffness.multiply(mode)) / + (stiffnessScale * (std::max)(1.0, vectorNorm(mode))); + EXPECT_LE(normalizedResidual, kRigidTolerance); + } + + const std::array q = { + 1.0, 1.0, 1.0, length, length, length, + 1.0, 1.0, 1.0, length, length, length}; + Matrix scaled{kElementDofCount, kElementDofCount}; + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = 0; column < kElementDofCount; ++column) { + scaled(row, column) = stiffness(row, column) / (q[row] * q[column]); + } + } + const auto eigenvalues = symmetricEigenvalues(scaled); + double maximumSingularValue = 0.0; + for (const double value : eigenvalues) { + maximumSingularValue = (std::max)(maximumSingularValue, std::abs(value)); + } + const auto positiveCount = std::count_if( + eigenvalues.begin(), eigenvalues.end(), [maximumSingularValue](double value) { + return std::abs(value) > kRigidTolerance * maximumSingularValue; + }); + EXPECT_EQ(positiveCount, 6); + for (const double value : eigenvalues) { + EXPECT_GE(value, -kRigidTolerance * maximumSingularValue); + } + + for (std::size_t component = 0; component < 6U; ++component) { + Vector deformation{kElementDofCount}; + deformation[6U + component] = 1.0; + EXPECT_GT(quadraticEnergy(stiffness, deformation), 0.0); + } +} + +TEST(EulerBeam3D, RotatedTransformPreservesWorkAndEnergy) { + const double inverseSqrtTwo = 1.0 / std::sqrt(2.0); + const std::array, 3> rotation = {{ + {{2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0}}, + {{-inverseSqrtTwo, inverseSqrtTwo, 0.0}}, + {{-inverseSqrtTwo / 3.0, -inverseSqrtTwo / 3.0, + 4.0 * inverseSqrtTwo / 3.0}}}}; + const Matrix transformation = transformationFromKnownRows(rotation); + + for (std::size_t row = 0; row < 3U; ++row) { + for (std::size_t column = 0; column < 3U; ++column) { + double dot = 0.0; + for (std::size_t component = 0; component < 3U; ++component) { + dot += rotation[row][component] * rotation[column][component]; + } + EXPECT_NEAR(dot, row == column ? 1.0 : 0.0, 1.0e-14); + } + } + EXPECT_NEAR(determinant(rotation), 1.0, 1.0e-14); + + const auto section = makeSection({-2.0, 2.0, 0.0}); + const auto material = makeMaterial(); + const auto beam = requireBeam( + makeNode({1.0, -2.0, 0.5}, 1U), + makeNode({3.0, 0.0, 1.5}, 2U), + section, + material); + const Matrix local = beam.localStiffness(); + const Matrix global = beam.globalStiffness(); + + Matrix expectedGlobal{kElementDofCount, kElementDofCount}; + const Matrix localTimesTransform = local.multiply(transformation); + for (std::size_t row = 0; row < kElementDofCount; ++row) { + for (std::size_t column = 0; column < kElementDofCount; ++column) { + for (std::size_t inner = 0; inner < kElementDofCount; ++inner) { + expectedGlobal(row, column) += + transformation(inner, row) * localTimesTransform(inner, column); + } + } + } + EXPECT_LE(normalizedMatrixError(global, expectedGlobal), kMatrixTolerance); + + Vector localDisplacement{kElementDofCount}; + for (std::size_t index = 0; index < localDisplacement.size(); ++index) { + localDisplacement[index] = 0.01 * static_cast(index + 1U) - 0.04; + } + const Vector globalDisplacement = transposeMultiply(transformation, localDisplacement); + const Vector localForce = local.multiply(localDisplacement); + const Vector globalForce = global.multiply(globalDisplacement); + const Vector expectedGlobalForce = transposeMultiply(transformation, localForce); + for (std::size_t index = 0; index < kElementDofCount; ++index) { + expectScaledNear(globalForce[index], expectedGlobalForce[index], kMatrixTolerance); + } + expectScaledNear( + quadraticEnergy(global, globalDisplacement), + quadraticEnergy(local, localDisplacement), + kMatrixTolerance); + + Vector globalVariation{kElementDofCount}; + for (std::size_t index = 0; index < globalVariation.size(); ++index) { + globalVariation[index] = 0.03 - 0.002 * static_cast(index); + } + const Vector localVariation = transformation.multiply(globalVariation); + expectScaledNear( + globalVariation.dot(globalForce), + localVariation.dot(localForce), + kMatrixTolerance); + + const BeamRecovery recovery = beam.recover(globalDisplacement); + EXPECT_NEAR( + recovery.gaussGeneralizedStrains[0U][0U], + (localDisplacement[6U] - localDisplacement[0U]) / 3.0, + 1.0e-14); +} + +TEST(EulerBeam3D, ConstantLineLoadMatchesAllSignedComponents) { + const double length = 4.0; + const ConstantLocalLineLoad load{2.5, -3.0, 5.5, -7.0}; + const Vector equivalent = alignedBeam(length).localEquivalentLoad(load); + const std::array expected = { + 5.0, -6.0, 11.0, -14.0, -22.0 / 3.0, -4.0, + 5.0, -6.0, 11.0, -14.0, 22.0 / 3.0, 4.0}; + ASSERT_EQ(equivalent.size(), expected.size()); + for (std::size_t index = 0; index < expected.size(); ++index) { + expectScaledNear(equivalent[index], expected[index], kMatrixTolerance); + } +} + +TEST(EulerBeam3D, AnalyticalAxialTorsionAndTwoPlaneBendingRecover) { + const double length = 3.0; + const auto section = makeSection(); + const auto material = makeMaterial(); + const auto beam = alignedBeam(length, section, material); + const Matrix stiffness = beam.localStiffness(); + const double shearModulus = + material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); + + const double axialForce = 1250.0; + const Vector axial = solveFixedFirstNode(stiffness, {axialForce, 0.0, 0.0, 0.0, 0.0, 0.0}); + expectScaledNear( + 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); + + const double torque = -870.0; + const Vector torsion = solveFixedFirstNode(stiffness, {0.0, 0.0, 0.0, torque, 0.0, 0.0}); + expectScaledNear( + 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); + + const double localYForce = 640.0; + const Vector localY = solveFixedFirstNode(stiffness, {0.0, localYForce, 0.0, 0.0, 0.0, 0.0}); + expectScaledNear( + localY[7U], + localYForce * length * length * length / + (3.0 * material.youngsModulus * section.i22), + kAnalyticalTolerance); + expectScaledNear( + 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( + localYRecovery.equilibriumEndActions[0U][5U], + -localYForce * length, + kMatrixTolerance); + expectScaledNear( + localYRecovery.endpointSectionResultants[0U][3U], + localYForce * length, + kMatrixTolerance); + 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( + localZ[8U], + localZForce * length * length * length / + (3.0 * material.youngsModulus * section.i11), + kAnalyticalTolerance); + expectScaledNear( + 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( + localZRecovery.equilibriumEndActions[0U][4U], + localZForce * length, + kMatrixTolerance); + expectScaledNear( + localZRecovery.endpointSectionResultants[0U][2U], + -localZForce * length, + kMatrixTolerance); + EXPECT_NEAR(localZRecovery.endpointSectionResultants[1U][2U], 0.0, 1.0e-8); +} + +TEST(EulerBeam3D, RejectsInvalidGeometryAndProperties) { + const Node origin = makeNode({0.0, 0.0, 0.0}, 1U); + const Node unitX = makeNode({1.0, 0.0, 0.0}, 2U); + const auto validSection = makeSection(); + const auto validMaterial = makeMaterial(); + + const auto expectFailure = [](const Result& result, + const std::string& code) { + ASSERT_FALSE(result.hasValue()); + EXPECT_EQ(result.status().failureCategory(), FailureCategory::model); + ASSERT_EQ(result.status().diagnostics().size(), 1U); + EXPECT_EQ(result.status().diagnostics()[0U].code, code); + }; + + expectFailure( + EulerBeam3D::create(origin, origin, validSection, validMaterial), + "invalid-beam-length"); + expectFailure( + EulerBeam3D::create( + origin, + makeNode({1.0e-12, 0.0, 0.0}, 2U), + validSection, + validMaterial), + "invalid-beam-length"); + + const double coordinate = 1048576.0; + const Node scaledFirst = makeNode({coordinate, 0.0, 0.0}, 1U); + const Node belowThreshold = makeNode({coordinate + coordinate * 0.5e-12, 0.0, 0.0}, 2U); + const Node aboveThreshold = makeNode({coordinate + coordinate * 2.0e-12, 0.0, 0.0}, 2U); + expectFailure( + EulerBeam3D::create(scaledFirst, belowThreshold, validSection, validMaterial), + "invalid-beam-length"); + EXPECT_TRUE(EulerBeam3D::create( + scaledFirst, aboveThreshold, validSection, validMaterial) + .hasValue()); + + auto parallelGuide = validSection; + parallelGuide.firstAxis = {1.0, 0.0, 0.0}; + expectFailure( + EulerBeam3D::create(origin, unitX, parallelGuide, validMaterial), + "invalid-beam-guide-vector"); + auto guideAtThreshold = validSection; + guideAtThreshold.firstAxis = {1.0, 1.0e-12, 0.0}; + expectFailure( + EulerBeam3D::create(origin, unitX, guideAtThreshold, validMaterial), + "invalid-beam-guide-vector"); + auto guideAboveThreshold = validSection; + guideAboveThreshold.firstAxis = {1.0, 2.0e-12, 0.0}; + EXPECT_TRUE(EulerBeam3D::create(origin, unitX, guideAboveThreshold, validMaterial).hasValue()); + + auto invalidMaterial = validMaterial; + invalidMaterial.youngsModulus = 0.0; + expectFailure( + EulerBeam3D::create(origin, unitX, validSection, invalidMaterial), + "invalid-beam-property"); + invalidMaterial = validMaterial; + invalidMaterial.poissonRatio = -2.0; + expectFailure( + EulerBeam3D::create(origin, unitX, validSection, invalidMaterial), + "invalid-beam-property"); + + for (std::size_t property = 0; property < 4U; ++property) { + auto invalidSection = validSection; + double* properties[] = { + &invalidSection.area, + &invalidSection.i11, + &invalidSection.i22, + &invalidSection.torsionalConstant}; + *properties[property] = 0.0; + expectFailure( + EulerBeam3D::create(origin, unitX, invalidSection, validMaterial), + "invalid-beam-property"); + } + + auto coupledSection = validSection; + coupledSection.i12 = 1.0e-9; + expectFailure( + EulerBeam3D::create(origin, unitX, coupledSection, validMaterial), + "unsupported-coupled-section"); +} + +TEST(EulerBeam3D, RecoversSectionPointAndDefaultCentroidS11) { + const double length = 2.0; + const double epsilon = 0.01; + const double kappaY = 0.02; + const double kappaZ = -0.03; + const auto material = makeMaterial(); + auto section = makeSection({0.0, 1.0, 0.0}, {{0.25, -0.5}, {-0.4, 0.3}}); + const auto beam = alignedBeam(length, section, material); + + Vector displacement{kElementDofCount}; + displacement[6U] = epsilon * length; + displacement[7U] = 0.5 * kappaZ * length * length; + displacement[8U] = -0.5 * kappaY * length * length; + displacement[10U] = kappaY * length; + displacement[11U] = kappaZ * length; + + const BeamRecovery recovery = beam.recover(displacement); + ASSERT_EQ(recovery.stressPoints.size(), 4U); + for (std::size_t gaussPoint = 0; gaussPoint < 2U; ++gaussPoint) { + for (std::size_t point = 0; point < section.sectionPoints.size(); ++point) { + const BeamStressPoint& stress = + recovery.stressPoints[gaussPoint * section.sectionPoints.size() + point]; + const double x1 = section.sectionPoints[point][0U]; + const double x2 = section.sectionPoints[point][1U]; + EXPECT_EQ(stress.gaussPoint, static_cast(gaussPoint + 1U)); + EXPECT_EQ(stress.sectionPoint, point + 1U); + EXPECT_DOUBLE_EQ(stress.x1, x1); + EXPECT_DOUBLE_EQ(stress.x2, x2); + expectScaledNear( + stress.s11, + material.youngsModulus * (epsilon + x2 * kappaY - x1 * kappaZ), + kMatrixTolerance); + } + } + + const auto defaultBeam = alignedBeam(length, makeSection(), material); + const BeamRecovery defaultRecovery = defaultBeam.recover(displacement); + ASSERT_EQ(defaultRecovery.stressPoints.size(), 2U); + for (std::size_t gaussPoint = 0; gaussPoint < 2U; ++gaussPoint) { + const BeamStressPoint& stress = defaultRecovery.stressPoints[gaussPoint]; + EXPECT_EQ(stress.gaussPoint, static_cast(gaussPoint + 1U)); + EXPECT_EQ(stress.sectionPoint, 0U); + EXPECT_DOUBLE_EQ(stress.x1, 0.0); + EXPECT_DOUBLE_EQ(stress.x2, 0.0); + EXPECT_EQ(stress.source, "fesa-default"); + expectScaledNear( + stress.s11, + material.youngsModulus * epsilon, + kMatrixTolerance); + } +} + +TEST(EulerBeam3D, ReproducesConstantStrainTwistAndCurvaturePatches) { + const double length = 2.8; + const double epsilon = -0.014; + const double twist = 0.023; + const double kappaY = -0.031; + const double kappaZ = 0.047; + const auto section = makeSection(); + const auto material = makeMaterial(); + const double shearModulus = + material.youngsModulus / (2.0 * (1.0 + material.poissonRatio)); + const auto beam = alignedBeam(length, section, material); + + Vector displacement{kElementDofCount}; + displacement[6U] = epsilon * length; + displacement[7U] = 0.5 * kappaZ * length * length; + displacement[8U] = -0.5 * kappaY * length * length; + displacement[9U] = twist * length; + displacement[10U] = kappaY * length; + displacement[11U] = kappaZ * length; + + const std::array expectedStrain = {epsilon, twist, kappaY, kappaZ}; + const std::array expectedResultant = { + material.youngsModulus * section.area * epsilon, + shearModulus * section.torsionalConstant * twist, + material.youngsModulus * section.i11 * kappaY, + material.youngsModulus * section.i22 * kappaZ}; + const BeamRecovery recovery = beam.recover(displacement); + for (std::size_t point = 0; point < 2U; ++point) { + for (std::size_t component = 0; component < 4U; ++component) { + expectScaledNear( + recovery.gaussGeneralizedStrains[point][component], + expectedStrain[component], + kMatrixTolerance); + expectScaledNear( + recovery.gaussGeneralizedResultants[point][component], + expectedResultant[component], + kMatrixTolerance); + expectScaledNear( + recovery.endpointSectionResultants[point][component], + expectedResultant[component], + kMatrixTolerance); + } + } + + const std::array endActionComponents = {0U, 3U, 4U, 5U}; + for (std::size_t component = 0; component < expectedResultant.size(); ++component) { + expectScaledNear( + recovery.equilibriumEndActions[0U][endActionComponents[component]], + -expectedResultant[component], + kMatrixTolerance); + expectScaledNear( + recovery.equilibriumEndActions[1U][endActionComponents[component]], + expectedResultant[component], + kMatrixTolerance); + } + EXPECT_NEAR(recovery.equilibriumEndActions[0U][1U], 0.0, 1.0e-8); + EXPECT_NEAR(recovery.equilibriumEndActions[0U][2U], 0.0, 1.0e-8); + EXPECT_NEAR(recovery.equilibriumEndActions[1U][1U], 0.0, 1.0e-8); + EXPECT_NEAR(recovery.equilibriumEndActions[1U][2U], 0.0, 1.0e-8); +} + +TEST(EulerBeam3D, OnePointNegativeControlHasRankFour) { + const double length = 3.7; + auto section = makeSection(); + section.area = 1.0; + section.i11 = 0.7; + section.i22 = 1.2; + section.torsionalConstant = 0.9; + const auto material = makeMaterial(4.0, 0.25); + + const Matrix onePoint = testOnlyOnePointStiffness(length, section, material); + const Matrix production = alignedBeam(length, section, material).localStiffness(); + EXPECT_EQ(symmetricRank(onePoint, kRigidTolerance), 4U); + EXPECT_EQ(symmetricRank(production, kRigidTolerance), 6U); +} + +} // namespace +} // namespace fesa