FESADev/docs/mitc4-shell-feast-eigen/numerical-review-evidence.md

# MITC4 독립 수치 검토 evidence

- Feature ID: `mitc4-shell-feast-eigen`
- Owner: `numerical-review-agent`
- Input contract revision: `r1`
- 기본 보고서: [numerical-review.md](numerical-review.md)
- 검토일 / web access date: `2026-09-06`

## NR-MFE-001 직접 계산 반례

관련 원문은 [requirements](requirements.md)의 `002`(기존 kernel 유지),
`008`(선택 질량의 유한 spectrum 보존), `014`(eigen 해석의 coordinate-invariance
evidence)와 [V5](formulation-verification.md#v5-rigid-mode-coordinate-및-density-불변식)다.
실제 [kernel](../../src/fesa/elements/mitc4_shell.cpp)의
`PhysicalTransformation20`(line 488), `DrillingTransformation4`(505),
`Stiffness`의 minimum tangent diagonal(633~653)을 읽기 전용으로 대조했다.

다음은 production 함수를 재호출한 expected value가 아니라 평판 다항식의 직접 적분이다.
Tangent axes `a=(0,1,0)`, `b=(-4/5,0,3/5)`, normal `d=(3/5,0,4/5)`를 둔다.
Node 순서는 `X1=-2a-b`, `X2=2a-b`, `X3=2a+b`, `X4=-2a+b`다.
`a×b=d`, 면적 8, 두께 1이고 constant director·positive Jacobian의 직사각형이다.
단일 element의 geometry-derived director와도 같다. `E=12, nu=0, rho=1`로 둔다.
모든 DOF를 고정하되 node1의 global `URZ` 하나만 free로 한다.

직사각형 `Lx=4,Ly=2`에서 `integral Nx² dA=1/6`,
`integral Ny² dA=2/3`, `|integral Nx Ny dA|=1/4`이다.
`D=Et³/[12(1-nu²)]=1`, `G=6`이며 MITC tying으로 각 nodal rotary
shear square integral은 `A/12`다. 따라서 각 node의 tangent rotary physical block은

\[
R=\begin{bmatrix}49/12&\pm1/8\\\pm1/8&23/6\end{bmatrix},
\quad j=\rho t^3A/48=1/6,\quad k_d=10^{-3}(23/6).
\]

Off-diagonal 부호는 corner의 `Nx Ny` 부호에 따라 달라지고 두 부호 모두 존재한다.
Global `eZ`의 tangent 성분은 `(0,3/5)`이고 drilling 성분은 `4/5`다.
Free mass는 `j(1-dZ²)=0.06>0`이므로 E2는 이 DOF를 삭제하지 않는다.
이 1차원 constrained pencil의 정확한 유한 eigenvalue는

\[
\lambda=\frac{(9/25)(23/6)+(16/25)k_d}{(9/25)(1/6)}
=23+(32/3)k_d.
\]

이제 전체 geometry를 global Z축으로 45° 회전한다. Free axis `eZ`와 모든 나머지
fixed subspace는 그대로이므로 동일 물리 BC가 계속 지원 component BC로 표현된다.
새 director는 `(3/(5 sqrt2),3/(5 sqrt2),4/5)`다. Least-aligned global-axis 규칙의
새 tangent와 회전 전 tangent 사이의 cosine/sine은 `-5/sqrt41,-4/sqrt41`이다.
모든 corner의 두 tangent diagonal 중 최소는 `937/246`으로 달라진다.
Physical `eZ^T Kphys eZ`와 질량은 불변이지만 `kd'=10^-3(937/246)`다.

| 값 | 회전 전 | 회전 후 |
| --- | ---: | ---: |
| `kd` | `0.0038333333333333336` | `0.003808943089430893` |
| `lambda` | `23.04088888888889` | `23.040628726287263` |

상대 차이는 `1.1291343961675893e-5`다. Tangent tie가 반대 axis를 선택해도
모든 corner/두 tangent diagonal의 minimum은 같다. 좌표변환으로 constraints를 바꾼
오류, input explicit director/MPC 기능, mass cutoff 또는 FEAST 오차를 이용한 반례가 아니다.
크기가 외부 tolerance보다 작다는 사실은 내부 coordinate-invariance 항등식의 근거가 아니다.

실행 명령은 아래 read-only here-script, `python -B -`, exit `0`이었다.

```powershell
@'
from math import sqrt
from fractions import Fraction as F
def dot(a,b): return sum(x*y for x,y in zip(a,b))
def cross(a,b): return (a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0])
def rot(v): return ((v[0]-v[1])/sqrt(2),(v[0]+v[1])/sqrt(2),v[2])
d=(.6,0,.8); a=(0.,1.,0.); b=cross(d,a)
Raa=F(49,12); Rbb=F(23,6); Rab=F(1,8); j=F(1,6)
k0=.001*min(Raa,Rbb)
dp=rot(d); ap0=(1-dp[0]**2,-dp[0]*dp[1],-dp[0]*dp[2])
ap=tuple(v/sqrt(dot(ap0,ap0)) for v in ap0)
u=dot(ap,rot(a)); v=dot(ap,rot(b))
diags=[float(Raa)*u*u+float(Rbb)*v*v+2*c*u*v for c in [-float(Rab),float(Rab)]]
diags += [float(Raa)*v*v+float(Rbb)*u*u-2*c*u*v for c in [-float(Rab),float(Rab)]]
k1=.001*min(diags)
l0=float(Rbb/j)+k0*.64/(float(j)*.36)
l1=float(Rbb/j)+k1*.64/(float(j)*.36)
print(float(k0),k1,l0,l1,abs(l1-l0)/l0)
assert abs(l1-l0)>1e-6
'@ | python -B -
```

이 계산은 기존 C++ kernel 실행 evidence가 아니다. 기존 승인 수식과 source의 계수
규칙이 갖는 mathematical counterexample이다. Runtime 회귀 test는 이후 Planning 소유다.

## 비충돌 수학 검토

| 대상 | 독립 검토 결과와 한계 |
| --- | --- |
| M1/M2 | `S=-[d]x`이므로 `S theta=theta×d`; global DOF `[u,theta]`와 20↔24 transformation의 부호가 일치한다. Shear factor는 mass에 들어가지 않는다. |
| M3~M5 | `C0/C1/C2` 차원은 `M/ML/ML²`; symmetric cross block은 transpose 관계다. Trilinear map의 J degree≤(2,2,2), integrand≤(4,4,4)이므로 3×3×3의 degree5 정확성이 충분하다. |
| Consistent positivity | Positive weights/J에서 quadratic form=0이면 모든 27점에서 Hq=0. Hq는 각 자연좌표 degree≤1이므로 다항식 자체가 0이다. z=0에서 bilinear 독립성으로 모든 ui=0, z 계수에서 `Si theta_i=0`. Null은 정확히 네 drilling이며 physical20은 PD다. |
| M6~M8 | Scalar moment row-sum→nodal block은 partition of unity와 일치한다. Thickness points가 서로 달라 strict Cauchy–Schwarz이고 central inertia>0, rank5다. Tiny nonzero director component를 0으로 만드는 근거는 없다. |
| M9/M10 | Jensen 차이의 1/2 double sum과 kinetic-energy의 추가 1/2를 구분한 식이 맞다. Total mass·spatial first moment·allocated thickness moments는 보존하지만 일반 orbital inertia나 consistent spectrum은 보존하지 않는다. |
| Geometry | 기존 `EvaluateGeometry`는 실제 J>0와 finite reciprocal/frame을 검사한다. 추가 27점에도 같은 predicate를 적용하면 그 점의 부적합을 거부한다. Stiffness 8점과 validation을 바꾸지 않는다. 기존 8점만 통과한 input의 mass failure는 감추지 않으며 전역 injectivity/warp 허용 확대를 주장하지 않는다. |
| E1/E2 | Length congruence의 단위는 맞다. 실제 global free space와 director null의 교집합으로 null dimension을 정하는 조건은 정확하다. Partial BC의 full-rank local mass는 허용되며 near-null은 conditioning failure와 구분해야 한다. |
| E3/E4 | 공통 director와 element별 positive nodal scalar kd로 `KZ=Z diag(kappa)`, `Kpz=0`, `Kzz>0`이 성립한다. Partial BC에 남는 kd를 제거하면 원 pencil을 바꾼다. 일반 Schur 예제의 5/2 대 단순 삭제의 3은 옳고 general dense fallback은 필요하지 않다. |
| E5/E8 | `Mp=LL^T`, `A=L^-1 Kp L^-T`, `phi=DCP L^-T x`는 exact congruence다. Identity B의 FEAST도 actual lumped 문제다. Mass normalization 단위, constrained exact zero 및 candidate commit 의미가 맞다. |
| V1~V4 | 원 assembled free `Kf y-lambda Mf y`를 복원 후 검산해야 null 복원 오류를 잡는다. Mixed-unit norm을 직접 섞지 않고 zero denominator branch를 정의한 식은 일관된다. Gram/energy/residual은 lowest-N completeness를 대신하지 않는다. |
| Rigid/density | 기존 tangent-projected rigid field는 continuum interpolation의 rigid motion을 재현한다. Massless drilling과 6 physical rigid modes는 다르다. Density scaling과 repeated-group projector 정의는 일관된다. Eigen coordinate suite의 적용 범위는 NR-MFE-001 해결 전 확정하지 않는다. |

## 독립 mass oracle

`rho=3,t=1/2`에서 아래 kernel-level geometry/director fixture를 사용했다.
이는 현재 input의 explicit normal 기능을 추가하는 제안이 아니다. Unit rational directors와
warped coordinates를 이용한 독립 moment 적분 사례다. 전체 patch input 생성이나 기존
geometry validator 실행을 대신하지 않는다.

```text
X = [(-1,-1,1/10), (1,-1,-1/10), (1,1,1/5), (-1,1,-1/10)]
d = [(0,0,1), (3/5,0,4/5), (2/3,1/3,2/3), (0,3/5,4/5)]
```

`fractions.Fraction`으로 coefficient dictionary `(xi power,eta power,zeta power)`를
만들고 M1의 trilinear map을 미분·determinant 전개했다. Monomial integral은 홀수 차수면
0, 모든 차수가 짝수면 `8/[(p+1)(q+1)(r+1)]`다. Production shape/J/mass 함수를
호출하지 않았다. 이 exact rational 결과를 별도의 3-point 및 2-point Gauss 평가와 대조했다.
실제 read-only here-script `python -B -`, exit `0`, wall time 약 `0.52 s`였다.

| 관찰 | 결과 |
| --- | --- |
| 실제 maximum integrand degree | `(4,4,4)` |
| 27점 J 범위 | `[0.12859919221615254, 0.27013414111718076]` |
| 3-point 최대 scalar entry absolute error | `1.1102230246251565e-16` |
| 2-point 최대 scalar entry absolute error | `0.0006018518518518257` |
| `C0_11, C1_11, C2_11` | `17609/28800`, `11003/1440000`, `88177/6912000` |
| `sum mi, sum hi, sum ji` | `7743/1600`, `739/12000`, `517/5120` |
| First moment, two independent rational expressions exactly equal | `(-79747/360000, -685951/2880000, 49651/320000)` |
| Minimum nodal `ji-hi²/mi` | `0.022318161338791092 > 0` |

Independent row-sum oracle는 다음과 같으며 C1의 강제 zero나 flat rotary 비례식을 잡는다.

```text
m = [298/225, 1759/1440, 6211/5760, 2333/1920]
h = [3191/192000, 1061/72000, 337/24000, 9323/576000]
j = [5969/216000, 8807/345600, 155527/6912000, 779/30720]
```

현재 수치값은 정식화의 독립 계산 evidence다. C++ test 통과, reference spectrum 통과 또는
양의 J가 전체 volume에서 보장된다는 실행 결과로 바꾸지 않는다.

## FEAST 독립 source 대조

Intel C 2026.0 URL edition의 [parameter 표](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2026-0/extended-eigensolver-input-parameters.html)를
2026-09-06 재확인했다. `fpm[49]=1`은 CSR interval의 total count를 `m0`로 반환하고 flag를
reset한다. Stochastic estimate라는 설명은 없다. 다른 standalone FEAST의 옵션 번호나
확률적 count 설명을 이 API에 대입하지 않는다. 값의 conditioning-independent 정확성은
문서가 보장하지 않으므로 이것만으로 implementation completeness를 승인하지 않는다.

[API의 epsout 정의](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2026-0/feast-scsrgv-feast-hcsrgv.html)는
trace 변화량을 endpoint absolute maximum으로 나눈다. 따라서 strict nonzero-width
finite interval이면 zero eigenvalue만 있다는 이유로 trace 분모가 0이 되지 않는다.
Size1, N=n, B=0 문제는 실제 backend test가 여전히 필요하다.
Count scratch와 solve m0를 분리한 E7은 output overwrite와 용량을 혼동하지 않는다.

[Output status 문서](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2026-0/extended-eigensolver-output-details.html)는
endpoint 근접 count failure와 넓은 구간 내 작은 cluster를 놓칠 수 있는 경우를 명시한다.
E6의 enclosure total n, monotone prefix/empty-prefix, safe bracket, multiplicity capacity,
returned count/Gram/full residual을 결합하는 방향은 이 한계를 숨기지 않는다.
Known-spectrum 검증 없는 `info=0` 또는 count만으로 lowest-N를 판정해서는 안 된다.
이 대조는 runtime 통과나 후속 numerical-control 값의 승인이 아니다.

## 재검토 범위

NR-MFE-001 때문에 bounded return을 먼저 수행한다. 아래 항목은 최종 numerical gate 전에
NR가 같은 승인 revision에서 구체적인 수치와 필수 case를 고정해야 하며 미완료로 남긴다.

- V1~V4의 symmetry/oracle/moment/centroid/Jensen/null/KZ/MZ/PZ/Cholesky/constraint/
  full residual/whitened residual/normalization/Gram/energy의 threshold와 zero branch.
- Positive block reciprocal-condition limit, thickness central-inertia cancellation,
  tiny nonzero director의 analytic rank 유지와 명시적 ill-conditioning failure.
- Operator-scale near-zero/negative/cluster classification, raw-value 보존과 endpoint ULP
  margin, contour/precision/refinement, finite count/search/endpoint/solve-retry budget.
- Deterministic initial integer generator/seed/QR sign·rank 규칙, thread/runtime boundary,
  n*m0/nnz/index overflow 및 count scratch/shifted factorization memory 검사.
- 필수 포트폴리오: affine M10와 위 curved moment, Jensen·rigid·density·source identity,
  비축정렬/부분 BC/near-null/null dimension, 일반 Schur의 오삭제 반례,
  실제 FEAST의 size1/N=n/zero/negative/repeated group>N/tight cluster/endpoint/
  insufficient subspace/count mismatch/Gram rank/overflow/allocation failure.
- 요구사항 정합 후 coordinate-invariance coverage와 NR-MFE-001 사례의 정확한 판정 역할.
  기존 K 수정, partial BC exclusion 또는 tolerance 완화로 임의 해결하지 않는다.

## 실행 evidence

다음 read-only 조사 명령은 exit `0`이었다. 원본 reference를 쓰거나 실행하지 않았다.

```powershell
Get-Content .codex/skills/fem-numerical-review/SKILL.md
rg -n 'PhysicalTransformation20|DrillingTransformation4|drilling_stiffness|EvaluateGeometry' src/fesa/elements/mitc4_shell.cpp
rg -n 'CreateShell|nodal_frames|directors' src/fesa/elements/element_factory.cpp src/fesa/model/shell_geometry.cpp
Get-Content 'reference/shellNL-eigen-meshfine/shellNL-eigen-meshfine.inp' -Tail 90
Get-Content 'reference/shellNL-eigen-meshfine/shellNL-eigen-meshfine eigenvalues.txt'
Get-FileHash 'reference/shellNL-eigen-meshfine/shellNL-eigen-meshfine.inp','reference/shellNL-eigen-meshfine/shellNL-eigen-meshfine eigenvalues.txt'
```

`python -B -` read-only inventory script도 exit `0`: `pathlib`와 keyword별 state로
nodes/elements/Set-2를 세고, 모든 nonempty TXT row에 fullmatch regex
`Frame:\s*Mode\s+(\d+):\s*Value\s*=\s*(\S+)\s+Freq\s*=\s*(\S+)\s*\(cycles/time\)`를
적용했다. `1..10` identity, finite positive Value, 104 unique boundary labels,
유일한 source step 및 input frequency data를 검사했다.
결과는 [reference-model](reference-model.md)의 inventory/hash와 같다.

필수 project/feature/skill reference 문서는 `Get-Content`와 scoped `rg`로 읽었다.
Web `open/find`로 위 공식 Intel 3개 page의 parameter/status/epsout을 확인했다.
Configured vault는 이번 자료의 readiness 조건이 아니므로 추가 조회를 요구하지 않았다.
External solver, production implementation, CMake/CTest, Harness, FEAST runtime 및
FESA-vs-reference 수치 비교는 단계 범위 밖이며 미실행이다.
Owner 경계 밖 파일·Hook state·reference artifact를 수정하거나 stage/commit하지 않았다.

Owned 3개 문서를 read-only `python -B -`로 metadata·상대 링크·trailing whitespace·
final newline·display-math delimiter 검사했다. Exit `0`, documents `3`, local links `35`,
issues `0`이다. 문서에 보존한 반례 코드를 추출해 같은 Python으로 재실행했고 exit `0`,
위 kd/lambda/relative difference가 재현되었다. `git diff --check` exit `0`; 타 owner의
`coordination.md` LF/CRLF warning만 출력되었다. 신규 owned 문서의 whitespace는 위 별도
검사로 확인했다. `git status --short`에서 owned 신규 3개와 Coordinator 변경을 구분했다.