diff --git a/PLAN.md b/PLAN.md index 6182af6..598df2e 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-08 degenerated-continuum displacement, covariant strain rows, and MITC shear tying. 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-09 material matrix, transform, and `2 x 2 x 2` integration scaffolding. 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; `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; `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; Step 7 rebuilt geometry/director/local-basis scaffolding, while strain tying, material/integration, stiffness/drilling, and patch benchmarks remain for steps 8 through 11. | P1R-07 validation output; future element tests | +| G4 - MITC4 element readiness | partial | MITC4 formulation was rewritten from local papers; Steps 7 and 8 rebuilt geometry/director/local-basis scaffolding, displacement interpolation, direct covariant strain rows, and MITC shear tying rows, while material/integration, stiffness/drilling, and patch benchmarks remain for steps 9 through 11. | P1R-07 and P1R-08 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 @@ -89,7 +89,7 @@ All milestones are intended to become one or more self-contained sprint contract | P1R-05 | completed | DOF generator | Rebuild six-DOF DofManager, constrained/free mapping, equation numbering, and full-vector reconstruction. | none | DOF mapping and reaction foundation tests | | P1R-06 | completed | results generator | Rebuild minimum results model and displacement CSV comparator. | none | U/RF schema tests and CSV comparator tests | | P1R-07 | completed | MITC4 generator | Implement MITC4 geometry, node order, tying points, directors, and local bases. | none | Shape/basis/diagnostic tests | -| P1R-08 | pending | MITC4 generator | Implement degenerated-continuum displacement, covariant strain rows, and MITC shear tying. | P1R-07 | Finite-difference and tying interpolation 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 | pending | 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-11 | pending | verification generator | Add MITC4 patch, locking-sensitivity, and benchmark tests. | P1R-10 | Membrane/bending/shear/twist/locking tests | diff --git a/PROGRESS.md b/PROGRESS.md index 6e2e1b1..9453556 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,36 @@ 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 7 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, and MITC4 geometry/director 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 8 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, and MITC4 displacement/strain/tying row 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-08 MITC4 covariant strain tying 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 MITC4 tests for global rotation to local `alpha`, `beta`, `gamma` mapping and verified drilling `gamma` does not enter the physical director increment. +- Added degenerated-continuum displacement interpolation tests at the midsurface and through-thickness points, including rotation-driven `q_k = -V2_k alpha_k + V1_k beta_k` behavior. +- Added direct covariant strain-row finite-difference tests in the documented order `[eps11, eps22, eps33, gamma23, gamma13, gamma12]`. +- Added MITC tying tests proving FESA `A/C` signs for `gamma13` and `B/D` signs for `gamma23`, and proving Gauss-point shear rows are interpolated from tying rows rather than using direct Gauss-point transverse shear. +- Implemented reusable MITC4 displacement derivative, direct strain evaluation, direct strain-row, tied strain-row, and row-evaluation helpers for later material/integration and stiffness steps. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding Step 8 tests; it failed as expected because the MITC4 displacement/strain/tying 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-09 material matrix, transform, and `2 x 2 x 2` integration scaffolding. +- The existing stiffness kernel still needs to be rebuilt against the new row APIs in steps 9 and 10; Step 8 intentionally did not integrate stiffness, add material transforms, or add S4R behavior. + ### 2026-05-04 - P1R-07 MITC4 geometry and directors completed Author: Codex diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index ef20ddd..e75d0d2 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -1301,6 +1301,64 @@ struct MITC4IntegrationBasis { } }; +using MITC4ElementDofVector = std::array; +using MITC4StrainVector = std::array; +using MITC4StrainRow = std::array; + +enum class MITC4StrainComponent { + Eps11 = 0, + Eps22 = 1, + Eps33 = 2, + Gamma23 = 3, + Gamma13 = 4, + Gamma12 = 5 +}; + +inline std::size_t strainComponentIndex(MITC4StrainComponent component) { + return static_cast(component); +} + +inline std::array mitc4StrainComponentLabels() { + return {"eps11", "eps22", "eps33", "gamma23", "gamma13", "gamma12"}; +} + +struct MITC4LocalRotations { + Real alpha = 0.0; + Real beta = 0.0; + Real gamma = 0.0; +}; + +struct MITC4DisplacementDerivatives { + ShapeData shape; + Vec3 displacement; + Vec3 du_dxi; + Vec3 du_deta; + Vec3 du_dzeta; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4StrainEvaluation { + MITC4StrainVector values{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4StrainRows { + std::array rows{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + inline Vec3 globalEX() { return {1.0, 0.0, 0.0}; } @@ -1317,6 +1375,10 @@ inline Diagnostic mitc4Diagnostic(std::string code, std::string message) { return makeDiagnostic(Severity::Error, std::move(code), std::move(message), "mitc4", "", 0); } +inline void appendDiagnostics(std::vector& target, const std::vector& source) { + target.insert(target.end(), source.begin(), source.end()); +} + inline MITC4MidsurfaceDerivatives mitc4MidsurfaceDerivatives(const std::array& coordinates, Real xi, Real eta) { MITC4MidsurfaceDerivatives result; result.shape = shapeFunctions(xi, eta); @@ -1436,6 +1498,132 @@ inline MITC4IntegrationBasis computeMITC4IntegrationBasis(const MITC4Geometry& g return result; } +inline Vec3 mitc4NodalTranslation(const MITC4ElementDofVector& values, std::size_t node) { + const std::size_t base = 6 * node; + return {values[base + 0], values[base + 1], values[base + 2]}; +} + +inline Vec3 mitc4NodalRotation(const MITC4ElementDofVector& values, std::size_t node) { + const std::size_t base = 6 * node; + return {values[base + 3], values[base + 4], values[base + 5]}; +} + +inline MITC4LocalRotations mitc4LocalRotations(const MITC4DirectorFrame& frame, const Vec3& global_rotation) { + return {dot(global_rotation, frame.v1), dot(global_rotation, frame.v2), dot(global_rotation, frame.vn)}; +} + +inline Vec3 mitc4DirectorIncrement(const MITC4DirectorFrame& frame, const Vec3& global_rotation) { + const MITC4LocalRotations rotations = mitc4LocalRotations(frame, global_rotation); + return (-rotations.alpha) * frame.v2 + rotations.beta * frame.v1; +} + +inline MITC4DisplacementDerivatives mitc4DisplacementDerivatives(const MITC4Geometry& geometry, const MITC4ElementDofVector& values, + Real xi, Real eta, Real zeta) { + MITC4DisplacementDerivatives result; + result.diagnostics = geometry.diagnostics; + result.shape = shapeFunctions(xi, eta); + for (std::size_t node = 0; node < 4; ++node) { + const Vec3 translation = mitc4NodalTranslation(values, node); + const Vec3 rotation = mitc4NodalRotation(values, node); + const Vec3 q = mitc4DirectorIncrement(geometry.nodal_frames[node], rotation); + const Real n = result.shape.n[node]; + const Real dn_dxi = result.shape.dr[node]; + const Real dn_deta = result.shape.ds[node]; + result.displacement = result.displacement + n * translation + (0.5 * zeta * geometry.thickness * n) * q; + result.du_dxi = result.du_dxi + dn_dxi * translation + (0.5 * zeta * geometry.thickness * dn_dxi) * q; + result.du_deta = result.du_deta + dn_deta * translation + (0.5 * zeta * geometry.thickness * dn_deta) * q; + result.du_dzeta = result.du_dzeta + (0.5 * geometry.thickness * n) * q; + } + return result; +} + +inline void assignMITC4CovariantStrain(MITC4StrainVector& values, const MITC4IntegrationBasis& basis, + const MITC4DisplacementDerivatives& derivatives) { + values[strainComponentIndex(MITC4StrainComponent::Eps11)] = dot(derivatives.du_dxi, basis.g1); + values[strainComponentIndex(MITC4StrainComponent::Eps22)] = dot(derivatives.du_deta, basis.g2); + values[strainComponentIndex(MITC4StrainComponent::Eps33)] = 0.0; + values[strainComponentIndex(MITC4StrainComponent::Gamma23)] = dot(derivatives.du_deta, basis.g3) + dot(derivatives.du_dzeta, basis.g2); + values[strainComponentIndex(MITC4StrainComponent::Gamma13)] = dot(derivatives.du_dxi, basis.g3) + dot(derivatives.du_dzeta, basis.g1); + values[strainComponentIndex(MITC4StrainComponent::Gamma12)] = dot(derivatives.du_dxi, basis.g2) + dot(derivatives.du_deta, basis.g1); +} + +inline MITC4StrainEvaluation mitc4DirectCovariantStrain(const MITC4Geometry& geometry, const MITC4ElementDofVector& values, + Real xi, Real eta, Real zeta) { + MITC4StrainEvaluation result; + const auto basis = computeMITC4IntegrationBasis(geometry, xi, eta, zeta); + appendDiagnostics(result.diagnostics, basis.diagnostics); + const auto derivatives = mitc4DisplacementDerivatives(geometry, values, xi, eta, zeta); + appendDiagnostics(result.diagnostics, derivatives.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + assignMITC4CovariantStrain(result.values, basis, derivatives); + return result; +} + +inline MITC4StrainRows mitc4DirectCovariantStrainRows(const MITC4Geometry& geometry, Real xi, Real eta, Real zeta) { + MITC4StrainRows result; + const auto basis = computeMITC4IntegrationBasis(geometry, xi, eta, zeta); + appendDiagnostics(result.diagnostics, basis.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + for (std::size_t dof = 0; dof < 24; ++dof) { + MITC4ElementDofVector unit{}; + unit.fill(0.0); + unit[dof] = 1.0; + const auto derivatives = mitc4DisplacementDerivatives(geometry, unit, xi, eta, zeta); + appendDiagnostics(result.diagnostics, derivatives.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + MITC4StrainVector values{}; + assignMITC4CovariantStrain(values, basis, derivatives); + for (std::size_t component = 0; component < 6; ++component) { + result.rows[component][dof] = values[component]; + } + } + return result; +} + +inline MITC4StrainEvaluation evaluateMITC4StrainRows(const MITC4StrainRows& rows, const MITC4ElementDofVector& values) { + MITC4StrainEvaluation result; + result.diagnostics = rows.diagnostics; + if (hasError(result.diagnostics)) { + return result; + } + for (std::size_t component = 0; component < 6; ++component) { + for (std::size_t dof = 0; dof < 24; ++dof) { + result.values[component] += rows.rows[component][dof] * values[dof]; + } + } + return result; +} + +inline MITC4StrainRows mitc4TiedCovariantStrainRows(const MITC4Geometry& geometry, Real xi, Real eta, Real zeta) { + MITC4StrainRows result = mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta); + const auto direct_a = mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, zeta); + const auto direct_b = mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, zeta); + const auto direct_c = mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, zeta); + const auto direct_d = mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, zeta); + appendDiagnostics(result.diagnostics, direct_a.diagnostics); + appendDiagnostics(result.diagnostics, direct_b.diagnostics); + appendDiagnostics(result.diagnostics, direct_c.diagnostics); + appendDiagnostics(result.diagnostics, direct_d.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + const std::size_t gamma23 = strainComponentIndex(MITC4StrainComponent::Gamma23); + const std::size_t gamma13 = strainComponentIndex(MITC4StrainComponent::Gamma13); + for (std::size_t dof = 0; dof < 24; ++dof) { + result.rows[gamma13][dof] = 0.5 * (1.0 - eta) * direct_a.rows[gamma13][dof] + + 0.5 * (1.0 + eta) * direct_c.rows[gamma13][dof]; + result.rows[gamma23][dof] = 0.5 * (1.0 - xi) * direct_b.rows[gamma23][dof] + + 0.5 * (1.0 + xi) * direct_d.rows[gamma23][dof]; + } + return result; +} + inline LocalBasis computeLocalBasis(const std::array& coordinates) { const MITC4Geometry geometry = buildMITC4Geometry(coordinates, 1.0); if (!geometry.ok()) { diff --git a/phases/1-linear-static-mitc4-rebaseline/index.json b/phases/1-linear-static-mitc4-rebaseline/index.json index b5994a3..755b2bc 100644 --- a/phases/1-linear-static-mitc4-rebaseline/index.json +++ b/phases/1-linear-static-mitc4-rebaseline/index.json @@ -10,7 +10,7 @@ { "step": 5, "name": "dof-manager-reaction-foundation", "status": "completed" }, { "step": 6, "name": "results-comparator-foundation", "status": "completed" }, { "step": 7, "name": "mitc4-geometry-directors", "status": "completed" }, - { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "pending" }, + { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "completed" }, { "step": 9, "name": "mitc4-material-integration", "status": "pending" }, { "step": 10, "name": "mitc4-stiffness-drilling", "status": "pending" }, { "step": 11, "name": "mitc4-patch-benchmark-tests", "status": "pending" }, diff --git a/tests/test_main.cpp b/tests/test_main.cpp index 80d8650..d04136e 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -1,5 +1,6 @@ #include "fesa/fesa.hpp" +#include #include #include #include @@ -140,6 +141,12 @@ void checkRightHandedOrthonormal(const fesa::Vec3& e1, const fesa::Vec3& e2, con FESA_CHECK_NEAR(fesa::dot(fesa::cross(e1, e2), e3), 1.0, 1.0e-14); } +std::array zeroElementDofs() { + std::array values{}; + values.fill(0.0); + return values; +} + fesa::Domain singleElementValidationDomain() { fesa::Domain domain; domain.nodes[1] = {1, {0, 0, 0}}; @@ -921,6 +928,148 @@ FESA_TEST(mitc4_geometry_reports_singular_geometry_and_thickness) { FESA_CHECK(fesa::containsDiagnostic(corner_basis.diagnostics, "FESA-MITC4-SINGULAR-JACOBIAN")); } +FESA_TEST(mitc4_rotation_transform_and_displacement_interpolation) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + FESA_CHECK(geometry.ok()); + const auto flat_rotation = fesa::mitc4LocalRotations(geometry.nodal_frames[0], {1.0, 2.0, 3.0}); + FESA_CHECK_NEAR(flat_rotation.alpha, 1.0, 1.0e-14); + FESA_CHECK_NEAR(flat_rotation.beta, 2.0, 1.0e-14); + FESA_CHECK_NEAR(flat_rotation.gamma, 3.0, 1.0e-14); + checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 0.0, 5.0}), {0, 0, 0}, 1.0e-14); + checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 2.0, 5.0}), {2, 0, 0}, 1.0e-14); + + std::array dofs = zeroElementDofs(); + for (std::size_t node = 0; node < 4; ++node) { + dofs[6 * node + 0] = 1.0; + dofs[6 * node + 1] = 0.5; + dofs[6 * node + 4] = 2.0; + dofs[6 * node + 5] = 99.0; + } + const auto mid = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 0.0); + FESA_CHECK(mid.ok()); + checkVecNear(mid.displacement, {1.0, 0.5, 0.0}, 1.0e-14); + const auto top = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 1.0); + FESA_CHECK(top.ok()); + checkVecNear(top.displacement, {1.2, 0.5, 0.0}, 1.0e-14); + + const std::array fallback_coords = {{{0, 0, 0}, {1, 0, 0}, {1, 0, -1}, {0, 0, -1}}}; + const auto fallback_geometry = fesa::buildMITC4Geometry(fallback_coords, 0.1); + FESA_CHECK(fallback_geometry.ok()); + const auto fallback_rotation = fesa::mitc4LocalRotations(fallback_geometry.nodal_frames[0], {2.0, 3.0, 5.0}); + FESA_CHECK_NEAR(fallback_rotation.alpha, -2.0, 1.0e-14); + FESA_CHECK_NEAR(fallback_rotation.beta, 5.0, 1.0e-14); + FESA_CHECK_NEAR(fallback_rotation.gamma, 3.0, 1.0e-14); +} + +FESA_TEST(mitc4_direct_covariant_strain_rows_match_finite_difference) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + FESA_CHECK(geometry.ok()); + constexpr fesa::Real xi = 0.2; + constexpr fesa::Real eta = -0.3; + constexpr fesa::Real zeta = 0.4; + const auto rows = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta); + FESA_CHECK(rows.ok()); + FESA_CHECK(fesa::mitc4StrainComponentLabels()[0] == "eps11"); + FESA_CHECK(fesa::mitc4StrainComponentLabels()[3] == "gamma23"); + const std::array checked_dofs = {0, 1, 2, 7, 9, 10, 11}; + const fesa::Real h = 1.0e-6; + for (std::size_t dof : checked_dofs) { + auto plus = zeroElementDofs(); + auto minus = zeroElementDofs(); + plus[dof] = h; + minus[dof] = -h; + const auto plus_strain = fesa::mitc4DirectCovariantStrain(geometry, plus, xi, eta, zeta); + const auto minus_strain = fesa::mitc4DirectCovariantStrain(geometry, minus, xi, eta, zeta); + FESA_CHECK(plus_strain.ok()); + FESA_CHECK(minus_strain.ok()); + for (std::size_t component = 0; component < 6; ++component) { + const fesa::Real finite_difference = (plus_strain.values[component] - minus_strain.values[component]) / (2.0 * h); + FESA_CHECK_NEAR(rows.rows[component][dof], finite_difference, 1.0e-8); + } + } + + const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13); + const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23); + for (std::size_t node = 0; node < 4; ++node) { + const std::size_t drilling_dof = 6 * node + 5; + for (std::size_t component = 0; component < 6; ++component) { + FESA_CHECK_NEAR(rows.rows[component][drilling_dof], 0.0, 1.0e-14); + } + } + FESA_CHECK(gamma13 == 4); + FESA_CHECK(gamma23 == 3); +} + +FESA_TEST(mitc4_mitc_tying_rows_use_fesa_ac_bd_signs) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + FESA_CHECK(geometry.ok()); + const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23); + const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13); + + const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, 0.0); + const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, 0.0); + const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, 0.0); + const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, 0.0); + const auto mitc_a = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, -1.0, 0.0); + const auto mitc_b = fesa::mitc4TiedCovariantStrainRows(geometry, -1.0, 0.0, 0.0); + const auto mitc_c = fesa::mitc4TiedCovariantStrainRows(geometry, 0.0, 1.0, 0.0); + const auto mitc_d = fesa::mitc4TiedCovariantStrainRows(geometry, 1.0, 0.0, 0.0); + FESA_CHECK(direct_a.ok()); + FESA_CHECK(direct_b.ok()); + FESA_CHECK(direct_c.ok()); + FESA_CHECK(direct_d.ok()); + FESA_CHECK(mitc_a.ok()); + FESA_CHECK(mitc_b.ok()); + FESA_CHECK(mitc_c.ok()); + FESA_CHECK(mitc_d.ok()); + for (std::size_t dof = 0; dof < 24; ++dof) { + FESA_CHECK_NEAR(mitc_a.rows[gamma13][dof], direct_a.rows[gamma13][dof], 1.0e-14); + FESA_CHECK_NEAR(mitc_c.rows[gamma13][dof], direct_c.rows[gamma13][dof], 1.0e-14); + FESA_CHECK_NEAR(mitc_b.rows[gamma23][dof], direct_b.rows[gamma23][dof], 1.0e-14); + FESA_CHECK_NEAR(mitc_d.rows[gamma23][dof], direct_d.rows[gamma23][dof], 1.0e-14); + } +} + +FESA_TEST(mitc4_mitc_gauss_shear_rows_are_interpolated_from_tying_rows) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + FESA_CHECK(geometry.ok()); + constexpr fesa::Real xi = 0.5; + constexpr fesa::Real eta = 0.25; + constexpr fesa::Real zeta = 0.0; + const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23); + const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13); + const auto direct_at_point = fesa::mitc4DirectCovariantStrainRows(geometry, xi, eta, zeta); + const auto direct_a = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, -1.0, zeta); + const auto direct_b = fesa::mitc4DirectCovariantStrainRows(geometry, -1.0, 0.0, zeta); + const auto direct_c = fesa::mitc4DirectCovariantStrainRows(geometry, 0.0, 1.0, zeta); + const auto direct_d = fesa::mitc4DirectCovariantStrainRows(geometry, 1.0, 0.0, zeta); + const auto tied = fesa::mitc4TiedCovariantStrainRows(geometry, xi, eta, zeta); + FESA_CHECK(direct_at_point.ok()); + FESA_CHECK(tied.ok()); + + const std::size_t node2_ry = 6 + 4; + const fesa::Real expected_gamma13 = + 0.5 * (1.0 - eta) * direct_a.rows[gamma13][node2_ry] + 0.5 * (1.0 + eta) * direct_c.rows[gamma13][node2_ry]; + FESA_CHECK_NEAR(tied.rows[gamma13][node2_ry], expected_gamma13, 1.0e-14); + FESA_CHECK(std::fabs(tied.rows[gamma13][node2_ry] - direct_at_point.rows[gamma13][node2_ry]) > 1.0e-4); + + const std::size_t node4_rx = 3 * 6 + 3; + const fesa::Real expected_gamma23 = + 0.5 * (1.0 - xi) * direct_b.rows[gamma23][node4_rx] + 0.5 * (1.0 + xi) * direct_d.rows[gamma23][node4_rx]; + FESA_CHECK_NEAR(tied.rows[gamma23][node4_rx], expected_gamma23, 1.0e-14); + FESA_CHECK(std::fabs(tied.rows[gamma23][node4_rx] - direct_at_point.rows[gamma23][node4_rx]) > 1.0e-4); + + for (std::size_t component : {std::size_t{0}, std::size_t{1}, std::size_t{2}, std::size_t{5}}) { + for (std::size_t dof = 0; dof < 24; ++dof) { + FESA_CHECK_NEAR(tied.rows[component][dof], direct_at_point.rows[component][dof], 1.0e-14); + } + } +} + 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];