diff --git a/PLAN.md b/PLAN.md index 044a335..9003ecf 100644 --- a/PLAN.md +++ b/PLAN.md @@ -13,7 +13,7 @@ Every new agent session must read this file together with `PROGRESS.md` before p - If an item becomes obsolete, move it to `PROGRESS.md` with a short reason instead of silently deleting it. ## Current Objective -Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-10 MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset. +Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebaseline`, starting with P1R-11 MITC4 patch, locking-sensitivity, and benchmark tests. The old `phases/1-linear-static-mitc4` path is historical and superseded by the paper-based MITC4 formulation reset. ## Required Reading For New Agents 1. `AGENTS.md` @@ -36,7 +36,7 @@ Continue the new Phase 1 rebaseline plan in `phases/1-linear-static-mitc4-rebase ## Active Phase Files - Active phase directory: `phases/1-linear-static-mitc4-rebaseline` - Execute with: `python scripts/execute.py 1-linear-static-mitc4-rebaseline` -- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step15.md` is the independent evaluator closeout. +- Step numbering is zero-based. `step0.md` is complete and recorded in `phases/1-linear-static-mitc4-rebaseline/step0-audit.md`; `step1.md` is complete and created the `quad_02_phase1.inp` normalized reference path; `step2.md` is complete and revalidated core harness guardrails; `step3.md` is complete and revalidated the Phase 1 parser/domain subset; `step4.md` is complete and strengthened validation/singular diagnostics; `step5.md` is complete and revalidated the DofManager/reaction foundation; `step6.md` is complete and revalidated the minimum result model plus displacement CSV comparator; `step7.md` is complete and revalidated MITC4 natural coordinates, tying points, center directors, and integration bases; `step8.md` is complete and revalidated degenerated-continuum displacement, direct covariant strain rows, and MITC shear tying rows; `step9.md` is complete and revalidated plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding; `step10.md` is complete and revalidated MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization; `step15.md` is the independent evaluator closeout. - Every step file contains a sprint contract with objective, required reading, scope, allowed files, explicit non-goals, tests to write first, reference artifacts, acceptance command, evaluator checklist, and handoff requirements. - Historical phase directory: `phases/1-linear-static-mitc4` - Historical phase status: blocked/superseded. Do not resume the old P1-15/P1-16 path unless the user explicitly requests recovery of that exact phase. @@ -76,7 +76,7 @@ Each gate should be satisfied before moving to the next implementation band unle | G1 - Build and validation | satisfied | Build system, test framework, and `scripts/validate_workspace.py` run real checks. | Validation command output | | G2 - Parser and domain | satisfied | Parser subset revalidated in step 3; validation and singular diagnostics revalidated in step 4. | Parser acceptance/rejection tests, validation negative tests, and validation output | | G3 - DOF/math/results infrastructure | partial | Core aliases, DOF mapping, validation harness, model diagnostic context, DofManager, sparse-connectivity inputs, full-vector reaction formula, result model metadata, and displacement CSV comparator were revalidated in steps 2, 5, and 6; assembly remains for step 12. | P1R-02, P1R-05, and P1R-06 validation output | -| G4 - MITC4 element readiness | partial | MITC4 formulation was rewritten from local papers; Steps 7 through 9 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, MITC shear tying rows, plane-stress material, convected-to-local transform, and `2 x 2 x 2` material integration scaffolding, while stiffness/drilling and patch benchmarks remain for steps 10 and 11. | P1R-07 through P1R-09 validation output; future element tests | +| G4 - MITC4 element readiness | partial | MITC4 formulation was rewritten from local papers; Steps 7 through 10 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, MITC shear tying rows, plane-stress material, convected-to-local transform, `2 x 2 x 2` material integration scaffolding, stiffness/internal force, six-DOF transform, and drilling stabilization, while patch/benchmark coverage remains for step 11. | P1R-07 through P1R-10 validation output; future element tests | | G5 - End-to-end solver | reopened | Linear static path must be revalidated through steps 13 and 14 after the MITC4 rebuild and `quad_02` compatibility path. | Future integration/reference regression output | ## Phase 1 Implementation Milestones @@ -91,7 +91,7 @@ All milestones are intended to become one or more self-contained sprint contract | P1R-07 | completed | MITC4 generator | Implement MITC4 geometry, node order, tying points, directors, and local bases. | none | Shape/basis/diagnostic tests | | P1R-08 | completed | MITC4 generator | Implement degenerated-continuum displacement, covariant strain rows, and MITC shear tying. | P1R-07 | Finite-difference and tying interpolation tests | | P1R-09 | completed | MITC4 generator | Implement material matrix, transform, and `2 x 2 x 2` integration scaffolding. | P1R-08 | Material/integration tests | -| P1R-10 | pending | MITC4 generator | Assemble MITC4 stiffness/internal force with six-DOF transform and drilling stabilization. | P1R-09, P1R-05 | Symmetry, rigid body, drilling sensitivity tests | +| P1R-10 | completed | MITC4 generator | Assemble MITC4 stiffness/internal force with six-DOF transform and drilling stabilization. | P1R-09, P1R-05 | Symmetry, rigid body, drilling sensitivity tests | | P1R-11 | pending | verification generator | Add MITC4 patch, locking-sensitivity, and benchmark tests. | P1R-10 | Membrane/bending/shear/twist/locking tests | | P1R-12 | pending | assembly generator | Rebuild assembly, solver adapter boundary, constrained solve, and full-vector RF recovery. | P1R-05, P1R-10 | Assembly and full-vector reaction tests | | P1R-13 | pending | analysis generator | Rebuild linear static workflow from input to U/RF result fields. | P1R-03, P1R-04, P1R-06, P1R-12 | End-to-end linear static tests | diff --git a/PROGRESS.md b/PROGRESS.md index 2c60639..0c30401 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,37 @@ Every new agent session must read this file together with `PLAN.md` before plann - Do not remove history unless the user explicitly asks for archival cleanup. ## Current Status -Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 9 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, and MITC4 material/transform/integration scaffolding have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. +Phase 1 has a new rebaseline phase definition in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 10 are complete. `quad_02_phase1.inp` is now the normalized Phase 1-compatible input path for the stored `quad_02` S4 reference pair, while the original `quad_02.inp` remains preserved unsupported provenance. Core numeric aliases, DOF mapping, validation harness, model diagnostic context, the Phase 1 parser/domain subset, validation/singular diagnostics, DofManager/reaction foundation, minimum result model metadata, displacement CSV comparator foundation, MITC4 geometry/director scaffolding, MITC4 displacement/strain/tying row scaffolding, MITC4 material/transform/integration scaffolding, and MITC4 stiffness/drilling/internal-force scaffolding have been revalidated. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. ## Completed Work +### 2026-05-04 - P1R-10 MITC4 stiffness and drilling completed +Author: Codex + +Changed files: +- `include/fesa/fesa.hpp` +- `tests/test_main.cpp` +- `phases/1-linear-static-mitc4-rebaseline/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Added tests and implementation for MITC4 element stiffness accumulation using tied strain rows, transformed material matrices, and `2 x 2 x 2` integration samples. +- Added a stiffness result API that exposes local physical stiffness, local drilling-stabilized stiffness, global stiffness, integration-point count, drilling reference diagonal, and drilling stiffness metadata. +- Replaced the old default drilling scale with the documented `drilling_stiffness_scale = 1.0e-3` and applied it to the minimum positive physical local stiffness diagonal before global transformation. +- Added diagnostics for invalid drilling scale and missing positive drilling reference diagonal. +- Added tests for stiffness symmetry, drilling scale sensitivity, rigid translation zero energy, documented drilling-only energy, and `f_int_e = K_e u_e` consistency. +- Routed `MITC4ElementKernel::stiffness` through the rebuilt formulation path so later assembly uses the Step 7 through Step 10 MITC4 helpers. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding Step 10 tests; it failed as expected because the MITC4 stiffness/drilling/internal-force APIs did not exist yet. +- After implementation, `python scripts/validate_workspace.py` configured CMake, built `fesa_core` and `fesa_tests`, and ran CTest successfully. +- CTest result: 1 test executable passed. + +Follow-up: +- Continue with P1R-11 MITC4 patch, locking-sensitivity, and benchmark tests. +- Step 10 intentionally did not tune the drilling scale to a reference case, add stress/resultant output, introduce reduced integration/hourglass logic, or rebuild sparse/global assembly beyond routing the existing kernel to the new element formulation. + ### 2026-05-04 - P1R-09 MITC4 material integration scaffolding completed Author: Codex diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index 001af04..58df44c 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -1885,189 +1885,217 @@ inline LocalBasis computeLocalBasis(const std::array& coordinates) { return {frame.v1, frame.v2, frame.vn}; } -struct NaturalDerivatives { - ShapeData shape; - std::array dx{}; - std::array dy{}; - Real det_j = 0.0; +struct ElementStiffnessOptions { + Real drilling_stiffness_scale = 1.0e-3; }; -inline NaturalDerivatives naturalDerivatives(const std::array, 4>& xy, Real r, Real s) { - NaturalDerivatives data; - data.shape = shapeFunctions(r, s); - Real j00 = 0.0; - Real j01 = 0.0; - Real j10 = 0.0; - Real j11 = 0.0; - for (std::size_t i = 0; i < 4; ++i) { - j00 += data.shape.dr[i] * xy[i][0]; - j01 += data.shape.ds[i] * xy[i][0]; - j10 += data.shape.dr[i] * xy[i][1]; - j11 += data.shape.ds[i] * xy[i][1]; +struct MITC4DrillingStabilizationResult { + DenseMatrix local_with_drilling; + Real reference_diagonal = 0.0; + Real drilling_stiffness = 0.0; + Real drilling_stiffness_scale = 0.0; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); } - data.det_j = j00 * j11 - j01 * j10; - if (std::fabs(data.det_j) < 1.0e-14) { - throw std::runtime_error("invalid shell element jacobian"); +}; + +struct MITC4ElementStiffnessResult { + DenseMatrix local_without_drilling; + DenseMatrix local_with_drilling; + DenseMatrix global; + std::size_t integration_point_count = 0; + Real drilling_reference_diagonal = 0.0; + Real drilling_stiffness = 0.0; + Real drilling_stiffness_scale = 0.0; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); } - const Real inv00 = j11 / data.det_j; - const Real inv01 = -j01 / data.det_j; - const Real inv10 = -j10 / data.det_j; - const Real inv11 = j00 / data.det_j; - for (std::size_t i = 0; i < 4; ++i) { - data.dx[i] = inv00 * data.shape.dr[i] + inv01 * data.shape.ds[i]; - data.dy[i] = inv10 * data.shape.dr[i] + inv11 * data.shape.ds[i]; +}; + +inline DenseMatrix mitc4LocalDofTransform(const MITC4Geometry& geometry) { + DenseMatrix transform(24, 24); + for (LocalIndex node = 0; node < 4; ++node) { + const LocalIndex base = 6 * node; + const auto& frame = geometry.nodal_frames[static_cast(node)]; + const std::array axes = {frame.v1, frame.v2, frame.vn}; + for (LocalIndex local_axis = 0; local_axis < 3; ++local_axis) { + const Vec3 axis = axes[static_cast(local_axis)]; + transform(base + local_axis, base + 0) = axis.x; + transform(base + local_axis, base + 1) = axis.y; + transform(base + local_axis, base + 2) = axis.z; + transform(base + 3 + local_axis, base + 3) = axis.x; + transform(base + 3 + local_axis, base + 4) = axis.y; + transform(base + 3 + local_axis, base + 5) = axis.z; + } } - return data; + return transform; } -struct ElementStiffnessOptions { - Real drilling_stiffness_scale = 1.0e-6; -}; +inline MITC4StrainRows mitc4TransformStrainRowsToLocalDofs(const MITC4StrainRows& global_rows, + const DenseMatrix& local_dof_transform) { + MITC4StrainRows local_rows; + local_rows.diagnostics = global_rows.diagnostics; + if (hasError(local_rows.diagnostics)) { + return local_rows; + } + for (std::size_t component = 0; component < 6; ++component) { + for (LocalIndex local_dof = 0; local_dof < 24; ++local_dof) { + Real value = 0.0; + for (LocalIndex global_dof = 0; global_dof < 24; ++global_dof) { + value += global_rows.rows[component][static_cast(global_dof)] * + local_dof_transform(local_dof, global_dof); + } + local_rows.rows[component][static_cast(local_dof)] = value; + } + } + return local_rows; +} + +inline void accumulateMITC4BtDB(DenseMatrix& stiffness, const MITC4StrainRows& rows, + const MITC4MaterialMatrix& material, Real factor) { + for (LocalIndex i = 0; i < 24; ++i) { + for (LocalIndex j = 0; j < 24; ++j) { + Real value = 0.0; + for (std::size_t a = 0; a < 6; ++a) { + for (std::size_t b = 0; b < 6; ++b) { + value += rows.rows[a][static_cast(i)] * material[a][b] * + rows.rows[b][static_cast(j)]; + } + } + stiffness.add(i, j, value * factor); + } + } +} + +inline Real mitc4MinimumPositivePhysicalLocalDiagonal(const DenseMatrix& local_without_drilling, + Real tolerance = 1.0e-12) { + Real minimum = std::numeric_limits::infinity(); + for (LocalIndex node = 0; node < 4; ++node) { + for (LocalIndex local_dof = 0; local_dof < 5; ++local_dof) { + const Real diagonal = local_without_drilling(6 * node + local_dof, 6 * node + local_dof); + if (isFinite(diagonal) && diagonal > tolerance) { + minimum = std::min(minimum, diagonal); + } + } + } + return minimum; +} + +inline MITC4DrillingStabilizationResult mitc4ApplyDrillingStabilization(const DenseMatrix& local_without_drilling, + Real drilling_stiffness_scale, + Real tolerance = 1.0e-12) { + MITC4DrillingStabilizationResult result; + result.local_with_drilling = local_without_drilling; + result.drilling_stiffness_scale = drilling_stiffness_scale; + if (!isFinite(drilling_stiffness_scale) || drilling_stiffness_scale < 0.0) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-DRILLING-SCALE", "MITC4 drilling stiffness scale must be non-negative and finite")); + return result; + } + + result.reference_diagonal = mitc4MinimumPositivePhysicalLocalDiagonal(local_without_drilling, tolerance); + if (!isFinite(result.reference_diagonal)) { + result.diagnostics.push_back(mitc4Diagnostic("FESA-MITC4-DRILLING-REFERENCE", + "MITC4 drilling stiffness reference diagonal is not positive")); + return result; + } + + result.drilling_stiffness = drilling_stiffness_scale * result.reference_diagonal; + for (LocalIndex node = 0; node < 4; ++node) { + const LocalIndex gamma = 6 * node + 5; + result.local_with_drilling.add(gamma, gamma, result.drilling_stiffness); + } + return result; +} + +inline DenseMatrix mitc4TransformLocalStiffnessToGlobal(const DenseMatrix& local_stiffness, + const DenseMatrix& local_dof_transform) { + DenseMatrix global(24, 24); + for (LocalIndex i = 0; i < 24; ++i) { + for (LocalIndex j = 0; j < 24; ++j) { + Real value = 0.0; + for (LocalIndex a = 0; a < 24; ++a) { + for (LocalIndex b = 0; b < 24; ++b) { + value += local_dof_transform(a, i) * local_stiffness(a, b) * local_dof_transform(b, j); + } + } + global(i, j) = value; + } + } + return global; +} + +inline MITC4ElementStiffnessResult mitc4ElementStiffness(const std::array& coordinates, + Real elastic_modulus, Real poisson_ratio, Real thickness, + ElementStiffnessOptions options = {}) { + MITC4ElementStiffnessResult result; + result.local_without_drilling = DenseMatrix(24, 24); + result.local_with_drilling = DenseMatrix(24, 24); + result.global = DenseMatrix(24, 24); + result.drilling_stiffness_scale = options.drilling_stiffness_scale; + + const MITC4Geometry geometry = buildMITC4Geometry(coordinates, thickness); + appendDiagnostics(result.diagnostics, geometry.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + + const MITC4MaterialIntegrationData integration = + mitc4BuildMaterialIntegrationData(geometry, elastic_modulus, poisson_ratio); + appendDiagnostics(result.diagnostics, integration.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + result.integration_point_count = integration.samples.size(); + + const DenseMatrix local_dof_transform = mitc4LocalDofTransform(geometry); + for (const MITC4MaterialIntegrationSample& sample : integration.samples) { + const MITC4StrainRows local_rows = + mitc4TransformStrainRowsToLocalDofs(sample.strain_rows, local_dof_transform); + appendDiagnostics(result.diagnostics, local_rows.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + const Real factor = std::fabs(sample.basis.jacobian) * sample.point.weight; + accumulateMITC4BtDB(result.local_without_drilling, local_rows, sample.convected_material, factor); + } + + const auto drilling = mitc4ApplyDrillingStabilization(result.local_without_drilling, options.drilling_stiffness_scale); + appendDiagnostics(result.diagnostics, drilling.diagnostics); + result.local_with_drilling = drilling.local_with_drilling; + result.drilling_reference_diagonal = drilling.reference_diagonal; + result.drilling_stiffness = drilling.drilling_stiffness; + result.drilling_stiffness_scale = drilling.drilling_stiffness_scale; + if (hasError(result.diagnostics)) { + return result; + } + + result.global = mitc4TransformLocalStiffnessToGlobal(result.local_with_drilling, local_dof_transform); + return result; +} + +inline std::vector mitc4ElementInternalForce(const MITC4ElementStiffnessResult& stiffness, + const std::vector& element_displacement) { + if (stiffness.global.rows() != static_cast(element_displacement.size())) { + throw std::runtime_error("MITC4 internal force size mismatch"); + } + return stiffness.global.multiply(element_displacement); +} class MITC4ElementKernel { public: DenseMatrix stiffness(const std::array& coordinates, Real elastic_modulus, Real poisson_ratio, Real thickness, ElementStiffnessOptions options = {}) const { - const LocalBasis basis = computeLocalBasis(coordinates); - std::array, 4> xy{}; - for (std::size_t i = 0; i < 4; ++i) { - const Vec3 relative = coordinates[i] - coordinates[0]; - xy[i] = {dot(relative, basis.e1), dot(relative, basis.e2)}; + const auto result = mitc4ElementStiffness(coordinates, elastic_modulus, poisson_ratio, thickness, options); + if (!result.ok()) { + throw std::runtime_error(result.diagnostics.empty() ? "invalid MITC4 stiffness" + : result.diagnostics.front().message); } - - DenseMatrix local(24, 24); - const Real membrane_factor = elastic_modulus * thickness / (1.0 - poisson_ratio * poisson_ratio); - const Real bending_factor = elastic_modulus * thickness * thickness * thickness / (12.0 * (1.0 - poisson_ratio * poisson_ratio)); - const Real shear_factor = (5.0 / 6.0) * elastic_modulus * thickness / (2.0 * (1.0 + poisson_ratio)); - const Real gauss = 1.0 / std::sqrt(3.0); - const std::array points = {-gauss, gauss}; - - for (Real r : points) { - for (Real s : points) { - const NaturalDerivatives d = naturalDerivatives(xy, r, s); - DenseMatrix b(8, 24); - addMembraneB(b, d); - addBendingB(b, d); - addMitcShearB(b, xy, r, s); - accumulateBtDB(local, b, membrane_factor, bending_factor, shear_factor, poisson_ratio, std::fabs(d.det_j)); - addDrilling(local, d, elastic_modulus * thickness * options.drilling_stiffness_scale, std::fabs(d.det_j)); - } - } - - return transformToGlobal(local, basis); - } - - private: - static void addMembraneB(DenseMatrix& b, const NaturalDerivatives& d) { - for (LocalIndex i = 0; i < 4; ++i) { - const LocalIndex c = 6 * i; - b(0, c + 0) = d.dx[static_cast(i)]; - b(1, c + 1) = d.dy[static_cast(i)]; - b(2, c + 0) = d.dy[static_cast(i)]; - b(2, c + 1) = d.dx[static_cast(i)]; - } - } - - static void addBendingB(DenseMatrix& b, const NaturalDerivatives& d) { - for (LocalIndex i = 0; i < 4; ++i) { - const LocalIndex c = 6 * i; - b(3, c + 4) = -d.dx[static_cast(i)]; - b(4, c + 3) = d.dy[static_cast(i)]; - b(5, c + 3) = d.dx[static_cast(i)]; - b(5, c + 4) = -d.dy[static_cast(i)]; - } - } - - static void addStandardShearRow(DenseMatrix& b, LocalIndex row, const NaturalDerivatives& d, bool gamma_xz, Real scale) { - for (LocalIndex i = 0; i < 4; ++i) { - const LocalIndex c = 6 * i; - if (gamma_xz) { - b(row, c + 2) += scale * d.dx[static_cast(i)]; - b(row, c + 4) += scale * d.shape.n[static_cast(i)]; - } else { - b(row, c + 2) += scale * d.dy[static_cast(i)]; - b(row, c + 3) -= scale * d.shape.n[static_cast(i)]; - } - } - } - - static void addMitcShearB(DenseMatrix& b, const std::array, 4>& xy, Real r, Real s) { - addStandardShearRow(b, 6, naturalDerivatives(xy, 0.0, -1.0), true, 0.5 * (1.0 - s)); - addStandardShearRow(b, 6, naturalDerivatives(xy, 0.0, 1.0), true, 0.5 * (1.0 + s)); - addStandardShearRow(b, 7, naturalDerivatives(xy, -1.0, 0.0), false, 0.5 * (1.0 - r)); - addStandardShearRow(b, 7, naturalDerivatives(xy, 1.0, 0.0), false, 0.5 * (1.0 + r)); - } - - static void accumulateBtDB(DenseMatrix& k, const DenseMatrix& b, Real membrane_factor, Real bending_factor, Real shear_factor, Real poisson_ratio, Real det_j) { - std::array, 8> d{}; - d[0][0] = membrane_factor; - d[0][1] = poisson_ratio * membrane_factor; - d[1][0] = poisson_ratio * membrane_factor; - d[1][1] = membrane_factor; - d[2][2] = membrane_factor * (1.0 - poisson_ratio) / 2.0; - d[3][3] = bending_factor; - d[3][4] = poisson_ratio * bending_factor; - d[4][3] = poisson_ratio * bending_factor; - d[4][4] = bending_factor; - d[5][5] = bending_factor * (1.0 - poisson_ratio) / 2.0; - d[6][6] = shear_factor; - d[7][7] = shear_factor; - for (LocalIndex i = 0; i < 24; ++i) { - for (LocalIndex j = 0; j < 24; ++j) { - Real value = 0.0; - for (LocalIndex row = 0; row < 8; ++row) { - for (LocalIndex col = 0; col < 8; ++col) { - value += b(row, i) * d[static_cast(row)][static_cast(col)] * b(col, j); - } - } - k.add(i, j, value * det_j); - } - } - } - - static void addDrilling(DenseMatrix& k, const NaturalDerivatives& d, Real drill_factor, Real det_j) { - std::array b{}; - for (LocalIndex i = 0; i < 4; ++i) { - const LocalIndex c = 6 * i; - b[static_cast(c + 0)] = -0.5 * d.dy[static_cast(i)]; - b[static_cast(c + 1)] = 0.5 * d.dx[static_cast(i)]; - b[static_cast(c + 5)] = -d.shape.n[static_cast(i)]; - } - for (LocalIndex i = 0; i < 24; ++i) { - for (LocalIndex j = 0; j < 24; ++j) { - k.add(i, j, b[static_cast(i)] * drill_factor * b[static_cast(j)] * det_j); - } - } - } - - static DenseMatrix transformToGlobal(const DenseMatrix& local, const LocalBasis& basis) { - DenseMatrix transform(24, 24); - const std::array axes = {basis.e1, basis.e2, basis.e3}; - for (LocalIndex node = 0; node < 4; ++node) { - for (LocalIndex local_axis = 0; local_axis < 3; ++local_axis) { - const Vec3 axis = axes[static_cast(local_axis)]; - const LocalIndex base = 6 * node; - transform(base + local_axis, base + 0) = axis.x; - transform(base + local_axis, base + 1) = axis.y; - transform(base + local_axis, base + 2) = axis.z; - transform(base + 3 + local_axis, base + 3) = axis.x; - transform(base + 3 + local_axis, base + 4) = axis.y; - transform(base + 3 + local_axis, base + 5) = axis.z; - } - } - DenseMatrix global(24, 24); - for (LocalIndex i = 0; i < 24; ++i) { - for (LocalIndex j = 0; j < 24; ++j) { - Real value = 0.0; - for (LocalIndex a = 0; a < 24; ++a) { - for (LocalIndex b = 0; b < 24; ++b) { - value += transform(a, i) * local(a, b) * transform(b, j); - } - } - global(i, j) = value; - } - } - return global; + return result.global; } }; diff --git a/phases/1-linear-static-mitc4-rebaseline/index.json b/phases/1-linear-static-mitc4-rebaseline/index.json index a4c2154..594ff69 100644 --- a/phases/1-linear-static-mitc4-rebaseline/index.json +++ b/phases/1-linear-static-mitc4-rebaseline/index.json @@ -12,7 +12,7 @@ { "step": 7, "name": "mitc4-geometry-directors", "status": "completed" }, { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "completed" }, { "step": 9, "name": "mitc4-material-integration", "status": "completed" }, - { "step": 10, "name": "mitc4-stiffness-drilling", "status": "pending" }, + { "step": 10, "name": "mitc4-stiffness-drilling", "status": "completed" }, { "step": 11, "name": "mitc4-patch-benchmark-tests", "status": "pending" }, { "step": 12, "name": "assembly-sparse-solver-path", "status": "pending" }, { "step": 13, "name": "linear-static-workflow", "status": "pending" }, diff --git a/tests/test_main.cpp b/tests/test_main.cpp index 3ab4c92..09c992a 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -147,6 +147,19 @@ std::array zeroElementDofs() { return values; } +std::vector zeroElementVector() { + return std::vector(24, 0.0); +} + +fesa::Real quadraticEnergy(const fesa::DenseMatrix& stiffness, const std::vector& values) { + const auto internal = stiffness.multiply(values); + fesa::Real energy = 0.0; + for (std::size_t i = 0; i < values.size(); ++i) { + energy += values[i] * internal[i]; + } + return energy; +} + fesa::Domain singleElementValidationDomain() { fesa::Domain domain; domain.nodes[1] = {1, {0, 0, 0}}; @@ -1219,6 +1232,106 @@ FESA_TEST(mitc4_material_integration_samples_carry_tied_rows_basis_and_transform } } +FESA_TEST(mitc4_element_stiffness_result_is_symmetric_and_uses_2x2x2_samples) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2); + FESA_CHECK(result.ok()); + FESA_CHECK(result.integration_point_count == 8); + FESA_CHECK(result.local_without_drilling.rows() == 24); + FESA_CHECK(result.local_with_drilling.rows() == 24); + FESA_CHECK(result.global.rows() == 24); + FESA_CHECK(result.global.cols() == 24); + for (fesa::LocalIndex i = 0; i < 24; ++i) { + for (fesa::LocalIndex j = 0; j < 24; ++j) { + FESA_CHECK_NEAR(result.local_without_drilling(i, j), result.local_without_drilling(j, i), 1.0e-8); + FESA_CHECK_NEAR(result.local_with_drilling(i, j), result.local_with_drilling(j, i), 1.0e-8); + FESA_CHECK_NEAR(result.global(i, j), result.global(j, i), 1.0e-8); + } + } +} + +FESA_TEST(mitc4_drilling_stiffness_uses_minimum_positive_physical_local_diagonal) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2); + FESA_CHECK(result.ok()); + fesa::Real expected_reference = std::numeric_limits::infinity(); + for (fesa::LocalIndex node = 0; node < 4; ++node) { + for (fesa::LocalIndex local_dof = 0; local_dof < 5; ++local_dof) { + const fesa::Real diagonal = result.local_without_drilling(6 * node + local_dof, 6 * node + local_dof); + if (std::isfinite(diagonal) && diagonal > 0.0) { + expected_reference = std::min(expected_reference, diagonal); + } + } + } + FESA_CHECK(std::isfinite(expected_reference)); + FESA_CHECK_NEAR(result.drilling_reference_diagonal, expected_reference, 1.0e-10); + FESA_CHECK_NEAR(result.drilling_stiffness_scale, 1.0e-3, 1.0e-15); + FESA_CHECK_NEAR(result.drilling_stiffness, 1.0e-3 * expected_reference, 1.0e-10); + for (fesa::LocalIndex node = 0; node < 4; ++node) { + const fesa::LocalIndex gamma = 6 * node + 5; + FESA_CHECK_NEAR(result.local_with_drilling(gamma, gamma) - result.local_without_drilling(gamma, gamma), + result.drilling_stiffness, 1.0e-12); + } + + fesa::ElementStiffnessOptions options; + options.drilling_stiffness_scale = 2.0e-3; + const auto scaled = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2, options); + FESA_CHECK(scaled.ok()); + FESA_CHECK_NEAR(scaled.drilling_reference_diagonal, result.drilling_reference_diagonal, 1.0e-10); + FESA_CHECK_NEAR(scaled.drilling_stiffness, 2.0 * result.drilling_stiffness, 1.0e-10); +} + +FESA_TEST(mitc4_drilling_reports_invalid_scale_and_missing_reference_diagonal) { + fesa::DenseMatrix zero_physical(24, 24); + auto missing_reference = fesa::mitc4ApplyDrillingStabilization(zero_physical, 1.0e-3); + FESA_CHECK(!missing_reference.ok()); + FESA_CHECK(fesa::containsDiagnostic(missing_reference.diagnostics, "FESA-MITC4-DRILLING-REFERENCE")); + + auto invalid_scale = fesa::mitc4ApplyDrillingStabilization(zero_physical, -1.0); + FESA_CHECK(!invalid_scale.ok()); + FESA_CHECK(fesa::containsDiagnostic(invalid_scale.diagnostics, "FESA-MITC4-DRILLING-SCALE")); +} + +FESA_TEST(mitc4_rigid_translation_energy_is_zero_and_drilling_energy_is_documented) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2); + FESA_CHECK(result.ok()); + + auto translate_x = zeroElementVector(); + auto translate_y = zeroElementVector(); + auto translate_z = zeroElementVector(); + for (int node = 0; node < 4; ++node) { + translate_x[static_cast(6 * node + 0)] = 1.0; + translate_y[static_cast(6 * node + 1)] = 1.0; + translate_z[static_cast(6 * node + 2)] = 1.0; + } + FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_x)) < 1.0e-8); + FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_y)) < 1.0e-8); + FESA_CHECK(std::fabs(quadraticEnergy(result.global, translate_z)) < 1.0e-8); + + auto drilling = zeroElementVector(); + for (int node = 0; node < 4; ++node) { + drilling[static_cast(6 * node + 5)] = 1.0; + } + FESA_CHECK_NEAR(quadraticEnergy(result.global, drilling), 4.0 * result.drilling_stiffness, 1.0e-8); +} + +FESA_TEST(mitc4_internal_force_is_stiffness_times_displacement_for_linear_phase1) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto result = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2); + FESA_CHECK(result.ok()); + std::vector displacement(24, 0.0); + for (std::size_t i = 0; i < displacement.size(); ++i) { + displacement[i] = 0.001 * static_cast(i + 1); + } + const auto expected = result.global.multiply(displacement); + const auto actual = fesa::mitc4ElementInternalForce(result, displacement); + FESA_CHECK(actual.size() == expected.size()); + for (std::size_t i = 0; i < expected.size(); ++i) { + FESA_CHECK_NEAR(actual[i], expected[i], 1.0e-12); + } +} + FESA_TEST(mitc4_shape_functions_and_stiffness_baseline) { auto shape = fesa::shapeFunctions(0.25, -0.5); const fesa::Real sum = shape.n[0] + shape.n[1] + shape.n[2] + shape.n[3];