From fa492f994d8cf80db81563346d90448b4e16013e Mon Sep 17 00:00:00 2001 From: NINI Date: Mon, 4 May 2026 13:38:20 +0900 Subject: [PATCH] feat: add mitc4 geometry directors --- PLAN.md | 8 +- PROGRESS.md | 28 ++- include/fesa/fesa.hpp | 228 +++++++++++++++++- .../index.json | 2 +- tests/test_main.cpp | 121 ++++++++++ 5 files changed, 373 insertions(+), 14 deletions(-) diff --git a/PLAN.md b/PLAN.md index 16d3677..6182af6 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-07 MITC4 geometry, node order, tying points, directors, and local bases. 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-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. ## 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; `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; `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 | reopened | MITC4 formulation was rewritten from local papers; element implementation must be rebuilt or revalidated through steps 7 through 11. | Revised `docs/MITC4_FORMULATION.md`, future element tests | +| 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 | | 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 @@ -88,7 +88,7 @@ All milestones are intended to become one or more self-contained sprint contract | P1R-04 | completed | validation generator | Rebuild validation and singular diagnostic coverage. | P1R-03 | Missing-reference and singular-prone negative tests | | 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 | pending | MITC4 generator | Implement MITC4 geometry, node order, tying points, directors, and local bases. | none | Shape/basis/diagnostic 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-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 | diff --git a/PROGRESS.md b/PROGRESS.md index 326afc7..6e2e1b1 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 6 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, and displacement CSV comparator foundation 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 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. ## Completed Work +### 2026-05-04 - P1R-07 MITC4 geometry and directors 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 shape-function derivatives, FESA/Abaqus S4 natural-coordinate node order, tying point labels `A`, `B`, `C`, and `D`, and edge-node mapping. +- Added Phase 1 MITC4 geometry structures for natural points, tying points, center midsurface derivatives, center-normal director policy, nodal director frames, and integration-point local Cartesian bases. +- Implemented the documented center director policy `Vn = normalize(G1_c x G2_c)`, nodal `V1 = normalize(EY x Vn)` with deterministic fallback axes, and right-handed `V2 = Vn x V1`. +- Added integration-basis construction from degenerated-continuum covariant basis vectors and diagnostics for invalid thickness, singular center normal, singular basis, and near-zero Jacobian. +- Replaced the legacy `computeLocalBasis()` internals with the new center-director geometry policy while leaving stiffness/strain/drilling reimplementation to later steps. + +Verification: +- First ran `python scripts/validate_workspace.py` after adding Step 7 tests; it failed as expected because the MITC4 geometry/director API 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-08 degenerated-continuum displacement, covariant strain rows, and MITC shear tying. +- The existing stiffness kernel still needs to be rebuilt against the new formulation in steps 8 through 10; Step 7 intentionally did not add stiffness, stress, strain, or S4R behavior. + ### 2026-05-04 - P1R-06 results model and displacement CSV comparator completed Author: Codex diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index 849fbbb..ef20ddd 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -198,6 +198,22 @@ inline Real norm(const Vec3& value) { return std::sqrt(dot(value, value)); } +inline bool isFinite(Real value) { + return std::isfinite(value); +} + +inline bool isFinite(const Vec3& value) { + return isFinite(value.x) && isFinite(value.y) && isFinite(value.z); +} + +inline std::optional normalizedIfValid(const Vec3& value, Real tolerance = 1.0e-12) { + const Real length = norm(value); + if (!isFinite(length) || length <= tolerance) { + return std::nullopt; + } + return (1.0 / length) * value; +} + inline Vec3 normalized(const Vec3& value) { const Real length = norm(value); if (length <= std::numeric_limits::epsilon()) { @@ -1222,15 +1238,211 @@ struct LocalBasis { Vec3 e3; }; +struct MITC4NaturalPoint { + Real xi = 0.0; + Real eta = 0.0; +}; + +struct MITC4TyingPoint { + std::string label; + MITC4NaturalPoint natural; + std::array edge_node_indices{}; +}; + +inline std::array mitc4NodeNaturalCoordinates() { + return {{{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}}}; +} + +inline std::array mitc4TyingPoints() { + return {{{"A", {0.0, -1.0}, {0, 1}}, + {"B", {-1.0, 0.0}, {0, 3}}, + {"C", {0.0, 1.0}, {3, 2}}, + {"D", {1.0, 0.0}, {1, 2}}}}; +} + +struct MITC4DirectorFrame { + Vec3 v1; + Vec3 v2; + Vec3 vn; +}; + +struct MITC4MidsurfaceDerivatives { + ShapeData shape; + Vec3 g1; + Vec3 g2; +}; + +struct MITC4Geometry { + std::array coordinates{}; + Real thickness = 0.0; + ShapeData center_shape; + Vec3 g1_center; + Vec3 g2_center; + Vec3 center_normal; + std::array nodal_frames{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4IntegrationBasis { + ShapeData shape; + Vec3 g1; + Vec3 g2; + Vec3 g3; + Real jacobian = 0.0; + LocalBasis local; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +inline Vec3 globalEX() { + return {1.0, 0.0, 0.0}; +} + +inline Vec3 globalEY() { + return {0.0, 1.0, 0.0}; +} + +inline Vec3 globalEZ() { + return {0.0, 0.0, 1.0}; +} + +inline Diagnostic mitc4Diagnostic(std::string code, std::string message) { + return makeDiagnostic(Severity::Error, std::move(code), std::move(message), "mitc4", "", 0); +} + +inline MITC4MidsurfaceDerivatives mitc4MidsurfaceDerivatives(const std::array& coordinates, Real xi, Real eta) { + MITC4MidsurfaceDerivatives result; + result.shape = shapeFunctions(xi, eta); + for (std::size_t i = 0; i < 4; ++i) { + result.g1 = result.g1 + result.shape.dr[i] * coordinates[i]; + result.g2 = result.g2 + result.shape.ds[i] * coordinates[i]; + } + return result; +} + +inline std::optional firstNormalizedCross(const std::array& axes, const Vec3& vector, Real tolerance) { + for (const Vec3& axis : axes) { + auto candidate = normalizedIfValid(cross(axis, vector), tolerance); + if (candidate) { + return candidate; + } + } + return std::nullopt; +} + +inline std::optional buildMITC4DirectorFrame(const Vec3& normal, Real tolerance) { + auto v1 = normalizedIfValid(cross(globalEY(), normal), tolerance); + if (!v1) { + v1 = firstNormalizedCross({globalEZ(), globalEX(), globalEY()}, normal, tolerance); + } + if (!v1) { + return std::nullopt; + } + auto v2 = normalizedIfValid(cross(normal, *v1), tolerance); + if (!v2) { + return std::nullopt; + } + return MITC4DirectorFrame{*v1, *v2, normal}; +} + +inline MITC4Geometry buildMITC4Geometry(const std::array& coordinates, Real thickness, Real tolerance = 1.0e-12) { + MITC4Geometry geometry; + geometry.coordinates = coordinates; + geometry.thickness = thickness; + if (!isFinite(thickness) || thickness <= tolerance) { + geometry.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-THICKNESS", "MITC4 shell thickness must be positive and finite")); + } + for (const Vec3& coordinate : coordinates) { + if (!isFinite(coordinate)) { + geometry.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-COORDINATE", "MITC4 element coordinates must be finite")); + break; + } + } + + const auto center = mitc4MidsurfaceDerivatives(coordinates, 0.0, 0.0); + geometry.center_shape = center.shape; + geometry.g1_center = center.g1; + geometry.g2_center = center.g2; + const auto normal = normalizedIfValid(cross(center.g1, center.g2), tolerance); + if (!normal) { + geometry.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-NORMAL", "MITC4 element center normal is near zero")); + return geometry; + } + geometry.center_normal = *normal; + + const auto frame = buildMITC4DirectorFrame(*normal, tolerance); + if (!frame) { + geometry.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 nodal director basis could not be constructed")); + return geometry; + } + geometry.nodal_frames.fill(*frame); + return geometry; +} + +inline MITC4IntegrationBasis computeMITC4IntegrationBasis(const MITC4Geometry& geometry, Real xi, Real eta, Real zeta, + Real tolerance = 1.0e-12) { + MITC4IntegrationBasis result; + result.diagnostics = geometry.diagnostics; + result.shape = shapeFunctions(xi, eta); + for (std::size_t i = 0; i < 4; ++i) { + const Vec3& coordinate = geometry.coordinates[i]; + const Vec3& normal = geometry.nodal_frames[i].vn; + result.g1 = result.g1 + result.shape.dr[i] * coordinate + + (0.5 * zeta * geometry.thickness * result.shape.dr[i]) * normal; + result.g2 = result.g2 + result.shape.ds[i] * coordinate + + (0.5 * zeta * geometry.thickness * result.shape.ds[i]) * normal; + result.g3 = result.g3 + (0.5 * geometry.thickness * result.shape.n[i]) * normal; + } + + result.jacobian = dot(cross(result.g1, result.g2), result.g3); + if (!isFinite(result.jacobian) || std::fabs(result.jacobian) <= tolerance) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 element Jacobian is near zero")); + } + + const auto e3 = normalizedIfValid(result.g3, tolerance); + if (!e3) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis normal is near zero")); + return result; + } + auto e1 = normalizedIfValid(cross(result.g2, *e3), tolerance); + if (!e1) { + e1 = firstNormalizedCross({globalEY(), globalEZ(), globalEX()}, *e3, tolerance); + } + if (!e1) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis tangent could not be constructed")); + return result; + } + const auto e2 = normalizedIfValid(cross(*e3, *e1), tolerance); + if (!e2) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-BASIS", "MITC4 integration basis is not right-handed")); + return result; + } + result.local = {*e1, *e2, *e3}; + return result; +} + inline LocalBasis computeLocalBasis(const std::array& coordinates) { - Vec3 v1 = 0.5 * ((coordinates[1] - coordinates[0]) + (coordinates[2] - coordinates[3])); - Vec3 v2 = 0.5 * ((coordinates[3] - coordinates[0]) + (coordinates[2] - coordinates[1])); - Vec3 e1 = normalized(v1); - v2 = v2 - dot(v2, e1) * e1; - Vec3 e2 = normalized(v2); - Vec3 e3 = normalized(cross(e1, e2)); - e2 = normalized(cross(e3, e1)); - return {e1, e2, e3}; + const MITC4Geometry geometry = buildMITC4Geometry(coordinates, 1.0); + if (!geometry.ok()) { + throw std::runtime_error("invalid MITC4 geometry"); + } + const MITC4DirectorFrame& frame = geometry.nodal_frames[0]; + return {frame.v1, frame.v2, frame.vn}; } struct NaturalDerivatives { diff --git a/phases/1-linear-static-mitc4-rebaseline/index.json b/phases/1-linear-static-mitc4-rebaseline/index.json index 4c64c6e..b5994a3 100644 --- a/phases/1-linear-static-mitc4-rebaseline/index.json +++ b/phases/1-linear-static-mitc4-rebaseline/index.json @@ -9,7 +9,7 @@ { "step": 4, "name": "validation-singular-diagnostics", "status": "completed" }, { "step": 5, "name": "dof-manager-reaction-foundation", "status": "completed" }, { "step": 6, "name": "results-comparator-foundation", "status": "completed" }, - { "step": 7, "name": "mitc4-geometry-directors", "status": "pending" }, + { "step": 7, "name": "mitc4-geometry-directors", "status": "completed" }, { "step": 8, "name": "mitc4-covariant-strain-tying", "status": "pending" }, { "step": 9, "name": "mitc4-material-integration", "status": "pending" }, { "step": 10, "name": "mitc4-stiffness-drilling", "status": "pending" }, diff --git a/tests/test_main.cpp b/tests/test_main.cpp index 857967e..80d8650 100644 --- a/tests/test_main.cpp +++ b/tests/test_main.cpp @@ -124,6 +124,22 @@ std::size_t diagnosticCount(const std::vector& diagnostics, co return count; } +void checkVecNear(const fesa::Vec3& actual, const fesa::Vec3& expected, fesa::Real tolerance) { + FESA_CHECK_NEAR(actual.x, expected.x, tolerance); + FESA_CHECK_NEAR(actual.y, expected.y, tolerance); + FESA_CHECK_NEAR(actual.z, expected.z, tolerance); +} + +void checkRightHandedOrthonormal(const fesa::Vec3& e1, const fesa::Vec3& e2, const fesa::Vec3& e3) { + FESA_CHECK_NEAR(fesa::norm(e1), 1.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::norm(e2), 1.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::norm(e3), 1.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::dot(e1, e2), 0.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::dot(e1, e3), 0.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::dot(e2, e3), 0.0, 1.0e-14); + FESA_CHECK_NEAR(fesa::dot(fesa::cross(e1, e2), e3), 1.0, 1.0e-14); +} + fesa::Domain singleElementValidationDomain() { fesa::Domain domain; domain.nodes[1] = {1, {0, 0, 0}}; @@ -800,6 +816,111 @@ FESA_TEST(displacement_comparator_reports_duplicate_actual_nodes) { FESA_CHECK(fesa::containsDiagnostic(compared.diagnostics, "FESA-COMPARE-DUPLICATE-ACTUAL")); } +FESA_TEST(mitc4_shape_functions_node_order_and_tying_points) { + auto center = fesa::shapeFunctions(0.0, 0.0); + const fesa::Real sum = center.n[0] + center.n[1] + center.n[2] + center.n[3]; + FESA_CHECK_NEAR(sum, 1.0, 1.0e-15); + FESA_CHECK_NEAR(center.dr[0], -0.25, 1.0e-15); + FESA_CHECK_NEAR(center.dr[1], 0.25, 1.0e-15); + FESA_CHECK_NEAR(center.dr[2], 0.25, 1.0e-15); + FESA_CHECK_NEAR(center.dr[3], -0.25, 1.0e-15); + FESA_CHECK_NEAR(center.ds[0], -0.25, 1.0e-15); + FESA_CHECK_NEAR(center.ds[1], -0.25, 1.0e-15); + FESA_CHECK_NEAR(center.ds[2], 0.25, 1.0e-15); + FESA_CHECK_NEAR(center.ds[3], 0.25, 1.0e-15); + + const auto node_points = fesa::mitc4NodeNaturalCoordinates(); + FESA_CHECK_NEAR(node_points[0].xi, -1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[0].eta, -1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[1].xi, 1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[1].eta, -1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[2].xi, 1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[2].eta, 1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[3].xi, -1.0, 1.0e-15); + FESA_CHECK_NEAR(node_points[3].eta, 1.0, 1.0e-15); + + for (std::size_t i = 0; i < 4; ++i) { + const auto corner = fesa::shapeFunctions(node_points[i].xi, node_points[i].eta); + FESA_CHECK_NEAR(corner.n[i], 1.0, 1.0e-15); + } + + const auto tying_points = fesa::mitc4TyingPoints(); + FESA_CHECK(tying_points[0].label == "A"); + FESA_CHECK_NEAR(tying_points[0].natural.xi, 0.0, 1.0e-15); + FESA_CHECK_NEAR(tying_points[0].natural.eta, -1.0, 1.0e-15); + FESA_CHECK((tying_points[0].edge_node_indices == std::array{0, 1})); + FESA_CHECK(tying_points[1].label == "B"); + FESA_CHECK_NEAR(tying_points[1].natural.xi, -1.0, 1.0e-15); + FESA_CHECK_NEAR(tying_points[1].natural.eta, 0.0, 1.0e-15); + FESA_CHECK((tying_points[1].edge_node_indices == std::array{0, 3})); + FESA_CHECK(tying_points[2].label == "C"); + FESA_CHECK_NEAR(tying_points[2].natural.xi, 0.0, 1.0e-15); + FESA_CHECK_NEAR(tying_points[2].natural.eta, 1.0, 1.0e-15); + FESA_CHECK((tying_points[2].edge_node_indices == std::array{3, 2})); + FESA_CHECK(tying_points[3].label == "D"); + FESA_CHECK_NEAR(tying_points[3].natural.xi, 1.0, 1.0e-15); + FESA_CHECK_NEAR(tying_points[3].natural.eta, 0.0, 1.0e-15); + FESA_CHECK((tying_points[3].edge_node_indices == std::array{1, 2})); +} + +FESA_TEST(mitc4_geometry_builds_flat_directors_and_integration_basis) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + FESA_CHECK(geometry.ok()); + checkVecNear(geometry.center_normal, {0, 0, 1}, 1.0e-14); + checkVecNear(geometry.g1_center, {0.5, 0, 0}, 1.0e-14); + checkVecNear(geometry.g2_center, {0, 0.5, 0}, 1.0e-14); + for (const auto& frame : geometry.nodal_frames) { + checkVecNear(frame.v1, {1, 0, 0}, 1.0e-14); + checkVecNear(frame.v2, {0, 1, 0}, 1.0e-14); + checkVecNear(frame.vn, {0, 0, 1}, 1.0e-14); + checkRightHandedOrthonormal(frame.v1, frame.v2, frame.vn); + } + + auto basis = fesa::computeMITC4IntegrationBasis(geometry, 0.0, 0.0, 0.0); + FESA_CHECK(basis.ok()); + checkVecNear(basis.g1, {0.5, 0, 0}, 1.0e-14); + checkVecNear(basis.g2, {0, 0.5, 0}, 1.0e-14); + checkVecNear(basis.g3, {0, 0, 0.1}, 1.0e-14); + checkVecNear(basis.local.e1, {1, 0, 0}, 1.0e-14); + checkVecNear(basis.local.e2, {0, 1, 0}, 1.0e-14); + checkVecNear(basis.local.e3, {0, 0, 1}, 1.0e-14); + checkRightHandedOrthonormal(basis.local.e1, basis.local.e2, basis.local.e3); + FESA_CHECK_NEAR(basis.jacobian, 0.025, 1.0e-14); +} + +FESA_TEST(mitc4_geometry_uses_deterministic_director_fallback_axis) { + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 0, -1}, {0, 0, -1}}}; + auto geometry = fesa::buildMITC4Geometry(coords, 0.1); + FESA_CHECK(geometry.ok()); + checkVecNear(geometry.center_normal, {0, 1, 0}, 1.0e-14); + for (const auto& frame : geometry.nodal_frames) { + checkVecNear(frame.v1, {-1, 0, 0}, 1.0e-14); + checkVecNear(frame.v2, {0, 0, 1}, 1.0e-14); + checkVecNear(frame.vn, {0, 1, 0}, 1.0e-14); + checkRightHandedOrthonormal(frame.v1, frame.v2, frame.vn); + } +} + +FESA_TEST(mitc4_geometry_reports_singular_geometry_and_thickness) { + const std::array flat = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + auto invalid_thickness = fesa::buildMITC4Geometry(flat, 0.0); + FESA_CHECK(!invalid_thickness.ok()); + FESA_CHECK(fesa::containsDiagnostic(invalid_thickness.diagnostics, "FESA-MITC4-THICKNESS")); + + const std::array collinear = {{{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}}}; + auto singular = fesa::buildMITC4Geometry(collinear, 0.1); + FESA_CHECK(!singular.ok()); + FESA_CHECK(fesa::containsDiagnostic(singular.diagnostics, "FESA-MITC4-SINGULAR-NORMAL")); + + const std::array collapsed_corner = {{{0, 0, 0}, {0, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + auto geometry = fesa::buildMITC4Geometry(collapsed_corner, 0.1); + FESA_CHECK(geometry.ok()); + auto corner_basis = fesa::computeMITC4IntegrationBasis(geometry, -1.0, -1.0, 0.0); + FESA_CHECK(!corner_basis.ok()); + FESA_CHECK(fesa::containsDiagnostic(corner_basis.diagnostics, "FESA-MITC4-SINGULAR-JACOBIAN")); +} + 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];