From 918e219c48c7d5e8642d9e2ddd823ddac82797b7 Mon Sep 17 00:00:00 2001 From: NINI Date: Tue, 5 May 2026 23:26:11 +0900 Subject: [PATCH] refactor: extract mitc4 material stiffness helpers --- CMakeLists.txt | 6 + PLAN.md | 6 +- PROGRESS.md | 41 +- include/fesa/Element/Element.hpp | 2 + include/fesa/Element/MITC4Kinematics.hpp | 21 +- .../fesa/Element/MITC4MaterialIntegration.hpp | 182 +++++++ include/fesa/Element/MITC4Stiffness.hpp | 228 ++++++++ .../Material/MITC4PlaneStressMaterial.hpp | 121 +++++ include/fesa/Material/Material.hpp | 1 + include/fesa/fesa.hpp | 486 +----------------- .../1-structure-alignment-refactor/index.json | 2 +- .../test_mitc4_stiffness_module_includes.cpp | 117 +++++ 12 files changed, 715 insertions(+), 498 deletions(-) create mode 100644 include/fesa/Element/MITC4MaterialIntegration.hpp create mode 100644 include/fesa/Element/MITC4Stiffness.hpp create mode 100644 include/fesa/Material/MITC4PlaneStressMaterial.hpp create mode 100644 tests/test_mitc4_stiffness_module_includes.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 105f2f8..d9b5697 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,6 +35,9 @@ target_compile_definitions(fesa_results_module_tests PRIVATE FESA_SOURCE_DIR="${ add_executable(fesa_element_module_tests tests/test_element_module_includes.cpp) target_link_libraries(fesa_element_module_tests PRIVATE fesa_core) +add_executable(fesa_mitc4_stiffness_module_tests tests/test_mitc4_stiffness_module_includes.cpp) +target_link_libraries(fesa_mitc4_stiffness_module_tests PRIVATE fesa_core) + if(MSVC) target_compile_options(fesa_core PRIVATE /W4 /permissive-) target_compile_options(fesa_tests PRIVATE /W4 /permissive-) @@ -43,6 +46,7 @@ if(MSVC) 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-) + target_compile_options(fesa_mitc4_stiffness_module_tests PRIVATE /W4 /permissive-) else() target_compile_options(fesa_core PRIVATE -Wall -Wextra -Wpedantic) target_compile_options(fesa_tests PRIVATE -Wall -Wextra -Wpedantic) @@ -51,6 +55,7 @@ else() 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) + target_compile_options(fesa_mitc4_stiffness_module_tests PRIVATE -Wall -Wextra -Wpedantic) endif() add_test(NAME fesa_tests COMMAND fesa_tests) @@ -59,3 +64,4 @@ 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) +add_test(NAME fesa_mitc4_stiffness_module_tests COMMAND fesa_mitc4_stiffness_module_tests) diff --git a/PLAN.md b/PLAN.md index 7116702..6a08fd1 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-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. +Execute the Phase 1 structure-alignment refactor in `phases/1-structure-alignment-refactor`, continuing with P1A-08 Assembly and Analysis workflow 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, P1A-06 extracted MITC4 geometry/strain helpers into Element, and P1A-07 extracted MITC4 material/stiffness helpers into Material and 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` 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. +- 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` is complete and extracted 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`. @@ -64,7 +64,7 @@ This phase is an architecture-preserving refactor. It must not change Phase 1 so | 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 | 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-07 | completed | 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 f8a8de2..1966034 100644 --- a/PROGRESS.md +++ b/PROGRESS.md @@ -13,10 +13,47 @@ 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-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. +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 remaining Assembly and Analysis workflow code still lives in `include/fesa/fesa.hpp` instead of the module directories documented in `docs/ARCHITECTURE.md`; P1A-00 through P1A-07 are complete, so the next step is P1A-08 Assembly and Analysis workflow 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-07 MITC4 material stiffness extraction completed +Author: Codex + +Changed files: +- `CMakeLists.txt` +- `include/fesa/Element/Element.hpp` +- `include/fesa/Element/MITC4Kinematics.hpp` +- `include/fesa/Element/MITC4MaterialIntegration.hpp` +- `include/fesa/Element/MITC4Stiffness.hpp` +- `include/fesa/Material/MITC4PlaneStressMaterial.hpp` +- `include/fesa/Material/Material.hpp` +- `include/fesa/fesa.hpp` +- `tests/test_mitc4_stiffness_module_includes.cpp` +- `phases/1-structure-alignment-refactor/index.json` +- `PLAN.md` +- `PROGRESS.md` + +Summary: +- Extracted MITC4 strain component ordering, strain vector alias, material matrix alias, plane-stress material matrix diagnostics, material-vector multiply, material-vector dot product, and material matrix transform into `include/fesa/Material/MITC4PlaneStressMaterial.hpp`. +- Extracted `2 x 2 x 2` Gauss integration points, covariant-to-local strain transform, tensor/vector conversion helpers, and MITC4 material integration sample/data construction into `include/fesa/Element/MITC4MaterialIntegration.hpp`. +- Extracted MITC4 stiffness options, local/global DOF transform, strain-row local transform, `B^T D B` accumulation, drilling stabilization, local-to-global stiffness transform, element stiffness, internal force, and `MITC4ElementKernel` into `include/fesa/Element/MITC4Stiffness.hpp`. +- Updated Element and Material facade headers so direct module includes expose the relocated MITC4 material/stiffness surface without including `fesa/fesa.hpp`. +- Reduced `include/fesa/fesa.hpp` to keep only the remaining Assembly/Analysis workflow and umbrella includes; MITC4 material/stiffness implementation bodies no longer live there. +- Preserved `drilling_stiffness_scale = 1.0e-3`, the minimum-positive-physical-local-diagonal drilling policy, `2 x 2 x 2` integration, stiffness symmetry, internal-force `K_e u_e`, and all existing MITC4 patch/locking characterization behavior. +- No parser, result schema, reference tolerance, S4R, reduced integration, hourglass, nonlinear, pressure-load, HDF5, MKL, or TBB behavior was added. +- Remaining large groups in `fesa.hpp` are Assembly helpers (`buildReducedSparsePattern`, `recoverFullReaction`, assembly/project functions) and Analysis workflow. + +Verification: +- First ran `python scripts\validate_workspace.py` after adding `fesa_mitc4_stiffness_module_tests`; it failed as expected because `fesa/Element/MITC4Stiffness.hpp` did not exist yet. +- 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`, `fesa_element_module_tests`, and `fesa_mitc4_stiffness_module_tests`, and ran CTest successfully. +- CTest result: 7 test executables passed. + +Follow-up: +- Continue with P1A-08 Assembly and Analysis workflow 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-06 MITC4 geometry strain extraction completed Author: Codex @@ -1132,7 +1169,7 @@ Verification: - `python scripts/validate_workspace.py` ran, but reported no configured validation commands. ## Known Blockers -- Phase 1 architecture is not yet accepted: production code is concentrated in `include/fesa/fesa.hpp` instead of the `src/` module layout documented in `docs/ARCHITECTURE.md`. +- Phase 1 architecture is not yet accepted: remaining Assembly and Analysis workflow code still lives in `include/fesa/fesa.hpp` instead of the module layout documented in `docs/ARCHITECTURE.md`. - No reaction-force reference artifact exists yet under `references/`. - The PRD target of three stored Phase 1 reference cases is not yet satisfied; only `quad_02_phase1` is an active stored displacement regression. - The current initial `quad_01.inp` reference contains `S4R`, `Part/Assembly/Instance`, `*Density`, and `NLGEOM=YES`, so it is not a Phase 1 parser acceptance case as-is. diff --git a/include/fesa/Element/Element.hpp b/include/fesa/Element/Element.hpp index b1b4a4f..43e4282 100644 --- a/include/fesa/Element/Element.hpp +++ b/include/fesa/Element/Element.hpp @@ -2,6 +2,8 @@ #include "fesa/Element/MITC4Geometry.hpp" #include "fesa/Element/MITC4Kinematics.hpp" +#include "fesa/Element/MITC4MaterialIntegration.hpp" +#include "fesa/Element/MITC4Stiffness.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/Element/MITC4Kinematics.hpp b/include/fesa/Element/MITC4Kinematics.hpp index e6d57a9..8862020 100644 --- a/include/fesa/Element/MITC4Kinematics.hpp +++ b/include/fesa/Element/MITC4Kinematics.hpp @@ -1,35 +1,16 @@ #pragma once #include "fesa/Element/MITC4Geometry.hpp" +#include "fesa/Material/MITC4PlaneStressMaterial.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; diff --git a/include/fesa/Element/MITC4MaterialIntegration.hpp b/include/fesa/Element/MITC4MaterialIntegration.hpp new file mode 100644 index 0000000..7ce33eb --- /dev/null +++ b/include/fesa/Element/MITC4MaterialIntegration.hpp @@ -0,0 +1,182 @@ +#pragma once + +#include "fesa/Element/MITC4Kinematics.hpp" +#include "fesa/Material/MITC4PlaneStressMaterial.hpp" + +#include +#include +#include + +namespace fesa { + +struct MITC4StrainTransform { + MITC4MaterialMatrix matrix{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4IntegrationPoint { + Real xi = 0.0; + Real eta = 0.0; + Real zeta = 0.0; + Real weight = 0.0; +}; + +struct MITC4MaterialIntegrationSample { + MITC4IntegrationPoint point; + MITC4IntegrationBasis basis; + MITC4StrainRows strain_rows; + MITC4MaterialMatrix local_material{}; + MITC4MaterialMatrix strain_transform{}; + MITC4MaterialMatrix convected_material{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4MaterialIntegrationData { + std::vector samples; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +inline std::array mitc4GaussQuadrature2x2x2() { + const Real gauss = 1.0 / std::sqrt(3.0); + const std::array points = {-gauss, gauss}; + std::array integration_points{}; + std::size_t index = 0; + for (Real xi : points) { + for (Real eta : points) { + for (Real zeta : points) { + integration_points[index++] = {xi, eta, zeta, 1.0}; + } + } + } + return integration_points; +} + +inline std::array, 3> mitc4TensorFromEngineeringComponent(std::size_t component) { + std::array, 3> tensor{}; + switch (static_cast(component)) { + case MITC4StrainComponent::Eps11: + tensor[0][0] = 1.0; + break; + case MITC4StrainComponent::Eps22: + tensor[1][1] = 1.0; + break; + case MITC4StrainComponent::Eps33: + tensor[2][2] = 1.0; + break; + case MITC4StrainComponent::Gamma23: + tensor[1][2] = 0.5; + tensor[2][1] = 0.5; + break; + case MITC4StrainComponent::Gamma13: + tensor[0][2] = 0.5; + tensor[2][0] = 0.5; + break; + case MITC4StrainComponent::Gamma12: + tensor[0][1] = 0.5; + tensor[1][0] = 0.5; + break; + } + return tensor; +} + +inline MITC4StrainVector mitc4EngineeringVectorFromTensor(const std::array, 3>& tensor) { + MITC4StrainVector vector{}; + vector[strainComponentIndex(MITC4StrainComponent::Eps11)] = tensor[0][0]; + vector[strainComponentIndex(MITC4StrainComponent::Eps22)] = tensor[1][1]; + vector[strainComponentIndex(MITC4StrainComponent::Eps33)] = tensor[2][2]; + vector[strainComponentIndex(MITC4StrainComponent::Gamma23)] = 2.0 * tensor[1][2]; + vector[strainComponentIndex(MITC4StrainComponent::Gamma13)] = 2.0 * tensor[0][2]; + vector[strainComponentIndex(MITC4StrainComponent::Gamma12)] = 2.0 * tensor[0][1]; + return vector; +} + +inline MITC4StrainTransform mitc4CovariantToLocalStrainTransform(const MITC4IntegrationBasis& basis, + Real tolerance = 1.0e-12) { + MITC4StrainTransform result; + result.diagnostics = basis.diagnostics; + const Real jacobian = dot(cross(basis.g1, basis.g2), basis.g3); + if (!isFinite(jacobian) || std::fabs(jacobian) <= tolerance) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 material transform Jacobian is near zero")); + return result; + } + if (hasError(result.diagnostics)) { + return result; + } + + const std::array contravariant = { + (1.0 / jacobian) * cross(basis.g2, basis.g3), + (1.0 / jacobian) * cross(basis.g3, basis.g1), + (1.0 / jacobian) * cross(basis.g1, basis.g2)}; + const std::array local = {basis.local.e1, basis.local.e2, basis.local.e3}; + std::array, 3> direction_cosines{}; + for (std::size_t local_axis = 0; local_axis < 3; ++local_axis) { + for (std::size_t convected_axis = 0; convected_axis < 3; ++convected_axis) { + direction_cosines[local_axis][convected_axis] = dot(local[local_axis], contravariant[convected_axis]); + } + } + + for (std::size_t column = 0; column < 6; ++column) { + const auto covariant_tensor = mitc4TensorFromEngineeringComponent(column); + std::array, 3> local_tensor{}; + for (std::size_t a = 0; a < 3; ++a) { + for (std::size_t b = 0; b < 3; ++b) { + for (std::size_t i = 0; i < 3; ++i) { + for (std::size_t j = 0; j < 3; ++j) { + local_tensor[a][b] += direction_cosines[a][i] * direction_cosines[b][j] * covariant_tensor[i][j]; + } + } + } + } + const auto local_vector = mitc4EngineeringVectorFromTensor(local_tensor); + for (std::size_t row = 0; row < 6; ++row) { + result.matrix[row][column] = local_vector[row]; + } + } + return result; +} + +inline MITC4MaterialIntegrationData mitc4BuildMaterialIntegrationData(const MITC4Geometry& geometry, + Real elastic_modulus, + Real poisson_ratio, + Real shear_correction = 5.0 / 6.0) { + MITC4MaterialIntegrationData data; + const auto material = mitc4PlaneStressMaterialMatrix(elastic_modulus, poisson_ratio, shear_correction); + appendDiagnostics(data.diagnostics, material.diagnostics); + if (hasError(data.diagnostics)) { + return data; + } + + for (const MITC4IntegrationPoint& point : mitc4GaussQuadrature2x2x2()) { + MITC4MaterialIntegrationSample sample; + sample.point = point; + sample.local_material = material.matrix; + sample.basis = computeMITC4IntegrationBasis(geometry, point.xi, point.eta, point.zeta); + appendDiagnostics(sample.diagnostics, sample.basis.diagnostics); + const auto transform = mitc4CovariantToLocalStrainTransform(sample.basis); + sample.strain_transform = transform.matrix; + appendDiagnostics(sample.diagnostics, transform.diagnostics); + sample.strain_rows = mitc4TiedCovariantStrainRows(geometry, point.xi, point.eta, point.zeta); + appendDiagnostics(sample.diagnostics, sample.strain_rows.diagnostics); + if (!hasError(sample.diagnostics)) { + sample.convected_material = mitc4TransformMaterialMatrix(sample.local_material, sample.strain_transform); + } + appendDiagnostics(data.diagnostics, sample.diagnostics); + data.samples.push_back(sample); + } + return data; +} + +} // namespace fesa diff --git a/include/fesa/Element/MITC4Stiffness.hpp b/include/fesa/Element/MITC4Stiffness.hpp new file mode 100644 index 0000000..0704464 --- /dev/null +++ b/include/fesa/Element/MITC4Stiffness.hpp @@ -0,0 +1,228 @@ +#pragma once + +#include "fesa/Element/MITC4MaterialIntegration.hpp" +#include "fesa/Math/Math.hpp" + +#include +#include +#include +#include +#include + +namespace fesa { + +struct ElementStiffnessOptions { + Real drilling_stiffness_scale = 1.0e-3; +}; + +struct MITC4DrillingStabilizationResult { + DenseMatrix local_with_drilling; + Real reference_diagonal = 0.0; + Real drilling_stiffness = 0.0; + Real drilling_stiffness_scale = 0.0; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +struct MITC4ElementStiffnessResult { + DenseMatrix local_without_drilling; + DenseMatrix local_with_drilling; + DenseMatrix global; + std::size_t integration_point_count = 0; + Real drilling_reference_diagonal = 0.0; + Real drilling_stiffness = 0.0; + Real drilling_stiffness_scale = 0.0; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +inline DenseMatrix mitc4LocalDofTransform(const MITC4Geometry& geometry) { + DenseMatrix transform(24, 24); + for (LocalIndex node = 0; node < 4; ++node) { + const LocalIndex base = 6 * node; + const auto& frame = geometry.nodal_frames[static_cast(node)]; + const std::array axes = {frame.v1, frame.v2, frame.vn}; + for (LocalIndex local_axis = 0; local_axis < 3; ++local_axis) { + const Vec3 axis = axes[static_cast(local_axis)]; + transform(base + local_axis, base + 0) = axis.x; + transform(base + local_axis, base + 1) = axis.y; + transform(base + local_axis, base + 2) = axis.z; + transform(base + 3 + local_axis, base + 3) = axis.x; + transform(base + 3 + local_axis, base + 4) = axis.y; + transform(base + 3 + local_axis, base + 5) = axis.z; + } + } + return transform; +} + +inline MITC4StrainRows mitc4TransformStrainRowsToLocalDofs(const MITC4StrainRows& global_rows, + const DenseMatrix& local_dof_transform) { + MITC4StrainRows local_rows; + local_rows.diagnostics = global_rows.diagnostics; + if (hasError(local_rows.diagnostics)) { + return local_rows; + } + for (std::size_t component = 0; component < 6; ++component) { + for (LocalIndex local_dof = 0; local_dof < 24; ++local_dof) { + Real value = 0.0; + for (LocalIndex global_dof = 0; global_dof < 24; ++global_dof) { + value += global_rows.rows[component][static_cast(global_dof)] * + local_dof_transform(local_dof, global_dof); + } + local_rows.rows[component][static_cast(local_dof)] = value; + } + } + return local_rows; +} + +inline void accumulateMITC4BtDB(DenseMatrix& stiffness, const MITC4StrainRows& rows, + const MITC4MaterialMatrix& material, Real factor) { + for (LocalIndex i = 0; i < 24; ++i) { + for (LocalIndex j = 0; j < 24; ++j) { + Real value = 0.0; + for (std::size_t a = 0; a < 6; ++a) { + for (std::size_t b = 0; b < 6; ++b) { + value += rows.rows[a][static_cast(i)] * material[a][b] * + rows.rows[b][static_cast(j)]; + } + } + stiffness.add(i, j, value * factor); + } + } +} + +inline Real mitc4MinimumPositivePhysicalLocalDiagonal(const DenseMatrix& local_without_drilling, + Real tolerance = 1.0e-12) { + Real minimum = std::numeric_limits::infinity(); + for (LocalIndex node = 0; node < 4; ++node) { + for (LocalIndex local_dof = 0; local_dof < 5; ++local_dof) { + const Real diagonal = local_without_drilling(6 * node + local_dof, 6 * node + local_dof); + if (isFinite(diagonal) && diagonal > tolerance) { + minimum = std::min(minimum, diagonal); + } + } + } + return minimum; +} + +inline MITC4DrillingStabilizationResult mitc4ApplyDrillingStabilization(const DenseMatrix& local_without_drilling, + Real drilling_stiffness_scale, + Real tolerance = 1.0e-12) { + MITC4DrillingStabilizationResult result; + result.local_with_drilling = local_without_drilling; + result.drilling_stiffness_scale = drilling_stiffness_scale; + if (!isFinite(drilling_stiffness_scale) || drilling_stiffness_scale < 0.0) { + result.diagnostics.push_back( + mitc4Diagnostic("FESA-MITC4-DRILLING-SCALE", "MITC4 drilling stiffness scale must be non-negative and finite")); + return result; + } + + result.reference_diagonal = mitc4MinimumPositivePhysicalLocalDiagonal(local_without_drilling, tolerance); + if (!isFinite(result.reference_diagonal)) { + result.diagnostics.push_back(mitc4Diagnostic("FESA-MITC4-DRILLING-REFERENCE", + "MITC4 drilling stiffness reference diagonal is not positive")); + return result; + } + + result.drilling_stiffness = drilling_stiffness_scale * result.reference_diagonal; + for (LocalIndex node = 0; node < 4; ++node) { + const LocalIndex gamma = 6 * node + 5; + result.local_with_drilling.add(gamma, gamma, result.drilling_stiffness); + } + return result; +} + +inline DenseMatrix mitc4TransformLocalStiffnessToGlobal(const DenseMatrix& local_stiffness, + const DenseMatrix& local_dof_transform) { + DenseMatrix global(24, 24); + for (LocalIndex i = 0; i < 24; ++i) { + for (LocalIndex j = 0; j < 24; ++j) { + Real value = 0.0; + for (LocalIndex a = 0; a < 24; ++a) { + for (LocalIndex b = 0; b < 24; ++b) { + value += local_dof_transform(a, i) * local_stiffness(a, b) * local_dof_transform(b, j); + } + } + global(i, j) = value; + } + } + return global; +} + +inline MITC4ElementStiffnessResult mitc4ElementStiffness(const std::array& coordinates, + Real elastic_modulus, Real poisson_ratio, Real thickness, + ElementStiffnessOptions options = {}) { + MITC4ElementStiffnessResult result; + result.local_without_drilling = DenseMatrix(24, 24); + result.local_with_drilling = DenseMatrix(24, 24); + result.global = DenseMatrix(24, 24); + result.drilling_stiffness_scale = options.drilling_stiffness_scale; + + const MITC4Geometry geometry = buildMITC4Geometry(coordinates, thickness); + appendDiagnostics(result.diagnostics, geometry.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + + const MITC4MaterialIntegrationData integration = + mitc4BuildMaterialIntegrationData(geometry, elastic_modulus, poisson_ratio); + appendDiagnostics(result.diagnostics, integration.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + result.integration_point_count = integration.samples.size(); + + const DenseMatrix local_dof_transform = mitc4LocalDofTransform(geometry); + for (const MITC4MaterialIntegrationSample& sample : integration.samples) { + const MITC4StrainRows local_rows = + mitc4TransformStrainRowsToLocalDofs(sample.strain_rows, local_dof_transform); + appendDiagnostics(result.diagnostics, local_rows.diagnostics); + if (hasError(result.diagnostics)) { + return result; + } + const Real factor = std::fabs(sample.basis.jacobian) * sample.point.weight; + accumulateMITC4BtDB(result.local_without_drilling, local_rows, sample.convected_material, factor); + } + + const auto drilling = mitc4ApplyDrillingStabilization(result.local_without_drilling, options.drilling_stiffness_scale); + appendDiagnostics(result.diagnostics, drilling.diagnostics); + result.local_with_drilling = drilling.local_with_drilling; + result.drilling_reference_diagonal = drilling.reference_diagonal; + result.drilling_stiffness = drilling.drilling_stiffness; + result.drilling_stiffness_scale = drilling.drilling_stiffness_scale; + if (hasError(result.diagnostics)) { + return result; + } + + result.global = mitc4TransformLocalStiffnessToGlobal(result.local_with_drilling, local_dof_transform); + return result; +} + +inline std::vector mitc4ElementInternalForce(const MITC4ElementStiffnessResult& stiffness, + const std::vector& element_displacement) { + if (stiffness.global.rows() != static_cast(element_displacement.size())) { + throw std::runtime_error("MITC4 internal force size mismatch"); + } + return stiffness.global.multiply(element_displacement); +} + +class MITC4ElementKernel { + public: + DenseMatrix stiffness(const std::array& coordinates, Real elastic_modulus, Real poisson_ratio, Real thickness, + ElementStiffnessOptions options = {}) const { + const auto result = mitc4ElementStiffness(coordinates, elastic_modulus, poisson_ratio, thickness, options); + if (!result.ok()) { + throw std::runtime_error(result.diagnostics.empty() ? "invalid MITC4 stiffness" + : result.diagnostics.front().message); + } + return result.global; + } +}; + +} // namespace fesa diff --git a/include/fesa/Material/MITC4PlaneStressMaterial.hpp b/include/fesa/Material/MITC4PlaneStressMaterial.hpp new file mode 100644 index 0000000..a0a7dce --- /dev/null +++ b/include/fesa/Material/MITC4PlaneStressMaterial.hpp @@ -0,0 +1,121 @@ +#pragma once + +#include "fesa/Math/Vector.hpp" +#include "fesa/Util/Diagnostics.hpp" + +#include +#include +#include +#include +#include + +namespace fesa { + +using MITC4StrainVector = 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 MITC4MaterialMatrixEvaluation { + MITC4MaterialMatrix matrix{}; + std::vector diagnostics; + + bool ok() const { + return !hasError(diagnostics); + } +}; + +inline Diagnostic mitc4MaterialDiagnostic(std::string code, std::string message) { + return makeDiagnostic(Severity::Error, std::move(code), std::move(message), "mitc4", "", 0); +} + +inline MITC4MaterialMatrixEvaluation mitc4PlaneStressMaterialMatrix(Real elastic_modulus, Real poisson_ratio, + Real shear_correction = 5.0 / 6.0, + Real tolerance = 1.0e-12) { + MITC4MaterialMatrixEvaluation result; + if (!isFinite(elastic_modulus) || elastic_modulus <= tolerance) { + result.diagnostics.push_back( + mitc4MaterialDiagnostic("FESA-MITC4-MATERIAL", "MITC4 elastic modulus must be positive and finite")); + } + if (!isFinite(poisson_ratio) || poisson_ratio <= -1.0 || poisson_ratio >= 0.5) { + result.diagnostics.push_back( + mitc4MaterialDiagnostic("FESA-MITC4-POISSON", "MITC4 isotropic Poisson ratio must satisfy -1 < nu < 0.5")); + } + if (!isFinite(shear_correction) || shear_correction <= tolerance) { + result.diagnostics.push_back(mitc4MaterialDiagnostic("FESA-MITC4-SHEAR-CORRECTION", + "MITC4 shear correction factor must be positive and finite")); + } + if (hasError(result.diagnostics)) { + return result; + } + + const Real scale = elastic_modulus / (1.0 - poisson_ratio * poisson_ratio); + const Real shear_modulus = elastic_modulus / (2.0 * (1.0 + poisson_ratio)); + const std::size_t eps11 = strainComponentIndex(MITC4StrainComponent::Eps11); + const std::size_t eps22 = strainComponentIndex(MITC4StrainComponent::Eps22); + const std::size_t gamma23 = strainComponentIndex(MITC4StrainComponent::Gamma23); + const std::size_t gamma13 = strainComponentIndex(MITC4StrainComponent::Gamma13); + const std::size_t gamma12 = strainComponentIndex(MITC4StrainComponent::Gamma12); + + result.matrix[eps11][eps11] = scale; + result.matrix[eps11][eps22] = poisson_ratio * scale; + result.matrix[eps22][eps11] = poisson_ratio * scale; + result.matrix[eps22][eps22] = scale; + result.matrix[gamma23][gamma23] = shear_correction * shear_modulus; + result.matrix[gamma13][gamma13] = shear_correction * shear_modulus; + result.matrix[gamma12][gamma12] = shear_modulus; + return result; +} + +inline MITC4StrainVector multiplyMITC4MaterialMatrix(const MITC4MaterialMatrix& matrix, + const MITC4StrainVector& vector) { + MITC4StrainVector result{}; + for (std::size_t row = 0; row < 6; ++row) { + for (std::size_t col = 0; col < 6; ++col) { + result[row] += matrix[row][col] * vector[col]; + } + } + return result; +} + +inline Real dotMITC4Vector(const MITC4StrainVector& a, const MITC4StrainVector& b) { + Real result = 0.0; + for (std::size_t i = 0; i < 6; ++i) { + result += a[i] * b[i]; + } + return result; +} + +inline MITC4MaterialMatrix mitc4TransformMaterialMatrix(const MITC4MaterialMatrix& local_material, + const MITC4MaterialMatrix& covariant_to_local) { + MITC4MaterialMatrix result{}; + for (std::size_t i = 0; i < 6; ++i) { + for (std::size_t j = 0; j < 6; ++j) { + Real value = 0.0; + for (std::size_t a = 0; a < 6; ++a) { + for (std::size_t b = 0; b < 6; ++b) { + value += covariant_to_local[a][i] * local_material[a][b] * covariant_to_local[b][j]; + } + } + result[i][j] = value; + } + } + return result; +} + +} // namespace fesa diff --git a/include/fesa/Material/Material.hpp b/include/fesa/Material/Material.hpp index 8359ce8..49713c3 100644 --- a/include/fesa/Material/Material.hpp +++ b/include/fesa/Material/Material.hpp @@ -1,5 +1,6 @@ #pragma once +#include "fesa/Material/MITC4PlaneStressMaterial.hpp" #include "fesa/ModuleInfo.hpp" namespace fesa::module { diff --git a/include/fesa/fesa.hpp b/include/fesa/fesa.hpp index d57d154..a0c0ad7 100644 --- a/include/fesa/fesa.hpp +++ b/include/fesa/fesa.hpp @@ -5,6 +5,7 @@ #include "fesa/Element/Element.hpp" #include "fesa/IO/IO.hpp" #include "fesa/Load/Load.hpp" +#include "fesa/Material/Material.hpp" #include "fesa/Math/Math.hpp" #include "fesa/ModuleInfo.hpp" #include "fesa/Property/Property.hpp" @@ -57,7 +58,9 @@ inline SparsePattern buildReducedSparsePattern(const Domain& domain, const DofMa return pattern; } -inline std::vector recoverFullReaction(const DenseMatrix& k_full, const std::vector& u_full, const std::vector& f_full) { +inline std::vector recoverFullReaction(const DenseMatrix& k_full, + const std::vector& u_full, + const std::vector& f_full) { if (k_full.rows() != k_full.cols() || static_cast(u_full.size()) != k_full.cols() || static_cast(f_full.size()) != k_full.rows()) { throw std::runtime_error("full reaction size mismatch"); @@ -69,471 +72,6 @@ inline std::vector recoverFullReaction(const DenseMatrix& k_full, const st return reaction; } -struct MITC4MaterialMatrixEvaluation { - MITC4MaterialMatrix matrix{}; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -struct MITC4StrainTransform { - MITC4MaterialMatrix matrix{}; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -struct MITC4IntegrationPoint { - Real xi = 0.0; - Real eta = 0.0; - Real zeta = 0.0; - Real weight = 0.0; -}; - -struct MITC4MaterialIntegrationSample { - MITC4IntegrationPoint point; - MITC4IntegrationBasis basis; - MITC4StrainRows strain_rows; - MITC4MaterialMatrix local_material{}; - MITC4MaterialMatrix strain_transform{}; - MITC4MaterialMatrix convected_material{}; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -struct MITC4MaterialIntegrationData { - std::vector samples; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -inline std::array mitc4GaussQuadrature2x2x2() { - const Real gauss = 1.0 / std::sqrt(3.0); - const std::array points = {-gauss, gauss}; - std::array integration_points{}; - std::size_t index = 0; - for (Real xi : points) { - for (Real eta : points) { - for (Real zeta : points) { - integration_points[index++] = {xi, eta, zeta, 1.0}; - } - } - } - return integration_points; -} - -inline MITC4MaterialMatrixEvaluation mitc4PlaneStressMaterialMatrix(Real elastic_modulus, Real poisson_ratio, - Real shear_correction = 5.0 / 6.0, - Real tolerance = 1.0e-12) { - MITC4MaterialMatrixEvaluation result; - if (!isFinite(elastic_modulus) || elastic_modulus <= tolerance) { - result.diagnostics.push_back( - mitc4Diagnostic("FESA-MITC4-MATERIAL", "MITC4 elastic modulus must be positive and finite")); - } - if (!isFinite(poisson_ratio) || poisson_ratio <= -1.0 || poisson_ratio >= 0.5) { - result.diagnostics.push_back( - mitc4Diagnostic("FESA-MITC4-POISSON", "MITC4 isotropic Poisson ratio must satisfy -1 < nu < 0.5")); - } - if (!isFinite(shear_correction) || shear_correction <= tolerance) { - result.diagnostics.push_back( - mitc4Diagnostic("FESA-MITC4-SHEAR-CORRECTION", "MITC4 shear correction factor must be positive and finite")); - } - if (hasError(result.diagnostics)) { - return result; - } - - const Real scale = elastic_modulus / (1.0 - poisson_ratio * poisson_ratio); - const Real shear_modulus = elastic_modulus / (2.0 * (1.0 + poisson_ratio)); - const std::size_t eps11 = strainComponentIndex(MITC4StrainComponent::Eps11); - const std::size_t eps22 = strainComponentIndex(MITC4StrainComponent::Eps22); - const std::size_t gamma23 = strainComponentIndex(MITC4StrainComponent::Gamma23); - const std::size_t gamma13 = strainComponentIndex(MITC4StrainComponent::Gamma13); - const std::size_t gamma12 = strainComponentIndex(MITC4StrainComponent::Gamma12); - - result.matrix[eps11][eps11] = scale; - result.matrix[eps11][eps22] = poisson_ratio * scale; - result.matrix[eps22][eps11] = poisson_ratio * scale; - result.matrix[eps22][eps22] = scale; - result.matrix[gamma23][gamma23] = shear_correction * shear_modulus; - result.matrix[gamma13][gamma13] = shear_correction * shear_modulus; - result.matrix[gamma12][gamma12] = shear_modulus; - return result; -} - -inline MITC4StrainVector multiplyMITC4MaterialMatrix(const MITC4MaterialMatrix& matrix, const MITC4StrainVector& vector) { - MITC4StrainVector result{}; - for (std::size_t row = 0; row < 6; ++row) { - for (std::size_t col = 0; col < 6; ++col) { - result[row] += matrix[row][col] * vector[col]; - } - } - return result; -} - -inline Real dotMITC4Vector(const MITC4StrainVector& a, const MITC4StrainVector& b) { - Real result = 0.0; - for (std::size_t i = 0; i < 6; ++i) { - result += a[i] * b[i]; - } - return result; -} - -inline MITC4MaterialMatrix mitc4TransformMaterialMatrix(const MITC4MaterialMatrix& local_material, - const MITC4MaterialMatrix& covariant_to_local) { - MITC4MaterialMatrix result{}; - for (std::size_t i = 0; i < 6; ++i) { - for (std::size_t j = 0; j < 6; ++j) { - Real value = 0.0; - for (std::size_t a = 0; a < 6; ++a) { - for (std::size_t b = 0; b < 6; ++b) { - value += covariant_to_local[a][i] * local_material[a][b] * covariant_to_local[b][j]; - } - } - result[i][j] = value; - } - } - return result; -} - -inline std::array, 3> mitc4TensorFromEngineeringComponent(std::size_t component) { - std::array, 3> tensor{}; - switch (static_cast(component)) { - case MITC4StrainComponent::Eps11: - tensor[0][0] = 1.0; - break; - case MITC4StrainComponent::Eps22: - tensor[1][1] = 1.0; - break; - case MITC4StrainComponent::Eps33: - tensor[2][2] = 1.0; - break; - case MITC4StrainComponent::Gamma23: - tensor[1][2] = 0.5; - tensor[2][1] = 0.5; - break; - case MITC4StrainComponent::Gamma13: - tensor[0][2] = 0.5; - tensor[2][0] = 0.5; - break; - case MITC4StrainComponent::Gamma12: - tensor[0][1] = 0.5; - tensor[1][0] = 0.5; - break; - } - return tensor; -} - -inline MITC4StrainVector mitc4EngineeringVectorFromTensor(const std::array, 3>& tensor) { - MITC4StrainVector vector{}; - vector[strainComponentIndex(MITC4StrainComponent::Eps11)] = tensor[0][0]; - vector[strainComponentIndex(MITC4StrainComponent::Eps22)] = tensor[1][1]; - vector[strainComponentIndex(MITC4StrainComponent::Eps33)] = tensor[2][2]; - vector[strainComponentIndex(MITC4StrainComponent::Gamma23)] = 2.0 * tensor[1][2]; - vector[strainComponentIndex(MITC4StrainComponent::Gamma13)] = 2.0 * tensor[0][2]; - vector[strainComponentIndex(MITC4StrainComponent::Gamma12)] = 2.0 * tensor[0][1]; - return vector; -} - -inline MITC4StrainTransform mitc4CovariantToLocalStrainTransform(const MITC4IntegrationBasis& basis, - Real tolerance = 1.0e-12) { - MITC4StrainTransform result; - result.diagnostics = basis.diagnostics; - const Real jacobian = dot(cross(basis.g1, basis.g2), basis.g3); - if (!isFinite(jacobian) || std::fabs(jacobian) <= tolerance) { - result.diagnostics.push_back( - mitc4Diagnostic("FESA-MITC4-SINGULAR-JACOBIAN", "MITC4 material transform Jacobian is near zero")); - return result; - } - if (hasError(result.diagnostics)) { - return result; - } - - const std::array contravariant = { - (1.0 / jacobian) * cross(basis.g2, basis.g3), - (1.0 / jacobian) * cross(basis.g3, basis.g1), - (1.0 / jacobian) * cross(basis.g1, basis.g2)}; - const std::array local = {basis.local.e1, basis.local.e2, basis.local.e3}; - std::array, 3> direction_cosines{}; - for (std::size_t local_axis = 0; local_axis < 3; ++local_axis) { - for (std::size_t convected_axis = 0; convected_axis < 3; ++convected_axis) { - direction_cosines[local_axis][convected_axis] = dot(local[local_axis], contravariant[convected_axis]); - } - } - - for (std::size_t column = 0; column < 6; ++column) { - const auto covariant_tensor = mitc4TensorFromEngineeringComponent(column); - std::array, 3> local_tensor{}; - for (std::size_t a = 0; a < 3; ++a) { - for (std::size_t b = 0; b < 3; ++b) { - for (std::size_t i = 0; i < 3; ++i) { - for (std::size_t j = 0; j < 3; ++j) { - local_tensor[a][b] += direction_cosines[a][i] * direction_cosines[b][j] * covariant_tensor[i][j]; - } - } - } - } - const auto local_vector = mitc4EngineeringVectorFromTensor(local_tensor); - for (std::size_t row = 0; row < 6; ++row) { - result.matrix[row][column] = local_vector[row]; - } - } - return result; -} - -inline MITC4MaterialIntegrationData mitc4BuildMaterialIntegrationData(const MITC4Geometry& geometry, Real elastic_modulus, - Real poisson_ratio, - Real shear_correction = 5.0 / 6.0) { - MITC4MaterialIntegrationData data; - const auto material = mitc4PlaneStressMaterialMatrix(elastic_modulus, poisson_ratio, shear_correction); - appendDiagnostics(data.diagnostics, material.diagnostics); - if (hasError(data.diagnostics)) { - return data; - } - - for (const MITC4IntegrationPoint& point : mitc4GaussQuadrature2x2x2()) { - MITC4MaterialIntegrationSample sample; - sample.point = point; - sample.local_material = material.matrix; - sample.basis = computeMITC4IntegrationBasis(geometry, point.xi, point.eta, point.zeta); - appendDiagnostics(sample.diagnostics, sample.basis.diagnostics); - const auto transform = mitc4CovariantToLocalStrainTransform(sample.basis); - sample.strain_transform = transform.matrix; - appendDiagnostics(sample.diagnostics, transform.diagnostics); - sample.strain_rows = mitc4TiedCovariantStrainRows(geometry, point.xi, point.eta, point.zeta); - appendDiagnostics(sample.diagnostics, sample.strain_rows.diagnostics); - if (!hasError(sample.diagnostics)) { - sample.convected_material = mitc4TransformMaterialMatrix(sample.local_material, sample.strain_transform); - } - appendDiagnostics(data.diagnostics, sample.diagnostics); - data.samples.push_back(sample); - } - return data; -} - -struct ElementStiffnessOptions { - Real drilling_stiffness_scale = 1.0e-3; -}; - -struct MITC4DrillingStabilizationResult { - DenseMatrix local_with_drilling; - Real reference_diagonal = 0.0; - Real drilling_stiffness = 0.0; - Real drilling_stiffness_scale = 0.0; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -struct MITC4ElementStiffnessResult { - DenseMatrix local_without_drilling; - DenseMatrix local_with_drilling; - DenseMatrix global; - std::size_t integration_point_count = 0; - Real drilling_reference_diagonal = 0.0; - Real drilling_stiffness = 0.0; - Real drilling_stiffness_scale = 0.0; - std::vector diagnostics; - - bool ok() const { - return !hasError(diagnostics); - } -}; - -inline DenseMatrix mitc4LocalDofTransform(const MITC4Geometry& geometry) { - DenseMatrix transform(24, 24); - for (LocalIndex node = 0; node < 4; ++node) { - const LocalIndex base = 6 * node; - const auto& frame = geometry.nodal_frames[static_cast(node)]; - const std::array axes = {frame.v1, frame.v2, frame.vn}; - for (LocalIndex local_axis = 0; local_axis < 3; ++local_axis) { - const Vec3 axis = axes[static_cast(local_axis)]; - transform(base + local_axis, base + 0) = axis.x; - transform(base + local_axis, base + 1) = axis.y; - transform(base + local_axis, base + 2) = axis.z; - transform(base + 3 + local_axis, base + 3) = axis.x; - transform(base + 3 + local_axis, base + 4) = axis.y; - transform(base + 3 + local_axis, base + 5) = axis.z; - } - } - return transform; -} - -inline MITC4StrainRows mitc4TransformStrainRowsToLocalDofs(const MITC4StrainRows& global_rows, - const DenseMatrix& local_dof_transform) { - MITC4StrainRows local_rows; - local_rows.diagnostics = global_rows.diagnostics; - if (hasError(local_rows.diagnostics)) { - return local_rows; - } - for (std::size_t component = 0; component < 6; ++component) { - for (LocalIndex local_dof = 0; local_dof < 24; ++local_dof) { - Real value = 0.0; - for (LocalIndex global_dof = 0; global_dof < 24; ++global_dof) { - value += global_rows.rows[component][static_cast(global_dof)] * - local_dof_transform(local_dof, global_dof); - } - local_rows.rows[component][static_cast(local_dof)] = value; - } - } - return local_rows; -} - -inline void accumulateMITC4BtDB(DenseMatrix& stiffness, const MITC4StrainRows& rows, - const MITC4MaterialMatrix& material, Real factor) { - for (LocalIndex i = 0; i < 24; ++i) { - for (LocalIndex j = 0; j < 24; ++j) { - Real value = 0.0; - for (std::size_t a = 0; a < 6; ++a) { - for (std::size_t b = 0; b < 6; ++b) { - value += rows.rows[a][static_cast(i)] * material[a][b] * - rows.rows[b][static_cast(j)]; - } - } - stiffness.add(i, j, value * factor); - } - } -} - -inline Real mitc4MinimumPositivePhysicalLocalDiagonal(const DenseMatrix& local_without_drilling, - Real tolerance = 1.0e-12) { - Real minimum = std::numeric_limits::infinity(); - for (LocalIndex node = 0; node < 4; ++node) { - for (LocalIndex local_dof = 0; local_dof < 5; ++local_dof) { - const Real diagonal = local_without_drilling(6 * node + local_dof, 6 * node + local_dof); - if (isFinite(diagonal) && diagonal > tolerance) { - minimum = std::min(minimum, diagonal); - } - } - } - return minimum; -} - -inline MITC4DrillingStabilizationResult mitc4ApplyDrillingStabilization(const DenseMatrix& local_without_drilling, - Real drilling_stiffness_scale, - Real tolerance = 1.0e-12) { - MITC4DrillingStabilizationResult result; - result.local_with_drilling = local_without_drilling; - result.drilling_stiffness_scale = drilling_stiffness_scale; - if (!isFinite(drilling_stiffness_scale) || drilling_stiffness_scale < 0.0) { - result.diagnostics.push_back( - mitc4Diagnostic("FESA-MITC4-DRILLING-SCALE", "MITC4 drilling stiffness scale must be non-negative and finite")); - return result; - } - - result.reference_diagonal = mitc4MinimumPositivePhysicalLocalDiagonal(local_without_drilling, tolerance); - if (!isFinite(result.reference_diagonal)) { - result.diagnostics.push_back(mitc4Diagnostic("FESA-MITC4-DRILLING-REFERENCE", - "MITC4 drilling stiffness reference diagonal is not positive")); - return result; - } - - result.drilling_stiffness = drilling_stiffness_scale * result.reference_diagonal; - for (LocalIndex node = 0; node < 4; ++node) { - const LocalIndex gamma = 6 * node + 5; - result.local_with_drilling.add(gamma, gamma, result.drilling_stiffness); - } - return result; -} - -inline DenseMatrix mitc4TransformLocalStiffnessToGlobal(const DenseMatrix& local_stiffness, - const DenseMatrix& local_dof_transform) { - DenseMatrix global(24, 24); - for (LocalIndex i = 0; i < 24; ++i) { - for (LocalIndex j = 0; j < 24; ++j) { - Real value = 0.0; - for (LocalIndex a = 0; a < 24; ++a) { - for (LocalIndex b = 0; b < 24; ++b) { - value += local_dof_transform(a, i) * local_stiffness(a, b) * local_dof_transform(b, j); - } - } - global(i, j) = value; - } - } - return global; -} - -inline MITC4ElementStiffnessResult mitc4ElementStiffness(const std::array& coordinates, - Real elastic_modulus, Real poisson_ratio, Real thickness, - ElementStiffnessOptions options = {}) { - MITC4ElementStiffnessResult result; - result.local_without_drilling = DenseMatrix(24, 24); - result.local_with_drilling = DenseMatrix(24, 24); - result.global = DenseMatrix(24, 24); - result.drilling_stiffness_scale = options.drilling_stiffness_scale; - - const MITC4Geometry geometry = buildMITC4Geometry(coordinates, thickness); - appendDiagnostics(result.diagnostics, geometry.diagnostics); - if (hasError(result.diagnostics)) { - return result; - } - - const MITC4MaterialIntegrationData integration = - mitc4BuildMaterialIntegrationData(geometry, elastic_modulus, poisson_ratio); - appendDiagnostics(result.diagnostics, integration.diagnostics); - if (hasError(result.diagnostics)) { - return result; - } - result.integration_point_count = integration.samples.size(); - - const DenseMatrix local_dof_transform = mitc4LocalDofTransform(geometry); - for (const MITC4MaterialIntegrationSample& sample : integration.samples) { - const MITC4StrainRows local_rows = - mitc4TransformStrainRowsToLocalDofs(sample.strain_rows, local_dof_transform); - appendDiagnostics(result.diagnostics, local_rows.diagnostics); - if (hasError(result.diagnostics)) { - return result; - } - const Real factor = std::fabs(sample.basis.jacobian) * sample.point.weight; - accumulateMITC4BtDB(result.local_without_drilling, local_rows, sample.convected_material, factor); - } - - const auto drilling = mitc4ApplyDrillingStabilization(result.local_without_drilling, options.drilling_stiffness_scale); - appendDiagnostics(result.diagnostics, drilling.diagnostics); - result.local_with_drilling = drilling.local_with_drilling; - result.drilling_reference_diagonal = drilling.reference_diagonal; - result.drilling_stiffness = drilling.drilling_stiffness; - result.drilling_stiffness_scale = drilling.drilling_stiffness_scale; - if (hasError(result.diagnostics)) { - return result; - } - - result.global = mitc4TransformLocalStiffnessToGlobal(result.local_with_drilling, local_dof_transform); - return result; -} - -inline std::vector mitc4ElementInternalForce(const MITC4ElementStiffnessResult& stiffness, - const std::vector& element_displacement) { - if (stiffness.global.rows() != static_cast(element_displacement.size())) { - throw std::runtime_error("MITC4 internal force size mismatch"); - } - return stiffness.global.multiply(element_displacement); -} - -class MITC4ElementKernel { - public: - DenseMatrix stiffness(const std::array& coordinates, Real elastic_modulus, Real poisson_ratio, Real thickness, - ElementStiffnessOptions options = {}) const { - const auto result = mitc4ElementStiffness(coordinates, elastic_modulus, poisson_ratio, thickness, options); - if (!result.ok()) { - throw std::runtime_error(result.diagnostics.empty() ? "invalid MITC4 stiffness" - : result.diagnostics.front().message); - } - return result.global; - } -}; - struct AssemblyResult { DenseMatrix k_full; std::vector f_full; @@ -568,12 +106,14 @@ inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dof for (const auto& [element_id, element] : domain.elements) { const ShellSection* section = shellSectionForElement(domain, element_id); if (section == nullptr) { - result.diagnostics.push_back({Severity::Error, "FESA-ASSEMBLY-MISSING-PROPERTY", "Element has no shell section", {}}); + result.diagnostics.push_back( + {Severity::Error, "FESA-ASSEMBLY-MISSING-PROPERTY", "Element has no shell section", {}}); continue; } const auto material_it = domain.materials.find(Domain::key(section->material)); if (material_it == domain.materials.end()) { - result.diagnostics.push_back({Severity::Error, "FESA-ASSEMBLY-MISSING-MATERIAL", "Element material is missing", {}}); + result.diagnostics.push_back( + {Severity::Error, "FESA-ASSEMBLY-MISSING-MATERIAL", "Element material is missing", {}}); continue; } std::array coordinates{}; @@ -581,7 +121,7 @@ inline AssemblyResult assembleSystem(const Domain& domain, const DofManager& dof coordinates[i] = domain.nodes.at(element.node_ids[i]).coordinates; } const auto stiffness = mitc4ElementStiffness(coordinates, material_it->second.elastic_modulus, - material_it->second.poisson_ratio, section->thickness, options); + material_it->second.poisson_ratio, section->thickness, options); result.diagnostics.insert(result.diagnostics.end(), stiffness.diagnostics.begin(), stiffness.diagnostics.end()); if (!stiffness.ok()) { continue; @@ -614,8 +154,7 @@ inline ReducedSystem projectToReducedSystem(const AssemblyResult& assembly, cons result.k = DenseMatrix(dofs.freeDofCount(), dofs.freeDofCount()); result.f = std::vector(static_cast(dofs.freeDofCount()), 0.0); result.free_full_indices = dofs.freeFullIndices(); - if (assembly.k_full.rows() != assembly.k_full.cols() || - assembly.k_full.rows() != dofs.fullDofCount() || + if (assembly.k_full.rows() != assembly.k_full.cols() || assembly.k_full.rows() != dofs.fullDofCount() || static_cast(assembly.f_full.size()) != dofs.fullDofCount()) { result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-ASSEMBLY-SIZE", "Full-space stiffness/load sizes do not match DofManager", "assembly")); @@ -628,7 +167,8 @@ inline ReducedSystem projectToReducedSystem(const AssemblyResult& assembly, cons } if (assembly.reduced_pattern.equation_count != dofs.freeDofCount() || assembly.reduced_pattern.nonzeroCount() == 0) { result.diagnostics.push_back(makeDiagnostic(Severity::Error, "FESA-SINGULAR-SPARSE-PATTERN", - "Reduced sparse pattern is empty or inconsistent with free DOFs", "assembly")); + "Reduced sparse pattern is empty or inconsistent with free DOFs", + "assembly")); return result; } for (LocalIndex i = 0; i < dofs.freeDofCount(); ++i) { @@ -656,6 +196,7 @@ struct AnalysisResult { class Analysis { public: virtual ~Analysis() = default; + AnalysisResult run(const Domain& domain) const { AnalysisResult result; initialize(domain, result); @@ -671,6 +212,7 @@ class Analysis { auto diagnostics = validateDomain(domain); result.diagnostics.insert(result.diagnostics.end(), diagnostics.begin(), diagnostics.end()); } + virtual void solve(const Domain& domain, AnalysisResult& result) const = 0; }; diff --git a/phases/1-structure-alignment-refactor/index.json b/phases/1-structure-alignment-refactor/index.json index 4afd2d2..dbaf8c9 100644 --- a/phases/1-structure-alignment-refactor/index.json +++ b/phases/1-structure-alignment-refactor/index.json @@ -9,7 +9,7 @@ { "step": 4, "name": "io-parser-extraction", "status": "completed" }, { "step": 5, "name": "results-reference-extraction", "status": "completed" }, { "step": 6, "name": "mitc4-geometry-strain-extraction", "status": "completed" }, - { "step": 7, "name": "mitc4-material-stiffness-extraction", "status": "pending" }, + { "step": 7, "name": "mitc4-material-stiffness-extraction", "status": "completed" }, { "step": 8, "name": "assembly-analysis-extraction", "status": "pending" }, { "step": 9, "name": "architecture-evaluator-closeout", "status": "pending" } ] diff --git a/tests/test_mitc4_stiffness_module_includes.cpp b/tests/test_mitc4_stiffness_module_includes.cpp new file mode 100644 index 0000000..b6894f9 --- /dev/null +++ b/tests/test_mitc4_stiffness_module_includes.cpp @@ -0,0 +1,117 @@ +#include "fesa/Element/Element.hpp" +#include "fesa/Element/MITC4Stiffness.hpp" +#include "fesa/Material/MITC4PlaneStressMaterial.hpp" +#include "fesa/Material/Material.hpp" + +#include +#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); + } +} + +fesa::Real matrixEnergy(const fesa::DenseMatrix& matrix, const std::vector& values) { + const auto product = matrix.multiply(values); + fesa::Real energy = 0.0; + for (std::size_t i = 0; i < values.size(); ++i) { + energy += values[i] * product[i]; + } + return energy; +} + +} // namespace + +int main() { + const auto law = fesa::mitc4PlaneStressMaterialMatrix(1000.0, 0.25); + check(law.ok(), "valid MITC4 plane-stress material should remain accepted"); + check(fesa::mitc4StrainComponentLabels()[0] == "eps11", "MITC4 strain labels changed"); + const auto eps11 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps11); + const auto eps22 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Eps22); + const auto gamma13 = fesa::strainComponentIndex(fesa::MITC4StrainComponent::Gamma13); + checkNear(law.matrix[eps11][eps11], 1066.6666666666667, 1.0e-10, "plane-stress D11 changed"); + checkNear(law.matrix[eps11][eps22], 266.6666666666667, 1.0e-10, "plane-stress D12 changed"); + checkNear(law.matrix[gamma13][gamma13], (5.0 / 6.0) * 400.0, 1.0e-10, + "plane-stress shear correction changed"); + + fesa::MITC4StrainVector strain{}; + strain[eps11] = 0.01; + strain[eps22] = -0.02; + strain[gamma13] = 0.03; + const auto stress = fesa::multiplyMITC4MaterialMatrix(law.matrix, strain); + check(fesa::dotMITC4Vector(strain, stress) > 0.0, "MITC4 material energy changed"); + + const auto invalid = fesa::mitc4PlaneStressMaterialMatrix(-1.0, 0.25); + check(!invalid.ok(), "invalid MITC4 material should remain diagnosed"); + check(fesa::containsDiagnostic(invalid.diagnostics, "FESA-MITC4-MATERIAL"), + "MITC4 invalid material diagnostic changed"); + + const auto points = fesa::mitc4GaussQuadrature2x2x2(); + check(points.size() == 8, "MITC4 integration point count changed"); + for (const auto& point : points) { + checkNear(point.weight, 1.0, 1.0e-15, "MITC4 integration weight 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 MITC4 geometry should remain valid"); + const auto integration = fesa::mitc4BuildMaterialIntegrationData(geometry, 1000.0, 0.25); + check(integration.ok(), "MITC4 material integration data should remain valid"); + check(integration.samples.size() == 8, "MITC4 material integration sample count changed"); + for (const auto& sample : integration.samples) { + check(sample.ok(), "MITC4 material integration sample should remain valid"); + check(sample.basis.ok(), "MITC4 material integration basis should remain valid"); + } + + const auto stiffness = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2); + check(stiffness.ok(), "MITC4 stiffness should remain valid"); + check(stiffness.global.rows() == 24 && stiffness.global.cols() == 24, "MITC4 global stiffness size changed"); + check(stiffness.integration_point_count == 8, "MITC4 stiffness integration count changed"); + checkNear(stiffness.drilling_stiffness_scale, 1.0e-3, 1.0e-15, "MITC4 drilling scale changed"); + check(stiffness.drilling_reference_diagonal > 0.0, "MITC4 drilling reference diagonal changed"); + check(stiffness.drilling_stiffness > 0.0, "MITC4 drilling stiffness changed"); + for (fesa::LocalIndex i = 0; i < 24; ++i) { + for (fesa::LocalIndex j = 0; j < 24; ++j) { + checkNear(stiffness.global(i, j), stiffness.global(j, i), 1.0e-8, "MITC4 stiffness symmetry changed"); + } + } + + std::vector displacement(24, 0.0); + displacement[0] = 0.01; + displacement[7] = -0.02; + displacement[14] = 0.03; + displacement[21] = 0.04; + const auto internal_force = fesa::mitc4ElementInternalForce(stiffness, displacement); + check(internal_force.size() == displacement.size(), "MITC4 internal force size changed"); + const auto expected_internal_force = stiffness.global.multiply(displacement); + for (std::size_t i = 0; i < internal_force.size(); ++i) { + checkNear(internal_force[i], expected_internal_force[i], 1.0e-12, "MITC4 internal force changed"); + } + check(matrixEnergy(stiffness.global, displacement) > 0.0, "MITC4 stiffness energy changed"); + + fesa::ElementStiffnessOptions options; + options.drilling_stiffness_scale = 2.0e-3; + const auto scaled = fesa::mitc4ElementStiffness(coords, 1000.0, 0.25, 0.2, options); + check(scaled.ok(), "scaled MITC4 stiffness should remain valid"); + checkNear(scaled.drilling_reference_diagonal, stiffness.drilling_reference_diagonal, 1.0e-10, + "MITC4 drilling reference changed with scale"); + checkNear(scaled.drilling_stiffness, 2.0 * stiffness.drilling_stiffness, 1.0e-10, + "MITC4 drilling scale response changed"); + + fesa::MITC4ElementKernel kernel; + const auto kernel_stiffness = kernel.stiffness(coords, 1000.0, 0.25, 0.2); + check(kernel_stiffness.rows() == 24 && kernel_stiffness.cols() == 24, "MITC4 kernel stiffness size changed"); + + return 0; +}