Compare commits

...

73 Commits

Author SHA1 Message Date
KOKO\Mimi 5855604318 docs(handoff): prepare internal release transition 2026-08-03 02:19:08 +09:00
KOKO\Mimi 6bc7cc3ada chore(beam-reference-qualification): step 3 output 2026-08-03 01:48:33 +09:00
KOKO\Mimi 8c448ffe45 docs(beam-reference-qualification): record dual-gate evidence 2026-08-03 01:48:10 +09:00
KOKO\Mimi ca0268ea5b chore(beam-reference-qualification): step 2 output 2026-08-03 01:46:00 +09:00
KOKO\Mimi 675379779f feat(beam-reference-qualification): correlate Abaqus beam results 2026-08-03 01:45:30 +09:00
KOKO\Mimi b68f6ee143 docs(beam-reference-qualification): define dual validation gates 2026-08-03 01:08:54 +09:00
KOKO\Mimi 02e2994c5b chore(beam-reference-qualification): step 1 output 2026-08-02 03:33:06 +09:00
KOKO\Mimi 4f96719663 feat(beam-reference-qualification): step 1 — reference-csv-adapters 2026-08-02 03:33:06 +09:00
KOKO\Mimi 2df90b95f8 chore(beam-reference-qualification): step 0 output 2026-08-02 03:22:08 +09:00
KOKO\Mimi 98c13c2ec2 feat(beam-reference-qualification): step 0 — comparison-metric-and-entity-matching 2026-08-02 03:22:08 +09:00
KOKO\Mimi 51939fba5b fix(result-contract-completion): enforce complete result provenance 2026-08-02 02:34:07 +09:00
KOKO\Mimi f23deb0ade chore(result-contract-completion): mark phase completed 2026-08-02 01:51:25 +09:00
KOKO\Mimi fb67007527 chore(result-contract-completion): step 2 output 2026-08-02 01:51:24 +09:00
KOKO\Mimi 4d04f3dbbe feat(result-contract-completion): step 2 — self-contained-hdf5 2026-08-02 01:51:24 +09:00
KOKO\Mimi 9d37465d92 chore(result-contract-completion): step 1 output 2026-08-02 01:28:37 +09:00
KOKO\Mimi e0a6a6fa70 feat(result-contract-completion): step 1 — complete-result-contract 2026-08-02 01:28:37 +09:00
KOKO\Mimi 218cfa9d50 chore(result-contract-completion): step 0 output 2026-08-02 01:14:43 +09:00
KOKO\Mimi 5fe57cb145 feat(result-contract-completion): step 0 — beam-element-end-recovery 2026-08-02 01:14:40 +09:00
KOKO\Mimi 5618636f0d docs(handoff): prepare result contract phase 2026-08-02 00:20:53 +09:00
KOKO\Mimi 7d52247d6c fix(assembly): stage benchmark TBB runtime 2026-08-01 23:42:47 +09:00
KOKO\Mimi ac8d24a2cf test(assembly): cover deterministic kernel failures 2026-08-01 23:40:56 +09:00
KOKO\Mimi 6ad9f3cd47 chore(deterministic-parallel-assembly): mark phase completed 2026-08-01 23:36:18 +09:00
KOKO\Mimi 6ee51a8502 chore(deterministic-parallel-assembly): step 2 output 2026-08-01 23:36:18 +09:00
KOKO\Mimi 127286cf2b feat(deterministic-parallel-assembly): step 2 — thread-count-determinism 2026-08-01 23:36:17 +09:00
KOKO\Mimi dc6baed8ed chore(deterministic-parallel-assembly): step 1 output 2026-08-01 23:21:27 +09:00
KOKO\Mimi 6c2e1f3ab1 feat(deterministic-parallel-assembly): step 1 — tbb-element-evaluation 2026-08-01 23:21:27 +09:00
KOKO\Mimi b9bb439766 chore(deterministic-parallel-assembly): step 0 output 2026-08-01 23:06:22 +09:00
KOKO\Mimi 40a7e6c3af feat(deterministic-parallel-assembly): step 0 — canonical-contribution-order 2026-08-01 23:06:21 +09:00
KOKO\Mimi 6ac474f19b docs(handoff): prepare parallel assembly phase 2026-08-01 20:53:18 +09:00
KOKO\Mimi cbb621bb28 fix(input): validate inactive Part records 2026-08-01 04:07:18 +09:00
KOKO\Mimi af886f3f90 fix(input): enforce strict Abaqus subset contract 2026-08-01 03:48:29 +09:00
KOKO\Mimi 9269847c83 chore(abaqus-subset-completion): mark phase completed 2026-08-01 03:23:25 +09:00
KOKO\Mimi eec4fe1f49 chore(abaqus-subset-completion): step 4 output 2026-08-01 03:23:25 +09:00
KOKO\Mimi 94c83a39c5 feat(abaqus-subset-completion): step 4 — step-bc-load-and-noop-directives 2026-08-01 03:23:25 +09:00
KOKO\Mimi af076c39a3 chore(abaqus-subset-completion): step 3 output 2026-08-01 02:59:40 +09:00
KOKO\Mimi a1fed69c47 feat(abaqus-subset-completion): step 3 — material-section-and-shear-defaults 2026-08-01 02:59:40 +09:00
KOKO\Mimi 8cd2d5b2ae chore(abaqus-subset-completion): step 2 output 2026-08-01 02:44:01 +09:00
KOKO\Mimi e6bc708be3 feat(abaqus-subset-completion): step 2 — single-instance-semantic-validation 2026-08-01 02:44:01 +09:00
KOKO\Mimi 3eeab2fbe4 chore(abaqus-subset-completion): step 1 output 2026-08-01 02:37:05 +09:00
KOKO\Mimi 3588aa2bb6 feat(abaqus-subset-completion): step 1 — part-and-assembly-set-resolution 2026-08-01 02:37:05 +09:00
KOKO\Mimi 91041f28c9 chore(abaqus-subset-completion): step 0 output 2026-08-01 02:24:24 +09:00
KOKO\Mimi 2b83b1128f feat(abaqus-subset-completion): step 0 — abaqus-input-contract 2026-08-01 02:24:24 +09:00
KOKO\Mimi 18361e4eb2 fix(results-and-pipeline): set CLI test runtime path 2026-08-01 01:26:06 +09:00
KOKO\Mimi e31ef999bd chore(results-and-pipeline): mark phase completed 2026-08-01 01:11:01 +09:00
KOKO\Mimi 93496c0157 chore(results-and-pipeline): step 3 output 2026-08-01 01:11:01 +09:00
KOKO\Mimi fe90528115 feat(results-and-pipeline): step 3 — cli-pipeline-integration 2026-08-01 01:11:01 +09:00
KOKO\Mimi 493219d3c2 chore(results-and-pipeline): step 2 output 2026-08-01 01:01:15 +09:00
KOKO\Mimi 11672df2c4 feat(results-and-pipeline): step 2 — linear-static-analysis 2026-08-01 01:01:15 +09:00
KOKO\Mimi d6306383a9 chore(results-and-pipeline): step 1 output 2026-08-01 00:51:34 +09:00
KOKO\Mimi 76a0c8ecde feat(results-and-pipeline): step 1 — minimal-hdf5-schema 2026-08-01 00:51:34 +09:00
KOKO\Mimi f8af2c1485 chore(results-and-pipeline): step 0 output 2026-08-01 00:14:54 +09:00
KOKO\Mimi 5b7de91988 feat(results-and-pipeline): step 0 — result-database 2026-08-01 00:14:53 +09:00
KOKO\Mimi 4b75b72968 modify handoff.md 2026-07-31 23:54:01 +09:00
KOKO\Mimi a02024929c docs(equation-and-linear-solve): align completed contracts 2026-07-31 16:37:34 +09:00
KOKO\Mimi 741fc9eaea fix(equation-and-linear-solve): address review findings 2026-07-31 16:36:07 +09:00
KOKO\Mimi 2dd17be5d0 chore(equation-and-linear-solve): mark phase completed 2026-07-31 16:16:28 +09:00
KOKO\Mimi 38096007f7 chore(equation-and-linear-solve): step 2 output 2026-07-31 16:16:28 +09:00
KOKO\Mimi e216660d51 feat(equation-and-linear-solve): step 2 — pardiso-linear-solver 2026-07-31 16:16:28 +09:00
KOKO\Mimi a78f36a172 chore(equation-and-linear-solve): step 1 output 2026-07-31 16:00:45 +09:00
KOKO\Mimi 1bf277cbf1 feat(equation-and-linear-solve): step 1 — essential-bc-elimination 2026-07-31 16:00:45 +09:00
KOKO\Mimi 30dcb05dd1 chore(equation-and-linear-solve): step 0 output 2026-07-31 15:45:03 +09:00
KOKO\Mimi 8342774c82 feat(equation-and-linear-solve): step 0 — symmetric-csr-assembly 2026-07-31 15:45:02 +09:00
KOKO\Mimi d393402174 modify handoff.md 2026-07-31 15:12:11 +09:00
KOKO\Mimi fa48a4d9d7 docs(fem-and-beam-kernel): align integration description 2026-07-31 02:29:09 +09:00
KOKO\Mimi 334f9e79ae chore(fem-and-beam-kernel): mark phase completed 2026-07-31 02:24:48 +09:00
KOKO\Mimi 7af3fba0d4 chore(fem-and-beam-kernel): step 3 output 2026-07-31 02:24:48 +09:00
KOKO\Mimi 93b2c3d43e feat(fem-and-beam-kernel): step 3 — timoshenko-stiffness-kernel 2026-07-31 02:24:48 +09:00
KOKO\Mimi 07f0bb3c13 chore(fem-and-beam-kernel): step 2 output 2026-07-31 02:08:27 +09:00
KOKO\Mimi 7ddd8e690e feat(fem-and-beam-kernel): step 2 — beam-local-frame 2026-07-31 02:08:27 +09:00
KOKO\Mimi 63a71b7e03 chore(fem-and-beam-kernel): step 1 output 2026-07-31 01:55:21 +09:00
KOKO\Mimi dc6b9f1bb4 feat(fem-and-beam-kernel): step 1 — dof-manager 2026-07-31 01:55:21 +09:00
KOKO\Mimi 7f28047d06 chore(fem-and-beam-kernel): step 0 output 2026-07-31 01:48:33 +09:00
KOKO\Mimi 80ea9759da feat(fem-and-beam-kernel): step 0 — quadrature-and-shape-functions 2026-07-31 01:48:33 +09:00
184 changed files with 19740 additions and 715 deletions
+1 -1
View File
@@ -60,7 +60,7 @@
- 전단강성이 생략되면 \(A_{sy}=A_{sz}=5A/6\)과 `SCF=0`을 Phase 1 기본값으로
적용한다.
- reference 비교는 metadata 없이 요청된 물리량과 CSV 경로를 명시한다. 현재
캔틸레버 샘플은 변위 반력을 비교하며 요소 내력과 도심 응력 비교 루틴은
캔틸레버 샘플은 변위, 반력 및 요소 단면력을 비교하며 도심 응력 비교 루틴은
synthetic CSV로 검증한다.
- 새 MSVC 빌드 경고를 추가하지 않는다.
- 변경은 요청 범위에 한정하고 Conventional Commits 형식의 메시지를 사용한다.
+46
View File
@@ -26,11 +26,29 @@ include(CTest)
include(cmake/FesaDependencies.cmake)
add_library(fesa_core STATIC
src/fesa/analysis/linear_static_analysis.cpp
src/fesa/analysis/run_solver.cpp
src/fesa/assembly/contribution.cpp
src/fesa/assembly/parallel_assembler.cpp
src/fesa/assembly/serial_assembler.cpp
src/fesa/constraints/essential_bc.cpp
src/fesa/core/version.cpp
src/fesa/elements/beam/beam3d2.cpp
src/fesa/fem/beam_frame.cpp
src/fesa/fem/dof_manager.cpp
src/fesa/fem/gauss_rule.cpp
src/fesa/fem/line2_shape.cpp
src/fesa/io/abaqus/active_input.cpp
src/fesa/io/abaqus/parser.cpp
src/fesa/io/abaqus/semantic_mapper.cpp
src/fesa/io/abaqus/set_resolver.cpp
src/fesa/io/hdf5/writer.cpp
src/fesa/model/domain.cpp
src/fesa/model/domain_builder.cpp
src/fesa/results/result_database.cpp
src/fesa/solvers/linear/pardiso_linear_solver.cpp
src/fesa/validation/comparison.cpp
src/fesa/validation/reference_csv.cpp
)
target_include_directories(fesa_core
@@ -40,6 +58,7 @@ target_include_directories(fesa_core
target_compile_features(fesa_core PUBLIC cxx_std_20)
target_compile_options(fesa_core PRIVATE /W4 /permissive- /EHsc)
target_link_libraries(fesa_core PRIVATE MKL::MKL TBB::tbb HDF5::HDF5)
add_executable(fesa
src/fesa/cli/main.cpp
@@ -49,6 +68,33 @@ 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
)
target_link_libraries(fesa_assembly_benchmark PRIVATE fesa_core)
target_compile_features(fesa_assembly_benchmark PRIVATE cxx_std_20)
target_compile_options(
fesa_assembly_benchmark PRIVATE /W4 /permissive- /EHsc
)
add_custom_command(
TARGET fesa_assembly_benchmark
POST_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_if_different
"$<TARGET_FILE:TBB::tbb>"
"$<TARGET_FILE_DIR:fesa_assembly_benchmark>"
)
add_subdirectory(tests)
endif()
+298
View File
@@ -0,0 +1,298 @@
# FESA Phase 1 Abaqus Input Subset
## 1. Status and scope
This document is the normative input contract for the FESA Phase 1 Abaqus
adapter. FESA accepts only the keywords, parameters, scopes, and data forms
defined here. It does not implement general Abaqus input syntax, and it never
silently ignores an unknown keyword or option.
Two mutually exclusive organizations are supported:
- a flat/orphan mesh whose mesh and set records are in global scope; or
- one or more `*PART` definitions followed by exactly one `*ASSEMBLY` that
contains exactly one untransformed `*INSTANCE` of the active Part.
Only the active Part is normalized into `Domain`. Names and external labels are
input identity; generated nonnegative FESA IDs and dense indices are separate.
FESA performs no unit conversion.
`tests/fixtures/abaqus/contract.tsv` is the executable valid/invalid fixture
matrix for this document. Its expected diagnostic code and line are part of
the contract. A diagnostic for a bad data row points to that row; a keyword,
parameter, or scope error points to the keyword row.
## 2. Common lexical rules
- Input is UTF-8. A UTF-8 BOM is accepted only at the start of the file.
- Blank lines and lines whose first non-whitespace characters are `**` are
ignored. A comment or blank line does not end the current keyword record.
- Keyword names, parameter names, flag parameters, and the enumerated values
named below are ASCII case-insensitive. Entity names are trimmed and then
matched exactly.
- Fields are comma-separated and surrounding ASCII whitespace is ignored.
- External node and element labels are positive signed 64-bit integers. DOF
numbers are integers 1 through 6. Real fields must be finite.
- A flag parameter has no `=` value. A valued parameter must have one nonempty
value. Duplicate parameters are `abaqus.syntax.duplicate_parameter`.
- Parameters not listed for a keyword are
`abaqus.syntax.unsupported_parameter`. Unknown keywords, including
`*INCLUDE`, are `abaqus.unsupported_keyword`.
- A non-comment data line without an open data-bearing keyword is
`abaqus.syntax.data_without_keyword`.
## 3. Scope and ordering
The parser maintains `global`, `part`, `assembly`, `instance`, and `step`
scope. `*STEP` is a global child scope: model-data records already opened in
global scope remain model data, while `*STATIC`, `*CLOAD`, `*RESTART`, and
`*OUTPUT` belong to the open Step.
The accepted ordering is:
```text
optional *HEADING and *PREPRINT
flat mesh, or one or more *PART blocks and one *ASSEMBLY block
global materials
optional global model-data *BOUNDARY
one *STEP
one *STATIC
optional *BOUNDARY and *CLOAD
optional no-op *RESTART and *OUTPUT
*END STEP
```
References are resolved after the complete deck is parsed. A material may
therefore follow the Part that uses it, and a nested set may refer to a set
declared later in the same scope.
## 4. Mesh and hierarchy keywords
### 4.1 `*NODE`
- Scope: flat global or Part; not Assembly, Instance, or Step.
- Parameters: none.
- Data: one or more `label, x, y, z` rows, with exactly four nonempty fields.
- Semantics: labels are unique in their mesh scope. Reuse in an inactive Part
is allowed because it is a different Part scope.
- Diagnostics: `abaqus.syntax.invalid_node_scope`,
`abaqus.semantic.invalid_node_data`, `abaqus.semantic.duplicate_node_label`.
### 4.2 `*ELEMENT`
- Scope: flat global or Part.
- Parameters: required `TYPE=B31`; optional `ELSET=<name>`; no others.
- Data: one or more `element_label, node_1_label, node_2_label` rows.
- Semantics: element labels are unique in their mesh scope. Both nodes must
exist in that scope. `ELSET=` adds every row to the named element set and
merges with an explicit set of the same name using sorted-unique membership.
- Diagnostics: `abaqus.syntax.invalid_element_scope`,
`abaqus.semantic.unsupported_element`,
`abaqus.semantic.invalid_element_data`,
`abaqus.semantic.duplicate_element_label`,
`abaqus.semantic.missing_node`.
### 4.3 Part delimiters
`*PART` is global-only, requires exactly `NAME=<name>`, and accepts no data.
Part names are unique. `*END PART` accepts no parameters or data and closes the
open Part. Nesting or a mismatched delimiter is invalid.
Diagnostics are `abaqus.syntax.invalid_part_scope`,
`abaqus.syntax.unexpected_end_part`, `abaqus.syntax.unclosed_part`, and
`abaqus.semantic.duplicate_part`.
### 4.4 Assembly delimiters
`*ASSEMBLY` is global-only, requires exactly `NAME=<name>`, and accepts no
data. Phase 1 accepts exactly one Assembly. `*END ASSEMBLY` has no parameters
or data and closes the open Assembly.
Diagnostics are `abaqus.syntax.invalid_assembly_scope`,
`abaqus.syntax.multiple_assemblies`,
`abaqus.syntax.unexpected_end_assembly`, and
`abaqus.syntax.unclosed_assembly`.
### 4.5 Instance delimiters
`*INSTANCE` is Assembly-only and requires exactly `NAME=<name>, PART=<name>`.
Phase 1 accepts exactly one Instance. No translation or rotation data and no
keyword record are allowed inside it. `*END INSTANCE` has no parameters or
data.
Diagnostics are `abaqus.syntax.invalid_instance_scope`,
`abaqus.syntax.instance_local_keyword`,
`abaqus.syntax.unexpected_end_instance`,
`abaqus.syntax.unclosed_instance`, `abaqus.semantic.instance_count`,
`abaqus.semantic.instance_transform`, and `abaqus.semantic.missing_part`.
A transform diagnostic points to the first transform data row.
Flat `*NODE`/`*ELEMENT` records combined with any Part/Assembly organization
are `abaqus.semantic.mixed_mesh_organization`.
## 5. Sets
`*NSET` and `*ELSET` are allowed in flat global, Part, or Assembly scope.
- Required parameter: respectively `NSET=<name>` or `ELSET=<name>`.
- Optional flag: `GENERATE`.
- Assembly scope additionally requires `INSTANCE=<active-instance-name>`.
`INSTANCE=` is forbidden in flat global and Part scope.
- Explicit data consists of comma-separated positive labels and/or names of
sets of the same kind and scope. Empty trailing fields are ignored.
- `GENERATE` data consists of exactly one `start, end, increment` row. All
values are positive integers, `start <= end`, and `(end-start)` is divisible
by `increment`.
- Repeated declarations of the same set merge. Nested references are resolved
independent of declaration order, cycles are rejected, and final membership
is deterministic sorted-unique.
- An Assembly set lifts active-Part local labels through its named Instance.
It cannot reference an inactive or unknown Instance.
Diagnostics are `abaqus.syntax.invalid_set_scope`,
`abaqus.semantic.invalid_generate`, `abaqus.semantic.set_cycle`,
`abaqus.semantic.missing_set_member`, and
`abaqus.semantic.wrong_instance`.
## 6. Material and Beam section
### 6.1 `*MATERIAL` and `*ELASTIC`
`*MATERIAL` is global-only, requires exactly `NAME=<name>`, and has no data.
Material names are unique. Its `*ELASTIC` child is global model data, has no
parameters, and has exactly one `young_modulus, poisson_ratio` row. Young's
modulus is finite and positive; Poisson's ratio is finite and satisfies
`-1 < nu < 0.5`. Temperature and field dependencies are not supported.
Diagnostics are `abaqus.syntax.invalid_material_scope`,
`abaqus.semantic.duplicate_material`,
`abaqus.semantic.elastic_without_material`,
`abaqus.semantic.missing_elastic`, and
`abaqus.semantic.invalid_elastic_data`.
### 6.2 `*BEAM GENERAL SECTION`
- Scope: flat global or Part.
- Parameters: required `SECTION=GENERAL`, `ELSET=<name>`, and
`MATERIAL=<name>`; no others.
- First data row: exactly `A, I_y, I_yz, I_z, J`. `A`, `I_y`, `I_z`, and `J`
are finite and positive; `I_yz` must be finite and exactly zero for Phase 1.
- Second data row: exactly three finite components of the local section-axis
reference direction. Model validation rejects a zero direction or one
parallel to an assigned element axis.
- The material and set may be declared later, but must resolve. Every active
B31 element has exactly one section assignment.
Diagnostics are `abaqus.syntax.invalid_section_scope`,
`abaqus.semantic.unsupported_section`,
`abaqus.semantic.invalid_section_data`,
`abaqus.semantic.missing_material`,
`abaqus.semantic.missing_element_set`,
`abaqus.semantic.missing_section`, and
`abaqus.semantic.duplicate_section_assignment`.
### 6.3 `*TRANSVERSE SHEAR STIFFNESS`
This optional record is in the same flat-global or Part scope as, and must
immediately follow, the affected `*BEAM GENERAL SECTION`. It has no parameters
and exactly one `K23, K13, SCF` data row. `K23` and `K13` are finite and
positive. Phase 1 accepts only numeric `SCF=0`; omitted/default `0.25`, nonzero
values, and the Abaqus `SCF` label are unsupported.
For isotropic `G=E/[2(1+nu)]`, FESA stores
`A_sy=K23/G` and `A_sz=K13/G` with source `input`. If this keyword is absent,
the semantic mapper stores `A_sy=A_sz=5A/6`, `SCF=0`, with source
`phase1_default`.
The data order follows the Abaqus 2024
[*TRANSVERSE SHEAR STIFFNESS* reference](https://docs.software.vt.edu/abaqusv2024/English/SIMACAEKEYRefMap/simakey-r-transverseshearstiffness.htm);
the restriction to numeric zero SCF and the effective-area mapping are FESA
Phase 1 decisions.
Diagnostics are `abaqus.syntax.invalid_transverse_shear_scope`,
`abaqus.semantic.orphan_transverse_shear`,
`abaqus.semantic.invalid_transverse_shear_data`, and
`abaqus.semantic.nonzero_scf`.
## 7. Linear static step, BC, and load
### 7.1 `*BOUNDARY`
`*BOUNDARY` is allowed as global model data before the Step or inside the sole
Step. It has no parameters. Each row is
`node-or-nset, first_dof[, last_dof[, value]]`. `last_dof` defaults to
`first_dof`; value defaults to zero. The inclusive DOF range is 1 through 6.
Repeated identical prescriptions are deduplicated; differing values for one
node/DOF are rejected.
Diagnostics are `abaqus.syntax.invalid_boundary_scope`,
`abaqus.semantic.invalid_boundary_data`,
`abaqus.semantic.invalid_dof`, `abaqus.semantic.invalid_dof_range`,
`abaqus.semantic.missing_node_target`, and
`abaqus.semantic.conflicting_boundary`.
### 7.2 `*CLOAD`
`*CLOAD` is Step-only and has no parameters. Each row is exactly
`node-or-nset, dof, magnitude`; DOF is 1 through 6 and magnitude is finite.
Loads on the same node/DOF are summed in input order after target resolution.
Diagnostics are `abaqus.syntax.invalid_cload_scope`,
`abaqus.semantic.invalid_cload_data`, `abaqus.semantic.invalid_dof`, and
`abaqus.semantic.missing_node_target`.
### 7.3 `*STEP`, `*STATIC`, and `*END STEP`
`*STEP` is global-only. Optional parameters are `NAME=<name>` and
`NLGEOM=NO`; omitted name becomes `Step-1`. Exactly one Step is required.
`NLGEOM=YES` and every other option are unsupported.
Exactly one `*STATIC` occurs inside the Step. It has no parameters and accepts
either no data row or one row of one through four finite positive values
`initial_increment[, time_period[, minimum_increment[, maximum_increment]]]`.
The values are accepted as load-step metadata; Phase 1 performs one linear
solve.
`*END STEP` has no parameters or data and closes the Step.
Diagnostics are `abaqus.syntax.invalid_step_scope`,
`abaqus.syntax.unexpected_end_step`, `abaqus.syntax.unclosed_step`,
`abaqus.semantic.step_count`, `abaqus.semantic.unsupported_step_option`,
`abaqus.semantic.missing_static`, and
`abaqus.semantic.invalid_static_data`.
## 8. Recognized no-op directives
These records are deliberately recognized and do not create Domain entities:
- `*HEADING`: global-only, no parameters, zero or more text data rows.
- `*PREPRINT`: global-only, optional `ECHO`, `MODEL`, `HISTORY`, and `CONTACT`
parameters, each with value `YES` or `NO`; no data.
- `*RESTART`: Step-only, optional `WRITE` flag and optional nonnegative integer
`FREQUENCY`; no data.
- `*OUTPUT`: Step-only, exactly one `FIELD` or `HISTORY` flag and optional
`VARIABLE=PRESELECT`; no data.
Diagnostics are the common unsupported-parameter diagnostic plus
`abaqus.syntax.invalid_heading_scope`,
`abaqus.syntax.invalid_preprint_scope`,
`abaqus.syntax.invalid_restart_scope`, and
`abaqus.syntax.invalid_output_scope`. There is no general no-op or
ignore-unknown path.
## 9. Fixture matrix contract
The tab-separated manifest columns are:
```text
case_id, outcome, fixture, expected_stage, expected_code, expected_line,
node_count, element_count, node_set_count, element_set_count,
prescribed_dof_count, nodal_load_count, shear_source, shear_area_y,
shear_area_z, checked_node_set, checked_element_set
```
For a valid case, the public `parse_deck()` and `map_deck_to_domain()` path must
produce a Domain matching every populated expected field. Set checks use
`name:local-label,...`. For an invalid case, that same public path must fail
and contain the exact stage, code, and source line in the manifest. Partial
decks and test-only semantic construction are not accepted.
+43 -2
View File
@@ -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: 일관 단위계와 결과 좌표계
@@ -257,7 +259,7 @@ flat `Domain`으로 정규화한다. 외부 entity는 `(instance name, part-loca
## ADR-015: 명시적 물리량 선택 기반 CSV 검증
**상태:** Accepted
**상태:** Superseded by ADR-017
**상황:** 현재 캔틸레버 reference에는 변위와 반력만 있고 per-model metadata는
요구하지 않는다. 요소 내력과 응력 비교 기능은 해당 CSV가 추가되기 전에 구현해야
@@ -297,3 +299,42 @@ toolset을 명시한다. CMake, CMake Presets, CTest 및 GoogleTest/GoogleMock
보장 대상이 아니다.
- 컴파일러 갱신에 따른 경고와 표준 라이브러리 동작은 전체 Debug/Release 검증에서
다시 확인해야 한다.
## ADR-017: FESA 정식화 적합성과 Abaqus 결과 상관성의 이중 gate
**상태:** Accepted
**상황:** FESA의 2절점 Timoshenko Beam은 선택적 감차적분과 `SCF=0`을 사용하고
Abaqus B31의 slenderness compensation을 구현하지 않는다. 따라서 `SCF=0.25`
Abaqus 결과에 기본 상대오차 \(10^{-5}\)를 적용하면 FESA 정식화 결함과 의도된
정식화 차이를 구분할 수 없다. 현재 캔틸레버 reference에는 변위, 반력 및 요소
단면력 CSV가 있다.
**결정:**
- FESA 정식화 적합성 gate는 해석해, 에너지, 강체 mode, 평형 및 엄격한 tolerance로
FESA 자체 정식화를 검증한다.
- Abaqus 결과 상관성 gate는 원본 Abaqus 입력과 CSV를 변경하지 않는다. FESA는
동일한 기하·재료·하중과 명시적 전단강성을 사용하되 `SCF=0`인 별도 입력을
production parser와 solver로 해석한다.
- 상관성 gate는 요청된 모든 entity와 component가 유일하게 매칭되고 유한한
component별 RMSE 및 Relative L2가 생성되면 evaluable이다. 서로 다른 정식화에
기본 상대오차 \(10^{-5}\) pass/fail을 적용하지 않는다.
- Relative L2의 reference norm이 영에 가까우면 해당 component의 characteristic
absolute scale norm을 분모 하한으로 사용한다.
- 병진과 회전, 힘과 모멘트처럼 단위가 다른 component를 하나의 RMSE 또는 norm에
혼합하지 않는다.
- Abaqus의 \((\mathbf t,\mathbf n_1,\mathbf n_2)\)와 FESA의
\((\mathbf e_x,\mathbf e_y,\mathbf e_z)\)를 각각 일치시킨 모델에서 Abaqus
`SF1,SF2,SF3,SM1,SM2,SM3`은 FESA \(N,V_y,V_z,T,M_y,M_z\) 순서로
`SF1,SF3,SF2,SM3,SM1,SM2`를 사용한다.
- 현재 캔틸레버 상관성 요청은 변위, 반력 및 요소 단면력을 명시한다. 요청하지 않은
응력은 통과로 보고하지 않는다.
**결과와 트레이드오프:**
- FESA 구현 회귀와 상용 solver와의 모델 상관성을 서로 오인하지 않는다.
- Abaqus 원본과 FESA 투영 입력을 함께 관리해야 하며 formulation 차이를
`docs/VALIDATION.md`에 기록해야 한다.
- 첫 상관성 보고서는 metric을 제시하지만 관측값에 맞춘 acceptance envelope를
만들지 않는다. 후속 envelope에는 해석적 또는 mesh study 근거가 필요하다.
+7 -6
View File
@@ -429,16 +429,17 @@ kernel을 추가한다.
- `tests/reference`: CSV 골든 결과와 FESA HDF5 결과 비교
- `reference/<model-id>`: 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를 직접 주입해 전체 파이프라인 결함을 숨기지
+270 -223
View File
@@ -1,82 +1,234 @@
# FESA Session Handoff
## 1. 문서 목적
## 1. 목적과 기준 문서
이 문서는 `domain-and-input-skeleton` 완료 후 새 세션에서
`fem-and-beam-kernel` Phase를 바로 시작하기 위한 인수인계 기록이다.
요구사항과 설계의 기준은 이 문서가 아니라 다음 파일이다.
이 문서는 `beam-reference-qualification` 완료 후 새 세션에서 마지막 Phase 1 단계인
`internal-release` 시작하기 위한 인수인계 기록이다. 과거 phase의 구현 역사를
반복하기보다 현재 기준선, 검증 결과, 남은 작업과 실행 순서를 제공한다.
- `/AGENTS.md`
- `/docs/PRD.md`
- `/docs/ARCHITECTURE.md`
- `/docs/ADR.md`
- `/docs/HARNESS.md`
- `/docs/superpowers/plans/2026-07-29-fesa-phase-1.md`
- `/phases/fem-and-beam-kernel/index.json`
- `/phases/fem-and-beam-kernel/step0.md`부터 `step3.md`
다음 문서를 우선 기준으로 사용한다.
내용이 충돌하면 위 기준 문서와 `phases/`의 현재 상태를 우선한다.
- `AGENTS.md`
- `docs/PRD.md`
- `docs/ARCHITECTURE.md`
- `docs/ADR.md`
- `docs/HARNESS.md`
- `docs/HDF5_SCHEMA.md`
- `docs/ABAQUS_INPUT_SUBSET.md`
- `docs/formulation/timoshenko-beam-3d.md`
- `docs/VALIDATION.md`
- `phases/internal-release/index.json`
- `phases/internal-release/step0.md`부터 `step4.md`
## 2. 현재 저장소 상태
내용이 충돌하면 `AGENTS.md`, PRD/ADR/아키텍처, 완료된 phase metadata와 실제
테스트를 우선한다. 특히 `internal-release`의 기존 Step 0과 Step 4에 남아 있는
reference 범위 설명은 4절의 최신 계약으로 바로잡은 뒤 release evidence에 사용한다.
2026-07-31 확인 기준:
## 2. Git과 Phase 상태
- 현재 브랜치: `dev`
- `dev`, `origin/dev`, `origin/HEAD` 기준 커밋:
`e7f0d53406bfdff0fbd610c2b595a8d214341353`
- 완료 Phase:
- `solver-bootstrap`
- `domain-and-input-skeleton`
- 다음 Phase: `fem-and-beam-kernel`
- 다음 Step: `0 - quadrature-and-shape-functions`
- `fem-and-beam-kernel`의 Step 0~3은 모두 `pending`
- 활성 제품 코드 blocker는 없다.
- `include/fesa/fem/`, `include/fesa/elements/beam/`,
`docs/formulation/timoshenko-beam-3d.md`는 아직 존재하지 않는다. 해당 Step에서
실제 구현과 함께 만든다.
2026-08-03 기준 구현 상태는 다음과 같다.
이 문서를 갱신하기 직전 작업 트리는 clean이었다. 문서 수정 자체가 아직 커밋되지 않은
경우 Harness 실행 전에 `docs/HANDOFF.md` 변경을 먼저 커밋하거나 별도로 정리해야 한다.
- 기준 브랜치: `dev`
- 이 문서 갱신 직전 구현 HEAD: `6bc7cc3adadded6f1a2a0ed45449b2264c0f7a4e`
- `feat-beam-reference-qualification``dev`에 fast-forward 병합된 뒤 로컬에서
삭제됐다.
- `internal-release`를 제외한 `phases/index.json`의 모든 phase가 `completed`다.
- 다음 Phase: `internal-release`
- 다음 Step: Step 0 `release-checklist`
- Step 0부터 Step 4까지 모두 `pending`이다.
- 이 문서의 갱신은 구현 기준선 다음의 docs-only commit이다.
## 3. 완료된 기반
새 세션에서는 기록된 hash를 강제로 맞추지 말고 로컬과 원격의 실제 상태를 먼저
확인한다. 이 HANDOFF commit이 push된 뒤에는 `dev``origin/dev`가 같은 commit을
가리켜야 한다.
### `solver-bootstrap`
```powershell
git switch dev
git status --short --branch
git rev-parse HEAD
git rev-parse origin/dev
git rev-list --left-right --count origin/dev...dev
```
- Visual Studio 2026 MSVC v145, Windows x64, C++20 CMake Preset
- `fesa_core``fesa` CLI 골격
- MKL, TBB, HDF5, GoogleTest package discovery와 dependency smoke test
- typed entity ID, `Vec3`, source location, diagnostic, status 값 타입
- Harness Python characterization test와 CMake/CTest 실행 계약
reset, rebase 또는 force push로 차이를 숨기지 않는다. 예상하지 못한 변경이 있으면
소유자를 확인하고 보존한다.
### `domain-and-input-skeleton`
## 3. 완료된 Phase 1 기준선
- solver semantic 값 타입:
- `Node`, `BeamElement`, `IsotropicElastic`, `BeamSection`
- `NodeSet`, `ElementSet`, `StepDefinition`
- `EntityOrigin`과 typed internal ID
- build 이후 const view만 제공하는 `Domain`
- `DomainBuilder`의 중복 ID/origin, 참조, 유한값, 물성, 길이, orientation,
경계조건 검증
- flat/orphan mesh와 Part/Assembly/단일 무변환 Instance scope를 보존하는 Abaqus
syntax parser
- 활성 Part만 선택해 flat `Domain`으로 정규화하는 semantic mapper
- 실제 Part/Instance/local label provenance 보존
- 생략된 전단강성에
\(A_{sy}=A_{sz}=5A/6\)과
`ShearPropertySource::phase1_default` 적용
- 여러 Instance, Instance transform, missing Part 등 최소 조직 오류의 source
diagnostic
- public parser/mapper 경로를 사용하는 flat/계층형 통합 테스트
현재 production 경로는 다음 기능을 제공한다.
현재 parser/mapper는 다음 Phase의 Beam kernel 테스트에 필요한 최소 semantic
`Domain`을 제공한다. 전체 Abaqus 부분집합, nested/`GENERATE` set, 완전한 scope
resolution, explicit transverse shear, no-op directive와 단일 step의 모든 option은
후속 `abaqus-subset-completion` Phase 책임이다. 이번 Phase에서 선행 구현하지 않는다.
- Abaqus `.inp` 제한 부분집합의 flat/orphan mesh 또는 좌표변환이 없는 단일
Part/Assembly/Instance 모델 파싱과 semantic validation
- 2절점 3D Isoparametric Timoshenko Beam, 등방성 선형 탄성, 일반 단면,
`*BOUNDARY`, `*CLOAD`, 단일 선형 정적 Step
- MSVC v145, Intel oneAPI MKL PARDISO와 TBB를 사용하는 Windows x64 해석 경로
- canonical contribution ordering에 기반한 결정론적 병렬 요소 평가와 조립
- 원래 full equation에서 변위, 반력과 평형을 복구하는 선형 정적 해석
- 두 요소 끝의 section strain, section force, centroid/recovery-point `Sxx` 회복
- node/element provenance, component와 좌표계 metadata를 포함한 완전한
`ResultDatabase`
- 입력 원본 없이 model, analysis와 결과를 재구성할 수 있는 HDF5 schema `2.0.0`
- 변위, 반력, 요소 단면력용 Abaqus CSV adapter와 component별 상관성 보고 CLI
## 4. 검증된 개발환경과 baseline
결과 계약과 Beam reference phase의 세부 완료 기록은 다음 파일에 있다.
새 PowerShell 세션에서 configure 또는 Harness 실행 전에 다음 환경 변수를 설정한다.
절대경로를 tracked CMake 파일이나 Preset에 넣지 않는다.
- `phases/result-contract-completion/index.json`
- `phases/beam-reference-qualification/index.json`
- `docs/HDF5_SCHEMA.md`
- `docs/VALIDATION.md`
`internal-release`에서는 위 solver 계약을 다시 설계하지 않는다. 설치, clean consumer,
scale measurement와 release evidence에 필요한 최소 변경만 수행한다.
## 4. 검증 상태와 남은 제한
### 4.1 이중 검증 gate
`docs/VALIDATION.md`가 검증 결과의 단일 상세 보고서다.
- Gate A — FESA 정식화 적합성: **PASS**
- Gate B — Abaqus 결과 상관성: **EVALUABLE**
Gate A는 해석해, strain energy, 강체 mode, 회전 불변성, 평형과 결정성을 엄격한
tolerance로 검증한다. Gate B는 Abaqus B31과 FESA 정식화의 차이를 인정하고
component별 RMSE와 Relative L2를 보고한다. Gate B의 `EVALUABLE`은 모든 요청
entity/component가 유일하게 매칭되고 metric이 유한하다는 뜻이며, 임의의 상관성
pass/fail threshold를 통과했다는 뜻은 아니다.
현재 cantilever reference의 주요 Relative L2는 다음과 같다.
| 결과 | Component | Relative L2 |
|---|---|---:|
| displacement | `Uz` | `0.0022349969291714038` |
| displacement | `Ry` | `0.0000010416581667076092` |
| reaction | `RFz` | `4.4393520322089023e-13` |
| reaction | `RMy` | `2.591919334788851e-13` |
| internal force | `Vz` | `2.7107472798172426e-13` |
| internal force | `My` | `0.082541022764834646` |
18개 전체 component의 count, RMSE와 Relative L2는 `docs/VALIDATION.md`를 사용한다.
### 4.2 Reference 모델과 축 매핑
- Abaqus provenance는 `reference/cantilever beam/cantilever beam.inp`와 세 CSV다.
- Abaqus 원본 모델의 `SCF=0.25`는 보존한다.
- FESA projection은 `reference/cantilever beam/cantilever beam fesa.inp`이며 다른
의미 입력을 유지하고 `SCF=0`을 사용한다.
- FESA에 Abaqus SCF 보정이나 결과 맞춤 계수를 추가하지 않는다.
요소력 CSV는 다음 순서로 매핑한다.
| FESA | Abaqus CSV |
|---|---|
| `N` | `SF1` |
| `Vy` | `SF3` |
| `Vz` | `SF2` |
| `T` | `SM3` |
| `My` | `SM1` |
| `Mz` | `SM2` |
### 4.3 아직 완료되지 않은 검증
- 실제 Abaqus 결과로 검증된 물리량은 변위, 반력과 요소 단면력이다.
- Abaqus 도심 응력 CSV는 제공되지 않았다. `StressCsv`는 synthetic fixture로 schema와
entity mapping만 검증했으며 응력은 **not yet Abaqus-qualified**다.
- synthetic internal-force fixture도 adapter 단위 검증에 계속 사용하지만, 요소
단면력 자체는 별도의 실제 Abaqus golden CSV와 상관성 비교가 완료됐다.
- 내부 배포용 Release build, install tree, clean consumer smoke test와 약 100,000 DOF
측정은 아직 실행되지 않았다. 증거 없이 완료 표시하지 않는다.
### 4.4 `internal-release` Step 문서의 계약 불일치
`phases/internal-release/step0.md``step4.md`에는 beam reference phase 이전의 문구가
남아 있다.
- “현재 Abaqus displacement/reaction comparison”은 변위·반력·요소 단면력
correlation으로 갱신해야 한다.
- “내력·응력은 synthetic coverage”라는 묶음 표현은 요소 단면력과 응력을 분리해야
한다. 요소 단면력은 real golden correlation과 synthetic adapter coverage가 모두
있고, 응력만 synthetic adapter coverage다.
새 세션은 Harness 실행 전에 이 두 Step 문서와 관련 checklist 문구를 현재
`docs/PRD.md` 8절 및 `docs/VALIDATION.md`와 일치시켜야 한다. 이 정렬은 검증 범위의
확장이 아니라 이미 완료된 증거를 정확히 기술하는 작업이다.
## 5. 다음 Phase: `internal-release`
Phase metadata는 `phases/internal-release/index.json`에 있다. Step은 순서대로 실행한다.
### Step 0 — `release-checklist`
- `docs/BUILDING.md`, `docs/INPUT_FORMAT.md`, `docs/RELEASE_CHECKLIST.md`를 작성한다.
- PRD 8절의 각 내부 배포 기준에 고유 checklist ID를 부여한다.
- 실제 target, preset, example과 evidence command를 CMake에서 재확인한다.
- 미실행 Release/install/benchmark 항목을 완료 표시하지 않는다.
- 4.4절의 stale reference 문구를 먼저 바로잡는다.
### Step 1 — `cmake-install-package`
- `cmake --install`로 내부 배포용 install tree를 만든다.
- CLI, `fesa_core`, public headers, CMake package config, example, schema/input/validation
문서와 runtime DLL inventory를 포함한다.
- install config에 absolute build path가 남는 실패 검사를 먼저 작성한다.
- installer, registry write, 외부 dependency 다운로드와 public ABI 약속은 범위 밖이다.
### Step 2 — `install-tree-smoke-test`
- source/build tree를 참조하지 않는 consumer configure/link test를 만든다.
- 설치된 CLI의 `--version`, example solve와 생성 HDF5 open을 검증한다.
- 개발 PATH의 `fesa.exe`나 source include fallback으로 결함을 숨기지 않는다.
### Step 3 — `phase1-scale-benchmark`
- 약 100,000 DOF Beam chain에서 generation/parsing, assembly, PARDISO solve,
recovery와 HDF5 write 시간을 분리해 측정한다.
- 재현 가능한 Windows memory metric, thread/solver 설정과 실제 DOF 수를 기록한다.
- finite result, 평형과 정상 종료는 검사하되 임의 시간·speedup 기준은 만들지 않는다.
### Step 4 — `release-evidence-gate`
- Debug와 Release configure/build/test를 새로 실행한다.
- test count가 0이 아닌지 확인한다.
- Harness pytest, reference, determinism, HDF5 inspection, install-tree smoke test와 scale
benchmark 증거를 checklist ID에 연결한다.
- 모든 수용 조건에 현재 실행 증거가 있을 때만 Step과 Phase를 `completed`로 바꾼다.
- 실패나 미실행 항목이 있으면 정확한 blocker를 남기고 release 완료를 선언하지 않는다.
## 6. 유지해야 할 경계
- C++20, Visual Studio 2026 MSVC v145와 Windows x64 기준을 유지한다.
- `core`, `model`, `fem`, `elements`에 Abaqus, MKL, TBB 또는 HDF5 API를 노출하지
않는다.
- solver semantic model과 result contract에 installer 또는 serialization 전용 타입을
추가하지 않는다.
- dependency는 개발 환경에 사전 설치된 버전을 사용하며 package 중 다운로드하지
않는다.
- install tree는 source/build tree의 절대경로에 의존하지 않아야 한다.
- Debug와 Release artifact/runtime을 혼합하지 않는다.
- FESA는 단위 변환을 수행하지 않는다.
- performance 수치를 correctness gate로 바꾸지 않는다.
- stress를 Abaqus-qualified로 표현하지 않는다.
- 테스트를 disable하거나 제외해 release evidence를 만들지 않는다.
- 사용자가 명시적으로 요청하지 않은 phase 실행에서는 `--push`를 사용하지 않는다.
## 7. 검증 기준선과 개발환경
### 7.1 마지막 확인 결과
2026-08-03의 beam reference qualification 완료 및 `dev` 병합 후 다음을 확인했다.
- `cmake --build --preset windows-debug`: 성공, 새 MSVC warning 없음
- `ctest --preset windows-debug --output-on-failure`: 68/68 통과
- `uv run --with pytest python -m pytest -v -rs`: 20/20 통과
- thread `{1,2,16}`에서 10회 반복한 조립/해석 결과: bitwise 동일
이 수치는 새 세션이 유지해야 할 Debug baseline이다. Release와 install-tree 결과로
확대 해석하지 않는다.
### 7.2 Package 설정
새 PowerShell 세션에서 configure 전에 현재 설치 위치를 확인한다. tracked preset이나
CMake 파일에 사용자별 절대경로를 넣지 않는다.
```powershell
$env:MKL_DIR = "C:\Program Files (x86)\Intel\oneAPI\2026.1\lib\cmake\mkl"
@@ -90,188 +242,83 @@ Test-Path "$env:HDF5_DIR\hdf5-config.cmake"
Test-Path "$env:GTest_DIR\GTestConfig.cmake"
```
2026-07-31 위 네 경로가 모두 존재함을 확인했다. 같은 환경에서 다음 baseline을
새로 검증했다.
- MSVC v145 Debug build 성공, 새 경고 없음
- CTest 12개 중 12개 성공
- Harness pytest 20개 중 20개 성공
- pytest가 실제로 20개를 수집했으므로 0-test 성공이 아님
검증 명령:
네 경로가 모두 유효한 같은 에서 preset을 실행한다. cache가 없거나 package 위치가
변경됐을 때만 `--fresh` configure를 사용한다.
```powershell
cmake --fresh --preset windows-debug
cmake --build --preset windows-debug
ctest --preset windows-debug --output-on-failure
uv run --with pytest python -m pytest -v -rs
```
CMake cache가 없거나 package 경로가 바뀐 경우에만 같은 환경 변수 세션에서 먼저
다음을 실행한다.
Release evidence는 `internal-release` Step에서 새로 생성한다.
```powershell
cmake --fresh --preset windows-debug
cmake --preset windows-release
cmake --build --preset windows-release
ctest --preset windows-release --output-on-failure
```
## 5. 다음 Phase 목표와 Step 순서
## 8. Harness 주의사항
`fem-and-beam-kernel`의 독립 deliverable은 실제 2절점 3D Timoshenko Beam의
local/global stiffness와 해석적 sanity test다. 네 Step을 순서대로 실행한다.
### Step 0 — `quadrature-and-shape-functions`
- `fem`에 1점/2점 1D Gauss rule, 2절점 선형 shape function과 derivative,
`length/2` Jacobian을 구현한다.
- partition of unity, endpoint interpolation, derivative sum zero, 적분 정확도,
invalid order/length를 실패 테스트로 먼저 고정한다.
- runtime registry나 Beam stiffness를 만들지 않는다.
### Step 1 — `dof-manager`
- 절점당 자유도 순서는
\(u_x,u_y,u_z,r_x,r_y,r_z\)의 6개다.
- `DofManager`가 full DOF, constrained/free equation numbering,
12개 요소 DOF mapping과 prescribed value를 포함한 full-vector reconstruction을
전담한다.
- external label이나 입력 선언 순서에 흔들리지 않는 deterministic numbering을
테스트로 고정한다.
- equation ID를 `Node``BeamElement`에 저장하지 않는다.
- sparse pattern, MPC, penalty DOF는 이 Step 범위가 아니다.
### Step 2 — `beam-local-frame`
- 두 절점과 mandatory orientation vector로 오른손 직교
`BeamFrame {ex, ey, ez}`를 만든다.
- Gram-Schmidt, 정규화, 직교성, determinant \(+1\), 회전/에너지 invariant를
테스트한다.
- zero length, zero orientation, 요소축과 평행한 orientation은 diagnostic으로
실패시킨다.
- `Matrix12`는 현재 코드에 없으므로 이 Step의 transformation 계약에 필요한 최소
고정 크기 타입만 만든다. 범용 동적 matrix 계층이나 registry를 만들지 않는다.
### Step 3 — `timoshenko-stiffness-kernel`
- production code보다 먼저
`/docs/formulation/timoshenko-beam-3d.md`를 작성한다.
- 문서에는 최소한 다음을 고정한다.
- 12개 local/global DOF와 component 순서
- 오른손 국부 좌표계 정의
- 변형률과 부호 규약
- \(G=E/[2(1+\nu)]\)와 constitutive matrix
- 자연좌표, shape derivative와 \(J=L/2\)
- 축·굽힘·비틀림 2점, 전단 1점 선택적 감차적분
- local/global transformation 규칙
- 실제 `Beam3D2Input`, `Beam3D2Contribution`, `compute_beam3d2` kernel을 구현한다.
- 대칭성, 강체운동 zero energy, 축/비틀림/단축·이축 굽힘 해석해,
shear-dominant 문제, 세장비 sweep과 회전 invariant를 실패 테스트로 먼저 만든다.
- 임시 가짜/닫힌형 stiffness로 파이프라인만 통과시키지 않는다.
## 6. 수치·정식화 계약
- 요소: 2절점 직선 3D Isoparametric Timoshenko Beam
- 절점당 6 DOF:
\(u_x,u_y,u_z,\theta_x,\theta_y,\theta_z\)
- 선형 shape function, \(\xi\in[-1,1]\)
- 축·굽힘·비틀림: 2점 Gauss 적분
- 전단: 1점 Gauss 적분
- 재료: 등방성 선형 탄성,
\(G=E/[2(1+\nu)]\)
- 단면:
\(A,I_y,I_z,J,A_{sy},A_{sz}\)
- 제한:
\(I_{yz}=0\), 도심/전단중심 일치, offset/warping 없음
- orientation은 필수이며 자동 추측하거나 임의 축으로 대체하지 않는다.
- tolerance는 테스트 이름이나 주석에 수학적/scale-aware 근거를 기록한다.
- FESA는 단위 변환을 하지 않는다.
정식화는 `docs/PRD.md` 3.2절, `docs/ARCHITECTURE.md` 8절,
`docs/ADR.md` ADR-005와 ADR-014를 우선한다.
## 7. 아키텍처 경계와 이번 Phase 제외 범위
- public header는 `include/fesa/`, 구현은 `src/fesa/`, 테스트는 `tests/`에 둔다.
- `core`, `model`, `fem`, `elements`는 Abaqus, MKL, TBB, HDF5 API에 의존하지
않는다.
- `fem`은 DOF, equation mapping, quadrature, shape function, Jacobian,
local/global mapping만 담당한다.
- `elements/beam`은 Beam local contribution과 kernel 계약을 담당한다.
- `Domain`은 읽기 전용 view로 사용하고 복제하거나 mutable accessor를 추가하지 않는다.
- 이번 Phase에서는 다음을 구현하지 않는다.
- sparse COO/CSR pattern과 조립
- essential BC elimination system
- MKL PARDISO solve
- oneTBB parallel assembly
- HDF5 writer와 result recovery
- 여러 요소 타입을 위한 registry/factory
- 비선형, warping, offset, \(I_{yz}\ne0\), 여러 Instance
## 8. 새 세션 시작 절차
먼저 이 문서와 Step 파일을 읽고 현재 상태를 재검증한다.
표준 실행 명령은 다음과 같다.
```powershell
git switch dev
git status --short --branch
git rev-parse HEAD
git rev-parse origin/dev
python scripts/execute.py internal-release
```
작업 트리가 clean이고 `dev`가 의도한 기준인지 확인한 뒤 4절의 package 환경 변수를
설정하고 baseline을 실행한다. 그 다음:
이전 beam reference phase에서 Harness child가 WindowsApps PowerShell을 시작하지 못해
`CreateProcessAsUserW` 오류가 반복됐다. 동일 저장소와 preset에서 직접 실행한 수용
명령은 성공했고, `stepN-output.json`에는 child 환경 blocker와 직접 실행 증거를
구분해 기록했다.
새 세션에서는 다음 원칙을 지킨다.
1. 먼저 표준 Harness 실행을 시도한다.
2. 같은 child shell 오류가 발생하면 중복 executor를 시작하지 말고 남은 process를
확인한다.
3. 직접 실행한 명령과 Harness 실행 결과를 혼동하지 않고 metadata에 사실대로 적는다.
4. 사용자 profile이나 drive root를 `--codex-add-dir`로 허용하지 않는다.
5. global Codex 설정을 임시 변경했다면 원래 값을 기록하고 종료 즉시 복원·재확인한다.
6. 사용자가 push를 명시하지 않은 phase 실행에는 `--push`를 추가하지 않는다.
## 9. 새 세션 시작 순서
1. `AGENTS.md`와 이 문서를 읽는다.
2. `docs/PRD.md` 8절, `docs/VALIDATION.md`, `phases/internal-release/index.json`
Step 0~4를 읽는다.
3. Git 상태와 `dev == origin/dev`를 확인한다.
4. 4.4절의 Step 0/4 reference 범위 문구를 최신 계약에 맞춘다.
5. Debug baseline을 새로 실행한다.
6. Step 0의 release 문서와 checklist 요구조건을 먼저 테스트 가능한 형태로 고정한다.
7. 다음 명령으로 Phase를 실행한다.
```powershell
python scripts/execute.py fem-and-beam-kernel
python scripts/execute.py internal-release
```
executor는 `feat-fem-and-beam-kernel` 브랜치를 생성하거나 checkout하고 Step 상태와
output metadata를 기록한다. 사용자가 명시적으로 요청하지 않은 한 `--push`
사용하지 않는다.
각 Step에서는 실패하는 검사 또는 미충족 evidence를 먼저 확인하고, 최소 변경으로
수용 조건을 만족시킨 뒤 focused test와 전체 test를 실행한다. Step output과 phase
metadata는 실제 명령 결과를 그대로 반영한다.
각 Step은 다음 순서를 지킨다.
## 10. `internal-release` 완료 조건
1. Step 파일의 필수 문서와 선행 구현을 모두 읽는다.
2. 성공 기준과 수학 invariant를 명시한다.
3. 실패 테스트를 먼저 작성하고 예상한 이유로 실패함을 확인한다.
4. 테스트를 통과시키는 최소 production code만 구현한다.
5. focused test, 전체 CTest, Harness pytest를 실행한다.
6. Step summary와 output metadata가 실제 결과와 일치하는지 확인한다.
7. Phase 종료 전 전체 diff를 아키텍처·정식화·테스트 기준으로 review한다.
- `phases/internal-release/index.json`의 Step 0~4가 모두 `completed`
- `phases/index.json``internal-release``completed`
- PRD 8절의 모든 기준이 고유 checklist ID와 실제 증거에 연결됨
- Debug/Release build에 새 MSVC warning이 없음
- Debug/Release CTest Harness pytest가 0개가 아닌 상태로 모두 통과
- `cmake --install` 결과에 요구 binary/library/header/document/example/runtime
inventory가 포함됨
- source/build tree를 숨긴 install consumer와 CLI/HDF5 smoke test 통과
- 약 100,000 DOF benchmark의 correctness, stage time, memory와 환경 기록 완료
- Gate A PASS와 Gate B EVALUABLE의 의미 및 Abaqus stress 제한을 release 문서가
정확히 유지함
- 실패하거나 미실행인 증거가 없는 경우에만 내부 배포 완료 선언
## 9. 알려진 Harness 실행 이력과 복구 주의사항
새 세션의 권장 첫 요청은 다음과 같다.
이전 Phase에서 production 코드와 무관한 두 executor 문제가 있었다.
- Step 2 child sandbox에서 MSBuild가 사용자 profile의 Microsoft SDK/FileTracker
경로를 읽거나 child compiler process를 실행하지 못했다. 같은 checkout의 root
PowerShell에서는 정확한 build와 focused/full CTest가 통과했다.
- Step 3 child Codex 실행은 2026-07-30 usage limit에 도달해
`2026-08-05 15:07` 이후 재시도하라는 응답을 냈다. 당시 Step은 root 세션에서
TDD와 수용 조건을 직접 수행해 완료했다.
새 세션에서는 이 상태가 여전히 유효하다고 단정하지 말고 Harness를 한 번 정상
실행해 확인한다. 같은 외부 문제가 반복되면:
1. 저장소 변경 전후 상태와 정확한 실패 stage를 보존한다.
2. 제품 코드 문제인지 child sandbox/quota 문제인지 root의 동일 명령으로 분리한다.
3. 수동 fallback 시에도 Step의 TDD, focused/full test, review를 생략하지 않는다.
4. 실패 output을 성공으로 위장하지 말고 이후 성공 증거와 복구 commit을 분리해
metadata와 Git 이력에 남긴다.
5. 광범위한 사용자 profile 경로를 `--codex-add-dir`로 허용하지 않는다.
## 10. 다음 Phase 완료 조건
- `phases/fem-and-beam-kernel/index.json`의 Step 0~3이 모두 `completed`
- `phases/index.json`에서 `fem-and-beam-kernel``completed`
- formulation 문서와 production kernel의 DOF, 좌표계, 부호, 적분 규칙 일치
- focused test와 전체 CTest 통과
- Harness pytest가 0개가 아닌 상태로 전체 통과
- 새 MSVC 경고 없음
- 강체운동, 해석해, shear-dominant, 세장비와 회전 invariant 검증 통과
- `fem`/`elements`에 Abaqus 또는 외부 library API 유입 없음
- 코드 review의 Critical/Important 항목 해결
- 사용자 선택 전 원격 push나 `dev` 병합을 수행하지 않음
새 세션의 권장 첫 요청:
> `docs/HANDOFF.md`와 `phases/fem-and-beam-kernel/step0.md`부터 `step3.md`를 읽고
> 현재 baseline을 확인한 뒤 `fem-and-beam-kernel` Phase를 시작해주세요.
> `docs/HANDOFF.md`와 `internal-release`의 index/step0~4를 읽고 현재 `dev`
> baseline과 beam reference 검증 범위를 확인해주세요. Step 0과 Step 4의 stale
> reference 문구를 PRD/VALIDATION에 맞춘 뒤 `internal-release` Phase를 시작해주세요.
+168
View File
@@ -0,0 +1,168 @@
# FESA HDF5 Schema 2.0.0
## 1. Scope and version compatibility
Schema `2.0.0` is the self-contained Phase 1 result contract. One file contains
the active normalized model, its single linear-static analysis definition and
solver settings, and every nodal and Beam result needed without the source
`.inp` file. FESA performs no unit conversion.
Schema `1.0.0` was the earlier minimal vertical slice. Version `2.0.0` changes
the required model and result objects, so it is a new major version rather than
an in-place change to `1.0.0`. The current writer and reader accept exactly
`2.0.0`; every other version fails with `hdf5.unsupported_schema`. A future
reader may explicitly add support for compatible minor versions, but must not
infer compatibility from a version prefix.
## 2. Common rules
- Root attributes are variable-length UTF-8 strings:
`schema_version="2.0.0"`, `fesa_version`, and
`unit_policy="consistent_input_units_no_conversion"`, `input_source`, and
`input_fingerprint`. `input_source` is the UTF-8 path supplied to the solve
request. `input_fingerprint` is `fnv1a64:` followed by the 16 lowercase
hexadecimal digits of FNV-1a 64 over the original input bytes; it is a
reproducibility identifier, not a cryptographic integrity guarantee.
- Integer datasets use the stated little-endian fixed-width type. Floating
datasets use IEEE 754 little-endian `float64`. Strings are variable-length
UTF-8.
- `dense_index` is a contiguous 0-based row index. Semantic `internal_id`
values are nonnegative and are not assumed to be dense or ordered.
- Flat/orphan mesh `part_name` and `instance_name` values are empty strings.
- Ragged arrays use an offset dataset of length `row_count + 1`. Offsets start
at zero, are nondecreasing, and the final offset equals the flattened row
count. Input order is preserved.
- Numeric field datasets carry UTF-8 `coordinate_system` and `components`
attributes where listed. `components` is a comma-separated ordered list.
- Step and frame group names are contiguous decimal indices beginning at zero.
Phase 1 requires exactly one analysis step and one result step with the same
name, and exactly one result frame in that step.
## 3. Required objects
### 3.1 Model
Let `N`, `E`, `M`, `S`, `NS`, and `ES` be the node, Beam element, material,
section, node-set, and element-set counts. Let `P` be the total number of
section recovery points, `NM` the total node-set membership count, and `EM` the
total element-set membership count.
| Path | Type | Rank and shape | Attributes / meaning |
|---|---|---|---|
| `/model/nodes/dense_index` | `uint64` | 1, `[N]` | contiguous row index |
| `/model/nodes/internal_id` | `int64` | 1, `[N]` | `NodeId` |
| `/model/nodes/part_name` | UTF-8 | 1, `[N]` | entity provenance |
| `/model/nodes/instance_name` | UTF-8 | 1, `[N]` | entity provenance |
| `/model/nodes/local_label` | `int64` | 1, `[N]` | external local label |
| `/model/nodes/coordinates` | `float64` | 2, `[N,3]` | `coordinate_system="global"`, `components="X,Y,Z"` |
| `/model/elements/dense_index` | `uint64` | 1, `[E]` | contiguous row index |
| `/model/elements/internal_id` | `int64` | 1, `[E]` | `ElementId` |
| `/model/elements/part_name` | UTF-8 | 1, `[E]` | entity provenance |
| `/model/elements/instance_name` | UTF-8 | 1, `[E]` | entity provenance |
| `/model/elements/local_label` | `int64` | 1, `[E]` | external local label |
| `/model/elements/connectivity` | `uint64` | 2, `[E,2]` | node `dense_index`, ordered end `-1,+1` |
| `/model/elements/material_id` | `int64` | 1, `[E]` | references material `internal_id` |
| `/model/elements/section_id` | `int64` | 1, `[E]` | references section `internal_id` |
| `/model/materials/internal_id` | `int64` | 1, `[M]` | `MaterialId` |
| `/model/materials/name` | UTF-8 | 1, `[M]` | material name |
| `/model/materials/young_modulus` | `float64` | 1, `[M]` | finite, positive |
| `/model/materials/poisson_ratio` | `float64` | 1, `[M]` | finite, `-1 < nu < 0.5` |
| `/model/sections/internal_id` | `int64` | 1, `[S]` | `SectionId` |
| `/model/sections/name` | UTF-8 | 1, `[S]` | section name |
| `/model/sections/area` | `float64` | 1, `[S]` | `A` |
| `/model/sections/moment_y` | `float64` | 1, `[S]` | `Iy` |
| `/model/sections/moment_z` | `float64` | 1, `[S]` | `Iz` |
| `/model/sections/torsion_constant` | `float64` | 1, `[S]` | `J` |
| `/model/sections/shear_area_y` | `float64` | 1, `[S]` | applied `Asy` |
| `/model/sections/shear_area_z` | `float64` | 1, `[S]` | applied `Asz` |
| `/model/sections/shear_source` | `uint8` | 1, `[S]` | `0=input`, `1=phase1_default` |
| `/model/sections/orientation` | `float64` | 2, `[S,3]` | `coordinate_system="global"`, `components="X,Y,Z"` |
| `/model/sections/recovery_point_offsets` | `uint64` | 1, `[S+1]` | offsets into `recovery_points` |
| `/model/sections/recovery_points` | `float64` | 2, `[P,2]` | `coordinate_system="element_local"`, `components="y,z"` |
| `/model/sets/node/names` | UTF-8 | 1, `[NS]` | exact node-set names |
| `/model/sets/node/member_offsets` | `uint64` | 1, `[NS+1]` | offsets into `members` |
| `/model/sets/node/members` | `int64` | 1, `[NM]` | `NodeId`, input set/member order |
| `/model/sets/element/names` | UTF-8 | 1, `[ES]` | exact element-set names |
| `/model/sets/element/member_offsets` | `uint64` | 1, `[ES+1]` | offsets into `members` |
| `/model/sets/element/members` | `int64` | 1, `[EM]` | `ElementId`, input set/member order |
### 3.2 Analysis
`/analysis/steps/0` has UTF-8 attribute `name`. Let `B` be the prescribed-DOF
count and `L` the nodal-load count.
| Path | Type | Rank and shape | Attributes / meaning |
|---|---|---|---|
| `/analysis/steps/0/boundary_conditions/node_ids` | `int64` | 1, `[B]` | target `NodeId` |
| `/analysis/steps/0/boundary_conditions/dofs` | `uint8` | 1, `[B]` | Abaqus/FESA DOF number 1 through 6 |
| `/analysis/steps/0/boundary_conditions/values` | `float64` | 1, `[B]` | prescribed value |
| `/analysis/steps/0/nodal_loads/node_ids` | `int64` | 1, `[L]` | target `NodeId` |
| `/analysis/steps/0/nodal_loads/values` | `float64` | 2, `[L,6]` | `coordinate_system="global"`, `components="Fx,Fy,Fz,Mx,My,Mz"` |
`/analysis/solver_settings` has the exact UTF-8 attributes used by the Phase 1
pipeline: `backend="mkl_pardiso"`,
`matrix_storage="symmetric_upper_csr"`,
`matrix_type="symmetric_positive_definite"`,
`constraint_method="essential_dof_elimination"`, and
`assembly="deterministic_serial"`. No unused future settings are stored.
### 3.3 Results
`/results/steps/0` has UTF-8 attribute `name`; frame group `0` has scalar
`float64` attribute `step_time`. Let `RN`, `RE`, `RP`, and `D` be the nodal
result, Beam result, flattened Beam recovery-point, and diagnostic counts.
| Path below `/results/steps/0/frames/0` | Type | Rank and shape | Attributes / meaning |
|---|---|---|---|
| `nodal/node_ids` | `int64` | 1, `[RN]` | `NodeId`; provenance is joined from `/model/nodes` |
| `nodal/displacement` | `float64` | 2, `[RN,6]` | `coordinate_system="global"`, `components="Ux,Uy,Uz,Rx,Ry,Rz"` |
| `nodal/reaction` | `float64` | 2, `[RN,6]` | `coordinate_system="global"`, `components="RFx,RFy,RFz,RMx,RMy,RMz"` |
| `element/beam/element_ids` | `int64` | 1, `[RE]` | `ElementId` |
| `element/beam/part_name` | UTF-8 | 1, `[RE]` | result provenance |
| `element/beam/instance_name` | UTF-8 | 1, `[RE]` | result provenance |
| `element/beam/local_label` | `int64` | 1, `[RE]` | result provenance |
| `element/beam/local_frame` | `float64` | 3, `[RE,3,3]` | `coordinate_system="global"`, `components="ex,ey,ez"`; last dimension is `X,Y,Z` |
| `element/beam/end_node_ids` | `int64` | 2, `[RE,2]` | ordered ends `-1,+1` |
| `element/beam/xi` | `float64` | 2, `[RE,2]` | `components="end_minus,end_plus"` |
| `element/beam/section_strain` | `float64` | 3, `[RE,2,6]` | `coordinate_system="element_local"`, `components="epsilon,gamma_y,gamma_z,kappa_x,kappa_y,kappa_z"` |
| `element/beam/section_force` | `float64` | 3, `[RE,2,6]` | `coordinate_system="element_local"`, `components="N,Vy,Vz,T,My,Mz"` |
| `element/beam/centroid_sigma_xx` | `float64` | 2, `[RE,2]` | `coordinate_system="element_local"`, `components="end_minus,end_plus"`, `quantity="sigma_xx"` |
| `element/beam/recovery_point_offsets` | `uint64` | 1, `[RE+1]` | offsets into recovery-point rows |
| `element/beam/recovery_point_sigma_xx` | `float64` | 2, `[RP,2]` | `coordinate_system="element_local"`, `components="end_minus,end_plus"`, `quantity="sigma_xx"`; point order comes from the referenced section |
| `diagnostics/stage` | `uint8` | 1, `[D]` | enum table below |
| `diagnostics/severity` | `uint8` | 1, `[D]` | `0=warning`, `1=error` |
| `diagnostics/code` | UTF-8 | 1, `[D]` | exact diagnostic code |
| `diagnostics/message` | UTF-8 | 1, `[D]` | exact diagnostic message |
| `diagnostics/has_source` | `uint8` | 1, `[D]` | `0=no source`, `1=source present` |
| `diagnostics/source_file` | UTF-8 | 1, `[D]` | empty when source is absent |
| `diagnostics/source_line` | `uint64` | 1, `[D]` | zero when source is absent |
| `diagnostics/source_column` | `uint64` | 1, `[D]` | zero when source is absent |
Diagnostic stage encoding follows the declaration order:
| Value | Stage |
|---:|---|
| 0 | `io` |
| 1 | `syntax` |
| 2 | `semantic` |
| 3 | `model` |
| 4 | `equation` |
| 5 | `solver` |
| 6 | `results` |
| 7 | `validation` |
## 4. Writer and reader contract
- The writer validates `ResultDatabase`, exact schema version, model/result ID
and provenance joins, Beam connectivity, recovery-point counts, and the
single-step name before creating the file.
- Required-object creation, write, flush, and close failures become
`DiagnosticStage::results` errors. A failed write is never reported as
success.
- The reader validates the exact version, required datatypes, ranks, shapes,
offsets, finite values, uniqueness, references, field metadata, and result
contracts. It does not return a partial database or partial snapshots.
- The public reader returns adapter-owned, read-only metadata, model, and
analysis snapshots plus the semantic `ResultDatabase`. No HDF5 object or
handle escapes the adapter, and the snapshots contain enough information to
reconstruct the Phase 1 model and run definition without the source deck.
+40 -13
View File
@@ -163,8 +163,8 @@ HDF5 결과는 다음 정보를 함께 갖는 자기완결형 파일이어야
- 여러 재료·단면과 중첩 집합
4. Reference 테스트
- Abaqus/Standard 2024 B31 결과
- 현재 캔틸레버의 변위 반력
- 요소 내력 및 요소 절점 단면 도심 응력 비교 계약의 synthetic CSV 검증
- 현재 캔틸레버의 변위, 반력 및 요소 단면력
- 요소 절점 단면 도심 응력 비교 계약의 synthetic CSV 검증
### 5.2 골든 데이터
@@ -174,8 +174,15 @@ Abaqus는 CI나 Harness에서 자동 실행하지 않는다. 별도 Abaqus 2024
비교 실행은 물리량과 해당 CSV 경로를 명시한다. 요청한 파일이 없으면 실패하고,
요청하지 않은 물리량은 통과로 보고하지 않는다. 현재 `reference/cantilever beam`
샘플은 변위 반력 비교한다. 요소 내력과 응력 CSV가 추가되기 전까지 해당
reader와 비교 kernel은 synthetic CSV로 검증한다.
샘플은 변위, 반력 및 요소 단면력을 비교한다. 요소 응력 CSV가 추가되기 전까지
해당 reader와 비교 kernel은 synthetic CSV로 검증한다.
FESA 정식화 적합성과 Abaqus 결과 상관성은 별도 gate로 운영한다. FESA 적합성
gate는 `SCF=0`, 선택적 감차적분 및 문서화된 FESA 정식화를 해석해와 physics
invariant로 엄격히 검증한다. Abaqus 상관성 gate는 `SCF=0.25`를 포함할 수 있는
원본 Abaqus 모델과 CSV를 보존하고, 동일한 기하·재료·하중에 `SCF=0`을 적용한
별도 FESA 입력을 production pipeline으로 해석해 결과 차이를 정량화한다. 상관성
gate는 서로 다른 정식화의 수치 일치를 주장하지 않는다.
CSV 식별 및 값 열:
@@ -185,17 +192,35 @@ CSV 식별 및 값 열:
`SF-SF1..SF-SF3`, `SM-SM1..SM-SM3`
- 요소 응력: `Part Instance Name`, `Element Label`, `Node Label`, `Sxx`
단일 Instance에서는 `Part Instance Name` 열을 생략할 수 있다. 내력은
`SF1,SF2,SF3,SM1,SM2,SM3`을 각각 \(N,V_y,V_z,T,M_y,M_z\)로 비교한다.
응력은 요소 절점의 단면 도심값 \(\sigma_{xx}=N/A\)를 비교한다.
단일 Instance에서는 `Part Instance Name` 열을 생략할 수 있다. Abaqus Beam의
단면축 \((\mathbf n_1,\mathbf n_2)\)를 FESA의 \((\mathbf e_y,\mathbf e_z)\)와
일치시킨 입력에서 요소 내력은 Abaqus CSV 순서를
`SF1,SF3,SF2,SM3,SM1,SM2`로 재배열해 FESA의
\(N,V_y,V_z,T,M_y,M_z\)와 비교한다. 응력은 요소 절점의 단면 도심값
\(\sigma_{xx}=N/A\)를 비교한다.
### 5.3 허용오차
- 단위·정식화 테스트는 정규화된 엄격한 tolerance를 사용한다.
- Abaqus 비교 기본 상대오차는 \(10^{-5}\)로 한다.
- 영에 가까운 결과는 특성 길이, 하중 및 응력에 기반한 절대오차를 함께 사용한다.
- formulation 또는 output 위치 차이로 별도 tolerance가 필요하면 comparison
test 설정과 `docs/VALIDATION.md`에 근거를 기록한다.
- FESA 단위·정식화 적합성 gate는 정규화된 엄격한 tolerance를 사용한다.
- 정식화가 일치하는 reference 비교 기본 상대오차는 \(10^{-5}\)로 한다.
- Abaqus B31과 FESA Beam의 정식화가 다른 상관성 gate는 component별 RMSE와
Relative L2를 보고한다. 물리량과 component가 다른 값을 하나의 norm으로
혼합하지 않는다.
- component \(c\)의 값 쌍을 \((F_{ic},A_{ic})\), characteristic absolute scale을
\(s_c\)라 하면
\[
\operatorname{RMSE}_c=
\sqrt{\frac{1}{n}\sum_i(F_{ic}-A_{ic})^2},\qquad
\operatorname{RelativeL2}_c=
\frac{\sqrt{\sum_i(F_{ic}-A_{ic})^2}}
{\max\left(\sqrt{\sum_iA_{ic}^2},\sqrt{n}s_c\right)}.
\]
- Abaqus 상관성 gate의 성공은 요청된 모든 entity/component가 매칭되고 유한한
metric이 생성됨을 뜻한다. 관측된 단일 샘플에 맞춘 임의 pass/fail tolerance는
두지 않는다. 이후 acceptance envelope를 추가하려면 해석적 또는 mesh study
근거와 함께 `docs/VALIDATION.md`에 사전 기록한다.
## 6. 개발 워크플로우
@@ -235,6 +260,8 @@ CSV 식별 및 값 열:
- 테스트 0개 수집이 아님을 확인
- 전체 입력-해석-출력 통합 테스트 통과
- physics sanity와 평형 잔차 기준 통과
- 현재 Abaqus 2024 변위·반력 골든 결과의 tolerance 통과
- FESA Beam 정식화 적합성 gate의 엄격한 tolerance 통과
- 현재 Abaqus 2024 변위·반력·요소 단면력과의 component별 RMSE 및 Relative L2
상관성 보고서 생성
- 요소 내력·도심 응력 CSV adapter와 비교 kernel의 synthetic 검증 통과
- HDF5 schema, 입력 부분집합, 정식화 및 검증 보고서 제공
+132
View File
@@ -0,0 +1,132 @@
# FESA Phase 1 검증 보고서
## 1. 검증 기준과 실행 증거
이 보고서는 2026-08-03에 `feat-beam-reference-qualification`에서 새 Debug 빌드와
테스트를 실행해 수집한 결과다. FESA Beam 검증은 다음 두 gate를 독립적으로
운영한다.
- Gate A — FESA 정식화 적합성: FESA가 채택한 선택적 감차적분 Timoshenko
정식화를 해석해, 에너지, 강체 mode와 평형으로 검증한다.
- Gate B — Abaqus 결과 상관성: Abaqus B31과 FESA의 정식화 차이를 인정하고
요청된 결과의 component별 RMSE와 Relative L2를 보고한다.
실행 결과는 다음과 같다.
| 명령 | 결과 |
|------|------|
| `cmake --build --preset windows-debug` | MSVC v145 Debug x64 빌드 성공, 새 경고 없음 |
| `ctest --preset windows-debug --output-on-failure` | 68/68 통과 |
| `uv run --with pytest python -m pytest -v -rs` | 20/20 통과 |
| `fesa solve``fesa-reference-compare` | 42개 요청 위치와 18개 component metric 생성, exit code 0 |
Harness child의 관리형 WindowsApps PowerShell은 `CreateProcessAsUserW` 오류로 명령을
시작하지 못했다. 위 명령은 동일 저장소와 preset에서 현재 세션이 직접 실행했으며,
executor 환경 문제를 테스트 성공으로 간주하지 않았다.
## 2. Gate A — FESA 정식화 적합성
### 2.1 요소 정식화와 회복
`tests/unit/elements/beam3d2_test.cpp`의 다음 검증이 모두 통과했다.
| 검증 | 증거 | 판정 기준 |
|------|------|-----------|
| 행렬 유한성·대칭성 | `Beam3D2.ProducesFiniteSymmetricLocalAndGlobalStiffness` | local/global 12×12 전 항 유한, 대칭 항 bitwise equality |
| 강체 mode | `RigidBody.SixIndependentModesHaveZeroStrainEnergy` | 3개 병진과 3개 회전 mode의 energy가 roundoff bound 이내 |
| 축·비틀림 | `Timoshenko.ReproducesAnalyticalAxialSubmatrix`, `ReproducesAnalyticalTorsionalSubmatrix` | 해석 stiffness 대비 `512*epsilon` 상대 규모 |
| 굽힘 y/z | `Timoshenko.ReproducesConstantCurvatureEnergyAboutLocalY/LocalZ` | 해석 strain energy 대비 `512*epsilon` 상대 규모 |
| 전단 y/z | `Timoshenko.ReproducesConstantShearEnergyInLocalY/LocalZ` | 해석 strain energy 대비 `512*epsilon` 상대 규모 |
| 세장비 | `Timoshenko.AvoidsShearLockingAcrossSlendernessSweep` | `L/h={2,10,100,1000}`, `4096*epsilon*(L/h)^2` 상대 규모 |
| 회전 불변성 | `Beam3D2.PreservesGlobalEnergyUnderRigidCoordinateRotation` | 회전 전후 global energy가 `4096*epsilon` 상대 규모 이내 |
| 단면력 회복 | `BeamRecovery.*`, `SectionForce.*` | 축력, 비틀림, 전단, My, Mz와 biaxial 부호를 양 끝에서 검증 |
### 2.2 Physics sanity와 결정성
| 검증 | 증거 | 결과 |
|------|------|------|
| full-system 평형 | `StaticEquilibrium.ReturnedFieldsSatisfyOriginalFullEquation` | `Ku-f-r` 각 DOF가 `1e-12` 이내 |
| 비영 지정 DOF | `ConstraintElimination.ShiftsNonzeroPrescribedValueAndPreservesOriginalSystem``LinearStaticAnalysis.SolvesAllConstrainedSystemWithoutPardiso` | reduced RHS 이동, full displacement와 reaction 검증 |
| 반력 | `Reaction.UsesOriginalFullEquilibriumEquation` | 원래 full equation에서 회복 |
| production pipeline | `MinimalCantileverPipeline.WritesReadableFiniteEquilibratedResults` | parser→solver→HDF5 결과 유한, 전역 반력 평형 `1e-12` 이내 |
| reference 평형 | `CantileverReference.CorrelatesAllAvailableAbaqusResults` | 절대 `1e-12`와 적용 하중 대비 상대 `2e-13` 중 큰 허용치 이내 |
| 병렬 결정성 | `ThreadCountDeterminism.AssemblyAndLinearStateAreBitwiseStable` | thread `{1,2,16}`에서 10회 반복, matrix/displacement/reaction bitwise 동일 |
Gate A disposition은 **PASS**다.
## 3. Gate B — Abaqus 2024 결과 상관성
### 3.1 모델과 비교 계약
- Abaqus provenance: `reference/cantilever beam/cantilever beam.inp`와 변위, 반력,
요소 단면력 CSV. 원본 모델의 `SCF=0.25`는 유지한다.
- FESA projection: `reference/cantilever beam/cantilever beam fesa.inp`. 기하, 재료,
단면, 명시적 전단강성, 경계조건과 하중은 같고 `SCF=0`만 적용한다.
- entity join: nodal 결과 11개씩은 `(Instance, Node Label)`, 요소 단면력 20개는
`(Instance, Element Label, End Node Label)`로 매칭한다.
- absolute scale: displacement `1e-10`, reaction `1e-8`, internal force `1e-8`.
- 서로 단위가 다른 component는 하나의 norm으로 합치지 않는다.
요소력 축 순서는 다음과 같다.
| FESA | Abaqus CSV |
|------|------------|
| `N` | `SF1` |
| `Vy` | `SF3` |
| `Vz` | `SF2` |
| `T` | `SM3` |
| `My` | `SM1` |
| `Mz` | `SM2` |
component \(c\)의 표본 수를 \(n\), 오차를 \(e_i=a_i-r_i\), absolute scale을
\(s_i\)라 하면 다음을 사용한다.
\[
\operatorname{RMSE}_c=\sqrt{\frac{1}{n}\sum_i e_i^2}
\]
\[
\operatorname{RelativeL2}_c=
\frac{\sqrt{\sum_i e_i^2}}
{\max\left(\sqrt{\sum_i r_i^2},\sqrt{\sum_i s_i^2}\right)}
\]
### 3.2 수집된 correlation metric
| 물리량 | component | count | RMSE | Relative L2 |
|--------|-----------|------:|-----:|------------:|
| displacement | Ux | 11 | 0 | 0 |
| displacement | Uy | 11 | 0 | 0 |
| displacement | Uz | 11 | 2.1970445023044601e-05 | 2.2349969291714038e-03 |
| displacement | Rx | 11 | 0 | 0 |
| displacement | Ry | 11 | 2.1672945177750166e-09 | 1.0416581667076092e-06 |
| displacement | Rz | 11 | 0 | 0 |
| reaction | RFx | 11 | 0 | 0 |
| reaction | RFy | 11 | 0 | 0 |
| reaction | RFz | 11 | 1.3385150002853335e-07 | 4.4393520322089023e-13 |
| reaction | RMx | 11 | 0 | 0 |
| reaction | RMy | 11 | 7.8149308366928907e-07 | 2.591919334788851e-13 |
| reaction | RMz | 11 | 0 | 0 |
| internal force | N | 20 | 0 | 0 |
| internal force | Vy | 20 | 0 | 0 |
| internal force | Vz | 20 | 2.7107472798172421e-07 | 2.7107472798172426e-13 |
| internal force | T | 20 | 0 | 0 |
| internal force | My | 20 | 474341.64902538335 | 0.082541022764834646 |
| internal force | Mz | 20 | 0 | 0 |
모든 요청 entity/component가 유일하게 매칭됐고 metric이 유한하므로 Gate B
disposition은 **EVALUABLE**이다. 이는 Abaqus B31과 FESA의 정식화가 동일하거나
임의 정확도 threshold를 통과했다는 뜻이 아니다. 현재 가장 큰 Relative L2는
`My=0.082541022764834646`, 다음은 `Uz=0.0022349969291714038`이다. pass/fail
envelope는 mesh 또는 해석적 연구 근거 없이 이 관측값에 맞춰 설정하지 않는다.
## 4. 검증 범위와 제한
- Abaqus 단면 도심 응력 CSV는 제공되지 않았다. `StressCsv`는 synthetic CSV schema와
`(Instance, Element, End Node)` 매핑만 검증하며, 응력은 **not yet
Abaqus-qualified**다.
- Gate B는 solver 간 correlation이며 Gate A의 해석적 정확도 검증을 대체하지 않는다.
- Abaqus golden 갱신에는 Abaqus 2024 환경과 수동 provenance 검토가 필요하다.
- FESA는 단위 변환을 수행하지 않으므로 입력과 CSV가 일관 단위계를 사용해야 한다.
현재 이중 gate 결론은 **Gate A PASS / Gate B EVALUABLE**이다.
+237
View File
@@ -0,0 +1,237 @@
# 2절점 3D Timoshenko Beam 정식화
## 1. 범위와 가정
FESA Phase 1 Beam은 소변형, 선형 탄성의 2절점 직선
isoparametric Timoshenko 요소다. 단면은 도심과 전단중심이 일치하는 주축 단면이며
\(I_{yz}=0\)이다. 단면 offset, warping, 기하·재료 비선형은 포함하지 않는다.
FESA는 단위 변환을 하지 않으므로 모든 입력은 하나의 일관 단위계를 사용해야 한다.
정식화의 유한요소 이산화와 수치적분은 Bathe,
*Finite Element Procedures*, 2판과 Hughes,
*The Finite Element Method: Linear Static and Dynamic Finite Element
Analysis*를 따른다. 3차원 isoparametric Beam의 국부 좌표계는 Bathe와
Bolourchi(1979)의 선형화된 부분을 사용한다. 선택적 감차적분은 Hughes, Taylor와
Kanoknukulchai(1977)의 원칙을 Beam 전단항에 적용한다. Plate 운동학은 사용하지
않는다.
## 2. 국부 좌표계와 자유도
절점 1에서 절점 2로 향하는 단위벡터를 \(\mathbf e_x\)로 둔다. 입력 orientation
벡터에서 \(\mathbf e_x\) 성분을 Gram-Schmidt로 제거하고 정규화한 벡터를
\(\mathbf e_y\)로 두며,
\[
\mathbf e_z=\mathbf e_x\times\mathbf e_y
\]
로 정의한다. 따라서 \((\mathbf e_x,\mathbf e_y,\mathbf e_z)\)는 오른손 직교
기저다. 길이가 0인 요소, 길이가 0인 orientation, 요소축과 평행한 orientation은
오류다.
국부 요소 자유도 벡터의 순서는 다음과 같다.
\[
\mathbf q_l =
\begin{bmatrix}
u_1&v_1&w_1&\theta_{x1}&\theta_{y1}&\theta_{z1}&
u_2&v_2&w_2&\theta_{x2}&\theta_{y2}&\theta_{z2}
\end{bmatrix}^{T}.
\]
\(u,v,w\)는 각각 국부 \(x,y,z\) 병진이고
\(\theta_x,\theta_y,\theta_z\)는 오른손 규칙의 회전 성분이다. 전역 자유도도 절점별
\([u_x,u_y,u_z,r_x,r_y,r_z]\) 순서를 사용한다.
기준축에서 떨어진 점 \((0,y,z)\)의 선형화된 변위장은
\[
U_x=u+z\theta_y-y\theta_z,\qquad
U_y=v-z\theta_x,\qquad
U_z=w+y\theta_x
\]
이다.
## 3. 변형률과 부호
일반화 변형률과 대응 단면력의 component 순서는
\[
\boldsymbol\varepsilon_s =
\begin{bmatrix}
\epsilon&\gamma_y&\gamma_z&\kappa_x&\kappa_y&\kappa_z
\end{bmatrix}^{T},
\qquad
\mathbf s =
\begin{bmatrix}
N&V_y&V_z&T&M_y&M_z
\end{bmatrix}^{T}
\]
로 고정한다. 위 변위장과 오른손 회전 부호로부터
\[
\begin{aligned}
\epsilon &= \frac{du}{dx},&
\gamma_y &= \frac{dv}{dx}-\theta_z,&
\gamma_z &= \frac{dw}{dx}+\theta_y,\\
\kappa_x &= \frac{d\theta_x}{dx},&
\kappa_y &= \frac{d\theta_y}{dx},&
\kappa_z &= \frac{d\theta_z}{dx}
\end{aligned}
\]
를 사용한다. 양의 \(N,V_y,V_z,T,M_y,M_z\)는 위 일반화 변형률과 양의 내부
일률 \(\delta\boldsymbol\varepsilon_s^T\mathbf s\)을 이루는 방향이다. 이 부호에서
단면 축응력은
\[
\sigma_{xx}=E\left(\epsilon+z\kappa_y-y\kappa_z\right)
\]
이고 \(M_y=\int_A z\sigma_{xx}\,dA\),
\(M_z=-\int_A y\sigma_{xx}\,dA\)이다.
## 4. 재료와 단면 constitutive matrix
등방성 선형 탄성의 전단계수는 입력 \(E,\nu\)로부터
\[
G=\frac{E}{2(1+\nu)}
\]
로 계산한다. 일반화 constitutive 관계는
\[
\mathbf s=\mathbf D\boldsymbol\varepsilon_s,\qquad
\mathbf D=
\operatorname{diag}
\left(EA,\;GA_{sy},\;GA_{sz},\;GJ,\;EI_y,\;EI_z\right).
\]
\(A,I_y,I_z,J,A_{sy},A_{sz}\)는 양의 값이어야 한다. 명시적인 전단강성이 없는
입력은 semantic mapper가 커널 호출 전에
\(A_{sy}=A_{sz}=5A/6\), `SCF=0`으로 정규화한다. 커널은 Abaqus의 slenderness
compensation을 추가하지 않는다.
## 5. 자연좌표, 형상함수와 Jacobian
자연좌표는 \(\xi\in[-1,1]\)이고 두 절점의 선형 형상함수는
\[
N_1(\xi)=\frac{1-\xi}{2},\qquad
N_2(\xi)=\frac{1+\xi}{2}.
\]
모든 여섯 component를 같은 형상함수로 보간한다. 직선 요소 길이를 \(L\)이라 하면
\[
J=\frac{dx}{d\xi}=\frac{L}{2},\qquad
\frac{dN_a}{dx}=\frac{1}{J}\frac{dN_a}{d\xi}
\]
이다.
절점 \(a\)의 자유도 순서
\([u_a,v_a,w_a,\theta_{xa},\theta_{ya},\theta_{za}]\)에 대한
strain-displacement block은
\[
\mathbf B_a(\xi)=
\begin{bmatrix}
N_{a,x}&0&0&0&0&0\\
0&N_{a,x}&0&0&0&-N_a\\
0&0&N_{a,x}&0&N_a&0\\
0&0&0&N_{a,x}&0&0\\
0&0&0&0&N_{a,x}&0\\
0&0&0&0&0&N_{a,x}
\end{bmatrix},
\qquad
\mathbf B=\begin{bmatrix}\mathbf B_1&\mathbf B_2\end{bmatrix}.
\]
## 6. 선택적 감차적분과 요소 강성
Constitutive matrix를 축·비틀림·굽힘 부분과 전단 부분으로 나눈다.
\[
\begin{aligned}
\mathbf D_{ab}&=
\operatorname{diag}(EA,0,0,GJ,EI_y,EI_z),\\
\mathbf D_s&=
\operatorname{diag}(0,GA_{sy},GA_{sz},0,0,0).
\end{aligned}
\]
국부 강성은
\[
\mathbf K_l =
\sum_{g=1}^{2}
\mathbf B(\xi_g)^T\mathbf D_{ab}\mathbf B(\xi_g)Jw_g
+
\sum_{g=1}^{1}
\mathbf B(\xi_g)^T\mathbf D_s\mathbf B(\xi_g)Jw_g
\]
로 직접 적분한다. 축·비틀림·굽힘은
\(\xi=\pm1/\sqrt{3},w=1\)인 2점 Gauss rule을 사용하고, 전단은
\(\xi=0,w=2\)인 1점 Gauss rule만 사용한다. 전단항을 2점 적분하지 않는다.
선형 회전장에서 1점 전단 적분은 요소 중앙의 평균 회전을 사용하므로 일정곡률
굽힘 mode의 불필요한 전단에너지를 제거한다.
`Matrix12``matrix[row][column]`의 12행 12열 고정 크기 저장을 사용한다. 적분은
각 Gauss point에서 \(\mathbf B^T\mathbf D\mathbf B\)의 상삼각 entry를 고정된
부동소수점 연산 순서로 누적한 뒤 하삼각에 복사한다.
## 7. 전역 변환
기존 `beam_transformation`이 만드는 \(\mathbf T\)는 각 병진·회전 3-vector를
전역 성분에서 국부 성분으로 투영한다.
\[
\mathbf q_l=\mathbf T\mathbf q_g.
\]
따라서 전역 강성과 에너지 관계는
\[
\mathbf K_g=\mathbf T^T\mathbf K_l\mathbf T,\qquad
\frac12\mathbf q_g^T\mathbf K_g\mathbf q_g
=\frac12\mathbf q_l^T\mathbf K_l\mathbf q_l
\]
이다.
## 8. 구현과 검증 invariant
커널은 위 적분식으로 local stiffness를 계산하고 변환식으로 global stiffness를
계산한다. 다음 invariant가 정식화 검증 기준이다.
- local/global stiffness는 roundoff 이내에서 대칭이다.
- 세 병진과 세 회전 강체 mode의 변형에너지는 0이다.
- 축 및 비틀림 submatrix는 각각 \(EA/L\), \(GJ/L\)을 재현한다.
- 일정곡률 단축·이축 굽힘 에너지는 각각
\(EI_y\kappa_y^2L/2\), \(EI_z\kappa_z^2L/2\)이다.
- 일정전단 mode의 에너지는 각각
\(GA_{sy}\gamma_y^2L/2\), \(GA_{sz}\gamma_z^2L/2\)이다.
- 일정곡률 mode의 전단에너지는 slenderness 변화와 무관하게 0이며 전체 에너지는
해석값과 일치한다.
- 요소와 자유도를 같은 강체 회전으로 회전하면 변형에너지가 보존된다.
## 9. 참고문헌과 FESA 차이
- K. J. Bathe, *Finite Element Procedures*, 2nd ed.
- T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic
Finite Element Analysis*.
- K. J. Bathe and S. Bolourchi, “Large Displacement Analysis of
Three-Dimensional Beam Structures,” 1979. FESA는 이 문헌의 선형화된 국부
좌표와 변환만 사용한다.
- T. J. R. Hughes, R. L. Taylor, and W. Kanoknukulchai, “A Simple and
Efficient Finite Element for Plate Bending,” 1977. FESA는 selective
integration 원칙만 사용한다.
- Abaqus 2024, “Choosing a Beam Element” 및 “BEAM GENERAL SECTION.”
FESA의 \(5A/6\) 생략값과 `SCF=0`은 Phase 1 명시적 계약이며 임의 일반 단면의
보편적인 Abaqus 기본값이 아니다. Abaqus B31 slenderness compensation도
구현하지 않는다.
@@ -0,0 +1,23 @@
#pragma once
#include <optional>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/model/domain.hpp>
#include <fesa/results/result_database.hpp>
namespace fesa {
struct AnalysisRunResult final {
bool succeeded;
std::optional<ResultDatabase> results;
std::vector<Diagnostic> diagnostics;
};
class LinearStaticAnalysis final {
public:
[[nodiscard]] AnalysisRunResult run(const Domain& domain) const;
};
} // namespace fesa
+17
View File
@@ -0,0 +1,17 @@
#pragma once
#include <filesystem>
#include <fesa/analysis/linear_static_analysis.hpp>
namespace fesa {
struct AnalysisRequest final {
std::filesystem::path input_path;
std::filesystem::path output_path;
};
[[nodiscard]] AnalysisRunResult run_solver(
const AnalysisRequest& request);
} // namespace fesa
+21
View File
@@ -0,0 +1,21 @@
#pragma once
#include <cstddef>
#include <fesa/assembly/equation_system.hpp>
#include <fesa/fem/dof_manager.hpp>
#include <fesa/model/domain.hpp>
namespace fesa {
struct AssemblyOptions final {
std::size_t max_threads;
std::size_t grain_size;
};
[[nodiscard]] EquationSystem assemble_parallel(
const Domain& domain,
const DofManager& dofs,
AssemblyOptions options);
} // namespace fesa
+28
View File
@@ -0,0 +1,28 @@
#pragma once
#include <cstddef>
#include <cstdint>
#include <span>
#include <vector>
#include <fesa/assembly/symmetric_csr.hpp>
#include <fesa/model/ids.hpp>
namespace fesa {
struct MatrixContribution final {
std::size_t row;
std::size_t column;
ElementId element;
std::uint16_t local_order;
double value;
};
[[nodiscard]] std::vector<MatrixContribution> canonicalize_contributions(
std::span<const MatrixContribution> contributions);
[[nodiscard]] SymmetricCsr merge_contributions(
std::size_t order,
std::span<const MatrixContribution> canonical);
} // namespace fesa
+14
View File
@@ -0,0 +1,14 @@
#pragma once
#include <vector>
#include <fesa/assembly/symmetric_csr.hpp>
namespace fesa {
struct EquationSystem final {
SymmetricCsr stiffness;
std::vector<double> force;
};
} // namespace fesa
@@ -0,0 +1,13 @@
#pragma once
#include <fesa/assembly/equation_system.hpp>
#include <fesa/fem/dof_manager.hpp>
#include <fesa/model/domain.hpp>
namespace fesa {
[[nodiscard]] EquationSystem assemble_serial(
const Domain& domain,
const DofManager& dofs);
} // namespace fesa
+18
View File
@@ -0,0 +1,18 @@
#pragma once
#include <cstddef>
#include <cstdint>
#include <vector>
namespace fesa {
// Zero-based CSR containing only the upper triangle. Column indices are
// strictly increasing within each row.
struct SymmetricCsr final {
std::size_t order;
std::vector<std::int32_t> row_offsets;
std::vector<std::int32_t> column_indices;
std::vector<double> values;
};
} // namespace fesa
+31
View File
@@ -0,0 +1,31 @@
#pragma once
#include <optional>
#include <span>
#include <vector>
#include <fesa/assembly/equation_system.hpp>
#include <fesa/core/diagnostic.hpp>
#include <fesa/fem/dof_manager.hpp>
namespace fesa {
struct ReducedSystem final {
SymmetricCsr stiffness;
std::vector<double> force;
};
struct ConstraintResult final {
std::optional<ReducedSystem> reduced_system;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] ConstraintResult eliminate_essential_bcs(
const EquationSystem& original,
const DofManager& dofs);
[[nodiscard]] std::vector<double> recover_reaction(
const EquationSystem& original,
std::span<const double> full_displacement);
} // namespace fesa
+52
View File
@@ -0,0 +1,52 @@
#pragma once
#include <array>
#include <optional>
#include <span>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/core/vec3.hpp>
#include <fesa/fem/beam_frame.hpp>
#include <fesa/model/beam_section.hpp>
#include <fesa/model/ids.hpp>
#include <fesa/model/material.hpp>
namespace fesa {
struct Beam3D2Input final {
std::array<Vec3, 2> coordinates;
std::array<NodeId, 2> node_ids;
IsotropicElastic material;
BeamSection section;
};
struct BeamSectionResult final {
double xi;
NodeId end_node;
std::array<double, 6> section_strain;
std::array<double, 6> section_force;
double centroid_sigma_xx;
std::vector<double> sigma_xx;
};
struct Beam3D2Contribution final {
Matrix12 local_stiffness;
Matrix12 global_stiffness;
BeamFrame frame;
};
struct BeamKernelResult final {
std::optional<Beam3D2Contribution> contribution;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] BeamKernelResult compute_beam3d2(
const Beam3D2Input& input);
[[nodiscard]] std::vector<BeamSectionResult> recover_beam3d2(
const Beam3D2Input& input,
std::span<const double, 12> element_displacement,
std::span<const std::array<double, 2>> recovery_points);
} // namespace fesa
+32
View File
@@ -0,0 +1,32 @@
#pragma once
#include <array>
#include <optional>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/core/vec3.hpp>
namespace fesa {
using Matrix12 = std::array<std::array<double, 12>, 12>;
struct BeamFrame final {
Vec3 ex;
Vec3 ey;
Vec3 ez;
};
struct BeamFrameResult final {
std::optional<BeamFrame> frame;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] BeamFrameResult make_beam_frame(
const Vec3& first,
const Vec3& second,
const Vec3& orientation);
[[nodiscard]] Matrix12 beam_transformation(const BeamFrame& frame);
} // namespace fesa
+58
View File
@@ -0,0 +1,58 @@
#pragma once
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <optional>
#include <span>
#include <vector>
#include <fesa/model/domain.hpp>
namespace fesa {
enum class NodeDof : std::uint8_t {
ux,
uy,
uz,
rx,
ry,
rz,
};
struct DofAddress final {
NodeId node;
NodeDof dof;
};
class DofManager final {
public:
[[nodiscard]] static DofManager build(const Domain& domain);
[[nodiscard]] std::size_t full_dof_count() const noexcept;
[[nodiscard]] std::size_t free_equation_count() const noexcept;
[[nodiscard]] std::optional<std::size_t> equation(
DofAddress address) const;
[[nodiscard]] std::optional<std::size_t> equation(
std::size_t full_dof) const;
[[nodiscard]] std::optional<double> prescribed_value(
DofAddress address) const;
[[nodiscard]] std::optional<double> prescribed_value(
std::size_t full_dof) const;
[[nodiscard]] std::size_t full_dof(DofAddress address) const;
[[nodiscard]] std::array<std::size_t, 12> element_full_dofs(
const BeamElement& element) const;
[[nodiscard]] std::vector<double> reconstruct_full(
std::span<const double> reduced) const;
private:
[[nodiscard]] std::size_t node_full_dof_base(NodeId node) const;
std::map<std::int64_t, std::size_t> node_full_dof_bases_;
std::vector<std::optional<std::size_t>> equations_;
std::vector<double> prescribed_values_;
std::size_t free_equation_count_ = 0;
};
} // namespace fesa
+14
View File
@@ -0,0 +1,14 @@
#pragma once
#include <span>
namespace fesa {
struct GaussPoint1D final {
double xi;
double weight;
};
[[nodiscard]] std::span<const GaussPoint1D> gauss_rule_1d(int order);
} // namespace fesa
+13
View File
@@ -0,0 +1,13 @@
#pragma once
#include <array>
namespace fesa {
[[nodiscard]] std::array<double, 2> line2_shape(double xi);
[[nodiscard]] std::array<double, 2> line2_shape_derivative(double xi);
[[nodiscard]] double line2_jacobian(double length);
} // namespace fesa
+28
View File
@@ -0,0 +1,28 @@
#pragma once
#include <optional>
#include <span>
#include <string>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/io/abaqus/deck_record.hpp>
namespace fesa {
struct ActiveInputView final {
bool flat = true;
std::string part_name;
std::string instance_name;
std::span<const DeckRecord> part_records;
std::span<const DeckRecord> assembly_records;
};
struct ActiveInputResult final {
std::optional<ActiveInputView> input;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] ActiveInputResult select_active_input(const ParsedDeck& deck);
} // namespace fesa
+2
View File
@@ -15,6 +15,7 @@ struct DeckRecord final {
std::map<std::string, std::string, std::less<>> parameters;
std::vector<std::vector<std::string>> data;
SourceLocation source;
std::vector<SourceLocation> data_sources;
};
struct ParsedPart final {
@@ -27,6 +28,7 @@ struct ParsedInstance final {
std::string name;
std::string part_name;
std::vector<std::vector<std::string>> transform_data;
std::vector<SourceLocation> transform_sources;
SourceLocation source;
};
+2
View File
@@ -2,6 +2,7 @@
#include <filesystem>
#include <optional>
#include <string>
#include <vector>
#include <fesa/core/diagnostic.hpp>
@@ -12,6 +13,7 @@ namespace fesa {
struct ParseDeckResult final {
std::optional<ParsedDeck> deck;
std::vector<Diagnostic> diagnostics;
std::string input_fingerprint;
};
[[nodiscard]] ParseDeckResult parse_deck(
+31
View File
@@ -0,0 +1,31 @@
#pragma once
#include <cstdint>
#include <string>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/io/abaqus/deck_record.hpp>
namespace fesa {
enum class ResolvedSetScope { global, part, assembly };
enum class ResolvedSetKind { node, element };
struct ResolvedSet final {
std::string scope_name;
std::string set_name;
std::vector<std::int64_t> sorted_unique_labels;
ResolvedSetScope scope = ResolvedSetScope::global;
ResolvedSetKind kind = ResolvedSetKind::node;
};
struct SetResolutionResult final {
std::vector<ResolvedSet> sets;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] SetResolutionResult resolve_sets(const ParsedDeck& deck);
} // namespace fesa
+102
View File
@@ -0,0 +1,102 @@
#pragma once
#include <array>
#include <cstdint>
#include <filesystem>
#include <optional>
#include <string>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/core/vec3.hpp>
#include <fesa/model/beam_section.hpp>
#include <fesa/model/domain.hpp>
#include <fesa/model/entity_origin.hpp>
#include <fesa/model/ids.hpp>
#include <fesa/results/result_database.hpp>
namespace fesa {
struct Hdf5NodeSnapshot final {
std::uint64_t dense_index;
NodeId id;
EntityOrigin origin;
Vec3 coordinates;
};
struct Hdf5ElementSnapshot final {
std::uint64_t dense_index;
ElementId id;
EntityOrigin origin;
std::array<std::uint64_t, 2> connectivity;
MaterialId material;
SectionId section;
};
struct Hdf5SectionSnapshot final {
SectionId id;
std::string name;
double area;
double iy;
double iz;
double torsion_j;
double shear_area_y;
double shear_area_z;
ShearPropertySource shear_source;
Vec3 orientation;
std::vector<std::array<double, 2>> recovery_points;
};
struct Hdf5ModelSnapshot final {
std::vector<Hdf5NodeSnapshot> nodes;
std::vector<Hdf5ElementSnapshot> elements;
std::vector<IsotropicElastic> materials;
std::vector<Hdf5SectionSnapshot> sections;
std::vector<NodeSet> node_sets;
std::vector<ElementSet> element_sets;
};
struct Hdf5MetadataSnapshot final {
std::string schema_version;
std::string fesa_version;
std::string unit_policy;
std::string input_source;
std::string input_fingerprint;
};
struct Hdf5InputIdentity final {
std::string source;
std::string fingerprint;
};
struct Hdf5SolverSettingsSnapshot final {
std::string backend;
std::string matrix_storage;
std::string matrix_type;
std::string constraint_method;
std::string assembly;
};
struct Hdf5AnalysisSnapshot final {
StepDefinition step;
Hdf5SolverSettingsSnapshot solver;
};
struct Hdf5ReadResult final {
std::optional<ResultDatabase> database;
std::vector<Diagnostic> diagnostics;
std::optional<Hdf5ModelSnapshot> model;
std::optional<Hdf5MetadataSnapshot> metadata;
std::optional<Hdf5AnalysisSnapshot> analysis;
};
[[nodiscard]] std::vector<Diagnostic> write_hdf5(
const std::filesystem::path& path,
const Domain& domain,
const ResultDatabase& database,
const Hdf5InputIdentity& input_identity);
[[nodiscard]] Hdf5ReadResult read_hdf5_results(
const std::filesystem::path& path);
} // namespace fesa
+2
View File
@@ -57,6 +57,8 @@ public:
[[nodiscard]] const Node& node(NodeId id) const;
[[nodiscard]] const Node& node(const EntityOrigin& origin) const;
[[nodiscard]] const IsotropicElastic& material(MaterialId id) const;
[[nodiscard]] const BeamSection& section(SectionId id) const;
private:
friend class DomainBuilder;
+77
View File
@@ -0,0 +1,77 @@
#pragma once
#include <array>
#include <string>
#include <string_view>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/core/status.hpp>
#include <fesa/elements/beam/beam3d2.hpp>
#include <fesa/model/entity_origin.hpp>
#include <fesa/model/ids.hpp>
namespace fesa {
enum class FieldCoordinateSystem { global, element_local };
struct NodalFrame final {
inline static constexpr FieldCoordinateSystem coordinate_system =
FieldCoordinateSystem::global;
inline static constexpr std::array<std::string_view, 6>
displacement_components{
"Ux", "Uy", "Uz", "Rx", "Ry", "Rz"};
inline static constexpr std::array<std::string_view, 6>
reaction_components{
"RFx", "RFy", "RFz", "RMx", "RMy", "RMz"};
std::vector<NodeId> node_ids;
std::vector<EntityOrigin> origins;
std::vector<std::array<double, 6>> displacement;
std::vector<std::array<double, 6>> reaction;
};
struct BeamElementFrame final {
inline static constexpr FieldCoordinateSystem coordinate_system =
FieldCoordinateSystem::element_local;
inline static constexpr std::array<std::string_view, 6>
section_strain_components{
"epsilon", "gamma_y", "gamma_z", "kappa_x", "kappa_y",
"kappa_z"};
inline static constexpr std::array<std::string_view, 6>
section_force_components{
"N", "Vy", "Vz", "T", "My", "Mz"};
inline static constexpr std::string_view axial_stress_component =
"sigma_xx";
ElementId element;
EntityOrigin origin;
BeamFrame local_frame;
std::array<BeamSectionResult, 2> end_results;
};
struct ElementFrame final {
std::vector<BeamElementFrame> beams;
};
struct ResultFrame final {
double step_time;
NodalFrame nodal;
ElementFrame element;
std::vector<Diagnostic> diagnostics;
};
struct ResultStep final {
std::string name;
std::vector<ResultFrame> frames;
};
struct ResultDatabase final {
std::string schema_version;
std::vector<ResultStep> steps;
};
[[nodiscard]] Status validate_result_database(
const ResultDatabase& database);
} // namespace fesa
@@ -0,0 +1,26 @@
#pragma once
#include <span>
#include <vector>
#include <fesa/assembly/symmetric_csr.hpp>
#include <fesa/core/diagnostic.hpp>
namespace fesa {
struct LinearSolveResult final {
std::vector<double> solution;
double relative_residual{};
std::vector<Diagnostic> diagnostics;
};
class LinearSolver {
public:
virtual ~LinearSolver() = default;
[[nodiscard]] virtual LinearSolveResult solve(
const SymmetricCsr& matrix,
std::span<const double> rhs) = 0;
};
} // namespace fesa
@@ -0,0 +1,22 @@
#pragma once
#include <span>
#include <fesa/solvers/linear/linear_solver.hpp>
namespace fesa {
class PardisoLinearSolver final : public LinearSolver {
public:
PardisoLinearSolver();
~PardisoLinearSolver() override;
PardisoLinearSolver(const PardisoLinearSolver&) = delete;
PardisoLinearSolver& operator=(
const PardisoLinearSolver&) = delete;
[[nodiscard]] LinearSolveResult solve(
const SymmetricCsr& matrix,
std::span<const double> rhs) override;
};
} // namespace fesa
+79
View File
@@ -0,0 +1,79 @@
#pragma once
#include <cstddef>
#include <cstdint>
#include <optional>
#include <span>
#include <string>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/results/result_database.hpp>
namespace fesa {
enum class ReferenceQuantity {
displacement,
reaction,
internal_force,
centroid_stress
};
struct Tolerance final {
double relative;
double absolute_scale;
};
struct ResultPosition final {
std::string instance_name;
std::int64_t entity_label;
std::optional<std::int64_t> end_node_label;
};
struct ComparisonSample final {
ReferenceQuantity quantity;
ResultPosition position;
std::vector<double> reference;
std::vector<double> actual;
Tolerance tolerance;
};
struct ComparisonReport final {
bool passed;
double maximum_normalized_error;
std::vector<Diagnostic> 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<ComponentCorrelationMetric> metrics;
std::vector<Diagnostic> failures;
};
struct ComparisonSampleMatch final {
std::optional<ComparisonSample> sample;
std::vector<Diagnostic> failures;
};
[[nodiscard]] ComparisonSampleMatch make_comparison_sample(
const ResultFrame& frame,
ReferenceQuantity quantity,
const ResultPosition& position,
std::span<const double> reference,
Tolerance tolerance);
[[nodiscard]] ComparisonReport compare_samples(
std::span<const ComparisonSample> samples);
[[nodiscard]] CorrelationReport correlate_samples(
std::span<const ComparisonSample> samples);
} // namespace fesa
+28
View File
@@ -0,0 +1,28 @@
#pragma once
#include <filesystem>
#include <string_view>
#include <vector>
#include <fesa/core/diagnostic.hpp>
#include <fesa/validation/comparison.hpp>
namespace fesa {
struct ReferenceRow final {
ReferenceQuantity quantity;
ResultPosition position;
std::vector<double> values;
};
struct ReferenceCsvReadResult final {
std::vector<ReferenceRow> rows;
std::vector<Diagnostic> diagnostics;
};
[[nodiscard]] ReferenceCsvReadResult read_reference_csv(
ReferenceQuantity quantity,
const std::filesystem::path& path,
std::string_view single_instance_name);
} // namespace fesa
+23 -6
View File
@@ -5,27 +5,44 @@
{
"step": 0,
"name": "abaqus-input-contract",
"status": "pending"
"status": "completed",
"summary": "Defined the normative Phase 1 Abaqus subset and expanded the public parser/mapper fixture matrix to 74 cases covering strict scopes, parameter forms, data ownership, exact mesh rows and references in active and inactive Parts, duplicate Parts, Assembly-only hierarchical targets, and source-accurate diagnostics.",
"started_at": "2026-08-01T02:16:39+0900",
"completed_at": "2026-08-01T02:24:23+0900"
},
{
"step": 1,
"name": "part-and-assembly-set-resolution",
"status": "pending"
"status": "completed",
"summary": "Added deterministic scope- and kind-aware Part/Assembly set resolution with explicit, generated, nested, forward, duplicate, empty, cycle, missing-member, and active-Instance validation.",
"started_at": "2026-08-01T02:24:24+0900",
"completed_at": "2026-08-01T02:37:05+0900"
},
{
"step": 2,
"name": "single-instance-semantic-validation",
"status": "pending"
"status": "completed",
"summary": "Added ActiveInputView selection and semantic validation for exclusive flat or single untransformed Instance input, preserving source diagnostics and excluding inactive Parts.",
"started_at": "2026-08-01T02:37:05+0900",
"completed_at": "2026-08-01T02:44:01+0900"
},
{
"step": 3,
"name": "material-section-and-shear-defaults",
"status": "pending"
"status": "completed",
"summary": "Added complete-deck global material and active-ELSET Beam section mapping with exact property diagnostics, explicit shear-area conversion, and preserved Phase 1 default shear source.",
"started_at": "2026-08-01T02:44:01+0900",
"completed_at": "2026-08-01T02:59:40+0900"
},
{
"step": 4,
"name": "step-bc-load-and-noop-directives",
"status": "pending"
"status": "completed",
"summary": "Completed single static Step mapping with deterministic BC canonicalization, accumulated nodal loads, explicit no-op directives, exact source diagnostics, and supplied cantilever normalization.",
"started_at": "2026-08-01T02:59:40+0900",
"completed_at": "2026-08-01T03:23:25+0900"
}
]
],
"created_at": "2026-08-01T01:59:08+0900",
"completed_at": "2026-08-01T03:23:25+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+18 -5
View File
@@ -5,22 +5,35 @@
{
"step": 0,
"name": "comparison-metric-and-entity-matching",
"status": "pending"
"status": "completed",
"summary": "Added CSV-independent normalized comparison metrics and diagnostics with ResultFrame origin and Beam end-node matching.",
"started_at": "2026-08-02T02:42:01+0900",
"completed_at": "2026-08-02T03:22:08+0900"
},
{
"step": 1,
"name": "reference-csv-adapters",
"status": "pending"
"status": "completed",
"summary": "Added strict four-schema Abaqus reference CSV adapters with canonical row mapping, synthetic internal-force/stress fixtures, and validation diagnostics.",
"started_at": "2026-08-02T03:22:08+0900",
"completed_at": "2026-08-02T03:33:06+0900"
},
{
"step": 2,
"name": "cantilever-reference-comparison",
"status": "pending"
"status": "completed",
"summary": "Added the SCF=0 FESA cantilever projection, production HDF5/CSV correlation CLI, component-wise RMSE and Relative L2, Abaqus-to-FESA beam-force axis mapping, and displacement/reaction/internal-force reference coverage.",
"started_at": "2026-08-02T03:33:07+0900",
"completed_at": "2026-08-03T01:45:42+0900"
},
{
"step": 3,
"name": "qualification-report",
"status": "pending"
"status": "completed",
"summary": "Recorded docs/VALIDATION.md with Gate A PASS, Gate B EVALUABLE, 68/68 CTest and 20/20 Harness pytest evidence, all 18 correlation metrics, determinism coverage, and the remaining Abaqus stress qualification gap.",
"started_at": "2026-08-03T01:45:43+0900",
"completed_at": "2026-08-03T01:48:16+0900"
}
]
],
"created_at": "2026-08-02T02:42:01+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -0,0 +1,7 @@
{
"step": 2,
"name": "cantilever-reference-comparison",
"exitCode": 0,
"stdout": "Harness child verification was blocked by the managed WindowsApps PowerShell sandbox. The same acceptance contract was completed directly: CMake Debug build passed; CTest passed 68/68; Harness pytest passed 20/20; the fresh solver and reference-comparison CLI produced finite component metrics for 11 displacement rows, 11 reaction rows, and 20 internal-force rows.",
"stderr": ""
}
+46 -18
View File
@@ -6,45 +6,73 @@
- `/docs/PRD.md`
- `/docs/ARCHITECTURE.md`
- `/docs/ADR.md`
- `/docs/formulation/timoshenko-beam-3d.md`
- `/reference/cantilever beam/cantilever beam.inp`
- `/reference/cantilever beam/cantilever beam displacements.csv`
- `/reference/cantilever beam/cantilever beam reactions.csv`
- `/reference/cantilever beam/cantilever beam elemental forces.csv`
- `/include/fesa/analysis/run_solver.hpp`
- `/include/fesa/io/hdf5/reader.hpp`
- `/include/fesa/io/hdf5/writer.hpp`
- `/include/fesa/validation/comparison.hpp`
- `/include/fesa/validation/reference_csv.hpp`
## 작업
제공된 계층형 캔틸레버를 production pipeline으로 해석하고 현재 존재하는 변위와
반력만 Abaqus 2024 결과와 비교한다.
FESA 정식화 적합성과 Abaqus 결과 상관성을 별도 gate로 검증한다.
- `tests/reference/cantilever_reference_test.cpp`와 reference compare CLI를 먼저
작성한다.
- comparison request는 Instance `Part-1-1`, relative tolerance `1e-5`,
displacement absolute scale `1e-10`, reaction absolute scale `1e-8`을 명시한다.
- HDF5 결과와 CSV를 public adapter로 읽어 `(Instance,Node Label)`로 join한다.
- 요청하지 않은 internal force/stress 파일을 검색하거나 pass로 보고하지 않는다.
- equilibrium과 finite result도 함께 assertion한다.
- Gate A는 기존 해석해, energy, rigid mode 및 equilibrium 테스트를 그대로 엄격히
통과시킨다. `SCF=0.25`를 kernel에 추가하거나 production parser가 무시하게 하지
않는다.
- 원본 `cantilever beam.inp`와 세 CSV는 Abaqus provenance로 보존한다. 동일한
기하·재료·하중과 전단강성을 사용하되 `SCF=0`
`reference/cantilever beam/cantilever beam fesa.inp`를 추가해 production
pipeline으로 해석한다.
- `tests/unit/validation/comparison_test.cpp`에 component별 RMSE와 Relative L2의
실패 테스트를 먼저 추가한다. `include/fesa/validation/comparison.hpp`에는
`ComponentCorrelationMetric``CorrelationReport`,
`correlate_samples(std::span<const ComparisonSample>)`를 공개한다.
- Relative L2는 reference L2 norm과 component별 absolute-scale norm 중 큰 값을
분모로 사용한다. 서로 다른 component를 하나의 norm으로 합치지 않는다.
- `tests/unit/validation/reference_csv_test.cpp`의 internal-force fixture 기대값을
Abaqus `SF1,SF3,SF2,SM3,SM1,SM2`에서 FESA
`N,Vy,Vz,T,My,Mz` 순서로 재배열하도록 먼저 변경하고 RED를 확인한다.
- `tests/reference/cantilever_reference_test.cpp`와 reference compare CLI는 Instance
`PART-1_1-1`, 변위, 반력 및 요소 단면력 CSV 경로와 각 물리량의 absolute scale을
명시한다.
- HDF5 결과와 CSV를 public adapter로 읽어 nodal 결과는
`(Instance,Node Label)`, 요소 단면력은
`(Instance,Element Label,End Node Label)`로 join한다.
- correlation CLI는 요청된 결과가 모두 매칭되고 metric이 유한할 때 성공하며
component별 `count`, `rmse`, `relative_l2`를 출력한다. 관측값을 이용한 임의
pass/fail tolerance를 적용하지 않는다.
- equilibrium과 finite result도 함께 assertion한다. 요청하지 않은 stress 파일을
검색하거나 pass로 보고하지 않는다.
## Acceptance Criteria
```powershell
cmake --build --preset windows-debug
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.inp" --output out\cantilever-beam.h5
.\out\build\windows-debug\Debug\fesa-reference-compare.exe --results out\cantilever-beam.h5 --instance Part-1-1 --displacements "reference\cantilever beam\cantilever beam displacements.csv" --reactions "reference\cantilever beam\cantilever beam reactions.csv" --relative-tolerance 1e-5 --displacement-absolute-scale 1e-10 --reaction-absolute-scale 1e-8
.\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
ctest --preset windows-debug --output-on-failure
```
## 검증 절차
1. reference test가 실제 오차를 보고하며 실패하는 것을 확인한다.
2. discrepancy마다 가장 작은 analytical test를 추가한 뒤 근거 있는 kernel만 수정한다.
3. tolerance를 넓혀 결함을 숨기지 않는다.
4. 전체 테스트와 최대 정규화 오차를 index summary에 기록한다.
1. RMSE, Relative L2 및 요소력 축 재배열 테스트가 각각 의도한 이유로 실패하는
RED를 확인한다.
2. 최소 comparison metric과 CSV adapter 변경으로 GREEN을 만든다.
3. reference test가 production solve/HDF5/CSV/correlation 경로를 실행하고 세
물리량의 component metric을 출력하는지 확인한다.
4. 전체 테스트와 component별 metric 요약을 index summary에 기록한다.
## 금지사항
- reference `.inp` 또는 CSV를 수정하지 마라. 이유: 원본 golden을 보존해야 한다.
- 미제공 내력/응력 Abaqus 검증을 통과했다고 주장하지 마라. 이유: 증거가 없다.
- 원본 Abaqus `.inp` 또는 CSV를 수정하지 마라. 이유: 원본 golden과 formulation
provenance를 보존해야 한다.
- Abaqus `SCF=0.25`를 FESA가 구현하거나 무시하지 마라. 이유: 승인된 FESA
정식화와 입력 계약을 바꾼다.
- 미제공 응력 Abaqus 검증을 통과했다고 주장하지 마라. 이유: 증거가 없다.
- test-only parser/solver 경로를 만들지 마라. 이유: production pipeline 검증이다.
@@ -0,0 +1,7 @@
{
"step": 3,
"name": "qualification-report",
"exitCode": 0,
"stdout": "Created docs/VALIDATION.md from fresh evidence. CMake Debug build passed, CTest passed 68/68, Harness pytest passed 20/20, and the focused Beam3D2/ThreadCountDeterminism/Reference tests passed 4/4. The report records Gate A PASS, Gate B EVALUABLE, all 18 finite correlation metrics, and Abaqus stress as not yet qualified.",
"stderr": ""
}
+8 -7
View File
@@ -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를 약화한다.
@@ -5,17 +5,28 @@
{
"step": 0,
"name": "canonical-contribution-order",
"status": "pending"
"status": "completed",
"summary": "Added canonical MatrixContribution ordering and deterministic upper-triangle CSR merge with bitwise permutation, cancellation, signed-zero, and serial-oracle coverage.",
"started_at": "2026-08-01T22:17:55+0900",
"completed_at": "2026-08-01T23:06:21+0900"
},
{
"step": 1,
"name": "tbb-element-evaluation",
"status": "pending"
"status": "completed",
"summary": "Added AssemblyOptions and oneTBB Beam element evaluation with element-local contributions, canonical serial merge, explicit execution limits, and bitwise serial parity tests.",
"started_at": "2026-08-01T23:06:22+0900",
"completed_at": "2026-08-01T23:21:27+0900"
},
{
"step": 2,
"name": "thread-count-determinism",
"status": "pending"
"status": "completed",
"summary": "Added 10-run bitwise CSR, RHS, displacement, and reaction checks at 1, 2, and 16 threads plus a production serial/parallel assembly benchmark.",
"started_at": "2026-08-01T23:21:27+0900",
"completed_at": "2026-08-01T23:36:17+0900"
}
]
],
"created_at": "2026-08-01T22:17:55+0900",
"completed_at": "2026-08-01T23:36:18+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+15 -4
View File
@@ -5,17 +5,28 @@
{
"step": 0,
"name": "symmetric-csr-assembly",
"status": "pending"
"status": "completed",
"started_at": "2026-07-31T15:19:35+0900",
"summary": "Added upper-triangle symmetric CSR with every-row diagonals and deterministic serial Beam stiffness/load assembly with stable floating-point merge tests.",
"completed_at": "2026-07-31T15:45:02+0900"
},
{
"step": 1,
"name": "essential-bc-elimination",
"status": "pending"
"status": "completed",
"started_at": "2026-07-31T15:45:03+0900",
"summary": "Added DofManager-owned essential BC elimination with zero/nonzero RHS shifts, all-constrained handling, full reconstruction, and finite original-system reaction recovery.",
"completed_at": "2026-07-31T16:00:44+0900"
},
{
"step": 2,
"name": "pardiso-linear-solver",
"status": "pending"
"status": "completed",
"started_at": "2026-07-31T16:00:45+0900",
"summary": "Added a noncopyable LP64 MKL PARDISO SPD solver with explicit diagnostic-aware RAII release, CSR validation, residual diagnostics, and regression tests.",
"completed_at": "2026-07-31T16:16:28+0900"
}
]
],
"created_at": "2026-07-31T15:19:35+0900",
"completed_at": "2026-07-31T16:16:28+0900"
}
File diff suppressed because one or more lines are too long
+2 -1
View File
@@ -33,8 +33,9 @@ struct EquationSystem final {
- sparsity pattern builder와 numeric contribution merge를 분리한다.
- `(row,column,element-origin,local-order)`의 안정된 순서로 합산한다.
- 연결되지 않은 자유도를 포함해 모든 CSR row에 diagonal entry를 보존한다.
- hand-calculated 2-element system, duplicate contribution, external ID 순서 변화,
CSR invariant를 실패 테스트로 먼저 작성한다.
비결합적 부동소수점 합산 순서와 CSR invariant를 실패 테스트로 먼저 작성한다.
## Acceptance Criteria
File diff suppressed because one or more lines are too long
+7 -8
View File
@@ -9,30 +9,29 @@
- `/include/fesa/assembly/symmetric_csr.hpp`
- `/include/fesa/assembly/equation_system.hpp`
- `/include/fesa/fem/dof_manager.hpp`
- `/include/fesa/model/step_definition.hpp`
## 작업
0과 비영 지정변위를 지원하는 essential-BC elimination과 full-vector 복원을 구현한다.
`DofManager`가 소유한 equation mapping과 지정변위로 essential-BC elimination
수행하고, 기존 `DofManager::reconstruct_full()`로 full vector를 복원한다.
```cpp
struct ReducedSystem final {
SymmetricCsr stiffness;
std::vector<double> force;
std::vector<std::size_t> free_to_full;
std::vector<double> prescribed_full;
};
[[nodiscard]] ConstraintResult eliminate_essential_bcs(
const EquationSystem& original,
const DofManager& dofs,
std::span<const PrescribedDof> prescribed);
const DofManager& dofs);
[[nodiscard]] std::vector<double> recover_reaction(
const EquationSystem& original,
std::span<const double> full_displacement);
```
- 작은 hand calculation으로 RHS shift, 0/비영 prescribed value, all constrained,
충돌 조건, \(r=Ku-f\) 반력 복원을 먼저 테스트한다.
\(r=Ku-f\) 반력 복원을 먼저 테스트한다.
- 별도 prescribed 인수나 full/free mapping 상태를 중복하지 않는다.
- 입력과 산술 결과의 NaN/Inf를 성공 결과로 반환하지 않는다.
## Acceptance Criteria
@@ -46,7 +45,7 @@ ctest --preset windows-debug --output-on-failure
1. 비영 지정값 테스트의 실패를 먼저 확인한다.
2. 원래 EquationSystem을 보존한 채 reduced system을 생성한다.
3. 복원 변위와 원래 평형식 반력을 assertion한다.
3. `DofManager` 복원 변위와 원래 평형식 반력을 assertion한다.
4. 전체 테스트와 index를 갱신한다.
## 금지사항
File diff suppressed because one or more lines are too long
@@ -32,6 +32,8 @@ class PardisoLinearSolver final : public LinearSolver {
public:
PardisoLinearSolver();
~PardisoLinearSolver() override;
PardisoLinearSolver(const PardisoLinearSolver&) = delete;
PardisoLinearSolver& operator=(const PardisoLinearSolver&) = delete;
[[nodiscard]] LinearSolveResult solve(
const SymmetricCsr&,
std::span<const double>) override;
@@ -42,6 +44,7 @@ public:
테스트한다.
- `mtype=2`, LP64 index, `iparm[34]=1`, matrix checker, analysis/factor/solve/release
phase를 사용한다.
- analysis/factor/solve/release의 모든 MKL 오류를 solver diagnostic으로 변환한다.
## Acceptance Criteria
+19 -5
View File
@@ -5,22 +5,36 @@
{
"step": 0,
"name": "quadrature-and-shape-functions",
"status": "pending"
"status": "completed",
"summary": "Added fixed 1D Gauss rules, Line2 shape primitives, and a validated length/2 Jacobian with invariant and invalid-input tests.",
"started_at": "2026-07-31T01:37:23+0900",
"completed_at": "2026-07-31T01:48:32+0900"
},
{
"step": 1,
"name": "dof-manager",
"status": "pending"
"status": "completed",
"summary": "Added deterministic six-DOF equation numbering, prescribed-value reconstruction, and 12-DOF Beam element mapping with validation tests.",
"started_at": "2026-07-31T01:48:33+0900",
"completed_at": "2026-07-31T01:55:21+0900"
},
{
"step": 2,
"name": "beam-local-frame",
"status": "pending"
"status": "completed",
"summary": "Added scale-aware right-handed Beam frames and fixed 12x12 global-to-local transformations with degeneracy, orthonormality, rotation, and energy-invariance tests.",
"started_at": "2026-07-31T01:55:21+0900",
"completed_at": "2026-07-31T02:08:27+0900"
},
{
"step": 3,
"name": "timoshenko-stiffness-kernel",
"status": "pending"
"status": "completed",
"summary": "Documented and implemented selectively integrated Beam3D2 local/global stiffness with rigid-body, analytical, slenderness, and rotation-invariance tests.",
"started_at": "2026-07-31T02:08:27+0900",
"completed_at": "2026-07-31T02:24:47+0900"
}
]
],
"created_at": "2026-07-31T01:37:23+0900",
"completed_at": "2026-07-31T02:24:48+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+14 -7
View File
@@ -12,31 +12,38 @@
},
{
"dir": "fem-and-beam-kernel",
"status": "pending"
"status": "completed",
"completed_at": "2026-07-31T02:24:48+0900"
},
{
"dir": "equation-and-linear-solve",
"status": "pending"
"status": "completed",
"completed_at": "2026-07-31T16:16:28+0900"
},
{
"dir": "results-and-pipeline",
"status": "pending"
"status": "completed",
"completed_at": "2026-08-01T01:11:01+0900"
},
{
"dir": "abaqus-subset-completion",
"status": "pending"
"status": "completed",
"completed_at": "2026-08-01T03:23:25+0900"
},
{
"dir": "deterministic-parallel-assembly",
"status": "pending"
"status": "completed",
"completed_at": "2026-08-01T23:36:18+0900"
},
{
"dir": "result-contract-completion",
"status": "pending"
"status": "completed",
"completed_at": "2026-08-02T01:51:25+0900"
},
{
"dir": "beam-reference-qualification",
"status": "pending"
"status": "completed",
"completed_at": "2026-08-03T01:48:16+0900"
},
{
"dir": "internal-release",
+15 -4
View File
@@ -5,17 +5,28 @@
{
"step": 0,
"name": "beam-element-end-recovery",
"status": "pending"
"status": "completed",
"summary": "Added actual Beam node IDs and signed two-end strain, resultant, centroid/recovery-point stress recovery using the stiffness frame and reduced-shear convention, with hand-calculated tests.",
"started_at": "2026-08-02T00:27:47+0900",
"completed_at": "2026-08-02T01:14:40+0900"
},
{
"step": 1,
"name": "complete-result-contract",
"status": "pending"
"status": "completed",
"summary": "Added nodal and Beam element provenance, explicit field metadata and validation, and production element-end recovery orchestration in LinearStaticAnalysis.",
"started_at": "2026-08-02T01:14:43+0900",
"completed_at": "2026-08-02T01:28:37+0900"
},
{
"step": 2,
"name": "self-contained-hdf5",
"status": "pending"
"status": "completed",
"summary": "Defined schema 2.0.0 and completed self-contained HDF5 model, analysis, nodal/Beam result, diagnostic writer/reader round trips with strict metadata and version validation.",
"started_at": "2026-08-02T01:28:37+0900",
"completed_at": "2026-08-02T01:51:24+0900"
}
]
],
"created_at": "2026-08-02T00:27:47+0900",
"completed_at": "2026-08-02T01:51:25+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+19 -5
View File
@@ -5,22 +5,36 @@
{
"step": 0,
"name": "result-database",
"status": "pending"
"status": "completed",
"started_at": "2026-08-01T00:01:45+0900",
"summary": "Added HDF5-independent nodal result aggregates and results-stage validation for sizes, duplicate nodes/steps/frames, and finite values with focused tests.",
"completed_at": "2026-08-01T00:14:53+0900"
},
{
"step": 1,
"name": "minimal-hdf5-schema",
"status": "pending"
"status": "completed",
"started_at": "2026-08-01T00:14:54+0900",
"summary": "Documented HDF5 schema 1.0.0 and added a move-only RAII writer/public reader round trip for model identity, connectivity, shear provenance, and nodal results.",
"completed_at": "2026-08-01T00:51:34+0900"
},
{
"step": 2,
"name": "linear-static-analysis",
"status": "pending"
"status": "completed",
"started_at": "2026-08-01T00:51:34+0900",
"summary": "Added LinearStaticAnalysis orchestration through serial assembly, BC elimination, PARDISO with an all-constrained bypass, reaction recovery, and deterministic nodal results with equilibrium tests.",
"completed_at": "2026-08-01T01:01:15+0900"
},
{
"step": 3,
"name": "cli-pipeline-integration",
"status": "pending"
"status": "completed",
"started_at": "2026-08-01T01:01:16+0900",
"summary": "Added public run_solver orchestration and the solve CLI with source-preserving diagnostics, plus end-to-end HDF5 pipeline tests.",
"completed_at": "2026-08-01T01:11:00+0900"
}
]
],
"created_at": "2026-08-01T00:01:45+0900",
"completed_at": "2026-08-01T01:11:01+0900"
}
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -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
1 Part Instance Name Node Label U-U1 U-U2 U-U3 UR-UR1 UR-UR2 UR-UR3
2 PART-1_1-1 1 0 0.000000E+00 0 0.000000E+00 -1.00E-32 -1.000000E-30 0 0.000000E+00 1.00E-31 1.000000E-29 0 0.000000E+00
3 PART-1_1-1 2 0 0.000000E+00 0 0.000000E+00 -2.87E-06 -2.900000E-04 0 0.000000E+00 5.43E-06 5.428570E-04 0 0.000000E+00
4 PART-1_1-1 3 0 0.000000E+00 0 0.000000E+00 -1.09E-05 -1.094290E-03 0 0.000000E+00 1.03E-05 1.028570E-03 0 0.000000E+00
5 PART-1_1-1 4 0 0.000000E+00 0 0.000000E+00 -2.35E-05 -2.355720E-03 0 0.000000E+00 1.46E-05 1.457140E-03 0 0.000000E+00
6 PART-1_1-1 5 0 0.000000E+00 0 0.000000E+00 -4.00E-05 -4.017140E-03 0 0.000000E+00 1.83E-05 1.828570E-03 0 0.000000E+00
7 PART-1_1-1 6 0 0.000000E+00 0 0.000000E+00 -6.01E-05 -6.021430E-03 0 0.000000E+00 2.14E-05 2.142860E-03 0 0.000000E+00
8 PART-1_1-1 7 0 0.000000E+00 0 0.000000E+00 -8.29E-05 -8.311430E-03 0 0.000000E+00 2.40E-05 2.400000E-03 0 0.000000E+00
9 PART-1_1-1 8 0 0.000000E+00 0 0.000000E+00 -1.08E-04 -1.083000E-02 0 0.000000E+00 2.60E-05 2.600000E-03 0 0.000000E+00
10 PART-1_1-1 9 0 0.000000E+00 0 0.000000E+00 -1.35E-04 -1.352000E-02 0 0.000000E+00 2.74E-05 2.742860E-03 0 0.000000E+00
11 PART-1_1-1 10 0 0.000000E+00 0 0.000000E+00 -1.63E-04 -1.632430E-02 0 0.000000E+00 2.83E-05 2.828570E-03 0 0.000000E+00
12 PART-1_1-1 11 0 0.000000E+00 0 0.000000E+00 -1.92E-04 -1.918570E-02 0 0.000000E+00 2.86E-05 2.857140E-03 0 0.000000E+00
@@ -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,,,,,
1 Part Instance Name Element Label Node Label SF-SF1 SF-SF2 SF-SF3 SM-SM1 SM-SM2 SM-SM3
2 PART-1_1-1 1 1 0.000000E+00 -1.000000E+06 0.000000E+00 9.500000E+06 0.000000E+00 0.000000E+00
3 PART-1_1-1 1 2 0.000000E+00 -1.000000E+06 0.000000E+00 9.000000E+06 0.000000E+00 0.000000E+00
4 PART-1_1-1 2 2 0.000000E+00 -1.000000E+06 0.000000E+00 9.000000E+06 0.000000E+00 0.000000E+00
5 PART-1_1-1 2 3 0.000000E+00 -1.000000E+06 0.000000E+00 8.000000E+06 0.000000E+00 0.000000E+00
6 PART-1_1-1 3 3 0.000000E+00 -1.000000E+06 0.000000E+00 8.000000E+06 0.000000E+00 0.000000E+00
7 PART-1_1-1 3 4 0.000000E+00 -1.000000E+06 0.000000E+00 7.000000E+06 0.000000E+00 0.000000E+00
8 PART-1_1-1 4 4 0.000000E+00 -1.000000E+06 0.000000E+00 7.000000E+06 0.000000E+00 0.000000E+00
9 PART-1_1-1 4 5 0.000000E+00 -1.000000E+06 0.000000E+00 6.000000E+06 0.000000E+00 0.000000E+00
10 PART-1_1-1 5 5 0.000000E+00 -1.000000E+06 0.000000E+00 6.000000E+06 0.000000E+00 0.000000E+00
11 PART-1_1-1 5 6 0.000000E+00 -1.000000E+06 0.000000E+00 5.000000E+06 0.000000E+00 0.000000E+00
12 PART-1_1-1 6 6 0.000000E+00 -1.000000E+06 0.000000E+00 5.000000E+06 0.000000E+00 0.000000E+00
13 PART-1_1-1 6 7 0.000000E+00 -1.000000E+06 0.000000E+00 4.000000E+06 0.000000E+00 0.000000E+00
14 PART-1_1-1 7 7 0.000000E+00 -1.000000E+06 0.000000E+00 4.000000E+06 0.000000E+00 0.000000E+00
15 PART-1_1-1 7 8 0.000000E+00 -1.000000E+06 0.000000E+00 3.000000E+06 0.000000E+00 0.000000E+00
16 PART-1_1-1 8 8 0.000000E+00 -1.000000E+06 0.000000E+00 3.000000E+06 0.000000E+00 0.000000E+00
17 PART-1_1-1 8 9 0.000000E+00 -1.000000E+06 0.000000E+00 2.000000E+06 0.000000E+00 0.000000E+00
18 PART-1_1-1 9 9 0.000000E+00 -1.000000E+06 0.000000E+00 2.000000E+06 0.000000E+00 0.000000E+00
19 PART-1_1-1 9 10 0.000000E+00 -1.000000E+06 0.000000E+00 1.000000E+06 0.000000E+00 0.000000E+00
20 PART-1_1-1 10 10 0.000000E+00 -1.000000E+06 0.000000E+00 1.000000E+06 0.000000E+00 0.000000E+00
21 PART-1_1-1 10 11 0.000000E+00 -1.000000E+06 0.000000E+00 5.000000E+05 0.000000E+00 0.000000E+00
@@ -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
@@ -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
1 Part Instance Name Node Label RF-RF1 RF-RF2 RF-RF3 RM-RM1 RM-RM2 RM-RM3
2 PART-1_1-1 1 0 0.000000E+00 0 0.000000E+00 1.00E+04 1.000000E+06 0 0.000000E+00 -1.00E+05 -1.000000E+07 0 0.000000E+00
3 PART-1_1-1 2 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
4 PART-1_1-1 3 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
5 PART-1_1-1 4 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
6 PART-1_1-1 5 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
7 PART-1_1-1 6 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
8 PART-1_1-1 7 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
9 PART-1_1-1 8 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
10 PART-1_1-1 9 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
11 PART-1_1-1 10 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
12 PART-1_1-1 11 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00 0 0.000000E+00
+25 -22
View File
@@ -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
**
@@ -0,0 +1,226 @@
#include <fesa/analysis/linear_static_analysis.hpp>
#include <algorithm>
#include <array>
#include <cstddef>
#include <exception>
#include <stdexcept>
#include <optional>
#include <string>
#include <utility>
#include <vector>
#include <fesa/assembly/equation_system.hpp>
#include <fesa/assembly/serial_assembler.hpp>
#include <fesa/constraints/essential_bc.hpp>
#include <fesa/elements/beam/beam3d2.hpp>
#include <fesa/fem/dof_manager.hpp>
#include <fesa/solvers/linear/pardiso_linear_solver.hpp>
namespace fesa {
namespace {
AnalysisRunResult failure(std::vector<Diagnostic> diagnostics) {
return {false, std::nullopt, std::move(diagnostics)};
}
AnalysisRunResult equation_failure(
std::string code,
std::string message) {
return failure({{
DiagnosticStage::equation,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
}});
}
AnalysisRunResult results_failure(
std::string code,
std::string message) {
return failure({{
DiagnosticStage::results,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
}});
}
NodalFrame build_nodal_frame(
const Domain& domain,
const DofManager& dofs,
const std::vector<double>& displacement,
const std::vector<double>& reaction) {
std::vector<NodeId> node_ids;
node_ids.reserve(domain.nodes().size());
for (const Node& node : domain.nodes()) {
node_ids.push_back(node.id);
}
std::ranges::sort(node_ids);
NodalFrame nodal;
nodal.node_ids = std::move(node_ids);
nodal.origins.reserve(nodal.node_ids.size());
nodal.displacement.reserve(nodal.node_ids.size());
nodal.reaction.reserve(nodal.node_ids.size());
for (const NodeId node_id : nodal.node_ids) {
nodal.origins.push_back(domain.node(node_id).origin);
std::array<double, 6> node_displacement{};
std::array<double, 6> node_reaction{};
for (std::size_t component = 0; component < 6; ++component) {
const std::size_t full_dof = dofs.full_dof({
node_id,
static_cast<NodeDof>(component),
});
node_displacement[component] = displacement[full_dof];
node_reaction[component] = reaction[full_dof];
}
nodal.displacement.push_back(node_displacement);
nodal.reaction.push_back(node_reaction);
}
return nodal;
}
ElementFrame build_element_frame(
const Domain& domain,
const DofManager& dofs,
const std::vector<double>& displacement) {
std::vector<const BeamElement*> elements;
elements.reserve(domain.beam_elements().size());
for (const BeamElement& element : domain.beam_elements()) {
elements.push_back(&element);
}
std::ranges::sort(
elements,
{},
[](const BeamElement* element) {
return element->id.value();
});
ElementFrame frame;
frame.beams.reserve(elements.size());
for (const BeamElement* element : elements) {
const Beam3D2Input input{
{
domain.node(element->nodes[0]).position,
domain.node(element->nodes[1]).position,
},
element->nodes,
domain.material(element->material),
domain.section(element->section),
};
const BeamKernelResult kernel = compute_beam3d2(input);
if (!kernel.contribution.has_value()) {
const std::string message = kernel.diagnostics.empty()
? "Beam recovery requires a valid element input."
: kernel.diagnostics.front().message;
throw std::runtime_error{message};
}
std::array<double, 12> element_displacement{};
const std::array<std::size_t, 12> full_dofs =
dofs.element_full_dofs(*element);
for (std::size_t local = 0; local < full_dofs.size(); ++local) {
element_displacement[local] = displacement[full_dofs[local]];
}
std::vector<BeamSectionResult> recovered = recover_beam3d2(
input,
element_displacement,
input.section.recovery_points);
if (recovered.size() != 2U) {
throw std::logic_error{
"Beam recovery must return exactly two end results."};
}
frame.beams.push_back({
element->id,
element->origin,
kernel.contribution->frame,
{
std::move(recovered[0]),
std::move(recovered[1]),
},
});
}
return frame;
}
} // namespace
AnalysisRunResult LinearStaticAnalysis::run(const Domain& domain) const {
const DofManager dofs = DofManager::build(domain);
std::optional<EquationSystem> original;
try {
original = assemble_serial(domain, dofs);
} catch (const std::exception& error) {
return equation_failure(
"equation.assembly_failed", error.what());
}
ConstraintResult constrained =
eliminate_essential_bcs(*original, dofs);
if (!constrained.reduced_system.has_value()) {
return failure(std::move(constrained.diagnostics));
}
std::vector<double> reduced_solution;
if (dofs.free_equation_count() != 0) {
PardisoLinearSolver solver;
LinearSolveResult solved = solver.solve(
constrained.reduced_system->stiffness,
constrained.reduced_system->force);
if (!solved.diagnostics.empty()) {
return failure(std::move(solved.diagnostics));
}
reduced_solution = std::move(solved.solution);
}
std::vector<double> displacement;
try {
displacement = dofs.reconstruct_full(reduced_solution);
} catch (const std::exception& error) {
return equation_failure(
"equation.solution_reconstruction_failed", error.what());
}
std::vector<double> reaction;
try {
reaction = recover_reaction(*original, displacement);
} catch (const std::exception& error) {
return equation_failure(
"equation.reaction_recovery_failed", error.what());
}
ElementFrame element;
try {
element = build_element_frame(domain, dofs, displacement);
} catch (const std::exception& error) {
return results_failure(
"results.element_recovery_failed", error.what());
}
ResultDatabase database{
"2.0.0",
{{
domain.step().name,
{{
1.0,
build_nodal_frame(
domain, dofs, displacement, reaction),
std::move(element),
{},
}},
}},
};
Status status = validate_result_database(database);
if (!status.succeeded) {
return failure(std::move(status.diagnostics));
}
return {true, std::move(database), {}};
}
} // namespace fesa
+51
View File
@@ -0,0 +1,51 @@
#include <fesa/analysis/run_solver.hpp>
#include <optional>
#include <string>
#include <utility>
#include <fesa/io/abaqus/parser.hpp>
#include <fesa/io/abaqus/semantic_mapper.hpp>
#include <fesa/io/hdf5/writer.hpp>
namespace fesa {
namespace {
std::string path_utf8(const std::filesystem::path& path) {
const std::u8string value = path.u8string();
return {reinterpret_cast<const char*>(value.data()), value.size()};
}
} // namespace
AnalysisRunResult run_solver(const AnalysisRequest& request) {
ParseDeckResult parsed = parse_deck(request.input_path);
if (!parsed.deck.has_value()) {
return {false, std::nullopt, std::move(parsed.diagnostics)};
}
DomainBuildResult mapped = map_deck_to_domain(*parsed.deck);
if (!mapped.domain.has_value()) {
return {false, std::nullopt, std::move(mapped.diagnostics)};
}
const Hdf5InputIdentity identity{
path_utf8(request.input_path),
parsed.input_fingerprint,
};
AnalysisRunResult run = LinearStaticAnalysis{}.run(*mapped.domain);
if (!run.succeeded || !run.results.has_value()) {
return run;
}
std::vector<Diagnostic> write_diagnostics = write_hdf5(
request.output_path, *mapped.domain, *run.results, identity);
if (!write_diagnostics.empty()) {
return {false, std::nullopt, std::move(write_diagnostics)};
}
return run;
}
} // namespace fesa
+132
View File
@@ -0,0 +1,132 @@
#include <fesa/assembly/contribution.hpp>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <limits>
#include <stdexcept>
#include <tuple>
#include <utility>
#include <vector>
namespace fesa {
namespace {
std::int32_t csr_index(const std::size_t value) {
if (value >
static_cast<std::size_t>(
std::numeric_limits<std::int32_t>::max())) {
throw std::overflow_error{
"Symmetric CSR exceeds the 32-bit index range."};
}
return static_cast<std::int32_t>(value);
}
bool contribution_less(
const MatrixContribution& left,
const MatrixContribution& right) {
return std::tuple{
left.row,
left.column,
left.element,
left.local_order} <
std::tuple{
right.row,
right.column,
right.element,
right.local_order};
}
} // namespace
std::vector<MatrixContribution> canonicalize_contributions(
const std::span<const MatrixContribution> contributions) {
std::vector<MatrixContribution> canonical{
contributions.begin(), contributions.end()};
std::ranges::sort(canonical, contribution_less);
return canonical;
}
SymmetricCsr merge_contributions(
const std::size_t order,
const std::span<const MatrixContribution> canonical) {
if (!std::ranges::is_sorted(canonical, contribution_less)) {
throw std::invalid_argument{
"Matrix contributions are not in canonical order."};
}
using Coordinate = std::pair<std::size_t, std::size_t>;
std::vector<Coordinate> coordinates;
coordinates.reserve(order + canonical.size());
for (std::size_t row = 0; row < order; ++row) {
coordinates.emplace_back(row, row);
}
for (const MatrixContribution& contribution : canonical) {
if (contribution.row > contribution.column ||
contribution.column >= order) {
throw std::invalid_argument{
"Matrix contribution is outside the upper triangle."};
}
coordinates.emplace_back(
contribution.row, contribution.column);
}
std::ranges::sort(coordinates);
coordinates.erase(
std::ranges::unique(coordinates).begin(),
coordinates.end());
SymmetricCsr matrix{
order,
std::vector<std::int32_t>(order + 1, 0),
{},
{},
};
matrix.column_indices.reserve(coordinates.size());
for (const auto [row, column] : coordinates) {
++matrix.row_offsets[row + 1];
matrix.column_indices.push_back(csr_index(column));
}
for (std::size_t row = 0; row < order; ++row) {
const std::size_t offset =
static_cast<std::size_t>(matrix.row_offsets[row]) +
static_cast<std::size_t>(matrix.row_offsets[row + 1]);
matrix.row_offsets[row + 1] = csr_index(offset);
}
matrix.values.resize(matrix.column_indices.size(), 0.0);
std::size_t contribution_index = 0;
while (contribution_index < canonical.size()) {
const MatrixContribution& first =
canonical[contribution_index];
double value = 0.0;
do {
value += canonical[contribution_index].value;
++contribution_index;
} while (
contribution_index < canonical.size() &&
canonical[contribution_index].row == first.row &&
canonical[contribution_index].column == first.column);
const auto row_begin =
matrix.column_indices.begin() + matrix.row_offsets[first.row];
const auto row_end =
matrix.column_indices.begin() +
matrix.row_offsets[first.row + 1];
const auto entry = std::lower_bound(
row_begin,
row_end,
csr_index(first.column));
if (entry == row_end || *entry != csr_index(first.column)) {
throw std::logic_error{
"Numeric contribution is absent from the CSR pattern."};
}
matrix.values[static_cast<std::size_t>(
std::distance(matrix.column_indices.begin(), entry))] = value;
}
return matrix;
}
} // namespace fesa
+209
View File
@@ -0,0 +1,209 @@
#include <fesa/assembly/assembler.hpp>
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <limits>
#include <numeric>
#include <optional>
#include <stdexcept>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
#include <fesa/assembly/contribution.hpp>
#include <fesa/elements/beam/beam3d2.hpp>
#include <oneapi/tbb/blocked_range.h>
#include <oneapi/tbb/parallel_for.h>
#include <oneapi/tbb/task_arena.h>
namespace fesa {
namespace {
struct ElementEvaluation final {
std::vector<MatrixContribution> contributions;
std::optional<std::string> error;
};
auto origin_key(const EntityOrigin& origin) {
return std::tie(
origin.instance_name,
origin.local_label,
origin.part_name);
}
std::string kernel_error_message(
const BeamElement& element,
const BeamKernelResult& result) {
std::string message =
"Beam element " + std::to_string(element.origin.local_label) +
" kernel failed";
for (const Diagnostic& diagnostic : result.diagnostics) {
message += ": " + diagnostic.code + " - " + diagnostic.message;
}
return message;
}
std::vector<std::size_t> canonical_element_order(const Domain& domain) {
std::vector<std::size_t> order(domain.beam_elements().size());
std::iota(order.begin(), order.end(), std::size_t{0});
std::ranges::sort(
order,
[&domain](const std::size_t left, const std::size_t right) {
return origin_key(domain.beam_elements()[left].origin) <
origin_key(domain.beam_elements()[right].origin);
});
return order;
}
std::vector<ElementId> canonical_contribution_element_ids(
const std::vector<std::size_t>& order) {
// Domain ElementIds are not ordered by input identity. These tie-break
// IDs encode the existing serial assembler's element-origin order.
std::vector<ElementId> ids(order.size(), ElementId{0});
for (std::size_t rank = 0; rank < order.size(); ++rank) {
ids[order[rank]] = ElementId{static_cast<std::int64_t>(rank)};
}
return ids;
}
ElementEvaluation evaluate_element(
const Domain& domain,
const DofManager& dofs,
const BeamElement& element,
const ElementId canonical_id) {
const BeamKernelResult result = compute_beam3d2({
{
domain.node(element.nodes[0]).position,
domain.node(element.nodes[1]).position,
},
element.nodes,
domain.material(element.material),
domain.section(element.section),
});
if (!result.contribution.has_value()) {
return {{}, kernel_error_message(element, result)};
}
ElementEvaluation evaluation;
evaluation.contributions.reserve(78);
const std::array<std::size_t, 12> full_dofs =
dofs.element_full_dofs(element);
for (std::size_t local_row = 0; local_row < full_dofs.size();
++local_row) {
for (std::size_t local_column = local_row;
local_column < full_dofs.size();
++local_column) {
evaluation.contributions.push_back({
std::min(
full_dofs[local_row],
full_dofs[local_column]),
std::max(
full_dofs[local_row],
full_dofs[local_column]),
canonical_id,
static_cast<std::uint16_t>(
local_row * full_dofs.size() + local_column),
result.contribution
->global_stiffness[local_row][local_column],
});
}
}
return evaluation;
}
std::vector<double> assemble_force(
const Domain& domain,
const DofManager& dofs) {
std::vector<double> force(dofs.full_dof_count(), 0.0);
for (const NodalLoad& load : domain.step().nodal_loads) {
for (std::size_t component = 0; component < load.values.size();
++component) {
const auto dof = static_cast<NodeDof>(component);
force[dofs.full_dof({load.node, dof})] +=
load.values[component];
}
}
return force;
}
void validate_options(const AssemblyOptions options) {
if (options.max_threads == 0) {
throw std::invalid_argument{
"Assembly max_threads must be greater than zero."};
}
if (options.max_threads >
static_cast<std::size_t>(std::numeric_limits<int>::max())) {
throw std::invalid_argument{
"Assembly max_threads exceeds the TBB task arena range."};
}
if (options.grain_size == 0) {
throw std::invalid_argument{
"Assembly grain_size must be greater than zero."};
}
}
} // namespace
EquationSystem assemble_parallel(
const Domain& domain,
const DofManager& dofs,
const AssemblyOptions options) {
validate_options(options);
const std::vector<std::size_t> element_order =
canonical_element_order(domain);
const std::vector<ElementId> element_ids =
canonical_contribution_element_ids(element_order);
std::vector<ElementEvaluation> evaluations(
domain.beam_elements().size());
oneapi::tbb::task_arena arena{
static_cast<int>(options.max_threads)};
arena.execute([&] {
oneapi::tbb::parallel_for(
oneapi::tbb::blocked_range<std::size_t>{
0,
evaluations.size(),
options.grain_size,
},
[&](const oneapi::tbb::blocked_range<std::size_t>& range) {
for (std::size_t index = range.begin();
index != range.end();
++index) {
ElementEvaluation local = evaluate_element(
domain,
dofs,
domain.beam_elements()[index],
element_ids[index]);
evaluations[index] = std::move(local);
}
});
});
std::vector<MatrixContribution> contributions;
contributions.reserve(domain.beam_elements().size() * 78);
for (const std::size_t element_index : element_order) {
ElementEvaluation& evaluation = evaluations[element_index];
if (evaluation.error.has_value()) {
throw std::runtime_error{std::move(*evaluation.error)};
}
contributions.insert(
contributions.end(),
std::make_move_iterator(evaluation.contributions.begin()),
std::make_move_iterator(evaluation.contributions.end()));
}
const std::vector<MatrixContribution> canonical =
canonicalize_contributions(contributions);
return {
merge_contributions(dofs.full_dof_count(), canonical),
assemble_force(domain, dofs),
};
}
} // namespace fesa
+245
View File
@@ -0,0 +1,245 @@
#include <fesa/assembly/serial_assembler.hpp>
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <limits>
#include <stdexcept>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
#include <fesa/elements/beam/beam3d2.hpp>
namespace fesa {
namespace {
struct NumericContribution final {
std::size_t row;
std::size_t column;
const EntityOrigin* element_origin;
std::size_t local_order;
double value;
};
std::int32_t csr_index(const std::size_t value) {
if (value >
static_cast<std::size_t>(
std::numeric_limits<std::int32_t>::max())) {
throw std::overflow_error{
"Symmetric CSR exceeds the 32-bit index range."};
}
return static_cast<std::int32_t>(value);
}
auto origin_key(const EntityOrigin& origin) {
return std::tie(
origin.instance_name,
origin.local_label,
origin.part_name);
}
SymmetricCsr build_sparsity_pattern(
const Domain& domain,
const DofManager& dofs) {
using Coordinate = std::pair<std::size_t, std::size_t>;
std::vector<Coordinate> coordinates;
coordinates.reserve(
dofs.full_dof_count() +
domain.beam_elements().size() * 78);
for (std::size_t full = 0; full < dofs.full_dof_count(); ++full) {
coordinates.emplace_back(full, full);
}
for (const BeamElement& element : domain.beam_elements()) {
const std::array<std::size_t, 12> full_dofs =
dofs.element_full_dofs(element);
for (std::size_t local_row = 0; local_row < full_dofs.size();
++local_row) {
for (std::size_t local_column = local_row;
local_column < full_dofs.size();
++local_column) {
coordinates.emplace_back(
std::min(
full_dofs[local_row],
full_dofs[local_column]),
std::max(
full_dofs[local_row],
full_dofs[local_column]));
}
}
}
std::ranges::sort(coordinates);
coordinates.erase(
std::ranges::unique(coordinates).begin(),
coordinates.end());
SymmetricCsr pattern{
dofs.full_dof_count(),
std::vector<std::int32_t>(dofs.full_dof_count() + 1, 0),
{},
{},
};
pattern.column_indices.reserve(coordinates.size());
for (const auto [row, column] : coordinates) {
if (row >= pattern.order || column >= pattern.order) {
throw std::invalid_argument{
"DofManager element mapping exceeds the full system order."};
}
++pattern.row_offsets[row + 1];
pattern.column_indices.push_back(csr_index(column));
}
for (std::size_t row = 0; row < pattern.order; ++row) {
const std::size_t offset =
static_cast<std::size_t>(pattern.row_offsets[row]) +
static_cast<std::size_t>(pattern.row_offsets[row + 1]);
pattern.row_offsets[row + 1] = csr_index(offset);
}
pattern.values.resize(pattern.column_indices.size(), 0.0);
return pattern;
}
std::runtime_error kernel_error(
const BeamElement& element,
const BeamKernelResult& result) {
std::string message =
"Beam element " + std::to_string(element.origin.local_label) +
" kernel failed";
for (const Diagnostic& diagnostic : result.diagnostics) {
message += ": " + diagnostic.code + " - " + diagnostic.message;
}
return std::runtime_error{std::move(message)};
}
std::vector<NumericContribution> collect_numeric_contributions(
const Domain& domain,
const DofManager& dofs) {
std::vector<NumericContribution> contributions;
contributions.reserve(domain.beam_elements().size() * 78);
for (const BeamElement& element : domain.beam_elements()) {
const BeamKernelResult result = compute_beam3d2({
{
domain.node(element.nodes[0]).position,
domain.node(element.nodes[1]).position,
},
element.nodes,
domain.material(element.material),
domain.section(element.section),
});
if (!result.contribution.has_value()) {
throw kernel_error(element, result);
}
const std::array<std::size_t, 12> full_dofs =
dofs.element_full_dofs(element);
for (std::size_t local_row = 0; local_row < full_dofs.size();
++local_row) {
for (std::size_t local_column = local_row;
local_column < full_dofs.size();
++local_column) {
contributions.push_back({
std::min(
full_dofs[local_row],
full_dofs[local_column]),
std::max(
full_dofs[local_row],
full_dofs[local_column]),
&element.origin,
local_row * full_dofs.size() + local_column,
result.contribution
->global_stiffness[local_row][local_column],
});
}
}
}
std::ranges::sort(
contributions,
[](const NumericContribution& left,
const NumericContribution& right) {
return std::tuple{
left.row,
left.column,
origin_key(*left.element_origin),
left.local_order} <
std::tuple{
right.row,
right.column,
origin_key(*right.element_origin),
right.local_order};
});
return contributions;
}
void merge_numeric_contributions(
SymmetricCsr& matrix,
const std::vector<NumericContribution>& contributions) {
std::size_t contribution_index = 0;
while (contribution_index < contributions.size()) {
const NumericContribution& first =
contributions[contribution_index];
double value = 0.0;
do {
value += contributions[contribution_index].value;
++contribution_index;
} while (
contribution_index < contributions.size() &&
contributions[contribution_index].row == first.row &&
contributions[contribution_index].column == first.column);
const auto row_begin =
matrix.column_indices.begin() +
matrix.row_offsets[first.row];
const auto row_end =
matrix.column_indices.begin() +
matrix.row_offsets[first.row + 1];
const auto entry = std::lower_bound(
row_begin,
row_end,
csr_index(first.column));
if (entry == row_end || *entry != csr_index(first.column)) {
throw std::logic_error{
"Numeric contribution is absent from the CSR pattern."};
}
matrix.values[static_cast<std::size_t>(
std::distance(matrix.column_indices.begin(), entry))] = value;
}
}
std::vector<double> assemble_force(
const Domain& domain,
const DofManager& dofs) {
std::vector<double> force(dofs.full_dof_count(), 0.0);
for (const NodalLoad& load : domain.step().nodal_loads) {
for (std::size_t component = 0; component < load.values.size();
++component) {
const auto dof = static_cast<NodeDof>(component);
force[dofs.full_dof({load.node, dof})] +=
load.values[component];
}
}
return force;
}
} // namespace
EquationSystem assemble_serial(
const Domain& domain,
const DofManager& dofs) {
SymmetricCsr stiffness = build_sparsity_pattern(domain, dofs);
const std::vector<NumericContribution> contributions =
collect_numeric_contributions(domain, dofs);
merge_numeric_contributions(stiffness, contributions);
return {
std::move(stiffness),
assemble_force(domain, dofs),
};
}
} // namespace fesa
+62 -1
View File
@@ -1,14 +1,75 @@
#include <fesa/analysis/run_solver.hpp>
#include <fesa/core/version.hpp>
#include <filesystem>
#include <iostream>
#include <string>
#include <string_view>
namespace {
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_usage() {
std::cerr << "Usage:\n"
<< " fesa solve <model.inp> --output <results.h5>\n"
<< " fesa --version\n";
}
} // namespace
int main(int argc, char* argv[]) {
if (argc == 2 && std::string_view{argv[1]} == "--version") {
std::cout << fesa::version() << '\n';
return 0;
}
std::cerr << "Usage: fesa --version\n";
if (argc == 5 && std::string_view{argv[1]} == "solve" &&
std::string_view{argv[3]} == "--output") {
const fesa::AnalysisRunResult result = fesa::run_solver({
std::filesystem::path{argv[2]},
std::filesystem::path{argv[4]},
});
if (result.succeeded) {
return 0;
}
for (const fesa::Diagnostic& diagnostic : result.diagnostics) {
print_diagnostic(diagnostic);
}
return 1;
}
print_usage();
return 1;
}
+254
View File
@@ -0,0 +1,254 @@
#include <fesa/constraints/essential_bc.hpp>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <optional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace fesa {
namespace {
std::optional<std::string> validate_system(
const EquationSystem& system) {
const SymmetricCsr& matrix = system.stiffness;
if (system.force.size() != matrix.order) {
return "Force size must equal the matrix order.";
}
if (matrix.order >
static_cast<std::size_t>(
std::numeric_limits<std::int32_t>::max())) {
return "Matrix order exceeds the 32-bit CSR index range.";
}
if (matrix.row_offsets.size() != matrix.order + 1 ||
matrix.row_offsets.empty() ||
matrix.row_offsets.front() != 0) {
return "CSR row offsets must contain order + 1 entries "
"starting at zero.";
}
if (matrix.column_indices.size() != matrix.values.size()) {
return "CSR column and value counts must match.";
}
if (!std::ranges::all_of(
matrix.values, [](const double value) {
return std::isfinite(value);
}) ||
!std::ranges::all_of(
system.force, [](const double value) {
return std::isfinite(value);
})) {
return "Matrix and force values must be finite.";
}
std::int32_t previous_offset = 0;
for (const std::int32_t offset : matrix.row_offsets) {
if (offset < previous_offset || offset < 0 ||
static_cast<std::size_t>(offset) >
matrix.column_indices.size()) {
return "CSR row offsets must be nondecreasing and in range.";
}
previous_offset = offset;
}
if (static_cast<std::size_t>(matrix.row_offsets.back()) !=
matrix.column_indices.size()) {
return "The final CSR row offset must equal the entry count.";
}
for (std::size_t row = 0; row < matrix.order; ++row) {
std::int32_t previous_column = -1;
const std::size_t begin =
static_cast<std::size_t>(matrix.row_offsets[row]);
const std::size_t end =
static_cast<std::size_t>(matrix.row_offsets[row + 1]);
for (std::size_t entry = begin; entry < end; ++entry) {
const std::int32_t column = matrix.column_indices[entry];
if (column < static_cast<std::int32_t>(row) ||
column >= static_cast<std::int32_t>(matrix.order) ||
column <= previous_column) {
return "CSR rows must contain sorted unique "
"upper-triangle columns.";
}
previous_column = column;
}
}
return std::nullopt;
}
Diagnostic equation_error(std::string code, std::string message) {
return {
DiagnosticStage::equation,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
};
}
std::int32_t csr_index(const std::size_t value) {
if (value >
static_cast<std::size_t>(
std::numeric_limits<std::int32_t>::max())) {
throw std::overflow_error{
"Reduced system exceeds the 32-bit CSR index range."};
}
return static_cast<std::int32_t>(value);
}
} // namespace
ConstraintResult eliminate_essential_bcs(
const EquationSystem& original,
const DofManager& dofs) {
ConstraintResult result;
if (const auto error = validate_system(original);
error.has_value()) {
result.diagnostics.push_back(equation_error(
"equation.invalid_system", *error));
return result;
}
if (original.stiffness.order != dofs.full_dof_count()) {
result.diagnostics.push_back(equation_error(
"equation.dof_count_mismatch",
"Matrix order must equal the DofManager full DOF count."));
return result;
}
const std::size_t full_count = original.stiffness.order;
std::vector<std::size_t> full_to_free(full_count, full_count);
std::vector<double> prescribed_full(full_count, 0.0);
for (std::size_t full = 0; full < full_count; ++full) {
const std::optional<std::size_t> equation =
dofs.equation(full);
if (equation.has_value()) {
full_to_free[full] = *equation;
} else {
prescribed_full[full] =
dofs.prescribed_value(full).value();
}
}
std::vector<double> reduced_force(
dofs.free_equation_count(), 0.0);
for (std::size_t full = 0; full < full_count; ++full) {
if (full_to_free[full] != full_count) {
reduced_force[full_to_free[full]] = original.force[full];
}
}
SymmetricCsr reduced{
dofs.free_equation_count(),
std::vector<std::int32_t>(
dofs.free_equation_count() + 1, 0),
{},
{},
};
reduced.column_indices.reserve(
original.stiffness.column_indices.size());
reduced.values.reserve(original.stiffness.values.size());
for (std::size_t full_row = 0; full_row < full_count; ++full_row) {
const std::size_t begin = static_cast<std::size_t>(
original.stiffness.row_offsets[full_row]);
const std::size_t end = static_cast<std::size_t>(
original.stiffness.row_offsets[full_row + 1]);
for (std::size_t entry = begin; entry < end; ++entry) {
const std::size_t full_column =
static_cast<std::size_t>(
original.stiffness.column_indices[entry]);
const double stiffness = original.stiffness.values[entry];
if (full_to_free[full_row] == full_count) {
if (full_to_free[full_column] != full_count) {
reduced_force[full_to_free[full_column]] -=
stiffness * prescribed_full[full_row];
}
} else if (full_to_free[full_column] == full_count) {
reduced_force[full_to_free[full_row]] -=
stiffness * prescribed_full[full_column];
} else {
reduced.column_indices.push_back(
csr_index(full_to_free[full_column]));
reduced.values.push_back(stiffness);
}
}
if (full_to_free[full_row] != full_count) {
reduced.row_offsets[full_to_free[full_row] + 1] =
csr_index(reduced.column_indices.size());
}
}
if (!std::ranges::all_of(
reduced_force, [](const double value) {
return std::isfinite(value);
})) {
result.diagnostics.push_back(equation_error(
"equation.nonfinite_result",
"Essential BC elimination produced a nonfinite force."));
return result;
}
result.reduced_system = ReducedSystem{
std::move(reduced),
std::move(reduced_force),
};
return result;
}
std::vector<double> recover_reaction(
const EquationSystem& original,
const std::span<const double> full_displacement) {
if (const auto error = validate_system(original);
error.has_value()) {
throw std::invalid_argument{*error};
}
if (full_displacement.size() != original.stiffness.order) {
throw std::invalid_argument{
"Full displacement size must equal the matrix order."};
}
if (!std::ranges::all_of(
full_displacement, [](const double value) {
return std::isfinite(value);
})) {
throw std::invalid_argument{
"Full displacement values must be finite."};
}
std::vector<double> reaction = original.force;
for (double& value : reaction) {
value = -value;
}
for (std::size_t row = 0;
row < original.stiffness.order;
++row) {
const std::size_t begin = static_cast<std::size_t>(
original.stiffness.row_offsets[row]);
const std::size_t end = static_cast<std::size_t>(
original.stiffness.row_offsets[row + 1]);
for (std::size_t entry = begin; entry < end; ++entry) {
const std::size_t column =
static_cast<std::size_t>(
original.stiffness.column_indices[entry]);
const double stiffness = original.stiffness.values[entry];
reaction[row] +=
stiffness * full_displacement[column];
if (column != row) {
reaction[column] +=
stiffness * full_displacement[row];
}
}
}
if (!std::ranges::all_of(
reaction, [](const double value) {
return std::isfinite(value);
})) {
throw std::invalid_argument{
"Reaction recovery produced a nonfinite value."};
}
return reaction;
}
} // namespace fesa
+349
View File
@@ -0,0 +1,349 @@
#include <fesa/elements/beam/beam3d2.hpp>
#include <array>
#include <cmath>
#include <cstddef>
#include <optional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
#include <fesa/fem/gauss_rule.hpp>
#include <fesa/fem/line2_shape.hpp>
namespace fesa {
namespace {
using StrainMatrix =
std::array<std::array<double, 12>, 6>;
BeamKernelResult error_result(std::string code, std::string message) {
BeamKernelResult result;
result.diagnostics.push_back({
DiagnosticStage::model,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
});
return result;
}
bool is_positive_finite(const double value) {
return std::isfinite(value) && value > 0.0;
}
std::array<double, 6> constitutive_values(
const Beam3D2Input& input) {
const double shear_modulus =
input.material.young / (2.0 * (1.0 + input.material.poisson));
return {
input.material.young * input.section.area,
shear_modulus * input.section.shear_area_y,
shear_modulus * input.section.shear_area_z,
shear_modulus * input.section.torsion_j,
input.material.young * input.section.iy,
input.material.young * input.section.iz,
};
}
std::optional<BeamKernelResult> validate_properties(
const Beam3D2Input& input) {
if (!std::isfinite(input.material.young) ||
!std::isfinite(input.material.poisson)) {
return error_result(
"model.nonfinite_value",
"Beam kernel requires finite elastic constants.");
}
if (input.material.young <= 0.0 ||
input.material.poisson <= -1.0 ||
input.material.poisson >= 0.5) {
return error_result(
"model.invalid_material",
"Beam kernel requires E > 0 and -1 < nu < 0.5.");
}
if (!is_positive_finite(input.section.area) ||
!is_positive_finite(input.section.iy) ||
!is_positive_finite(input.section.iz) ||
!is_positive_finite(input.section.torsion_j) ||
!is_positive_finite(input.section.shear_area_y) ||
!is_positive_finite(input.section.shear_area_z)) {
return error_result(
"model.invalid_section",
"Beam kernel requires positive finite section properties.");
}
return std::nullopt;
}
StrainMatrix strain_matrix(
const double xi,
const double jacobian) {
StrainMatrix strain{};
const auto shape = line2_shape(xi);
const auto natural_derivative = line2_shape_derivative(xi);
for (std::size_t node = 0; node < shape.size(); ++node) {
const std::size_t offset = node * 6;
const double derivative = natural_derivative[node] / jacobian;
strain[0][offset] = derivative;
strain[1][offset + 1] = derivative;
strain[1][offset + 5] = -shape[node];
strain[2][offset + 2] = derivative;
strain[2][offset + 4] = shape[node];
strain[3][offset + 3] = derivative;
strain[4][offset + 4] = derivative;
strain[5][offset + 5] = derivative;
}
return strain;
}
template <std::size_t ComponentCount>
void integrate_components(
Matrix12& stiffness,
const std::span<const GaussPoint1D> rule,
const std::array<std::size_t, ComponentCount>& components,
const std::array<double, 6>& constitutive,
const double jacobian) {
for (const GaussPoint1D& point : rule) {
const StrainMatrix strain = strain_matrix(point.xi, jacobian);
const double integration_weight = jacobian * point.weight;
for (std::size_t row = 0; row < stiffness.size(); ++row) {
for (std::size_t column = row;
column < stiffness[row].size();
++column) {
double entry = 0.0;
for (const std::size_t component : components) {
entry +=
strain[component][row] *
constitutive[component] *
strain[component][column];
}
stiffness[row][column] +=
entry * integration_weight;
}
}
}
}
void mirror_upper_triangle(Matrix12& matrix) {
for (std::size_t row = 0; row < matrix.size(); ++row) {
for (std::size_t column = row + 1;
column < matrix[row].size();
++column) {
matrix[column][row] = matrix[row][column];
}
}
}
Matrix12 transform_stiffness(
const Matrix12& local,
const Matrix12& transformation) {
Matrix12 global{};
for (std::size_t row = 0; row < global.size(); ++row) {
for (std::size_t column = row;
column < global[row].size();
++column) {
double entry = 0.0;
for (std::size_t local_row = 0;
local_row < local.size();
++local_row) {
for (std::size_t local_column = 0;
local_column < local[local_row].size();
++local_column) {
entry +=
transformation[local_row][row] *
local[local_row][local_column] *
transformation[local_column][column];
}
}
global[row][column] = entry;
}
}
mirror_upper_triangle(global);
return global;
}
std::array<double, 12> transform_displacement(
const Matrix12& transformation,
const std::span<const double, 12> global_displacement) {
std::array<double, 12> local_displacement{};
for (std::size_t row = 0; row < transformation.size(); ++row) {
for (std::size_t column = 0;
column < transformation[row].size();
++column) {
local_displacement[row] +=
transformation[row][column] * global_displacement[column];
}
}
return local_displacement;
}
std::array<double, 6> evaluate_strain(
const StrainMatrix& strain_matrix,
const std::array<double, 12>& local_displacement) {
std::array<double, 6> strain{};
for (std::size_t component = 0; component < strain.size(); ++component) {
for (std::size_t dof = 0; dof < local_displacement.size(); ++dof) {
strain[component] +=
strain_matrix[component][dof] * local_displacement[dof];
}
}
return strain;
}
bool is_finite(const Matrix12& matrix) {
for (const auto& row : matrix) {
for (const double value : row) {
if (!std::isfinite(value)) {
return false;
}
}
}
return true;
}
} // namespace
BeamKernelResult compute_beam3d2(const Beam3D2Input& input) {
if (const auto validation = validate_properties(input);
validation.has_value()) {
return *validation;
}
BeamFrameResult frame_result = make_beam_frame(
input.coordinates[0],
input.coordinates[1],
input.section.orientation);
if (!frame_result.frame.has_value()) {
return {std::nullopt, std::move(frame_result.diagnostics)};
}
const Vec3 axis{
input.coordinates[1].x - input.coordinates[0].x,
input.coordinates[1].y - input.coordinates[0].y,
input.coordinates[1].z - input.coordinates[0].z,
};
const double length = std::hypot(axis.x, axis.y, axis.z);
const double jacobian = length / 2.0;
if (!std::isfinite(jacobian) || jacobian <= 0.0) {
return error_result(
"model.zero_length_element",
"Beam kernel requires a representable positive Jacobian.");
}
const std::array<double, 6> constitutive =
constitutive_values(input);
for (const double value : constitutive) {
if (!is_positive_finite(value)) {
return error_result(
"model.nonfinite_value",
"Beam constitutive stiffness is not finite and positive.");
}
}
Matrix12 local{};
integrate_components(
local,
gauss_rule_1d(2),
std::array<std::size_t, 4>{0, 3, 4, 5},
constitutive,
jacobian);
integrate_components(
local,
gauss_rule_1d(1),
std::array<std::size_t, 2>{1, 2},
constitutive,
jacobian);
mirror_upper_triangle(local);
const Matrix12 transformation =
beam_transformation(*frame_result.frame);
Matrix12 global = transform_stiffness(local, transformation);
if (!is_finite(local) || !is_finite(global)) {
return error_result(
"model.nonfinite_value",
"Beam stiffness contains a nonfinite value.");
}
return {
Beam3D2Contribution{
std::move(local),
std::move(global),
*frame_result.frame,
},
{},
};
}
std::vector<BeamSectionResult> recover_beam3d2(
const Beam3D2Input& input,
const std::span<const double, 12> element_displacement,
const std::span<const std::array<double, 2>> recovery_points) {
const BeamKernelResult kernel = compute_beam3d2(input);
if (!kernel.contribution.has_value()) {
const std::string message = kernel.diagnostics.empty()
? "Beam recovery requires a valid Beam3D2 input."
: kernel.diagnostics.front().message;
throw std::invalid_argument{message};
}
const Vec3 axis{
input.coordinates[1].x - input.coordinates[0].x,
input.coordinates[1].y - input.coordinates[0].y,
input.coordinates[1].z - input.coordinates[0].z,
};
const double jacobian = std::hypot(axis.x, axis.y, axis.z) / 2.0;
const Matrix12 transformation =
beam_transformation(kernel.contribution->frame);
const std::array<double, 12> local_displacement =
transform_displacement(transformation, element_displacement);
const std::array<double, 6> center_strain = evaluate_strain(
strain_matrix(0.0, jacobian),
local_displacement);
const std::array<double, 6> constitutive =
constitutive_values(input);
std::vector<BeamSectionResult> results;
results.reserve(input.node_ids.size());
for (std::size_t end = 0; end < input.node_ids.size(); ++end) {
const double xi = end == 0 ? -1.0 : 1.0;
std::array<double, 6> section_strain = evaluate_strain(
strain_matrix(xi, jacobian),
local_displacement);
section_strain[1] = center_strain[1];
section_strain[2] = center_strain[2];
std::array<double, 6> section_force{};
for (std::size_t component = 0;
component < section_force.size();
++component) {
section_force[component] =
constitutive[component] * section_strain[component];
}
std::vector<double> sigma_xx;
sigma_xx.reserve(recovery_points.size());
for (const auto& point : recovery_points) {
const double y = point[0];
const double z = point[1];
sigma_xx.push_back(
input.material.young *
(section_strain[0] + z * section_strain[4] -
y * section_strain[5]));
}
results.push_back({
xi,
input.node_ids[end],
section_strain,
section_force,
section_force[0] / input.section.area,
std::move(sigma_xx),
});
}
return results;
}
} // namespace fesa
+146
View File
@@ -0,0 +1,146 @@
#include <fesa/fem/beam_frame.hpp>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <utility>
namespace fesa {
namespace {
constexpr double kParallelTolerance =
64.0 * std::numeric_limits<double>::epsilon();
double length(const Vec3 value) {
return std::hypot(value.x, value.y, value.z);
}
double max_abs_component(const Vec3 value) {
return std::max(
std::abs(value.x),
std::max(std::abs(value.y), std::abs(value.z)));
}
double dot(const Vec3 first, const Vec3 second) {
return first.x * second.x + first.y * second.y + first.z * second.z;
}
Vec3 cross(const Vec3 first, const Vec3 second) {
return {
first.y * second.z - first.z * second.y,
first.z * second.x - first.x * second.z,
first.x * second.y - first.y * second.x,
};
}
Vec3 normalized(const Vec3 value) {
const double scale = max_abs_component(value);
const Vec3 scaled{
value.x / scale,
value.y / scale,
value.z / scale,
};
const double scaled_length = length(scaled);
return {
scaled.x / scaled_length,
scaled.y / scaled_length,
scaled.z / scaled_length,
};
}
BeamFrameResult error_result(std::string code, std::string message) {
BeamFrameResult result;
result.diagnostics.push_back({
DiagnosticStage::model,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
});
return result;
}
} // namespace
BeamFrameResult make_beam_frame(
const Vec3& first,
const Vec3& second,
const Vec3& orientation) {
if (!is_finite(first) || !is_finite(second)) {
return error_result(
"model.nonfinite_value",
"Beam frame requires finite node coordinates.");
}
const Vec3 axis{
second.x - first.x,
second.y - first.y,
second.z - first.z,
};
if (!is_finite(axis)) {
return error_result(
"model.nonfinite_value",
"Beam frame axis is not representable as a finite vector.");
}
if (max_abs_component(axis) == 0.0) {
return error_result(
"model.zero_length_element",
"Beam frame requires distinct node coordinates.");
}
if (!is_finite(orientation)) {
return error_result(
"model.nonfinite_value",
"Beam frame requires a finite orientation vector.");
}
if (max_abs_component(orientation) == 0.0) {
return error_result(
"model.invalid_orientation",
"Beam frame requires a nonzero orientation vector.");
}
const Vec3 ex = normalized(axis);
const Vec3 orientation_unit = normalized(orientation);
const double axial_projection = dot(orientation_unit, ex);
const Vec3 transverse{
orientation_unit.x - axial_projection * ex.x,
orientation_unit.y - axial_projection * ex.y,
orientation_unit.z - axial_projection * ex.z,
};
const double transverse_length = length(transverse);
if (transverse_length <= kParallelTolerance) {
return error_result(
"model.invalid_orientation",
"Beam orientation is parallel to the element axis.");
}
const Vec3 ey = normalized(transverse);
const Vec3 ez = cross(ex, ey);
return {
BeamFrame{ex, ey, ez},
{},
};
}
Matrix12 beam_transformation(const BeamFrame& frame) {
Matrix12 transformation{};
const std::array<Vec3, 3> basis{frame.ex, frame.ey, frame.ez};
constexpr std::array<std::size_t, 4> block_offsets{0, 3, 6, 9};
for (const std::size_t offset : block_offsets) {
for (std::size_t row = 0; row < basis.size(); ++row) {
transformation[offset + row][offset] = basis[row].x;
transformation[offset + row][offset + 1] = basis[row].y;
transformation[offset + row][offset + 2] = basis[row].z;
}
}
return transformation;
}
} // namespace fesa
+129
View File
@@ -0,0 +1,129 @@
#include <fesa/fem/dof_manager.hpp>
#include <algorithm>
#include <stdexcept>
#include <vector>
namespace fesa {
namespace {
constexpr std::size_t dofs_per_node = 6;
} // namespace
DofManager DofManager::build(const Domain& domain) {
DofManager manager;
std::vector<NodeId> node_ids;
node_ids.reserve(domain.nodes().size());
for (const Node& node : domain.nodes()) {
node_ids.push_back(node.id);
}
std::ranges::sort(node_ids);
for (std::size_t index = 0; index < node_ids.size(); ++index) {
manager.node_full_dof_bases_.emplace(
node_ids[index].value(), index * dofs_per_node);
}
const std::size_t full_dof_count = node_ids.size() * dofs_per_node;
manager.equations_.resize(full_dof_count);
manager.prescribed_values_.resize(full_dof_count);
std::vector<bool> constrained(full_dof_count, false);
for (const PrescribedDof& prescribed : domain.step().prescribed_dofs) {
const NodeDof dof =
static_cast<NodeDof>(prescribed.dof - std::uint8_t{1});
const std::size_t full =
manager.full_dof({prescribed.node, dof});
constrained[full] = true;
manager.prescribed_values_[full] = prescribed.value;
}
for (std::size_t full = 0; full < full_dof_count; ++full) {
if (!constrained[full]) {
manager.equations_[full] = manager.free_equation_count_;
++manager.free_equation_count_;
}
}
return manager;
}
std::size_t DofManager::full_dof_count() const noexcept {
return equations_.size();
}
std::size_t DofManager::free_equation_count() const noexcept {
return free_equation_count_;
}
std::optional<std::size_t> DofManager::equation(
const DofAddress address) const {
return equation(full_dof(address));
}
std::optional<std::size_t> DofManager::equation(
const std::size_t full_dof) const {
return equations_.at(full_dof);
}
std::optional<double> DofManager::prescribed_value(
const DofAddress address) const {
return prescribed_value(full_dof(address));
}
std::optional<double> DofManager::prescribed_value(
const std::size_t full_dof) const {
if (equations_.at(full_dof).has_value()) {
return std::nullopt;
}
return prescribed_values_.at(full_dof);
}
std::array<std::size_t, 12> DofManager::element_full_dofs(
const BeamElement& element) const {
std::array<std::size_t, 12> full_dofs{};
for (std::size_t node = 0; node < element.nodes.size(); ++node) {
const std::size_t base = node_full_dof_base(element.nodes[node]);
for (std::size_t component = 0; component < dofs_per_node;
++component) {
full_dofs[node * dofs_per_node + component] = base + component;
}
}
return full_dofs;
}
std::vector<double> DofManager::reconstruct_full(
const std::span<const double> reduced) const {
if (reduced.size() != free_equation_count_) {
throw std::invalid_argument{
"Reduced vector size must equal the free equation count."};
}
std::vector<double> full = prescribed_values_;
for (std::size_t index = 0; index < equations_.size(); ++index) {
if (equations_[index].has_value()) {
full[index] = reduced[*equations_[index]];
}
}
return full;
}
std::size_t DofManager::full_dof(const DofAddress address) const {
const auto component = static_cast<std::uint8_t>(address.dof);
if (component >= dofs_per_node) {
throw std::invalid_argument{"Node DOF must be one of ux through rz."};
}
return node_full_dof_base(address.node) + component;
}
std::size_t DofManager::node_full_dof_base(const NodeId node) const {
const auto found = node_full_dof_bases_.find(node.value());
if (found == node_full_dof_bases_.end()) {
throw std::out_of_range{"Node ID is not present in the DofManager."};
}
return found->second;
}
} // namespace fesa
+33
View File
@@ -0,0 +1,33 @@
#include <fesa/fem/gauss_rule.hpp>
#include <array>
#include <stdexcept>
namespace fesa {
namespace {
constexpr std::array<GaussPoint1D, 1> order_one{{
{0.0, 2.0},
}};
constexpr double inverse_sqrt_three = 0.57735026918962576451;
constexpr std::array<GaussPoint1D, 2> order_two{{
{-inverse_sqrt_three, 1.0},
{inverse_sqrt_three, 1.0},
}};
} // namespace
std::span<const GaussPoint1D> gauss_rule_1d(const int order) {
switch (order) {
case 1:
return order_one;
case 2:
return order_two;
default:
throw std::invalid_argument{"Gauss rule order must be 1 or 2."};
}
}
} // namespace fesa
+31
View File
@@ -0,0 +1,31 @@
#include <fesa/fem/line2_shape.hpp>
#include <cmath>
#include <stdexcept>
namespace fesa {
std::array<double, 2> line2_shape(const double xi) {
return {0.5 * (1.0 - xi), 0.5 * (1.0 + xi)};
}
std::array<double, 2> line2_shape_derivative(const double) {
return {-0.5, 0.5};
}
double line2_jacobian(const double length) {
if (!std::isfinite(length) || length <= 0.0) {
throw std::invalid_argument{
"Line2 Jacobian requires a positive finite length."};
}
const double jacobian = length / 2.0;
if (jacobian == 0.0) {
throw std::invalid_argument{
"Line2 Jacobian must be representable as a positive double."};
}
return jacobian;
}
} // namespace fesa
+149
View File
@@ -0,0 +1,149 @@
#include <fesa/io/abaqus/active_input.hpp>
#include <algorithm>
#include <set>
#include <string>
#include <string_view>
#include <utility>
namespace fesa {
namespace {
const std::string* parameter(
const DeckRecord& record,
const std::string_view name) {
const auto found = record.parameters.find(name);
return found == record.parameters.end() ? nullptr : &found->second;
}
bool is_flat_model_record(const DeckRecord& record) {
return record.keyword == "NODE" || record.keyword == "ELEMENT" ||
record.keyword == "NSET" || record.keyword == "ELSET" ||
record.keyword == "BEAM GENERAL SECTION" ||
record.keyword == "TRANSVERSE SHEAR STIFFNESS";
}
ActiveInputResult failure(
std::string code,
std::string message,
const SourceLocation& source) {
return {
std::nullopt,
{{
DiagnosticStage::semantic,
Severity::error,
std::move(code),
std::move(message),
source,
}},
};
}
} // namespace
ActiveInputResult select_active_input(const ParsedDeck& deck) {
const bool hierarchical =
!deck.parts.empty() || deck.assembly.has_value();
if (!hierarchical) {
return {
ActiveInputView{
true,
{},
{},
deck.global_records,
{},
},
{},
};
}
const auto flat_record =
std::ranges::find_if(deck.global_records, is_flat_model_record);
if (flat_record != deck.global_records.end()) {
return failure(
"abaqus.semantic.mixed_mesh_organization",
"Flat model records cannot be mixed with Part/Assembly input.",
flat_record->source);
}
std::set<std::string, std::less<>> part_names;
for (const ParsedPart& part : deck.parts) {
if (!part_names.insert(part.name).second) {
return failure(
"abaqus.semantic.duplicate_part",
"Part name '" + part.name + "' is defined more than once.",
part.source);
}
}
if (!deck.assembly.has_value()) {
const SourceLocation source =
deck.parts.empty() ? SourceLocation{} : deck.parts.front().source;
return failure(
"abaqus.semantic.assembly_count",
"Hierarchical Phase 1 input requires exactly one Assembly.",
source);
}
const ParsedAssembly& assembly = *deck.assembly;
if (assembly.instances.size() != 1U) {
const SourceLocation& source = assembly.instances.size() > 1U
? assembly.instances[1].source
: assembly.source;
return failure(
"abaqus.semantic.instance_count",
"Phase 1 requires exactly one Instance.",
source);
}
const ParsedInstance& instance = assembly.instances.front();
if (!instance.transform_data.empty()) {
const SourceLocation& source = instance.transform_sources.empty()
? instance.source
: instance.transform_sources.front();
return failure(
"abaqus.semantic.instance_transform",
"Instance translation and rotation data are unsupported.",
source);
}
const auto part =
std::ranges::find(deck.parts, instance.part_name, &ParsedPart::name);
if (part == deck.parts.end()) {
return failure(
"abaqus.semantic.missing_part",
"Instance '" + instance.name + "' references missing Part '" +
instance.part_name + "'.",
instance.source);
}
for (const DeckRecord& record : assembly.records) {
if (record.keyword != "NSET" && record.keyword != "ELSET") {
continue;
}
const std::string* record_instance = parameter(record, "INSTANCE");
if (record_instance == nullptr || *record_instance != instance.name) {
const std::string* name = parameter(
record, record.keyword == "NSET" ? "NSET" : "ELSET");
return failure(
"abaqus.semantic.wrong_instance",
"Assembly set '" +
(name == nullptr ? std::string{} : *name) +
"' must reference the active Instance.",
record.source);
}
}
return {
ActiveInputView{
false,
part->name,
instance.name,
part->records,
assembly.records,
},
{},
};
}
} // namespace fesa
+481 -7
View File
@@ -2,9 +2,16 @@
#include <algorithm>
#include <array>
#include <charconv>
#include <cstddef>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <initializer_list>
#include <iterator>
#include <optional>
#include <set>
#include <sstream>
#include <string>
#include <string_view>
#include <utility>
@@ -13,14 +20,19 @@
namespace fesa {
namespace {
enum class Scope { global, part, assembly, instance };
enum class Scope { global, part, assembly, instance, step };
struct KeywordLine final {
std::string keyword;
std::map<std::string, std::string, std::less<>> parameters;
std::set<std::string, std::less<>> flag_parameters;
SourceLocation source;
};
const std::string* parameter(
const KeywordLine& keyword,
std::string_view name);
std::string_view trim(const std::string_view value) {
constexpr std::string_view whitespace{" \t\f\v\r\n"};
const std::size_t first = value.find_first_not_of(whitespace);
@@ -41,6 +53,20 @@ std::string uppercase_ascii(std::string value) {
return value;
}
std::string input_fingerprint(const std::string_view bytes) {
std::uint64_t fingerprint = 14695981039346656037ULL;
for (const char byte : bytes) {
fingerprint ^=
static_cast<std::uint8_t>(static_cast<unsigned char>(byte));
fingerprint *= 1099511628211ULL;
}
std::ostringstream encoded;
encoded << "fnv1a64:" << std::hex << std::setfill('0')
<< std::setw(16) << fingerprint;
return encoded.str();
}
std::vector<std::string> split_fields(const std::string_view value) {
std::vector<std::string> fields;
std::size_t first = 0;
@@ -93,17 +119,202 @@ bool is_supported_record(const std::string_view keyword) {
std::string_view{"ELASTIC"},
std::string_view{"ELEMENT"},
std::string_view{"ELSET"},
std::string_view{"END STEP"},
std::string_view{"HEADING"},
std::string_view{"MATERIAL"},
std::string_view{"NODE"},
std::string_view{"NSET"},
std::string_view{"OUTPUT"},
std::string_view{"PREPRINT"},
std::string_view{"RESTART"},
std::string_view{"STATIC"},
std::string_view{"STEP"},
std::string_view{"TRANSVERSE SHEAR STIFFNESS"},
};
return std::ranges::find(supported, keyword) != supported.end();
}
bool is_known_keyword(const std::string_view keyword) {
constexpr std::array scope_keywords{
std::string_view{"ASSEMBLY"},
std::string_view{"END ASSEMBLY"},
std::string_view{"END INSTANCE"},
std::string_view{"END PART"},
std::string_view{"INSTANCE"},
std::string_view{"PART"},
};
return is_supported_record(keyword) ||
std::ranges::find(scope_keywords, keyword) != scope_keywords.end();
}
ParseDeckResult invalid_record_scope(const KeywordLine& keyword) {
std::string code = "abaqus.syntax.invalid_step_scope";
if (keyword.keyword == "NODE") {
code = "abaqus.syntax.invalid_node_scope";
} else if (keyword.keyword == "ELEMENT") {
code = "abaqus.syntax.invalid_element_scope";
} else if (keyword.keyword == "NSET" || keyword.keyword == "ELSET") {
code = "abaqus.syntax.invalid_set_scope";
} else if (keyword.keyword == "MATERIAL" ||
keyword.keyword == "ELASTIC") {
code = "abaqus.syntax.invalid_material_scope";
} else if (keyword.keyword == "BEAM GENERAL SECTION") {
code = "abaqus.syntax.invalid_section_scope";
} else if (keyword.keyword == "TRANSVERSE SHEAR STIFFNESS") {
code = "abaqus.syntax.invalid_transverse_shear_scope";
} else if (keyword.keyword == "BOUNDARY") {
code = "abaqus.syntax.invalid_boundary_scope";
} else if (keyword.keyword == "CLOAD") {
code = "abaqus.syntax.invalid_cload_scope";
} else if (keyword.keyword == "HEADING") {
code = "abaqus.syntax.invalid_heading_scope";
} else if (keyword.keyword == "PREPRINT") {
code = "abaqus.syntax.invalid_preprint_scope";
} else if (keyword.keyword == "RESTART") {
code = "abaqus.syntax.invalid_restart_scope";
} else if (keyword.keyword == "OUTPUT") {
code = "abaqus.syntax.invalid_output_scope";
} else if (keyword.keyword == "PART") {
code = "abaqus.syntax.invalid_part_scope";
} else if (keyword.keyword == "ASSEMBLY") {
code = "abaqus.syntax.invalid_assembly_scope";
} else if (keyword.keyword == "INSTANCE") {
code = "abaqus.syntax.invalid_instance_scope";
}
return syntax_failure(
std::move(code),
"*" + keyword.keyword + " is invalid in the current input scope.",
keyword.source);
}
bool is_allowed_parameter(
const std::string_view name,
const std::initializer_list<std::string_view> allowed) {
return std::ranges::find(allowed, name) != allowed.end();
}
std::optional<ParseDeckResult> validate_parameter_names(
const KeywordLine& keyword,
const std::initializer_list<std::string_view> allowed) {
for (const auto& [name, value] : keyword.parameters) {
static_cast<void>(value);
if (!is_allowed_parameter(name, allowed)) {
return syntax_failure(
"abaqus.syntax.unsupported_parameter",
"Unsupported parameter '" + name + "' on *" +
keyword.keyword + ".",
keyword.source);
}
}
return std::nullopt;
}
std::optional<ParseDeckResult> unsupported_parameter_value(
const KeywordLine& keyword,
std::string_view name);
std::optional<ParseDeckResult> validate_parameter_schema(
const KeywordLine& keyword,
const std::initializer_list<std::string_view> valued,
const std::initializer_list<std::string_view> flags = {}) {
for (const auto& [name, value] : keyword.parameters) {
if (is_allowed_parameter(name, valued)) {
if (keyword.flag_parameters.contains(name) || value.empty()) {
return syntax_failure(
"abaqus.syntax.invalid_parameter",
"Parameter '" + name + "' on *" + keyword.keyword +
" requires a nonempty value.",
keyword.source);
}
} else if (is_allowed_parameter(name, flags)) {
if (!keyword.flag_parameters.contains(name)) {
return unsupported_parameter_value(keyword, name);
}
} else {
return syntax_failure(
"abaqus.syntax.unsupported_parameter",
"Unsupported parameter '" + name + "' on *" +
keyword.keyword + ".",
keyword.source);
}
}
return std::nullopt;
}
std::optional<ParseDeckResult> unsupported_parameter_value(
const KeywordLine& keyword,
const std::string_view name) {
return syntax_failure(
"abaqus.syntax.unsupported_parameter",
"Unsupported value for parameter '" + std::string{name} +
"' on *" + keyword.keyword + ".",
keyword.source);
}
std::optional<ParseDeckResult> validate_noop_parameters(
const KeywordLine& keyword) {
if (keyword.keyword == "HEADING") {
return validate_parameter_schema(keyword, {});
}
if (keyword.keyword == "PREPRINT") {
if (auto error = validate_parameter_schema(
keyword, {"ECHO", "MODEL", "HISTORY", "CONTACT"})) {
return error;
}
for (const auto& [name, value] : keyword.parameters) {
const std::string normalized = uppercase_ascii(value);
if (normalized != "YES" && normalized != "NO") {
return unsupported_parameter_value(keyword, name);
}
}
return std::nullopt;
}
if (keyword.keyword == "RESTART") {
if (auto error = validate_parameter_schema(
keyword, {"FREQUENCY"}, {"WRITE"})) {
return error;
}
if (const std::string* write = parameter(keyword, "WRITE");
write != nullptr && !write->empty()) {
return unsupported_parameter_value(keyword, "WRITE");
}
if (const std::string* frequency = parameter(keyword, "FREQUENCY");
frequency != nullptr) {
std::int64_t value = 0;
const auto parsed = std::from_chars(
frequency->data(), frequency->data() + frequency->size(), value);
if (frequency->empty() || parsed.ec != std::errc{} ||
parsed.ptr != frequency->data() + frequency->size() ||
value < 0) {
return unsupported_parameter_value(keyword, "FREQUENCY");
}
}
return std::nullopt;
}
if (keyword.keyword == "OUTPUT") {
if (auto error = validate_parameter_schema(
keyword, {"VARIABLE"}, {"FIELD", "HISTORY"})) {
return error;
}
const bool field = keyword.parameters.contains("FIELD");
const bool history = keyword.parameters.contains("HISTORY");
if (field == history) {
return unsupported_parameter_value(keyword, "FIELD/HISTORY");
}
const std::string* flag = parameter(
keyword, field ? std::string_view{"FIELD"}
: std::string_view{"HISTORY"});
if (flag == nullptr || !flag->empty()) {
return unsupported_parameter_value(
keyword, field ? "FIELD" : "HISTORY");
}
if (const std::string* variable = parameter(keyword, "VARIABLE");
variable != nullptr && uppercase_ascii(*variable) != "PRESELECT") {
return unsupported_parameter_value(keyword, "VARIABLE");
}
return std::nullopt;
}
return std::nullopt;
}
std::optional<KeywordLine> parse_keyword_line(
const std::string_view line,
const SourceLocation& source,
@@ -120,6 +331,7 @@ std::optional<KeywordLine> parse_keyword_line(
KeywordLine parsed{
uppercase_ascii(fields[0]),
{},
{},
source,
};
for (std::size_t index = 1; index < fields.size(); ++index) {
@@ -152,6 +364,9 @@ std::optional<KeywordLine> parse_keyword_line(
source);
return std::nullopt;
}
if (equals == std::string_view::npos) {
parsed.flag_parameters.insert(key);
}
}
return parsed;
}
@@ -176,20 +391,34 @@ ParseDeckResult missing_parameter(
} // namespace
ParseDeckResult parse_deck(const std::filesystem::path& path) {
std::ifstream input{path, std::ios::binary};
if (!input) {
std::ifstream file{path, std::ios::binary};
if (!file) {
return failure(
DiagnosticStage::io,
"abaqus.io.open_failed",
"Unable to open Abaqus input file.",
SourceLocation{path, 0U, 0U});
}
const std::string source_bytes{
std::istreambuf_iterator<char>{file},
std::istreambuf_iterator<char>{},
};
if (file.bad()) {
return failure(
DiagnosticStage::io,
"abaqus.io.read_failed",
"Failed while reading Abaqus input file.",
SourceLocation{path, 0U, 0U});
}
std::istringstream input{source_bytes};
ParsedDeck deck;
Scope scope = Scope::global;
std::optional<ParsedPart> current_part;
std::optional<ParsedAssembly> current_assembly;
std::optional<ParsedInstance> current_instance;
std::optional<SourceLocation> current_step_source;
bool completed_step = false;
DeckRecord* current_record = nullptr;
std::string line;
@@ -214,8 +443,18 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
const std::vector<std::string> fields = split_fields(content);
if (scope == Scope::instance) {
current_instance->transform_data.push_back(fields);
current_instance->transform_sources.push_back(SourceLocation{
path,
line_number,
first_nonspace + 1U,
});
} else if (current_record != nullptr) {
current_record->data.push_back(fields);
current_record->data_sources.push_back(SourceLocation{
path,
line_number,
first_nonspace + 1U,
});
} else {
return syntax_failure(
"abaqus.syntax.data_without_keyword",
@@ -239,6 +478,69 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
KeywordLine& keyword = *parsed;
current_record = nullptr;
if (keyword.keyword == "STEP") {
if (scope != Scope::global) {
return syntax_failure(
"abaqus.syntax.invalid_step_scope",
"*STEP is only valid in global input scope.",
source);
}
if (auto error = validate_parameter_schema(
keyword, {"NAME", "NLGEOM"})) {
return std::move(*error);
}
if (const std::string* name = parameter(keyword, "NAME");
name != nullptr && name->empty()) {
return syntax_failure(
"abaqus.syntax.invalid_parameter",
"*STEP parameter NAME requires a nonempty value.",
source);
}
if (const std::string* nlgeom = parameter(keyword, "NLGEOM");
nlgeom != nullptr && nlgeom->empty()) {
return syntax_failure(
"abaqus.syntax.invalid_parameter",
"*STEP parameter NLGEOM requires a value.",
source);
}
deck.global_records.push_back({
std::move(keyword.keyword),
std::move(keyword.parameters),
{},
source,
});
current_step_source = source;
scope = Scope::step;
continue;
}
if (keyword.keyword == "END STEP") {
if (auto error = validate_parameter_names(keyword, {})) {
return std::move(*error);
}
if (scope != Scope::step || !current_step_source.has_value()) {
return syntax_failure(
"abaqus.syntax.unexpected_end_step",
"*END STEP does not match an open *STEP.",
source);
}
deck.global_records.push_back({
std::move(keyword.keyword),
std::move(keyword.parameters),
{},
source,
});
current_step_source.reset();
completed_step = true;
scope = Scope::global;
continue;
}
if (scope == Scope::global && completed_step &&
is_known_keyword(keyword.keyword)) {
return invalid_record_scope(keyword);
}
if (keyword.keyword == "PART") {
if (scope != Scope::global) {
return syntax_failure(
@@ -246,6 +548,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*PART is only valid in global input scope.",
source);
}
if (auto error = validate_parameter_schema(keyword, {"NAME"})) {
return std::move(*error);
}
const std::string* name = parameter(keyword, "NAME");
if (name == nullptr || name->empty()) {
return missing_parameter(keyword, "NAME");
@@ -262,6 +567,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*END PART does not match an open *PART.",
source);
}
if (auto error = validate_parameter_schema(keyword, {})) {
return std::move(*error);
}
deck.parts.push_back(std::move(*current_part));
current_part.reset();
scope = Scope::global;
@@ -275,6 +583,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*ASSEMBLY is only valid in global input scope.",
source);
}
if (auto error = validate_parameter_schema(keyword, {"NAME"})) {
return std::move(*error);
}
if (deck.assembly.has_value() ||
current_assembly.has_value()) {
return syntax_failure(
@@ -299,6 +610,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*END ASSEMBLY does not match an open *ASSEMBLY.",
source);
}
if (auto error = validate_parameter_schema(keyword, {})) {
return std::move(*error);
}
deck.assembly = std::move(*current_assembly);
current_assembly.reset();
scope = Scope::global;
@@ -313,6 +627,10 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*INSTANCE is only valid in an open *ASSEMBLY.",
source);
}
if (auto error =
validate_parameter_schema(keyword, {"NAME", "PART"})) {
return std::move(*error);
}
const std::string* name = parameter(keyword, "NAME");
if (name == nullptr || name->empty()) {
return missing_parameter(keyword, "NAME");
@@ -322,7 +640,7 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
return missing_parameter(keyword, "PART");
}
current_instance =
ParsedInstance{*name, *part_name, {}, source};
ParsedInstance{*name, *part_name, {}, {}, source};
scope = Scope::instance;
continue;
}
@@ -336,6 +654,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"*END INSTANCE does not match an open *INSTANCE.",
source);
}
if (auto error = validate_parameter_schema(keyword, {})) {
return std::move(*error);
}
current_assembly->instances.push_back(
std::move(*current_instance));
current_instance.reset();
@@ -343,12 +664,149 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
continue;
}
if (keyword.keyword == "HEADING" ||
keyword.keyword == "PREPRINT") {
if (scope != Scope::global) {
return syntax_failure(
keyword.keyword == "HEADING"
? "abaqus.syntax.invalid_heading_scope"
: "abaqus.syntax.invalid_preprint_scope",
"*" + keyword.keyword +
" is only valid in global input scope.",
source);
}
if (auto error = validate_noop_parameters(keyword)) {
return std::move(*error);
}
deck.global_records.push_back({
std::move(keyword.keyword),
std::move(keyword.parameters),
{},
source,
});
if (deck.global_records.back().keyword == "HEADING") {
current_record = &deck.global_records.back();
}
continue;
}
if (keyword.keyword == "RESTART" || keyword.keyword == "OUTPUT") {
if (scope != Scope::step) {
return syntax_failure(
keyword.keyword == "RESTART"
? "abaqus.syntax.invalid_restart_scope"
: "abaqus.syntax.invalid_output_scope",
"*" + keyword.keyword +
" is only valid inside *STEP.",
source);
}
if (auto error = validate_noop_parameters(keyword)) {
return std::move(*error);
}
deck.global_records.push_back({
std::move(keyword.keyword),
std::move(keyword.parameters),
{},
source,
});
continue;
}
if (keyword.keyword == "STATIC") {
if (scope != Scope::step) {
return syntax_failure(
"abaqus.syntax.invalid_step_scope",
"*STATIC is only valid inside *STEP.",
source);
}
if (auto error = validate_parameter_names(keyword, {})) {
return std::move(*error);
}
} else if (keyword.keyword == "CLOAD") {
if (scope != Scope::step) {
return syntax_failure(
"abaqus.syntax.invalid_cload_scope",
"*CLOAD is only valid inside *STEP.",
source);
}
if (auto error = validate_parameter_names(keyword, {})) {
return std::move(*error);
}
} else if (keyword.keyword == "BOUNDARY") {
if (scope != Scope::global && scope != Scope::step) {
return syntax_failure(
"abaqus.syntax.invalid_boundary_scope",
"*BOUNDARY is valid only in global or Step scope.",
source);
}
if (auto error = validate_parameter_names(keyword, {})) {
return std::move(*error);
}
}
if (!is_supported_record(keyword.keyword)) {
return syntax_failure(
"abaqus.unsupported_keyword",
"Unsupported Abaqus keyword *" + keyword.keyword + ".",
source);
}
if (scope == Scope::step && keyword.keyword != "STATIC" &&
keyword.keyword != "BOUNDARY" && keyword.keyword != "CLOAD") {
return invalid_record_scope(keyword);
}
if (scope != Scope::instance) {
const bool mesh_scope = scope == Scope::global || scope == Scope::part;
const bool set_scope = mesh_scope || scope == Scope::assembly;
if ((keyword.keyword == "NODE" || keyword.keyword == "ELEMENT" ||
keyword.keyword == "BEAM GENERAL SECTION" ||
keyword.keyword == "TRANSVERSE SHEAR STIFFNESS") &&
!mesh_scope) {
return invalid_record_scope(keyword);
}
if ((keyword.keyword == "NSET" || keyword.keyword == "ELSET") &&
!set_scope) {
return invalid_record_scope(keyword);
}
if ((keyword.keyword == "MATERIAL" ||
keyword.keyword == "ELASTIC") &&
scope != Scope::global) {
return invalid_record_scope(keyword);
}
}
if (scope != Scope::instance) {
std::optional<ParseDeckResult> parameter_error;
if (keyword.keyword == "NODE" || keyword.keyword == "ELASTIC" ||
keyword.keyword == "TRANSVERSE SHEAR STIFFNESS") {
parameter_error = validate_parameter_schema(keyword, {});
} else if (keyword.keyword == "ELEMENT") {
parameter_error = validate_parameter_schema(
keyword, {"TYPE", "ELSET"});
} else if (keyword.keyword == "NSET") {
parameter_error = scope == Scope::assembly
? validate_parameter_schema(
keyword,
{"NSET", "INSTANCE"},
{"GENERATE"})
: validate_parameter_schema(
keyword, {"NSET"}, {"GENERATE"});
} else if (keyword.keyword == "ELSET") {
parameter_error = scope == Scope::assembly
? validate_parameter_schema(
keyword,
{"ELSET", "INSTANCE"},
{"GENERATE"})
: validate_parameter_schema(
keyword, {"ELSET"}, {"GENERATE"});
} else if (keyword.keyword == "MATERIAL") {
parameter_error = validate_parameter_schema(keyword, {"NAME"});
} else if (keyword.keyword == "BEAM GENERAL SECTION") {
parameter_error = validate_parameter_schema(
keyword, {"SECTION", "ELSET", "MATERIAL"});
}
if (parameter_error.has_value()) {
return std::move(*parameter_error);
}
}
DeckRecord next_record{
std::move(keyword.keyword),
@@ -359,7 +817,9 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
switch (scope) {
case Scope::global:
deck.global_records.push_back(std::move(next_record));
if (deck.global_records.back().keyword != "MATERIAL") {
current_record = &deck.global_records.back();
}
break;
case Scope::part:
current_part->records.push_back(std::move(next_record));
@@ -375,6 +835,10 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"abaqus.syntax.instance_local_keyword",
"Keyword records inside *INSTANCE are unsupported.",
source);
case Scope::step:
deck.global_records.push_back(std::move(next_record));
current_record = &deck.global_records.back();
break;
}
}
@@ -404,8 +868,18 @@ ParseDeckResult parse_deck(const std::filesystem::path& path) {
"Abaqus *ASSEMBLY is not closed by *END ASSEMBLY.",
current_assembly->source);
}
if (current_step_source.has_value()) {
return syntax_failure(
"abaqus.syntax.unclosed_step",
"Abaqus *STEP is not closed by *END STEP.",
*current_step_source);
}
return {std::move(deck), {}};
return {
std::move(deck),
{},
input_fingerprint(source_bytes),
};
}
} // namespace fesa
File diff suppressed because it is too large Load Diff
+411
View File
@@ -0,0 +1,411 @@
#include <fesa/io/abaqus/set_resolver.hpp>
#include <algorithm>
#include <charconv>
#include <compare>
#include <cstddef>
#include <cstdint>
#include <map>
#include <set>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
namespace fesa {
namespace {
struct SetNamespace final {
ResolvedSetScope scope;
std::string scope_name;
ResolvedSetKind kind;
auto operator<=>(const SetNamespace&) const = default;
};
struct SetKey final {
SetNamespace name_space;
std::string set_name;
auto operator<=>(const SetKey&) const = default;
};
struct RawMember final {
std::string text;
SourceLocation source;
};
struct RawSet final {
std::vector<RawMember> members;
};
enum class VisitState { unvisited, visiting, resolved, failed };
const std::string* parameter(
const DeckRecord& record,
const std::string_view name) {
const auto found = record.parameters.find(name);
return found == record.parameters.end() ? nullptr : &found->second;
}
bool has_parameter(
const DeckRecord& record,
const std::string_view name) {
return record.parameters.contains(name);
}
SourceLocation data_source(
const DeckRecord& record,
const std::size_t row) {
return row < record.data_sources.size()
? record.data_sources[row]
: record.source;
}
bool parse_positive_label(
const std::string_view text,
std::int64_t& value) {
const auto parsed =
std::from_chars(text.data(), text.data() + text.size(), value);
return parsed.ec == std::errc{} &&
parsed.ptr == text.data() + text.size() && value > 0;
}
class SetResolver final {
public:
explicit SetResolver(const ParsedDeck& deck) : deck_{deck} {}
[[nodiscard]] SetResolutionResult resolve() {
collect_scope(
deck_.global_records,
{ResolvedSetScope::global, "global", ResolvedSetKind::node});
for (const ParsedPart& part : deck_.parts) {
collect_scope(
part.records,
{ResolvedSetScope::part,
part.name,
ResolvedSetKind::node});
}
collect_assembly();
std::vector<ResolvedSet> sets;
sets.reserve(raw_sets_.size());
for (const auto& [key, raw_set] : raw_sets_) {
static_cast<void>(raw_set);
if (!resolve_set(key)) {
continue;
}
sets.push_back({
key.name_space.scope_name,
key.set_name,
resolved_sets_.at(key),
key.name_space.scope,
key.name_space.kind,
});
}
if (!diagnostics_.empty()) {
sets.clear();
}
return {std::move(sets), std::move(diagnostics_)};
}
private:
void add_error(
std::string code,
std::string message,
const SourceLocation& source) {
diagnostics_.push_back({
DiagnosticStage::semantic,
Severity::error,
std::move(code),
std::move(message),
source,
});
}
void collect_scope(
const std::vector<DeckRecord>& records,
SetNamespace name_space) {
collect_entities(records, name_space);
collect_set_records(records, std::move(name_space), nullptr);
}
void collect_entities(
const std::vector<DeckRecord>& records,
const SetNamespace& base_namespace) {
for (const DeckRecord& record : records) {
ResolvedSetKind kind;
if (record.keyword == "NODE") {
kind = ResolvedSetKind::node;
} else if (record.keyword == "ELEMENT") {
kind = ResolvedSetKind::element;
} else {
continue;
}
SetNamespace entity_namespace = base_namespace;
entity_namespace.kind = kind;
std::set<std::int64_t>& labels = entities_[entity_namespace];
for (std::size_t row = 0; row < record.data.size(); ++row) {
if (record.data[row].empty()) {
continue;
}
std::int64_t label = 0;
if (!parse_positive_label(record.data[row][0], label)) {
continue;
}
labels.insert(label);
if (kind != ResolvedSetKind::element) {
continue;
}
const std::string* set_name = parameter(record, "ELSET");
if (set_name == nullptr || set_name->empty()) {
continue;
}
raw_sets_[{entity_namespace, *set_name}].members.push_back({
std::to_string(label),
data_source(record, row),
});
}
}
}
void collect_set_records(
const std::vector<DeckRecord>& records,
const SetNamespace& base_namespace,
const std::string* active_instance) {
for (const DeckRecord& record : records) {
SetNamespace set_namespace = base_namespace;
std::string_view name_parameter;
if (record.keyword == "NSET") {
set_namespace.kind = ResolvedSetKind::node;
name_parameter = "NSET";
} else if (record.keyword == "ELSET") {
set_namespace.kind = ResolvedSetKind::element;
name_parameter = "ELSET";
} else {
continue;
}
const std::string* set_name = parameter(record, name_parameter);
if (set_name == nullptr || set_name->empty()) {
add_error(
"abaqus.semantic.missing_parameter",
"*" + record.keyword + " requires parameter " +
std::string{name_parameter} + ".",
record.source);
continue;
}
if (set_namespace.scope == ResolvedSetScope::assembly) {
const std::string* instance = parameter(record, "INSTANCE");
if (active_instance == nullptr || instance == nullptr ||
*instance != *active_instance) {
add_error(
"abaqus.semantic.wrong_instance",
"Assembly set '" + *set_name +
"' must reference the active Instance.",
record.source);
continue;
}
}
RawSet& raw_set = raw_sets_[{set_namespace, *set_name}];
if (has_parameter(record, "GENERATE")) {
collect_generate(record, raw_set);
} else {
collect_explicit(record, raw_set);
}
}
}
void collect_explicit(
const DeckRecord& record,
RawSet& raw_set) {
for (std::size_t row = 0; row < record.data.size(); ++row) {
for (const std::string& field : record.data[row]) {
if (field.empty()) {
continue;
}
raw_set.members.push_back({field, data_source(record, row)});
}
}
}
void collect_generate(
const DeckRecord& record,
RawSet& raw_set) {
const SourceLocation source =
record.data.empty() ? record.source : data_source(record, 0U);
if (record.data.size() != 1U || record.data[0].size() != 3U ||
std::ranges::any_of(
record.data[0],
[](const std::string& field) { return field.empty(); })) {
add_error(
"abaqus.semantic.invalid_generate",
"*" + record.keyword +
", GENERATE requires exactly start, end, increment.",
source);
return;
}
std::int64_t start = 0;
std::int64_t end = 0;
std::int64_t increment = 0;
if (!parse_positive_label(record.data[0][0], start) ||
!parse_positive_label(record.data[0][1], end) ||
!parse_positive_label(record.data[0][2], increment) ||
start > end || (end - start) % increment != 0) {
add_error(
"abaqus.semantic.invalid_generate",
"Invalid *" + record.keyword + " generate range.",
source);
return;
}
for (std::int64_t label = start;; label += increment) {
raw_set.members.push_back({std::to_string(label), source});
if (label == end) {
break;
}
}
}
void collect_assembly() {
if (!deck_.assembly.has_value()) {
return;
}
const ParsedAssembly& assembly = *deck_.assembly;
const ParsedInstance* active_instance =
assembly.instances.size() == 1U
? &assembly.instances.front()
: nullptr;
SetNamespace assembly_namespace{
ResolvedSetScope::assembly,
assembly.name,
ResolvedSetKind::node,
};
if (active_instance != nullptr) {
const auto part = std::ranges::find(
deck_.parts, active_instance->part_name, &ParsedPart::name);
if (part != deck_.parts.end()) {
for (const ResolvedSetKind kind : {
ResolvedSetKind::node,
ResolvedSetKind::element}) {
const SetNamespace part_namespace{
ResolvedSetScope::part,
part->name,
kind,
};
SetNamespace lifted_namespace = assembly_namespace;
lifted_namespace.kind = kind;
const auto labels = entities_.find(part_namespace);
if (labels != entities_.end()) {
entities_[lifted_namespace] = labels->second;
}
}
}
}
const std::string* instance_name =
active_instance == nullptr ? nullptr : &active_instance->name;
collect_set_records(
assembly.records, assembly_namespace, instance_name);
}
bool resolve_set(const SetKey& key) {
VisitState& state = states_[key];
if (state == VisitState::resolved) {
return true;
}
if (state == VisitState::failed) {
return false;
}
state = VisitState::visiting;
bool succeeded = true;
std::vector<std::int64_t> labels;
for (const RawMember& member : raw_sets_.at(key).members) {
std::int64_t label = 0;
if (parse_positive_label(member.text, label)) {
const auto entity_namespace = entities_.find(key.name_space);
if (entity_namespace == entities_.end() ||
!entity_namespace->second.contains(label)) {
add_error(
"abaqus.semantic.missing_set_member",
"Set '" + key.set_name +
"' references missing entity label " +
std::to_string(label) + ".",
member.source);
succeeded = false;
continue;
}
labels.push_back(label);
continue;
}
const SetKey nested_key{key.name_space, member.text};
const auto nested = raw_sets_.find(nested_key);
if (nested == raw_sets_.end()) {
add_error(
"abaqus.semantic.missing_set_member",
"Set '" + key.set_name + "' references missing set '" +
member.text + "'.",
member.source);
succeeded = false;
continue;
}
if (states_[nested_key] == VisitState::visiting) {
add_error(
"abaqus.semantic.set_cycle",
"Set '" + key.set_name +
"' closes a nested set reference cycle through '" +
member.text + "'.",
member.source);
succeeded = false;
continue;
}
if (!resolve_set(nested_key)) {
succeeded = false;
continue;
}
const std::vector<std::int64_t>& nested_labels =
resolved_sets_.at(nested_key);
labels.insert(
labels.end(), nested_labels.begin(), nested_labels.end());
}
if (!succeeded) {
state = VisitState::failed;
return false;
}
std::ranges::sort(labels);
labels.erase(std::ranges::unique(labels).begin(), labels.end());
resolved_sets_[key] = std::move(labels);
state = VisitState::resolved;
return true;
}
const ParsedDeck& deck_;
std::map<SetNamespace, std::set<std::int64_t>> entities_;
std::map<SetKey, RawSet> raw_sets_;
std::map<SetKey, VisitState> states_;
std::map<SetKey, std::vector<std::int64_t>> resolved_sets_;
std::vector<Diagnostic> diagnostics_;
};
} // namespace
SetResolutionResult resolve_sets(const ParsedDeck& deck) {
return SetResolver{deck}.resolve();
}
} // namespace fesa
File diff suppressed because it is too large Load Diff
+18
View File
@@ -55,6 +55,24 @@ const Node& Domain::node(const EntityOrigin& origin) const {
return nodes_[found->second];
}
const IsotropicElastic& Domain::material(const MaterialId id) const {
const auto found = material_indices_.find(id.value());
if (found == material_indices_.end()) {
throw std::out_of_range{
"Material ID is not present in the Domain."};
}
return materials_[found->second];
}
const BeamSection& Domain::section(const SectionId id) const {
const auto found = section_indices_.find(id.value());
if (found == section_indices_.end()) {
throw std::out_of_range{
"Section ID is not present in the Domain."};
}
return sections_[found->second];
}
Domain::OriginKey Domain::origin_key(const EntityOrigin& origin) {
return {origin.instance_name, origin.local_label};
}
+201
View File
@@ -0,0 +1,201 @@
#include <fesa/results/result_database.hpp>
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <set>
#include <string>
#include <utility>
namespace fesa {
namespace {
void add_error(
std::vector<Diagnostic>& diagnostics,
std::string code,
std::string message) {
diagnostics.push_back({
DiagnosticStage::results,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
});
}
bool is_finite(const std::array<double, 6>& field) {
return std::ranges::all_of(
field,
[](const double component) {
return std::isfinite(component);
});
}
void validate_nodal_frame(
const NodalFrame& nodal,
std::vector<Diagnostic>& diagnostics) {
if (
nodal.displacement.size() != nodal.node_ids.size() ||
nodal.reaction.size() != nodal.node_ids.size()) {
add_error(
diagnostics,
"results.nodal_size_mismatch",
"Nodal IDs, displacement, and reaction fields must have "
"matching sizes.");
}
if (nodal.origins.size() != nodal.node_ids.size()) {
add_error(
diagnostics,
"results.nodal_origin_size_mismatch",
"Nodal IDs and origins must have matching sizes.");
}
std::set<std::int64_t> node_ids;
for (const NodeId node_id : nodal.node_ids) {
if (!node_ids.insert(node_id.value()).second) {
add_error(
diagnostics,
"results.duplicate_node_id",
"Nodal frame contains duplicate node ID " +
std::to_string(node_id.value()) + ".");
}
}
for (const auto& displacement : nodal.displacement) {
if (!is_finite(displacement)) {
add_error(
diagnostics,
"results.nonfinite_value",
"Nodal displacement contains a nonfinite component.");
}
}
for (const auto& reaction : nodal.reaction) {
if (!is_finite(reaction)) {
add_error(
diagnostics,
"results.nonfinite_value",
"Nodal reaction contains a nonfinite component.");
}
}
}
void validate_beam_section_result(
const BeamSectionResult& result,
std::vector<Diagnostic>& diagnostics) {
if (
!std::isfinite(result.xi) ||
!is_finite(result.section_strain) ||
!is_finite(result.section_force) ||
!std::isfinite(result.centroid_sigma_xx) ||
!std::ranges::all_of(
result.sigma_xx,
[](const double value) { return std::isfinite(value); })) {
add_error(
diagnostics,
"results.nonfinite_value",
"Beam section result contains a nonfinite value.");
}
}
void validate_element_frame(
const ElementFrame& element,
const NodalFrame& nodal,
std::vector<Diagnostic>& diagnostics) {
std::set<std::int64_t> nodal_ids;
for (const NodeId node_id : nodal.node_ids) {
nodal_ids.insert(node_id.value());
}
std::set<std::int64_t> element_ids;
for (const BeamElementFrame& beam : element.beams) {
if (!element_ids.insert(beam.element.value()).second) {
add_error(
diagnostics,
"results.duplicate_element_id",
"Element frame contains duplicate element ID " +
std::to_string(beam.element.value()) + ".");
}
if (
!is_finite(beam.local_frame.ex) ||
!is_finite(beam.local_frame.ey) ||
!is_finite(beam.local_frame.ez)) {
add_error(
diagnostics,
"results.nonfinite_value",
"Beam local frame contains a nonfinite component.");
}
const BeamSectionResult& first = beam.end_results[0];
const BeamSectionResult& second = beam.end_results[1];
if (first.end_node == second.end_node) {
add_error(
diagnostics,
"results.duplicate_beam_end_node",
"Beam element result contains the same node at both ends.");
}
if (
first.xi != -1.0 || second.xi != 1.0 ||
!nodal_ids.contains(first.end_node.value()) ||
!nodal_ids.contains(second.end_node.value())) {
add_error(
diagnostics,
"results.invalid_beam_connectivity",
"Beam end results must follow (-1, +1) connectivity and "
"reference nodes in the nodal frame.");
}
if (first.sigma_xx.size() != second.sigma_xx.size()) {
add_error(
diagnostics,
"results.recovery_point_count_mismatch",
"Beam end results must have matching recovery-point "
"counts.");
}
validate_beam_section_result(first, diagnostics);
validate_beam_section_result(second, diagnostics);
}
}
} // namespace
Status validate_result_database(const ResultDatabase& database) {
std::vector<Diagnostic> diagnostics;
std::set<std::string> step_names;
for (const ResultStep& step : database.steps) {
if (!step_names.insert(step.name).second) {
add_error(
diagnostics,
"results.duplicate_step_name",
"Result database contains duplicate step name '" +
step.name + "'.");
}
std::set<double> frame_times;
for (const ResultFrame& frame : step.frames) {
if (!std::isfinite(frame.step_time)) {
add_error(
diagnostics,
"results.nonfinite_value",
"Result frame has a nonfinite step time.");
} else if (!frame_times.insert(frame.step_time).second) {
add_error(
diagnostics,
"results.duplicate_frame_time",
"Result step '" + step.name +
"' contains duplicate frame time " +
std::to_string(frame.step_time) + ".");
}
validate_nodal_frame(frame.nodal, diagnostics);
validate_element_frame(
frame.element, frame.nodal, diagnostics);
}
}
return {diagnostics.empty(), std::move(diagnostics)};
}
} // namespace fesa
@@ -0,0 +1,330 @@
#include <fesa/solvers/linear/pardiso_linear_solver.hpp>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <optional>
#include <string>
#include <string_view>
#include <type_traits>
#include <utility>
#include <vector>
#include <mkl.h>
namespace fesa {
namespace {
static_assert(
std::is_same_v<MKL_INT, std::int32_t>,
"FESA requires the oneMKL LP64 interface.");
Diagnostic solver_error(std::string code, std::string message) {
return {
DiagnosticStage::solver,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
};
}
std::optional<std::string> validate_matrix(
const SymmetricCsr& matrix) {
if (matrix.order == 0) {
return "PARDISO requires a nonempty reduced system.";
}
if (matrix.order >
static_cast<std::size_t>(
std::numeric_limits<MKL_INT>::max())) {
return "Matrix order exceeds the oneMKL LP64 index range.";
}
if (matrix.row_offsets.size() != matrix.order + 1 ||
matrix.row_offsets.front() != 0) {
return "CSR row offsets must contain order + 1 entries "
"starting at zero.";
}
if (matrix.column_indices.size() != matrix.values.size()) {
return "CSR column and value counts must match.";
}
MKL_INT previous_offset = 0;
for (const MKL_INT offset : matrix.row_offsets) {
if (offset < previous_offset || offset < 0 ||
static_cast<std::size_t>(offset) >
matrix.column_indices.size()) {
return "CSR row offsets must be nondecreasing and in range.";
}
previous_offset = offset;
}
if (static_cast<std::size_t>(matrix.row_offsets.back()) !=
matrix.column_indices.size()) {
return "The final CSR row offset must equal the entry count.";
}
for (std::size_t row = 0; row < matrix.order; ++row) {
MKL_INT previous_column = -1;
bool has_diagonal = false;
const std::size_t begin =
static_cast<std::size_t>(matrix.row_offsets[row]);
const std::size_t end =
static_cast<std::size_t>(matrix.row_offsets[row + 1]);
for (std::size_t entry = begin; entry < end; ++entry) {
const MKL_INT column = matrix.column_indices[entry];
if (column < static_cast<MKL_INT>(row) ||
column >= static_cast<MKL_INT>(matrix.order) ||
column <= previous_column) {
return "CSR rows must contain sorted unique "
"upper-triangle columns.";
}
if (!std::isfinite(matrix.values[entry])) {
return "CSR values must be finite.";
}
has_diagonal =
has_diagonal || column == static_cast<MKL_INT>(row);
previous_column = column;
}
if (!has_diagonal) {
return "Every CSR row must contain its diagonal entry.";
}
}
return std::nullopt;
}
class PardisoSession final {
public:
PardisoSession() {
pardisoinit(handles_.data(), &matrix_type_, parameters_.data());
parameters_[26] = 1;
parameters_[34] = 1;
}
PardisoSession(const PardisoSession&) = delete;
PardisoSession& operator=(const PardisoSession&) = delete;
~PardisoSession() noexcept {
static_cast<void>(release());
}
MKL_INT release() noexcept {
if (!active_) {
return 0;
}
active_ = false;
constexpr MKL_INT release_all = -1;
MKL_INT error = 0;
pardiso(
handles_.data(),
&max_factorizations_,
&matrix_number_,
&matrix_type_,
&release_all,
&order_,
matrix_->values.data(),
matrix_->row_offsets.data(),
matrix_->column_indices.data(),
permutation_.data(),
&right_hand_side_count_,
parameters_.data(),
&message_level_,
right_hand_side_->data(),
solution_->data(),
&error);
return error;
}
MKL_INT execute(
const MKL_INT phase,
const SymmetricCsr& matrix,
std::vector<double>& right_hand_side,
std::vector<double>& solution) {
order_ = static_cast<MKL_INT>(matrix.order);
matrix_ = &matrix;
right_hand_side_ = &right_hand_side;
solution_ = &solution;
permutation_.resize(matrix.order);
active_ = true;
MKL_INT error = 0;
pardiso(
handles_.data(),
&max_factorizations_,
&matrix_number_,
&matrix_type_,
&phase,
&order_,
matrix.values.data(),
matrix.row_offsets.data(),
matrix.column_indices.data(),
permutation_.data(),
&right_hand_side_count_,
parameters_.data(),
&message_level_,
right_hand_side.data(),
solution.data(),
&error);
return error;
}
private:
std::array<void*, 64> handles_{};
std::array<MKL_INT, 64> parameters_{};
std::vector<MKL_INT> permutation_;
const SymmetricCsr* matrix_{};
std::vector<double>* right_hand_side_{};
std::vector<double>* solution_{};
MKL_INT order_{};
MKL_INT max_factorizations_{1};
MKL_INT matrix_number_{1};
MKL_INT matrix_type_{2};
MKL_INT right_hand_side_count_{1};
MKL_INT message_level_{};
bool active_{};
};
double relative_residual(
const SymmetricCsr& matrix,
const std::span<const double> rhs,
const std::span<const double> solution) {
std::vector<double> residual(rhs.begin(), rhs.end());
for (double& value : residual) {
value = -value;
}
for (std::size_t row = 0; row < matrix.order; ++row) {
const std::size_t begin =
static_cast<std::size_t>(matrix.row_offsets[row]);
const std::size_t end =
static_cast<std::size_t>(matrix.row_offsets[row + 1]);
for (std::size_t entry = begin; entry < end; ++entry) {
const std::size_t column =
static_cast<std::size_t>(matrix.column_indices[entry]);
const double value = matrix.values[entry];
residual[row] += value * solution[column];
if (column != row) {
residual[column] += value * solution[row];
}
}
}
double residual_squared = 0.0;
double rhs_squared = 0.0;
for (std::size_t index = 0; index < rhs.size(); ++index) {
residual_squared += residual[index] * residual[index];
rhs_squared += rhs[index] * rhs[index];
}
const double residual_norm = std::sqrt(residual_squared);
const double rhs_norm = std::sqrt(rhs_squared);
return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm;
}
std::string pardiso_failure(
const std::string_view phase,
const MKL_INT error) {
return "PARDISO " + std::string{phase} +
" failed with error " + std::to_string(error) + ".";
}
} // namespace
PardisoLinearSolver::PardisoLinearSolver() = default;
PardisoLinearSolver::~PardisoLinearSolver() = default;
LinearSolveResult PardisoLinearSolver::solve(
const SymmetricCsr& matrix,
const std::span<const double> rhs) {
LinearSolveResult result;
if (const auto error = validate_matrix(matrix);
error.has_value()) {
result.diagnostics.push_back(
solver_error("solver.invalid_csr", *error));
return result;
}
if (rhs.size() != matrix.order) {
result.diagnostics.push_back(solver_error(
"solver.dimension_mismatch",
"Right-hand side size must equal the matrix order."));
return result;
}
if (!std::ranges::all_of(rhs, [](const double value) {
return std::isfinite(value);
})) {
result.diagnostics.push_back(solver_error(
"solver.invalid_rhs",
"Right-hand side values must be finite."));
return result;
}
std::vector<double> right_hand_side(rhs.begin(), rhs.end());
std::vector<double> solution(matrix.order, 0.0);
{
PardisoSession session;
constexpr MKL_INT analyze = 11;
MKL_INT error = session.execute(
analyze, matrix, right_hand_side, solution);
if (error != 0) {
result.diagnostics.push_back(solver_error(
"solver.analysis_failed",
pardiso_failure("analysis", error)));
}
if (error == 0) {
constexpr MKL_INT factorize = 22;
error = session.execute(
factorize, matrix, right_hand_side, solution);
if (error != 0) {
result.diagnostics.push_back(solver_error(
"solver.factorization_failed",
pardiso_failure("factorization", error)));
}
}
if (error == 0) {
constexpr MKL_INT solve_system = 33;
error = session.execute(
solve_system, matrix, right_hand_side, solution);
if (error != 0) {
result.diagnostics.push_back(solver_error(
"solver.solve_failed",
pardiso_failure("solve", error)));
}
}
if (const MKL_INT release_error = session.release();
release_error != 0) {
result.diagnostics.push_back(solver_error(
"solver.release_failed",
pardiso_failure("release", release_error)));
}
}
if (!result.diagnostics.empty()) {
return result;
}
if (!std::ranges::all_of(solution, [](const double value) {
return std::isfinite(value);
})) {
result.diagnostics.push_back(solver_error(
"solver.nonfinite_solution",
"PARDISO produced a nonfinite solution."));
return result;
}
result.relative_residual =
relative_residual(matrix, rhs, solution);
if (!std::isfinite(result.relative_residual)) {
result.diagnostics.push_back(solver_error(
"solver.nonfinite_residual",
"The independently computed relative residual is nonfinite."));
return result;
}
result.solution = std::move(solution);
return result;
}
} // namespace fesa
+591
View File
@@ -0,0 +1,591 @@
#include <fesa/validation/comparison.hpp>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <limits>
#include <map>
#include <optional>
#include <set>
#include <sstream>
#include <string>
#include <string_view>
#include <tuple>
#include <utility>
#include <vector>
namespace fesa {
namespace {
using PositionKey = std::tuple<
ReferenceQuantity,
std::string,
std::int64_t,
std::optional<std::int64_t>>;
void add_failure(
std::vector<Diagnostic>& failures,
std::string code,
std::string message) {
failures.push_back({
DiagnosticStage::validation,
Severity::error,
std::move(code),
std::move(message),
std::nullopt,
});
}
std::string_view quantity_name(const ReferenceQuantity quantity) {
switch (quantity) {
case ReferenceQuantity::displacement:
return "displacement";
case ReferenceQuantity::reaction:
return "reaction";
case ReferenceQuantity::internal_force:
return "internal_force";
case ReferenceQuantity::centroid_stress:
return "centroid_stress";
}
return "unknown";
}
std::size_t expected_component_count(
const ReferenceQuantity quantity) {
switch (quantity) {
case ReferenceQuantity::displacement:
case ReferenceQuantity::reaction:
case ReferenceQuantity::internal_force:
return 6;
case ReferenceQuantity::centroid_stress:
return 1;
}
return 0;
}
std::string component_name(
const ReferenceQuantity quantity,
const std::size_t index) {
switch (quantity) {
case ReferenceQuantity::displacement:
return std::string{NodalFrame::displacement_components[index]};
case ReferenceQuantity::reaction:
return std::string{NodalFrame::reaction_components[index]};
case ReferenceQuantity::internal_force:
return std::string{
BeamElementFrame::section_force_components[index]};
case ReferenceQuantity::centroid_stress:
return std::string{BeamElementFrame::axial_stress_component};
}
return "unknown";
}
std::string number_text(const double value) {
std::ostringstream stream;
stream << std::setprecision(std::numeric_limits<double>::max_digits10)
<< value;
return stream.str();
}
std::string tolerance_text(const Tolerance tolerance) {
return " relative_tolerance=" + number_text(tolerance.relative) +
" absolute_scale=" + number_text(tolerance.absolute_scale);
}
std::string position_text(
const ReferenceQuantity quantity,
const ResultPosition& position) {
std::string text = "quantity=" + std::string{quantity_name(quantity)} +
" instance=" + position.instance_name +
" entity=" + std::to_string(position.entity_label);
if (position.end_node_label.has_value()) {
text += " end_node=" +
std::to_string(*position.end_node_label);
} else {
text += " end_node=n/a";
}
return text;
}
std::string scalar_failure_text(
const ComparisonSample& sample,
const std::size_t component,
const double normalized_error) {
return position_text(sample.quantity, sample.position) +
" component=" + component_name(sample.quantity, component) +
" reference=" + number_text(sample.reference[component]) +
" actual=" + number_text(sample.actual[component]) +
" normalized_error=" + number_text(normalized_error) +
tolerance_text(sample.tolerance);
}
std::string unevaluable_failure_text(
const ReferenceQuantity quantity,
const ResultPosition& position,
const Tolerance tolerance,
const std::string_view reason) {
return position_text(quantity, position) +
" component=n/a reference=n/a actual=n/a "
"normalized_error=inf" + tolerance_text(tolerance) +
" reason=" + std::string{reason};
}
bool origin_matches(
const EntityOrigin& origin,
const std::string& instance_name,
const std::int64_t local_label) {
return origin.instance_name == instance_name &&
origin.local_label == local_label;
}
ComparisonSampleMatch matching_failure(
const ReferenceQuantity quantity,
const ResultPosition& position,
const Tolerance tolerance,
std::string code,
const std::string_view reason) {
std::vector<Diagnostic> failures;
add_failure(
failures,
std::move(code),
unevaluable_failure_text(
quantity, position, tolerance, reason));
return {std::nullopt, std::move(failures)};
}
std::vector<double> as_vector(const std::array<double, 6>& values) {
return {values.begin(), values.end()};
}
} // namespace
ComparisonSampleMatch make_comparison_sample(
const ResultFrame& frame,
const ReferenceQuantity quantity,
const ResultPosition& position,
const std::span<const double> reference,
const Tolerance tolerance) {
std::vector<double> actual;
if (
quantity == ReferenceQuantity::displacement ||
quantity == ReferenceQuantity::reaction) {
if (position.end_node_label.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.invalid_result_position",
"nodal_position_has_end_node");
}
std::optional<std::size_t> matched_index;
for (
std::size_t index = 0;
index < frame.nodal.origins.size();
++index) {
if (origin_matches(
frame.nodal.origins[index],
position.instance_name,
position.entity_label)) {
if (matched_index.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"ambiguous_result_origin");
}
matched_index = index;
}
}
if (!matched_index.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"unknown_result_origin");
}
const auto& field =
quantity == ReferenceQuantity::displacement
? frame.nodal.displacement
: frame.nodal.reaction;
if (*matched_index >= field.size()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.component_count_mismatch",
"missing_actual_components");
}
actual = as_vector(field[*matched_index]);
} else {
if (!position.end_node_label.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.invalid_element_node_pair",
"missing_end_node");
}
const BeamElementFrame* matched_beam = nullptr;
for (const BeamElementFrame& beam : frame.element.beams) {
if (origin_matches(
beam.origin,
position.instance_name,
position.entity_label)) {
if (matched_beam != nullptr) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"ambiguous_result_origin");
}
matched_beam = &beam;
}
}
if (matched_beam == nullptr) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"unknown_result_origin");
}
std::optional<NodeId> end_node;
for (
std::size_t index = 0;
index < frame.nodal.origins.size() &&
index < frame.nodal.node_ids.size();
++index) {
if (origin_matches(
frame.nodal.origins[index],
position.instance_name,
*position.end_node_label)) {
if (end_node.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"ambiguous_end_node_origin");
}
end_node = frame.nodal.node_ids[index];
}
}
if (!end_node.has_value()) {
return matching_failure(
quantity,
position,
tolerance,
"validation.unknown_result_origin",
"unknown_end_node_origin");
}
const BeamSectionResult* matched_end = nullptr;
for (const BeamSectionResult& end : matched_beam->end_results) {
if (end.end_node == *end_node) {
matched_end = &end;
break;
}
}
if (matched_end == nullptr) {
return matching_failure(
quantity,
position,
tolerance,
"validation.invalid_element_node_pair",
"node_is_not_element_end");
}
if (quantity == ReferenceQuantity::internal_force) {
actual = as_vector(matched_end->section_force);
} else {
actual = {matched_end->centroid_sigma_xx};
}
}
ComparisonSample matched{
quantity,
position,
{reference.begin(), reference.end()},
std::move(actual),
tolerance,
};
return {std::move(matched), {}};
}
ComparisonReport compare_samples(
const std::span<const ComparisonSample> samples) {
ComparisonReport report{true, 0.0, {}};
std::set<PositionKey> positions;
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"));
report.maximum_normalized_error =
std::numeric_limits<double>::infinity();
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()));
report.maximum_normalized_error =
std::numeric_limits<double>::infinity();
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"));
report.maximum_normalized_error =
std::numeric_limits<double>::infinity();
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<double>::infinity()));
report.maximum_normalized_error =
std::numeric_limits<double>::infinity();
continue;
}
const double denominator =
sample.tolerance.absolute_scale +
sample.tolerance.relative * std::abs(reference);
if (!(denominator > 0.0) || !std::isfinite(denominator)) {
add_failure(
report.failures,
"validation.invalid_tolerance",
scalar_failure_text(
sample,
component,
std::numeric_limits<double>::infinity()));
report.maximum_normalized_error =
std::numeric_limits<double>::infinity();
continue;
}
const double normalized_error =
std::abs(actual - reference) / denominator;
report.maximum_normalized_error = std::max(
report.maximum_normalized_error,
normalized_error);
if (!(normalized_error <= 1.0)) {
add_failure(
report.failures,
"validation.tolerance_exceeded",
scalar_failure_text(
sample, component, normalized_error));
}
}
}
report.passed = report.failures.empty();
return report;
}
CorrelationReport correlate_samples(
const std::span<const ComparisonSample> 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<PositionKey> positions;
std::map<std::pair<ReferenceQuantity, std::size_t>, 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<double>::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<double>(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
@@ -0,0 +1,305 @@
#include <fesa/io/hdf5/writer.hpp>
#include <fesa/validation/comparison.hpp>
#include <fesa/validation/reference_csv.hpp>
#include <charconv>
#include <cmath>
#include <filesystem>
#include <iomanip>
#include <iostream>
#include <limits>
#include <optional>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
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<fesa::Diagnostic>& diagnostics) {
for (const fesa::Diagnostic& diagnostic : diagnostics) {
print_diagnostic(diagnostic);
}
}
void print_usage() {
std::cerr
<< "Usage:\n"
<< " fesa-reference-compare --results <results.h5> "
"--instance <name> --displacements <displacements.csv> "
"--reactions <reactions.csv> "
"--internal-forces <internal-forces.csv> "
"--displacement-absolute-scale <value> "
"--reaction-absolute-scale <value> "
"--internal-force-absolute-scale <value>\n";
}
std::optional<double> 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<ComparisonRequest> parse_request(
const int argc,
char* argv[]) {
if (argc != 17) {
return std::nullopt;
}
std::optional<std::filesystem::path> results;
std::optional<std::string> instance;
std::optional<std::filesystem::path> displacements;
std::optional<std::filesystem::path> reactions;
std::optional<std::filesystem::path> internal_forces;
std::optional<double> displacement_absolute_scale;
std::optional<double> reaction_absolute_scale;
std::optional<double> 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<fesa::ReferenceRow>& rows,
const fesa::Tolerance tolerance,
std::vector<fesa::ComparisonSample>& 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<ComparisonRequest> 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<fesa::ComparisonSample> 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<double>::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;
}

Some files were not shown because too many files have changed in this diff Show More