diff --git a/docs/business-rules.md b/docs/business-rules.md index db73b5d..98cadb0 100644 --- a/docs/business-rules.md +++ b/docs/business-rules.md @@ -53,6 +53,32 @@ procedure의 검증을 모두 통과한 candidate만 accepted state와 history 시점부터 적용한다. 반력은 열저장까지 포함한 full residual의 constrained projection이다. 중앙면 양쪽 flux trace는 서로 다른 위치이며 평균하지 않는다. +## MITC4 drilling의 선택과 인공 모멘트 + +현재 모델 전체는 한 코드 설정으로 직접 초기 director 스프링 또는 기존 상대 비틀림을 +선택한다. 기본값은 직접 스프링이며 선형·비선형 정적·정적 pressure·frequency·HHT에 +같이 적용한다. Beam·Hex8·DS4에는 적용하지 않는다. Runtime·section별 선택은 현재 없다. + +직접 방식에서 `d0`는 초기 global 단위 director, `phi`는 전체 global rotation vector다. +`gamma=d0·phi`, `Udrill=0.5*kd*sum(gamma²)`, `frot=kd*gamma*d0`, +`Krot=kd*d0*d0^T`다. 병진·절점 간 coupling·drilling force의2차 방향미분은 exact zero다. +Kd는 BC 전 초기 physical rotary 대각8개의 overflow 방지 평균에 `1e-3`을 곱해 고정한다. +Zero 대각도 평균에 포함하며 음수/nonfinite·전부 zero·nonpositive kd는 실패다. +Floor·density 치환·trial 재조정은 없다. 유한 강체 회전의 인공 에너지·모멘트는 허용한다. + +직접 HHT hidden 식은 실제 관성의 normalized null `n(phi)`에 대해 +`sum_e kd*gamma_ei*(n·d_ei0)-n·Fext_rot`이며 같은 trial의 고차 미분을 사용한다. +공유 node의 외력 항은 한 번만 뺀다. 일정 drilling 접선이 null의 고차 미분을 없애지는 않는다. +기존 방식의 계산은 polar1·관성3차 근사와 기존 판정을 보존하며 정확도 실패를 자동 수정하지 않는다. + +내부 spatial `Mdrill=sum[x cross fdrill_trans+Jl(phi)^(-T)*fdrill_rot]`를 source 순서로 +누적한다. 선형은 reference 위치/작은 회전, 비선형/HHT는 현재 위치/기존 SO(3) 변환이다. +기존 raw 평형은 그대로이고 보정 평형은 `[raw_force,raw_moment-Mdrill]`다. +직접 방식의 기존 전역 gate만 보정값을 사용하고 기존 방식은 raw를 사용한다. +기존 scale·exact-zero·허용치는 유지하며 인공 모멘트를 분모에 더하지 않는다. +비선형 정적의 전역 수치 gate는 기존 reference 비교기에 있으며 procedure에 새 gate를 만들지 않는다. +Marker 없는 옛 결과의 방식은 미상이다. 선형/HHT 변위 비교기에는 새 평형 수치 family를 만들지 않는다. + ## 외부 수치 비교 먼저 source identity·row 집합·component·schema·단위·좌표·위치와 finite 여부를 검사한다. diff --git a/docs/contracts.md b/docs/contracts.md index 868b542..392ac9c 100644 --- a/docs/contracts.md +++ b/docs/contracts.md @@ -95,6 +95,40 @@ Schema·필수 field·dtype·shape·inventory·finite 불일치는 읽기나 비 누락을 0·NaN·기본 marker로 채워 이용하지 않는다. Mode shape는 응답 변위가 아니며 frequency와 raw lambda도 다른 양이다. 작은 음의 lambda를 abs/clamp하여 raw 값을 바꾸지 않는다. +## MITC4 drilling 추가 진단 + +기존 feature/procedure/schema 번호를 유지한다. 새 MITC4 구조 결과는 실제 계산한 state의 +방식을 저장하며 writer가 현재 코드 기본값으로 추정하지 않는다. Beam·Hex8·thermal의 +metadata·dataset·report 계약에는 아래 필드를 추가하지 않는다. + +| 위치 | 필드 | 형식과 의미 | +| --- | --- | --- | +| `/metadata` attribute | `mitc4_drilling_formulation` | UTF-8 scalar; `direct-initial-director-spring` 또는 `legacy-corotated-twist` | +| `/metadata` attribute | `drilling_diagnostics_version` | U64LE scalar, 1 | +| `/steps/Step-1/frames//global` | `drilling_moment` | F64LE `[3]`, global spatial 내부 복원모멘트 XYZ, force×length | +| 같은 위치 | `equilibrium_with_drilling` | F64LE `[6]`, raw force XYZ와 raw moment에서 drilling moment를 차감한 XYZ | + +기존 `global/equilibrium`이 있는 선형/비선형 정적·HHT의 모든 frame에 두 dataset을 +추가하며 HHT frame0도 포함한다. Eigen에는 두 metadata attribute만 추가하고 가상 평형 +frame을 만들지 않는다. Raw equilibrium의 기존 값·단위·의미는 유지한다. + +일반 reader의 두 marker가 모두 없는 옛 파일은 옛 계약으로 읽고 진단은 미상/null이다. +한 marker만 있거나 알 수 없는 방식/version·필수 dataset 누락·dtype/shape/finite 오류는 +실패다. 모든 applicable frame에서 보정 force는 raw와 같고 보정 moment는 binary64의 +`raw-Mdrill`과 정확히 같아야 한다. 작은 값에 absolute tolerance를 더해 허용하지 않는다. +새 candidate self-check는 actual state로 marker 필요 여부를 판단하므로 두 marker를 +모두 삭제해도 옛 파일로 오인해 통과시키지 않는다. Write·close·readback·최종 교체 실패는 +이전 최종 HDF5를 보존한다. + +기존 MITC4 JSON report는 schema/version·기존 키를 유지하고 `drilling_diagnostics`를 +추가한다. 옛 markerless 결과는 `null`이다. 새 값은 `formulation` 문자열, `version:1`, +`frames` 배열이다. 각 frame은 `storage_step:"Step-1"`, 정수 `frame_index`, +숫자 배열 `raw_equilibrium[6]`, `drilling_moment[3]`, `equilibrium_with_drilling[6]`을 +가지며 storage frame 번호 순서로 나열한다. 기존 C++ typed report도 같은 정보를 보존한다. +Eigen typed 진단의 frames는 빈 배열이며 별도 eigen JSON API를 새로 만들지 않는다. +비선형 정적의 기존 moment family만 marker에 따라 판정값을 고른다. 선형/HHT 비교기는 +진단을 사전검사하며 기존 family·U blocking/UR warning·수치 verdict를 유지한다. + ## C++ consumer와 비교 도구 공개 interface는 `include/fesa`의 backend-neutral 값과 `Status`/`Result`를 사용한다. diff --git a/docs/engineering-notes.md b/docs/engineering-notes.md index cd5bc26..f7c8160 100644 --- a/docs/engineering-notes.md +++ b/docs/engineering-notes.md @@ -1,5 +1,34 @@ # 재현된 함정과 대응 +## 작은 직접 스프링에서 중간 곱이 먼저 사라진 경우 + +양의 subnormal kd에 먼저 `0.5*kd`를 계산하면 energy의 AD 미분에서 force/tangent가 +사라질 수 있다. 직접 force와 일정 tangent는 명시적 선형식으로 계산하고 energy는 +`0.5*gamma*(kd*gamma)`로 평가한다. 양의 작은 scale을 floor로 올리는 대응은 하지 않는다. +독립 tiny-scale oracle로 energy·force·tangent를 각각 검사했다. + +## 관성 null과 drilling primitive가 달라진 경우 + +Cubic 관성의 실제 null은 정확 SO(3)의 closed null과 다르다. 직접 스프링의 hidden +미분에 후자를 사용하면 실제 관성과 다른 제약을 풀게 된다. 실제 cubic Exp/Jacobian의 +normalized solve와 그 미분을 사용한다. 별도 physical force의 null 투영도 정확히0이라고 +가정할 수 없어 `n^T(full force-external)`로 canonical drilling scalar를 재구성하지 않는다. +Full mass/null·triad·physical residual gate는 그대로 남으며 이 선택이 기존 D3 실패를 해결하지는 않는다. + +## 빈 null과 미평가 scalar를 혼동한 경우 + +방향 응답을 요청하지 않은 assembly의 canonical scalar가 비어 있다는 이유만으로 +all-rotation-fixed/all-fixed 초기화를 거부했던 회귀를 재현했다. Scalar는 실제 admissible +null이 소비할 때 요구한다. 필요한 scalar의 누락/Inf는 계속 실패이며 미평가를0으로 채우지 않는다. +방식이 다른 cache도 typed configuration identity로 구별한다. + +## 새 candidate를 옛 markerless 파일로 오인한 경우 + +호환 reader만으로 새 파일을 검사하면 marker 두 개를 모두 지운 candidate가 옛 파일로 +읽힐 수 있다. Writer self-check는 actual state의 진단 존재를 별도로 요구한다. +양 marker/양 dataset 삭제와 write/readback close 실패 주입에서 이전 final bytes 보존을 확인했다. +일반 옛 파일 reader의 null 진단은 유지하며 형식 검사를 위한 tolerance floor를 추가하지 않는다. + ## 긴 임시 경로가 출력 충돌처럼 보인 경우 Integration 검증에서 TMP 아래 lease 경로가 260/261자가 되어 `output-busy`가 발생했다. diff --git a/docs/hht-implicit-dynamics/README.md b/docs/hht-implicit-dynamics/README.md index 03faa99..d6209b4 100644 --- a/docs/hht-implicit-dynamics/README.md +++ b/docs/hht-implicit-dynamics/README.md @@ -2,6 +2,22 @@ ## 범위와 현재 상태 +2026-10-20 현재 drilling 기본값은 `gamma_i=d_i0·phi_i`의 직접 작은 스프링이며 +기존 상대 비틀림도 모델 전체의 코드 설정으로 선택한다. 직접 경로는 drilling용 +polar/swing/twist를 생략하고 physical 회전·관성·고정 null 기저·시간적분은 유지한다. +직접 hidden scalar는 실제 cubic 관성이 계산한 normalized null `n_i(phi)`에 대해 +`sum_e kd*gamma_ei*(n_i·d_i0)-n_i·Fext_rot`다. 그 같은 식의 gradient·Hessian· +필요한 3차 방향미분을 사용한다. 초기 director를 현재 null로 고정하지 않는다. + +현재 source `6bbbb2b`를 전체 실행한 Debug Native/CTest는 각각1035개 중 +1024 PASS·기존10 FAIL·1 INCONCLUSIVE SKIP이며 별도 Python164개는 PASS다. +정식729node·101frame Direct Release 해석과 외부9node U의101 families 비교는 +CLI0/비교0, 최대 정규화 RMS0.2282326248391787%로 통과했고 UR 경고200개를 유지했다. +기존 방식은 별도 Debug의 관련 검사와 실제 소형3frame HHT 출력을 확인했다. +729node 기존 방식의 외부 U 재해석은 이번에 하지 않았다. 과거 결과를 현재 PASS로 승계하지 않는다. +기존10개 물리 실패·N7 미판정·정식 자료 누락이 남으므로 전체 검증 완료는 아니다. +아래 2026-10-19 수용 기록은 이전 수렴 실험의 결정이며 현재 전체 PASS를 뜻하지 않는다. + 2026-10-19 후속 사용자 승인으로 이번 수렴 규칙 실험은 **Abaqus 변위 검증 통과·실패10개 보류로 조건부 완료**한다. [수용 범위](../tracking/decisions/0011-tolerance-experiment-acceptance.md)를 따르며, 아래 미해결 물리 검사· N7와 실제 FAIL/SKIP은 보존한다. 전체 상태729절점·101프레임 진단과 관측 시간 비교는 @@ -76,7 +92,7 @@ Schema identity는 `(hht-implicit-dynamics,hht-implicit-dynamics,2)`다. `reduced-coordinate-HHT-with-physical-rate-reconstruction`, `mass_formulation=nodal-lumped-kinetic`, `acceleration_convention=physical-endpoint-and-auxiliary-separated`다. -모든 result dataset version도 2이며 A의 quantity_semantics는 +기존 schema2 result dataset의 version은2이며 A의 quantity_semantics는 `full-chart-HHT-filter-combination-not-predictor-driver`다. 잘못된 marker·혼합 version·mass· inventory·nonfinite는 dynamic-output-schema 실패다. Historical v1 reader branch는 과거 full-coordinate 결과만 읽으며 marker를 보충하거나 v2로 재표시하지 않는다. @@ -92,6 +108,13 @@ inventory·nonfinite는 dynamic-output-schema 실패다. Historical v1 reader br | nodal/reaction | [N,6]; support-on-structure spatial RF/RM, moment는 Jl^-T 변환 | | element/shell/local_frame, generalized_strain, section_resultant, stress | [E,4,3,3], [E,4,8], [E,4,8], [E,4,3,3]; stress는 bottom/middle/top × S11/S22/S12 | | global/energy, global/equilibrium | [5] kinetic/physical/drill/work/defect, [6] 관성 포함 force/moment balance | +| global/drilling_moment, global/equilibrium_with_drilling | [3] 내부 spatial drilling moment, [6] raw force와 raw moment에서 drilling moment를 차감한 진단; frame0 포함 | + +새 MITC4 출력에는 실제 방식의 UTF-8 `mitc4_drilling_formulation`과 U64LE +`drilling_diagnostics_version=1`을 metadata에 추가한다. 기존 schema2와 raw equilibrium은 +유지한다. 직접 방식의 기존 전역 gate만 보정 moment를 사용하고 기존 방식은 raw 판정이다. +Marker 없는 옛 결과는 방식 미상이며 새 진단을 추정하지 않는다. 변위 비교기는 새 진단의 +모든 frame을 사전검사하지만 새 평형 수치 family나 verdict를 만들지 않는다. Row는 model source inventory와 결합한다. Generalized moment와 spatial moment, physical a와 A, drilling force와 physical stress를 혼합하지 않는다. 영 BC는 q/v/a/A exact zero다. diff --git a/docs/hht-implicit-dynamics/formulation-fixed-null-basis.md b/docs/hht-implicit-dynamics/formulation-fixed-null-basis.md index 7c9600f..75437e3 100644 --- a/docs/hht-implicit-dynamics/formulation-fixed-null-basis.md +++ b/docs/hht-implicit-dynamics/formulation-fixed-null-basis.md @@ -16,12 +16,18 @@ component selector를 아래 고정 ON 기저로 대체한다. 기존 ## 1. 적용 가정과 물리 계약 Total Lagrangian MITC4, 큰 변위·유한회전·작은 탄성 변형률, 등방 선형 탄성, -objective drilling potential, `nodal-lumped-kinetic`, 기존 공간 적분·두께 적분과 +선택한 drilling potential, `nodal-lumped-kinetic`, 기존 공간 적분·두께 적분과 fixed-initial total rotation vector를 유지한다. 재료 history update는 없다. 관성·내력·공간 각속도·단면력은 [기존 물리 정식화](formulation.md)와 [관성 상세](formulation-inertia.md)가 소유한다. 좌표 변경은 constitutive law, mass measure, drilling stiffness 또는 recovery formula를 바꾸지 않는다. +2026-10-20 기본 직접 스프링은 별도로 승인한 numerical potential 변경이다. +고정 기저·partial-BC 교집합·graph/rate 복원·regularity 계약은 유지하며, hidden primitive만 +직접 힘과 실제 cubic 관성 null에서 일관되게 계산한다. 기존 objective twist를 선택하면 +그 현행 계산과 판정을 보존한다. 현재 null을 초기 director로 고정하거나 partial BC를 +full clamp로 바꾸어 새 스프링의 rank/regularity 실패를 피하지 않는다. + Full node-major ordering은 `q_i=(u_x,u_y,u_z,phi_X,phi_Y,phi_Z)`다. `u`: length, `phi`: rad, `v=dot q`, `a=ddot q`이며 `v_phi`와 `a_phi`는 공간 각속도·각가속도가 아니다. Zero fixed BC를 제거한 공간에서 diff --git a/docs/hht-implicit-dynamics/formulation-inertia.md b/docs/hht-implicit-dynamics/formulation-inertia.md index fcfdbbb..b367d16 100644 --- a/docs/hht-implicit-dynamics/formulation-inertia.md +++ b/docs/hht-implicit-dynamics/formulation-inertia.md @@ -10,6 +10,14 @@ ## 1. 회전과 관성의 미분 +2026-10-20 drilling 선택 변경은 아래 관성의 수학 유도나 현행 cubic 근사 차수를 바꾸지 +않는다. 기본 직접 스프링의 hidden 미분은 실제 cubic 관성 null을 사용한다. 아래 정확 +SO(3) 유도와 현재 근사 구현의 D3/null 항등식 검증 상태를 구분한다. +직접 스프링은 유한 강체 회전에서 인공 내부 spatial moment `Mdrill`을 허용하므로, +아래 objective total energy의 무일 가정은 직접 drilling에 적용하지 않는다. +Raw 관성 포함 평형을 보존하고 직접 방식의 기존 전역 gate에는 moment에서 `Mdrill`을 +차감한다. 기존 방식은 raw 판정이며 physical mass·momentum·a/A·history 정의는 그대로다. + `P=[phi]_×`, `theta=|phi|`에 대해 \[ diff --git a/docs/hht-implicit-dynamics/formulation-reduced-coordinate.md b/docs/hht-implicit-dynamics/formulation-reduced-coordinate.md index 8aa87b4..2de3b57 100644 --- a/docs/hht-implicit-dynamics/formulation-reduced-coordinate.md +++ b/docs/hht-implicit-dynamics/formulation-reduced-coordinate.md @@ -15,7 +15,16 @@ physics와 graph/rate 계약은 유지하되, identity-coordinate 전용 설명 그대로 적용하지 않는다. 실제 augmented scaling·초기화·CA/N7의 새 basis 경계는 해당 상세 문서를 따른다. -기존 Total Lagrangian MITC4, objective numerical drilling potential, isotropic linear elasticity, nodal lumped kinetic measure, 공간 적분과 fixed-initial global chart를 유지한다. 좌표는 node-major `q=(u_x,u_y,u_z,phi_X,phi_Y,phi_Z)`, `v=dot q`, `a=ddot q`; 병진 단위는 length, 회전은 rad다. 질량 추가, damping, stiffness modification, stress averaging은 없다. 공개 load subset은 기존 constant global force, 구속은 fixed zero다. 아래 smooth F(t)의 시간미분은 C08a 내부 manufactured oracle에만 해당한다. +기존 Total Lagrangian MITC4, 선택한 numerical drilling potential, isotropic linear elasticity, nodal lumped kinetic measure, 공간 적분과 fixed-initial global chart를 유지한다. 좌표는 node-major `q=(u_x,u_y,u_z,phi_X,phi_Y,phi_Z)`, `v=dot q`, `a=ddot q`; 병진 단위는 length, 회전은 rad다. 질량 추가, damping, stress averaging은 없다. 공개 load subset은 기존 constant global force, 구속은 fixed zero다. 아래 smooth F(t)의 시간미분은 C08a 내부 manufactured oracle에만 해당한다. + +2026-10-20 기본 직접 스프링의 primitive는 +`psi_ei=kd*(d_i0·phi_i)*(n_actual(phi_i)·d_i0)`다. `n_actual`은 기존 +`EvaluateDynamicRotation`의 cubic Exp/Jacobian에서 계산한 `normalize(J_cubic^-1*d_cubic)`를 +사용하며, 아래 정확 SO(3)의 closed 식으로 대체하지 않는다. Shared node에서 요소별 +primitive를 source 순서로 합하고 `n_actual·Fext_rot`을 한 번만 뺀다. +Gradient·Hessian·3차 방향미분도 같은 primitive와 실제 null에서 유도한다. +Metric은 `sum_e kd*grad(gamma_ei)*grad(gamma_ei)^T`다. Force의 일정 접선을 이유로 +null의 고차 미분을 지우지 않는다. 다음 `k*S*gamma`는 기존 상대 비틀림 방식의 유도다. Zero constrained 성분을 정확히 제거한 free 공간에서 @@ -32,7 +41,14 @@ S_i=1/\|J_l^{-1}d_i\|. \tag{RC-02} \] -이는 기존 `drilling_constraint`의 scalar이며 raw twist gamma가 아니다. Physical force의 null work가 0인 정확 identity를 사용하므로 energy의 4차 tensor는 필요 없다. Translation force의 null work는 0이다. Shared node external term은 한번만 추가한다. +이는 기존 `drilling_constraint`의 scalar이며 raw twist gamma가 아니다. 이 유도는 physical force의 null work가 0인 정확 identity를 사용하므로 energy의 4차 tensor는 필요 없다. Translation force의 null work는 0이다. Shared node external term은 한번만 추가한다. + +현재 cubic 관성에서는 별도로 계산한 physical force의 null 투영이 이 정확 identity를 +항상 만족하지 않는다. 직접 방식의 hidden scalar를 `n^T(fphysical+fdrill-Fext)`로 +재구성하면 primitive 미분과 다른 식이 되므로 canonical scalar를 그대로 사용한다. +Full mass/null·physical residual·triad acceptance는 별도로 유지하며 근사 오차를 숨기지 않는다. +빈 admissible null inventory에는 scalar를 요구하지 않고, 필요한 null이 있을 때 미평가 +scalar를0으로 치환하거나 사용할 수 없는 directional 응답을 재사용하지 않는다. ## 2. Fixed graph와 부분 구속 diff --git a/docs/linear-static-mitc4-shell/README.md b/docs/linear-static-mitc4-shell/README.md index 3ffc967..67161d4 100644 --- a/docs/linear-static-mitc4-shell/README.md +++ b/docs/linear-static-mitc4-shell/README.md @@ -11,6 +11,13 @@ Composite·offset·material nonlinearity·mixed beam-shell은 이 기능 범위 ## 채택한 수식과 이유 +2026-10-20부터 모델 전체의 drilling 방식은 코드의 단일 설정으로 선택한다. +기본 `kDirectSpring`은 초기 global 단위 director `d0`와 전체 회전 `phi`에 대해 +`gamma=d0·phi`, `Udrill=0.5*kd*sum(gamma²)`를 사용한다. +`kLegacyCorotatedTwist`는 기존 상대 비틀림 방식이다. 선형 정적과 원형상 고유치의 +drilling block은 두 선택에서 같으며, 차이를 만들기 위해 physical 강성을 바꾸지 않는다. +직접 방식은 유한 강체 회전의 인공 에너지·모멘트를 허용한다. + 각 절점의 6 global DOF에서 director 접평면의 5 physical DOF를 분리한다. Bilinear geometry와 canonical edge-midpoint tying으로 공변 횡전단장을 보간한다. 이는 직접 shear interpolation의 locking을 줄이기 위한 MITC 선택이다. @@ -48,6 +55,11 @@ Key는 frame·instance·source node·component이며 row-set·finite·unit/locat 그 밖은 relative error `<=5%`, family `RMS(error)/S<=1%`를 동시에 요구한다. `S=0`은 exact zero다. UR 수치 실패는 warning이며 잘못된 identity를 허용하는 의미는 아니다. 내부 free residual·전역 force/moment balance는 normalized `1e-10` 한계다. +직접 방식의 기존 전역 gate는 raw moment에서 별도 내부 spatial drilling moment를 +차감한 값을 사용한다. 기존 방식은 raw 판정을 유지하며 분모·허용치는 같다. +기존 `global/equilibrium`은 raw 값이다. 새 파일에는 실제 방식과 진단 version, +`global/drilling_moment[3]`, `global/equilibrium_with_drilling[6]`을 추가한다. +선형 변위 비교기에 새 평형 수치 family를 추가하지 않고 진단의 형식·관계만 사전검사한다. ## 구현과 데이터 흐름 @@ -61,6 +73,11 @@ Procedure는 해석 상태, writer는 serialization을 소유하며 element에 e ## 검증 의미와 한계 +2026-10-20 source `6bbbb2b`의 두 방식 별도 실행에서 정식 S4 입력의 생성·U 비교를 +통과했고 실제 선형 변위도 정확히 같았다. 직접 방식의 인공 torque·보정 부호·실제 +불평형 거부 검사를 포함했다. 전체 Debug 검사는 Native/CTest 각각 +1035개 중1024 PASS·기존10 FAIL·1 INCONCLUSIVE SKIP으로 전체 수치 검증 통과는 아니다. + 과거 S4 reference는 147 blocking U rows와 평형·대칭·constitutive·physical energy 검사를 통과했다. 그때 중심 U3는 약 `-2.35504418e-5`, normalized free residual은 약 `1.08e-14`였다. 이는 당시 실행 evidence이며 이번 문서 정리에 따른 현재 build/test 통과 선언이 아니다. diff --git a/docs/linear-static-mitc4-shell/formulation.md b/docs/linear-static-mitc4-shell/formulation.md index ecb787b..3aa2574 100644 --- a/docs/linear-static-mitc4-shell/formulation.md +++ b/docs/linear-static-mitc4-shell/formulation.md @@ -5,7 +5,14 @@ > > 적용 경계: 물리 정식화와 당시 수치 stabilization 유도를 구분한다. 현재 공변 trace-average > scale·frozen lifecycle은 [고유치 질량·stabilization](../mitc4-shell-feast-eigen2/formulation-mass-stabilization.md)을 -> 함께 따른다. 이전 minimum-diagonal 식을 현재 기준으로 재도입하지 않으며 S3 partial-BC 보류를 유지한다. +> 함께 따른다. 이전 minimum-diagonal 식을 현재 기준으로 재도입하지 않는다. +> +> 2026-10-20 채택: 단일 코드 설정의 기본값은 초기 director 직접 스프링이며 기존 상대 +> 비틀림도 선택한다. 선형/eigen의 drilling block과 trace-average scale은 두 방식에서 같다. +> 직접 방식의 기존 전역 moment gate는 raw에서 별도 내부 drilling moment를 차감한다. +> 아래 초기 문서의 drilling 전용 출력 금지와 objective nonlinear 후보 설명은 당시 범위다. +> 현재는 실제 방식 metadata와 raw/인공/보정 평형을 별도로 저장하며, 직접 방식의 유한 +> 강체 회전 인공 에너지를 허용한다. 물리 MITC strain·constitutive·적분 유도는 유지한다. ## Metadata diff --git a/docs/mitc4-shell-feast-eigen2/README.md b/docs/mitc4-shell-feast-eigen2/README.md index d56b8aa..0569305 100644 --- a/docs/mitc4-shell-feast-eigen2/README.md +++ b/docs/mitc4-shell-feast-eigen2/README.md @@ -15,6 +15,9 @@ Mass는 `3×3×3` 적분의 scalar thickness moments를 row-sum하여 nodal `6× 3 병진·2 director 접평면 회전 관성과 warped geometry cross coupling을 보존한다. Numerical drill은 analytic mass null이며 인공 drilling inertia를 넣지 않는다. Stabilization scale은 초기 physical rotary diagonal 8개 trace 평균의 `10^-3`이다. +2026-10-20부터 모델 전체의 코드 설정은 직접 초기 director 스프링을 기본으로 하며 +기존 상대 비틀림도 선택할 수 있다. 원형상 eigen의 K/M과 drilling block은 두 방식에서 +같다. Artificial drilling mass를 추가하거나 서로 다른 spectrum을 만들지 않는다. Global 실제 BC를 먼저 적용하고 analytic null과 free 공간의 교집합만 제거한다. Partial rotational BC를 full clamp로 대체하지 않고 finite constrained spectrum을 유지한다. Positive nodal mass space를 whitening하고 sparse symmetric 문제를 FEAST에 전달한다. @@ -41,6 +44,9 @@ HDF5는 `schema_version=0`, `eigen_schema_version=1`, procedure `linear-frequenc Constrained 성분은 exact zero이고 sign은 최대 scaled global component의 stable rule로 정한다. Raw 작은 음수 lambda를 clamp하지 않는다. Frequency는 lambda>=0 mode subset만 생성한다. `verification`에는 전체 c개 residual·Gram·null 검사와 backend call/control evidence를 보존한다. +실제 선택은 `/metadata`의 `mitc4_drilling_formulation`과 U64LE +`drilling_diagnostics_version=1`에 기록한다. Eigen에는 평형 frame이나 가상 drilling +moment dataset을 만들지 않는다. Typed 비교 report의 진단 frames도 빈 배열이다. | 원본 | 비교 계약 | | --- | --- | @@ -69,10 +75,16 @@ Domain→local K/M→ordered assembly→reduced solve→global phi 복원→all- ## 실제 evidence와 미해결 범위 +2026-10-20 source `6bbbb2b`의 두 방식 별도 실행에서 정식10-mode 생성·비교를 통과했고 +실제 eigen table은 정확히 같았다. 전체 Debug Native/CTest는 각각1035개 중 +1024 PASS·기존10 FAIL·1 INCONCLUSIVE SKIP이다. 해당 Value 비교 통과와 전체 물리 검증을 구분한다. + 2026-09-10 사용자 승인 첫 구현에서 전체 417/417 PASS와 원본 10 Value의 모드별5% 통과 기록이 있다. 최대 외부 상대차는 약 `4.261331%`다. 이전 RMS1% 불일치 이력은 소급 PASS로 고치지 않는다. 이 기록은 당시 revision의 evidence이며 이번 문서 정리는 새 build/test를 실행하지 않았다. Rational mass oracle·Jensen energy·좌표 재계산 공변성·density scaling·all-c 실패 주입이 보완 근거다. -S3 public static 평형 충돌과 Abaqus Value 정확도 개선은 후속 보류다. +S3 partial-BC의 public static은 직접 방식에서 인공 모멘트 차감 후 기존 gate를 통과한다. +기존 방식의 raw gate 실패는 명시적 검사로 보존했다. 이는 partial BC를 full clamp로 +바꾸거나 Value 허용치를 완화한 해결이 아니다. Abaqus Value 정확도 개선은 후속 보류다. 실제 512회 vendor 호출 또는 실제 HDF5 close-error의 모든 조합을 주입한 것은 아니다. Retained shell/shellNL 회귀와 내부 검증 의무를 유지하며 `reference/temp` 제외 예제를 재사용하지 않는다. diff --git a/docs/mitc4-shell-feast-eigen2/formulation-mass-stabilization.md b/docs/mitc4-shell-feast-eigen2/formulation-mass-stabilization.md index d554302..8acdfbc 100644 --- a/docs/mitc4-shell-feast-eigen2/formulation-mass-stabilization.md +++ b/docs/mitc4-shell-feast-eigen2/formulation-mass-stabilization.md @@ -136,7 +136,14 @@ spectrum·static response의 공변성을 검증한다. 다른 BC를 같은 문 여섯 물리 rigid modes는 양의 mass를 갖는다. Pure numerical drill `[0;di]`는 mass=0, stiffness>0이다. Linear UR에 full spatial Omega를 넣는 field와 이 physical 대표값을 구분한다. -## S2. 비선형 objective potential과 고정 scale +## S2. 비선형 drilling 선택과 고정 scale + +2026-10-20 현재 기본값은 초기 global 단위 director `d_i0`를 사용한 직접 스프링이다. +`gamma_i=d_i0·phi_i`, `Pi_d=0.5*kd*sum(gamma_i²)`, +`r_i=kd*gamma_i*d_i0`, `K_ii=kd*d_i0*d_i0^T`이며 병진·절점 간 coupling은0이다. +Drilling polar/swing/atan2를 평가하지 않고 유한 강체 회전의 인공 모멘트를 허용한다. +S1의 초기 scale과 mass null은 유지한다. 원형상 선형/eigen block은 두 방식에서 같다. +다음 objective 유도는 선택 가능한 기존 상대 비틀림 방식에 한정한다. 기존 [nonlinear 정식화 §10](../nonlinear-static-newton-raphson/formulation.md)의 polar corotation, total nodal rotation, minimal swing 및 objective `gamma_i=atan2(s_i,c_i)`를 @@ -148,7 +155,9 @@ corotation, total nodal rotation, minimal swing 및 objective `gamma_i=atan2(s_i K_d=k_d\left(D_\gamma^TD_\gamma+\sum_i\gamma_iD^2\gamma_i\right). \tag{S2} \] -Polar/swing/atan2의 미분과 branch validation을 생략하지 않는다. Trial마다 scale을 바꾸지 않는다. +기존 방식에서는 polar/swing/atan2의 미분과 branch validation을 생략하지 않는다. +현행 구현의 polar1과 미분 정확도 실패를 보존하며 과거 polar12로 복원하지 않는다. +두 방식 모두 trial마다 scale을 바꾸지 않는다. 일반 nonlinear Hessian 전체의 SPD는 주장하지 않는다. 초기에도 Dgamma는 corotation coupling을 포함하므로 linear Td와 동일시하지 않는다. 기존 physical 선형극한, directional residual/Hessian, objective rigid orbit·state rollback을 검증한다. 새 eigen에는 이 nonlinear state를 전달하지 않는다. @@ -173,3 +182,7 @@ Static 대응 load는 node1의 tangent moment `m=(5/3)M0 b`, 다른 load=0이다 `d·m=0,mZ=M0`; `theta1Z=M0/Kff`, full residual `r=Kq-F`다. 회전 뒤 `m'=Qm`을 사용한다. Pure global eZ moment는 drilling 성분 때문에 기존 load 계약에 부합하지 않는다. Eigen의 constrained residual을 static reaction으로 출력하지 않는다. + +2026-10-20 현재 이 static 사례의 직접 방식은 별도 인공 spatial 모멘트를 차감한 기존 +gate를 통과하고, 기존 방식은 raw gate에서 계속 거부된다. S3의 K/M/lambda target이나 +partial BC를 바꾼 결과가 아니며 Value 비교 예외를 정적 평형에 적용하지 않는다. diff --git a/docs/nonlinear-static-newton-raphson/README.md b/docs/nonlinear-static-newton-raphson/README.md index e80d2ae..463ef56 100644 --- a/docs/nonlinear-static-newton-raphson/README.md +++ b/docs/nonlinear-static-newton-raphson/README.md @@ -2,6 +2,11 @@ ## 목적과 경계 +2026-10-20 현재 drilling 기본값은 초기 director의 직접 작은 스프링이다. +모델 전체의 한 코드 설정으로 기존 상대 비틀림도 선택하며 정적 pressure에도 같은 +선택을 적용한다. 직접 방식은 drilling의 polar/swing/twist 계산을 생략한다. +아래 polar 근사 실험의 실패는 보존한 기존 방식의 근거이며 직접 스프링 식을 뜻하지 않는다. + 2026-10-18 후속 사용자 지시의 [local 근사 실험](../abaqus-tolerance/local-kernel-experiment.md)에서 공유 nonlinear MITC4의 polar 반복을 12→1회로 변경했다. 후속 사용자 지시로 관련 접선 테스트 오차를 1e-3으로 완화한 뒤에도 drilling constraint의 gradient 항등식 검사가 @@ -20,8 +25,12 @@ 초기 director/frame에서 같은 trial의 기하·strain·internal force·consistent tangent를 계산한다. Rotation-vector generalized moment와 global spatial moment는 일의 짝이 다르다. 반력의 공개 moment 복구에서 chart 변환을 적용하며 성분명 변경으로 대체하지 않는다. -Physical MITC strain·constitutive response와 objective numerical drilling potential을 분리한다. +Physical MITC strain·constitutive response와 선택한 numerical drilling potential을 분리한다. Drilling scale은 초기 상태에서 정하고 trial마다 다시 tuning하지 않는다. +직접 방식은 `gamma_i=d_i0·phi_i`, `f_i=kd*gamma_i*d_i0`, +`K_ii=kd*d_i0*d_i0^T`다. 병진과 절점 간 drilling coupling은 exact zero다. +유한 강체 회전에도 인공 복원모멘트가 생기는 것을 허용한다. 기존 방식의 objective +유도와 polar1 구현을 이 직접 방식에 적용하지 않는다. 평형 residual은 internal minus external이며 follower load가 있으면 tangent는 `Kint-dFext/dq`다. Unsymmetric tangent를 대칭화하지 않는다. @@ -58,6 +67,11 @@ Accepted free residual은 matching external scale이 nonzero이면 infinity-norm exact-zero이면 typed absolute `<=1e-3` force 또는 force-length다. Factor-one 전역 force/moment balance는 각각 normalized `<=1e-8`이다. 이는 기존 comparator의 승인 계약이며 두 양을 하나의 `1e-10` 기준으로 요약하지 않는다. Linear 정적·HHT·eigen의 별도 gate는 바꾸지 않는다. +이 전역 수치 판정은 기존 reference 비교기의 책임이다. Procedure에 새 전역 gate를 +만들지 않는다. 비교기는 직접 방식 marker일 때 기존 moment family에 +`raw_moment-Mdrill`을, 기존 방식·marker 없는 옛 파일에는 raw 값을 사용한다. +`Mdrill`은 현재 위치의 병진 모멘트팔과 `Jl(phi)^(-T)`로 변환한 내부 drilling +모멘트의 source 순서 합이다. 분모에 인공 모멘트를 더하지 않는다. ## 구현과 상태 수명 @@ -73,8 +87,14 @@ Domain→non-owning model→stable full/free DOF→local trial response→ordere ## 검증에서 배운 점과 남은 경계 +2026-10-20 source `6bbbb2b`의 기본 직접 방식과 별도 기존 방식에서 정식 coarse +nonlinear 입력의 생성·기존 비교를 통과했다. 두 방식의 저장·reader·원자적 실패 보존도 +검사했다. 보존한 기존 방식의 gradient 항등식은 차이 약0.00329135로 계속 실패하며, +전체 Native/CTest 각각1035개는1024 PASS·10 FAIL·1 INCONCLUSIVE SKIP이다. +이는 현재 실행 결과이며 아래 과거 기록이나 coarse U 통과로 모든 비선형 검증을 대체하지 않는다. + 기존 nonlinear reference와 rollback·finite rotation·tangent·equilibrium 검증 기록이 있다. -이는 새 revision의 build/test 통과를 의미하지 않으며 이번 작업은 문서 정리만 수행한다. +과거 기록만으로 새 revision의 build/test 통과를 의미하지 않는다. Reaction은 단일 global 합만 맞춰 local reaction 정확성을 대체할 수 없다. 큰 자세 변화와 원점 이동에서 spatial moment balance를 함께 확인해야 한다. 외부 U 통과는 warning UR 또는 review RF/RM의 수치 일치를 뜻하지 않는다. diff --git a/docs/nonlinear-static-newton-raphson/formulation.md b/docs/nonlinear-static-newton-raphson/formulation.md index 60ded28..3a5c87c 100644 --- a/docs/nonlinear-static-newton-raphson/formulation.md +++ b/docs/nonlinear-static-newton-raphson/formulation.md @@ -5,7 +5,40 @@ > > 적용 경계: 물리 정식화와 당시 수치 stabilization 유도를 구분한다. 현재 공변 trace-average > scale·frozen lifecycle은 [고유치 질량·stabilization](../mitc4-shell-feast-eigen2/formulation-mass-stabilization.md)을 -> 함께 따른다. 이전 minimum-diagonal 식을 현재 기준으로 재도입하지 않으며 S3 partial-BC 보류를 유지한다. +> 함께 따른다. 이전 minimum-diagonal 식을 현재 기준으로 재도입하지 않는다. +> S3 partial-BC는 직접 방식의 보정 gate와 기존 방식의 raw gate를 구분한다. + +## 2026-10-20 drilling 채택 범위 + +기본값은 초기 global 단위 director `d_i0`와 total global rotation vector `phi_i`의 +직접 스프링이다. 아래 §10의 corotation·minimal swing·atan2 유도는 선택 가능한 +기존 상대 비틀림 방식의 원문으로 보존한다. 기존 구현은 기준 source +`52e321e15c70377d4318aace3f3bc64b3ea1a65c`의 polar1을 유지하며, 원문의 정확한 +객관성·미분 항등식이 이 근사 구현에서 모두 통과했다는 뜻은 아니다. + +```text +gamma_i = d_i0 · phi_i +kd = 1e-3 * average(initial physical K20 tangent-rotation diagonals, 8 entries) +Udrill = 0.5 * kd * sum_i gamma_i^2 +f_drill,i,rot = kd * gamma_i * d_i0 +K_drill,ii,rot,rot = kd * d_i0 * transpose(d_i0) +D² f_drill[v,v] = 0 +``` + +병진 force/tangent와 서로 다른 절점의 coupling은 exact zero다. Global URZ 또는 +현재 director를 `d_i0` 대신 사용하지 않는다. Scale은 BC 전 초기 physical K20에서 +기존 결정적 overflow 방지 평균으로 한 번 계산하며 zero 대각도 8개 평균에 포함한다. +음수/nonfinite 대각, 모두 zero인 대각, nonfinite/nonpositive kd는 실패다. +Clamp·floor·density 치환이나 trial별 재조정은 없다. + +직접 경로는 drilling의 polar/swing/twist/atan2를 평가하지 않는다. Physical strain· +geometry·관성 계산은 유지한다. 초기 director와 전체 회전 기준이므로 유한 강체 회전의 +인공 에너지·복원모멘트를 허용하며 아래 기존 방식의 self-equilibrated 성질을 요구하지 않는다. +내부 generalized drilling 힘은 현재 위치의 모멘트팔과 `Jl(phi)^(-T)` 변환으로 spatial +`Mdrill`을 만들고 source 순서로 누적한다. Raw 평형은 보존하고 직접 방식의 기존 moment +판정에만 `raw_moment-Mdrill`을 사용한다. 기존 방식은 raw 판정을 유지한다. +비선형 정적 procedure에는 원래 진단만 있었으며 수치 gate는 기존 reference 비교기에 있다. +새 procedure gate나 인공 모멘트를 포함한 정규화 분모를 만들지 않는다. ## 1. Metadata diff --git a/docs/operations.md b/docs/operations.md index 20ffee6..7ba5f24 100644 --- a/docs/operations.md +++ b/docs/operations.md @@ -72,6 +72,16 @@ foreach ($test_name in 'fesa_unit_tests','fesa_integration_tests','fesa_referenc 생성된 경로에서도 실제 최장 lease 경로를 확인한다. 너무 길면 승인된 짧은 작업 위치를 사용한다. +## MITC4 drilling 방식 변경 + +`src/fesa/elements/mitc4_drilling_internal.h`의 `kProductionDrillingMode` 한 곳을 바꾼다. +기본값은 `Mitc4DrillingMode::kDirectSpring`, 기존 방식은 +`Mitc4DrillingMode::kLegacyCorotatedTwist`다. 변경 후 solver·CLI·검사를 모두 rebuild한다. +같은 모델의 모든 MITC4 정적·pressure·frequency·HHT에 적용되며 CLI/INP/환경변수나 +section별 override는 현재 없다. 시험용 명시적 alternate mode 경로는 production 설정이 아니다. +서로 다른 mode의 실행물·DLL·입력 sibling 결과를 별도 위치에 두고 실제 metadata marker로 +선택을 확인한다. 선형 변위와 원형상 eigen값이 두 방식에서 같아도 정상이다. + ## CMake 보조 빌드 같은 VS x64 환경에서 로컬 CONFIG package 환경변수를 설정한 뒤 실행한다. diff --git a/docs/shell-surface-pressure-load/README.md b/docs/shell-surface-pressure-load/README.md index f589571..ca8002e 100644 --- a/docs/shell-surface-pressure-load/README.md +++ b/docs/shell-surface-pressure-load/README.md @@ -41,6 +41,11 @@ Follower 여부는 procedure로 결정하며 keyword 옵션을 조용히 무시 결과는 기존 linear schema0 또는 nonlinear schema1의 displacement/reaction·shell recovery다. 별도의 가상 pressure stress 또는 nodal averaged stress를 추가하지 않는다. Nonlinear 비교는 exact factor-one accepted frame만 사용한다. +MITC4의 직접 스프링 기본값과 기존 상대 비틀림 선택은 pressure에도 모델 전체로 적용한다. +Pressure force/tangent 정의는 그대로다. 새 방식 marker가 있는 nonlinear 비교의 기존 +moment family만 직접 방식의 `raw-Mdrill` 또는 기존 방식의 raw를 선택한다. +기존 scale·`1e-8` 기준과 U/UR 정책은 유지한다. 정식 pressure bundle 부재로 +실제 외부 pressure 비교는 아직 미완료이며 내부 pressure 검사가 이를 대신하지 않는다. | 원본 bundle | 비교 의미 | | --- | --- | diff --git a/docs/solver-core/README.md b/docs/solver-core/README.md index a08195d..f0b9db6 100644 --- a/docs/solver-core/README.md +++ b/docs/solver-core/README.md @@ -45,7 +45,7 @@ Named exception **`mitc4-eigen-first-implementation`**(2026-09-10 승인)은 eig Configure/build/test의 compiler·구성·build tree를 함께 맞추고 discovery0건·다른 compiler를 실패로 다룬다. MSVC runtime DLL staging을 확인하며 개발자 PATH의 우연한 성공을 배포 검증으로 삼지 않는다. Dense allocation 크기 곱 overflow는 할당 전에 검사한다. Header는 self-contained `.h`/guard, 파일·namespace snake_case, 타입·함수 PascalCase, 변수 snake_case, class member trailing `_`, Google 기반2-space format을 따른다. Vendor target에 FESA warning 정책을 강제하지 않는다. -현재 공백은 shell HHT의 최신 전체 검증·시간수렴, thermal 새 필수 reference, eigen S3 partial-BC 정적 충돌, 성능 목표 미달이다. 각 기능에 구체 한계를 보존한다. 기존1×1 PARDISO 반복의 heap failure는 SpMV 없는 baseline에서도 재현됐고 vendor 원인은 미확정이다. 기본 정책의 해당 크기에만 local1-thread/TLS 복원을 적용했으며 보조100회 성공을 라이브러리 전체 해결로 확대하지 않는다. +현재 공백은 shell HHT의 전체 물리·시간정확도, thermal 새 필수 reference, 기존 drilling 방식의 S3 partial-BC raw 정적 gate 한계, 성능 목표 미달이다. S3의 기본 직접 방식은 별도 인공 모멘트 차감으로 기존 gate를 통과하며 기존 방식의 raw 거부를 유지한다. 각 기능에 구체 한계를 보존한다. 기존1×1 PARDISO 반복의 heap failure는 SpMV 없는 baseline에서도 재현됐고 vendor 원인은 미확정이다. 기본 정책의 해당 크기에만 local1-thread/TLS 복원을 적용했으며 보조100회 성공을 라이브러리 전체 해결로 확대하지 않는다. ## 빌드와 직접 검사 diff --git a/docs/tracking/decisions/0012-mitc4-drilling-options.md b/docs/tracking/decisions/0012-mitc4-drilling-options.md new file mode 100644 index 0000000..2a7254e --- /dev/null +++ b/docs/tracking/decisions/0012-mitc4-drilling-options.md @@ -0,0 +1,39 @@ +# MITC4 drilling의 직접 스프링 기본값과 기존 방식 선택 + +## 상황 + +기존 비선형 상대 비틀림은 drilling의 polar·minimal swing·atan2와 고차 미분을 계산한다. +사용자는 이 계산을 생략하는 작은 직접 스프링을 기본값으로 요청했고, 초기 director와 +전체 회전 기준에서 유한 강체 회전의 인공 에너지·복원모멘트가 생기는 것을 허용했다. +기존 polar1·관성3차 근사와 정확도 실패는 별도 보존 대상이다. + +## 결정 + +모델 전체의 compile-time 설정 하나로 직접 초기 director 스프링 또는 기존 상대 비틀림을 +선택하며 기본값은 직접 방식이다. `gamma=d0·phi`, `U=0.5*kd*sum(gamma²)`를 사용하고 +초기 physical rotary 대각8개 평균의 `1e-3` scale을 고정한다. +선형·비선형 정적·pressure·frequency·HHT에 같이 적용한다. +직접 HHT hidden 미분은 실제 cubic 관성이 계산한 null과 같은 scalar에서 도출한다. + +인공 내부 spatial 모멘트를 별도로 기록한다. 기존 raw 평형을 보존하고 직접 방식의 기존 +moment gate는 raw에서 인공 모멘트를 차감한다. 기존 방식은 raw 판정·현재 계산을 유지한다. +기존 feature/procedure/schema와 외부 U/UR 기준을 바꾸지 않고 추가 marker와 진단으로 구분한다. + +## 검토한 대안 + +기존 상대 비틀림만 유지하면 drilling 회전 계산 생략 요구를 충족하지 못한다. +직접 스프링에 상대 회전 보정을 다시 넣으면 허용한 비객관성을 제거하는 대신 그 계산을 +다시 도입하게 된다. 현재는 모델·section별 runtime 입력을 늘리는 대신 한 코드 설정을 택했다. +Global URZ에만 스프링을 놓으면 기울어진 초기 director의 drilling을 나타내지 못하므로 쓰지 않는다. + +## 결과와 한계 + +직접 force의 병진·절점 간 coupling과2차 방향미분은0이지만 HHT null의 고차 미분은 필요하다. +인공 에너지를 physical stress/resultant/energy에 혼합하거나 spring 질량을 추가하지 않는다. +Raw 평형만 보고 직접 방식의 유한 강체 회전을 불평형으로 거부해서는 안 되며, 보정 분모에 +인공 모멘트를 더해 실제 불평형을 숨겨서도 안 된다. 옛 markerless 파일의 방식은 미상이다. +현재 선형/eigen block이 두 방식에서 같은 것은 정상이다. + +이 선택은 기존 근사 실험의 D3·Legacy constraint/C08·mass-null/triad 실패나 비선형 시간정확도를 +해결하는 결정이 아니다. 정식 Direct HHT U 통과와 전체 Debug의10 FAIL/1 SKIP은 함께 남는다. +공식 성능 개선이나 전체 solver correctness의 인증도 포함하지 않는다. diff --git a/docs/tracking/decisions/index.md b/docs/tracking/decisions/index.md index 6ac7aed..8ec304d 100644 --- a/docs/tracking/decisions/index.md +++ b/docs/tracking/decisions/index.md @@ -13,3 +13,4 @@ | `0009-hex8-numerics.md` | Mean-volume 강성·reference row-sum mass·high+low 상태를 사용한다. | | `0010-thermal-discretization.md` | 다섯 두께 온도점·consistent capacity·Backward Euler를 사용한다. | | `0011-tolerance-experiment-acceptance.md` | 이번 MITC4 HHT 실험은 변위 검증 통과·실패10개 보류로 수용하며 실제 FAIL/SKIP과 후속 검증을 보존한다. | +| `0012-mitc4-drilling-options.md` | 초기 director 직접 스프링을 기본으로 하고 기존 상대 비틀림을 코드로 선택한다. 인공 spatial 모멘트와 raw/보정 평형을 구별한다. | diff --git a/docs/tracking/findings.md b/docs/tracking/findings.md index 9d1051f..e86ec1d 100644 --- a/docs/tracking/findings.md +++ b/docs/tracking/findings.md @@ -1,5 +1,28 @@ # 해결되지 않은 문제 +## 2026-10-20 현재 실행과 미완료 경계 + +직접 스프링 기본 source `6bbbb2b`의 전체 Native/CTest는 각각1035개 중 +1024 PASS·10 FAIL·1 INCONCLUSIVE SKIP이다. 기존 실패를 ID뿐 아니라 현재 raw 수치와 +발생 단계로 대조했다. 새 직접 스프링 수식·저장 연결 검사는 통과했으며 아래 근사 문제를 +이번 변경에서 자동 복원하거나 허용치를 더 완화하지 않았다. + +- 관성 D3 두 검사의 첫 차이는 약0.22761592214, 명시 Legacy drilling gradient 차이는 + 약0.00329135061318이다. Legacy C08a hidden ratio는0.00123371144037로1e-3을 초과한다. +- C08 controller4개는 첫 격자(h=1/64,1/16,1/256,1/64)의 inc1에서100회 correction을 + 소진했다. 시간수렴·차수 판정 이전 실패이며 이전 통과로 대신하지 않는다. +- LargeDirect는 accepted inc1/t=.25/iteration5의 mass-null 실패, UniformTranslation의 + 유한회전 subrun은 inc2/t=.02의 triad/rotation mismatch가 유지됐다. +- N7 finest 비교 하나는 arithmetic/solve resolution의 INCONCLUSIVE SKIP이다. + +정식 직접 HHT의729nodes·101frames와9node U101 families는 이번에 실제 통과했지만 +비선형 q/v/physical a의 시간정확도와 위 독립 물리 실패를 해결하지는 않는다. +기존 방식의 대형 HHT 외부 U 재해석은 미실행이다. 별도 관련162검사 중7 FAIL은 기존5개와 +Direct 기본값 전용 검사를 잘못 포함한2개로 분리했다. 두 선택 오류도 실제 FAIL로 보존했다. +정식 fine/fine_02/pressure·새 thermal 자료는 없으며 `reference/temp`를 쓰지 않았다. +이 미완료는 승인된 근사 정확도 조사와 외부 자료 제공이 필요한 범위이므로, 이번 drilling +선택 구현을 전체 solver 검증 통과로 확대하지 않는다. + ## 이번 실험의 보류 결정 2026-10-19 사용자는 Abaqus U 비교를 이번 실험의 검증 통과로 수용하고 아래 실패10개를 @@ -74,14 +97,17 @@ - 다음 접근: 새 reference의 identity·시점·hash와 blocking 계약을 확정해 별도로 검증한다. 음수 온도 clamp·원본 보정·RMS 완화로 기존 기록을 통과로 바꾸지 않는다. -## Eigen partial-BC 정적 평형 충돌 +## 기존 drilling 방식의 partial-BC raw 정적 gate 한계 -- 조건·관측: MITC4 eigen의 S3 partial rotational BC 사례와 public static 평형 사이의 - 충돌이 기능 기록에 보류 상태로 남아 있다. -- 영향: 모드별 Value 5% 승인 예외가 이 정적 평형 문제나 Abaqus 알고리즘 동등성을 해결하지 않는다. -- 지금 해결하지 않는 이유: 수식·공개 정적 평형 계약의 별도 조사이며 문서만으로 고칠 수 없다. -- 다음 접근: 동일 source BC에서 static과 eigen의 null/finite 공간·힘/모멘트 의미를 독립 대조한다. - 부분 회전구속을 full clamp로 대체하지 않는다. +- 조건·관측: 보존한 S3 partial rotational BC의 public static 사례는 현재 기본 직접 방식에서 + raw 인공 모멘트를 별도로 복구·차감한 뒤 기존1e-10 gate를 통과한다. 기존 방식의 raw + gate는 같은 사례에서 계속 거부되며 명시 Legacy 검사와 별도 build의 실제 진단에 보존했다. +- 영향: 직접 방식의 해당 충돌은 해결했지만 기존 방식의 raw 수용 범위를 넓힌 것은 아니다. + Eigen Value5% 예외나 모든 partial-BC 비선형 모델의 성공을 뜻하지 않는다. +- 지금 남겨 두는 이유: 기존 방식의 수치 판정 보존이 승인된 계약이다. 그 gate를 자동으로 + 보정값으로 바꾸거나 부분 회전구속을 full clamp로 대체하지 않는다. +- 다음 접근: 기존 방식의 수용 범위를 변경하려면 같은 BC에서 null/finite 공간과 + generalized/spatial 힘·모멘트의 계약을 별도로 판정한다. ## 성능 목표 미달과 과거 회귀 diff --git a/docs/tracking/status.md b/docs/tracking/status.md index 3a4ae86..a5677da 100644 --- a/docs/tracking/status.md +++ b/docs/tracking/status.md @@ -1,5 +1,38 @@ # 구현과 검증 상태 +## 2026-10-20 MITC4 drilling 선택 + +초기 director 직접 스프링을 기본으로 하고 기존 상대 비틀림을 코드 한 곳에서 선택하는 +구현을 선형·비선형 정적·pressure·frequency·HHT와 저장·reader·비교기에 연결했다. +직접 방식의 비객관성은 허용하며 raw/인공/보정 spatial 평형을 구별한다. +아래 실행 근거의 제품 source는 `6bbbb2ba8c8078c7a637886201750753da0af5c0`이다. + +| 현재 실제 실행 | 결과와 captured exit | +| --- | --- | +| Native x64 Debug 전체 | unit915:904 PASS/10 FAIL/1 SKIP, exit1; integration74와 reference46은 모두 PASS, exit0/0 | +| CMake x64 Debug 전체 CTest | 같은1035 identity,1024 PASS/10 FAIL/1 SKIP, exit8 | +| 별도 Python 전체 |164 PASS, exit0; pytest9.1.1/NumPy2.4.6/h5py3.16.0의 격리 Python3.11 환경 | +| 기존 방식 별도 Debug 관련 검사 |162 unique=155 PASS/7 FAIL. 기존 물리5개와 Direct 기본값 전용 선택 오류2개를 분리하며 실제 FAIL은 보존 | +| 두 방식의 정식 linear/coarse-NL/eigen |실제 생성·기존 비교 모두 PASS. 현재 선형 변위와 eigen table은 두 방식에서 정확히 동일 | +| 직접 방식 정식 HHT |Release CLI0/비교0,729nodes·101frames, 외부9node U101 families PASS; 최대 정규화 RMS0.2282326248391787%, UR 경고200개 | +| 기존 방식 실제 HHT |별도 Debug 소형3frame CLI/readback PASS.729node 외부 U 재해석은 미실행 | + +두 Debug build와 별도 기존 방식 Debug/공식 Direct Release build는 exit0이며 실제 +MSVC19.51.36260.0 x64/v145, toolset directory14.51.36231과 probe false/OFF를 확인했다. +Compile source·등록/실행 identity·CL/link·staged/loaded DLL·원본 hash를 대조했다. +정식 원본34파일은 불변이고 실제10 HDF5/128평형 frame의 추가 marker·스키마·대수 관계를 검사했다. +실행 명령/exit·raw/XML·현재 실패 원인·선별/제외 목록과 hash는 +`.build/drilling-options-20261020/T4/`에 보존했다. 이 자료가 없는 checkout의 현재 실행 근거로 간주하지 않는다. + +Native/CTest의10 FAIL은 기존 D3·명시 Legacy constraint/C08·후속 mass-null/triad 원인과 +현재 수치/실패 단계를 대조했다. N7 하나는 해상도 INCONCLUSIVE SKIP이다. Legacy의 선택 오류2개는 +Direct mode literal 또는 production Direct 수용을 전제한 검사이며 기본 Direct 전체에서는 PASS했다. +반력·영운동·실제 Legacy marker는 별도 HHT 출력으로 확인했고, 오류 선택의 FAIL을 PASS로 바꾸지 않았다. +정식 fine/fine_02/pressure·새 thermal 자료 누락과 기존 NT11 수치 불일치도 남는다. +구현과 관련 실행은 확보했으나 전체 물리/시간정확도 검증 완료는 아니다. 성능 인증은 하지 않았다. + +## 2026-10-19 수렴 실험 기록 + 2026-10-19 후속 사용자 결정으로 이번 MITC4 HHT 실험은 **변위 검증 통과·실패10개 보류로 조건부 완료**한다. [수용 결정](decisions/0011-tolerance-experiment-acceptance.md)을 따르며 아래 전체 검사 raw 통계와 N7·자료 부재·미해결 물리 문제는 보존한다. 전체 물리/시간정확도 검증 통과로 확대하지 않는다. @@ -60,7 +93,7 @@ Fresh Refactor·통합 Native·통합 CMake 실행물에서 같은 scratch 입 2. 구성된 전체 실행 결과를 바탕으로 MITC4 local 미분/물리 실패와 실제 후보 후속 gate 실패를 분리한다. 후보 비선형 q/v/physical a 시간수렴·N7 및 후보 원본 비교를 확보한다. 전체 검사 실행 사실과 미해결 수치 검증을 구분하고, 자료가 필요한 비교는 확보 뒤 수행한다. -3. Eigen partial-BC 정적 평형 충돌과 thermal 외부 불일치의 영향을 각각 판정한다. +3. 기존 drilling 방식의 partial-BC raw 정적 gate 한계와 thermal 외부 불일치의 영향을 각각 판정한다. 기준 완화·질량 변경·원본 보정으로 해결하지 않는다. 4. 같은 물리·허용치·출력·정밀도를 보존한 상태에서 성능 회귀 원인을 조사하고 최종 source의 Release 시간과 수치 결과를 다시 측정한다. 검증되지 않은 가속을 완료로 기록하지 않는다. diff --git a/include/fesa/analysis/analysis_state.h b/include/fesa/analysis/analysis_state.h index 0222ae2..e1baff6 100644 --- a/include/fesa/analysis/analysis_state.h +++ b/include/fesa/analysis/analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/status.h" @@ -72,6 +73,8 @@ class AnalysisState { double PhysicalStrainEnergy() const noexcept; /// @brief Returns global force and moment equilibrium components. const std::array& Equilibrium() const noexcept; + /// @brief Returns actual MITC4 mode and artificial spatial moment evidence. + const std::optional& Mitc4Drilling() const noexcept; /// @brief Returns normalized shell verification metrics. const std::array& VerificationMetrics() const noexcept; @@ -96,6 +99,7 @@ class AnalysisState { double physical_strain_energy_{0.0}; std::array equilibrium_{}; std::array verification_metrics_{}; + std::optional mitc4_drilling_; }; } // namespace fesa diff --git a/include/fesa/analysis/dynamic_analysis_state.h b/include/fesa/analysis/dynamic_analysis_state.h index 96bca40..d96f83c 100644 --- a/include/fesa/analysis/dynamic_analysis_state.h +++ b/include/fesa/analysis/dynamic_analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/source_identity.h" @@ -54,6 +55,7 @@ struct DynamicAcceptedFrameState { double external_work = 0.0; double energy_defect = 0.0; std::array global_equilibrium{}; + std::optional mitc4_drilling; std::array linear_momentum{}; std::array angular_momentum{}; std::array linear_momentum_rate{}; diff --git a/include/fesa/analysis/eigen_analysis_state.h b/include/fesa/analysis/eigen_analysis_state.h index d5aecd0..cbe376c 100644 --- a/include/fesa/analysis/eigen_analysis_state.h +++ b/include/fesa/analysis/eigen_analysis_state.h @@ -8,6 +8,7 @@ #include #include "fesa/math/matrix.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/solvers/eigen/eigen_solver.h" namespace fesa { @@ -53,6 +54,7 @@ struct EigenAnalysisState { explicit EigenAnalysisState(FeastEigenResult backend_result) : backend{std::move(backend_result)} {} + std::optional drilling_mode; std::size_t full_dof_count{0U}; std::size_t free_dof_count{0U}; std::size_t mass_null_dimension{0U}; diff --git a/include/fesa/analysis/nonlinear_analysis_state.h b/include/fesa/analysis/nonlinear_analysis_state.h index 7efee97..c0e3127 100644 --- a/include/fesa/analysis/nonlinear_analysis_state.h +++ b/include/fesa/analysis/nonlinear_analysis_state.h @@ -3,6 +3,7 @@ #include #include +#include #include #include "fesa/core/source_identity.h" @@ -46,6 +47,7 @@ struct NonlinearTrialState { std::vector shell_results; double physical_strain_energy; std::array global_equilibrium; + std::optional mitc4_drilling; }; /// @brief Stores one authoritative nonlinear increment frame. @@ -69,6 +71,7 @@ struct AcceptedFrameState { std::vector shell_results; double physical_strain_energy; std::array global_equilibrium; + std::optional mitc4_drilling; }; /// @brief Owns committed nonlinear state and validated accepted frames. diff --git a/include/fesa/elements/element.h b/include/fesa/elements/element.h index 5430172..6a1ec5d 100644 --- a/include/fesa/elements/element.h +++ b/include/fesa/elements/element.h @@ -1,11 +1,14 @@ #ifndef FESA_ELEMENTS_ELEMENT_H_ #define FESA_ELEMENTS_ELEMENT_H_ +#include +#include #include #include #include "fesa/core/status.h" #include "fesa/elements/element_dof_layout.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/math/matrix.h" #include "fesa/math/vector.h" #include "fesa/results/result_records.h" @@ -30,6 +33,8 @@ struct BeamElementResultRows { struct ShellElementResultRows { std::vector rows; double physical_strain_energy{0.0}; + std::optional drilling_mode; + std::array drilling_moment{}; }; /// @brief Selects the physical recovery shape of one runtime element. diff --git a/include/fesa/elements/element_capabilities.h b/include/fesa/elements/element_capabilities.h index d6c610d..c68ed6a 100644 --- a/include/fesa/elements/element_capabilities.h +++ b/include/fesa/elements/element_capabilities.h @@ -178,6 +178,8 @@ class ShellReferenceElement : public Element { virtual Result Stiffness() const = 0; /// @copydoc Mitc4Shell::Mass virtual Result Mass(double density) const = 0; + /// @brief Returns the actual immutable reference-kernel drilling formulation. + virtual Mitc4DrillingMode DrillingMode() const noexcept = 0; }; /// @brief Borrows the same linear MITC4 kernel for stiffness and recovery. @@ -190,6 +192,7 @@ class Mitc4ReferenceAdapter final : public ShellReferenceElement { Result Recover(const Vector& u) const override; Result Stiffness() const override; Result Mass(double density) const override; + Mitc4DrillingMode DrillingMode() const noexcept override; private: const Mitc4Shell& kernel_; diff --git a/include/fesa/elements/mitc4_drilling_mode.h b/include/fesa/elements/mitc4_drilling_mode.h new file mode 100644 index 0000000..8df7201 --- /dev/null +++ b/include/fesa/elements/mitc4_drilling_mode.h @@ -0,0 +1,26 @@ +#ifndef FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ +#define FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ + +namespace fesa { + +/// @brief Identifies the immutable numerical MITC4 drilling formulation. +enum class Mitc4DrillingMode { kDirectSpring, kLegacyCorotatedTwist }; + +/// @brief Returns the stable storage name of a supported drilling formulation. +constexpr const char* Mitc4DrillingModeName(Mitc4DrillingMode mode) noexcept { + switch (mode) { + case Mitc4DrillingMode::kDirectSpring: + return "direct-initial-director-spring"; + case Mitc4DrillingMode::kLegacyCorotatedTwist: + return "legacy-corotated-twist"; + } + return ""; +} + +namespace mitc4_drilling_internal { +struct TestAccess; +} // namespace mitc4_drilling_internal + +} // namespace fesa + +#endif // FESA_ELEMENTS_MITC4_DRILLING_MODE_H_ diff --git a/include/fesa/elements/mitc4_shell.h b/include/fesa/elements/mitc4_shell.h index 36e954d..657008e 100644 --- a/include/fesa/elements/mitc4_shell.h +++ b/include/fesa/elements/mitc4_shell.h @@ -8,6 +8,7 @@ #include "fesa/core/status.h" #include "fesa/elements/element.h" #include "fesa/elements/mitc4_shell_definition.h" +#include "fesa/elements/mitc4_drilling_mode.h" #include "fesa/materials/isotropic_linear_elastic_material.h" #include "fesa/math/matrix.h" #include "fesa/math/vector.h" @@ -92,6 +93,11 @@ class Mitc4Shell final : public Element { std::array, 4> initial_directors, const ShellSection& section, const LinearElasticMaterial& material); + /// @brief Returns this immutable kernel's production drilling formulation. + [[nodiscard]] Mitc4DrillingMode DrillingMode() const noexcept { + return drilling_mode_; + } + /// @brief Returns the factory-bound MITC4 DOF order. const ElementDofLayout& DofLayout() const noexcept override; @@ -163,6 +169,7 @@ class Mitc4Shell final : public Element { private: friend class ElementFactory; + friend struct mitc4_drilling_internal::TestAccess; /// @brief Binds Domain identity after the numerical candidate is valid. void BindRuntime(ElementDofLayout layout, EntityIndex element_index, @@ -203,6 +210,7 @@ class Mitc4Shell final : public Element { double thickness_; double youngs_modulus_; double poisson_ratio_; + Mitc4DrillingMode drilling_mode_; SourceLocation source_location_; std::string identity_; ElementDofLayout dof_layout_; diff --git a/include/fesa/elements/nonlinear_mitc4_shell.h b/include/fesa/elements/nonlinear_mitc4_shell.h index 4359400..e880fdc 100644 --- a/include/fesa/elements/nonlinear_mitc4_shell.h +++ b/include/fesa/elements/nonlinear_mitc4_shell.h @@ -54,6 +54,8 @@ struct NonlinearMitc4InternalEquilibrium { Matrix drilling_metric{0U, 0U}; double drilling_energy = 0.0; Vector force_directional_second{0U}; + /// @brief Stores D2(f_drill)[v,v] separately; direct springs give exact zero. + Vector drilling_force_directional_second{0U}; Vector drilling_constraint{0U}; Matrix drilling_constraint_gradient{0U, 0U}; std::array drilling_constraint_hessian{ @@ -104,6 +106,11 @@ class NonlinearMitc4Shell { std::array initial_frames, const ShellSection& section, const LinearElasticMaterial& material); + /// @brief Returns this immutable kernel's numerical drilling formulation. + [[nodiscard]] Mitc4DrillingMode DrillingMode() const noexcept { + return drilling_mode_; + } + /// @brief Evaluates trial physical kinematics and fixed recovery rows. /// @param trial Global nodal translation and fixed-global rotation vectors. /// @return Physical recovery or a fail-closed trial diagnostic. @@ -134,6 +141,7 @@ class NonlinearMitc4Shell { private: friend struct nonlinear_mitc4_internal::ConfigurationSnapshot; + friend struct mitc4_drilling_internal::TestAccess; /// @brief Retains validated fixed-reference data without analysis ownership. NonlinearMitc4Shell(std::array coordinates, @@ -155,6 +163,7 @@ class NonlinearMitc4Shell { double youngs_modulus_; double poisson_ratio_; double drilling_stiffness_; + Mitc4DrillingMode drilling_mode_; SourceLocation source_location_; std::string identity_; }; diff --git a/include/fesa/results/result_records.h b/include/fesa/results/result_records.h index c1a5848..b867cab 100644 --- a/include/fesa/results/result_records.h +++ b/include/fesa/results/result_records.h @@ -4,9 +4,11 @@ #include #include #include +#include #include #include "fesa/core/source_identity.h" +#include "fesa/elements/mitc4_drilling_mode.h" namespace fesa { @@ -70,6 +72,15 @@ struct ShellResultRow { std::array stress; }; +/// @brief Separates artificial internal spatial moment from raw equilibrium. +/// @note The actual immutable kernel mode accompanies every production frame. +struct Mitc4DrillingEvidence { + Mitc4DrillingMode mode; + std::array drilling_moment{}; + // Raw force balance and raw moment balance minus positive internal drilling. + std::array equilibrium_with_drilling{}; +}; + /// @brief Carries a complete shell result candidate for atomic validation. struct ShellStateCandidate { std::vector rows; @@ -79,6 +90,7 @@ struct ShellStateCandidate { // [FREE_RESIDUAL_NORMALIZED,FORCE_BALANCE_NORMALIZED, // MOMENT_BALANCE_NORMALIZED]. std::array verification_metrics{}; + std::optional mitc4_drilling; }; } // namespace fesa diff --git a/msvc/Fesa.Common.props b/msvc/Fesa.Common.props index 34bf8c1..f8a2f99 100644 --- a/msvc/Fesa.Common.props +++ b/msvc/Fesa.Common.props @@ -45,7 +45,7 @@ WIN32;_WINDOWS;%(PreprocessorDefinitions) NDEBUG;WIN32;_WINDOWS;%(PreprocessorDefinitions) $(FesaRoot)include;%(AdditionalIncludeDirectories) - $(FesaRoot)include;$(FesaRoot)src\fesa;%(AdditionalIncludeDirectories) + $(FesaRoot)include;$(FesaRoot)src\fesa;%(AdditionalIncludeDirectories) $(FesaExternalOptions) /MP2 %(AdditionalOptions) TurnOffAllWarnings diff --git a/scripts/verification/compare_hht_dynamic.py b/scripts/verification/compare_hht_dynamic.py index 074830a..98114f0 100644 --- a/scripts/verification/compare_hht_dynamic.py +++ b/scripts/verification/compare_hht_dynamic.py @@ -197,6 +197,38 @@ def dataset(file: h5py.File, path: str, shape: tuple[int, ...]) -> h5py.Dataset: return result +def drilling_diagnostics(file: h5py.File) -> dict[str, Any] | None: + """Validate additive diagnostics without introducing a balance verdict.""" + metadata = file["metadata"] + mode_present = "mitc4_drilling_formulation" in metadata.attrs + version_present = "drilling_diagnostics_version" in metadata.attrs + require(mode_present == version_present, "partial-drilling-markers", "Both drilling markers are required.") + if not mode_present: + return None + mode = string_attribute(metadata, "mitc4_drilling_formulation") + require(mode in ("direct-initial-director-spring", "legacy-corotated-twist"), + "unknown-drilling-formulation", mode) + numeric_attribute(metadata, "drilling_diagnostics_version", " tuple[list[np.ndarray], dict[str, Any]]: """Precheck every result node/frame/component, then project by source key.""" @@ -319,9 +351,10 @@ def read_results(path: Path, nodes: dict[int, tuple[float, ...]], if frame == 0: require((values == 0.0).all(), "nonzero-initial-displacement", "Frame zero must have q0=0.") projected.append(values[selection]) + drilling = drilling_diagnostics(file) return projected, {"source_input_identity": source_identity, "source_node_count": NODE_COUNT, "accepted_frame_count": FRAME_COUNT, "checked_result_scalars": NODE_COUNT * FRAME_COUNT * 6, - "projected_node_labels": list(NODE_LABELS)} + "projected_node_labels": list(NODE_LABELS), "drilling_diagnostics": drilling} def evaluate_family(increment: int, quantity: str, result: list[float], @@ -374,6 +407,7 @@ def evaluate_family(increment: int, quantity: str, result: list[float], def compare(input_path: Path, reference_path: Path, results_path: Path) -> dict[str, Any]: report: dict[str, Any] = { + "drilling_diagnostics": None, "schema": "hht-reference-comparison-v1", "case": CASE, "policy": dict(POLICY), "identities": {}, "precheck": {"passed": False}, "families": [], "blocking_failures": [], "warnings": [], "blocking_passed": False, "verdict": "failed", @@ -386,6 +420,7 @@ def compare(input_path: Path, reference_path: Path, results_path: Path) -> dict[ nodes, content_identity = input_inventory(input_path) reference = read_reference(reference_path) result, evidence = read_results(results_path, nodes, content_identity) + report["drilling_diagnostics"] = evidence.pop("drilling_diagnostics") families = [] for frame in range(FRAME_COUNT): for quantity, first in (("U", 0), ("UR", 3)): diff --git a/src/fesa/analysis/analysis_state.cpp b/src/fesa/analysis/analysis_state.cpp index 786be82..dabb8a5 100644 --- a/src/fesa/analysis/analysis_state.cpp +++ b/src/fesa/analysis/analysis_state.cpp @@ -1,5 +1,7 @@ #include "fesa/analysis/analysis_state.h" +#include "results/mitc4_drilling_evidence.h" + #include #include #include @@ -167,9 +169,18 @@ Status AnalysisState::CommitShellResults( "Shell energy, equilibrium, and normalized metrics must be finite."); } + if (!results_internal::ValidDrillingEvidence(candidate.mitc4_drilling, + candidate.equilibrium) || + (mitc4_drilling_ && !results_internal::SameDrillingInventory( + mitc4_drilling_, candidate.mitc4_drilling))) + return ShellCandidateFailure("invalid-shell-drilling-evidence", + identity_.step_name, + "Drilling evidence must be finite, consistent, and retain its mode."); + shell_results_ = std::move(candidate.rows); physical_strain_energy_ = candidate.physical_strain_energy; equilibrium_ = candidate.equilibrium; + mitc4_drilling_ = candidate.mitc4_drilling; verification_metrics_ = candidate.verification_metrics; return Status::Ok(); } @@ -187,6 +198,11 @@ const std::array& AnalysisState::Equilibrium() const noexcept { return equilibrium_; } +const std::optional& AnalysisState::Mitc4Drilling() + const noexcept { + return mitc4_drilling_; +} + const std::array& AnalysisState::VerificationMetrics() const noexcept { return verification_metrics_; diff --git a/src/fesa/analysis/dynamic_analysis_state.cpp b/src/fesa/analysis/dynamic_analysis_state.cpp index 3e374c0..06bb1f7 100644 --- a/src/fesa/analysis/dynamic_analysis_state.cpp +++ b/src/fesa/analysis/dynamic_analysis_state.cpp @@ -1,5 +1,7 @@ #include "fesa/analysis/dynamic_analysis_state.h" +#include "results/mitc4_drilling_evidence.h" + #include #include #include @@ -273,6 +275,12 @@ Status DynamicAnalysisState::CommitAcceptedFrame( return Failure(FailureCategory::kSolver, "Frame physical evidence is nonfinite or has invalid size."); } + if (!results_internal::ValidDrillingEvidence(frame.mitc4_drilling, + frame.global_equilibrium) || + (!accepted_frames_.empty() && !results_internal::SameDrillingInventory( + frame.mitc4_drilling, accepted_frames_.back().mitc4_drilling))) + return Failure(FailureCategory::kSolver, + "Accepted drilling evidence must be finite, consistent, and retain its mode."); for (const std::size_t dof : constrained_dofs_) { if (frame.generalized_coordinates[dof] != 0.0 || frame.velocity[dof] != 0.0 || frame.acceleration[dof] != 0.0 || diff --git a/src/fesa/analysis/hht_dynamic_algebra.cpp b/src/fesa/analysis/hht_dynamic_algebra.cpp index 09d6803..9311989 100644 --- a/src/fesa/analysis/hht_dynamic_algebra.cpp +++ b/src/fesa/analysis/hht_dynamic_algebra.cpp @@ -547,6 +547,18 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const SparseMatrix& gradient, DynamicAcceptedFrameState* frame, const bool require_acceleration) { + const bool direct = frame->mitc4_drilling && + frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring; + if (direct) { + for (const std::size_t node : nodes) { + if (node >= assembly.hidden_constraint.Size()) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must cover every admissible null node."); + if (!std::isfinite(assembly.hidden_constraint[node])) + return Failure(domain, "dynamic-invalid-algebra-state", + "Canonical hidden values must remain finite."); + } + } Vector psi{nodes.size()}; Vector psi_scale{nodes.size()}; Vector velocity_defect{nodes.size()}; @@ -575,9 +587,10 @@ Status HiddenEvidence(const Domain& domain, const DynamicAssembly& assembly, const double g = frame->internal_force[full] - frame->external_force[full]; const double n = rotation.drill_null[c]; - psi[row] += n * g; + if (!direct) psi[row] += n * g; psi_scale[row] += std::abs(n) * frame->residual_scale[full]; } + if (direct) psi[row] = assembly.hidden_constraint[node]; velocity_defect[row] += velocity_correction; AddCompensated(assembly.hidden_second[node], &acceleration_defect[row], &acceleration_correction); diff --git a/src/fesa/analysis/hht_dynamic_assembly.cpp b/src/fesa/analysis/hht_dynamic_assembly.cpp index f4a5e1a..8faf890 100644 --- a/src/fesa/analysis/hht_dynamic_assembly.cpp +++ b/src/fesa/analysis/hht_dynamic_assembly.cpp @@ -18,6 +18,7 @@ #include "fesa/elements/finite_rotation.h" #include "fesa/loads/load.h" #include "fesa/model/source_target_resolver.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa::hht_internal { @@ -498,6 +499,10 @@ Result EvaluateAssembly( internal[e] = std::move(f.Value()); inertia[e] = std::move(m.Value()); }); + if (shells.empty()) + return Result::Failure(Failure("HHT requires a shell kernel.")); + frame.mitc4_drilling = + Mitc4DrillingEvidence{shells.front().internal_kernel->DrillingMode()}; frame.inertia_force = Vector(count); frame.internal_force = Vector(count); frame.residual = Vector(count); @@ -515,6 +520,7 @@ Result EvaluateAssembly( std::vector> constraint_gradients(domain.Nodes().size()); std::vector external_done(domain.Nodes().size(), false); + Vector hidden_constraint(directional_third ? domain.Nodes().size() : 0U); Vector third(count); Vector third_scale(count); std::array, 6> entries; @@ -532,9 +538,31 @@ Result EvaluateAssembly( return Result::Failure(*failures[e]); const auto& f = *internal[e]; const auto& m = *inertia[e]; + if (shells[e].internal_kernel->DrillingMode() != frame.mitc4_drilling->mode || + f.drilling_force.Size() != 24U) + return Result::Failure(Failure("Invalid actual drilling mode or force.")); + for (std::size_t node = 0U; node < 4U; ++node) { + const auto local = 6U * node; + const auto global = shells[e].scatter[local]; + const Vector3 force = Triple(f.drilling_force, local); + const auto spatial = SpatialMomentSO3( + Triple(frame.generalized_coordinates, global + 3U), + Triple(f.drilling_force, local + 3U)); + if (!spatial.HasValue()) + return Result::Failure(spatial.GetStatus()); + const Vector3 position = Vector3(domain.Nodes()[global / 6U].coordinates) + + Triple(frame.generalized_coordinates, global); + const Vector3 moment = position.Cross(force) + spatial.Value(); + for (std::size_t c = 0U; c < 3U; ++c) + frame.mitc4_drilling->drilling_moment[c] += moment[c]; + } + for (const double value : frame.mitc4_drilling->drilling_moment) + if (!std::isfinite(value)) + return Result::Failure(Failure("Drilling spatial moment overflowed.")); if (directional_third) { for (std::size_t local = 0U; local < 4U; ++local) { const auto node = shells[e].scatter[6U * local] / 6U; + hidden_constraint[node] += f.drilling_constraint[local]; auto& gradient = constraint_gradients[node]; auto& hessian = constraint_hessians[node]; for (std::size_t i = 0U; i < 24U; ++i) @@ -555,6 +583,7 @@ Result EvaluateAssembly( Triple(frame.external_force, first), {}, {}); if (!external.HasValue()) return Result::Failure(external.GetStatus()); + hidden_constraint[node] += external.Value().value; for (std::size_t i = 0U; i < 3U; ++i) { gradient[first + i] += external.Value().gradient[i]; for (std::size_t j = 0U; j < 3U; ++j) @@ -740,7 +769,8 @@ Result EvaluateAssembly( std::move(third), std::move(third_scale), std::move(rotations), momentum_scale, angular_scale, std::move(hidden_second), std::move(hidden_second_scale), std::move(hidden_gradient), - std::move(trial_constraints), std::move(field_flux)}); + std::move(trial_constraints), std::move(field_flux), + std::move(hidden_constraint)}); } Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs, @@ -776,10 +806,20 @@ Status CompleteEquilibrium(const Domain& domain, const DofManager& dofs, frame->global_equilibrium[c] = force_error[c]; frame->global_equilibrium[3U + c] = moment_error[c]; } + Vector3 acceptance_moment = moment_error; + if (frame->mitc4_drilling) { + results_internal::CompleteDrillingEvidence(frame->global_equilibrium, + &*frame->mitc4_drilling); + if (!results_internal::ValidDrillingEvidence(frame->mitc4_drilling, + frame->global_equilibrium)) + return Failure("Invalid dynamic drilling equilibrium evidence."); + if (frame->mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) + acceptance_moment = moment_error - Vector3(frame->mitc4_drilling->drilling_moment); + } if (!force_error.IsFinite() || !moment_error.IsFinite() || !std::isfinite(force_scale) || !std::isfinite(moment_scale) || Ratio(force_error.Norm(), force_scale) > 1e-8 || - Ratio(moment_error.Norm(), moment_scale) > 1e-8) + Ratio(acceptance_moment.Norm(), moment_scale) > 1e-8) return Failure("Dynamic physical force or moment balance exceeds 1e-8."); return Status::Ok(); } diff --git a/src/fesa/analysis/hht_dynamic_internal.h b/src/fesa/analysis/hht_dynamic_internal.h index e18dd62..5c5642c 100644 --- a/src/fesa/analysis/hht_dynamic_internal.h +++ b/src/fesa/analysis/hht_dynamic_internal.h @@ -97,6 +97,9 @@ struct DynamicAssembly { std::optional hidden_gradient; std::optional trial_constraints; std::array, 2> field_flux; + // Source-ordered scalar primitive, including the external term exactly once. + // Empty when directional primitives were not evaluated; never a zero fallback. + Vector hidden_constraint{0U}; }; /// @brief Builds checked shell evaluators in stable source order. diff --git a/src/fesa/analysis/hht_reduced_operators.cpp b/src/fesa/analysis/hht_reduced_operators.cpp index be4781f..e94d382 100644 --- a/src/fesa/analysis/hht_reduced_operators.cpp +++ b/src/fesa/analysis/hht_reduced_operators.cpp @@ -82,7 +82,8 @@ const hht_internal::DirectionalShellConstraintData* MatchingConstraints( nonlinear_mitc4_internal::ConfigurationSnapshot::Read(*shell.internal_kernel); if (record.configuration.source.file != configuration.source.file || record.configuration.source.line != configuration.source.line || - record.configuration.identity != configuration.identity) + record.configuration.identity != configuration.identity || + record.configuration.drilling_mode != configuration.drilling_mode) return nullptr; for (std::size_t i = 0U; i < configuration.values.size(); ++i) if (!std::isfinite(record.configuration.values[i]) || diff --git a/src/fesa/analysis/linear_frequency_procedure.cpp b/src/fesa/analysis/linear_frequency_procedure.cpp index d351098..3c9a0c9 100644 --- a/src/fesa/analysis/linear_frequency_procedure.cpp +++ b/src/fesa/analysis/linear_frequency_procedure.cpp @@ -217,6 +217,13 @@ Result ExecuteCandidate(const Domain& domain, if (!backend.HasValue()) return Result::Failure(backend.GetStatus()); EigenAnalysisState state(std::move(backend.Value())); + for (const auto& reference : references) { + if (state.drilling_mode && *state.drilling_mode != reference->DrillingMode()) + return Result::Failure(Failure( + "invalid-shell-drilling-evidence", "mode", + "All frequency kernels must retain one actual drilling mode.")); + state.drilling_mode = reference->DrillingMode(); + } state.full_dof_count = dofs.FullDofCount(); state.free_dof_count = reduction.Value().free_count; state.mass_null_dimension = reduction.Value().null_count; diff --git a/src/fesa/analysis/nonlinear_analysis_state.cpp b/src/fesa/analysis/nonlinear_analysis_state.cpp index 23c3ee4..d96f3f0 100644 --- a/src/fesa/analysis/nonlinear_analysis_state.cpp +++ b/src/fesa/analysis/nonlinear_analysis_state.cpp @@ -12,6 +12,7 @@ #include "fesa/elements/element_definition.h" #include "math/finite_values.h" #include "results/shell_result_recovery.h" +#include "results/mitc4_drilling_evidence.h" #include "analysis/shell_state_inventory.h" namespace fesa { @@ -110,6 +111,13 @@ Status ValidateAcceptedFrame(const AcceptedFrameState& frame, "Accepted frame physical result inventory is invalid."); } + if (!results_internal::ValidDrillingEvidence(frame.mitc4_drilling, + frame.global_equilibrium) || + (accepted_frame_count != 0U && !results_internal::SameDrillingInventory( + frame.mitc4_drilling, committed.mitc4_drilling))) + return StateFailure(FailureCategory::kSolver, "invalid-shell-drilling-evidence", {}, + "Accepted drilling evidence must be finite and retain its mode."); + const double gauss = 1.0 / std::sqrt(3.0); const std::array expected_locations{ShellMidsurfaceLocation::kGp1, @@ -243,7 +251,7 @@ Status NonlinearAnalysisState::CommitAcceptedFrame( frame.reaction, frame.shell_results, frame.physical_strain_energy, - frame.global_equilibrium}; + frame.global_equilibrium, frame.mitc4_drilling}; std::vector accepted_candidate = accepted_frames_; accepted_candidate.push_back(frame); diff --git a/src/fesa/analysis/nonlinear_static_procedure.cpp b/src/fesa/analysis/nonlinear_static_procedure.cpp index 4412f22..a47760e 100644 --- a/src/fesa/analysis/nonlinear_static_procedure.cpp +++ b/src/fesa/analysis/nonlinear_static_procedure.cpp @@ -33,6 +33,7 @@ #include "assembly/ordered_contribution_internal.h" #include "math/finite_values.h" #include "elements/shell_kernel_values.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa { namespace { @@ -57,6 +58,7 @@ struct InternalAssembly { SparseMatrix tangent; std::vector shell_rows; double physical_strain_energy; + Mitc4DrillingEvidence mitc4_drilling; }; Status ProcedureFailure(const FailureCategory category, const std::string& code, @@ -175,6 +177,10 @@ Result BuildBaseline(const Domain& domain) { Result> BuildRuntimeShells(const Domain& domain, const DofManager& dofs) { std::vector shells; + if (domain.Elements().Empty()) + return Result>::Failure(ProcedureFailure( + FailureCategory::kModel, "unsupported-nonlinear-element", {}, "ELEMENT", + {}, "The nonlinear static procedure requires a MITC4 shell element.")); shells.reserve(domain.Elements().Size()); for (std::size_t index = 0U; index < domain.Elements().Size(); ++index) { const auto* definition = @@ -320,7 +326,7 @@ NonlinearMitc4TrialDofs ExtractElementTrial( Result AssembleInternal( const std::vector& shells, const DofManager& dofs, - const Vector& coordinates, const ParallelFor& parallel_for) { + const Domain& domain, const Vector& coordinates, const ParallelFor& parallel_for) { std::vector> local( shells.size()); std::vector> failures(shells.size()); @@ -341,6 +347,7 @@ Result AssembleInternal( std::vector rows; rows.reserve(shells.size() * 4U); double energy = 0.0; + Mitc4DrillingEvidence drilling{shells.front().kernel->DrillingMode()}; for (std::size_t index = 0U; index < shells.size(); ++index) { if (failures[index].has_value()) { return Result::Failure(*failures[index]); @@ -354,6 +361,28 @@ Result AssembleInternal( "NONLINEAR_STATIC_PROCEDURE", std::to_string(index), "A nonlinear shell internal residual or tangent is nonfinite.")); } + if (shells[index].kernel->DrillingMode() != drilling.mode || + local[index]->drilling_force.Size() != kShellDofCount || + !IsFinite(local[index]->drilling_force)) + return Result::Failure(ProcedureFailure( + FailureCategory::kSolver, "invalid-shell-drilling-evidence", {}, + "NONLINEAR_STATIC_PROCEDURE", {}, "Drilling mode or force is invalid.")); + const auto element_trial = ExtractElementTrial(shells[index], coordinates); + for (std::size_t node = 0U; node < 4U; ++node) { + const auto first = 6U * node; + const auto global_node = shells[index].scatter[first] / 6U; + const auto& values = local[index]->drilling_force; + const Vector3 drilling_force{values[first], values[first + 1U], values[first + 2U]}; + const auto spatial = SpatialMomentSO3(element_trial.rotation_vector[node], + {values[first + 3U], values[first + 4U], values[first + 5U]}); + if (!spatial.HasValue()) + return Result::Failure(spatial.GetStatus()); + const Vector3 position = Vector3(domain.Nodes()[global_node].coordinates) + + element_trial.displacement[node]; + const Vector3 moment = position.Cross(drilling_force) + spatial.Value(); + for (std::size_t c = 0U; c < 3U; ++c) + drilling.drilling_moment[c] += moment[c]; + } for (std::size_t row = 0U; row < kShellDofCount; ++row) { force[shells[index].scatter[row]] += local[index]->internal_residual[row]; assembly_internal::AppendLocalMatrixRow( @@ -373,7 +402,8 @@ Result AssembleInternal( } energy += local[index]->physical_strain_energy; } - if (!IsFinite(force) || !std::isfinite(energy)) { + if (!IsFinite(force) || !std::isfinite(energy) || + !IsFinite(drilling.drilling_moment)) { return Result::Failure(ProcedureFailure( FailureCategory::kSolver, "nonfinite-newton-data", {}, "NONLINEAR_STATIC_PROCEDURE", {}, @@ -387,7 +417,7 @@ Result AssembleInternal( return Result::Failure(tangent.GetStatus()); } return Result::Success( - {std::move(force), std::move(tangent.Value()), std::move(rows), energy}); + {std::move(force), std::move(tangent.Value()), std::move(rows), energy, drilling}); } Result Difference(const Vector& left, const Vector& right) { @@ -699,7 +729,7 @@ AcceptedFrameState AcceptedFrame(const std::uint32_t frame_index, trial.reaction, trial.shell_results, trial.physical_strain_energy, - trial.global_equilibrium}; + trial.global_equilibrium, trial.mitc4_drilling}; } } // namespace @@ -807,7 +837,7 @@ Result ExecuteNonlinearStaticCore( NewtonIterationStage::kKinematicsEvaluated; Result internal = AssembleInternal( - shells.Value(), dofs.Value(), + shells.Value(), dofs.Value(), domain, trial.Value().generalized_coordinates, parallel_for); if (!internal.HasValue()) { return Result>::Failure( @@ -879,6 +909,9 @@ Result ExecuteNonlinearStaticCore( NewtonIterationStage::kExternalAssembled)); } trial.Value().global_equilibrium = global_equilibrium.Value(); + trial.Value().mitc4_drilling = internal.Value().mitc4_drilling; + results_internal::CompleteDrillingEvidence( + trial.Value().global_equilibrium, &*trial.Value().mitc4_drilling); Result partition = constraints.Partition(tangent.Value(), dofs.Value()); diff --git a/src/fesa/elements/element_capabilities.cpp b/src/fesa/elements/element_capabilities.cpp index 772fc0c..f780864 100644 --- a/src/fesa/elements/element_capabilities.cpp +++ b/src/fesa/elements/element_capabilities.cpp @@ -95,6 +95,10 @@ Result Ds4ThermalAdapter::Recover( Mitc4ReferenceAdapter::Mitc4ReferenceAdapter(const Mitc4Shell& kernel) : kernel_{kernel} {} +Mitc4DrillingMode Mitc4ReferenceAdapter::DrillingMode() const noexcept { + return kernel_.DrillingMode(); +} + const ElementDofLayout& Mitc4ReferenceAdapter::DofLayout() const noexcept { return kernel_.DofLayout(); } diff --git a/src/fesa/elements/mitc4_drilling_internal.h b/src/fesa/elements/mitc4_drilling_internal.h new file mode 100644 index 0000000..820e303 --- /dev/null +++ b/src/fesa/elements/mitc4_drilling_internal.h @@ -0,0 +1,28 @@ +#ifndef FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ +#define FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ + +#include "fesa/elements/nonlinear_mitc4_shell.h" + +namespace fesa::mitc4_drilling_internal { + +// The sole production selection. Change this constant and rebuild all targets. +inline constexpr Mitc4DrillingMode kProductionDrillingMode = + Mitc4DrillingMode::kDirectSpring; + +/// @brief Provides explicit alternate formulations only for internal tests. +/// @note Production factories use the single compile-time selection above. +struct TestAccess { + static Mitc4Shell WithMode(Mitc4Shell shell, Mitc4DrillingMode mode) { + shell.drilling_mode_ = mode; + return shell; + } + static NonlinearMitc4Shell WithMode(NonlinearMitc4Shell shell, + Mitc4DrillingMode mode) { + shell.drilling_mode_ = mode; + return shell; + } +}; + +} // namespace fesa::mitc4_drilling_internal + +#endif // FESA_SRC_FESA_ELEMENTS_MITC4_DRILLING_INTERNAL_H_ diff --git a/src/fesa/elements/mitc4_shell.cpp b/src/fesa/elements/mitc4_shell.cpp index 6aea36e..ef7f002 100644 --- a/src/fesa/elements/mitc4_shell.cpp +++ b/src/fesa/elements/mitc4_shell.cpp @@ -12,6 +12,7 @@ #include "math/finite_values.h" #include "math/vector3_operations.h" #include "elements/shell_kernel_values.h" +#include "elements/mitc4_drilling_internal.h" namespace fesa { @@ -957,6 +958,22 @@ Result Mitc4Shell::Recover( ShellElementResultRows rows{}; rows.rows.reserve(kNodeCount); rows.physical_strain_energy = recovered.Value().strain_energy; + const auto stiffness = Stiffness(); + if (!stiffness.HasValue()) + return Result::Failure(stiffness.GetStatus()); + const Vector drilling_force = + stiffness.Value().drilling_global24.Multiply(element_displacement); + rows.drilling_mode = DrillingMode(); + Vector3 drilling_moment; + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const std::size_t first = 6U * node; + const Vector3 force{drilling_force[first], drilling_force[first + 1U], + drilling_force[first + 2U]}; + const Vector3 moment{drilling_force[first + 3U], drilling_force[first + 4U], + drilling_force[first + 5U]}; + drilling_moment = drilling_moment + coordinates_[node].Cross(force) + moment; + } + rows.drilling_moment = drilling_moment.Components(); for (std::size_t point = 0U; point < recovered.Value().points.size(); ++point) { const auto& physical_point = recovered.Value().points[point]; @@ -1002,6 +1019,7 @@ Mitc4Shell::Mitc4Shell(std::array coordinates, thickness_{thickness}, youngs_modulus_{youngs_modulus}, poisson_ratio_{poisson_ratio}, + drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode}, source_location_{std::move(source_location)}, identity_{std::move(identity)} {} diff --git a/src/fesa/elements/nonlinear_mitc4_configuration_internal.h b/src/fesa/elements/nonlinear_mitc4_configuration_internal.h index 618eebc..304c50a 100644 --- a/src/fesa/elements/nonlinear_mitc4_configuration_internal.h +++ b/src/fesa/elements/nonlinear_mitc4_configuration_internal.h @@ -12,6 +12,7 @@ namespace fesa::nonlinear_mitc4_internal { /// key. struct ConfigurationSnapshot { std::array values{}; + Mitc4DrillingMode drilling_mode = Mitc4DrillingMode::kDirectSpring; SourceLocation source{}; std::string identity; diff --git a/src/fesa/elements/nonlinear_mitc4_shell.cpp b/src/fesa/elements/nonlinear_mitc4_shell.cpp index d6711cb..e8d41e9 100644 --- a/src/fesa/elements/nonlinear_mitc4_shell.cpp +++ b/src/fesa/elements/nonlinear_mitc4_shell.cpp @@ -12,6 +12,7 @@ #include "fesa/elements/finite_rotation.h" #include "elements/nonlinear_mitc4_configuration_internal.h" +#include "elements/mitc4_drilling_internal.h" #include "math/dense_blas_internal.h" #include "math/jet_arithmetic.h" #include "math/jet_linear_solve.h" @@ -39,6 +40,7 @@ nonlinear_mitc4_internal::ConfigurationSnapshot::Read( result.values[index++] = element.youngs_modulus_; result.values[index++] = element.poisson_ratio_; result.values[index] = element.drilling_stiffness_; + result.drilling_mode = element.drilling_mode_; result.source = element.source_location_; result.identity = element.identity_; return result; @@ -1227,8 +1229,91 @@ Jet DrillSpinRate(const JetVector3& phi, const Vector3& director) { return Reciprocal(Sqrt(Dot(preimage, preimage))); } +/// @brief Uses precisely the current inertia kernel's cubic rotation maps. +/// @note This is deliberately distinct from the legacy analytic null. Keep the +/// cubic Exp/Jacobian order synchronized with EvaluateDynamicRotation; it is +/// not an approximation to be silently replaced by the exact exponential +/// inverse. +JetVector3 InertiaDrillNull(const JetVector3& phi, const Vector3& director) { + const JetMatrix3 skew = Skew(phi); + JetMatrix3 rotation = Identity(); + JetMatrix3 jacobian = Identity(); + JetMatrix3 term = Identity(); + double rotation_factorial = 1.0; + double jacobian_factorial = 1.0; + for (std::size_t order = 1U; order <= 3U; ++order) { + term = Multiply(term, skew); + rotation_factorial *= static_cast(order); + jacobian_factorial *= static_cast(order + 1U); + rotation = Add(rotation, Scale(1.0 / rotation_factorial, term)); + jacobian = Add(jacobian, Scale(1.0 / jacobian_factorial, term)); + } + const JetVector3 rhs = Multiply(rotation, ConstantVector(director)); + Matrix primal{3U, 3U}; + for (std::size_t row = 0U; row < 3U; ++row) + for (std::size_t column = 0U; column < 3U; ++column) + primal(row, column) = jacobian[row][column].value; + const auto factor = dense_blas_internal::DenseFactor::Lu(primal); + const std::array channels{1U, kGlobalDofCount + 1U, + kHessianSize + 1U + kGlobalDofCount, + kGlobalDofCount}; + JetVector3 preimage{}; + if (!factor.HasValue() || + !math_internal::SolveJetAction(jacobian, rhs, factor.Value(), channels, + Constant, ReadChannel, WriteChannel, + preimage)) { + preimage.fill(Constant(std::numeric_limits::quiet_NaN())); + return preimage; + } + return Scale(Reciprocal(Sqrt(Dot(preimage, preimage))), preimage); +} + +/// @brief Evaluates the fixed initial-director spring without corotation. +Jet DirectDrillingEnergy( + const std::array& initial_directors, + const JetTrial& trial, double drilling_stiffness, Matrix* metric, + std::array* constraints) { + Jet energy{}; + double* tangent = energy.second.Data(); + for (std::size_t node = 0U; node < kNodeCount; ++node) { + const auto d0 = ConstantVector(initial_directors[node]); + const Jet gamma = Dot(d0, trial.rotation[node]); + // Retain the validated spring scale through subnormal intermediates. + // Fill its linear derivatives explicitly, without differentiating a + // floating half-energy that can round kd before doubling it again. + energy.value += 0.5 * gamma.value * (drilling_stiffness * gamma.value); + energy.directional_first += + drilling_stiffness * gamma.value * gamma.directional_first; + energy.directional_second += + drilling_stiffness * gamma.directional_first * gamma.directional_first; + for (std::size_t i = 0U; i < 3U; ++i) { + const auto row = 6U * node + 3U + i; + energy.first[row] = + drilling_stiffness * gamma.value * initial_directors[node][i]; + energy.mixed_first[row] = drilling_stiffness * gamma.directional_first * + initial_directors[node][i]; + } + // D2 f_drill[v,v] is identically zero; mixed_second retains exact zeros. + if (constraints != nullptr) { + const auto null = + InertiaDrillNull(trial.rotation[node], initial_directors[node]); + (*constraints)[node] = drilling_stiffness * gamma * Dot(null, d0); + } + for (std::size_t i = 0U; i < 3U; ++i) + for (std::size_t j = 0U; j < 3U; ++j) { + const auto row = 6U * node + 3U + i; + const auto column = 6U * node + 3U + j; + const double value = drilling_stiffness * initial_directors[node][i] * + initial_directors[node][j]; + tangent[row * kGlobalDofCount + column] = value; + (*metric)(row, column) = value; + } + } + return energy; +} + /// @brief Evaluates the objective numerical drilling potential and variations. -Jet DrillingEnergy(const std::array& coordinates, +Jet LegacyDrillingEnergy(const std::array& coordinates, const std::array& initial_directors, const JetTrial& trial, const double thickness, const double drilling_stiffness, bool& valid, Matrix* metric, @@ -1323,6 +1408,7 @@ NonlinearMitc4Shell::NonlinearMitc4Shell( youngs_modulus_{youngs_modulus}, poisson_ratio_{poisson_ratio}, drilling_stiffness_{drilling_stiffness}, + drilling_mode_{mitc4_drilling_internal::kProductionDrillingMode}, source_location_{std::move(source_location)}, identity_{std::move(identity)} {} @@ -1640,7 +1726,10 @@ NonlinearMitc4Shell::EvaluateExternalDrillingConstraint( jet_phi[component] = Variable(phi[component], component, phi_dot[component]); } - const JetVector3 null_direction = DrillNull(jet_phi, initial_directors_[node]); + const JetVector3 null_direction = + drilling_mode_ == Mitc4DrillingMode::kDirectSpring + ? InertiaDrillNull(jet_phi, initial_directors_[node]) + : DrillNull(jet_phi, initial_directors_[node]); const Jet scalar = -1.0 * Dot(null_direction, ConstantVector(force)); const Jet scalar_time = -1.0 * Dot(null_direction, ConstantVector(force_rate)); @@ -1741,10 +1830,17 @@ Result EvaluateInternalResidualAndTangent( element.poisson_ratio_, valid); Matrix drilling_metric{kGlobalDofCount, kGlobalDofCount}; std::array constraints{}; - const Jet drilling = DrillingEnergy( - element.coordinates_, element.initial_directors_, jet_trial, - element.thickness_, element.drilling_stiffness_, valid, &drilling_metric, - coordinate_direction != nullptr ? &constraints : nullptr); + const Jet drilling = + element.drilling_mode_ == Mitc4DrillingMode::kDirectSpring + ? DirectDrillingEnergy( + element.initial_directors_, jet_trial, + element.drilling_stiffness_, &drilling_metric, + coordinate_direction != nullptr ? &constraints : nullptr) + : LegacyDrillingEnergy( + element.coordinates_, element.initial_directors_, jet_trial, + element.thickness_, element.drilling_stiffness_, valid, + &drilling_metric, + coordinate_direction != nullptr ? &constraints : nullptr); const Jet internal_energy = physical + drilling; if (!valid || !IsFinite(internal_energy)) { return Failure( @@ -1815,10 +1911,13 @@ Result EvaluateInternalResidualAndTangent( equilibrium.drilling_metric = std::move(drilling_metric); equilibrium.drilling_energy = drilling.value; equilibrium.force_directional_second = Vector{kGlobalDofCount}; + equilibrium.drilling_force_directional_second = Vector{kGlobalDofCount}; for (std::size_t index = 0U; index < kGlobalDofCount; ++index) { equilibrium.internal_force[index] = equilibrium.internal_residual[index]; equilibrium.physical_force[index] = physical.first[index]; equilibrium.drilling_force[index] = drilling.first[index]; + equilibrium.drilling_force_directional_second[index] = + drilling.mixed_second[index]; equilibrium.force_directional_second[index] = internal_energy.mixed_second[index]; for (std::size_t column = 0U; column < kGlobalDofCount; ++column) { diff --git a/src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp b/src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp index db70de4..67bc522 100644 --- a/src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/dynamic_hdf5_results_writer.cpp @@ -18,6 +18,10 @@ #include "io/hdf5/hdf5_primitives.h" #include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "io/hdf5/mitc4_file_hooks.h" + namespace fesa { namespace { @@ -849,13 +853,13 @@ void CheckTables(const hid_t file, const DynamicAnalysisState& state, Status SelfCheck(const std::filesystem::path& path, const Domain& domain, const DynamicAnalysisState& state, const hdf5_internal::WriterModelData& model, - const std::vector& diagnostics) { + const std::vector& diagnostics, herr_t (*close_reader)(hid_t)) { try { Hdf5Handle file{ hdf5_internal::RequireHdf5Id( H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), "Unable to reopen dynamic candidate."), - H5Fclose}; + close_reader}; Hdf5Handle metadata{hdf5_internal::RequireHdf5Id( H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), "The dynamic metadata group is missing."), @@ -902,9 +906,16 @@ Status SelfCheck(const std::filesystem::path& path, const Domain& domain, "Unable to inspect accepted frame group."); Check(information.nlinks == state.AcceptedFrames().size(), "Unexpected frame group inventory."); + const auto drilling = hdf5_internal::RequireDrillingDiagnostics( + file.Get(), state.AcceptedFrames().front().mitc4_drilling->mode); + hdf5_internal::RequireDrilling(drilling.frames.size() == state.AcceptedFrames().size()); + std::size_t drilling_index = 0U; for (const auto& frame : state.AcceptedFrames()) { + hdf5_internal::CheckDrillingFrame(drilling.frames[drilling_index++], frame.frame_index, + frame.global_equilibrium, *frame.mitc4_drilling); WriteFrame(file.Get(), domain, frame, true); } + hdf5_internal::RequireHdf5(file.CloseChecked(), "Unable to close dynamic readback."); return Status::Ok(); } catch (const std::exception& error) { return Failure("dynamic-output-schema", error.what()); @@ -913,10 +924,10 @@ Status SelfCheck(const std::filesystem::path& path, const Domain& domain, } // namespace -Status DynamicHdf5ResultsWriter::Write( +Status hdf5_internal::WriteDynamicFile( const std::filesystem::path& output_path, const Domain& domain, const DynamicAnalysisState& state, - const std::vector& diagnostics) { + const std::vector& diagnostics, const Mitc4FileHooks& hooks) { if (output_path.empty() || output_path.filename().empty() || !state.IsInitialized() || state.AcceptedFrames().empty() || domain.Dynamic() == nullptr) @@ -924,6 +935,11 @@ Status DynamicHdf5ResultsWriter::Write( "Dynamic output requires a named path, initialized state, " "and accepted frames."); try { + for (const auto& frame : state.AcceptedFrames()) { + hdf5_internal::ValidateDrillingEvidence(frame.global_equilibrium, frame.mitc4_drilling); + hdf5_internal::RequireDrilling(frame.mitc4_drilling->mode == + state.AcceptedFrames().front().mitc4_drilling->mode); + } hdf5_internal::Hdf5ErrorSilencer silence; const auto& step = *domain.Dynamic(); Check(state.AcceptedFrames().size() == @@ -965,24 +981,28 @@ Status DynamicHdf5ResultsWriter::Write( H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT), "Unable to create dynamic candidate."), - H5Fclose}; + hooks.candidate_close}; (void)hdf5_internal::CreateGroup(file.Get(), "/model/shell"); hdf5_internal::WriteModel(file.Get(), domain, model); WriteMetadata(file.Get(), domain, state.AcceptedFrames().front()); + hdf5_internal::WriteDrillingMetadata(file.Get(), state.AcceptedFrames().front().mitc4_drilling->mode); WriteDensity(file.Get(), domain); (void)hdf5_internal::CreateGroup(file.Get(), "/steps"); (void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1"); WriteAcceptedFrames(file.Get(), state.AcceptedFrames()); for (const auto& frame : state.AcceptedFrames()) { WriteFrame(file.Get(), domain, frame); + hdf5_internal::WriteDrillingFrame(file.Get(), + "/steps/Step-1/frames/" + std::to_string(frame.frame_index), *frame.mitc4_drilling); } WriteDiagnostics(file.Get(), diagnostics); hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), "Unable to flush dynamic candidate."); hdf5_internal::RequireHdf5(file.CloseChecked(), "Unable to close dynamic candidate."); + if (hooks.before_self_check) hooks.before_self_check(candidate_path); const Status check = - SelfCheck(candidate_path, domain, state, model, diagnostics); + SelfCheck(candidate_path, domain, state, model, diagnostics, hooks.reader_close); const Status finalization = hdf5_internal::AtomicFinalizeValidatedCandidate( candidate_path, output_path, check); if (finalization.IsOk()) cleanup.Release(); @@ -994,4 +1014,10 @@ Status DynamicHdf5ResultsWriter::Write( } } +Status DynamicHdf5ResultsWriter::Write( + const std::filesystem::path& output_path, const Domain& domain, + const DynamicAnalysisState& state, const std::vector& diagnostics) { + return hdf5_internal::WriteDynamicFile(output_path, domain, state, diagnostics); +} + } // namespace fesa diff --git a/src/fesa/io/hdf5/eigen_hdf5_results_writer.cpp b/src/fesa/io/hdf5/eigen_hdf5_results_writer.cpp index 60688de..dc6ee25 100644 --- a/src/fesa/io/hdf5/eigen_hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/eigen_hdf5_results_writer.cpp @@ -22,6 +22,10 @@ #include "io/hdf5/hdf5_raii.h" #include "math/finite_values.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "io/hdf5/mitc4_file_hooks.h" + namespace fesa { namespace { @@ -875,13 +879,14 @@ void WriteDiagnostics(hid_t file, const std::vector& diagnostics) { } // namespace -Status EigenHdf5ResultsWriter::Write( +Status hdf5_internal::WriteEigenFile( const std::filesystem::path& output_path, const Domain& domain, const EigenAnalysisState& state, - const std::vector& diagnostics) { + const std::vector& diagnostics, const Mitc4FileHooks& hooks) { const Status valid = Validate(domain, state, diagnostics); if (!valid.IsOk()) return valid; try { + hdf5_internal::RequireDrilling(state.drilling_mode.has_value()); hdf5_internal::Hdf5ErrorSilencer silence; const auto candidate_path = hdf5_internal::MakeCandidatePath(output_path); hdf5_internal::CandidateFileGuard cleanup{candidate_path}; @@ -890,7 +895,7 @@ Status EigenHdf5ResultsWriter::Write( H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT), "Unable to create eigen candidate."), - H5Fclose}; + hooks.candidate_close}; hdf5_internal::WriterModelData model{ {}, state.constraint_mask, state.prescribed_displacement}; (void)hdf5_internal::CreateGroup(file.Get(), "/model"); @@ -899,6 +904,7 @@ Status EigenHdf5ResultsWriter::Write( WriteDensities(file.Get(), domain); WriteDiagnostics(file.Get(), diagnostics); WriteEigenMetadata(file.Get()); + hdf5_internal::WriteDrillingMetadata(file.Get(), *state.drilling_mode); WriteStep(file.Get(), domain, state); (void)hdf5_internal::CreateGroup(file.Get(), kEigenRoot); WriteModes(file.Get(), state, domain.Frequency()->source_step_name_raw); @@ -940,8 +946,9 @@ Status EigenHdf5ResultsWriter::Write( "Unable to flush eigen candidate."); hdf5_internal::RequireHdf5(file.CloseChecked(), "Unable to close eigen candidate."); + if (hooks.before_self_check) hooks.before_self_check(candidate_path); const Status checked = hdf5_internal::SelfCheckEigenFile( - candidate_path, domain, state, diagnostics); + candidate_path, domain, state, diagnostics, hooks.reader_close); const Status result = hdf5_internal::AtomicFinalizeValidatedCandidate( candidate_path, output_path, checked); if (result.IsOk()) cleanup.Release(); @@ -953,4 +960,10 @@ Status EigenHdf5ResultsWriter::Write( } } +Status EigenHdf5ResultsWriter::Write( + const std::filesystem::path& output_path, const Domain& domain, + const EigenAnalysisState& state, const std::vector& diagnostics) { + return hdf5_internal::WriteEigenFile(output_path, domain, state, diagnostics); +} + } // namespace fesa diff --git a/src/fesa/io/hdf5/eigen_hdf5_self_check.cpp b/src/fesa/io/hdf5/eigen_hdf5_self_check.cpp index b41907e..9a923a2 100644 --- a/src/fesa/io/hdf5/eigen_hdf5_self_check.cpp +++ b/src/fesa/io/hdf5/eigen_hdf5_self_check.cpp @@ -15,6 +15,8 @@ #include "io/hdf5/hdf5_primitives.h" #include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + namespace fesa::hdf5_internal { namespace { @@ -726,13 +728,21 @@ bool CheckModalRows(hid_t file, const Domain& domain, Status SelfCheckEigenFile(const std::filesystem::path& path, const Domain& domain, const EigenAnalysisState& state, - const std::vector& diagnostics) { + const std::vector& diagnostics, + herr_t (*close_reader)(hid_t)) { try { Hdf5ErrorSilencer silence; Hdf5Handle file{RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), "Unable to reopen eigen candidate."), - H5Fclose}; + close_reader}; + const auto drilling = ReadDrillingDiagnostics(file.Get()); + if (state.drilling_mode) { + RequireDrilling(drilling.has_value() && drilling->formulation == + Mitc4DrillingModeName(*state.drilling_mode) && drilling->frames.empty()); + } else { + RequireDrilling(!drilling.has_value()); + } for (const char* required : {"/metadata", "/model/nodes", "/model/elements", "/model/shell/densities", "/diagnostics", "/steps/Step-1/eigen/modes", diff --git a/src/fesa/io/hdf5/eigen_hdf5_self_check.h b/src/fesa/io/hdf5/eigen_hdf5_self_check.h index da7458e..2cbd068 100644 --- a/src/fesa/io/hdf5/eigen_hdf5_self_check.h +++ b/src/fesa/io/hdf5/eigen_hdf5_self_check.h @@ -1,6 +1,8 @@ #ifndef FESA_IO_HDF5_EIGEN_HDF5_SELF_CHECK_H_ #define FESA_IO_HDF5_EIGEN_HDF5_SELF_CHECK_H_ +#include + #include #include @@ -17,7 +19,8 @@ namespace fesa::hdf5_internal { /// @brief Reopens and validates a closed authoritative eigen HDF5 candidate. Status SelfCheckEigenFile(const std::filesystem::path& path, const Domain& domain, const EigenAnalysisState& state, - const std::vector& diagnostics); + const std::vector& diagnostics, + herr_t (*close_reader)(hid_t) = H5Fclose); } // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/hdf5_results_writer.cpp b/src/fesa/io/hdf5/hdf5_results_writer.cpp index 077ca81..4640766 100644 --- a/src/fesa/io/hdf5/hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/hdf5_results_writer.cpp @@ -12,12 +12,16 @@ #include "io/hdf5/hdf5_result_writer.h" #include "io/hdf5/hdf5_self_check.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "io/hdf5/mitc4_file_hooks.h" + namespace fesa { -Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, +Status hdf5_internal::WriteStaticFile(const std::filesystem::path& output_path, const Domain& domain, const AnalysisState& state, - const std::vector& diagnostics) { + const std::vector& diagnostics, const Mitc4FileHooks& hooks) { auto plan_result = hdf5_internal::BuildWriterPlan(output_path, domain, state, diagnostics); if (!plan_result.HasValue()) { @@ -26,6 +30,8 @@ Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, hdf5_internal::WriterPlan plan = std::move(plan_result.Value()); try { + if (hdf5_internal::IsShellDomain(domain)) + hdf5_internal::ValidateDrillingEvidence(state.Equilibrium(), state.Mitc4Drilling()); hdf5_internal::Hdf5ErrorSilencer silence_backend_errors; const std::filesystem::path candidate_path = hdf5_internal::MakeCandidatePath(output_path); @@ -35,7 +41,7 @@ Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT), "Unable to create the temporary HDF5 file."), - H5Fclose}; + hooks.candidate_close}; (void)hdf5_internal::CreateGroup(file.Get(), "/model"); (void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1/frames/0/nodal"); @@ -50,14 +56,19 @@ Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, } hdf5_internal::WriteModel(file.Get(), domain, plan.model_data); hdf5_internal::WriteResults(file.Get(), domain, state, diagnostics, plan); + if (hdf5_internal::IsShellDomain(domain)) { + hdf5_internal::WriteDrillingMetadata(file.Get(), state.Mitc4Drilling()->mode); + hdf5_internal::WriteDrillingFrame(file.Get(), hdf5_internal::kStepRoot, *state.Mitc4Drilling()); + } hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), "Unable to flush the temporary HDF5 file."); hdf5_internal::RequireHdf5( file.CloseChecked(), "Unable to close the temporary HDF5 file after writing."); + if (hooks.before_self_check) hooks.before_self_check(candidate_path); const Status self_check = hdf5_internal::SelfCheck( - candidate_path, domain, state, diagnostics.size(), plan.model_data); + candidate_path, domain, state, diagnostics.size(), plan.model_data, hooks.reader_close); const Status atomic = hdf5_internal::AtomicFinalizeValidatedCandidate( candidate_path, output_path, self_check); if (atomic.IsOk()) { @@ -71,4 +82,10 @@ Status Hdf5ResultsWriter::Write(const std::filesystem::path& output_path, } } +Status Hdf5ResultsWriter::Write( + const std::filesystem::path& output_path, const Domain& domain, + const AnalysisState& state, const std::vector& diagnostics) { + return hdf5_internal::WriteStaticFile(output_path, domain, state, diagnostics); +} + } // namespace fesa diff --git a/src/fesa/io/hdf5/hdf5_self_check.cpp b/src/fesa/io/hdf5/hdf5_self_check.cpp index 2779079..95e8973 100644 --- a/src/fesa/io/hdf5/hdf5_self_check.cpp +++ b/src/fesa/io/hdf5/hdf5_self_check.cpp @@ -13,6 +13,8 @@ #include "io/hdf5/hdf5_raii.h" #include "io/hdf5/hdf5_primitives.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + namespace fesa::hdf5_internal { Hdf5Handle OpenDatasetForCheck(const hid_t file, const char* path) { @@ -516,18 +518,30 @@ void RequireSchemaV1ShellIdentity(const hid_t file, const char* path, void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, const std::size_t diagnostic_count, - const WriterModelData& model_data) { + const WriterModelData& model_data, + herr_t (*close_reader)(hid_t)) { if (H5Fis_hdf5(path.string().c_str()) <= 0) { throw Hdf5Failure{"The temporary output is not an HDF5 file."}; } Hdf5Handle file{ RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), "Unable to reopen the temporary HDF5 file read-only."), - H5Fclose}; + close_reader}; Hdf5Handle metadata{ RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), "The HDF5 metadata group is missing."), H5Gclose}; + if (IsShellDomain(domain)) { + const auto drilling = ReadDrillingDiagnostics(file.Get()); + if (state.Mitc4Drilling()) { + ValidateDrillingEvidence(state.Equilibrium(), state.Mitc4Drilling()); + RequireDrilling(drilling.has_value() && drilling->formulation == + Mitc4DrillingModeName(state.Mitc4Drilling()->mode) && drilling->frames.size() == 1U); + CheckDrillingFrame(drilling->frames.front(), 0U, state.Equilibrium(), *state.Mitc4Drilling()); + } else { + RequireDrilling(!drilling.has_value()); + } + } RequireUint64Attribute(metadata.Get(), "schema_version", 0U); RequireStringAttribute(metadata.Get(), "solver_version", std::string{SolverVersion()}); @@ -751,18 +765,36 @@ void SelfCheckFile(const std::filesystem::path& path, const Domain& domain, void NonlinearSelfCheckFile(const std::filesystem::path& path, const Domain& domain, const NonlinearAnalysisState& state, - const std::size_t diagnostic_count) { + const std::size_t diagnostic_count, + herr_t (*close_reader)(hid_t)) { if (H5Fis_hdf5(path.string().c_str()) <= 0) { throw Hdf5Failure{"The nonlinear candidate is not an HDF5 file."}; } Hdf5Handle file{ RequireHdf5Id(H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), "Unable to reopen nonlinear candidate read-only."), - H5Fclose}; + close_reader}; Hdf5Handle metadata{ RequireHdf5Id(H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), "The nonlinear metadata group is missing."), H5Gclose}; + RequireDrilling(!state.AcceptedFrames().empty()); + const auto drilling = ReadDrillingDiagnostics(file.Get()); + const auto& expected_drilling = state.AcceptedFrames().front().mitc4_drilling; + if (expected_drilling) { + RequireDrilling(drilling.has_value() && drilling->formulation == + Mitc4DrillingModeName(expected_drilling->mode) && + drilling->frames.size() == state.AcceptedFrames().size()); + for (std::size_t i = 0; i < drilling->frames.size(); ++i) { + const auto& frame = state.AcceptedFrames()[i]; + ValidateDrillingEvidence(frame.global_equilibrium, frame.mitc4_drilling); + RequireDrilling(frame.mitc4_drilling->mode == expected_drilling->mode); + CheckDrillingFrame(drilling->frames[i], frame.frame_index, frame.global_equilibrium, *frame.mitc4_drilling); + } + } else { + RequireDrilling(!drilling.has_value()); + for (const auto& frame : state.AcceptedFrames()) RequireDrilling(!frame.mitc4_drilling); + } RequireUint64Attribute(metadata.Get(), "schema_version", 1U); RequireStringAttribute(metadata.Get(), "feature_id", "nonlinear-static-newton-raphson"); @@ -919,9 +951,10 @@ void NonlinearSelfCheckFile(const std::filesystem::path& path, Status SelfCheck(const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, const std::size_t diagnostic_count, - const WriterModelData& model_data) { + const WriterModelData& model_data, + herr_t (*close_reader)(hid_t)) { try { - SelfCheckFile(path, domain, state, diagnostic_count, model_data); + SelfCheckFile(path, domain, state, diagnostic_count, model_data, close_reader); return Status::Ok(); } catch (const Hdf5Failure& failure) { return OutputFailure("hdf5-write-failure", failure.what()); @@ -933,9 +966,10 @@ Status SelfCheck(const std::filesystem::path& path, const Domain& domain, Status NonlinearSelfCheck(const std::filesystem::path& path, const Domain& domain, const NonlinearAnalysisState& state, - const std::size_t diagnostic_count) { + const std::size_t diagnostic_count, + herr_t (*close_reader)(hid_t)) { try { - NonlinearSelfCheckFile(path, domain, state, diagnostic_count); + NonlinearSelfCheckFile(path, domain, state, diagnostic_count, close_reader); return Status::Ok(); } catch (const Hdf5Failure& failure) { return OutputFailure("hdf5-write-failure", failure.what()); diff --git a/src/fesa/io/hdf5/hdf5_self_check.h b/src/fesa/io/hdf5/hdf5_self_check.h index be6fb94..7369eda 100644 --- a/src/fesa/io/hdf5/hdf5_self_check.h +++ b/src/fesa/io/hdf5/hdf5_self_check.h @@ -15,13 +15,15 @@ namespace fesa::hdf5_internal { /// @brief Reopens and validates the complete closed schema-v0 candidate. Status SelfCheck(const std::filesystem::path& path, const Domain& domain, const AnalysisState& state, std::size_t diagnostic_count, - const WriterModelData& model_data); + const WriterModelData& model_data, + herr_t (*close_reader)(hid_t) = H5Fclose); /// @brief Reopens and validates the complete closed schema-v1 candidate. Status NonlinearSelfCheck(const std::filesystem::path& path, const Domain& domain, const NonlinearAnalysisState& state, - std::size_t diagnostic_count); + std::size_t diagnostic_count, + herr_t (*close_reader)(hid_t) = H5Fclose); } // namespace fesa::hdf5_internal diff --git a/src/fesa/io/hdf5/mitc4_drilling_diagnostics.h b/src/fesa/io/hdf5/mitc4_drilling_diagnostics.h new file mode 100644 index 0000000..5a0f826 --- /dev/null +++ b/src/fesa/io/hdf5/mitc4_drilling_diagnostics.h @@ -0,0 +1,172 @@ +#ifndef FESA_IO_HDF5_MITC4_DRILLING_DIAGNOSTICS_H_ +#define FESA_IO_HDF5_MITC4_DRILLING_DIAGNOSTICS_H_ + +#include +#include +#include +#include +#include +#include + +#include "fesa/results/result_records.h" +#include "io/hdf5/hdf5_primitives.h" + +namespace fesa::hdf5_internal { + +struct StoredDrillingFrame { + std::string storage_step; + std::uint64_t frame_index; + std::array raw_equilibrium; + std::array drilling_moment; + std::array equilibrium_with_drilling; +}; + +struct StoredDrillingDiagnostics { + std::string formulation; + std::uint64_t version; + std::vector frames; +}; + +inline void RequireDrilling(bool valid) { + if (!valid) throw Hdf5Failure{"Invalid MITC4 drilling diagnostics contract."}; +} + +inline void ValidateDrillingEvidence( + const std::array& raw, + const std::optional& evidence) { + RequireDrilling(evidence.has_value()); + RequireDrilling(Mitc4DrillingModeName(evidence->mode)[0] != '\0' && + IsFinite(raw) && IsFinite(evidence->drilling_moment) && + IsFinite(evidence->equilibrium_with_drilling)); + for (std::size_t i = 0; i < 6U; ++i) { + RequireDrilling(evidence->equilibrium_with_drilling[i] == + (i < 3U ? raw[i] : raw[i] - evidence->drilling_moment[i - 3U])); + } +} + +inline void WriteDrillingMetadata(hid_t file, Mitc4DrillingMode mode) { + RequireDrilling(Mitc4DrillingModeName(mode)[0] != '\0'); + Hdf5Handle metadata{RequireHdf5Id(H5Gopen2(file, "/metadata", H5P_DEFAULT), + "Missing drilling metadata group."), H5Gclose}; + WriteStringAttribute(metadata.Get(), "mitc4_drilling_formulation", + Mitc4DrillingModeName(mode)); + WriteUint64Attribute(metadata.Get(), "drilling_diagnostics_version", 1U); +} + +inline void WriteDrillingFrame(hid_t file, const std::string& root, + const Mitc4DrillingEvidence& evidence) { + WriteModelDoubleDataset(file, root + "/global/drilling_moment", {3U}, + evidence.drilling_moment.data(), 3U, "MDRILL_X,MDRILL_Y,MDRILL_Z", + "force*length,force*length,force*length", "global-cartesian", "global-origin"); + WriteModelDoubleDataset(file, root + "/global/equilibrium_with_drilling", {6U}, + evidence.equilibrium_with_drilling.data(), 6U, + "FORCE_1,FORCE_2,FORCE_3,MOMENT_1,MOMENT_2,MOMENT_3", + "force,force,force,force*length,force*length,force*length", + "global-cartesian", "global-origin"); +} + +inline std::string ReadDrillingString(hid_t object, const char* name) { + Hdf5Handle attribute{RequireHdf5Id(H5Aopen(object, name, H5P_DEFAULT), + "Missing drilling string."), H5Aclose}; + Hdf5Handle type{H5Aget_type(attribute.Get()), H5Tclose}; + Hdf5Handle space{H5Aget_space(attribute.Get()), H5Sclose}; + RequireDrilling(H5Sget_simple_extent_type(space.Get()) == H5S_SCALAR && + H5Tget_class(type.Get()) == H5T_STRING && + H5Tis_variable_str(type.Get()) > 0 && + H5Tget_cset(type.Get()) == H5T_CSET_UTF8); + char* value = nullptr; + RequireHdf5(H5Aread(attribute.Get(), type.Get(), &value), "Cannot read drilling string."); + const std::string result = value == nullptr ? "" : value; + if (value != nullptr) H5free_memory(value); + RequireDrilling(IsValidUtf8(result)); + return result; +} + +template +inline std::array ReadDrillingValues(hid_t file, const std::string& path) { + Hdf5Handle dataset{RequireHdf5Id(H5Dopen2(file, path.c_str(), H5P_DEFAULT), + "Missing drilling dataset."), H5Dclose}; + Hdf5Handle type{H5Dget_type(dataset.Get()), H5Tclose}; + Hdf5Handle space{H5Dget_space(dataset.Get()), H5Sclose}; + hsize_t dimension = 0U; + RequireDrilling(H5Tequal(type.Get(), H5T_IEEE_F64LE) > 0 && + H5Sget_simple_extent_ndims(space.Get()) == 1 && + H5Sget_simple_extent_dims(space.Get(), &dimension, nullptr) >= 0 && + dimension == N); + std::array values{}; + RequireHdf5(H5Dread(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, values.data()), "Cannot read drilling dataset."); + RequireDrilling(IsFinite(values)); + return values; +} + +// Markerless historical files retain their old contract; no mode is inferred. +inline std::optional ReadDrillingDiagnostics(hid_t file) { + Hdf5Handle metadata{RequireHdf5Id(H5Gopen2(file, "/metadata", H5P_DEFAULT), + "Missing drilling metadata group."), H5Gclose}; + const htri_t mode = H5Aexists(metadata.Get(), "mitc4_drilling_formulation"); + const htri_t version = H5Aexists(metadata.Get(), "drilling_diagnostics_version"); + RequireDrilling(mode >= 0 && version >= 0 && (mode > 0) == (version > 0)); + if (mode == 0) return std::nullopt; + StoredDrillingDiagnostics result{}; + result.formulation = ReadDrillingString(metadata.Get(), "mitc4_drilling_formulation"); + RequireDrilling(result.formulation == Mitc4DrillingModeName(Mitc4DrillingMode::kDirectSpring) || + result.formulation == Mitc4DrillingModeName(Mitc4DrillingMode::kLegacyCorotatedTwist)); + Hdf5Handle attribute{H5Aopen(metadata.Get(), "drilling_diagnostics_version", H5P_DEFAULT), H5Aclose}; + Hdf5Handle type{H5Aget_type(attribute.Get()), H5Tclose}; + Hdf5Handle space{H5Aget_space(attribute.Get()), H5Sclose}; + RequireDrilling(H5Tequal(type.Get(), H5T_STD_U64LE) > 0 && + H5Sget_simple_extent_type(space.Get()) == H5S_SCALAR); + RequireHdf5(H5Aread(attribute.Get(), H5T_NATIVE_UINT64, &result.version), "Cannot read drilling version."); + RequireDrilling(result.version == 1U); + if (H5Lexists(file, "/steps/Step-1/frames", H5P_DEFAULT) > 0) { + Hdf5Handle frames{H5Gopen2(file, "/steps/Step-1/frames", H5P_DEFAULT), H5Gclose}; + H5G_info_t info{}; + RequireHdf5(H5Gget_info(frames.Get(), &info), "Cannot inspect drilling frames."); + for (hsize_t i = 0; i < info.nlinks; ++i) { + const auto size = H5Lget_name_by_idx(frames.Get(), ".", H5_INDEX_NAME, H5_ITER_INC, + i, nullptr, 0U, H5P_DEFAULT); + RequireDrilling(size > 0); + std::vector name(static_cast(size) + 1U); + RequireDrilling(H5Lget_name_by_idx(frames.Get(), ".", H5_INDEX_NAME, H5_ITER_INC, + i, name.data(), name.size(), H5P_DEFAULT) == size); + const std::string token{name.data()}; + const std::string root = "/steps/Step-1/frames/" + token + "/global/"; + // Each new MITC4 static/dynamic frame has an existing raw equilibrium. + StoredDrillingFrame frame{}; + frame.storage_step = "Step-1"; + RequireDrilling(!token.empty() && token.find_first_not_of("0123456789") == std::string::npos); + frame.frame_index = std::stoull(token); + RequireDrilling(std::to_string(frame.frame_index) == token); + frame.raw_equilibrium = ReadDrillingValues<6>(file, root + "equilibrium"); + frame.drilling_moment = ReadDrillingValues<3>(file, root + "drilling_moment"); + frame.equilibrium_with_drilling = ReadDrillingValues<6>(file, root + "equilibrium_with_drilling"); + for (std::size_t j = 0; j < 6U; ++j) + RequireDrilling(frame.equilibrium_with_drilling[j] == (j < 3U ? frame.raw_equilibrium[j] : + frame.raw_equilibrium[j] - frame.drilling_moment[j - 3U])); + result.frames.push_back(frame); + } + std::sort(result.frames.begin(), result.frames.end(), + [](const auto& left, const auto& right) { return left.frame_index < right.frame_index; }); + } + return result; +} + +inline StoredDrillingDiagnostics RequireDrillingDiagnostics(hid_t file, Mitc4DrillingMode mode) { + const auto result = ReadDrillingDiagnostics(file); + RequireDrilling(result.has_value() && result->formulation == Mitc4DrillingModeName(mode)); + return *result; +} + +inline void CheckDrillingFrame(const StoredDrillingFrame& stored, + std::uint64_t index, + const std::array& raw, + const Mitc4DrillingEvidence& evidence) { + RequireDrilling(stored.frame_index == index && stored.raw_equilibrium == raw && + stored.drilling_moment == evidence.drilling_moment && + stored.equilibrium_with_drilling == evidence.equilibrium_with_drilling); +} + +} // namespace fesa::hdf5_internal + +#endif // FESA_IO_HDF5_MITC4_DRILLING_DIAGNOSTICS_H_ diff --git a/src/fesa/io/hdf5/mitc4_file_hooks.h b/src/fesa/io/hdf5/mitc4_file_hooks.h new file mode 100644 index 0000000..4e629b7 --- /dev/null +++ b/src/fesa/io/hdf5/mitc4_file_hooks.h @@ -0,0 +1,42 @@ +#ifndef FESA_IO_HDF5_MITC4_FILE_HOOKS_H_ +#define FESA_IO_HDF5_MITC4_FILE_HOOKS_H_ + +#include + +#include +#include +#include + +#include "fesa/core/status.h" + +namespace fesa { +class Domain; +class AnalysisState; +class NonlinearAnalysisState; +class DynamicAnalysisState; +struct EigenAnalysisState; +struct Diagnostic; + +namespace hdf5_internal { + +// Private lifecycle fault seams. Production uses real closers and no callback. +// Injected closers must release the handle even when reporting failure. +struct Mitc4FileHooks { + herr_t (*candidate_close)(hid_t) = H5Fclose; + herr_t (*reader_close)(hid_t) = H5Fclose; + std::function before_self_check; +}; + +Status WriteStaticFile(const std::filesystem::path&, const Domain&, + const AnalysisState&, const std::vector&, const Mitc4FileHooks& = {}); +Status WriteNonlinearFile(const std::filesystem::path&, const Domain&, + const NonlinearAnalysisState&, const std::vector&, const Mitc4FileHooks& = {}); +Status WriteDynamicFile(const std::filesystem::path&, const Domain&, + const DynamicAnalysisState&, const std::vector&, const Mitc4FileHooks& = {}); +Status WriteEigenFile(const std::filesystem::path&, const Domain&, + const EigenAnalysisState&, const std::vector&, const Mitc4FileHooks& = {}); + +} // namespace hdf5_internal +} // namespace fesa + +#endif // FESA_IO_HDF5_MITC4_FILE_HOOKS_H_ diff --git a/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp b/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp index cadd201..481dcb2 100644 --- a/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp +++ b/src/fesa/io/hdf5/nonlinear_hdf5_results_writer.cpp @@ -22,6 +22,10 @@ #include "math/finite_values.h" #include "results/shell_result_recovery.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "io/hdf5/mitc4_file_hooks.h" + namespace fesa { namespace { @@ -955,10 +959,10 @@ Status ValidateDiagnostics(const std::vector& diagnostics) { } // namespace -Status NonlinearHdf5ResultsWriter::Write( +Status hdf5_internal::WriteNonlinearFile( const std::filesystem::path& output_path, const Domain& domain, const NonlinearAnalysisState& state, - const std::vector& diagnostics) { + const std::vector& diagnostics, const Mitc4FileHooks& hooks) { if (output_path.empty() || output_path.filename().empty()) { return Failure("invalid-output-path", "The HDF5 output path must name a file."); } @@ -981,6 +985,11 @@ Status NonlinearHdf5ResultsWriter::Write( } try { + for (const auto& frame : state.AcceptedFrames()) { + hdf5_internal::ValidateDrillingEvidence(frame.global_equilibrium, frame.mitc4_drilling); + hdf5_internal::RequireDrilling(frame.mitc4_drilling->mode == + state.AcceptedFrames().front().mitc4_drilling->mode); + } hdf5_internal::Hdf5ErrorSilencer silence_backend_errors; const std::filesystem::path candidate_path = hdf5_internal::MakeCandidatePath(output_path); @@ -989,25 +998,29 @@ Status NonlinearHdf5ResultsWriter::Write( H5Fcreate(candidate_path.string().c_str(), H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT), "Unable to create nonlinear temporary HDF5 file."), - H5Fclose}; + hooks.candidate_close}; (void)hdf5_internal::CreateGroup(file.Get(), "/model"); (void)hdf5_internal::CreateGroup(file.Get(), "/model/shell"); (void)hdf5_internal::CreateGroup(file.Get(), "/steps/Step-1"); WriteMetadata(file.Get(), domain); + hdf5_internal::WriteDrillingMetadata(file.Get(), state.AcceptedFrames().front().mitc4_drilling->mode); WriteNodes(file.Get(), domain); WriteElements(file.Get(), domain); WriteModelData(file.Get(), domain, model_data); WriteAcceptedFrameTable(file.Get(), state.AcceptedFrames()); for (const auto& frame : state.AcceptedFrames()) { WriteFrame(file.Get(), domain, frame); + hdf5_internal::WriteDrillingFrame(file.Get(), + "/steps/Step-1/frames/" + std::to_string(frame.frame_index), *frame.mitc4_drilling); } WriteDiagnostics(file.Get(), diagnostics); hdf5_internal::RequireHdf5(H5Fflush(file.Get(), H5F_SCOPE_GLOBAL), "Unable to flush nonlinear temporary HDF5 file."); hdf5_internal::RequireHdf5( file.CloseChecked(), "Unable to close nonlinear temporary HDF5 file."); + if (hooks.before_self_check) hooks.before_self_check(candidate_path); const Status self_check = hdf5_internal::NonlinearSelfCheck( - candidate_path, domain, state, diagnostics.size()); + candidate_path, domain, state, diagnostics.size(), hooks.reader_close); const Status atomic = hdf5_internal::AtomicFinalizeValidatedCandidate( candidate_path, output_path, self_check); if (atomic.IsOk()) { @@ -1021,4 +1034,10 @@ Status NonlinearHdf5ResultsWriter::Write( } } +Status NonlinearHdf5ResultsWriter::Write( + const std::filesystem::path& output_path, const Domain& domain, + const NonlinearAnalysisState& state, const std::vector& diagnostics) { + return hdf5_internal::WriteNonlinearFile(output_path, domain, state, diagnostics); +} + } // namespace fesa diff --git a/src/fesa/results/analysis_state_commit.cpp b/src/fesa/results/analysis_state_commit.cpp index df6ed96..6058e83 100644 --- a/src/fesa/results/analysis_state_commit.cpp +++ b/src/fesa/results/analysis_state_commit.cpp @@ -48,7 +48,8 @@ Status CommitRecoveryCandidate(RecoveryCandidate candidate, !candidate.shell.rows.empty() || candidate.shell.physical_strain_energy != 0.0 || candidate.shell.equilibrium != zero_equilibrium || - candidate.shell.verification_metrics != zero_metrics; + candidate.shell.verification_metrics != zero_metrics || + candidate.shell.mitc4_drilling.has_value(); const bool has_complete_shell_inventory = !candidate.expected_shell_elements.empty() && !candidate.shell.rows.empty(); diff --git a/src/fesa/results/global_equilibrium_recovery.cpp b/src/fesa/results/global_equilibrium_recovery.cpp index 3f36bca..cd44b7b 100644 --- a/src/fesa/results/global_equilibrium_recovery.cpp +++ b/src/fesa/results/global_equilibrium_recovery.cpp @@ -10,6 +10,7 @@ #include #include "fesa/math/vector3.h" +#include "results/mitc4_drilling_evidence.h" namespace fesa::results_internal { namespace { @@ -243,11 +244,24 @@ Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs, candidate.shell.equilibrium[component] = force_balance[component]; candidate.shell.equilibrium[3U + component] = moment_balance[component]; } + Vector3 acceptance_moment = moment_balance; + if (candidate.shell.mitc4_drilling) { + CompleteDrillingEvidence(candidate.shell.equilibrium, + &*candidate.shell.mitc4_drilling); + if (!ValidDrillingEvidence(candidate.shell.mitc4_drilling, + candidate.shell.equilibrium)) + return RecoveryFailure("invalid-shell-drilling-evidence", + {domain.SourcePath(), 0U}, "global-equilibrium", + "Shell drilling diagnostics must be finite and consistent."); + if (candidate.shell.mitc4_drilling->mode == Mitc4DrillingMode::kDirectSpring) + acceptance_moment = moment_balance - + Vector3(candidate.shell.mitc4_drilling->drilling_moment); + } const double force_metric = NormalizedBalance(force_balance.Components(), (std::max)(applied_force_scale, reaction_force_scale)); const double moment_metric = NormalizedBalance( - moment_balance.Components(), + acceptance_moment.Components(), (std::max)(applied_moment_scale, reaction_moment_scale)); candidate.shell.verification_metrics = {global.normalized_free_residual, force_metric, moment_metric}; diff --git a/src/fesa/results/mitc4_drilling_evidence.h b/src/fesa/results/mitc4_drilling_evidence.h new file mode 100644 index 0000000..29ed8b3 --- /dev/null +++ b/src/fesa/results/mitc4_drilling_evidence.h @@ -0,0 +1,54 @@ +#ifndef FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ +#define FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ + +#include +#include +#include +#include + +#include "fesa/results/result_records.h" + +namespace fesa::results_internal { + +inline bool ValidDrillingMode(const Mitc4DrillingMode mode) { + return mode == Mitc4DrillingMode::kDirectSpring || + mode == Mitc4DrillingMode::kLegacyCorotatedTwist; +} + +/// @brief Validates exact raw-minus-internal diagnostics without a tolerance. +inline bool ValidDrillingEvidence( + const std::optional& evidence, + const std::array& raw) { + if (!evidence) return true; // Older manually constructed fixtures. + if (!ValidDrillingMode(evidence->mode)) return false; + for (std::size_t i = 0U; i < 6U; ++i) { + const double expected = i < 3U ? raw[i] + : raw[i] - evidence->drilling_moment[i - 3U]; + if (!std::isfinite(raw[i]) || !std::isfinite(expected) || + !std::isfinite(evidence->equilibrium_with_drilling[i]) || + evidence->equilibrium_with_drilling[i] != expected) + return false; + } + for (const double value : evidence->drilling_moment) + if (!std::isfinite(value)) return false; + return true; +} + +inline bool SameDrillingInventory( + const std::optional& left, + const std::optional& right) { + return left.has_value() == right.has_value() && + (!left || left->mode == right->mode); +} + +/// @brief Keeps force balance raw and subtracts positive artificial moment once. +inline void CompleteDrillingEvidence( + const std::array& raw, Mitc4DrillingEvidence* evidence) { + evidence->equilibrium_with_drilling = raw; + for (std::size_t i = 0U; i < 3U; ++i) + evidence->equilibrium_with_drilling[3U + i] -= evidence->drilling_moment[i]; +} + +} // namespace fesa::results_internal + +#endif // FESA_RESULTS_MITC4_DRILLING_EVIDENCE_H_ diff --git a/src/fesa/results/shell_result_recovery.cpp b/src/fesa/results/shell_result_recovery.cpp index b7c00da..b443b0f 100644 --- a/src/fesa/results/shell_result_recovery.cpp +++ b/src/fesa/results/shell_result_recovery.cpp @@ -4,6 +4,8 @@ #include #include +#include "results/mitc4_drilling_evidence.h" + namespace fesa::results_internal { bool IsFiniteShellRow(const ShellResultRow& row) { @@ -29,6 +31,27 @@ Status ShellResultRecovery(const SourceLocation& location, const EntityIndex expected_element, const ShellElementResultRows& rows, RecoveryCandidate& candidate) { + if ((rows.drilling_mode && !ValidDrillingMode(*rows.drilling_mode)) || + !IsFiniteArray(rows.drilling_moment) || + (candidate.has_shell_results && + (candidate.shell.mitc4_drilling.has_value() != rows.drilling_mode.has_value() || + (rows.drilling_mode && + candidate.shell.mitc4_drilling->mode != *rows.drilling_mode)))) { + return RecoveryFailure("invalid-shell-drilling-evidence", location, + source_id.source_label_text, + "Shell recovery requires one finite actual drilling mode."); + } + if (rows.drilling_mode) { + if (!candidate.shell.mitc4_drilling) + candidate.shell.mitc4_drilling = Mitc4DrillingEvidence{*rows.drilling_mode}; + for (std::size_t i = 0U; i < 3U; ++i) { + candidate.shell.mitc4_drilling->drilling_moment[i] += rows.drilling_moment[i]; + if (!std::isfinite(candidate.shell.mitc4_drilling->drilling_moment[i])) + return RecoveryFailure("nonfinite-recovery-value", location, + source_id.source_label_text, + "Source-ordered drilling moment overflowed."); + } + } const double accumulated_energy = candidate.shell.physical_strain_energy + rows.physical_strain_energy; if (!std::isfinite(rows.physical_strain_energy) || diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9c4a638..1c81290 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -162,6 +162,8 @@ add_executable( reference/thermal_reference_cases_test.cpp ) +target_include_directories(fesa_reference_tests PRIVATE ${PROJECT_SOURCE_DIR}/src/fesa) + target_link_libraries( fesa_reference_tests PRIVATE diff --git a/tests/integration/analysis/linear_frequency_procedure_test.cpp b/tests/integration/analysis/linear_frequency_procedure_test.cpp index cd15396..fab87a8 100644 --- a/tests/integration/analysis/linear_frequency_procedure_test.cpp +++ b/tests/integration/analysis/linear_frequency_procedure_test.cpp @@ -18,6 +18,9 @@ #include "fesa/analysis/prepared_analysis.h" #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" +#include "fesa/elements/element_factory.h" +#include "fesa/results/result_recovery.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/fem/dof_manager.h" #include "fesa/model/domain.h" #include "fesa/model/model_types.h" @@ -132,6 +135,8 @@ TEST(LinearFrequencyProcedureExecutor, fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value()); ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); const auto& state = result.Value(); + ASSERT_TRUE(state.drilling_mode); + EXPECT_EQ(*state.drilling_mode, fesa::Mitc4DrillingMode::kDirectSpring); ASSERT_EQ(state.modes.size(), 1U); EXPECT_EQ(state.mass_null_dimension, 0U); EXPECT_EQ(state.finite_dimension, 1U); @@ -189,10 +194,10 @@ TEST(LinearFrequencyProcedureExecutor, } } -// E2-C09 is a pencil oracle. The public static contract conflict is recorded -// separately; direct covariance does not certify spatial moment equilibrium. +// E2-C09 keeps raw spatial imbalance visible while the direct-spring gate +// subtracts only the separately recovered artificial internal moment. TEST(LinearFrequencyProcedureExecutor, - PartialBoundaryStaticRejectsUnbalancedDrillingMoment) { + PartialBoundaryStaticSeparatesArtificialMomentAndPreservesLegacyGate) { const auto domain = PanelDomain(false, false, true); ASSERT_TRUE(domain.HasValue()); const fesa::SerialParallelFor serial; @@ -200,9 +205,32 @@ TEST(LinearFrequencyProcedureExecutor, const auto result = fesa::LinearStaticProcedureExecutor{serial, solver}.Execute( domain.Value()); - ASSERT_FALSE(result.HasValue()); - EXPECT_NE(StatusText(result.GetStatus()) - .find("global-equilibrium-tolerance-failure"), + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + const auto& state = result.Value().state; + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_NE(state.Equilibrium()[5U], 0.0); + for (std::size_t c = 0U; c < 3U; ++c) { + EXPECT_NEAR(state.Equilibrium()[3U + c], state.Mitc4Drilling()->drilling_moment[c], 1e-14); + EXPECT_NEAR(state.Mitc4Drilling()->equilibrium_with_drilling[3U + c], 0.0, 1e-14); + } + auto model = fesa::AnalysisModel::Create(domain.Value()); + ASSERT_TRUE(model.HasValue()); + auto dofs = fesa::DofManager::Create(model.Value()); + ASSERT_TRUE(dofs.HasValue()); + auto element = fesa::ElementFactory{}.Create(domain.Value().Elements()[0U], domain.Value()); + ASSERT_TRUE(element.HasValue()); + auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + *shell = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell, fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + fesa::ElementView view{std::cref(*shell)}; + auto stiffness = fesa::SparseAssembler::Assemble(view, dofs.Value(), serial); + ASSERT_TRUE(stiffness.HasValue()); + auto legacy_state = fesa::AnalysisState::Create(dofs.Value(), {"Step-1", 0U}); + const auto legacy = fesa::ResultRecovery::Recover(model.Value(), view, dofs.Value(), + stiffness.Value(), state.Displacement(), state.ExternalForce(), legacy_state); + ASSERT_FALSE(legacy.IsOk()); + EXPECT_NE(StatusText(legacy).find("global-equilibrium-tolerance-failure"), std::string::npos); } diff --git a/tests/integration/analysis/nonlinear_static_procedure_test.cpp b/tests/integration/analysis/nonlinear_static_procedure_test.cpp index 8f41c34..cf4bcd6 100644 --- a/tests/integration/analysis/nonlinear_static_procedure_test.cpp +++ b/tests/integration/analysis/nonlinear_static_procedure_test.cpp @@ -16,6 +16,13 @@ #include #include "fesa/elements/finite_rotation.h" +#include "fesa/analysis/prepared_analysis.h" +#include "fesa/analysis/nonlinear_static_procedure.h" +#include "fesa/analysis/newton_iteration_log.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/mitc4_shell.h" +#include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" #include "fesa/math/vector3.h" namespace { @@ -343,3 +350,57 @@ TEST(NonlinearStaticProcedure, } } // namespace + +TEST(Mitc4DrillingEvidence, NonlinearAndPressureFramesRetainSignedTorqueOracle) { + for (const bool pressure : {false, true}) { + SCOPED_TRACE(pressure); + TempDirectory directory{"drilling"}; + auto deck = AllConstrainedNonlinearShellDeck(); + const auto boundary = deck.find("All, 1, 6"); + ASSERT_NE(boundary, std::string::npos); + deck.replace(boundary, std::string{"All, 1, 6"}.size(), "All, 1, 5\nAll, 6, 6, 0.125"); + if (pressure) { + const auto assembly_end = deck.find("*End Assembly"); + deck.insert(assembly_end, "*Surface, TYPE=ELEMENT, NAME=Pressure\nShell-1.1, SPOS\n"); + deck.insert(deck.find("*End Step"), "*Dsload\nPressure, P, 4.0\n"); + } + const auto path = directory.Path() / "drilling.inp"; + WriteText(path, deck); + auto prepared = fesa::PrepareAnalysisInput(path); + ASSERT_TRUE(prepared.HasValue()); + const auto& domain = prepared.Value().domain; + auto element = fesa::ElementFactory{}.Create(domain.Elements()[0U], domain); + ASSERT_TRUE(element.HasValue()); + const auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + const auto stiffness = shell->Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double trace = 0.0; + for (std::size_t node = 0U; node < 4U; ++node) + for (std::size_t c = 3U; c < 5U; ++c) + trace += stiffness.Value().physical_local20(5U * node + c, 5U * node + c); + const double kd = 1e-3 * trace / 8.0; + fesa::SerialParallelFor serial; + fesa::MklPardisoTangentSolver tangent; + fesa::NewtonIterationLogWriter log; + ASSERT_TRUE(log.Open(directory.Path() / "drilling.log").IsOk()); + auto result = fesa::NonlinearStaticProcedureExecutor{serial, tangent}.Execute(domain, log); + ASSERT_TRUE(log.Close().IsOk()); + ASSERT_TRUE(result.HasValue()) << (result.GetStatus().Diagnostics().empty() ? "" : result.GetStatus().Diagnostics()[0U].message); + ASSERT_EQ(result.Value().state.AcceptedFrames().size(), 10U); + for (const auto& frame : result.Value().state.AcceptedFrames()) { + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_EQ(frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + const double torque = 4.0 * kd * 0.125 * frame.load_factor; + EXPECT_NEAR(frame.global_equilibrium[5U], torque, 1e-13); + EXPECT_NEAR(frame.mitc4_drilling->drilling_moment[2U], torque, 1e-13); + for (const double value : frame.mitc4_drilling->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-13); + EXPECT_NEAR(frame.physical_strain_energy, 0.0, 1e-20); + if (pressure) EXPECT_NE(frame.external_force[2U], 0.0); + } + ASSERT_TRUE(result.Value().state.CommittedState().mitc4_drilling); + EXPECT_EQ(result.Value().state.CommittedState().mitc4_drilling->drilling_moment, + result.Value().state.AcceptedFrames().back().mitc4_drilling->drilling_moment); + } +} diff --git a/tests/integration/app/fesa_application_test.cpp b/tests/integration/app/fesa_application_test.cpp index efdbb4f..04942be 100644 --- a/tests/integration/app/fesa_application_test.cpp +++ b/tests/integration/app/fesa_application_test.cpp @@ -17,6 +17,8 @@ #include #include +#include "io/hdf5/mitc4_drilling_diagnostics.h" + namespace { constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; @@ -525,4 +527,40 @@ TEST(Mitc4ShellCli, WritesShellHdf5ThroughTheExistingApplicationRoute) { const auto result = RunApplication({input.string()}); ASSERT_EQ(result.exit_code, 0) << result.standard_error; ExpectShellHdf5Identity(output, input); + const auto file = OpenFile(output); + const auto drilling = fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()); + ASSERT_TRUE(drilling); + EXPECT_EQ(drilling->formulation, "direct-initial-director-spring"); + EXPECT_EQ(drilling->version, 1U); + ASSERT_EQ(drilling->frames.size(), 1U); + EXPECT_EQ(drilling->frames.front().frame_index, 0U); +} + + +TEST(Mitc4ShellCli, DrillingActualHhtPipelineWritesInitialAndAcceptedDiagnostics) { + TempDirectory directory{"shell-drilling-hht"}; + const auto input = directory.Path() / "shell.inp"; + const auto output = directory.Path() / "shell_results.h5"; + std::string deck = AllConstrainedShellDeck(); + const auto material = deck.find("*Material, name=Steel\n"); + ASSERT_NE(material, std::string::npos); + deck.insert(material + std::string{"*Material, name=Steel\n"}.size(), "*Density\n2.\n"); + const auto step = deck.find("*Step,"); + ASSERT_NE(step, std::string::npos); + deck.replace(step, std::string::npos, + "*Step, name=Load, nlgeom=YES\n*Dynamic, Direct\n0.25, 0.5\n*Cload\nAll, 3, 1.\n*End Step\n"); + WriteText(input, deck); + const auto result = RunApplication({input.string()}); + ASSERT_EQ(result.exit_code, 0) << result.standard_error; + const auto file = OpenFile(output); + const auto drilling = fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()); + ASSERT_TRUE(drilling); + EXPECT_EQ(drilling->formulation, "direct-initial-director-spring"); + ASSERT_EQ(drilling->frames.size(), 3U); + for (std::size_t i = 0; i < 3U; ++i) { + EXPECT_EQ(drilling->frames[i].frame_index, i); + EXPECT_EQ(drilling->frames[i].raw_equilibrium, (std::array{})); + EXPECT_EQ(drilling->frames[i].drilling_moment, (std::array{})); + EXPECT_EQ(drilling->frames[i].equilibrium_with_drilling, (std::array{})); + } } diff --git a/tests/reference/eigen_reference_comparison.cpp b/tests/reference/eigen_reference_comparison.cpp index 6a80dd5..d0920dd 100644 --- a/tests/reference/eigen_reference_comparison.cpp +++ b/tests/reference/eigen_reference_comparison.cpp @@ -20,6 +20,8 @@ #include "fesa/io/abaqus/domain_mapper.h" #include "fesa/io/abaqus/input_reader.h" +#include "mitc4_drilling_reader.h" + namespace fesa::test { namespace { @@ -385,6 +387,9 @@ Result EigenReferenceComparison::Compare( } EigenReferenceComparisonReport report{}; + report.drilling_diagnostics = ReadMitc4DrillingReport(file.Get()); + if (report.drilling_diagnostics && !report.drilling_diagnostics->frames.empty()) + Fail("schema-mismatch", "Eigen drilling diagnostics must have no equilibrium frames."); report.case_id = reference_case.case_id; report.source_step_name = frequency.source_step_name_raw; report.requested_mode_count = requested; diff --git a/tests/reference/eigen_reference_comparison.h b/tests/reference/eigen_reference_comparison.h index cb6c002..35e53b4 100644 --- a/tests/reference/eigen_reference_comparison.h +++ b/tests/reference/eigen_reference_comparison.h @@ -7,6 +7,7 @@ #include #include "fesa/core/status.h" +#include "mitc4_drilling_report.h" #include "reference_tolerance_policy.h" namespace fesa::test { @@ -36,6 +37,7 @@ struct EigenReferenceComparisonReport { bool passed; std::string policy_id; double relative_rms_tolerance; + std::optional drilling_diagnostics; }; /// @brief Compares modal eigenvalues without changing the approved input or diff --git a/tests/reference/eigen_reference_comparison_test.cpp b/tests/reference/eigen_reference_comparison_test.cpp index 8b3dc6b..3e0a5ff 100644 --- a/tests/reference/eigen_reference_comparison_test.cpp +++ b/tests/reference/eigen_reference_comparison_test.cpp @@ -123,6 +123,25 @@ TEST(EigenReferenceComparison, EXPECT_EQ(ok.Value().policy_id, "mitc4-eigen-first-implementation"); EXPECT_DOUBLE_EQ(ok.Value().relative_rms_tolerance, 0.05); EXPECT_EQ(ok.Value().families.size(), 2U); + EXPECT_FALSE(ok.Value().drilling_diagnostics); + for (const std::string mode : {"direct-initial-director-spring", "legacy-corotated-twist"}) { + f = H5Fopen(h5.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT); + m = H5Gopen2(f, "/metadata", H5P_DEFAULT); + if (H5Aexists(m, "mitc4_drilling_formulation") > 0) H5Adelete(m, "mitc4_drilling_formulation"); + if (H5Aexists(m, "drilling_diagnostics_version") > 0) H5Adelete(m, "drilling_diagnostics_version"); + Str(m, "mitc4_drilling_formulation", mode); + U(m, "drilling_diagnostics_version", 1U); + H5Gclose(m); + H5Fclose(f); + const auto compared = fesa::test::EigenReferenceComparison::Compare(reference_case); + ASSERT_TRUE(compared.HasValue()); + ASSERT_TRUE(compared.Value().drilling_diagnostics); + EXPECT_EQ(compared.Value().drilling_diagnostics->formulation, mode); + EXPECT_EQ(compared.Value().drilling_diagnostics->version, 1U); + EXPECT_TRUE(compared.Value().drilling_diagnostics->frames.empty()); + EXPECT_EQ(compared.Value().families.size(), 2U); + } + // Only synthetic FESA output is changed; each TXT mode remains an independent // scalar family, including at both inclusive five-percent boundaries. for (const auto& sample : diff --git a/tests/reference/mitc4_drilling_reader.h b/tests/reference/mitc4_drilling_reader.h new file mode 100644 index 0000000..878f9ac --- /dev/null +++ b/tests/reference/mitc4_drilling_reader.h @@ -0,0 +1,21 @@ +#ifndef FESA_TESTS_REFERENCE_MITC4_DRILLING_READER_H_ +#define FESA_TESTS_REFERENCE_MITC4_DRILLING_READER_H_ + +#include "mitc4_drilling_report.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +namespace fesa::test { + +inline std::optional ReadMitc4DrillingReport(hid_t file) { + const auto stored = hdf5_internal::ReadDrillingDiagnostics(file); + if (!stored) return std::nullopt; + Mitc4DrillingDiagnosticsReport result{stored->formulation, stored->version, {}}; + for (const auto& frame : stored->frames) + result.frames.push_back({frame.storage_step, frame.frame_index, frame.raw_equilibrium, + frame.drilling_moment, frame.equilibrium_with_drilling}); + return result; +} + +} // namespace fesa::test + +#endif // FESA_TESTS_REFERENCE_MITC4_DRILLING_READER_H_ diff --git a/tests/reference/mitc4_drilling_report.h b/tests/reference/mitc4_drilling_report.h new file mode 100644 index 0000000..660a518 --- /dev/null +++ b/tests/reference/mitc4_drilling_report.h @@ -0,0 +1,58 @@ +#ifndef FESA_TESTS_REFERENCE_MITC4_DRILLING_REPORT_H_ +#define FESA_TESTS_REFERENCE_MITC4_DRILLING_REPORT_H_ + +#include +#include +#include +#include +#include +#include +#include + +namespace fesa::test { + +struct Mitc4DrillingFrameReport { + std::string storage_step; + std::uint64_t frame_index; + std::array raw_equilibrium; + std::array drilling_moment; + std::array equilibrium_with_drilling; +}; + +struct Mitc4DrillingDiagnosticsReport { + std::string formulation; + std::uint64_t version; + std::vector frames; +}; + +inline void WriteDrillingDiagnosticsJson( + std::ostream& stream, const std::optional& report) { + if (!report) { stream << "null"; return; } + stream << "{\"formulation\":\"" << report->formulation + << "\",\"version\":" << report->version << ",\"frames\":["; + const auto numbers = [&](const auto& values) { + stream << '['; + for (std::size_t i = 0; i < values.size(); ++i) { + if (i != 0U) stream << ','; + stream << values[i]; + } + stream << ']'; + }; + for (std::size_t i = 0; i < report->frames.size(); ++i) { + if (i != 0U) stream << ','; + const auto& frame = report->frames[i]; + stream << "{\"storage_step\":\"Step-1\",\"frame_index\":" << frame.frame_index + << ",\"raw_equilibrium\":"; + numbers(frame.raw_equilibrium); + stream << ",\"drilling_moment\":"; + numbers(frame.drilling_moment); + stream << ",\"equilibrium_with_drilling\":"; + numbers(frame.equilibrium_with_drilling); + stream << '}'; + } + stream << "]}"; +} + +} // namespace fesa::test + +#endif // FESA_TESTS_REFERENCE_MITC4_DRILLING_REPORT_H_ diff --git a/tests/reference/mitc4_reference_comparison.cpp b/tests/reference/mitc4_reference_comparison.cpp index bc8100c..f90181a 100644 --- a/tests/reference/mitc4_reference_comparison.cpp +++ b/tests/reference/mitc4_reference_comparison.cpp @@ -26,6 +26,8 @@ #include "reference_tolerance_policy.h" +#include "mitc4_drilling_reader.h" + namespace fesa::test { namespace { @@ -573,6 +575,7 @@ std::vector ReadDisplacement(const hid_t file, } struct Hdf5Projection { + std::optional drilling_diagnostics; std::vector rows; std::string source_element_type; std::string internal_formulation; @@ -619,7 +622,7 @@ Hdf5Projection ReadHdf5(const Mitc4ReferenceCase& reference_case) { values.begin() + static_cast(row * kComponents.size()), kComponents.size(), rows[row].values.begin()); } - return {std::move(rows), reference_case.expected_source_element_type, + return {ReadMitc4DrillingReport(file.Get()), std::move(rows), reference_case.expected_source_element_type, kInternalFormulation, kIntegrationRule}; } @@ -739,6 +742,7 @@ Result Mitc4ReferenceComparison::Compare( } Mitc4ComparisonReport report{}; + report.drilling_diagnostics = std::move(hdf5.drilling_diagnostics); report.case_id = reference_case.case_id; report.source_element_type = std::move(hdf5.source_element_type); report.internal_formulation = std::move(hdf5.internal_formulation); @@ -889,7 +893,9 @@ Status Mitc4ReferenceComparison::WriteDeterministicJson( } stream.imbue(std::locale::classic()); stream << std::setprecision(std::numeric_limits::max_digits10); - stream << "{\"case_id\":\"" << JsonEscape(report.case_id) + stream << "{\"drilling_diagnostics\":"; + WriteDrillingDiagnosticsJson(stream, report.drilling_diagnostics); + stream << ",\"case_id\":\"" << JsonEscape(report.case_id) << "\",\"source_element_type\":\"" << JsonEscape(report.source_element_type) << "\",\"internal_formulation\":\"" diff --git a/tests/reference/mitc4_reference_comparison.h b/tests/reference/mitc4_reference_comparison.h index 2fb07b8..e3753e8 100644 --- a/tests/reference/mitc4_reference_comparison.h +++ b/tests/reference/mitc4_reference_comparison.h @@ -8,6 +8,7 @@ #include #include "fesa/core/status.h" +#include "mitc4_drilling_report.h" namespace fesa::test { @@ -77,6 +78,7 @@ struct Mitc4ComparisonReport { std::vector warnings; std::size_t worst_row; bool passed; + std::optional drilling_diagnostics; }; class Mitc4ReferenceComparison { diff --git a/tests/reference/mitc4_reference_comparison_test.cpp b/tests/reference/mitc4_reference_comparison_test.cpp index a65397f..0a35450 100644 --- a/tests/reference/mitc4_reference_comparison_test.cpp +++ b/tests/reference/mitc4_reference_comparison_test.cpp @@ -22,6 +22,8 @@ #error FESA_TEST_BINARY_DIR must identify the CMake binary root. #endif +#include "io/hdf5/mitc4_drilling_diagnostics.h" + namespace { constexpr const char* kInstanceName = "PART-1-1"; @@ -771,3 +773,56 @@ TEST(Mitc4ReferenceComparison, RequiresOnlyDeclaredInputCsvAndHdf5) { } } // namespace + + +TEST(Mitc4ReferenceComparison, DrillingOldMarkerlessReportIsNull) { + ContractFixture fixture{"drilling-old"}; + const auto compared = fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()); + ASSERT_TRUE(compared.HasValue()); + EXPECT_FALSE(compared.Value().drilling_diagnostics); + const auto json = fixture.Root() / "drilling-old.json"; + ASSERT_TRUE(fesa::test::Mitc4ReferenceComparison::WriteDeterministicJson(compared.Value(), json).IsOk()); + EXPECT_NE(ReadBytes(json).find("\"drilling_diagnostics\":null"), std::string::npos); +} + +TEST(Mitc4ReferenceComparison, DrillingNewMarkersPreserveAllFramesWithoutNewVerdict) { + for (const auto mode : {fesa::Mitc4DrillingMode::kDirectSpring, + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist}) { + ContractFixture fixture{"drilling-new"}; + { + Hdf5Handle file{H5Fopen(fixture.Results().string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::WriteDrillingMetadata(file.Get(), mode); + // Nonzero raw balance is diagnostic only for the linear displacement comparator. + for (const int index : {0, 2, 10}) { + const std::string root = "/steps/Step-1/frames/" + std::to_string(index); + (void)fesa::hdf5_internal::CreateGroup(file.Get(), (root + "/global").c_str()); + const std::array raw{1., 2., 3., 4., 5., 6.}; + const fesa::Mitc4DrillingEvidence evidence{mode, {4., 5., 6.}, {1., 2., 3., 0., 0., 0.}}; + fesa::hdf5_internal::WriteModelDoubleDataset(file.Get(), root + "/global/equilibrium", {6U}, + raw.data(), 6U, "", "", "", ""); + fesa::hdf5_internal::WriteDrillingFrame(file.Get(), root, evidence); + } + } + const auto compared = fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()); + ASSERT_TRUE(compared.HasValue()); + EXPECT_TRUE(compared.Value().passed); + EXPECT_EQ(compared.Value().metrics.size(), 2U); + const auto& diagnostics = compared.Value().drilling_diagnostics; + ASSERT_TRUE(diagnostics); + EXPECT_EQ(diagnostics->formulation, fesa::Mitc4DrillingModeName(mode)); + ASSERT_EQ(diagnostics->frames.size(), 3U); + EXPECT_EQ(diagnostics->frames[0].frame_index, 0U); + EXPECT_EQ(diagnostics->frames[1].frame_index, 2U); + EXPECT_EQ(diagnostics->frames[2].frame_index, 10U); + const auto json = fixture.Root() / "drilling-new.json"; + ASSERT_TRUE(fesa::test::Mitc4ReferenceComparison::WriteDeterministicJson(compared.Value(), json).IsOk()); + const auto bytes = ReadBytes(json); + EXPECT_NE(bytes.find("\"drilling_diagnostics\":{\"formulation\":\"" + + diagnostics->formulation + "\",\"version\":1,\"frames\":["), std::string::npos); + { + Hdf5Handle file{H5Fopen(fixture.Results().string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/2/global/drilling_moment", H5P_DEFAULT), 0); + } + EXPECT_FALSE(fesa::test::Mitc4ReferenceComparison::Compare(fixture.ReferenceCase()).HasValue()); + } +} diff --git a/tests/reference/nonlinear_static_reference_cases_test.cpp b/tests/reference/nonlinear_static_reference_cases_test.cpp index 0bf9c0d..04e7668 100644 --- a/tests/reference/nonlinear_static_reference_cases_test.cpp +++ b/tests/reference/nonlinear_static_reference_cases_test.cpp @@ -21,6 +21,8 @@ #error FESA_TEST_BINARY_DIR must identify the CMake binary root. #endif +#include "io/hdf5/mitc4_drilling_diagnostics.h" + namespace { std::string ReadBytes(const std::filesystem::path& path) { @@ -409,3 +411,75 @@ TEST(NonlinearStaticReference, EXPECT_FALSE(nonzero.normalized_metric.has_value()) << identity; } } + + +TEST(NonlinearStaticReference, DrillingSyntheticMomentSelectsOnlyExistingInternalFamily) { + const auto reference_case = MakeCase("shellnl-s4-final"); + const auto copied = fesa::test::CopyInputToExecutionDirectory(reference_case); + ASSERT_TRUE(copied.HasValue()); + const auto generated = fesa::test::RunFesaApplicationOnCopiedInput(reference_case, copied.Value()); + ASSERT_TRUE(generated.HasValue()); + const auto& path = generated.Value().generated_results_hdf5; + const auto baseline = fesa::test::CompareHdf5ToReferenceCsv(reference_case, path); + ASSERT_TRUE(baseline.HasValue()); + ASSERT_TRUE(baseline.Value().drilling_diagnostics); + const auto final_index = baseline.Value().drilling_diagnostics->frames.back().frame_index; + const std::string root = "/steps/Step-1/frames/" + std::to_string(final_index) + "/global/"; + const auto moment_gate = [](const auto& report) -> const auto& { + const auto found = std::find_if(report.internal_equilibrium.begin(), report.internal_equilibrium.end(), + [](const auto& family) { return family.family_identity == "coarse.global_moment_balance"; }); + if (found == report.internal_equilibrium.end()) throw std::runtime_error{"Missing existing moment family."}; + return *found; + }; + // This mutates only a generated working output into an algebra fixture. + // Its comparator verdict is not physics reference evidence. + for (const std::string mode : {"direct-initial-director-spring", "legacy-corotated-twist", "markerless"}) { + { + fesa::hdf5_internal::Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + if (mode != "markerless") { + fesa::hdf5_internal::WriteStringAttribute(metadata.Get(), "mitc4_drilling_formulation", mode); + fesa::hdf5_internal::WriteUint64Attribute(metadata.Get(), "drilling_diagnostics_version", 1U); + } + for (const auto& entry : std::vector>>{ + {"equilibrium", {0.,0.,0.,1e6,0.,0.}}, {"drilling_moment", {1e6,0.,0.}}, + {"equilibrium_with_drilling", {0.,0.,0.,0.,0.,0.}}}) { + fesa::hdf5_internal::Hdf5Handle dataset{H5Dopen2(file.Get(), (root + entry.first).c_str(), H5P_DEFAULT), H5Dclose}; + ASSERT_GE(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, entry.second.data()), 0); + } + } + const auto compared = fesa::test::CompareHdf5ToReferenceCsv(reference_case, path); + ASSERT_TRUE(compared.HasValue()); + EXPECT_TRUE(compared.Value().external_blocking_passed); + EXPECT_EQ(compared.Value().families.size(), baseline.Value().families.size()); + EXPECT_EQ(compared.Value().internal_equilibrium.size(), baseline.Value().internal_equilibrium.size()); + const auto& gate = moment_gate(compared.Value()); + EXPECT_EQ(gate.scale, moment_gate(baseline.Value()).scale); + EXPECT_EQ(gate.numerator, mode == "direct-initial-director-spring" ? 0.0 : 1e6); + EXPECT_EQ(gate.passed, mode == "direct-initial-director-spring"); + EXPECT_EQ(compared.Value().drilling_diagnostics.has_value(), mode != "markerless"); + } + // Restore new metadata, retain exact raw-minus-M algebra, and expose a real imbalance. + { + fesa::hdf5_internal::Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::WriteDrillingMetadata(file.Get(), fesa::Mitc4DrillingMode::kDirectSpring); + for (const auto& entry : std::vector>>{ + {"equilibrium", {0.,0.,0.,2e6,0.,0.}}, {"equilibrium_with_drilling", {0.,0.,0.,1e6,0.,0.}}}) { + fesa::hdf5_internal::Hdf5Handle dataset{H5Dopen2(file.Get(), (root + entry.first).c_str(), H5P_DEFAULT), H5Dclose}; + ASSERT_GE(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, + H5P_DEFAULT, entry.second.data()), 0); + } + } + const auto imbalanced = fesa::test::CompareHdf5ToReferenceCsv(reference_case, path); + ASSERT_TRUE(imbalanced.HasValue()); + EXPECT_FALSE(moment_gate(imbalanced.Value()).passed); + EXPECT_EQ(moment_gate(imbalanced.Value()).scale, moment_gate(baseline.Value()).scale); + { + fesa::hdf5_internal::Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/0/global/drilling_moment", H5P_DEFAULT), 0); + } + EXPECT_FALSE(fesa::test::CompareHdf5ToReferenceCsv(reference_case, path).HasValue()); +} diff --git a/tests/reference/nonlinear_static_reference_comparison.cpp b/tests/reference/nonlinear_static_reference_comparison.cpp index 874e701..b55dae6 100644 --- a/tests/reference/nonlinear_static_reference_comparison.cpp +++ b/tests/reference/nonlinear_static_reference_comparison.cpp @@ -37,6 +37,8 @@ #error FESA_TEST_BINARY_DIR must identify the CMake binary root. #endif +#include "mitc4_drilling_reader.h" + namespace fesa::test { namespace { @@ -896,6 +898,8 @@ std::vector ReadConstraintMask(const hid_t file, } struct Hdf5Projection { + std::uint64_t selected_frame_index; + std::optional drilling_diagnostics; std::filesystem::path copied_input_path; std::vector nodes; std::vector displacement; @@ -972,7 +976,7 @@ Hdf5Projection ReadNonlinearHdf5( const std::uint64_t frame_index = ReadFactorOneFrame(file.Get()); const std::string root = "/steps/Step-1/frames/" + std::to_string(frame_index) + "/nodal/"; - return {copied_input, + return {frame_index, ReadMitc4DrillingReport(file.Get()), copied_input, nodes, ReadNodalDataset(file.Get(), root + "displacement", "UX,UY,UZ,URX,URY,URZ", contract.node_count), @@ -1320,6 +1324,16 @@ void AppendInternalEquilibrium(NonlinearReferenceComparisonReport* report, applied_moment_scale += applied_resultant.Norm(); reaction_moment_scale += constrained_resultant.Norm(); } + const std::array* moment_balance = &projection.global_equilibrium; + if (projection.drilling_diagnostics && + projection.drilling_diagnostics->formulation == "direct-initial-director-spring") { + const auto& frames = projection.drilling_diagnostics->frames; + const auto selected = std::find_if(frames.begin(), frames.end(), [&](const auto& frame) { + return frame.frame_index == projection.selected_frame_index; + }); + if (selected == frames.end()) Fail("schema-mismatch", "Missing factor-one drilling diagnostics."); + moment_balance = &selected->equilibrium_with_drilling; + } report->constrained_reaction_resultant = { reaction_force.Components(), reaction_moment.Components()}; const std::string prefix = ContractFor(reference_case).family_prefix; @@ -1329,9 +1343,9 @@ void AppendInternalEquilibrium(NonlinearReferenceComparisonReport* report, projection.global_equilibrium[2U]}, applied_force_scale, reaction_force_scale); AppendGlobalBalance(report, prefix + "global_moment_balance", - {projection.global_equilibrium[3U], - projection.global_equilibrium[4U], - projection.global_equilibrium[5U]}, + {(*moment_balance)[3U], + (*moment_balance)[4U], + (*moment_balance)[5U]}, applied_moment_scale, reaction_moment_scale); AppendFreeResidual(report, prefix + "accepted_free_force_residual", projection, false); @@ -1385,6 +1399,7 @@ void CompareLinearCase(NonlinearReferenceComparisonReport& report, : diagnostics.front().message; Fail("schema-mismatch", message); } + report.drilling_diagnostics = legacy.Value().drilling_diagnostics; for (const auto& metric : legacy.Value().metrics) { ReferenceToleranceFamilyDecision tolerance{}; tolerance.identity = metric.blocking ? "linear.translation" : "linear.rotation"; @@ -1418,6 +1433,7 @@ void CompareNonlinearCase(NonlinearReferenceComparisonReport& report, report.generated_results_hdf5); report.copied_input_path = hdf5.copied_input_path; report.generated_log_path = DerivedLogPath(hdf5.copied_input_path); + report.drilling_diagnostics = hdf5.drilling_diagnostics; RequireExactNodeLabels(hdf5.nodes, contract.node_count); const auto displacement_rows = ReadCsvWideRows( reference_case.displacement_csv_path, true, kDisplacementComponents, @@ -1475,6 +1491,7 @@ void ComparePressureCase(NonlinearReferenceComparisonReport& report, ReadNonlinearHdf5(reference_case, report.generated_results_hdf5); report.copied_input_path = hdf5.copied_input_path; report.generated_log_path = DerivedLogPath(hdf5.copied_input_path); + report.drilling_diagnostics = hdf5.drilling_diagnostics; RequireExactNodeLabels(hdf5.nodes, contract.node_count); RequirePressureInstance(hdf5.nodes); const auto displacement_rows = diff --git a/tests/reference/nonlinear_static_reference_comparison.h b/tests/reference/nonlinear_static_reference_comparison.h index a0ac0aa..f5de95d 100644 --- a/tests/reference/nonlinear_static_reference_comparison.h +++ b/tests/reference/nonlinear_static_reference_comparison.h @@ -10,6 +10,7 @@ #include #include "fesa/core/status.h" +#include "mitc4_drilling_report.h" #include "reference_tolerance_policy.h" namespace fesa::test { @@ -114,6 +115,7 @@ struct NonlinearReferenceComparisonReport { ConstrainedReactionResultantReport constrained_reaction_resultant; bool external_blocking_passed; bool internal_blocking_passed; + std::optional drilling_diagnostics; }; enum class NonlinearReferenceClassification { diff --git a/tests/test_compare_hht_dynamic.py b/tests/test_compare_hht_dynamic.py index 3bbe9a7..6328712 100644 --- a/tests/test_compare_hht_dynamic.py +++ b/tests/test_compare_hht_dynamic.py @@ -349,3 +349,88 @@ def test_cli_writes_json_and_exit_two_for_missing_results(tmp_path, artifacts): evidence = json.loads(report.read_text()) assert not evidence["precheck"]["passed"] assert evidence["families"] == [] + + +DRILLING_MODES = ("direct-initial-director-spring", "legacy-corotated-twist") + + +def drilling_result(tmp_path, artifacts, mode=DRILLING_MODES[0]): + path = changed_result(tmp_path, artifacts) + with h5py.File(path, "r+") as file: + file["metadata"].attrs["mitc4_drilling_formulation"] = mode + file["metadata"].attrs["drilling_diagnostics_version"] = np.uint64(1) + for frame in range(101): + root = file.require_group(f"steps/Step-1/frames/{frame}/global") + raw = np.array([1., 2., 3., 4., 5., 6.], dtype="mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) next.mitc4_drilling->mode = static_cast(99); + if (kind == 3) next.mitc4_drilling->equilibrium_with_drilling[3U] = 1.0; + if (kind == 4) next.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::quiet_NaN(); + EXPECT_FALSE(fixture.state.CommitAcceptedFrame(next).IsOk()) << kind; + ASSERT_EQ(fixture.state.AcceptedFrames().size(), 1U); + EXPECT_EQ(fixture.state.CommittedState().frame_index, 0U); + EXPECT_EQ(fixture.state.CommittedState().mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + } + auto valid = MakeFrame(*fixture.domain, 1U); + valid.mitc4_drilling = initial.mitc4_drilling; + EXPECT_TRUE(fixture.state.CommitAcceptedFrame(valid).IsOk()); +} diff --git a/tests/unit/analysis/hht_dynamic_algebra_test.cpp b/tests/unit/analysis/hht_dynamic_algebra_test.cpp index a6da082..06e1c87 100644 --- a/tests/unit/analysis/hht_dynamic_algebra_test.cpp +++ b/tests/unit/analysis/hht_dynamic_algebra_test.cpp @@ -28,7 +28,8 @@ struct Fixture { Fixture MakeFixture(const std::size_t free_nodes = 1U, const bool partial_rotation = false, const bool translations_free = true, - const double youngs_modulus = 1000.0) { + const double youngs_modulus = 1000.0, + const bool all_rotations_constrained = false) { fesa::ModelDefinition definition{}; definition.source_path = "hht-algebra.inp"; definition.nodes = { @@ -66,7 +67,10 @@ Fixture MakeFixture(const std::size_t free_nodes = 1U, step.boundaries.push_back( {name, 1, 6, 0.0, {definition.source_path, 10U}}); } else { - if (partial_rotation) + if (all_rotations_constrained) + step.boundaries.push_back( + {name, 4, 6, 0.0, {definition.source_path, 10U}}); + else if (partial_rotation) step.boundaries.push_back( {name, 4, 5, 0.0, {definition.source_path, 10U}}); if (!translations_free) @@ -705,3 +709,57 @@ TEST(HhtEndpointFields, MeasuresActualFullChartAcrossNonlinearRootAndPartialBoun } } // namespace + +TEST(HhtDynamicAlgebra, NonDirectionalAssemblyRetainsValidEmptyNullBoundary) { + for (const std::size_t free_nodes : {4U, 0U}) { + SCOPED_TRACE(free_nodes); + const auto fixture = MakeFixture(free_nodes, false, true, 1000.0, true); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, + fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + auto trial = EmptyTrial(fixture); + trial.external_force[0U] = 0.4; + fesa::Vector load_scale(24U); + load_scale[0U] = 0.4; + fesa::SerialParallelFor serial; + const auto parameters = fesa::MakeHhtParameters(-0.05); + ASSERT_TRUE(parameters.HasValue()); + const auto assembly = fesa::hht_internal::EvaluateAssembly( + *fixture.domain, fixture.dofs, shells.Value(), trial, load_scale, + 0.1, parameters.Value(), serial, false); + ASSERT_TRUE(assembly.HasValue()); + ASSERT_TRUE(assembly.Value().frame.mitc4_drilling); + EXPECT_EQ(assembly.Value().frame.mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(assembly.Value().hidden_constraint.Size(), 0U); + EXPECT_FALSE(assembly.Value().hidden_gradient.has_value()); + const auto nodes = fesa::hht_internal::NullNodes(fixture.dofs, + assembly.Value().rotations); + ASSERT_TRUE(nodes.empty()); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + const auto acceleration = fesa::hht_internal::InitializeAcceleration( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, positive, tangent); + EXPECT_TRUE(acceleration.HasValue()) + << (acceleration.GetStatus().Diagnostics().empty() ? "" : + acceleration.GetStatus().Diagnostics()[0U].message); + auto frame = assembly.Value().frame; + if (acceleration.HasValue()) { + frame.acceleration = acceleration.Value(); + // rho*t*area/4 = .2 per node; the .4 force gives acceleration 2. + EXPECT_NEAR(acceleration.Value()[0U], free_nodes == 0U ? 0.0 : 2.0, 1e-12); + for (const auto full : fixture.dofs.ConstrainedDofs()) + EXPECT_DOUBLE_EQ(acceleration.Value()[full], 0.0); + } + const auto complete = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, nodes, + positive, tangent, &frame); + EXPECT_TRUE(complete.IsOk()) + << (complete.Diagnostics().empty() ? "" : complete.Diagnostics()[0U].message); + if (complete.IsOk()) { + EXPECT_EQ(frame.hidden_constraint.Size(), 0U); + EXPECT_EQ(frame.hidden_velocity_defect.Size(), 0U); + EXPECT_EQ(frame.hidden_acceleration_defect.Size(), 0U); + } + } +} diff --git a/tests/unit/analysis/hht_dynamic_procedure_test.cpp b/tests/unit/analysis/hht_dynamic_procedure_test.cpp index 75594e9..4b9ca29 100644 --- a/tests/unit/analysis/hht_dynamic_procedure_test.cpp +++ b/tests/unit/analysis/hht_dynamic_procedure_test.cpp @@ -216,6 +216,12 @@ void ExpectSameDynamicHistory(const RigidTranslationRun& left, EXPECT_DOUBLE_EQ(lhs.drilling_energy, rhs.drilling_energy); EXPECT_DOUBLE_EQ(lhs.external_work, rhs.external_work); EXPECT_DOUBLE_EQ(lhs.energy_defect, rhs.energy_defect); + ASSERT_TRUE(lhs.mitc4_drilling); + ASSERT_TRUE(rhs.mitc4_drilling); + EXPECT_EQ(lhs.mitc4_drilling->mode, rhs.mitc4_drilling->mode); + EXPECT_EQ(lhs.mitc4_drilling->drilling_moment, rhs.mitc4_drilling->drilling_moment); + EXPECT_EQ(lhs.mitc4_drilling->equilibrium_with_drilling, + rhs.mitc4_drilling->equilibrium_with_drilling); } } @@ -357,6 +363,10 @@ TEST(HhtDynamicProcedure, ClampedStepLoadHasPhysicalReactionAndZeroMotion) { EXPECT_EQ(frame.time_integration_formulation, fesa::DynamicTimeIntegrationFormulation::kReducedCoordinateHht); EXPECT_EQ(frame.iteration_count, 0U); + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_EQ(frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(frame.mitc4_drilling->drilling_moment, (std::array{})); + EXPECT_EQ(frame.mitc4_drilling->equilibrium_with_drilling, frame.global_equilibrium); EXPECT_DOUBLE_EQ(frame.reaction[2U], -4.0); EXPECT_DOUBLE_EQ(frame.external_force[2U], 4.0); EXPECT_DOUBLE_EQ(frame.kinetic_energy, 0.0); diff --git a/tests/unit/analysis/hht_reduced_shell_test.cpp b/tests/unit/analysis/hht_reduced_shell_test.cpp index 752198a..10577c7 100644 --- a/tests/unit/analysis/hht_reduced_shell_test.cpp +++ b/tests/unit/analysis/hht_reduced_shell_test.cpp @@ -22,10 +22,13 @@ #include "analysis/hht_reduced_internal.h" #include "math/dense_blas_internal.h" +#include "fesa/elements/mitc4_shell.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/analysis/analysis_model.h" #include "fesa/analysis/hht_time_integration.h" #include "fesa/assembly/parallel_for.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" namespace { @@ -625,6 +628,10 @@ TEST(HhtReducedOperators, [](Assembly& a) { a.trial_constraints->shells[0U].configuration.identity += "x"; }, + [](Assembly& a) { + a.trial_constraints->shells[0U].configuration.drilling_mode = + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + }, [](Assembly& a) { a.trial_constraints->shells[0U].q[0U] = -0.0; }, [](Assembly& a) { a.trial_constraints->shells[0U].v[0U] = -0.0; }, [](Assembly& a) { @@ -935,11 +942,15 @@ TEST(HhtReducedOperators, NonfiniteMatchingTimeKeyMissPreservesFreshSuccess) { } } -TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { +void CheckShellConstraintScatter(const fesa::Mitc4DrillingMode mode) { const auto fixture = MakeFixture(false); const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); ASSERT_TRUE(shells.HasValue()); + SCOPED_TRACE(static_cast(mode)); + for (const auto& shell : shells.Value()) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, mode); fesa::DynamicAcceptedFrameState frame{}; for (auto* values : {&frame.generalized_coordinates, &frame.velocity, &frame.acceleration, &frame.algorithmic_acceleration, @@ -980,7 +991,10 @@ TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { assembly.Value().frame.residual_scale[6U * node + 3U + component]; EXPECT_DOUBLE_EQ(row.value_scale, expected_scale); EXPECT_FALSE(row.gradient.empty()); - EXPECT_FALSE(row.hessian.empty()); + if (mode == fesa::Mitc4DrillingMode::kLegacyCorotatedTwist) + EXPECT_FALSE(row.hessian.empty()); + else + EXPECT_TRUE(row.hessian.empty()); EXPECT_LT(row.gradient.size(), fixture.dofs.FullDofCount()); for (const auto& entry : row.hessian) { EXPECT_LT(entry.row, fixture.dofs.FullDofCount()); @@ -989,6 +1003,14 @@ TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { } } +TEST(HhtReducedOperators, ScatterShellDerivativesByConstraintRow) { + CheckShellConstraintScatter(fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + +TEST(HhtReducedOperators, DirectSpringScattersZeroInitialHessianByConstraintRow) { + CheckShellConstraintScatter(fesa::Mitc4DrillingMode::kDirectSpring); +} + TEST(HhtReducedEmbedding, RejectsNonfiniteGraphResidualBeforeAcceptance) { const auto fixture = MakeFixture(false); const auto chart = fesa::hht_reduced_internal::BuildFixedChart( @@ -1184,11 +1206,15 @@ TEST(HhtReducedEmbedding, EXPECT_DOUBLE_EQ(rates.Value().z_acceleration[0U], 1.0); } -TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { +void CheckActualShellRatesWithMixedTimeLoad(const fesa::Mitc4DrillingMode mode) { const auto fixture = MakeFixture(false, true); const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); ASSERT_TRUE(shells.HasValue()); + SCOPED_TRACE(static_cast(mode)); + for (const auto& shell : shells.Value()) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, mode); const fesa::Vector zero(fixture.dofs.FullDofCount()); const auto initial = Assemble(fixture, shells.Value(), zero, zero, zero, zero, true); @@ -1470,13 +1496,26 @@ TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { ++row) { auto& omitted = omitted_external_gradient.rows[row].gradient; const auto& no_force = no_force_operators.Value().rows[row].gradient; - ASSERT_EQ(omitted.size(), no_force.size()); - for (std::size_t entry = 0U; entry < omitted.size(); ++entry) { - ASSERT_EQ(omitted[entry].dof, no_force[entry].dof); - saw_external_position_gradient = - saw_external_position_gradient || - omitted[entry].value != no_force[entry].value; - omitted[entry].value = no_force[entry].value; + if (mode == fesa::Mitc4DrillingMode::kLegacyCorotatedTwist) { + ASSERT_EQ(omitted.size(), no_force.size()); + for (std::size_t entry = 0U; entry < omitted.size(); ++entry) { + ASSERT_EQ(omitted[entry].dof, no_force[entry].dof); + saw_external_position_gradient = + saw_external_position_gradient || + omitted[entry].value != no_force[entry].value; + omitted[entry].value = no_force[entry].value; + } + } else { + // Direct internal/external derivatives can cancel exact sparse entries. + // Compare the full union before replacing only the external contribution. + fesa::Vector with_force(fixture.dofs.FullDofCount()); + fesa::Vector without_force(fixture.dofs.FullDofCount()); + for (const auto& entry : omitted) with_force[entry.dof] = entry.value; + for (const auto& entry : no_force) without_force[entry.dof] = entry.value; + for (std::size_t dof = 0U; dof < with_force.Size(); ++dof) + saw_external_position_gradient = saw_external_position_gradient || + with_force[dof] != without_force[dof]; + omitted = no_force; } } ASSERT_TRUE(saw_external_position_gradient); @@ -1560,6 +1599,14 @@ TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { EXPECT_GT(omitted_second_time_hidden, 1.0e-8); } +TEST(HhtReducedEmbedding, ReconstructsActualShellRatesWithMixedTimeLoad) { + CheckActualShellRatesWithMixedTimeLoad(fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + +TEST(HhtReducedEmbedding, DirectSpringReconstructsActualShellRatesWithMixedTimeLoad) { + CheckActualShellRatesWithMixedTimeLoad(fesa::Mitc4DrillingMode::kDirectSpring); +} + TEST(HhtReducedChannelNewton, RequiresEveryTranslationGateAndThaws) { fesa::hht_reduced_internal::FixedChart chart{}; chart.full_dof_count = 6U; @@ -1737,3 +1784,178 @@ TEST(HhtReducedOperators, NullMetricKeepsFiniteLeftAssociatedProducts) { } } // namespace + +TEST(Mitc4DrillingEvidence, HhtActualMomentMatchesClosedFormAndSourceOrder) { + const auto fixture = MakeFixture(false, true); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + for (std::size_t node = 0U; node < fixture.domain->Nodes().size(); ++node) { + q[6U * node + 0U] = 0.03 * static_cast(node); + q[6U * node + 3U] = 0.13 + 0.01 * static_cast(node); + q[6U * node + 4U] = -0.09; + q[6U * node + 5U] = 0.21 - 0.01 * static_cast(node); + } + const auto serial = Assemble(fixture, shells.Value(), q, zero, zero, zero, true); + ASSERT_TRUE(serial.HasValue()); + ASSERT_TRUE(serial.Value().frame.mitc4_drilling); + fesa::Vector3 expected; + for (std::size_t element = 0U; element < fixture.domain->ShellElements().Size(); ++element) { + const auto& definition = fixture.domain->ShellElements()[element]; + std::array nodes{}; + std::array, 4> directors{}; + for (std::size_t i = 0U; i < 4U; ++i) { + nodes[i] = &fixture.domain->Nodes()[definition.node_indices[i]]; + directors[i] = fixture.domain->ShellNodeInitialFrames()[definition.node_indices[i]].director; + } + const auto* section = dynamic_cast( + &fixture.domain->Properties()[definition.section_index]); + const auto* material = dynamic_cast( + &fixture.domain->Materials()[definition.material_index]); + ASSERT_NE(section, nullptr); + ASSERT_NE(material, nullptr); + const auto kernel = fesa::Mitc4Shell::Create(nodes, directors, *section, *material); + ASSERT_TRUE(kernel.HasValue()); + const auto stiffness = kernel.Value().Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double trace = 0.0; + for (std::size_t i = 0U; i < 4U; ++i) + for (std::size_t c = 3U; c < 5U; ++c) + trace += stiffness.Value().physical_local20(5U * i + c, 5U * i + c); + const double kd = 1e-3 * trace / 8.0; + for (std::size_t i = 0U; i < 4U; ++i) { + const auto node = definition.node_indices[i]; + const fesa::Vector3 phi{q[6U * node + 3U], q[6U * node + 4U], q[6U * node + 5U]}; + const fesa::Vector3 director{directors[i]}; + const double angle = phi.Norm(); + const double coefficient = (1.0 - 0.5 * angle / std::tan(0.5 * angle)) / (angle * angle); + // Independent closed inverse-transpose Rodrigues identity; no dense solve. + expected = expected + kd * phi.Dot(director) * + (director + 0.5 * phi.Cross(director) + coefficient * phi.Cross(phi.Cross(director))); + } + } + for (std::size_t c = 0U; c < 3U; ++c) + EXPECT_NEAR(serial.Value().frame.mitc4_drilling->drilling_moment[c], expected[c], 2e-14); + EXPECT_EQ(serial.Value().frame.mitc4_drilling->mode, fesa::Mitc4DrillingMode::kDirectSpring); + auto frame = serial.Value().frame; + const auto parameters = fesa::MakeHhtParameters(-0.05); + ASSERT_TRUE(parameters.HasValue()); + fesa::TbbParallelFor parallel; + const auto tbb = fesa::hht_internal::EvaluateAssembly(*fixture.domain, fixture.dofs, + shells.Value(), frame, zero, 0.1, parameters.Value(), parallel, true); + ASSERT_TRUE(tbb.HasValue()); + ASSERT_TRUE(tbb.Value().frame.mitc4_drilling); + for (std::size_t c = 0U; c < 3U; ++c) + ExpectSameBits(tbb.Value().frame.mitc4_drilling->drilling_moment[c], + serial.Value().frame.mitc4_drilling->drilling_moment[c]); +} + +TEST(Mitc4DrillingEvidence, HhtCorrectedGateRejectsOppositeSignAndRealImbalance) { + const auto fixture = MakeFixture(false); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + for (std::size_t node = 0U; node < 4U; ++node) q[6U * node + 5U] = 0.125; + const auto assembly = Assemble(fixture, shells.Value(), q, zero, zero, zero, false); + ASSERT_TRUE(assembly.HasValue()); + EXPECT_EQ(assembly.Value().hidden_constraint.Size(), 0U); + const auto directional = Assemble(fixture, shells.Value(), q, zero, zero, zero, true); + ASSERT_TRUE(directional.HasValue()); + ASSERT_EQ(directional.Value().hidden_constraint.Size(), 4U); + for (std::size_t node = 0U; node < 4U; ++node) + EXPECT_GT(directional.Value().hidden_constraint[node], 0.0); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + auto unavailable_frame = assembly.Value().frame; + const auto unavailable = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), 2.0, + fesa::hht_internal::NullNodes(fixture.dofs, assembly.Value().rotations), + positive, tangent, &unavailable_frame); + EXPECT_FALSE(unavailable.IsOk()); + EXPECT_EQ(unavailable_frame.hidden_constraint.Size(), 0U); + auto frame = assembly.Value().frame; + // A support supplying the complete internal force gives exact static balance. + frame.reaction = frame.internal_force; + ASSERT_TRUE(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs, + assembly.Value(), &frame).IsOk()); + ASSERT_TRUE(frame.mitc4_drilling); + EXPECT_GT(frame.global_equilibrium[5U], 0.0); + EXPECT_NEAR(frame.global_equilibrium[5U], frame.mitc4_drilling->drilling_moment[2U], 1e-14); + for (const double value : frame.mitc4_drilling->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-14); + for (const int kind : {0, 1, 2, 3}) { + auto invalid = frame; + if (kind == 0) invalid.mitc4_drilling->drilling_moment[2U] *= -1.0; + if (kind == 1) invalid.reaction[0U] += 1.0; + if (kind == 2) invalid.reaction[4U] += 1.0; + if (kind == 3) invalid.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + EXPECT_FALSE(fesa::hht_internal::CompleteEquilibrium(*fixture.domain, fixture.dofs, + assembly.Value(), &invalid).IsOk()) << kind; + } +} + +TEST(Mitc4DrillingEvidence, DirectHiddenCertificateUsesCanonicalPrimitive) { + const auto fixture = MakeFixture(false); + const auto shells = fesa::hht_internal::BuildRuntimeShells(*fixture.domain, fixture.dofs); + ASSERT_TRUE(shells.HasValue()); + const fesa::Vector zero(fixture.dofs.FullDofCount()); + auto q = zero; + auto external = zero; + fesa::NonlinearMitc4TrialDofs trial{}; + for (std::size_t node = 0U; node < 4U; ++node) { + q[6U * node] = 0.03 * static_cast(node); + trial.displacement[node] = {q[6U * node], 0.0, 0.0}; + trial.rotation_vector[node] = {0.00005, 0.000035, 0.00002}; + for (std::size_t c = 0U; c < 3U; ++c) + q[6U * node + 3U + c] = trial.rotation_vector[node][c]; + } + const auto local = shells.Value()[0U].internal->EvaluateTrial(trial, &zero); + ASSERT_TRUE(local.HasValue()); + external = local.Value().drilling_force; + const auto assembly = Assemble(fixture, shells.Value(), q, zero, zero, external, true); + ASSERT_TRUE(assembly.HasValue()); + const auto length = fesa::hht_internal::ReferenceLength(*fixture.domain); + ASSERT_TRUE(length.HasValue()); + const auto null_nodes = fesa::hht_internal::NullNodes(fixture.dofs, assembly.Value().rotations); + ASSERT_EQ(null_nodes.size(), 4U); + fesa::MklPardisoSolver positive; + fesa::MklPardisoTangentSolver tangent; + auto frame = assembly.Value().frame; + const auto status = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, assembly.Value(), length.Value(), null_nodes, + positive, tangent, &frame); + ASSERT_TRUE(status.IsOk()) << (status.Diagnostics().empty() ? "" : status.Diagnostics()[0U].message); + const auto chart = fesa::hht_reduced_internal::BuildFixedChart(fixture.dofs, assembly.Value().rotations); + ASSERT_TRUE(chart.HasValue()); + const auto canonical = fesa::hht_reduced_internal::EvaluateConstraintOperators( + *fixture.domain, fixture.dofs, shells.Value(), assembly.Value(), chart.Value(), zero, zero); + ASSERT_TRUE(canonical.HasValue()); + for (std::size_t row = 0U; row < 4U; ++row) { + EXPECT_DOUBLE_EQ(frame.hidden_constraint[row], canonical.Value().rows[row].value); + EXPECT_NEAR(frame.hidden_constraint[row], 0.0, 1e-17); + EXPECT_DOUBLE_EQ(frame.hidden_velocity_defect[row], 0.0); + EXPECT_DOUBLE_EQ(frame.hidden_acceleration_defect[row], 0.0); + } + const auto accepted = frame; + for (const bool missing : {true, false}) { + auto invalid = assembly.Value(); + if (missing) + invalid.hidden_constraint = fesa::Vector(0U); + else + invalid.hidden_constraint[0U] = std::numeric_limits::infinity(); + const auto rejected = fesa::hht_internal::CompleteConstraintEvidence( + *fixture.domain, fixture.dofs, invalid, length.Value(), null_nodes, + positive, tangent, &frame); + EXPECT_FALSE(rejected.IsOk()); + for (const auto& pair : { + std::make_pair(&frame.hidden_constraint, &accepted.hidden_constraint), + std::make_pair(&frame.hidden_velocity_defect, &accepted.hidden_velocity_defect), + std::make_pair(&frame.hidden_acceleration_defect, &accepted.hidden_acceleration_defect)}) { + ASSERT_EQ(pair.first->Size(), pair.second->Size()); + for (std::size_t row = 0U; row < pair.first->Size(); ++row) + ExpectSameBits((*pair.first)[row], (*pair.second)[row]); + } + } +} diff --git a/tests/unit/analysis/hht_shell_physics_test.cpp b/tests/unit/analysis/hht_shell_physics_test.cpp index abac0f0..733ec91 100644 --- a/tests/unit/analysis/hht_shell_physics_test.cpp +++ b/tests/unit/analysis/hht_shell_physics_test.cpp @@ -28,6 +28,7 @@ #include "fesa/analysis/newton_iteration_log.h" #include "fesa/assembly/parallel_for.h" #include "fesa/elements/finite_rotation.h" +#include "elements/mitc4_drilling_internal.h" #include "fesa/solvers/linear/mkl_pardiso_solver.h" #include "fesa/solvers/linear/mkl_pardiso_tangent_solver.h" @@ -120,6 +121,14 @@ struct Fixture { double length; }; +// The finite rigid-orbit and free-spin C08 fixtures prescribe the historical +// objective potential and its hidden velocity equation explicitly. +void UseLegacyDrilling(const Fixture& fixture) { + for (const auto& shell : fixture.shells) + *shell.internal_kernel = fesa::mitc4_drilling_internal::TestAccess::WithMode( + *shell.internal_kernel, fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + Fixture MakeFixture(const std::size_t divisions, const bool lifecycle_translation = false) { fesa::ModelDefinition definition{}; @@ -643,6 +652,7 @@ TEST(HhtShellPhysics, ManufacturedHiddenCountsExternalHessianOnce) { TEST(HhtShellPhysics, C08aRigidOrbitMatchesIndependentVolumeForcePowerAndMomentum) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); const auto parameters = Require(fesa::MakeHhtParameters(-0.05)); const auto initial = Evaluate(fixture, EmptyFrame(fixture), 1.0 / 16.0, parameters, false); @@ -3078,6 +3088,7 @@ std::array AccelerationSplit(const Frame& basis, TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); const auto grid = RunPlate(fixture, 1.0 / 64.0, true); ASSERT_EQ(grid.history.size(), 65U); const auto exact = RigidOrbit(fixture, 0.0); @@ -3115,6 +3126,7 @@ TEST(HhtShellPhysics, ReducedManufacturedControllerKeepsExactInitialOrbit) { TEST(HhtShellPhysics, C08aRigidManufacturedOrbitReducedCoordinateTimeConvergence) { const auto fixture = MakeFixture(1U); + UseLegacyDrilling(fixture); std::vector grids; std::vector> errors; for (const unsigned denominator : {16U, 32U, 64U, 128U, 256U, 512U, 1024U}) { @@ -3297,6 +3309,7 @@ TEST(HhtShellPhysics, C08bInitialPhysicalSpectrumProbe) { TEST(HhtShellPhysics, PredictorStartCompletesPreviouslyFailingLumpedFreePlateGrid) { const auto fixture = MakeFixture(2U); + UseLegacyDrilling(fixture); const auto evidence = RunPlate(fixture, 1.0 / 256.0, false); ASSERT_EQ(evidence.history.size(), 65U); EXPECT_DOUBLE_EQ(evidence.history.back().time, kEndTime); @@ -3307,6 +3320,7 @@ TEST(HhtShellPhysics, TEST(HhtShellPhysics, C08bFreeFiniteRotationPlateSelfConvergenceAndPhysicalInvariants) { const auto fixture = MakeFixture(2U); + UseLegacyDrilling(fixture); std::vector grids; for (const unsigned denominator : {64U, 128U, 256U, 512U, 1024U, 32768U, 65536U, 131072U}) { diff --git a/tests/unit/analysis/nonlinear_analysis_state_test.cpp b/tests/unit/analysis/nonlinear_analysis_state_test.cpp index 4fc9ad3..7c2aa06 100644 --- a/tests/unit/analysis/nonlinear_analysis_state_test.cpp +++ b/tests/unit/analysis/nonlinear_analysis_state_test.cpp @@ -454,3 +454,28 @@ TEST(NonlinearAnalysisState, StoresNoUnapprovedVerificationMetrics) { EXPECT_FALSE((HasVerificationMetrics::value)); EXPECT_FALSE((HasVerificationMetrics::value)); } + +TEST(NonlinearAnalysisState, DrillingEvidenceRetainsModeAndRejectsPartialCommit) { + const auto fixture = MakeFixture(); + auto result = fesa::NonlinearAnalysisState::Create(*fixture.domain, *fixture.dofs, + MakeBaseline(*fixture.domain)); + ASSERT_TRUE(result.HasValue()); + auto& state = result.Value(); + auto first = MakeFrame(state.CommittedState(), 0U, 1U, 0.5, 1U); + first.mitc4_drilling = fesa::Mitc4DrillingEvidence{fesa::Mitc4DrillingMode::kDirectSpring}; + ASSERT_TRUE(state.CommitAcceptedFrame(first).IsOk()); + for (const int kind : {0, 1, 2, 3, 4}) { + auto next = MakeFrame(state.CommittedState(), 1U, 2U, 1.0, 1U); + next.mitc4_drilling = first.mitc4_drilling; + if (kind == 0) next.mitc4_drilling.reset(); + if (kind == 1) next.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) next.mitc4_drilling->mode = static_cast(99); + if (kind == 3) next.mitc4_drilling->equilibrium_with_drilling[3U] = 1.0; + if (kind == 4) next.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::infinity(); + EXPECT_FALSE(state.CommitAcceptedFrame(next).IsOk()) << kind; + ASSERT_EQ(state.AcceptedFrames().size(), 1U); + EXPECT_DOUBLE_EQ(state.CommittedState().load_factor, 0.5); + EXPECT_EQ(state.CommittedState().mitc4_drilling->mode, + fesa::Mitc4DrillingMode::kDirectSpring); + } +} diff --git a/tests/unit/elements/element_capabilities_test.cpp b/tests/unit/elements/element_capabilities_test.cpp index 585d749..7e6b556 100644 --- a/tests/unit/elements/element_capabilities_test.cpp +++ b/tests/unit/elements/element_capabilities_test.cpp @@ -133,6 +133,7 @@ TEST(ElementCapabilities, ShellTrialPreservesCompleteDirectionalResponseAndFailu Same(a.drilling_metric, b.drilling_metric); Same(a.drilling_energy, b.drilling_energy); Same(a.force_directional_second, b.force_directional_second); + Same(a.drilling_force_directional_second, b.drilling_force_directional_second); Same(a.drilling_constraint, b.drilling_constraint); Same(a.drilling_constraint_gradient, b.drilling_constraint_gradient); for (std::size_t n = 0; n < 4; ++n) { diff --git a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp index 5256a2b..20bcb00 100644 --- a/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp +++ b/tests/unit/elements/nonlinear_mitc4_shell_tangent_test.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -13,6 +14,8 @@ #include "fesa/elements/nonlinear_mitc4_inertia.h" #include "fesa/elements/nonlinear_mitc4_shell.h" #include "fesa/model/model_types.h" +#include "elements/mitc4_drilling_internal.h" +#include "elements/nonlinear_mitc4_configuration_internal.h" namespace { @@ -44,7 +47,7 @@ std::array Frames() { {3U, {0.0, 0.0, 1.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}}}}; } -fesa::NonlinearMitc4Shell Shell(const std::array& nodes) { +fesa::NonlinearMitc4Shell ProductionShell(const std::array& nodes) { const fesa::ShellSection section{ "Section", 0.2, 0U, {"nonlinear-mitc4-tangent.inp", 20U}}; const fesa::LinearElasticMaterial material{ @@ -57,6 +60,11 @@ fesa::NonlinearMitc4Shell Shell(const std::array& nodes) { return shell.Value(); } +fesa::NonlinearMitc4Shell LegacyShell(const std::array& nodes) { + return fesa::mitc4_drilling_internal::TestAccess::WithMode( + ProductionShell(nodes), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); +} + fesa::NonlinearMitc4TrialDofs Trial() { fesa::NonlinearMitc4TrialDofs trial{}; trial.displacement[0U] = {-0.02, 0.01, 0.03}; @@ -200,7 +208,7 @@ void ExpectRelative(const double actual, const double expected) { // NSNR-ELEM-TAN-001 TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto recovery = shell.RecoverState(trial); const auto result = fesa::EvaluateInternalResidualAndTangent(shell, trial); @@ -230,7 +238,7 @@ TEST(NonlinearMitc4Residual, EvaluatesResidualTangentAndRecoveryAtOneTrial) { // NSNR-ELEM-TAN-002 and NSNR-ELEM-TAN-003 TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto direction = WorkConjugateDirection(nodes); const double length = ReferenceLength(nodes); @@ -295,7 +303,7 @@ TEST(NonlinearMitc4Tangent, MatchesTheApprovedDirectionalSchedule) { TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { const auto nodes = Nodes(); const auto equilibrium = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), {}); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), {}); ASSERT_TRUE(equilibrium.HasValue()); const double length = ReferenceLength(nodes); @@ -326,9 +334,9 @@ TEST(NonlinearMitc4Tangent, PreservesRigidTranslationNullspaceAndSymmetry) { TEST(NonlinearMitc4Residual, DoesNotAcceptExternalLoadShares) { const auto nodes = Nodes(); const auto first = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial()); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial()); const auto second = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), Trial()); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), Trial()); ASSERT_TRUE(first.HasValue()); ASSERT_TRUE(second.HasValue()); for (std::size_t index = 0U; index < 24U; ++index) { @@ -357,7 +365,7 @@ TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) { trial.rotation_vector[node] = rotation; } const auto result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial); ASSERT_TRUE(result.HasValue()); EXPECT_LE(result.Value().internal_residual.Norm(), kNumericalTolerance); EXPECT_LE(result.Value().physical_strain_energy, kNumericalTolerance); @@ -366,7 +374,7 @@ TEST(NonlinearMitc4Objectivity, FiniteRigidMotionHasZeroInternalForce) { // NSNR-ELEM-TAN-008 TEST(NonlinearMitc4Patch, RetainsPhysicalRecoveryIdentities) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); constexpr double kMembraneStrain = 1.0e-2; const double green_strain = kMembraneStrain + 0.5 * kMembraneStrain * kMembraneStrain; @@ -409,7 +417,7 @@ TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) { trial.rotation_vector[2U] = {std::numeric_limits::quiet_NaN(), 0.0, 0.0}; const auto result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), trial); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), trial); ASSERT_FALSE(result.HasValue()); EXPECT_EQ(result.GetStatus().Diagnostics()[0U].code, "invalid-nonlinear-shell-trial"); @@ -417,13 +425,13 @@ TEST(NonlinearMitc4Residual, RejectsNonfiniteTrialData) { auto overflow = Trial(); overflow.displacement[1U] = {1.0e300, 0.0, 0.0}; const auto overflow_result = - fesa::EvaluateInternalResidualAndTangent(Shell(nodes), overflow); + fesa::EvaluateInternalResidualAndTangent(LegacyShell(nodes), overflow); EXPECT_FALSE(overflow_result.HasValue()); } void CheckThirdSweep(const fesa::NonlinearMitc4TrialDofs& trial) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); auto direction = Direction(); double norm = 0.0; const auto unscaled = ToVector(direction); @@ -488,7 +496,7 @@ TEST(NonlinearMitc4Residual, nodes[1].coordinates = {4.23076923, 2.30769231, 0.0}; nodes[2].coordinates = {4.23076923, 2.69230769, 0.0}; nodes[3].coordinates = {3.84615385, 2.69230769, 0.0}; - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto evaluated = fesa::EvaluateInternalResidualAndTangent(shell, {}); ASSERT_TRUE(evaluated.HasValue()); for (std::size_t i = 0; i < 24U; ++i) @@ -566,7 +574,7 @@ TEST(NonlinearMitc4Residual, for (bool warped : {false, true}) { auto nodes = Nodes(); if (warped) nodes[2U].coordinates[2U] = 0.1; - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); auto trial = Trial(); trial.rotation_vector[0U] = {.6, -.2, .1}; const auto direction = Direction(); @@ -636,7 +644,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExposesExactDrillingConstraintThirdDirectionalGradient) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const auto trial = Trial(); const auto direction = Direction(); const auto velocity = ToVector(direction); @@ -679,7 +687,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExternalDrillingConstraintIncludesForceTimeAndCurvatureTerms) { const auto nodes = Nodes(); - const auto shell = Shell(nodes); + const auto shell = LegacyShell(nodes); const fesa::Vector3 phi{.6, -.2, .1}; const fesa::Vector3 phi_dot{.2, -.1, .3}; const fesa::Vector3 force{1.7, -.8, .4}; @@ -742,7 +750,7 @@ TEST(NonlinearMitc4Residual, TEST(NonlinearMitc4Residual, ExternalDrillingConstraintRejectsInvalidRotationAndHasZeroForceValue) { - const auto shell = Shell(Nodes()); + const auto shell = LegacyShell(Nodes()); const auto zero = shell.EvaluateExternalDrillingConstraint( 0U, {.2, -.1, .3}, {.1, .2, -.3}, {}, {}, {}); ASSERT_TRUE(zero.HasValue()); @@ -838,4 +846,470 @@ TEST(NonlinearMitc4Residual, TinyRotationPreservesPhysicalForceAndRecoveryEnergy } } +// Independent fixed-initial-director spring oracle: a total-chart rigid spin +// carries numerical energy even though the physical shell is undeformed. +TEST(Mitc4DrillingSelection, ProductionUsesDirectInitialDirectorSpring) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + const fesa::ShellSection section{"Section", 0.2, 0U, {"spring.inp", 1U}}; + const fesa::LinearElasticMaterial material{ + "Material", 120.0, 0.25, {"spring.inp", 2U}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + const auto stiffness = linear.Value().Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double rotary_trace = 0.0; + for (std::size_t node = 0; node < 4; ++node) + for (std::size_t c = 3; c < 5; ++c) + rotary_trace += + stiffness.Value().physical_local20(5 * node + c, 5 * node + c); + const double kd = 1e-3 * rotary_trace / 8.; + fesa::NonlinearMitc4TrialDofs trial{}; + const fesa::Vector3 phi{0., 0., 0.4}; + const auto rotation = fesa::ExpSkew(phi); + for (std::size_t node = 0; node < 4; ++node) { + const fesa::Vector3 x{nodes[node].coordinates}; + trial.rotation_vector[node] = phi; + fesa::Vector3 rotated{ + rotation(0, 0) * x[0] + rotation(0, 1) * x[1] + rotation(0, 2) * x[2], + rotation(1, 0) * x[0] + rotation(1, 1) * x[1] + rotation(1, 2) * x[2], + rotation(2, 0) * x[0] + rotation(2, 1) * x[1] + rotation(2, 2) * x[2]}; + trial.displacement[node] = rotated - x; + } + const auto actual = fesa::EvaluateInternalResidualAndTangent(shell, trial); + ASSERT_TRUE(actual.HasValue()); + EXPECT_NEAR(actual.Value().physical_strain_energy, 0., 1e-25); + EXPECT_NEAR(actual.Value().drilling_energy, .5 * kd * 4 * .4 * .4, 1e-14); + for (std::size_t row = 0; row < 24; ++row) { + const double force = row % 6 == 5 ? kd * .4 : 0.; + EXPECT_NEAR(actual.Value().drilling_force[row], force, 1e-14); + for (std::size_t col = 0; col < 24; ++col) { + const double tangent = row == col && row % 6 == 5 ? kd : 0.; + EXPECT_NEAR(actual.Value().drilling_tangent(row, col), tangent, 1e-14); + } + } +} + +TEST(Mitc4DrillingSelection, + DirectSpringUsesEachInitialDirectorAndConstantBlocks) { + for (const bool warped : {false, true}) { + auto nodes = Nodes(); + if (warped) nodes[2].coordinates[2] = .1; + auto frames = Frames(); + const auto transform = [](const fesa::Matrix& r, const Coordinates& x) { + Coordinates y{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) y[i] += r(i, j) * x[j]; + return y; + }; + const auto global = fesa::ExpSkew({.4, -.3, .2}); + std::array directors{}; + for (std::size_t n = 0; n < 4; ++n) { + const auto r = + fesa::ExpSkew({0., warped ? .03 * static_cast(n) : 0., 0.}); + frames[n].director = transform(global, transform(r, frames[n].director)); + frames[n].tangent_a = + transform(global, transform(r, frames[n].tangent_a)); + frames[n].tangent_b = + transform(global, transform(r, frames[n].tangent_b)); + nodes[n].coordinates = transform(global, nodes[n].coordinates); + directors[n] = frames[n].director; + } + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + const fesa::LinearElasticMaterial material{ + "Material", 120., .25, {"spring.inp", 2}}; + const auto shell = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), frames, section, material); + const auto linear = fesa::Mitc4Shell::Create(NodePointers(nodes), directors, + section, material); + ASSERT_TRUE(shell.HasValue()); + ASSERT_TRUE(linear.HasValue()); + const auto k = linear.Value().Stiffness(); + ASSERT_TRUE(k.HasValue()); + double trace = 0.; + for (std::size_t n = 0; n < 4; ++n) + for (std::size_t c = 3; c < 5; ++c) + trace += k.Value().physical_local20(5 * n + c, 5 * n + c); + const double kd = 1e-3 * trace / 8.; + auto trial = Trial(); + trial.rotation_vector[0] = {.65, -.4, .3}; + const auto direction = ToVector(Direction()); + const auto response = fesa::EvaluateInternalResidualAndTangent( + shell.Value(), trial, &direction); + const auto shifted = fesa::EvaluateInternalResidualAndTangent( + shell.Value(), AddScaled(trial, Direction(), .2), &direction); + ASSERT_TRUE(response.HasValue()); + ASSERT_TRUE(shifted.HasValue()); + double energy = 0.; + for (std::size_t n = 0; n < 4; ++n) { + const double gamma = + fesa::Vector3{directors[n]}.Dot(trial.rotation_vector[n]); + energy += .5 * kd * gamma * gamma; + const fesa::Vector3 d{directors[n]}; + const auto rotation = fesa::EvaluateDynamicRotation( + trial.rotation_vector[n], d, Direction().rotation_vector[n], {}); + ASSERT_TRUE(rotation.HasValue()); + EXPECT_NEAR(response.Value().drilling_constraint[n], + kd * gamma * rotation.Value().drill_null.Dot(d), 1e-14); + for (std::size_t c = 0; c < 3; ++c) { + double ds = 0.; + for (std::size_t j = 0; j < 3; ++j) + ds += d[j] * rotation.Value().null_gradient(j, c); + EXPECT_NEAR( + response.Value().drilling_constraint_gradient(n, 6 * n + 3 + c), + kd * (d[c] * rotation.Value().drill_null.Dot(d) + gamma * ds), + 1e-14); + } + for (std::size_t c = 0; c < 3; ++c) + EXPECT_NEAR(response.Value().drilling_force[6 * n + 3 + c], + kd * gamma * directors[n][c], 1e-14); + } + EXPECT_NEAR(response.Value().drilling_energy, energy, 1e-14); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_DOUBLE_EQ(response.Value().drilling_force_directional_second[i], + 0.); + if (i % 6 < 3) EXPECT_DOUBLE_EQ(response.Value().drilling_force[i], 0.); + for (std::size_t j = 0; j < 24; ++j) { + const double expected = + i / 6 == j / 6 && i % 6 >= 3 && j % 6 >= 3 + ? kd * directors[i / 6][i % 6 - 3] * directors[i / 6][j % 6 - 3] + : 0.; + EXPECT_NEAR(response.Value().drilling_tangent(i, j), expected, 1e-14); + EXPECT_NEAR(response.Value().drilling_metric(i, j), expected, 1e-14); + EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j), + shifted.Value().drilling_tangent(i, j)); + if (expected == 0.) + EXPECT_DOUBLE_EQ(response.Value().drilling_tangent(i, j), 0.); + } + } + } +} + +TEST(Mitc4DrillingSelection, + DirectConstraintUsesPhysicalNullAndItsHighDerivatives) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + auto trial = Trial(); + trial.rotation_vector[0] = {.65, -.4, .3}; + const auto direction = Direction(); + const auto velocity = ToVector(direction); + const auto response = + fesa::EvaluateInternalResidualAndTangent(shell, trial, &velocity); + ASSERT_TRUE(response.HasValue()); + const auto zero = + fesa::EvaluateInternalResidualAndTangent(shell, {}, &velocity); + ASSERT_TRUE(zero.HasValue()); + const double kd = zero.Value().drilling_tangent(5, 5); + // Independent inertia API provides n, Dn and D2n, never a twist angle. + const fesa::Vector3 d{0., 0., 1.}; + const auto scalar_gradient = [&](const fesa::Vector3& phi, std::size_t c) { + const auto r = fesa::EvaluateDynamicRotation(phi, d, {}, {}); + EXPECT_TRUE(r.HasValue()); + double projection_derivative = 0.; + for (std::size_t k = 0; k < 3; ++k) + projection_derivative += r.Value().null_gradient(k, c) * d[k]; + return kd * (d[c] * r.Value().drill_null.Dot(d) + + d.Dot(phi) * projection_derivative); + }; + constexpr double h = .002; + for (std::size_t n = 0; n < 4; ++n) { + const auto phi = trial.rotation_vector[n]; + const auto v = direction.rotation_vector[n]; + const auto r = fesa::EvaluateDynamicRotation(phi, d, v, {}); + ASSERT_TRUE(r.HasValue()); + const double gamma = d.Dot(phi); + const double gamma_dot = d.Dot(v); + double projection_dot = 0.; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + projection_dot += d[i] * r.Value().null_gradient(i, j) * v[j]; + const double expected_second = + kd * (2. * gamma_dot * projection_dot + + gamma * r.Value().null_directional_second.Dot(d)); + double actual_second = 0.; + EXPECT_NEAR(response.Value().drilling_constraint[n], + kd * gamma * r.Value().drill_null.Dot(d), 1e-14); + for (std::size_t i = 0; i < 24; ++i) { + const double gradient = + i / 6 == n && i % 6 >= 3 ? scalar_gradient(phi, i % 6 - 3) : 0.; + EXPECT_NEAR(response.Value().drilling_constraint_gradient(n, i), gradient, + 1e-14); + for (std::size_t j = 0; j < 24; ++j) + actual_second += velocity[i] * + response.Value().drilling_constraint_hessian[n](i, j) * + velocity[j]; + if (i / 6 == n && i % 6 >= 3) { + const auto c = i % 6 - 3; + const double fd = + (scalar_gradient(phi + h * v, c) - 2. * scalar_gradient(phi, c) + + scalar_gradient(phi - h * v, c)) / + (h * h); + const double actual = + response.Value() + .drilling_constraint_third_directional_gradient[n][i]; + EXPECT_NEAR(actual, fd, + 2e-5 * (std::abs(actual) + std::abs(fd) + kd * 1e-4)); + } else { + EXPECT_DOUBLE_EQ( + response.Value() + .drilling_constraint_third_directional_gradient[n][i], + 0.); + } + } + EXPECT_NEAR(actual_second, expected_second, 1e-14); + } +} + +TEST(Mitc4DrillingSelection, DirectExternalConstraintUsesActualInertiaNull) { + const auto nodes = Nodes(); + const auto shell = ProductionShell(nodes); + const fesa::Vector3 phi{.65, -.4, .3}, velocity{.2, .1, -.15}; + const fesa::Vector3 force{.7, -.3, .2}, rate{-.1, .2, .4}, + acceleration{.5, -.2, .1}; + const auto actual = shell.EvaluateExternalDrillingConstraint( + 0, phi, velocity, force, rate, acceleration); + const auto rotation = + fesa::EvaluateDynamicRotation(phi, {0., 0., 1.}, velocity, {}); + ASSERT_TRUE(actual.HasValue()); + ASSERT_TRUE(rotation.HasValue()); + const double h = .002; + const auto plus = shell.EvaluateExternalDrillingConstraint( + 0, phi + h * velocity, velocity, + force + h * rate + .5 * h * h * acceleration, rate + h * acceleration, + acceleration); + const auto minus = shell.EvaluateExternalDrillingConstraint( + 0, phi - h * velocity, velocity, + force - h * rate + .5 * h * h * acceleration, rate - h * acceleration, + acceleration); + ASSERT_TRUE(plus.HasValue()); + ASSERT_TRUE(minus.HasValue()); + for (std::size_t c = 0; c < 3; ++c) { + EXPECT_NEAR( + actual.Value().gradient_directional_first[c], + (plus.Value().gradient[c] - minus.Value().gradient[c]) / (2. * h), + 1e-7); + EXPECT_NEAR(actual.Value().gradient_directional_second[c], + (plus.Value().gradient[c] - 2. * actual.Value().gradient[c] + + minus.Value().gradient[c]) / + (h * h), + 2e-7); + std::array basis{}; + basis[c] = 1.; + constexpr double coordinate_h = 1e-4; + const auto positive = shell.EvaluateExternalDrillingConstraint( + 0, phi + coordinate_h * fesa::Vector3{basis}, velocity, force, rate, + acceleration); + const auto negative = shell.EvaluateExternalDrillingConstraint( + 0, phi - coordinate_h * fesa::Vector3{basis}, velocity, force, rate, + acceleration); + ASSERT_TRUE(positive.HasValue()); + ASSERT_TRUE(negative.HasValue()); + for (std::size_t j = 0; j < 3; ++j) + EXPECT_NEAR( + actual.Value().hessian(j, c), + (positive.Value().gradient[j] - negative.Value().gradient[j]) / + (2. * coordinate_h), + 1e-7); + } + const auto& n = rotation.Value(); + EXPECT_NEAR(actual.Value().value, -n.drill_null.Dot(force), 1e-14); + fesa::Vector3 n_dot{}; + std::array values{}; + for (std::size_t i = 0; i < 3; ++i) + for (std::size_t j = 0; j < 3; ++j) + values[i] += n.null_gradient(i, j) * velocity[j]; + n_dot = fesa::Vector3{values}; + EXPECT_NEAR(actual.Value().directional_first, + -n_dot.Dot(force) - n.drill_null.Dot(rate), 1e-14); + EXPECT_NEAR(actual.Value().directional_second, + -n.null_directional_second.Dot(force) - 2. * n_dot.Dot(rate) - + n.drill_null.Dot(acceleration), + 1e-14); + for (std::size_t c = 0; c < 3; ++c) { + double expected = 0.; + for (std::size_t i = 0; i < 3; ++i) + expected -= n.null_gradient(i, c) * force[i]; + EXPECT_NEAR(actual.Value().gradient[c], expected, 1e-14); + } +} + +TEST(Mitc4DrillingSelection, + ModesAreImmutableReferenceInputsAndPreservePhysicalResponse) { + const auto nodes = Nodes(); + const auto direct = ProductionShell(nodes); + const auto legacy = LegacyShell(nodes); + EXPECT_EQ(direct.DrillingMode(), fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_EQ(legacy.DrillingMode(), + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + EXPECT_STREQ(fesa::Mitc4DrillingModeName(direct.DrillingMode()), + "direct-initial-director-spring"); + EXPECT_STREQ(fesa::Mitc4DrillingModeName(legacy.DrillingMode()), + "legacy-corotated-twist"); + const auto a = + fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(direct); + const auto b = + fesa::nonlinear_mitc4_internal::ConfigurationSnapshot::Read(legacy); + EXPECT_EQ(a.values, b.values); + EXPECT_NE(a.drilling_mode, b.drilling_mode); + const auto direction = ToVector(Direction()); + const auto r = + fesa::EvaluateInternalResidualAndTangent(direct, Trial(), &direction); + const auto l = + fesa::EvaluateInternalResidualAndTangent(legacy, Trial(), &direction); + ASSERT_TRUE(r.HasValue()); + ASSERT_TRUE(l.HasValue()); + EXPECT_DOUBLE_EQ(r.Value().physical_strain_energy, + l.Value().physical_strain_energy); + EXPECT_NE(r.Value().drilling_energy, l.Value().drilling_energy); + for (std::size_t i = 0; i < 24; ++i) { + EXPECT_DOUBLE_EQ(r.Value().physical_force[i], l.Value().physical_force[i]); + for (std::size_t j = 0; j < 24; ++j) + EXPECT_NEAR( + r.Value().internal_tangent(i, j) - r.Value().drilling_tangent(i, j), + l.Value().internal_tangent(i, j) - l.Value().drilling_tangent(i, j), + 1e-14); + } + for (std::size_t gp = 0; gp < 4; ++gp) + EXPECT_EQ(r.Value().recovery.gauss_points[gp].section_resultant, + l.Value().recovery.gauss_points[gp].section_resultant); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + const fesa::LinearElasticMaterial material{ + "Material", 120., .25, {"spring.inp", 2}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + EXPECT_EQ(linear.Value().DrillingMode(), direct.DrillingMode()); + const auto alternate = fesa::mitc4_drilling_internal::TestAccess::WithMode( + linear.Value(), fesa::Mitc4DrillingMode::kLegacyCorotatedTwist); + const auto k1 = linear.Value().Stiffness(); + const auto k2 = alternate.Stiffness(); + ASSERT_TRUE(k1.HasValue()); + ASSERT_TRUE(k2.HasValue()); + for (std::size_t i = 0; i < 24; ++i) + for (std::size_t j = 0; j < 24; ++j) + EXPECT_DOUBLE_EQ(k1.Value().stabilized_global24(i, j), + k2.Value().stabilized_global24(i, j)); +} + +TEST(Mitc4DrillingSelection, DirectRejectsInvalidScaleAndNonfiniteTrial) { + const auto nodes = Nodes(); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + for (double youngs_modulus : {1e-320, (std::numeric_limits::max)()}) { + const fesa::LinearElasticMaterial material{ + "Material", youngs_modulus, .25, {"spring.inp", 2}}; + const auto invalid = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), Frames(), section, material); + ASSERT_FALSE(invalid.HasValue()); + EXPECT_EQ(invalid.GetStatus().Diagnostics()[0].code, + youngs_modulus < 1.0 ? "invalid-nonlinear-shell-drilling" + : "invalid-shell-material"); + } + const auto shell = ProductionShell(nodes); + for (double value : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity()}) { + auto trial = Trial(); + trial.rotation_vector[0] = {value, 0., 0.}; + EXPECT_FALSE( + fesa::EvaluateInternalResidualAndTangent(shell, trial).HasValue()); + } + const auto valid = fesa::EvaluateInternalResidualAndTangent(shell, {}); + ASSERT_TRUE(valid.HasValue()); + EXPECT_DOUBLE_EQ(valid.Value().drilling_energy, 0.); +} + +TEST(Mitc4DrillingSelection, DirectRetainsPositiveSubnormalSpringScale) { + const auto nodes = Nodes(); + const fesa::ShellSection section{"Section", .2, 0, {"spring.inp", 1}}; + unsigned accepted = 0; + for (double youngs_modulus : + {1e-315, 1e-316, 1e-317, 1e-318, 2e-318, 3e-318, 1e-319}) { + const fesa::LinearElasticMaterial material{ + "Material", youngs_modulus, .25, {"spring.inp", 2}}; + const auto linear = fesa::Mitc4Shell::Create( + NodePointers(nodes), + {{{0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}, {0., 0., 1.}}}, section, + material); + ASSERT_TRUE(linear.HasValue()); + const auto stiffness = linear.Value().Stiffness(); + const auto shell = fesa::NonlinearMitc4Shell::Create( + NodePointers(nodes), Frames(), section, material); + if (!stiffness.HasValue()) { + EXPECT_FALSE(shell.HasValue()); + continue; + } + ASSERT_TRUE(shell.HasValue()); + ++accepted; + fesa::NonlinearMitc4TrialDofs trial{}; + trial.rotation_vector[0] = {0., 0., 1.}; + const auto actual = + fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial); + ASSERT_TRUE(actual.HasValue()); + // Scale is already validated. No later half-scale rounding may change K. + const double kd = stiffness.Value().drilling_stiffness; + EXPECT_EQ(actual.Value().drilling_tangent(5, 5), kd); + EXPECT_EQ(actual.Value().drilling_force[5], kd); + trial.rotation_vector[0] = {0., 0., 2.}; + const auto double_angle = + fesa::EvaluateInternalResidualAndTangent(shell.Value(), trial); + ASSERT_TRUE(double_angle.HasValue()); + // U(2) = 2*kd is representable even for the least positive spring scale. + EXPECT_EQ(double_angle.Value().drilling_energy, 2. * kd); + EXPECT_EQ(double_angle.Value().drilling_force[5], 2. * kd); + } + EXPECT_GT(accepted, 0U); +} + +TEST(Mitc4DrillingSelection, LegacyHeadFrozenEvidence) { + const auto nodes = Nodes(); + const auto shell = LegacyShell(nodes); + const auto direction = ToVector(Direction()); + const auto response = + fesa::EvaluateInternalResidualAndTangent(shell, Trial(), &direction); + ASSERT_TRUE(response.HasValue()); + const auto& r = response.Value(); + std::array hashes{}; + hashes.fill(14695981039346656037ULL); + const auto add = [&](std::size_t channel, double value) { + std::uint64_t bits{}; + std::memcpy(&bits, &value, sizeof(bits)); + for (unsigned b = 0; b < 8; ++b) { + hashes[channel] ^= (bits >> (8 * b)) & 255ULL; + hashes[channel] *= 1099511628211ULL; + } + }; + add(0, r.drilling_energy); + for (std::size_t i = 0; i < 24; ++i) { + add(1, r.drilling_force[i]); + add(4, r.force_directional_second[i]); + for (std::size_t j = 0; j < 24; ++j) { + add(2, r.drilling_tangent(i, j)); + add(3, r.drilling_metric(i, j)); + } + } + for (std::size_t n = 0; n < 4; ++n) { + add(5, r.drilling_constraint[n]); + for (std::size_t i = 0; i < 24; ++i) { + add(5, r.drilling_constraint_gradient(n, i)); + add(7, r.drilling_constraint_third_directional_gradient[n][i]); + for (std::size_t j = 0; j < 24; ++j) + add(6, r.drilling_constraint_hessian[n](i, j)); + } + } + // Captured before implementation from HEAD + // 52e321e15c70377d4318aace3f3bc64b3ea1a65c, native MSVC x64 Debug. Entire + // energy/force/tangent/metric/directional/psi inventories are pinned + // independently, including the historical polar1 constraint inconsistency. Do + // not regenerate these from the selected mode. + constexpr std::array frozen{ + 16395275265220000061ULL, 12405324468590059764ULL, 79034762627498281ULL, + 7305417493946306304ULL, 7191018337857616990ULL, 4845391414006668707ULL, + 11613163243225266145ULL, 3872161451216566265ULL}; + EXPECT_EQ(hashes, frozen); +} + } // namespace diff --git a/tests/unit/io/hdf5/dynamic_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/dynamic_hdf5_results_writer_test.cpp index 0de145c..2776a7b 100644 --- a/tests/unit/io/hdf5/dynamic_hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/dynamic_hdf5_results_writer_test.cpp @@ -23,6 +23,10 @@ #include "fesa/fem/dof_manager.h" #include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "mitc4_fault_test_helpers.h" + namespace { using fesa::hdf5_internal::Hdf5Handle; @@ -145,6 +149,8 @@ fesa::DynamicAcceptedFrameState MakeFrame(const fesa::Domain& domain, {fesa::ShellSectionPosition::kMiddle, 0.0, {}}, {fesa::ShellSectionPosition::kTop, 1.0, {}}}}}); } + frame.mitc4_drilling = fesa::Mitc4DrillingEvidence{ + fesa::Mitc4DrillingMode::kDirectSpring, {}, frame.global_equilibrium}; return frame; } @@ -670,3 +676,58 @@ TEST(DynamicHdf5ResultsWriter, FailedAtomicReplacementPreservesPreviousBytes) { } } // namespace + + +TEST(DynamicHdf5ResultsWriter, DrillingModesIncludeInitialAndEveryAcceptedFrame) { + for (const auto mode : {fesa::Mitc4DrillingMode::kDirectSpring, + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist}) { + auto fixture = MakeFixture(); + for (std::uint32_t i = 0; i < 3U; ++i) { + auto frame = MakeFrame(*fixture.domain, i); + frame.mitc4_drilling->mode = mode; + ASSERT_TRUE(fixture.state.CommitAcceptedFrame(frame).IsOk()); + } + TempFile path{"drilling-roundtrip"}; + ASSERT_TRUE(fesa::DynamicHdf5ResultsWriter{}.Write(path.Path(), *fixture.domain, fixture.state, {}).IsOk()); + Hdf5Handle file{H5Fopen(path.Path().string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), H5Fclose}; + const auto stored = fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()); + ASSERT_TRUE(stored); + EXPECT_EQ(stored->formulation, fesa::Mitc4DrillingModeName(mode)); + ASSERT_EQ(stored->frames.size(), 3U); + for (std::size_t i = 0; i < 3U; ++i) { + EXPECT_EQ(stored->frames[i].frame_index, i); + EXPECT_EQ(stored->frames[i].raw_equilibrium, fixture.state.AcceptedFrames()[i].global_equilibrium); + EXPECT_EQ(stored->frames[i].equilibrium_with_drilling, + fixture.state.AcceptedFrames()[i].mitc4_drilling->equilibrium_with_drilling); + } + } +} + +TEST(DynamicHdf5ResultsWriter, DrillingMissingActualStatePreservesPriorFinal) { + auto fixture = MakeFixture(); + for (std::uint32_t i = 0; i < 3U; ++i) { + auto frame = MakeFrame(*fixture.domain, i); + frame.mitc4_drilling.reset(); + ASSERT_TRUE(fixture.state.CommitAcceptedFrame(frame).IsOk()); + } + TempFile path{"drilling-missing"}; + std::ofstream(path.Path(), std::ios::binary) << "prior"; + EXPECT_FALSE(fesa::DynamicHdf5ResultsWriter{}.Write(path.Path(), *fixture.domain, fixture.state, {}).IsOk()); + std::ifstream stream{path.Path(), std::ios::binary}; + EXPECT_EQ((std::string{std::istreambuf_iterator{stream}, {}}), "prior"); +} + + +TEST(DynamicHdf5ResultsWriter, DrillingLifecycleFaultsPreserveFinalBytes) { + auto fixture = MakeFixture(); + CommitHistory(fixture); + TempFile path{"drilling-faults"}; + for (const std::string fault : {"candidate-close", "reader-close", "markers", "datasets"}) { + SCOPED_TRACE(fault); + std::ofstream(path.Path(), std::ios::binary | std::ios::trunc) << "prior"; + EXPECT_FALSE(fesa::hdf5_internal::WriteDynamicFile(path.Path(), *fixture.domain, fixture.state, {}, + fesa::test::Mitc4LifecycleFault(fault)).IsOk()); + std::ifstream stream{path.Path(), std::ios::binary}; + EXPECT_EQ((std::string{std::istreambuf_iterator{stream}, {}}), "prior"); + } +} diff --git a/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp index 926d667..e3510d1 100644 --- a/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/eigen_hdf5_results_writer_test.cpp @@ -19,6 +19,10 @@ #include "fesa/model/domain.h" #include "io/hdf5/eigen_hdf5_self_check.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "mitc4_fault_test_helpers.h" + namespace { fesa::ModelDefinition Definition(const std::filesystem::path& source) { @@ -81,6 +85,10 @@ TEST(EigenHdf5ResultsWriter, WritesModalSchemaAndRejectsTamperedState) { const hid_t file = H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); ASSERT_GE(file, 0); + const auto drilling = fesa::hdf5_internal::ReadDrillingDiagnostics(file); + ASSERT_TRUE(drilling); + EXPECT_EQ(drilling->formulation, fesa::Mitc4DrillingModeName(*state.drilling_mode)); + EXPECT_TRUE(drilling->frames.empty()); EXPECT_GT(H5Lexists(file, "/steps/Step-1/eigen/modes", H5P_DEFAULT), 0); const hid_t controls = H5Dopen2(file, "/steps/Step-1/eigen/backend/base_fpm", H5P_DEFAULT); @@ -105,7 +113,8 @@ TEST(EigenHdf5ResultsWriter, WritesModalSchemaAndRejectsTamperedState) { [](auto& candidate) { candidate.constraint_mask[0] = 1U; }, [](auto& candidate) { candidate.nodal_reduction[0].rcond_inf = 1e-15; }, [](auto& candidate) { candidate.backend.base_fpm[5U] = 1; }, - [](auto& candidate) { candidate.eta_gram = 1e-5; }}; + [](auto& candidate) { candidate.eta_gram = 1e-5; }, + [](auto& candidate) { candidate.drilling_mode.reset(); }}; for (const auto& mutate : mutations) { auto candidate = state; mutate(candidate); @@ -133,6 +142,27 @@ TEST(EigenHdf5ResultsWriter, WritesModalSchemaAndRejectsTamperedState) { locked_final.close(); ASSERT_TRUE( fesa::hdf5_internal::SelfCheckEigenFile(path, domain, state, {}).IsOk()); + for (const auto mode : {fesa::Mitc4DrillingMode::kDirectSpring, + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist}) { + state.drilling_mode = mode; + ASSERT_TRUE(writer.Write(path, domain, state, {}).IsOk()); + fesa::hdf5_internal::Hdf5Handle stored_file{ + H5Fopen(path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT), H5Fclose}; + const auto stored = fesa::hdf5_internal::ReadDrillingDiagnostics(stored_file.Get()); + ASSERT_TRUE(stored); + EXPECT_EQ(stored->formulation, fesa::Mitc4DrillingModeName(mode)); + EXPECT_TRUE(stored->frames.empty()); + } + { + fesa::hdf5_internal::Hdf5Handle stored_file{ + H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::Hdf5Handle metadata{ + H5Gopen2(stored_file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + } + EXPECT_FALSE(fesa::hdf5_internal::SelfCheckEigenFile(path, domain, state, {}).IsOk()); + ASSERT_TRUE(writer.Write(path, domain, state, {}).IsOk()); const hid_t corrupt_file = H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT); ASSERT_GE(corrupt_file, 0); @@ -153,3 +183,22 @@ TEST(EigenHdf5ResultsWriter, WritesModalSchemaAndRejectsTamperedState) { } } // namespace + + +TEST(EigenHdf5ResultsWriter, DrillingLifecycleFaultsPreserveFinalBytes) { + const auto path = std::filesystem::temp_directory_path() / "fesa-eigen-drilling-fault.h5"; + auto domain = fesa::Domain::Create(Definition(path.parent_path() / "shell.inp")); + ASSERT_TRUE(domain.HasValue()); + auto state = fesa::LinearFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(state.HasValue()); + for (const std::string fault : {"candidate-close", "reader-close", "markers"}) { + SCOPED_TRACE(fault); + std::ofstream(path, std::ios::binary | std::ios::trunc) << "prior"; + EXPECT_FALSE(fesa::hdf5_internal::WriteEigenFile(path, domain.Value(), state.Value(), {}, + fesa::test::Mitc4LifecycleFault(fault)).IsOk()); + std::ifstream stream{path, std::ios::binary}; + EXPECT_EQ((std::string{std::istreambuf_iterator{stream}, {}}), "prior"); + } + std::error_code ignored; + std::filesystem::remove(path, ignored); +} diff --git a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp index 5250cd4..aeb0ab9 100644 --- a/tests/unit/io/hdf5/hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/hdf5_results_writer_test.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -27,6 +28,13 @@ #include "fesa/math/vector3.h" #include "fesa/model/domain.h" +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "mitc4_fault_test_helpers.h" + +#include "io/hdf5/hdf5_result_writer.h" +#include "io/hdf5/hdf5_self_check.h" + namespace { constexpr const char* kStepRoot = "/steps/Step-1/frames/0"; @@ -241,7 +249,8 @@ fesa::ModelDefinition MakeShellDefinition(const std::filesystem::path& source) { return definition; } -WriterFixture MakeShellFixture(const std::filesystem::path& source) { +WriterFixture MakeShellFixture(const std::filesystem::path& source, + std::optional mode = fesa::Mitc4DrillingMode::kDirectSpring) { auto domain_result = fesa::Domain::Create(MakeShellDefinition(source)); if (!domain_result.HasValue()) { throw std::runtime_error{ @@ -303,6 +312,8 @@ WriterFixture MakeShellFixture(const std::filesystem::path& source) { candidate.physical_strain_energy = 123.5; candidate.equilibrium = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; candidate.verification_metrics = {1.0e-13, 2.0e-13, 3.0e-13}; + if (mode) candidate.mitc4_drilling = fesa::Mitc4DrillingEvidence{ + *mode, {0.5, 1.0, 1.5}, {1.0, 2.0, 3.0, 3.5, 4.0, 4.5}}; const fesa::Status commit = state->CommitShellResults({0U}, std::move(candidate)); if (!commit.IsOk()) { @@ -1343,3 +1354,149 @@ TEST(Hdf5ResultsWriter, EXPECT_EQ(ReadBytes(final), valid_bytes); EXPECT_EQ(EntryCount(directory.Path()), 1U); } + + +TEST(Hdf5ResultsWriter, DrillingActualModesRoundTripWithoutDefaultInference) { + for (const auto mode : {fesa::Mitc4DrillingMode::kDirectSpring, + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist}) { + TempDirectory directory{"drilling-roundtrip"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp", mode); + const auto output = directory.Path() / "results.h5"; + ASSERT_TRUE(fesa::Hdf5ResultsWriter{}.Write(output, *fixture.domain, *fixture.state, {}).IsOk()); + const auto file = OpenFile(output); + const auto stored = fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()); + ASSERT_TRUE(stored); + EXPECT_EQ(stored->formulation, fesa::Mitc4DrillingModeName(mode)); + ASSERT_EQ(stored->frames.size(), 1U); + EXPECT_EQ(stored->frames[0].raw_equilibrium, fixture.state->Equilibrium()); + EXPECT_EQ(stored->frames[0].drilling_moment, fixture.state->Mitc4Drilling()->drilling_moment); + EXPECT_EQ(stored->frames[0].equilibrium_with_drilling, + fixture.state->Mitc4Drilling()->equilibrium_with_drilling); + } +} + +TEST(Hdf5ResultsWriter, DrillingMissingActualStatePreservesFinal) { + TempDirectory directory{"drilling-missing"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp", std::nullopt); + const auto output = directory.Path() / "results.h5"; + const std::vector sentinel{'p','r','i','o','r'}; + WriteBytes(output, sentinel); + EXPECT_FALSE(fesa::Hdf5ResultsWriter{}.Write(output, *fixture.domain, *fixture.state, {}).IsOk()); + EXPECT_EQ(ReadBytes(output), sentinel); + EXPECT_EQ(EntryCount(directory.Path()), 1U); +} + +TEST(Hdf5ResultsWriter, DrillingReadbackCannotTreatRemovedMarkersAsHistorical) { + for (const bool remove_markers : {true, false}) { + TempDirectory directory{"drilling-readback"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp"); + const auto output = directory.Path() / "results.h5"; + ASSERT_TRUE(fesa::Hdf5ResultsWriter{}.Write(output, *fixture.domain, *fixture.state, {}).IsOk()); + { + fesa::hdf5_internal::Hdf5Handle file{ + H5Fopen(output.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + if (remove_markers) { + fesa::hdf5_internal::Hdf5Handle metadata{ + H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + EXPECT_FALSE(fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get())); + } else { + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/0/global/drilling_moment", H5P_DEFAULT), 0); + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/0/global/equilibrium_with_drilling", H5P_DEFAULT), 0); + } + } + const auto plan = fesa::hdf5_internal::BuildWriterPlan(output, *fixture.domain, *fixture.state, {}); + ASSERT_TRUE(plan.HasValue()); + EXPECT_FALSE(fesa::hdf5_internal::SelfCheck(output, *fixture.domain, *fixture.state, 0U, + plan.Value().model_data).IsOk()); + } +} + + +TEST(Hdf5ResultsWriter, DrillingReaderRejectsMalformedAndInconsistentEvidence) { + for (const std::string mutation : {"partial-mode", "partial-version", "mode", "version", "version-type", + "mode-type", "shape", "type", "endian", "nan", "inf", "algebra", "force"}) { + SCOPED_TRACE(mutation); + TempDirectory directory{"drilling-invalid"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp"); + const auto output = directory.Path() / "results.h5"; + ASSERT_TRUE(fesa::Hdf5ResultsWriter{}.Write(output, *fixture.domain, *fixture.state, {}).IsOk()); + fesa::hdf5_internal::Hdf5ErrorSilencer silence; + fesa::hdf5_internal::Hdf5Handle file{ + H5Fopen(output.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::Hdf5Handle metadata{ + H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + const std::string path = "/steps/Step-1/frames/0/global/drilling_moment"; + if (mutation == "partial-mode" || mutation == "partial-version") { + ASSERT_GE(H5Adelete(metadata.Get(), mutation == "partial-mode" ? + "mitc4_drilling_formulation" : "drilling_diagnostics_version"), 0); + } else if (mutation == "mode" || mutation == "mode-type") { + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + if (mutation == "mode") fesa::hdf5_internal::WriteStringAttribute(metadata.Get(), "mitc4_drilling_formulation", "unknown"); + else fesa::hdf5_internal::WriteUint64Attribute(metadata.Get(), "mitc4_drilling_formulation", 1U); + } else if (mutation == "version" || mutation == "version-type") { + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + if (mutation == "version") fesa::hdf5_internal::WriteUint64Attribute(metadata.Get(), "drilling_diagnostics_version", 2U); + else fesa::hdf5_internal::WriteDoubleAttribute(metadata.Get(), "drilling_diagnostics_version", 1.0); + } else if (mutation == "shape" || mutation == "type" || mutation == "endian") { + ASSERT_GE(H5Ldelete(file.Get(), path.c_str(), H5P_DEFAULT), 0); + const hsize_t size = mutation == "shape" ? 2U : 3U; + fesa::hdf5_internal::Hdf5Handle space{H5Screate_simple(1, &size, nullptr), H5Sclose}; + fesa::hdf5_internal::Hdf5Handle dataset{H5Dcreate2(file.Get(), path.c_str(), + mutation == "type" ? H5T_IEEE_F32LE : (mutation == "endian" ? H5T_IEEE_F64BE : H5T_IEEE_F64LE), + space.Get(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT), H5Dclose}; + ASSERT_GE(dataset.Get(), 0); + } else { + const std::string target = mutation == "inf" ? "/steps/Step-1/frames/0/global/equilibrium" : + (mutation == "nan" ? path : "/steps/Step-1/frames/0/global/equilibrium_with_drilling"); + fesa::hdf5_internal::Hdf5Handle dataset{H5Dopen2(file.Get(), target.c_str(), H5P_DEFAULT), H5Dclose}; + std::array values{1., 2., 3., 3.5, 4., 4.5}; + if (mutation == "nan") values[0] = std::numeric_limits::quiet_NaN(); + else if (mutation == "inf") values[0] = std::numeric_limits::infinity(); + else values[mutation == "force" ? 0U : 3U] = 1e-100; + ASSERT_GE(H5Dwrite(dataset.Get(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, values.data()), 0); + } + EXPECT_THROW((void)fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()), + fesa::hdf5_internal::Hdf5Failure); + } +} + + +TEST(Hdf5ResultsWriter, DrillingHistoricalSelfCheckDoesNotInferMode) { + TempDirectory directory{"drilling-historical"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp"); + const auto output = directory.Path() / "results.h5"; + ASSERT_TRUE(fesa::Hdf5ResultsWriter{}.Write(output, *fixture.domain, *fixture.state, {}).IsOk()); + { + fesa::hdf5_internal::Hdf5Handle file{ + H5Fopen(output.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + fesa::hdf5_internal::Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/0/global/drilling_moment", H5P_DEFAULT), 0); + ASSERT_GE(H5Ldelete(file.Get(), "/steps/Step-1/frames/0/global/equilibrium_with_drilling", H5P_DEFAULT), 0); + } + auto historical = MakeShellFixture(directory.Path() / "shell.inp", std::nullopt); + const auto plan = fesa::hdf5_internal::BuildWriterPlan(output, *historical.domain, *historical.state, {}); + ASSERT_TRUE(plan.HasValue()); + EXPECT_TRUE(fesa::hdf5_internal::SelfCheck(output, *historical.domain, *historical.state, 0U, + plan.Value().model_data).IsOk()); + EXPECT_FALSE(historical.state->Mitc4Drilling()); +} + + +TEST(Hdf5ResultsWriter, DrillingLifecycleFaultsPreserveFinalBytes) { + TempDirectory directory{"drilling-faults"}; + auto fixture = MakeShellFixture(directory.Path() / "shell.inp"); + const auto path = directory.Path() / "results.h5"; + const std::vector prior{'p','r','i','o','r'}; + for (const std::string fault : {"candidate-close", "reader-close", "markers", "datasets"}) { + SCOPED_TRACE(fault); + WriteBytes(path, prior); + EXPECT_FALSE(fesa::hdf5_internal::WriteStaticFile(path, *fixture.domain, *fixture.state, {}, + fesa::test::Mitc4LifecycleFault(fault)).IsOk()); + EXPECT_EQ(ReadBytes(path), prior); + EXPECT_EQ(EntryCount(directory.Path()), 1U); + } +} diff --git a/tests/unit/io/hdf5/mitc4_fault_test_helpers.h b/tests/unit/io/hdf5/mitc4_fault_test_helpers.h new file mode 100644 index 0000000..6479fd3 --- /dev/null +++ b/tests/unit/io/hdf5/mitc4_fault_test_helpers.h @@ -0,0 +1,42 @@ +#ifndef FESA_TESTS_UNIT_IO_HDF5_MITC4_FAULT_TEST_HELPERS_H_ +#define FESA_TESTS_UNIT_IO_HDF5_MITC4_FAULT_TEST_HELPERS_H_ + +#include + +#include "io/hdf5/mitc4_file_hooks.h" +#include "io/hdf5/hdf5_raii.h" + +namespace fesa::test { + +inline herr_t CloseMitc4FileAndReportFailure(hid_t file) { + (void)H5Fclose(file); + return -1; +} + +inline hdf5_internal::Mitc4FileHooks Mitc4LifecycleFault(const std::string& fault) { + hdf5_internal::Mitc4FileHooks hooks; + if (fault == "candidate-close") hooks.candidate_close = CloseMitc4FileAndReportFailure; + if (fault == "reader-close") hooks.reader_close = CloseMitc4FileAndReportFailure; + if (fault == "markers" || fault == "datasets") { + hooks.before_self_check = [fault](const std::filesystem::path& path) { + hdf5_internal::Hdf5Handle file{hdf5_internal::RequireHdf5Id( + H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), "Cannot corrupt test candidate."), H5Fclose}; + if (fault == "markers") { + hdf5_internal::Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + hdf5_internal::RequireHdf5(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), "Cannot remove mode."); + hdf5_internal::RequireHdf5(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), "Cannot remove version."); + } else { + hdf5_internal::RequireHdf5(H5Ldelete(file.Get(), + "/steps/Step-1/frames/0/global/drilling_moment", H5P_DEFAULT), "Cannot remove moment."); + hdf5_internal::RequireHdf5(H5Ldelete(file.Get(), + "/steps/Step-1/frames/0/global/equilibrium_with_drilling", H5P_DEFAULT), "Cannot remove corrected balance."); + } + hdf5_internal::RequireHdf5(file.CloseChecked(), "Cannot close corrupted candidate."); + }; + } + return hooks; +} + +} // namespace fesa::test + +#endif // FESA_TESTS_UNIT_IO_HDF5_MITC4_FAULT_TEST_HELPERS_H_ diff --git a/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp b/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp index ee8d73b..c357821 100644 --- a/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp +++ b/tests/unit/io/hdf5/nonlinear_hdf5_results_writer_test.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -25,6 +26,11 @@ #include "fesa/model/domain.h" #include "fesa/results/nonlinear_results_writer.h" #include "io/hdf5/hdf5_raii.h" +#include "io/hdf5/hdf5_self_check.h" + +#include "io/hdf5/mitc4_drilling_diagnostics.h" + +#include "mitc4_fault_test_helpers.h" namespace { @@ -184,7 +190,8 @@ std::vector MakeShellRows(const double seed) { fesa::AcceptedFrameState MakeFrame(const fesa::NonlinearTrialState& committed, const std::uint32_t frame_index, const std::uint32_t increment_index, - const double load_factor) { + const double load_factor, + std::optional mode = fesa::Mitc4DrillingMode::kDirectSpring) { fesa::AcceptedFrameState frame{ frame_index, increment_index, @@ -223,10 +230,13 @@ fesa::AcceptedFrameState MakeFrame(const fesa::NonlinearTrialState& committed, frame.reaction[offset + component] = frame.residual[offset + component]; } } + if (mode) frame.mitc4_drilling = fesa::Mitc4DrillingEvidence{ + *mode, {}, frame.global_equilibrium}; return frame; } -Fixture MakeFixture(const std::filesystem::path& source) { +Fixture MakeFixture(const std::filesystem::path& source, + std::optional mode = fesa::Mitc4DrillingMode::kDirectSpring) { auto domain_result = fesa::Domain::Create(MakeDefinition(source)); if (!domain_result.HasValue()) { throw std::runtime_error{"Unable to construct nonlinear HDF5 Domain."}; @@ -248,10 +258,10 @@ Fixture MakeFixture(const std::filesystem::path& source) { } fesa::NonlinearAnalysisState state = std::move(state_result.Value()); if (!state.CommitAcceptedFrame( - MakeFrame(state.CommittedState(), 0U, 1U, 0.5)) + MakeFrame(state.CommittedState(), 0U, 1U, 0.5, mode)) .IsOk() || !state.CommitAcceptedFrame( - MakeFrame(state.CommittedState(), 1U, 2U, 1.0)) + MakeFrame(state.CommittedState(), 1U, 2U, 1.0, mode)) .IsOk()) { throw std::runtime_error{"Unable to commit nonlinear HDF5 test frames."}; } @@ -615,3 +625,55 @@ TEST(NonlinearHdf5ResultsWriter, CandidateFinalizationFailurePreservesPriorFinal EXPECT_EQ(status.Diagnostics()[0U].code, "hdf5-finalization-failure"); EXPECT_EQ(ReadBytes(output), sentinel); } + + +TEST(NonlinearHdf5ResultsWriter, DrillingModesAndEveryAcceptedFrameRoundTrip) { + for (const auto mode : {fesa::Mitc4DrillingMode::kDirectSpring, + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist}) { + TempDirectory directory{"drilling"}; + auto fixture = MakeFixture(directory.Path() / "shell.inp", mode); + const auto path = directory.Path() / "results.h5"; + ASSERT_TRUE(fesa::NonlinearHdf5ResultsWriter{}.Write(path, *fixture.domain, fixture.state, {}).IsOk()); + { + const auto file = OpenFile(path); + const auto stored = fesa::hdf5_internal::ReadDrillingDiagnostics(file.Get()); + ASSERT_TRUE(stored); + EXPECT_EQ(stored->formulation, fesa::Mitc4DrillingModeName(mode)); + ASSERT_EQ(stored->frames.size(), fixture.state.AcceptedFrames().size()); + EXPECT_EQ(stored->frames.front().frame_index, 0U); + EXPECT_EQ(stored->frames.back().raw_equilibrium, fixture.state.AcceptedFrames().back().global_equilibrium); + } + { + Hdf5Handle file{H5Fopen(path.string().c_str(), H5F_ACC_RDWR, H5P_DEFAULT), H5Fclose}; + Hdf5Handle metadata{H5Gopen2(file.Get(), "/metadata", H5P_DEFAULT), H5Gclose}; + ASSERT_GE(H5Adelete(metadata.Get(), "mitc4_drilling_formulation"), 0); + ASSERT_GE(H5Adelete(metadata.Get(), "drilling_diagnostics_version"), 0); + } + EXPECT_FALSE(fesa::hdf5_internal::NonlinearSelfCheck(path, *fixture.domain, fixture.state, 0U).IsOk()); + } +} + +TEST(NonlinearHdf5ResultsWriter, DrillingMissingActualStatePreservesPriorFinal) { + TempDirectory directory{"drilling-missing"}; + auto fixture = MakeFixture(directory.Path() / "shell.inp", std::nullopt); + const auto path = directory.Path() / "results.h5"; + std::ofstream(path, std::ios::binary) << "prior"; + EXPECT_FALSE(fesa::NonlinearHdf5ResultsWriter{}.Write(path, *fixture.domain, fixture.state, {}).IsOk()); + const auto bytes = ReadBytes(path); + EXPECT_EQ(std::string(bytes.begin(), bytes.end()), "prior"); +} + + +TEST(NonlinearHdf5ResultsWriter, DrillingLifecycleFaultsPreserveFinalBytes) { + TempDirectory directory{"drilling-faults"}; + auto fixture = MakeFixture(directory.Path() / "shell.inp"); + const auto path = directory.Path() / "results.h5"; + for (const std::string fault : {"candidate-close", "reader-close", "markers", "datasets"}) { + SCOPED_TRACE(fault); + std::ofstream(path, std::ios::binary | std::ios::trunc) << "prior"; + EXPECT_FALSE(fesa::hdf5_internal::WriteNonlinearFile(path, *fixture.domain, fixture.state, {}, + fesa::test::Mitc4LifecycleFault(fault)).IsOk()); + const auto bytes = ReadBytes(path); + EXPECT_EQ(std::string(bytes.begin(), bytes.end()), "prior"); + } +} diff --git a/tests/unit/results/result_recovery_components_test.cpp b/tests/unit/results/result_recovery_components_test.cpp index 1f9a057..8c2785d 100644 --- a/tests/unit/results/result_recovery_components_test.cpp +++ b/tests/unit/results/result_recovery_components_test.cpp @@ -17,9 +17,12 @@ #include "fesa/assembly/parallel_for.h" #include "fesa/assembly/sparse_assembler.h" #include "fesa/elements/element.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/mitc4_shell.h" #include "fesa/fem/dof_manager.h" #include "fesa/math/vector3.h" #include "fesa/model/domain.h" +#include "fesa/results/result_recovery.h" #include "results/analysis_state_commit.h" #include "results/beam_result_recovery.h" #include "results/global_equilibrium_recovery.h" @@ -627,3 +630,99 @@ TEST(ResultRecoveryComponents, ExpectStatusCode(status, "nonfinite-recovery-value"); ExpectAnalysisStateEqual(state, prior); } + +TEST(Mitc4DrillingEvidence, LinearDirectSpringAcceptsKnownArtificialTorque) { + auto definition = MakeShellDefinition(); + for (std::size_t node = 0U; node < 4U; ++node) { + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 1, 5, 0.0, {definition.source_path, 51U}}); + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 6, 6, 0.125, {definition.source_path, 52U}}); + } + const auto fixture = MakeShellFixture(std::move(definition)); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; node < 4U; ++node) + state.Displacement()[6U * node + 5U] = 0.125; + const auto status = fesa::ResultRecovery::Recover( + *fixture.model, *fixture.dofs, *fixture.stiffness, state); + ASSERT_TRUE(status.IsOk()) << (status.Diagnostics().empty() ? "" : status.Diagnostics()[0U].message); +} + +TEST(Mitc4DrillingEvidence, LinearTorqueOracleAndSignedGateKeepRealImbalance) { + auto definition = MakeShellDefinition(); + for (std::size_t node = 0U; node < 4U; ++node) { + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 1, 5, 0.0, {definition.source_path, 51U}}); + definition.steps[0U].boundaries.push_back( + {std::to_string(node + 1U), 6, 6, 0.125, {definition.source_path, 52U}}); + } + const auto fixture = MakeShellFixture(std::move(definition)); + auto element = fesa::ElementFactory{}.Create(fixture.domain->Elements()[0U], + *fixture.domain); + ASSERT_TRUE(element.HasValue()); + const auto* shell = dynamic_cast(element.Value().get()); + ASSERT_NE(shell, nullptr); + const auto stiffness = shell->Stiffness(); + ASSERT_TRUE(stiffness.HasValue()); + double rotary_trace = 0.0; + for (std::size_t node = 0U; node < 4U; ++node) + for (std::size_t c = 3U; c < 5U; ++c) + rotary_trace += stiffness.Value().physical_local20(5U * node + c, 5U * node + c); + const double expected_torque = 4.0 * 0.125 * (1e-3 * rotary_trace / 8.0); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + for (std::size_t node = 0U; node < 4U; ++node) + state.Displacement()[6U * node + 5U] = 0.125; + ASSERT_TRUE(fesa::ResultRecovery::Recover(*fixture.model, *fixture.dofs, + *fixture.stiffness, state).IsOk()); + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_EQ(state.Mitc4Drilling()->mode, fesa::Mitc4DrillingMode::kDirectSpring); + EXPECT_NEAR(state.Equilibrium()[5U], expected_torque, 1e-14); + EXPECT_NEAR(state.Mitc4Drilling()->drilling_moment[2U], expected_torque, 1e-14); + EXPECT_DOUBLE_EQ(state.Mitc4Drilling()->drilling_moment[0U], 0.0); + EXPECT_DOUBLE_EQ(state.Mitc4Drilling()->drilling_moment[1U], 0.0); + for (double value : state.Mitc4Drilling()->equilibrium_with_drilling) + EXPECT_NEAR(value, 0.0, 1e-14); + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), 0.0); + + auto global = fesa::results_internal::GlobalEquilibriumRecovery( + *fixture.domain, *fixture.dofs, *fixture.stiffness, state.Displacement(), + state.ExternalForce()); + ASSERT_TRUE(global.HasValue()); + for (const int kind : {0, 1, 2, 3}) { + auto changed = global.Value(); + fesa::results_internal::RecoveryCandidate candidate{24U}; + candidate.shell.mitc4_drilling = *state.Mitc4Drilling(); + if (kind == 0) candidate.shell.mitc4_drilling->drilling_moment[2U] *= -1.0; + if (kind == 1) changed.residual[0U] += 1.0; + if (kind == 2) changed.residual[4U] += 1.0; + if (kind == 3) candidate.shell.mitc4_drilling->mode = + fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + const auto rejected = fesa::results_internal::RecoverShellGlobalEvidence( + *fixture.domain, *fixture.dofs, state.ExternalForce(), changed, candidate); + EXPECT_FALSE(rejected.IsOk()) << kind; + } +} + +TEST(Mitc4DrillingEvidence, StaticCommitRejectsInvalidEvidenceAtomically) { + const auto fixture = MakeShellFixture(MakeShellDefinition()); + auto state = fesa::AnalysisState::Create(*fixture.dofs, {"Step-1", 0U}); + fesa::ShellStateCandidate candidate{}; + candidate.rows = MakeShellRows(0U, 1.0).rows; + candidate.mitc4_drilling = fesa::Mitc4DrillingEvidence{ + fesa::Mitc4DrillingMode::kDirectSpring, {0.0, 0.0, 2.0}, {0, 0, 0, 0, 0, -2}}; + ASSERT_TRUE(state.CommitShellResults({0U}, candidate).IsOk()); + for (const int kind : {0, 1, 2, 3, 4}) { + auto invalid = candidate; + invalid.physical_strain_energy = 77.0; + if (kind == 0) invalid.mitc4_drilling->mode = static_cast(99); + if (kind == 1) invalid.mitc4_drilling->mode = fesa::Mitc4DrillingMode::kLegacyCorotatedTwist; + if (kind == 2) invalid.mitc4_drilling->drilling_moment[0U] = std::numeric_limits::infinity(); + if (kind == 3) invalid.mitc4_drilling->equilibrium_with_drilling[5U] = 2.0; + if (kind == 4) invalid.mitc4_drilling.reset(); + EXPECT_FALSE(state.CommitShellResults({0U}, invalid).IsOk()) << kind; + EXPECT_DOUBLE_EQ(state.PhysicalStrainEnergy(), 0.0); + ASSERT_TRUE(state.Mitc4Drilling()); + EXPECT_EQ(state.Mitc4Drilling()->mode, candidate.mitc4_drilling->mode); + EXPECT_EQ(state.Mitc4Drilling()->drilling_moment, candidate.mitc4_drilling->drilling_moment); + } +}