From 675379779fe07ba329dc91ca768de8c997b2ac01 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Mon, 3 Aug 2026 01:45:30 +0900 Subject: [PATCH] feat(beam-reference-qualification): correlate Abaqus beam results --- AGENTS.md | 2 +- CMakeLists.txt | 10 + docs/ADR.md | 4 +- docs/ARCHITECTURE.md | 13 +- include/fesa/validation/comparison.hpp | 18 ++ phases/beam-reference-qualification/step2.md | 2 +- phases/beam-reference-qualification/step3.md | 15 +- .../cantilever beam displacements.csv | 24 +- .../cantilever beam elemental forces.csv | 21 ++ .../{Job-1.inp => cantilever beam fesa.inp} | 51 +-- .../cantilever beam reactions.csv | 24 +- reference/cantilever beam/cantilever beam.inp | 47 +-- src/fesa/validation/comparison.cpp | 149 +++++++++ .../validation/reference_compare_main.cpp | 305 ++++++++++++++++++ src/fesa/validation/reference_csv.cpp | 14 +- tests/CMakeLists.txt | 8 + .../io/minimal_deck_to_domain_test.cpp | 4 +- tests/reference/CMakeLists.txt | 38 +++ tests/reference/cantilever_reference_test.cpp | 241 ++++++++++++++ tests/unit/validation/comparison_test.cpp | 60 ++++ tests/unit/validation/reference_csv_test.cpp | 31 +- 21 files changed, 962 insertions(+), 119 deletions(-) create mode 100644 reference/cantilever beam/cantilever beam elemental forces.csv rename reference/cantilever beam/{Job-1.inp => cantilever beam fesa.inp} (63%) create mode 100644 src/fesa/validation/reference_compare_main.cpp create mode 100644 tests/reference/CMakeLists.txt create mode 100644 tests/reference/cantilever_reference_test.cpp diff --git a/AGENTS.md b/AGENTS.md index 90e3a6a..8dedd79 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -60,7 +60,7 @@ - 전단강성이 생략되면 \(A_{sy}=A_{sz}=5A/6\)과 `SCF=0`을 Phase 1 기본값으로 적용한다. - reference 비교는 metadata 없이 요청된 물리량과 CSV 경로를 명시한다. 현재 - 캔틸레버 샘플은 변위와 반력을 비교하며 요소 내력과 도심 응력 비교 루틴은 + 캔틸레버 샘플은 변위, 반력 및 요소 단면력을 비교하며 도심 응력 비교 루틴은 synthetic CSV로 검증한다. - 새 MSVC 빌드 경고를 추가하지 않는다. - 변경은 요청 범위에 한정하고 Conventional Commits 형식의 메시지를 사용한다. diff --git a/CMakeLists.txt b/CMakeLists.txt index 5254573..03ab119 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,6 +68,16 @@ target_link_libraries(fesa PRIVATE fesa_core) target_compile_features(fesa PRIVATE cxx_std_20) target_compile_options(fesa PRIVATE /W4 /permissive- /EHsc) +add_executable(fesa-reference-compare + src/fesa/validation/reference_compare_main.cpp +) + +target_link_libraries(fesa-reference-compare PRIVATE fesa_core) +target_compile_features(fesa-reference-compare PRIVATE cxx_std_20) +target_compile_options( + fesa-reference-compare PRIVATE /W4 /permissive- /EHsc +) + if(BUILD_TESTING) add_executable(fesa_assembly_benchmark tests/performance/assembly_benchmark.cpp diff --git a/docs/ADR.md b/docs/ADR.md index 8da0769..d9cc7a6 100644 --- a/docs/ADR.md +++ b/docs/ADR.md @@ -164,7 +164,7 @@ HDF5 파일에 저장하고 root에 schema version을 기록한다. ## ADR-010: 검증 계층별 허용오차 -**상태:** Accepted +**상태:** Accepted for threshold tests; cross-formulation correlation superseded by ADR-017 **상황:** 모든 물리량에 하나의 상대오차를 적용하면 영에 가까운 값이나 서로 다른 규모의 결과를 올바르게 판정할 수 없다. @@ -177,6 +177,8 @@ HDF5 파일에 저장하고 root에 schema version을 기록한다. - 영에 가까운 값과 큰 값 모두 의미 있게 비교할 수 있다. - 모델별 예외 tolerance에는 문서화된 수치 근거가 필요하다. - 단일 tolerance보다 comparison request와 helper가 복잡해진다. +- Abaqus와 FESA의 정식화가 다른 상관성 비교에는 ADR-017의 component별 RMSE와 + Relative L2 계약을 적용한다. ## ADR-011: 일관 단위계와 결과 좌표계 diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 98524f7..c961e3f 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -429,16 +429,17 @@ kernel을 추가한다. - `tests/reference`: CSV 골든 결과와 FESA HDF5 결과 비교 - `reference/`: Abaqus 입력과 현재 사용할 수 있는 결과 CSV -reference comparison request가 비교할 물리량과 CSV 경로, 상대 tolerance 및 -물리량별 절대 scale을 명시한다. 요청한 CSV가 없으면 실패하며 요청하지 않은 결과를 -통과로 표시하지 않는다. 현재 캔틸레버는 변위와 반력만 요청하고, 요소 내력과 -단면 도심 응력 adapter는 synthetic CSV로 검증한다. +reference comparison request가 비교할 물리량과 CSV 경로 및 물리량별 절대 scale을 +명시한다. 요청한 CSV가 없으면 실패하며 요청하지 않은 결과를 통과로 표시하지 +않는다. 현재 캔틸레버는 변위, 반력 및 요소 단면력을 요청하고, 단면 도심 응력 +adapter는 synthetic CSV로 검증한다. 요소 내력 CSV의 `(Instance, Element Label, Node Label)` 위치에서 -`SF1,SF2,SF3,SM1,SM2,SM3`을 \(N,V_y,V_z,T,M_y,M_z\)로 매핑한다. 응력 CSV의 +`SF1,SF3,SF2,SM3,SM1,SM2`를 \(N,V_y,V_z,T,M_y,M_z\)로 매핑한다. 응력 CSV의 같은 위치에 있는 `Sxx`는 단면 도심값 \(N/A\)와 비교한다. 단일 Instance에서는 Instance 열 생략을 허용하되 comparison request가 제공한 Instance 이름으로 -보완한다. +보완한다. Abaqus와 FESA의 정식화가 다른 상관성 비교는 component별 RMSE와 +Relative L2를 보고하며 관측값으로 만든 pass/fail tolerance를 적용하지 않는다. reference helper는 반드시 public parser와 analysis 경로로 FESA 결과를 생성한다. 테스트 전용 경로로 Domain이나 matrix를 직접 주입해 전체 파이프라인 결함을 숨기지 diff --git a/include/fesa/validation/comparison.hpp b/include/fesa/validation/comparison.hpp index 963504c..01b48c1 100644 --- a/include/fesa/validation/comparison.hpp +++ b/include/fesa/validation/comparison.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include #include @@ -43,6 +44,20 @@ struct ComparisonReport final { std::vector failures; }; +struct ComponentCorrelationMetric final { + ReferenceQuantity quantity; + std::size_t component_index; + std::size_t value_count; + double root_mean_square_error; + double relative_l2_error; +}; + +struct CorrelationReport final { + bool evaluable; + std::vector metrics; + std::vector failures; +}; + struct ComparisonSampleMatch final { std::optional sample; std::vector failures; @@ -58,4 +73,7 @@ struct ComparisonSampleMatch final { [[nodiscard]] ComparisonReport compare_samples( std::span samples); +[[nodiscard]] CorrelationReport correlate_samples( + std::span samples); + } // namespace fesa diff --git a/phases/beam-reference-qualification/step2.md b/phases/beam-reference-qualification/step2.md index 8622057..31658f7 100644 --- a/phases/beam-reference-qualification/step2.md +++ b/phases/beam-reference-qualification/step2.md @@ -52,7 +52,7 @@ FESA 정식화 적합성과 Abaqus 결과 상관성을 별도 gate로 검증한 ```powershell cmake --build --preset windows-debug -ctest --preset windows-debug -R "ValidationComparison|ReferenceCsv" --output-on-failure +ctest --preset windows-debug -R "CorrelationMetric|ReferenceCsv|InternalForceCsv" --output-on-failure ctest --preset windows-debug -R CantileverReference --output-on-failure .\out\build\windows-debug\Debug\fesa.exe solve "reference\cantilever beam\cantilever beam fesa.inp" --output out\cantilever-beam.h5 .\out\build\windows-debug\Debug\fesa-reference-compare.exe --results out\cantilever-beam.h5 --instance PART-1_1-1 --displacements "reference\cantilever beam\cantilever beam displacements.csv" --reactions "reference\cantilever beam\cantilever beam reactions.csv" --internal-forces "reference\cantilever beam\cantilever beam elemental forces.csv" --displacement-absolute-scale 1e-10 --reaction-absolute-scale 1e-8 --internal-force-absolute-scale 1e-8 diff --git a/phases/beam-reference-qualification/step3.md b/phases/beam-reference-qualification/step3.md index af1cc4d..6419f09 100644 --- a/phases/beam-reference-qualification/step3.md +++ b/phases/beam-reference-qualification/step3.md @@ -20,9 +20,10 @@ rotated frame, slenderness sweep - physics: equilibrium, symmetry, reaction, nonzero prescribed DOF - determinism: tested thread counts와 repeated runs -- Abaqus: 현재 cantilever displacement/reaction의 tolerance와 maximum error -- contract-only: synthetic internal-force/stress CSV schema와 component mapping -- 미제공 Abaqus internal-force/stress는 `not yet Abaqus-qualified`라고 명시한다. +- Abaqus correlation: 현재 cantilever displacement/reaction/internal-force의 + component별 RMSE와 Relative L2 +- contract-only: synthetic stress CSV schema와 component mapping +- 미제공 Abaqus stress는 `not yet Abaqus-qualified`라고 명시한다. 보고서 수치를 새 test output에서 수집하며 수동 추정값을 쓰지 않는다. @@ -34,18 +35,18 @@ ctest --preset windows-debug --output-on-failure ctest --preset windows-debug -R "Reference|Beam3D2|Determinism" --output-on-failure ``` -모든 실행이 통과하고 보고서의 test 이름, tolerance, 최대오차와 disposition이 실제 -출력과 일치해야 한다. +모든 실행이 통과하고 보고서의 test 이름, component별 RMSE·Relative L2와 +disposition이 실제 출력과 일치해야 한다. ## 검증 절차 1. 전체 suite를 새로 실행한다. 2. 결과를 benchmark/quantity별 표에 기록한다. -3. synthetic coverage와 Abaqus-backed qualification을 명확히 분리한다. +3. FESA 정식화 gate, Abaqus correlation 및 synthetic coverage를 명확히 분리한다. 4. index summary에 보고서 경로와 test counts를 기록한다. ## 금지사항 - 실행하지 않은 결과를 보고서에 쓰지 마라. 이유: 검증 증거가 아니다. -- Abaqus 내력/응력 qualification을 추론하지 마라. 이유: CSV가 아직 없다. +- Abaqus 응력 qualification을 추론하지 마라. 이유: CSV가 아직 없다. - 실패 테스트를 제외하거나 disable하지 마라. 이유: release gate를 약화한다. diff --git a/reference/cantilever beam/cantilever beam displacements.csv b/reference/cantilever beam/cantilever beam displacements.csv index 2eb9c64..24ae5db 100644 --- a/reference/cantilever beam/cantilever beam displacements.csv +++ b/reference/cantilever beam/cantilever beam displacements.csv @@ -1,12 +1,12 @@ - Node Label, U-U1, U-U2, U-U3, UR-UR1, UR-UR2, UR-UR3 -1,0,0,-1.00E-32,0,1.00E-31,0 -2,0,0,-2.87E-06,0,5.43E-06,0 -3,0,0,-1.09E-05,0,1.03E-05,0 -4,0,0,-2.35E-05,0,1.46E-05,0 -5,0,0,-4.00E-05,0,1.83E-05,0 -6,0,0,-6.01E-05,0,2.14E-05,0 -7,0,0,-8.29E-05,0,2.40E-05,0 -8,0,0,-1.08E-04,0,2.60E-05,0 -9,0,0,-1.35E-04,0,2.74E-05,0 -10,0,0,-1.63E-04,0,2.83E-05,0 -11,0,0,-1.92E-04,0,2.86E-05,0 +Part Instance Name, Node Label, U-U1, U-U2, U-U3, UR-UR1, UR-UR2, UR-UR3 +PART-1_1-1,1,0.000000E+00,0.000000E+00,-1.000000E-30,0.000000E+00,1.000000E-29,0.000000E+00 +PART-1_1-1,2,0.000000E+00,0.000000E+00,-2.900000E-04,0.000000E+00,5.428570E-04,0.000000E+00 +PART-1_1-1,3,0.000000E+00,0.000000E+00,-1.094290E-03,0.000000E+00,1.028570E-03,0.000000E+00 +PART-1_1-1,4,0.000000E+00,0.000000E+00,-2.355720E-03,0.000000E+00,1.457140E-03,0.000000E+00 +PART-1_1-1,5,0.000000E+00,0.000000E+00,-4.017140E-03,0.000000E+00,1.828570E-03,0.000000E+00 +PART-1_1-1,6,0.000000E+00,0.000000E+00,-6.021430E-03,0.000000E+00,2.142860E-03,0.000000E+00 +PART-1_1-1,7,0.000000E+00,0.000000E+00,-8.311430E-03,0.000000E+00,2.400000E-03,0.000000E+00 +PART-1_1-1,8,0.000000E+00,0.000000E+00,-1.083000E-02,0.000000E+00,2.600000E-03,0.000000E+00 +PART-1_1-1,9,0.000000E+00,0.000000E+00,-1.352000E-02,0.000000E+00,2.742860E-03,0.000000E+00 +PART-1_1-1,10,0.000000E+00,0.000000E+00,-1.632430E-02,0.000000E+00,2.828570E-03,0.000000E+00 +PART-1_1-1,11,0.000000E+00,0.000000E+00,-1.918570E-02,0.000000E+00,2.857140E-03,0.000000E+00 diff --git a/reference/cantilever beam/cantilever beam elemental forces.csv b/reference/cantilever beam/cantilever beam elemental forces.csv new file mode 100644 index 0000000..9720e64 --- /dev/null +++ b/reference/cantilever beam/cantilever beam elemental forces.csv @@ -0,0 +1,21 @@ +Part Instance Name, Element Label, Node Label, SF-SF1, SF-SF2, SF-SF3, SM-SM1, SM-SM2, SM-SM3,,,,, +PART-1_1-1,1,1,0.000000E+00,-1.000000E+06,0.000000E+00,9.500000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,1,2,0.000000E+00,-1.000000E+06,0.000000E+00,9.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,2,2,0.000000E+00,-1.000000E+06,0.000000E+00,9.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,2,3,0.000000E+00,-1.000000E+06,0.000000E+00,8.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,3,3,0.000000E+00,-1.000000E+06,0.000000E+00,8.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,3,4,0.000000E+00,-1.000000E+06,0.000000E+00,7.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,4,4,0.000000E+00,-1.000000E+06,0.000000E+00,7.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,4,5,0.000000E+00,-1.000000E+06,0.000000E+00,6.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,5,5,0.000000E+00,-1.000000E+06,0.000000E+00,6.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,5,6,0.000000E+00,-1.000000E+06,0.000000E+00,5.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,6,6,0.000000E+00,-1.000000E+06,0.000000E+00,5.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,6,7,0.000000E+00,-1.000000E+06,0.000000E+00,4.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,7,7,0.000000E+00,-1.000000E+06,0.000000E+00,4.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,7,8,0.000000E+00,-1.000000E+06,0.000000E+00,3.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,8,8,0.000000E+00,-1.000000E+06,0.000000E+00,3.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,8,9,0.000000E+00,-1.000000E+06,0.000000E+00,2.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,9,9,0.000000E+00,-1.000000E+06,0.000000E+00,2.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,9,10,0.000000E+00,-1.000000E+06,0.000000E+00,1.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,10,10,0.000000E+00,-1.000000E+06,0.000000E+00,1.000000E+06,0.000000E+00,0.000000E+00,,,,, +PART-1_1-1,10,11,0.000000E+00,-1.000000E+06,0.000000E+00,5.000000E+05,0.000000E+00,0.000000E+00,,,,, diff --git a/reference/cantilever beam/Job-1.inp b/reference/cantilever beam/cantilever beam fesa.inp similarity index 63% rename from reference/cantilever beam/Job-1.inp rename to reference/cantilever beam/cantilever beam fesa.inp index 700fab7..8b0f975 100644 --- a/reference/cantilever beam/Job-1.inp +++ b/reference/cantilever beam/cantilever beam fesa.inp @@ -1,10 +1,8 @@ *Heading -** Job name: Job-1 Model name: Model-1 -** Generated by: Abaqus/CAE Learning Edition 2024 +** FESA projection of the Abaqus 2024 cantilever reference model. +** Geometry, material, section, boundary conditions, and load are unchanged. +** SCF is set to zero because FESA does not use Abaqus slenderness compensation. *Preprint, echo=NO, model=NO, history=NO, contact=NO -** -** PARTS -** *Part, name=PART-1_1 *Node 1, 0., 0., 0. @@ -31,37 +29,23 @@ 10, 10, 11 *Elset, elset=Set-1, generate 1, 10, 1 -*Elset, elset=Set-2, generate - 1, 10, 1 -** Section: Section-1 Profile: Profile-1 *Beam General Section, elset=Set-1, material=Material-1, section=GENERAL 1., 0.0833333, 0., 0.0833333, 0.140833 0.,1.,0. +*Transverse Shear Stiffness +6.73077e+10, 6.73077e+10, 0. *End Part -** -** -** ASSEMBLY -** *Assembly, name=Assembly -** *Instance, name=PART-1_1-1, part=PART-1_1 *End Instance -** *Nset, nset=Set-3, instance=PART-1_1-1 1, *Nset, nset=Set-4, instance=PART-1_1-1 11, *End Assembly -** -** MATERIALS -** *Material, name=Material-1 *Elastic 2.1e+11, 0.3 -** -** BOUNDARY CONDITIONS -** -** Name: BC-1 Type: Displacement/Rotation *Boundary Set-3, 1, 1 Set-3, 2, 2 @@ -69,34 +53,9 @@ Set-3, 3, 3 Set-3, 4, 4 Set-3, 5, 5 Set-3, 6, 6 -** ---------------------------------------------------------------- -** -** STEP: Step-1 -** *Step, name=Step-1, nlgeom=NO *Static 1., 1., 1e-05, 1. -** -** LOADS -** -** Name: Load-1 Type: Concentrated force *Cload Set-4, 3, -1e+06 -** -** OUTPUT REQUESTS -** -*Restart, write, frequency=0 -** -** FIELD OUTPUT: F-Output-1 -** -*Output, field -*Node Output -CF, RF, TF, U -*Element Output, directions=YES -LE, NFORC, NFORCSO, PE, PEEQ, PEMAG, S, SF -*Contact Output, variable=PRESELECT -** -** HISTORY OUTPUT: H-Output-1 -** -*Output, history, variable=PRESELECT *End Step diff --git a/reference/cantilever beam/cantilever beam reactions.csv b/reference/cantilever beam/cantilever beam reactions.csv index 245e7cc..83584a1 100644 --- a/reference/cantilever beam/cantilever beam reactions.csv +++ b/reference/cantilever beam/cantilever beam reactions.csv @@ -1,12 +1,12 @@ - Node Label, RF-RF1, RF-RF2, RF-RF3, RM-RM1, RM-RM2, RM-RM3 -1,0,0,1.00E+04,0,-1.00E+05,0 -2,0,0,0,0,0,0 -3,0,0,0,0,0,0 -4,0,0,0,0,0,0 -5,0,0,0,0,0,0 -6,0,0,0,0,0,0 -7,0,0,0,0,0,0 -8,0,0,0,0,0,0 -9,0,0,0,0,0,0 -10,0,0,0,0,0,0 -11,0,0,0,0,0,0 +Part Instance Name, Node Label, RF-RF1, RF-RF2, RF-RF3, RM-RM1, RM-RM2, RM-RM3 +PART-1_1-1,1,0.000000E+00,0.000000E+00,1.000000E+06,0.000000E+00,-1.000000E+07,0.000000E+00 +PART-1_1-1,2,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,3,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,4,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,5,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,6,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,7,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,8,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,9,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,10,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +PART-1_1-1,11,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 diff --git a/reference/cantilever beam/cantilever beam.inp b/reference/cantilever beam/cantilever beam.inp index 60d62b7..9883654 100644 --- a/reference/cantilever beam/cantilever beam.inp +++ b/reference/cantilever beam/cantilever beam.inp @@ -5,7 +5,7 @@ ** ** PARTS ** -*Part, name=Part-1 +*Part, name=PART-1_1 *Node 1, 0., 0., 0. 2, 1., 0., 0. @@ -29,18 +29,16 @@ 8, 8, 9 9, 9, 10 10, 10, 11 -*Nset, nset=Set-1, generate - 1, 11, 1 *Elset, elset=Set-1, generate 1, 10, 1 -*Nset, nset=Set-2, generate - 1, 11, 1 *Elset, elset=Set-2, generate 1, 10, 1 ** Section: Section-1 Profile: Profile-1 *Beam General Section, elset=Set-1, material=Material-1, section=GENERAL 1., 0.0833333, 0., 0.0833333, 0.140833 0.,1.,0. +*Transverse Shear Stiffness +6.73077e+10, 6.73077e+10, 0.25 *End Part ** ** @@ -48,13 +46,13 @@ ** *Assembly, name=Assembly ** -*Instance, name=Part-1-1, part=Part-1 +*Instance, name=PART-1_1-1, part=PART-1_1 *End Instance ** -*Nset, nset=Set-1, instance=Part-1-1 - 11, -*Nset, nset=Set-2, instance=Part-1-1 +*Nset, nset=Set-3, instance=PART-1_1-1 1, +*Nset, nset=Set-4, instance=PART-1_1-1 + 11, *End Assembly ** ** MATERIALS @@ -62,6 +60,17 @@ *Material, name=Material-1 *Elastic 2.1e+11, 0.3 +** +** BOUNDARY CONDITIONS +** +** Name: BC-1 Type: Displacement/Rotation +*Boundary +Set-3, 1, 1 +Set-3, 2, 2 +Set-3, 3, 3 +Set-3, 4, 4 +Set-3, 5, 5 +Set-3, 6, 6 ** ---------------------------------------------------------------- ** ** STEP: Step-1 @@ -70,22 +79,11 @@ *Static 1., 1., 1e-05, 1. ** -** BOUNDARY CONDITIONS -** -** Name: BC-1 Type: Displacement/Rotation -*Boundary -Set-2, 1, 1 -Set-2, 2, 2 -Set-2, 3, 3 -Set-2, 4, 4 -Set-2, 5, 5 -Set-2, 6, 6 -** ** LOADS ** ** Name: Load-1 Type: Concentrated force *Cload -Set-1, 3, -10000. +Set-4, 3, -1e+06 ** ** OUTPUT REQUESTS ** @@ -93,7 +91,12 @@ Set-1, 3, -10000. ** ** FIELD OUTPUT: F-Output-1 ** -*Output, field, variable=PRESELECT +*Output, field +*Node Output +CF, RF, TF, U +*Element Output, directions=YES +LE, NFORC, NFORCSO, PE, PEEQ, PEMAG, S, SF +*Contact Output, variable=PRESELECT ** ** HISTORY OUTPUT: H-Output-1 ** diff --git a/src/fesa/validation/comparison.cpp b/src/fesa/validation/comparison.cpp index 75b2fe6..be1d0bb 100644 --- a/src/fesa/validation/comparison.cpp +++ b/src/fesa/validation/comparison.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -439,4 +440,152 @@ ComparisonReport compare_samples( return report; } +CorrelationReport correlate_samples( + const std::span samples) { + struct Accumulator final { + double squared_error{}; + double squared_reference{}; + double squared_absolute_scale{}; + std::size_t value_count{}; + }; + + CorrelationReport report{true, {}, {}}; + if (samples.empty()) { + add_failure( + report.failures, + "validation.empty_comparison", + "No comparison samples were provided for correlation."); + report.evaluable = false; + return report; + } + + std::set positions; + std::map, Accumulator> + accumulators; + + for (const ComparisonSample& sample : samples) { + const PositionKey key{ + sample.quantity, + sample.position.instance_name, + sample.position.entity_label, + sample.position.end_node_label, + }; + if (!positions.insert(key).second) { + add_failure( + report.failures, + "validation.duplicate_result_position", + unevaluable_failure_text( + sample.quantity, + sample.position, + sample.tolerance, + "duplicate_result_position")); + continue; + } + + const std::size_t expected = + expected_component_count(sample.quantity); + if ( + sample.reference.size() != expected || + sample.actual.size() != expected) { + add_failure( + report.failures, + "validation.component_count_mismatch", + unevaluable_failure_text( + sample.quantity, + sample.position, + sample.tolerance, + "component_count_mismatch") + + " expected=" + std::to_string(expected) + + " reference_count=" + + std::to_string(sample.reference.size()) + + " actual_count=" + + std::to_string(sample.actual.size())); + continue; + } + + if ( + !std::isfinite(sample.tolerance.relative) || + !std::isfinite(sample.tolerance.absolute_scale) || + sample.tolerance.relative < 0.0 || + sample.tolerance.absolute_scale < 0.0) { + add_failure( + report.failures, + "validation.invalid_tolerance", + unevaluable_failure_text( + sample.quantity, + sample.position, + sample.tolerance, + "invalid_tolerance")); + continue; + } + + for (std::size_t component = 0; component < expected; ++component) { + const double reference = sample.reference[component]; + const double actual = sample.actual[component]; + if (!std::isfinite(reference) || !std::isfinite(actual)) { + add_failure( + report.failures, + "validation.nonfinite_comparison_value", + scalar_failure_text( + sample, + component, + std::numeric_limits::infinity())); + continue; + } + + const double error = actual - reference; + Accumulator& accumulator = + accumulators[{sample.quantity, component}]; + accumulator.squared_error += error * error; + accumulator.squared_reference += reference * reference; + accumulator.squared_absolute_scale += + sample.tolerance.absolute_scale * + sample.tolerance.absolute_scale; + ++accumulator.value_count; + } + } + + for (const auto& [key, accumulator] : accumulators) { + const auto [quantity, component] = key; + const double error_l2 = std::sqrt(accumulator.squared_error); + const double root_mean_square_error = + error_l2 / + std::sqrt(static_cast(accumulator.value_count)); + const double reference_l2 = + std::sqrt(accumulator.squared_reference); + const double absolute_scale_l2 = + std::sqrt(accumulator.squared_absolute_scale); + const double denominator = + std::max(reference_l2, absolute_scale_l2); + const double relative_l2_error = denominator > 0.0 + ? error_l2 / denominator + : error_l2 == 0.0 + ? 0.0 + : std::numeric_limits< + double>::infinity(); + + report.metrics.push_back({ + quantity, + component, + accumulator.value_count, + root_mean_square_error, + relative_l2_error, + }); + if ( + !std::isfinite(root_mean_square_error) || + !std::isfinite(relative_l2_error)) { + add_failure( + report.failures, + "validation.nonfinite_correlation_metric", + "quantity=" + std::string{quantity_name(quantity)} + + " component=" + component_name(quantity, component) + + " rmse=" + number_text(root_mean_square_error) + + " relative_l2=" + number_text(relative_l2_error)); + } + } + + report.evaluable = report.failures.empty(); + return report; +} + } // namespace fesa diff --git a/src/fesa/validation/reference_compare_main.cpp b/src/fesa/validation/reference_compare_main.cpp new file mode 100644 index 0000000..ab654df --- /dev/null +++ b/src/fesa/validation/reference_compare_main.cpp @@ -0,0 +1,305 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +struct ComparisonRequest final { + std::filesystem::path results; + std::string instance; + std::filesystem::path displacements; + std::filesystem::path reactions; + std::filesystem::path internal_forces; + double displacement_absolute_scale; + double reaction_absolute_scale; + double internal_force_absolute_scale; +}; + +std::string_view quantity_name(const fesa::ReferenceQuantity quantity) { + switch (quantity) { + case fesa::ReferenceQuantity::displacement: + return "displacement"; + case fesa::ReferenceQuantity::reaction: + return "reaction"; + case fesa::ReferenceQuantity::internal_force: + return "internal_force"; + case fesa::ReferenceQuantity::centroid_stress: + return "centroid_stress"; + } + return "unknown"; +} + +std::string_view component_name( + const fesa::ReferenceQuantity quantity, + const std::size_t component) { + switch (quantity) { + case fesa::ReferenceQuantity::displacement: + return fesa::NodalFrame::displacement_components[component]; + case fesa::ReferenceQuantity::reaction: + return fesa::NodalFrame::reaction_components[component]; + case fesa::ReferenceQuantity::internal_force: + return fesa::BeamElementFrame::section_force_components[component]; + case fesa::ReferenceQuantity::centroid_stress: + return fesa::BeamElementFrame::axial_stress_component; + } + return "unknown"; +} + +std::string_view stage_name(const fesa::DiagnosticStage stage) { + switch (stage) { + case fesa::DiagnosticStage::io: + return "io"; + case fesa::DiagnosticStage::syntax: + return "syntax"; + case fesa::DiagnosticStage::semantic: + return "semantic"; + case fesa::DiagnosticStage::model: + return "model"; + case fesa::DiagnosticStage::equation: + return "equation"; + case fesa::DiagnosticStage::solver: + return "solver"; + case fesa::DiagnosticStage::results: + return "results"; + case fesa::DiagnosticStage::validation: + return "validation"; + } + return "unknown"; +} + +void print_diagnostic(const fesa::Diagnostic& diagnostic) { + std::cerr << stage_name(diagnostic.stage) << " [" << diagnostic.code + << "]"; + if (diagnostic.source.has_value()) { + const fesa::SourceLocation& source = *diagnostic.source; + std::cerr << ' ' << source.file.string() << ':' << source.line << ':' + << source.column; + } + std::cerr << ": " << diagnostic.message << '\n'; +} + +void print_diagnostics(const std::vector& diagnostics) { + for (const fesa::Diagnostic& diagnostic : diagnostics) { + print_diagnostic(diagnostic); + } +} + +void print_usage() { + std::cerr + << "Usage:\n" + << " fesa-reference-compare --results " + "--instance --displacements " + "--reactions " + "--internal-forces " + "--displacement-absolute-scale " + "--reaction-absolute-scale " + "--internal-force-absolute-scale \n"; +} + +std::optional parse_number(const std::string_view text) { + double value = 0.0; + const auto parsed = std::from_chars( + text.data(), text.data() + text.size(), value); + if ( + parsed.ec != std::errc{} || + parsed.ptr != text.data() + text.size() || + !std::isfinite(value)) { + return std::nullopt; + } + return value; +} + +std::optional parse_request( + const int argc, + char* argv[]) { + if (argc != 17) { + return std::nullopt; + } + + std::optional results; + std::optional instance; + std::optional displacements; + std::optional reactions; + std::optional internal_forces; + std::optional displacement_absolute_scale; + std::optional reaction_absolute_scale; + std::optional internal_force_absolute_scale; + + for (int index = 1; index < argc; index += 2) { + const std::string_view option{argv[index]}; + const std::string_view value{argv[index + 1]}; + if (option == "--results" && !results.has_value()) { + results = std::filesystem::path{value}; + } else if (option == "--instance" && !instance.has_value()) { + instance = value; + } else if ( + option == "--displacements" && + !displacements.has_value()) { + displacements = std::filesystem::path{value}; + } else if (option == "--reactions" && !reactions.has_value()) { + reactions = std::filesystem::path{value}; + } else if ( + option == "--internal-forces" && + !internal_forces.has_value()) { + internal_forces = std::filesystem::path{value}; + } else if ( + option == "--displacement-absolute-scale" && + !displacement_absolute_scale.has_value()) { + displacement_absolute_scale = parse_number(value); + } else if ( + option == "--reaction-absolute-scale" && + !reaction_absolute_scale.has_value()) { + reaction_absolute_scale = parse_number(value); + } else if ( + option == "--internal-force-absolute-scale" && + !internal_force_absolute_scale.has_value()) { + internal_force_absolute_scale = parse_number(value); + } else { + return std::nullopt; + } + } + + if ( + !results.has_value() || results->empty() || + !instance.has_value() || instance->empty() || + !displacements.has_value() || displacements->empty() || + !reactions.has_value() || reactions->empty() || + !internal_forces.has_value() || internal_forces->empty() || + !displacement_absolute_scale.has_value() || + !reaction_absolute_scale.has_value() || + !internal_force_absolute_scale.has_value() || + *displacement_absolute_scale < 0.0 || + *reaction_absolute_scale < 0.0 || + *internal_force_absolute_scale < 0.0) { + return std::nullopt; + } + + return ComparisonRequest{ + std::move(*results), + std::move(*instance), + std::move(*displacements), + std::move(*reactions), + std::move(*internal_forces), + *displacement_absolute_scale, + *reaction_absolute_scale, + *internal_force_absolute_scale, + }; +} + +bool append_samples( + const fesa::ResultFrame& frame, + const std::vector& rows, + const fesa::Tolerance tolerance, + std::vector& samples) { + bool matched = true; + for (const fesa::ReferenceRow& row : rows) { + fesa::ComparisonSampleMatch match = fesa::make_comparison_sample( + frame, + row.quantity, + row.position, + row.values, + tolerance); + if (!match.sample.has_value()) { + print_diagnostics(match.failures); + matched = false; + continue; + } + samples.push_back(std::move(*match.sample)); + } + return matched; +} + +} // namespace + +int main(int argc, char* argv[]) { + const std::optional request = + parse_request(argc, argv); + if (!request.has_value()) { + print_usage(); + return 1; + } + + const fesa::Hdf5ReadResult results = + fesa::read_hdf5_results(request->results); + if (!results.database.has_value()) { + print_diagnostics(results.diagnostics); + return 1; + } + + const fesa::ReferenceCsvReadResult displacements = + fesa::read_reference_csv( + fesa::ReferenceQuantity::displacement, + request->displacements, + request->instance); + const fesa::ReferenceCsvReadResult reactions = + fesa::read_reference_csv( + fesa::ReferenceQuantity::reaction, + request->reactions, + request->instance); + const fesa::ReferenceCsvReadResult internal_forces = + fesa::read_reference_csv( + fesa::ReferenceQuantity::internal_force, + request->internal_forces, + request->instance); + if ( + !displacements.diagnostics.empty() || + !reactions.diagnostics.empty() || + !internal_forces.diagnostics.empty()) { + print_diagnostics(displacements.diagnostics); + print_diagnostics(reactions.diagnostics); + print_diagnostics(internal_forces.diagnostics); + return 1; + } + + const fesa::ResultFrame& frame = + results.database->steps.front().frames.front(); + std::vector samples; + samples.reserve( + displacements.rows.size() + reactions.rows.size() + + internal_forces.rows.size()); + const bool displacements_matched = append_samples( + frame, + displacements.rows, + {0.0, request->displacement_absolute_scale}, + samples); + const bool reactions_matched = append_samples( + frame, + reactions.rows, + {0.0, request->reaction_absolute_scale}, + samples); + const bool internal_forces_matched = append_samples( + frame, + internal_forces.rows, + {0.0, request->internal_force_absolute_scale}, + samples); + if ( + !displacements_matched || !reactions_matched || + !internal_forces_matched) { + return 1; + } + + const fesa::CorrelationReport report = fesa::correlate_samples(samples); + std::cout << std::setprecision(std::numeric_limits::max_digits10); + for (const fesa::ComponentCorrelationMetric& metric : report.metrics) { + std::cout << "metric quantity=" << quantity_name(metric.quantity) + << " component=" + << component_name(metric.quantity, metric.component_index) + << " count=" << metric.value_count + << " rmse=" << metric.root_mean_square_error + << " relative_l2=" << metric.relative_l2_error << '\n'; + } + print_diagnostics(report.failures); + return report.evaluable ? 0 : 1; +} diff --git a/src/fesa/validation/reference_csv.cpp b/src/fesa/validation/reference_csv.cpp index a17e0ae..e42e453 100644 --- a/src/fesa/validation/reference_csv.cpp +++ b/src/fesa/validation/reference_csv.cpp @@ -64,6 +64,12 @@ std::vector split_fields(const std::string_view line) { return fields; } +void remove_trailing_empty_fields(std::vector& fields) { + while (!fields.empty() && fields.back().empty()) { + fields.pop_back(); + } +} + ReferenceCsvReadResult failure( const DiagnosticStage stage, std::string code, @@ -196,7 +202,7 @@ std::vector value_columns( case ReferenceQuantity::reaction: return {"RF-RF1", "RF-RF2", "RF-RF3", "RM-RM1", "RM-RM2", "RM-RM3"}; case ReferenceQuantity::internal_force: - return {"SF-SF1", "SF-SF2", "SF-SF3", "SM-SM1", "SM-SM2", "SM-SM3"}; + return {"SF-SF1", "SF-SF3", "SF-SF2", "SM-SM3", "SM-SM1", "SM-SM2"}; case ReferenceQuantity::centroid_stress: return {"Sxx"}; } @@ -283,7 +289,8 @@ ReferenceCsvReadResult read_reference_csv( 1U); } - const std::vector headers = split_fields(line); + std::vector headers = split_fields(line); + remove_trailing_empty_fields(headers); ColumnIndices column_indices; for (std::size_t index = 0U; index < headers.size(); ++index) { if (headers[index].empty() || @@ -340,7 +347,8 @@ ReferenceCsvReadResult read_reference_csv( if (trim(line).empty()) { continue; } - const std::vector fields = split_fields(line); + std::vector fields = split_fields(line); + remove_trailing_empty_fields(fields); if (fields.size() != headers.size()) { return validation_failure( "validation.reference_csv_invalid_row", diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 56f1971..cb10e38 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -773,6 +773,12 @@ add_test( --gtest_filter=ComparisonMetric.* ) +add_test( + NAME CorrelationMetric + COMMAND "$" + --gtest_filter=CorrelationMetric.* +) + add_test( NAME EntityMatching COMMAND "$" @@ -796,3 +802,5 @@ add_test( COMMAND "$" --gtest_filter=StressCsv.* ) + +add_subdirectory(reference) diff --git a/tests/integration/io/minimal_deck_to_domain_test.cpp b/tests/integration/io/minimal_deck_to_domain_test.cpp index 2e8016a..9569fbf 100644 --- a/tests/integration/io/minimal_deck_to_domain_test.cpp +++ b/tests/integration/io/minimal_deck_to_domain_test.cpp @@ -702,7 +702,7 @@ TEST(Cload, ReportsInvalidDofAtTheDataRow) { TEST(SuppliedCantilever, NormalizesReferenceModelThroughPublicParserAndMapper) { const std::filesystem::path path = std::filesystem::path{FESA_TEST_SOURCE_DIR}.parent_path() / - "reference" / "cantilever beam" / "cantilever beam.inp"; + "reference" / "cantilever beam" / "cantilever beam fesa.inp"; const auto result = parse_and_map(path); @@ -723,7 +723,7 @@ TEST(SuppliedCantilever, NormalizesReferenceModelThroughPublicParserAndMapper) { EXPECT_EQ(domain.step().nodal_loads.front().node, node->id); EXPECT_EQ( domain.step().nodal_loads.front().values, - (std::array{0.0, 0.0, -10000.0, 0.0, 0.0, 0.0})); + (std::array{0.0, 0.0, -1.0e6, 0.0, 0.0, 0.0})); } } // namespace diff --git a/tests/reference/CMakeLists.txt b/tests/reference/CMakeLists.txt new file mode 100644 index 0000000..24a3d04 --- /dev/null +++ b/tests/reference/CMakeLists.txt @@ -0,0 +1,38 @@ +add_executable(fesa_cantilever_reference_tests + cantilever_reference_test.cpp +) + +target_compile_features( + fesa_cantilever_reference_tests PRIVATE cxx_std_20 +) +target_compile_options( + fesa_cantilever_reference_tests PRIVATE /W4 /permissive- /EHsc +) +target_compile_definitions( + fesa_cantilever_reference_tests + PRIVATE + FESA_REFERENCE_COMPARE_PATH="$" + FESA_REPOSITORY_ROOT="${CMAKE_SOURCE_DIR}" + FESA_TEST_BINARY_DIR="${CMAKE_BINARY_DIR}" +) +target_link_libraries( + fesa_cantilever_reference_tests + PRIVATE + fesa_core + GTest::gtest_main +) +add_dependencies( + fesa_cantilever_reference_tests + fesa-reference-compare +) + +add_test( + NAME CantileverReference + COMMAND "$" +) + +set_property( + TEST CantileverReference + PROPERTY ENVIRONMENT_MODIFICATION + ${FESA_DEPENDENCY_RUNTIME_MODIFICATIONS} +) diff --git a/tests/reference/cantilever_reference_test.cpp b/tests/reference/cantilever_reference_test.cpp new file mode 100644 index 0000000..18875ff --- /dev/null +++ b/tests/reference/cantilever_reference_test.cpp @@ -0,0 +1,241 @@ +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +constexpr double kEquilibriumRelativeTolerance = 2.0e-13; +constexpr double kEquilibriumAbsoluteTolerance = 1.0e-12; +constexpr std::string_view kInstanceName = "PART-1_1-1"; + +class TemporaryPath final { +public: + explicit TemporaryPath(std::filesystem::path path) + : path_{std::move(path)} { + std::error_code error; + std::filesystem::create_directories(path_.parent_path(), error); + if (error) { + throw std::runtime_error{ + "Failed to create reference test directory."}; + } + std::filesystem::remove(path_, error); + } + + ~TemporaryPath() { + std::error_code error; + std::filesystem::remove(path_, error); + } + + TemporaryPath(const TemporaryPath&) = delete; + TemporaryPath& operator=(const TemporaryPath&) = delete; + + [[nodiscard]] const std::filesystem::path& path() const noexcept { + return path_; + } + +private: + std::filesystem::path path_; +}; + +std::filesystem::path reference_path(const std::string_view name) { + return std::filesystem::path{FESA_REPOSITORY_ROOT} / "reference" / + "cantilever beam" / name; +} + +std::filesystem::path test_output_path(const std::string_view name) { + return std::filesystem::path{FESA_TEST_BINARY_DIR} / "testing" / name; +} + +std::string quote(const std::filesystem::path& path) { + return '"' + path.string() + '"'; +} + +std::string read_text(const std::filesystem::path& path) { + std::ifstream input{path, std::ios::binary}; + return { + std::istreambuf_iterator{input}, + std::istreambuf_iterator{}, + }; +} + +const fesa::Vec3& coordinates_for( + const fesa::Hdf5ModelSnapshot& model, + const fesa::NodeId node_id) { + const auto found = std::ranges::find_if( + model.nodes, + [node_id](const fesa::Hdf5NodeSnapshot& node) { + return node.id == node_id; + }); + if (found == model.nodes.end()) { + throw std::runtime_error{ + "Result or load references an unknown node ID."}; + } + return found->coordinates; +} + +void add_nodal_resultant( + std::array& force, + std::array& moment, + const fesa::Vec3& position, + const std::array& values) { + force[0] += values[0]; + force[1] += values[1]; + force[2] += values[2]; + moment[0] += values[3] + position.y * values[2] - + position.z * values[1]; + moment[1] += values[4] + position.z * values[0] - + position.x * values[2]; + moment[2] += values[5] + position.x * values[1] - + position.y * values[0]; +} + +void expect_finite(const fesa::ResultFrame& frame) { + for (const auto& displacement : frame.nodal.displacement) { + for (const double value : displacement) { + EXPECT_TRUE(std::isfinite(value)); + } + } + for (const auto& reaction : frame.nodal.reaction) { + for (const double value : reaction) { + EXPECT_TRUE(std::isfinite(value)); + } + } + for (const fesa::BeamElementFrame& beam : frame.element.beams) { + EXPECT_TRUE(fesa::is_finite(beam.local_frame.ex)); + EXPECT_TRUE(fesa::is_finite(beam.local_frame.ey)); + EXPECT_TRUE(fesa::is_finite(beam.local_frame.ez)); + for (const fesa::BeamSectionResult& end : beam.end_results) { + EXPECT_TRUE(std::isfinite(end.xi)); + EXPECT_TRUE(std::isfinite(end.centroid_sigma_xx)); + for (const double value : end.section_strain) { + EXPECT_TRUE(std::isfinite(value)); + } + for (const double value : end.section_force) { + EXPECT_TRUE(std::isfinite(value)); + } + for (const double value : end.sigma_xx) { + EXPECT_TRUE(std::isfinite(value)); + } + } + } +} + +TEST(CantileverReference, CorrelatesAllAvailableAbaqusResults) { + const TemporaryPath results{ + test_output_path("cantilever-reference.h5")}; + const TemporaryPath standard_output{ + test_output_path("cantilever-reference-compare.stdout.txt")}; + const TemporaryPath error_output{ + test_output_path("cantilever-reference-compare.stderr.txt")}; + + const fesa::AnalysisRunResult run = fesa::run_solver({ + reference_path("cantilever beam fesa.inp"), + results.path(), + }); + ASSERT_TRUE(run.succeeded); + ASSERT_TRUE(run.diagnostics.empty()); + + const fesa::Hdf5ReadResult read = + fesa::read_hdf5_results(results.path()); + ASSERT_TRUE(read.diagnostics.empty()); + ASSERT_TRUE(read.database.has_value()); + ASSERT_TRUE(read.model.has_value()); + ASSERT_TRUE(read.analysis.has_value()); + ASSERT_EQ(read.database->steps.size(), 1U); + ASSERT_EQ(read.database->steps[0].frames.size(), 1U); + + const fesa::ResultFrame& frame = + read.database->steps[0].frames[0]; + expect_finite(frame); + + std::array total_force{}; + std::array total_moment{}; + std::array applied_force{}; + std::array applied_moment{}; + ASSERT_EQ(frame.nodal.node_ids.size(), frame.nodal.reaction.size()); + for (std::size_t index = 0; index < frame.nodal.node_ids.size(); ++index) { + add_nodal_resultant( + total_force, + total_moment, + coordinates_for(*read.model, frame.nodal.node_ids[index]), + frame.nodal.reaction[index]); + } + for (const fesa::NodalLoad& load : read.analysis->step.nodal_loads) { + const fesa::Vec3& position = + coordinates_for(*read.model, load.node); + add_nodal_resultant( + total_force, + total_moment, + position, + load.values); + add_nodal_resultant( + applied_force, + applied_moment, + position, + load.values); + } + for (std::size_t component = 0; component < total_force.size(); ++component) { + const double tolerance = std::max( + kEquilibriumAbsoluteTolerance, + kEquilibriumRelativeTolerance * + std::abs(applied_force[component])); + EXPECT_NEAR(total_force[component], 0.0, tolerance); + } + for (std::size_t component = 0; component < total_moment.size(); ++component) { + const double tolerance = std::max( + kEquilibriumAbsoluteTolerance, + kEquilibriumRelativeTolerance * + std::abs(applied_moment[component])); + EXPECT_NEAR(total_moment[component], 0.0, tolerance); + } + + const std::string command = + '"' + quote(std::filesystem::path{FESA_REFERENCE_COMPARE_PATH}) + + " --results " + quote(results.path()) + + " --instance " + std::string{kInstanceName} + + " --displacements " + + quote(reference_path("cantilever beam displacements.csv")) + + " --reactions " + + quote(reference_path("cantilever beam reactions.csv")) + + " --internal-forces " + + quote(reference_path("cantilever beam elemental forces.csv")) + + std::string{" --displacement-absolute-scale 1e-10"} + + " --reaction-absolute-scale 1e-8" + + " --internal-force-absolute-scale 1e-8" + + " 1>" + quote(standard_output.path()) + + " 2>" + quote(error_output.path()) + '"'; + const int exit_code = std::system(command.c_str()); + + const std::string standard_text = read_text(standard_output.path()); + const std::string error_text = read_text(error_output.path()); + + EXPECT_EQ(exit_code, 0) + << standard_text << error_text; + EXPECT_NE( + standard_text.find("quantity=displacement"), + std::string::npos); + EXPECT_NE( + standard_text.find("quantity=reaction"), + std::string::npos); + EXPECT_NE( + standard_text.find("quantity=internal_force"), + std::string::npos); + EXPECT_NE(standard_text.find("rmse="), std::string::npos); + EXPECT_NE(standard_text.find("relative_l2="), std::string::npos); +} + +} // namespace diff --git a/tests/unit/validation/comparison_test.cpp b/tests/unit/validation/comparison_test.cpp index b1c7820..0b19854 100644 --- a/tests/unit/validation/comparison_test.cpp +++ b/tests/unit/validation/comparison_test.cpp @@ -240,6 +240,66 @@ TEST(ComparisonMetric, RejectsComponentCountMismatch) { report.failures, "validation.component_count_mismatch")); } +TEST(CorrelationMetric, ComputesComponentWiseRmseAndRelativeL2) { + const std::vector inputs{ + sample( + fesa::ReferenceQuantity::displacement, + {"Part-1-1", 101, std::nullopt}, + {3.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 1.0e-9}), + sample( + fesa::ReferenceQuantity::displacement, + {"Part-1-1", 102, std::nullopt}, + {4.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 1.0e-9}), + }; + + const auto report = fesa::correlate_samples(inputs); + + ASSERT_TRUE(report.evaluable); + ASSERT_TRUE(report.failures.empty()); + ASSERT_EQ(report.metrics.size(), 6U); + EXPECT_EQ(report.metrics[0].quantity, fesa::ReferenceQuantity::displacement); + EXPECT_EQ(report.metrics[0].component_index, 0U); + EXPECT_EQ(report.metrics[0].value_count, 2U); + EXPECT_DOUBLE_EQ( + report.metrics[0].root_mean_square_error, + std::sqrt(12.5)); + EXPECT_DOUBLE_EQ(report.metrics[0].relative_l2_error, 1.0); + EXPECT_DOUBLE_EQ(report.metrics[1].root_mean_square_error, 0.0); + EXPECT_DOUBLE_EQ(report.metrics[1].relative_l2_error, 0.0); +} + +TEST(CorrelationMetric, UsesAbsoluteScaleNormForNearZeroReference) { + const std::vector inputs{ + sample( + fesa::ReferenceQuantity::reaction, + {"Part-1-1", 101, std::nullopt}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {3.0e-9, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 1.0e-9}), + sample( + fesa::ReferenceQuantity::reaction, + {"Part-1-1", 102, std::nullopt}, + {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, + {4.0e-9, 0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 1.0e-9}), + }; + + const auto report = fesa::correlate_samples(inputs); + + ASSERT_TRUE(report.evaluable); + ASSERT_EQ(report.metrics.size(), 6U); + EXPECT_DOUBLE_EQ( + report.metrics[0].root_mean_square_error, + std::sqrt(12.5) * 1.0e-9); + EXPECT_DOUBLE_EQ( + report.metrics[0].relative_l2_error, + 5.0 / std::sqrt(2.0)); +} + TEST(EntityMatching, MatchesNodalResultByInstanceAndExternalLabel) { const auto frame = result_frame(); const std::array reference{0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; diff --git a/tests/unit/validation/reference_csv_test.cpp b/tests/unit/validation/reference_csv_test.cpp index 6b896cb..7d6ce61 100644 --- a/tests/unit/validation/reference_csv_test.cpp +++ b/tests/unit/validation/reference_csv_test.cpp @@ -72,14 +72,14 @@ TEST(ReferenceCsv, ReadsSuppliedDisplacementsWithWhitespaceHeader) { EXPECT_EQ( result.rows.front().quantity, fesa::ReferenceQuantity::displacement); - EXPECT_EQ(result.rows.front().position.instance_name, "Part-1-1"); + EXPECT_EQ(result.rows.front().position.instance_name, "PART-1_1-1"); EXPECT_EQ(result.rows.front().position.entity_label, 1); EXPECT_FALSE(result.rows.front().position.end_node_label.has_value()); EXPECT_EQ( result.rows.front().values, - (std::vector{0.0, 0.0, -1.0e-32, 0.0, 1.0e-31, 0.0})); + (std::vector{0.0, 0.0, -1.0e-30, 0.0, 1.0e-29, 0.0})); EXPECT_EQ(result.rows.back().position.entity_label, 11); - EXPECT_EQ(result.rows.back().values[2], -1.92e-4); + EXPECT_EQ(result.rows.back().values[2], -1.91857e-2); } TEST(ReferenceCsv, ReadsSuppliedReactionsWithWhitespaceHeader) { @@ -90,12 +90,12 @@ TEST(ReferenceCsv, ReadsSuppliedReactionsWithWhitespaceHeader) { ASSERT_TRUE(result.diagnostics.empty()); ASSERT_EQ(result.rows.size(), 11U); - EXPECT_EQ(result.rows.front().position.instance_name, "Part-1-1"); + EXPECT_EQ(result.rows.front().position.instance_name, "PART-1_1-1"); EXPECT_EQ(result.rows.front().position.entity_label, 1); EXPECT_FALSE(result.rows.front().position.end_node_label.has_value()); EXPECT_EQ( result.rows.front().values, - (std::vector{0.0, 0.0, 1.0e4, 0.0, -1.0e5, 0.0})); + (std::vector{0.0, 0.0, 1.0e6, 0.0, -1.0e7, 0.0})); } TEST(InternalForceCsv, MapsAllSixComponentsAndElementEndPosition) { @@ -115,7 +115,26 @@ TEST(InternalForceCsv, MapsAllSixComponentsAndElementEndPosition) { EXPECT_EQ(*result.rows.front().position.end_node_label, 101); EXPECT_EQ( result.rows.front().values, - (std::vector{1.25, -2.5, 3.75, -4.0, 5.5, -6.25})); + (std::vector{1.25, 3.75, -2.5, -6.25, -4.0, 5.5})); +} + +TEST(InternalForceCsv, AcceptsTrailingEmptyAbaqusExportColumns) { + const TemporaryCsv input{ + "fesa-reference-trailing-empty-columns.csv", + "Part Instance Name,Element Label,Node Label,SF-SF1,SF-SF2," + "SF-SF3,SM-SM1,SM-SM2,SM-SM3,,,\n" + "PART-1_1-1,1,1,1,2,3,4,5,6,,,\n"}; + + const auto result = fesa::read_reference_csv( + fesa::ReferenceQuantity::internal_force, + input.path(), + "unused-request-name"); + + ASSERT_TRUE(result.diagnostics.empty()); + ASSERT_EQ(result.rows.size(), 1U); + EXPECT_EQ( + result.rows.front().values, + (std::vector{1.0, 3.0, 2.0, 6.0, 4.0, 5.0})); } TEST(StressCsv, FillsOmittedInstanceAndReadsCentroidStress) {