From 150653c3c762c80f77d19246cad6f1a95ec71d45 Mon Sep 17 00:00:00 2001 From: NINI Date: Tue, 5 May 2026 23:11:23 +0900 Subject: [PATCH] refactor: extract mitc4 geometry strain helpers --- CMakeLists.txt | 6 + PLAN.md | 6 +- PROGRESS.md | 36 +- include/fesa/Element/Element.hpp | 2 + include/fesa/Element/MITC4Geometry.hpp | 265 +++++++++++ include/fesa/Element/MITC4Kinematics.hpp | 205 +++++++++ include/fesa/fesa.hpp | 430 +----------------- .../1-structure-alignment-refactor/index.json | 2 +- tests/test_element_module_includes.cpp | 173 +++++++ 9 files changed, 691 insertions(+), 434 deletions(-) create mode 100644 include/fesa/Element/MITC4Geometry.hpp create mode 100644 include/fesa/Element/MITC4Kinematics.hpp create mode 100644 tests/test_element_module_includes.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5487834..105f2f8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,6 +32,9 @@ add_executable(fesa_results_module_tests tests/test_results_module_includes.cpp) target_link_libraries(fesa_results_module_tests PRIVATE fesa_core) target_compile_definitions(fesa_results_module_tests PRIVATE FESA_SOURCE_DIR="${CMAKE_CURRENT_SOURCE_DIR}") +add_executable(fesa_element_module_tests tests/test_element_module_includes.cpp) +target_link_libraries(fesa_element_module_tests PRIVATE fesa_core) + if(MSVC) target_compile_options(fesa_core PRIVATE /W4 /permissive-) target_compile_options(fesa_tests PRIVATE /W4 /permissive-) @@ -39,6 +42,7 @@ if(MSVC) target_compile_options(fesa_math_module_tests PRIVATE /W4 /permissive-) target_compile_options(fesa_io_module_tests PRIVATE /W4 /permissive-) target_compile_options(fesa_results_module_tests PRIVATE /W4 /permissive-) + target_compile_options(fesa_element_module_tests PRIVATE /W4 /permissive-) else() target_compile_options(fesa_core PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_tests PRIVATE -Wall -Wextra -Wpedantic) @@ -46,6 +50,7 @@ else() target_compile_options(fesa_math_module_tests PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_io_module_tests PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_results_module_tests PRIVATE -Wall -Wextra -Wpedantic) + target_compile_options(fesa_element_module_tests PRIVATE -Wall -Wextra -Wpedantic) endif() add_test(NAME fesa_tests COMMAND fesa_tests) @@ -53,3 +58,4 @@ add_test(NAME fesa_core_module_tests COMMAND fesa_core_module_tests) add_test(NAME fesa_math_module_tests COMMAND fesa_math_module_tests) add_test(NAME fesa_io_module_tests COMMAND fesa_io_module_tests) add_test(NAME fesa_results_module_tests COMMAND fesa_results_module_tests) +add_test(NAME fesa_element_module_tests COMMAND fesa_element_module_tests) diff --git a/PLAN.md b/PLAN.md index 67d9de7..7116702 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 -Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignment-refactor`, continuing with P1A-06 MITC4 geometry/strain extraction. P1A-00 completed the architecture drift audit, P1A-01 created the module scaffold, P1A-02 extracted Core/Util plus Phase 1 Boundary/Load/Property model ownership, P1A-03 extracted Math primitives plus the solver adapter boundary, P1A-04 extracted the Abaqus Phase 1 parser into IO, and P1A-05 extracted Results/reference comparison without changing solver behavior. This phase must align the current monolithic `include/fesa/fesa.hpp` implementation with the module ownership model in `docs/ARCHITECTURE.md` without changing solver behavior. Product-level Phase 1 reference gaps R-010 and R-013 remain open and must not be hidden by the refactor. +Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignment-refactor`, continuing with P1A-07 MITC4 material/stiffness extraction. P1A-00 completed the architecture drift audit, P1A-01 created the module scaffold, P1A-02 extracted Core/Util plus Phase 1 Boundary/Load/Property model ownership, P1A-03 extracted Math primitives plus the solver adapter boundary, P1A-04 extracted the Abaqus Phase 1 parser into IO, P1A-05 extracted Results/reference comparison, and P1A-06 extracted MITC4 geometry/strain helpers into Element without changing solver behavior. This phase must align the current monolithic `include/fesa/fesa.hpp` implementation with the module ownership model in `docs/ARCHITECTURE.md` without changing solver behavior. Product-level Phase 1 reference gaps R-010 and R-013 remain open and must not be hidden by the refactor. ## Required Reading For New Agents 1. `AGENTS.md` @@ -37,7 +37,7 @@ Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignmen ## Phase Files - Active phase directory: `phases/1-structure-alignment-refactor` - Execute with: `python scripts/execute.py 1-structure-alignment-refactor` -- Step numbering is zero-based. `step0.md` is complete and wrote `phases/1-structure-alignment-refactor/step0-architecture-map.md`; `step1.md` is complete and created module scaffold headers, source directories, CMake source discovery, and umbrella compatibility smoke coverage; `step2.md` is complete and extracted Core/Util domain, diagnostics, DofManager ownership, AnalysisModel/AnalysisState, and Phase 1 Boundary/Load/Property model ownership; `step3.md` is complete and extracted Math primitives, sparse pattern data, dense matrix support, and solver adapter boundary; `step4.md` is complete and extracted the Abaqus parser into IO; `step5.md` is complete and extracted Results/reference comparison code; `step6.md` extracts MITC4 geometry/strain helpers; `step7.md` extracts MITC4 material/stiffness helpers; `step8.md` extracts Assembly and Analysis workflow; `step9.md` is the independent architecture evaluator closeout. +- Step numbering is zero-based. `step0.md` is complete and wrote `phases/1-structure-alignment-refactor/step0-architecture-map.md`; `step1.md` is complete and created module scaffold headers, source directories, CMake source discovery, and umbrella compatibility smoke coverage; `step2.md` is complete and extracted Core/Util domain, diagnostics, DofManager ownership, AnalysisModel/AnalysisState, and Phase 1 Boundary/Load/Property model ownership; `step3.md` is complete and extracted Math primitives, sparse pattern data, dense matrix support, and solver adapter boundary; `step4.md` is complete and extracted the Abaqus parser into IO; `step5.md` is complete and extracted Results/reference comparison code; `step6.md` is complete and extracted MITC4 geometry/strain helpers; `step7.md` extracts MITC4 material/stiffness helpers; `step8.md` extracts Assembly and Analysis workflow; `step9.md` is the independent architecture evaluator closeout. - Completed phase directory: `phases/1-linear-static-mitc4-rebaseline` - Historical execution command: `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; `step10.md` is complete and revalidated MITC4 stiffness, internal force, six-DOF transform, and drilling stabilization; `step11.md` is complete and added MITC4 membrane, bending, shear, twist, drilling-sensitivity, and thin-cantilever locking-sensitivity tests; `step12.md` is complete and revalidated full-space assembly, reduced projection, deterministic sparse-pattern scaffold, solver adapter injection, and full-vector internal/reaction force state; `step13.md` is complete and revalidated active AnalysisModel construction plus input-to-AnalysisState-to-U/RF result workflow; `step14.md` is complete and added the first stored Abaqus displacement regression for `quad_02_phase1`; `step15.md` is complete and recorded the independent evaluator closeout in `phases/1-linear-static-mitc4-rebaseline/step15-evaluator-report.md`. @@ -63,7 +63,7 @@ This phase is an architecture-preserving refactor. It must not change Phase 1 so | P1A-03 | completed | generator | Extract Math and solver adapter boundaries. | P1A-02 | Linear solver interface remains adapter-ready; int64 paths unchanged | | P1A-04 | completed | generator | Extract Abaqus parser into IO. | P1A-02 | Parser subset and unsupported-feature diagnostics unchanged | | P1A-05 | completed | generator | Extract Results model, writer boundary, CSV loader, and reference comparator. | P1A-02, P1A-04 | `U`/`RF` schema and `quad_02_phase1` regression unchanged | -| P1A-06 | pending | generator | Extract MITC4 geometry, director, strain, and tying helpers into Element. | P1A-03 | Geometry/strain tests and formulation signs unchanged | +| P1A-06 | completed | generator | Extract MITC4 geometry, director, strain, and tying helpers into Element. | P1A-03 | Geometry/strain tests and formulation signs unchanged | | P1A-07 | pending | generator | Extract MITC4 material, integration, stiffness, drilling, and internal-force helpers. | P1A-06 | Patch, drilling, stiffness, and locking-sensitivity tests unchanged | | P1A-08 | pending | generator | Extract Assembly and Analysis workflow. | P1A-02, P1A-03, P1A-05, P1A-07 | Full-vector RF, solver injection, and end-to-end reference regression unchanged | | P1A-09 | pending | evaluator | Independently evaluate final architecture alignment. | P1A-08 | `src/` ownership matches `ARCHITECTURE.md`; umbrella header is facade only | diff --git a/PROGRESS.md b/PROGRESS.md index 2b2a23f..f8a8de2 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,44 @@ 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 completed rebaseline execution path in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 15 are complete, and P1R-15 recorded a pass-with-documented-gaps evaluator closeout. The follow-up architecture refactor phase in `phases/1-structure-alignment-refactor` is underway because the current production implementation is concentrated in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00 through P1A-05 are complete, so the next step is P1A-06 MITC4 geometry/strain extraction. `quad_02_phase1.inp` is 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, MITC4 stiffness/drilling/internal-force scaffolding, MITC4 patch/locking-sensitivity tests, full-space assembly, reduced projection, sparse-pattern scaffold, solver adapter injection, full-vector internal/reaction force state, active AnalysisModel construction, input-to-AnalysisState-to-U/RF result workflow, and the first stored Abaqus displacement regression have been revalidated. Full PRD Phase 1 completion still depends on the open architecture/reference gaps R-014, R-010, and R-013. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. +Phase 1 has a completed rebaseline execution path in `phases/1-linear-static-mitc4-rebaseline`. Steps 0 through 15 are complete, and P1R-15 recorded a pass-with-documented-gaps evaluator closeout. The follow-up architecture refactor phase in `phases/1-structure-alignment-refactor` is underway because the current production implementation is concentrated in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00 through P1A-06 are complete, so the next step is P1A-07 MITC4 material/stiffness extraction. `quad_02_phase1.inp` is 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, MITC4 stiffness/drilling/internal-force scaffolding, MITC4 patch/locking-sensitivity tests, full-space assembly, reduced projection, sparse-pattern scaffold, solver adapter injection, full-vector internal/reaction force state, active AnalysisModel construction, input-to-AnalysisState-to-U/RF result workflow, and the first stored Abaqus displacement regression have been revalidated. Full PRD Phase 1 completion still depends on the open architecture/reference gaps R-014, R-010, and R-013. The old `phases/1-linear-static-mitc4` path is historical and superseded after the MITC4 formulation reset. ## Completed Work +### 2026-05-05 - P1A-06 MITC4 geometry strain extraction completed +Author: Codex + +Changed files: +- `CMakeLists.txt` +- `include/fesa/Element/Element.hpp` +- `include/fesa/Element/MITC4Geometry.hpp` +- `include/fesa/Element/MITC4Kinematics.hpp` +- `include/fesa/fesa.hpp` +- `tests/test_element_module_includes.cpp` +- `phases/1-structure-alignment-refactor/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Extracted `ShapeData`, `shapeFunctions`, `LocalBasis`, MITC4 natural-coordinate and tying-point helpers, director frames, geometry construction, integration basis construction, MITC4 diagnostics, and `computeLocalBasis` into `include/fesa/Element/MITC4Geometry.hpp`. +- Extracted MITC4 element DOF/strain vector aliases, strain component ordering helpers, local rotation mapping, director increment, displacement derivatives, direct covariant strain evaluation/rows, row evaluation, and MITC transverse shear tying rows into `include/fesa/Element/MITC4Kinematics.hpp`. +- Updated `include/fesa/Element/Element.hpp` to expose the MITC4 geometry/kinematics layer through the Element module. +- Updated `include/fesa/fesa.hpp` to include the Element facade and removed the relocated MITC4 geometry/strain definitions from the umbrella body. +- Added `fesa_element_module_tests`, a direct Element include smoke/regression test that does not include `fesa/fesa.hpp`. +- Preserved FESA/Abaqus S4 node order, A/B/C/D tying labels and signs, director fallback policy, invalid thickness/singular geometry diagnostics, direct covariant strain rows, and MITC shear interpolation behavior. +- No material law, stiffness integration, drilling stabilization, assembly, parser, results, or analysis behavior was changed. +- Remaining large groups in `fesa.hpp` are MITC4 material transform/integration/stiffness/internal-force helpers, Assembly helpers (`buildReducedSparsePattern`, `recoverFullReaction`, assembly/project functions), and Analysis workflow. + +Verification: +- First ran `python scripts\validate_workspace.py` after adding the direct Element include test; it failed as expected because `fesa/Element/Element.hpp` did not yet expose MITC4 geometry/strain symbols. +- After extraction, `python scripts\validate_workspace.py` configured CMake, built `fesa_core`, `fesa_tests`, `fesa_core_module_tests`, `fesa_math_module_tests`, `fesa_io_module_tests`, `fesa_results_module_tests`, and `fesa_element_module_tests`, and ran CTest successfully. +- CTest result: 6 test executables passed. + +Follow-up: +- Continue with P1A-07 MITC4 material/stiffness extraction. +- Keep R-014 open until P1A-09 independently accepts the final architecture alignment. +- Keep R-010 and R-013 open; this refactor does not add Abaqus RF artifacts or additional stored reference cases. + ### 2026-05-05 - P1A-05 Results reference extraction completed Author: Codex diff --git a/include/fesa/Element/Element.hpp b/include/fesa/Element/Element.hpp index 26adae0..b1b4a4f 100644 --- a/include/fesa/Element/Element.hpp +++ b/include/fesa/Element/Element.hpp @@ -1,5 +1,7 @@ #pragma once +#include "fesa/Element/MITC4Geometry.hpp" +#include "fesa/Element/MITC4Kinematics.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/Element/MITC4Geometry.hpp b/include/fesa/Element/MITC4Geometry.hpp new file mode 100644 index 0000000..78ddf8e --- /dev/null +++ b/include/fesa/Element/MITC4Geometry.hpp @@ -0,0 +1,265 @@ +#pragma once + +#include "fesa/Math/Math.hpp" +#include "fesa/Util/Util.hpp" + +#include +#include +#include +#include +#include + +namespace fesa { + +struct ShapeData { + std::array n{}; + std::array dr{}; + std::array ds{}; +}; + +inline ShapeData shapeFunctions(Real r, Real s) { + return { { + 0.25 * (1.0 - r) * (1.0 - s), + 0.25 * (1.0 + r) * (1.0 - s), + 0.25 * (1.0 + r) * (1.0 + s), + 0.25 * (1.0 - r) * (1.0 + s), + }, + { + -0.25 * (1.0 - s), + 0.25 * (1.0 - s), + 0.25 * (1.0 + s), + -0.25 * (1.0 + s), + }, + { + -0.25 * (1.0 - r), + -0.25 * (1.0 + r), + 0.25 * (1.0 + r), + 0.25 * (1.0 - r), + } }; +} + +struct LocalBasis { + Vec3 e1; + Vec3 e2; + 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 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); + 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) { + 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}; +} + +} // namespace fesa diff --git a/include/fesa/Element/MITC4Kinematics.hpp b/include/fesa/Element/MITC4Kinematics.hpp new file mode 100644 index 0000000..e6d57a9 --- /dev/null +++ b/include/fesa/Element/MITC4Kinematics.hpp @@ -0,0 +1,205 @@ +#pragma once + +#include "fesa/Element/MITC4Geometry.hpp" + +#include +#include +#include +#include + +namespace fesa { + +using MITC4ElementDofVector = std::array; +using MITC4StrainVector = std::array; +using MITC4StrainRow = std::array; +using MITC4MaterialMatrix = std::array, 6>; + +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 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; +} + +} // namespace fesa diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index f3b654b..d57d154 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -2,6 +2,7 @@ #include "fesa/Boundary/Boundary.hpp" #include "fesa/Core/Core.hpp" +#include "fesa/Element/Element.hpp" #include "fesa/IO/IO.hpp" #include "fesa/Load/Load.hpp" #include "fesa/Math/Math.hpp" @@ -68,161 +69,6 @@ inline std::vector recoverFullReaction(const DenseMatrix& k_full, const st return reaction; } -struct ShapeData { - std::array n{}; - std::array dr{}; - std::array ds{}; -}; - -inline ShapeData shapeFunctions(Real r, Real s) { - return {{ - 0.25 * (1.0 - r) * (1.0 - s), - 0.25 * (1.0 + r) * (1.0 - s), - 0.25 * (1.0 + r) * (1.0 + s), - 0.25 * (1.0 - r) * (1.0 + s), - }, - { - -0.25 * (1.0 - s), - 0.25 * (1.0 - s), - 0.25 * (1.0 + s), - -0.25 * (1.0 + s), - }, - { - -0.25 * (1.0 - r), - -0.25 * (1.0 + r), - 0.25 * (1.0 + r), - 0.25 * (1.0 - r), - }}; -} - -struct LocalBasis { - Vec3 e1; - Vec3 e2; - 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); - } -}; - -using MITC4ElementDofVector = std::array; -using MITC4StrainVector = std::array; -using MITC4StrainRow = std::array; -using MITC4MaterialMatrix = std::array, 6>; - -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); - } -}; - struct MITC4MaterialMatrixEvaluation { MITC4MaterialMatrix matrix{}; std::vector diagnostics; @@ -271,271 +117,6 @@ struct MITC4MaterialIntegrationData { } }; -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 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); - 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 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 std::array mitc4GaussQuadrature2x2x2() { const Real gauss = 1.0 / std::sqrt(3.0); const std::array points = {-gauss, gauss}; @@ -739,15 +320,6 @@ inline MITC4MaterialIntegrationData mitc4BuildMaterialIntegrationData(const MITC return data; } -inline LocalBasis computeLocalBasis(const std::array& coordinates) { - 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 ElementStiffnessOptions { Real drilling_stiffness_scale = 1.0e-3; }; diff --git a/phases/1-structure-alignment-refactor/index.json b/phases/1-structure-alignment-refactor/index.json index 94fe961..4afd2d2 100644 --- a/phases/1-structure-alignment-refactor/index.json +++ b/phases/1-structure-alignment-refactor/index.json @@ -8,7 +8,7 @@ { "step": 3, "name": "math-solver-extraction", "status": "completed" }, { "step": 4, "name": "io-parser-extraction", "status": "completed" }, { "step": 5, "name": "results-reference-extraction", "status": "completed" }, - { "step": 6, "name": "mitc4-geometry-strain-extraction", "status": "pending" }, + { "step": 6, "name": "mitc4-geometry-strain-extraction", "status": "completed" }, { "step": 7, "name": "mitc4-material-stiffness-extraction", "status": "pending" }, { "step": 8, "name": "assembly-analysis-extraction", "status": "pending" }, { "step": 9, "name": "architecture-evaluator-closeout", "status": "pending" } diff --git a/tests/test_element_module_includes.cpp b/tests/test_element_module_includes.cpp new file mode 100644 index 0000000..92c318d --- /dev/null +++ b/tests/test_element_module_includes.cpp @@ -0,0 +1,173 @@ +#include "fesa/Element/Element.hpp" + +#include +#include +#include + +namespace { + +void check(bool value, const char* message) { + if (!value) { + throw std::runtime_error(message); + } +} + +void checkNear(fesa::Real actual, fesa::Real expected, fesa::Real tolerance, const char* message) { + if (std::fabs(actual - expected) > tolerance) { + throw std::runtime_error(message); + } +} + +void checkVecNear(const fesa::Vec3& actual, const fesa::Vec3& expected, fesa::Real tolerance, const char* message) { + checkNear(actual.x, expected.x, tolerance, message); + checkNear(actual.y, expected.y, tolerance, message); + checkNear(actual.z, expected.z, tolerance, message); +} + +void checkRightHandedOrthonormal(const fesa::Vec3& e1, const fesa::Vec3& e2, const fesa::Vec3& e3) { + checkNear(fesa::norm(e1), 1.0, 1.0e-14, "e1 norm changed"); + checkNear(fesa::norm(e2), 1.0, 1.0e-14, "e2 norm changed"); + checkNear(fesa::norm(e3), 1.0, 1.0e-14, "e3 norm changed"); + checkNear(fesa::dot(e1, e2), 0.0, 1.0e-14, "e1/e2 orthogonality changed"); + checkNear(fesa::dot(e1, e3), 0.0, 1.0e-14, "e1/e3 orthogonality changed"); + checkNear(fesa::dot(e2, e3), 0.0, 1.0e-14, "e2/e3 orthogonality changed"); + checkNear(fesa::dot(fesa::cross(e1, e2), e3), 1.0, 1.0e-14, "basis handedness changed"); +} + +fesa::MITC4ElementDofVector zeroElementDofs() { + fesa::MITC4ElementDofVector values{}; + values.fill(0.0); + return values; +} + +} // namespace + +int main() { + const auto center = fesa::shapeFunctions(0.0, 0.0); + checkNear(center.n[0] + center.n[1] + center.n[2] + center.n[3], 1.0, 1.0e-15, "shape sum changed"); + checkNear(center.dr[0], -0.25, 1.0e-15, "center dN1/dxi changed"); + checkNear(center.ds[2], 0.25, 1.0e-15, "center dN3/deta changed"); + + const auto nodes = fesa::mitc4NodeNaturalCoordinates(); + checkNear(nodes[0].xi, -1.0, 1.0e-15, "node 1 xi changed"); + checkNear(nodes[0].eta, -1.0, 1.0e-15, "node 1 eta changed"); + checkNear(nodes[2].xi, 1.0, 1.0e-15, "node 3 xi changed"); + checkNear(nodes[2].eta, 1.0, 1.0e-15, "node 3 eta changed"); + for (std::size_t i = 0; i < 4; ++i) { + const auto corner = fesa::shapeFunctions(nodes[i].xi, nodes[i].eta); + checkNear(corner.n[i], 1.0, 1.0e-15, "corner interpolation changed"); + } + + const auto tying = fesa::mitc4TyingPoints(); + check(tying[0].label == "A" && tying[0].natural.xi == 0.0 && tying[0].natural.eta == -1.0, + "tying point A changed"); + check((tying[0].edge_node_indices == std::array{0, 1}), "tying edge A changed"); + check(tying[1].label == "B" && tying[1].natural.xi == -1.0 && tying[1].natural.eta == 0.0, + "tying point B changed"); + check((tying[1].edge_node_indices == std::array{0, 3}), "tying edge B changed"); + check(tying[2].label == "C" && tying[2].natural.xi == 0.0 && tying[2].natural.eta == 1.0, + "tying point C changed"); + check((tying[2].edge_node_indices == std::array{3, 2}), "tying edge C changed"); + check(tying[3].label == "D" && tying[3].natural.xi == 1.0 && tying[3].natural.eta == 0.0, + "tying point D changed"); + check((tying[3].edge_node_indices == std::array{1, 2}), "tying edge D changed"); + + const std::array coords = {{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + const auto geometry = fesa::buildMITC4Geometry(coords, 0.2); + check(geometry.ok(), "flat geometry should remain valid"); + checkVecNear(geometry.center_normal, {0, 0, 1}, 1.0e-14, "center normal changed"); + for (const auto& frame : geometry.nodal_frames) { + checkVecNear(frame.v1, {1, 0, 0}, 1.0e-14, "director v1 changed"); + checkVecNear(frame.v2, {0, 1, 0}, 1.0e-14, "director v2 changed"); + checkVecNear(frame.vn, {0, 0, 1}, 1.0e-14, "director normal changed"); + checkRightHandedOrthonormal(frame.v1, frame.v2, frame.vn); + } + + const auto basis = fesa::computeMITC4IntegrationBasis(geometry, 0.0, 0.0, 0.0); + check(basis.ok(), "flat integration basis should remain valid"); + checkVecNear(basis.g1, {0.5, 0, 0}, 1.0e-14, "g1 changed"); + checkVecNear(basis.g2, {0, 0.5, 0}, 1.0e-14, "g2 changed"); + checkVecNear(basis.g3, {0, 0, 0.1}, 1.0e-14, "g3 changed"); + checkRightHandedOrthonormal(basis.local.e1, basis.local.e2, basis.local.e3); + checkNear(basis.jacobian, 0.025, 1.0e-14, "Jacobian changed"); + + const auto invalid_thickness = fesa::buildMITC4Geometry(coords, 0.0); + check(!invalid_thickness.ok(), "invalid thickness should remain diagnosed"); + check(fesa::containsDiagnostic(invalid_thickness.diagnostics, "FESA-MITC4-THICKNESS"), "thickness diagnostic changed"); + const std::array collinear = {{{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}}}; + const auto singular = fesa::buildMITC4Geometry(collinear, 0.1); + check(!singular.ok(), "singular normal should remain diagnosed"); + check(fesa::containsDiagnostic(singular.diagnostics, "FESA-MITC4-SINGULAR-NORMAL"), "normal diagnostic changed"); + + const auto rotations = fesa::mitc4LocalRotations(geometry.nodal_frames[0], {1.0, 2.0, 3.0}); + checkNear(rotations.alpha, 1.0, 1.0e-14, "alpha rotation changed"); + checkNear(rotations.beta, 2.0, 1.0e-14, "beta rotation changed"); + checkNear(rotations.gamma, 3.0, 1.0e-14, "gamma rotation changed"); + checkVecNear(fesa::mitc4DirectorIncrement(geometry.nodal_frames[0], {0.0, 2.0, 5.0}), {2, 0, 0}, 1.0e-14, + "director increment changed"); + + auto 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); + check(mid.ok(), "midsurface displacement interpolation changed"); + checkVecNear(mid.displacement, {1.0, 0.5, 0.0}, 1.0e-14, "midsurface displacement value changed"); + const auto top = fesa::mitc4DisplacementDerivatives(geometry, dofs, 0.0, 0.0, 1.0); + check(top.ok(), "top displacement interpolation changed"); + checkVecNear(top.displacement, {1.2, 0.5, 0.0}, 1.0e-14, "top displacement value changed"); + + 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); + check(rows.ok(), "direct covariant strain rows should remain valid"); + check(fesa::mitc4StrainComponentLabels()[0] == "eps11", "strain label eps11 changed"); + check(fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23) == 3, "gamma23 index changed"); + check(fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13) == 4, "gamma13 index changed"); + + const fesa::Real h = 1.0e-6; + const std::array checked_dofs = {0, 7, 10, 11}; + 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); + check(plus_strain.ok() && minus_strain.ok(), "direct strain perturbation changed"); + 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); + checkNear(rows.rows[component][dof], finite_difference, 1.0e-8, "direct strain row finite difference changed"); + } + } + + const auto gamma23 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma23); + const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13); + const auto direct_at_point = fesa::mitc4DirectCovariantStrainRows(geometry, 0.5, 0.25, 0.0); + 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 tied = fesa::mitc4TiedCovariantStrainRows(geometry, 0.5, 0.25, 0.0); + check(direct_at_point.ok() && tied.ok(), "MITC tied rows should remain valid"); + + const std::size_t node2_ry = 6 + 4; + const fesa::Real expected_gamma13 = + 0.5 * (1.0 - 0.25) * direct_a.rows[gamma13][node2_ry] + 0.5 * (1.0 + 0.25) * direct_c.rows[gamma13][node2_ry]; + checkNear(tied.rows[gamma13][node2_ry], expected_gamma13, 1.0e-14, "MITC gamma13 tying interpolation changed"); + check(std::fabs(tied.rows[gamma13][node2_ry] - direct_at_point.rows[gamma13][node2_ry]) > 1.0e-4, + "MITC gamma13 should not be direct Gauss shear"); + + const std::size_t node4_rx = 3 * 6 + 3; + const fesa::Real expected_gamma23 = + 0.5 * (1.0 - 0.5) * direct_b.rows[gamma23][node4_rx] + 0.5 * (1.0 + 0.5) * direct_d.rows[gamma23][node4_rx]; + checkNear(tied.rows[gamma23][node4_rx], expected_gamma23, 1.0e-14, "MITC gamma23 tying interpolation changed"); + check(std::fabs(tied.rows[gamma23][node4_rx] - direct_at_point.rows[gamma23][node4_rx]) > 1.0e-4, + "MITC gamma23 should not be direct Gauss shear"); + + return 0; +}