# 3D Isoparametric Euler–Bernoulli Beam Formulation

## Metadata

- feature_id: `3d-isoparametric-euler-beam`
- source_requirement: `docs/linear-static-3d-euler-beam/requirements.md`
- source_research: `docs/linear-static-3d-euler-beam/research.md`
- approved_design: `docs/superpowers/specs/2026-08-08-linear-static-3d-euler-beam-design.md`
- status: `ready-for-numerical-review`
- owner_agent: `formulation-agent`
- date: `2026-08-09`

## 1. 목적과 정식화 요약

이 문서는 선형 정적 구조해석 파이프라인을 검증하기 위한 2절점 3차원
Euler–Bernoulli beam 요소의 수학 계약을 정의한다. 요소는 축 변형, Saint-Venant
비틀림, 두 주축 방향 굽힘을 포함하며 절점마다 3개의 병진 자유도와 3개의 회전
자유도를 가진다.

요소 계산의 핵심 흐름은 다음과 같다.

$$
\boxed{
\mathbf X_1,\mathbf X_2,\mathbf a
\rightarrow L,\mathbf R,\mathbf T
\rightarrow \mathbf B(\xi),\mathbf D
\rightarrow \mathbf K_e^l
\rightarrow \mathbf K_e^g
\rightarrow \mathbf K\mathbf d=\mathbf F
\rightarrow \boldsymbol\varepsilon,\mathbf s,\mathbf f_e
}
$$

행렬식으로 쓰면

$$
\begin{aligned}
\boldsymbol\varepsilon(\xi)&=\mathbf B(\xi)\mathbf d_e^l,\\
\mathbf s(\xi)&=\mathbf D\boldsymbol\varepsilon(\xi),\\
\mathbf K_e^l&=\int_0^L\mathbf B^T\mathbf D\mathbf B\,dx,\\
\mathbf K_e^g&=\mathbf T^T\mathbf K_e^l\mathbf T,\\
\mathbf r_e^g&=\mathbf K_e^g\mathbf d_e^g-\mathbf f_e^g.
\end{aligned}
$$

성분형은

$$
\varepsilon_i=\sum_{a=1}^{12}B_{ia}d_a^l,\qquad
s_i=\sum_{j=1}^{4}D_{ij}\varepsilon_j,
$$

$$
K_{ab}^{l}=\int_0^L\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}D_{ij}B_{jb}\,dx,\qquad
r_a^g=\sum_{b=1}^{12}K_{ab}^{g}d_b^g-f_a^g.
$$

[[Beam and Frame Finite Elements]]는 Euler–Bernoulli beam의 횡변위와 회전
자유도에 cubic Hermite 보간이 필요하고, 공간 frame 요소가 축·굽힘·비틀림과
좌표 변환을 결합한다는 이론적 근거를 제공한다. [[Isoparametric Finite Elements]]는
자연좌표, Jacobian, 미분 변환 및 수치적분의 공통 틀을 제공한다.

## 2. Scope and Assumptions

### 2.1 포함 범위

- 2절점 직선 prismatic 공간 beam 요소
- Abaqus `TYPE=B33`에만 연결되는 shear-free cubic Euler–Bernoulli 의미
- small displacement 및 small rotation
- 선형 정적 해석
- 균질 선형 탄성 재료
- 축 변형, 두 방향 Euler–Bernoulli 굽힘, Saint-Venant 비틀림
- 절점당 자유도 $[u,v,w,\theta_x,\theta_y,\theta_z]$
- 단면 도심을 지나는 국부 $x$축
- 단면 주축과 국부 $y,z$축의 일치: $I_{yz}=0$
- 보존력과 일관 등가절점하중
- formulation 검증용 상수 국부 line-load 등가절점벡터; Abaqus `*DLOAD` 입력은 제외
- 자연좌표 $\xi\in[-1,1]$을 사용하는 1차원 등매개 기하 mapping

### 2.2 제외 범위

- Timoshenko 전단 변형과 shear correction factor
- Abaqus `TYPE=B31`의 Euler–Bernoulli 요소 재해석
- 기하비선형, 재료비선형, P–$\Delta$, geometric stiffness
- 단면 warping, bimoment 및 warping 자유도
- shear-center/centroid offset, rigid offset 및 end release
- 곡선 beam, taper 또는 요소 내부에서 변하는 $E,G,A,I_y,I_z,J$
- 비주축 단면의 $I_{yz}\ne0$ 굽힘 결합
- 질량, 감쇠, 동적 해석
- 열변형, 초기변형률, 초기응력
- follower load

### 2.3 “Isoparametric”의 정확한 의미

이 요소에서 두 절점의 위치와 축변위 및 비틀림 회전은 같은 2절점 Lagrange
함수로 보간한다. Euler–Bernoulli 굽힘은 횡변위와 그 기울기를 동시에 절점에서
보간해야 하므로 cubic Hermite 함수를 사용한다. 따라서 이 문서에서
“isoparametric”은 **자연좌표 기반의 기하 mapping과 적분 틀**을 뜻하며, 여섯
운동학 변수 모두에 동일한 Lagrange 함수를 사용한다는 뜻은 아니다.

### 2.4 단위

일반 입력에는 특정 단위계를 고정하지 않고 일관 단위계를 사용한다. 승인된
`reference/cantilever beam/` 비교만 외부 계약에 따라 SI를 사용하며, 이 사실은
일반 `.inp` 입력의 단위계를 SI로 추론한다는 뜻이 아니다.

| 양 | 차원 |
| --- | --- |
| $x,L,u,v,w$ | length |
| $\theta_x,\theta_y,\theta_z$ | dimensionless (radian) |
| $E,G$ | force / length$^2$ |
| $A$ | length$^2$ |
| $I_y,I_z,J$ | length$^4$ |
| $N$ | force |
| $T,M_y,M_z$ | force $\times$ length |

## 3. 좌표계와 부호 규약

### 3.1 국부 직교좌표계

전역 절점 위치를 $\mathbf X_1,\mathbf X_2\in\mathbb R^3$라 하고 Abaqus first
section axis `n1`으로 제공되는 guide vector를 $\mathbf a=\mathbf n_1$이라 한다.
요소 tangent $\mathbf t$를 국부 $x$축으로, 투영·정규화한 `n1`을 국부 $y$축으로,
$\mathbf t\times\mathbf n_1$ 방향을 국부 $z$축으로 둔다.

$$
\begin{aligned}
\Delta\mathbf X&=\mathbf X_2-\mathbf X_1,\\
L&=\|\Delta\mathbf X\|,\\
\mathbf e_x&=\frac{\Delta\mathbf X}{L},\\
\widetilde{\mathbf e}_y&=\mathbf a-(\mathbf a\cdot\mathbf e_x)\mathbf e_x,\\
\mathbf e_y&=\frac{\widetilde{\mathbf e}_y}
{\|\widetilde{\mathbf e}_y\|},\\
\mathbf e_z&=\mathbf e_x\times\mathbf e_y.
\end{aligned}
$$

성분형은

$$
\Delta X_I=X_{2I}-X_{1I},\qquad
L=\sqrt{\sum_{I=1}^{3}(\Delta X_I)^2},\qquad
e_{xI}=\frac{\Delta X_I}{L},
$$

$$
\widetilde e_{yI}=a_I-\left(\sum_{J=1}^{3}a_Je_{xJ}\right)e_{xI},
\qquad
e_{yI}=\frac{\widetilde e_{yI}}
{\sqrt{\sum_{J=1}^{3}\widetilde e_{yJ}^2}},
$$

$$
e_{zI}=\sum_{J=1}^{3}\sum_{K=1}^{3}\epsilon_{IJK}e_{xJ}e_{yK}.
$$

좌표 크기와 guide vector 크기를 반영한 승인된 validation은 다음과 같다.

$$
L>10^{-12}\max(1,\|\mathbf X_1\|,\|\mathbf X_2\|),\qquad
\|\widetilde{\mathbf e}_y\|>
10^{-12}\max(1,\|\mathbf a\|),\qquad
\mathbf e_i\cdot\mathbf e_j=\delta_{ij},\qquad
\mathbf e_x\times\mathbf e_y=\mathbf e_z.
$$

zero guide vector, guide vector가 요소 축과 평행한 경우, 또는 zero/near-zero length는
국부 $y,z$축과 $I_y,I_z$의 의미를 정할 수 없으므로 model validation 오류로 처리한다.
임의의 fallback 축을 조용히 선택하지 않는다.

### 3.2 변위와 회전 부호

- $u,v,w$: 각각 국부 $+x,+y,+z$ 방향 병진변위
- $\theta_x,\theta_y,\theta_z$: 각각 국부축에 대한 오른손 법칙의 양의 회전
- Euler–Bernoulli 무전단 조건:

$$
\boxed{\theta_z=\frac{dv}{dx},\qquad
\theta_y=-\frac{dw}{dx}}
$$

두 식의 부호가 다르다는 점이 $w$–$\theta_y$ 굽힘 강성행렬의 $6EI/L^2$
항 부호를 결정한다.

### 3.3 단면력 부호

단면 일반화 변형률과 단면력을 다음 순서로 둔다.

$$
\boldsymbol\varepsilon=
\begin{bmatrix}\varepsilon_0&\kappa_x&\kappa_y&\kappa_z\end{bmatrix}^T,
\qquad
\mathbf s=
\begin{bmatrix}N&T&M_y&M_z\end{bmatrix}^T.
$$

$N$은 인장일 때 양수이고 $T,M_y,M_z$는 해당 국부축 오른손 법칙을 따른다.

## 4. Primary Variables and DOFs

### 4.1 연속장

약형의 독립장은 중심선 축변위 $u(x)$, 비틀림 회전 $\theta_x(x)$, 두 횡변위
$v(x),w(x)$이다. 단면 굽힘 회전은 횡변위 기울기에서 유도한다.

$$
u,\theta_x\in H^1(0,L),\qquad v,w\in H^2(0,L).
$$

### 4.2 절점 자유도 벡터

국부 요소 자유도 순서는 다음으로 고정한다.

$$
\boxed{
\mathbf d_e^l=
\begin{bmatrix}
u_1&v_1&w_1&\theta_{x1}&\theta_{y1}&\theta_{z1}&
u_2&v_2&w_2&\theta_{x2}&\theta_{y2}&\theta_{z2}
\end{bmatrix}^T
}
$$

성분 번호는 다음과 같다.

| $a$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $d_a^l$ | $u_1$ | $v_1$ | $w_1$ | $\theta_{x1}$ | $\theta_{y1}$ | $\theta_{z1}$ | $u_2$ | $v_2$ | $w_2$ | $\theta_{x2}$ | $\theta_{y2}$ | $\theta_{z2}$ |

전역 nodal component 순서는 각 node에서 `[UX, UY, UZ, URX, URY, URZ]`이고,
$\mathbf d_e^l=\mathbf T\mathbf d_e^g$는 이 node/component block 순서를 유지한 채
국부 $[u,v,w,\theta_x,\theta_y,\theta_z]$로 회전한다.

요소는 equation id를 소유하지 않는다. 전역 equation numbering과 constrained/free
mapping은 프로젝트 아키텍처에 따라 DofManager가 소유해야 한다.

## 5. Strong Form and Boundary Conditions

### 5.1 단면력–변위 관계

상수 단면과 선형 탄성에서

$$
\begin{aligned}
N&=EA\,u',\\
T&=GJ\,\theta_x',\\
M_y&=EI_y\,\theta_y'=-EI_y\,w'',\\
M_z&=EI_z\,\theta_z'=EI_z\,v''.
\end{aligned}
$$

여기서 $(\cdot)'=d(\cdot)/dx$이다.

### 5.2 미분방정식

국부 분포하중을 $p_x,p_y,p_z$ 및 단위길이당 분포 비틀림 모멘트를 $m_x$라
하면 요소 내부 평형식은

$$
\boxed{
\begin{aligned}
-N'&=p_x,\\
-T'&=m_x,\\
M_z''&=p_y,\\
-M_y''&=p_z.
\end{aligned}}
$$

변위로 직접 쓰면

$$
\boxed{
\begin{aligned}
-(EAu')'&=p_x,\\
-(GJ\theta_x')'&=m_x,\\
(EI_zv'')''&=p_y,\\
(EI_yw'')''&=p_z.
\end{aligned}}
$$

상수 물성이고 요소 내부 분포하중이 없으면

$$
u''=0,\qquad \theta_x''=0,\qquad v''''=0,\qquad w''''=0.
$$

따라서 선형 축/비틀림 보간과 cubic Hermite 굽힘 보간이 각각 homogeneous
strong form을 정확히 만족한다.

### 5.3 경계조건과 자연경계 부호

경계 외향 법선의 1차원 부호를 $n=-1$ at $x=0$, $n=+1$ at $x=L$로 둔다.
각 끝에서는 다음 essential/natural 쌍 중 하나를 지정한다.

| 장 | essential | natural |
| --- | --- | --- |
| 축 | $u=\bar u$ | $F_x=nN$ |
| 비틀림 | $\theta_x=\bar\theta_x$ | $M_x=nT$ |
| $x$–$y$ 평면 굽힘 | $v=\bar v$ | $F_y=-nM_z'$ |
| $x$–$y$ 평면 회전 | $\theta_z=\bar\theta_z$ | $M_z^{\mathrm{end}}=nM_z$ |
| $x$–$z$ 평면 굽힘 | $w=\bar w$ | $F_z=nM_y'$ |
| $x$–$z$ 평면 회전 | $\theta_y=\bar\theta_y$ | $M_y^{\mathrm{end}}=nM_y$ |

굽힘 방정식은 4차이므로 각 굽힘 평면에서 변위/회전 또는 전단력/모멘트의
두 경계조건이 필요하다.

## 6. Weak or Variational Form

허용 시험함수 $\delta u,\delta\theta_x,\delta v,\delta w$는 essential
boundary에서 0이다. 두 번의 부분적분을 적용하면 다음 약형을 얻는다.

$$
\boxed{
\begin{aligned}
&\int_0^L
\left[
EA\,\delta u' u'
+GJ\,\delta\theta_x'\theta_x'
+EI_y\,\delta w''w''
+EI_z\,\delta v''v''
\right]dx\\
&=\int_0^L
\left[
\delta u\,p_x
+\delta\theta_x\,m_x
+\delta v\,p_y
+\delta w\,p_z
\right]dx
+\delta\mathbf d_{\Gamma}^T\bar{\mathbf f}_{\Gamma}.
\end{aligned}}
$$

경계 가상일의 성분형은

$$
\delta\mathbf d_{\Gamma}^T\bar{\mathbf f}_{\Gamma}
=\sum_{x\in\{0,L\}}
\left(
\delta u F_x+\delta v F_y+\delta w F_z
+\delta\theta_xM_x+\delta\theta_yM_y^{\mathrm{end}}
+\delta\theta_zM_z^{\mathrm{end}}
\right).
$$

일반화 변형률과 단면력을 사용하면 내부 가상일은 더 간단히

$$
\delta W_{\mathrm{int}}^e
=\int_0^L\delta\boldsymbol\varepsilon^T\mathbf s\,dx
=\int_0^L\delta\boldsymbol\varepsilon^T
\mathbf D\boldsymbol\varepsilon\,dx
$$

이고, 성분형은

$$
\delta W_{\mathrm{int}}^e
=\int_0^L
\left(
N\,\delta\varepsilon_0
+T\,\delta\kappa_x
+M_y\,\delta\kappa_y
+M_z\,\delta\kappa_z
\right)dx.
$$

선형 탄성에서는 동일한 식을 전체 잠재에너지의 정류조건으로 쓸 수 있다.

$$
\Pi_e(\mathbf d_e^l)
=\frac12(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l
-(\mathbf d_e^l)^T\mathbf f_e^l,
\qquad
\frac{\partial\Pi_e}{\partial\mathbf d_e^l}=\mathbf0.
$$

## 7. Mapping and Natural Coordinates

### 7.1 2절점 등매개 기하 mapping

자연좌표 $\xi\in[-1,1]$에서

$$
N_1^L(\xi)=\frac{1-\xi}{2},\qquad
N_2^L(\xi)=\frac{1+\xi}{2}.
$$

기하 위치는

$$
\mathbf X(\xi)=N_1^L(\xi)\mathbf X_1+N_2^L(\xi)\mathbf X_2.
$$

국부 축좌표 $x\in[0,L]$는

$$
x(\xi)=\frac{L}{2}(1+\xi),\qquad
r=\frac{x}{L}=\frac{1+\xi}{2}.
$$

Jacobian과 미분 변환은

$$
J_x=\frac{dx}{d\xi}=\frac{L}{2}>0,
\qquad
dx=J_xd\xi,
$$

$$
\frac{d}{dx}=\frac{1}{J_x}\frac{d}{d\xi}
=\frac{2}{L}\frac{d}{d\xi},\qquad
\frac{d^2}{dx^2}=\frac{4}{L^2}\frac{d^2}{d\xi^2}.
$$

직선 2절점 요소에서는 $J_x$가 상수이다. $J_x\le J_{\mathrm{tol}}$이면 zero
또는 near-zero length 요소이므로 계산을 중단한다.

## 8. Discretization

### 8.1 축변위와 비틀림 보간

$$
u(\xi)=N_1^Lu_1+N_2^Lu_2,
\qquad
\theta_x(\xi)=N_1^L\theta_{x1}+N_2^L\theta_{x2}.
$$

성분형 보간 조건은

$$
N_i^L(\xi_j)=\delta_{ij},\qquad
N_1^L+N_2^L=1.
$$

### 8.2 Cubic Hermite 함수

$r=x/L$을 사용하면

$$
\begin{aligned}
H_1(r)&=1-3r^2+2r^3,\\
H_2(r)&=L(r-2r^2+r^3),\\
H_3(r)&=3r^2-2r^3,\\
H_4(r)&=L(-r^2+r^3).
\end{aligned}
$$

자연좌표 성분형은

$$
\begin{aligned}
H_1(\xi)&=\frac14(2-3\xi+\xi^3),\\
H_2(\xi)&=\frac{L}{8}(1-\xi-\xi^2+\xi^3),\\
H_3(\xi)&=\frac14(2+3\xi-\xi^3),\\
H_4(\xi)&=\frac{L}{8}(-1-\xi+\xi^2+\xi^3).
\end{aligned}
$$

횡변위 보간은 회전 부호를 반영하여

$$
\boxed{
\begin{aligned}
v(x)&=H_1v_1+H_2\theta_{z1}+H_3v_2+H_4\theta_{z2},\\
w(x)&=H_1w_1-H_2\theta_{y1}+H_3w_2-H_4\theta_{y2}.
\end{aligned}}
$$

1차 미분은

$$
\begin{aligned}
H_1'&=\frac{6(r^2-r)}{L},&
H_2'&=1-4r+3r^2,\\
H_3'&=\frac{6(r-r^2)}{L},&
H_4'&=-2r+3r^2.
\end{aligned}
$$

2차 미분은

$$
\begin{aligned}
H_1''&=\frac{-6+12r}{L^2},&
H_2''&=\frac{-4+6r}{L},\\
H_3''&=\frac{6-12r}{L^2},&
H_4''&=\frac{-2+6r}{L}.
\end{aligned}
$$

### 8.3 Kronecker 및 기울기 조건

$$
\begin{array}{c|cccc|cccc}
&H_1&H_2&H_3&H_4&H_1'&H_2'&H_3'&H_4'\\ \hline
x=0&1&0&0&0&0&1&0&0\\
x=L&0&0&1&0&0&0&0&1
\end{array}
$$

따라서

$$
v(0)=v_1,\quad v'(0)=\theta_{z1},\quad
v(L)=v_2,\quad v'(L)=\theta_{z2},
$$

$$
w(0)=w_1,\quad w'(0)=-\theta_{y1},\quad
w(L)=w_2,\quad w'(L)=-\theta_{y2}.
$$

또한 $H_1+H_3=1$이다. Hermite 보간의 회전 함수 $H_2,H_4$는 길이
차원을 가지므로 네 함수를 단순히 모두 더하는 partition-of-unity 검사는 의미가
없다. 대신 상수 병진과 선형 rigid rotation field의 정확 재현을 검사한다.

### 8.4 6성분 운동학 보간행렬

중심선 운동학 벡터를

$$
\mathbf q(x)=
\begin{bmatrix}u&v&w&\theta_x&\theta_y&\theta_z\end{bmatrix}^T
=\mathbf N_q(x)\mathbf d_e^l
$$

로 두면

$$
\mathbf N_q=
\begin{bmatrix}
N_1^L&0&0&0&0&0&N_2^L&0&0&0&0&0\\
0&H_1&0&0&0&H_2&0&H_3&0&0&0&H_4\\
0&0&H_1&0&-H_2&0&0&0&H_3&0&-H_4&0\\
0&0&0&N_1^L&0&0&0&0&0&N_2^L&0&0\\
0&0&-H_1'&0&H_2'&0&0&0&-H_3'&0&H_4'&0\\
0&H_1'&0&0&0&H_2'&0&H_3'&0&0&0&H_4'
\end{bmatrix}_{6\times12}.
$$

성분형은

$$
q_i(x)=\sum_{a=1}^{12}N_{q,ia}(x)d_a^l.
$$

## 9. Kinematics

### 9.1 3차원 단면 운동학

도심에서 국부 단면좌표가 $(y,z)$인 점의 강체 병진·회전 부분은

$$
\begin{bmatrix}U_x\\U_y\\U_z\end{bmatrix}
=
\begin{bmatrix}u\\v\\w\end{bmatrix}
+
\begin{bmatrix}\theta_x\\\theta_y\\\theta_z\end{bmatrix}
\times
\begin{bmatrix}0\\y\\z\end{bmatrix}
=
\begin{bmatrix}
u+z\theta_y-y\theta_z\\
v-z\theta_x\\
w+y\theta_x
\end{bmatrix}.
$$

Euler–Bernoulli 조건을 대입한 축방향 변위와 변형률은

$$
U_x=u-zw'-yv',
$$

$$
\boxed{
\varepsilon_{xx}(x,y,z)
=\frac{\partial U_x}{\partial x}
=\varepsilon_0+z\kappa_y-y\kappa_z
}
$$

이다.

위 식은 이 요소가 보유하는 6개 단면운동 자유도만 나타내며 비원형 단면의
Saint-Venant warping 변위장은 포함하지 않는다. 축 및 굽힘 변형률은 이 운동학에서
유도하고, 비틀림은 별도의 1차원 에너지/단면력 관계 $T=GJ\kappa_x$로 표현한다.

### 9.2 일반화 변형률

$$
\boxed{
\boldsymbol\varepsilon=
\begin{bmatrix}
\varepsilon_0\\\kappa_x\\\kappa_y\\\kappa_z
\end{bmatrix}
=
\begin{bmatrix}
u'\\\theta_x'\\\theta_y'\\\theta_z'
\end{bmatrix}
=
\begin{bmatrix}
u'\\\theta_x'\\-w''\\v''
\end{bmatrix}}
$$

### 9.3 Strain–displacement matrix

$$
\boldsymbol\varepsilon(x)=\mathbf B(x)\mathbf d_e^l
$$

이며 $r=x/L$에서 성분행렬은

$$
\boxed{
\mathbf B=
\begin{bmatrix}
-\frac1L&0&0&0&0&0&\frac1L&0&0&0&0&0\\
0&0&0&-\frac1L&0&0&0&0&0&\frac1L&0&0\\
0&0&\frac{6-12r}{L^2}&0&\frac{-4+6r}{L}&0&
0&0&\frac{-6+12r}{L^2}&0&\frac{-2+6r}{L}&0\\
0&\frac{-6+12r}{L^2}&0&0&0&\frac{-4+6r}{L}&
0&\frac{6-12r}{L^2}&0&0&0&\frac{-2+6r}{L}
\end{bmatrix}_{4\times12}}
$$

즉 비영 성분은

$$
\begin{aligned}
\varepsilon_0&=\frac{u_2-u_1}{L},\\
\kappa_x&=\frac{\theta_{x2}-\theta_{x1}}{L},\\
\kappa_y&=\frac{6-12r}{L^2}w_1
+\frac{-4+6r}{L}\theta_{y1}
+\frac{-6+12r}{L^2}w_2
+\frac{-2+6r}{L}\theta_{y2},\\
\kappa_z&=\frac{-6+12r}{L^2}v_1
+\frac{-4+6r}{L}\theta_{z1}
+\frac{6-12r}{L^2}v_2
+\frac{-2+6r}{L}\theta_{z2}.
\end{aligned}
$$

## 10. Constitutive Contract

승인된 `*BEAM GENERAL SECTION, SECTION=GENERAL` 수학 mapping은 입력 순서
`A, I11, I12, I22, J`에 대해

$$
\boxed{I_y=I11,\qquad I_z=I22,\qquad I12=0}
$$

이다. 여기서 section axis 1은 `n1 -> local y`, section axis 2는
`t x n1 -> local z`이므로 $I11$은 local $y$축 굽힘의 $I_y$, $I22$는 local
$z$축 굽힘의 $I_z$에 대응한다. $I12\ne0$을 대각 행렬로 근사하지 않는다.

### 10.1 단면 적분

도심축 조건은

$$
\int_A y\,dA=0,\qquad \int_A z\,dA=0
$$

이고 주축 조건은

$$
I_{yz}=\int_A yz\,dA=0.
$$

단면 특성은

$$
A=\int_A dA,\qquad
I_y=\int_A z^2dA,\qquad
I_z=\int_A y^2dA.
$$

축응력과 단면력 정의는

$$
\sigma_{xx}=E\varepsilon_{xx},
$$

$$
N=\int_A\sigma_{xx}dA,\qquad
M_y=\int_Az\sigma_{xx}dA,\qquad
M_z=-\int_Ay\sigma_{xx}dA.
$$

이를 적분하면

$$
N=EA\varepsilon_0,\qquad
M_y=EI_y\kappa_y,\qquad
M_z=EI_z\kappa_z.
$$

Saint-Venant 비틀림은

$$
T=GJ\kappa_x
$$

로 둔다. 여기서 $J$는 polar moment와 항상 같은 값이 아니라, 일반 단면에 대한
Saint-Venant torsional constant이다.

### 10.2 일반화 탄성행렬

$$
\boxed{
\mathbf s=\mathbf D\boldsymbol\varepsilon,
\qquad
\mathbf D=
\begin{bmatrix}
EA&0&0&0\\
0&GJ&0&0\\
0&0&EI_y&0\\
0&0&0&EI_z
\end{bmatrix}_{4\times4}}
$$

성분형은

$$
s_i=\sum_{j=1}^{4}D_{ij}\varepsilon_j,
$$

$$
N=EA\varepsilon_0,\quad
T=GJ\kappa_x,\quad
M_y=EI_y\kappa_y,\quad
M_z=EI_z\kappa_z.
$$

필수 유효성 조건은

$$
E>0,\quad G>0,\quad A>0,\quad J>0,\quad I_y>0,\quad I_z>0.
$$

V0의 균질 등방성 재료는 입력의 정확한 $E,\nu$로부터

$$
G=\frac{E}{2(1+\nu)}
$$

를 계산하고 위의 $E>0$, $G>0$ 조건을 적용한다. 독립 $G$ 입력은 승인된 V0 입력
계약이 아니다.

### 10.3 비주축 단면에 대한 참고식

현재 범위 밖이지만 $I_{yz}\ne0$이면 굽힘 블록은

$$
\begin{bmatrix}M_y\\M_z\end{bmatrix}
=E
\begin{bmatrix}
I_y&-I_{yz}\\
-I_{yz}&I_z
\end{bmatrix}
\begin{bmatrix}\kappa_y\\\kappa_z\end{bmatrix}
$$

가 된다. 현재 구현 계약에서 이 항을 0으로 조용히 버리면 안 되며 입력을 거부하거나
향후 별도 기능으로 확장해야 한다.

## 11. Element Equations

### 11.1 내부력, 잔차 및 강성

$$
\mathbf f_{\mathrm{int},e}^l
=\int_0^L\mathbf B^T\mathbf s\,dx
=\int_0^L\mathbf B^T\mathbf D\mathbf B\,dx\,\mathbf d_e^l
=\mathbf K_e^l\mathbf d_e^l.
$$

성분형은

$$
f_{\mathrm{int},a}^l
=\int_0^L\sum_{i=1}^{4}B_{ia}s_i\,dx,
$$

$$
K_{ab}^l
=\int_0^L\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}D_{ij}B_{jb}\,dx.
$$

요소 잔차는

$$
\boxed{\mathbf r_e^l=\mathbf K_e^l\mathbf d_e^l-\mathbf f_{\mathrm{ext},e}^l}
$$

이다. 선형 문제이므로 tangent와 stiffness가 같고 변위와 무관하다.

### 11.2 부분 강성행렬

축 자유도 $[u_1,u_2]$:

$$
\mathbf K_a=\frac{EA}{L}
\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.
$$

비틀림 자유도 $[\theta_{x1},\theta_{x2}]$:

$$
\mathbf K_t=\frac{GJ}{L}
\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.
$$

$v$–$\theta_z$ 자유도 $[v_1,\theta_{z1},v_2,\theta_{z2}]$:

$$
\mathbf K_{bz}=\frac{EI_z}{L^3}
\begin{bmatrix}
12&6L&-12&6L\\
6L&4L^2&-6L&2L^2\\
-12&-6L&12&-6L\\
6L&2L^2&-6L&4L^2
\end{bmatrix}.
$$

$w$–$\theta_y$ 자유도 $[w_1,\theta_{y1},w_2,\theta_{y2}]$:

$$
\mathbf K_{by}=\frac{EI_y}{L^3}
\begin{bmatrix}
12&-6L&-12&-6L\\
-6L&4L^2&6L&2L^2\\
-12&6L&12&6L\\
-6L&2L^2&6L&4L^2
\end{bmatrix}.
$$

대표 적분 성분은

$$
K_{v_1v_1}=EI_z\int_0^L(H_1'')^2dx=\frac{12EI_z}{L^3},
$$

$$
K_{v_1\theta_{z1}}=EI_z\int_0^LH_1''H_2''dx
=\frac{6EI_z}{L^2},
$$

$$
K_{w_1\theta_{y1}}=EI_y\int_0^L(-H_1'')H_2''dx
=-\frac{6EI_y}{L^2}.
$$

마지막 식의 음수가 $\theta_y=-w'$에서 직접 발생한다.

### 11.3 전체 12×12 국부 강성행렬

다음 약어를 둔다.

$$
\begin{gathered}
a=\frac{EA}{L},\qquad t=\frac{GJ}{L},\\
b_y=\frac{12EI_y}{L^3},\quad c_y=\frac{6EI_y}{L^2},
\quad d_y=\frac{4EI_y}{L},\quad e_y=\frac{2EI_y}{L},\\
b_z=\frac{12EI_z}{L^3},\quad c_z=\frac{6EI_z}{L^2},
\quad d_z=\frac{4EI_z}{L},\quad e_z=\frac{2EI_z}{L}.
\end{gathered}
$$

그러면 4.2절의 자유도 순서에 대한 국부 강성행렬은

$$
\boxed{
\mathbf K_e^l=
\begin{bmatrix}
a&0&0&0&0&0&-a&0&0&0&0&0\\
0&b_z&0&0&0&c_z&0&-b_z&0&0&0&c_z\\
0&0&b_y&0&-c_y&0&0&0&-b_y&0&-c_y&0\\
0&0&0&t&0&0&0&0&0&-t&0&0\\
0&0&-c_y&0&d_y&0&0&0&c_y&0&e_y&0\\
0&c_z&0&0&0&d_z&0&-c_z&0&0&0&e_z\\
-a&0&0&0&0&0&a&0&0&0&0&0\\
0&-b_z&0&0&0&-c_z&0&b_z&0&0&0&-c_z\\
0&0&-b_y&0&c_y&0&0&0&b_y&0&c_y&0\\
0&0&0&-t&0&0&0&0&0&t&0&0\\
0&0&-c_y&0&e_y&0&0&0&c_y&0&d_y&0\\
0&c_z&0&0&0&e_z&0&-c_z&0&0&0&d_z
\end{bmatrix}}
$$

이어야 한다.

### 11.4 국부 요소 끝력의 성분형

분포하중이 없는 경우

$$
\mathbf f_{\mathrm{int},e}^l=
\begin{bmatrix}
F_{x1}&F_{y1}&F_{z1}&M_{x1}&M_{y1}&M_{z1}&
F_{x2}&F_{y2}&F_{z2}&M_{x2}&M_{y2}&M_{z2}
\end{bmatrix}^T
=\mathbf K_e^l\mathbf d_e^l.
$$

각 성분은

$$
\begin{aligned}
F_{x1}&=a(u_1-u_2),\\
F_{y1}&=b_z(v_1-v_2)+c_z(\theta_{z1}+\theta_{z2}),\\
F_{z1}&=b_y(w_1-w_2)-c_y(\theta_{y1}+\theta_{y2}),\\
M_{x1}&=t(\theta_{x1}-\theta_{x2}),\\
M_{y1}&=-c_yw_1+d_y\theta_{y1}+c_yw_2+e_y\theta_{y2},\\
M_{z1}&=c_zv_1+d_z\theta_{z1}-c_zv_2+e_z\theta_{z2},\\[2mm]
F_{x2}&=-a(u_1-u_2),\\
F_{y2}&=-b_z(v_1-v_2)-c_z(\theta_{z1}+\theta_{z2}),\\
F_{z2}&=-b_y(w_1-w_2)+c_y(\theta_{y1}+\theta_{y2}),\\
M_{x2}&=-t(\theta_{x1}-\theta_{x2}),\\
M_{y2}&=-c_yw_1+e_y\theta_{y1}+c_yw_2+d_y\theta_{y2},\\
M_{z2}&=c_zv_1+e_z\theta_{z1}-c_zv_2+d_z\theta_{z2}.
\end{aligned}
$$

### 11.5 차원과 기본 행렬 특성

$$
\mathbf B:4\times12,\quad
\mathbf D:4\times4,\quad
\mathbf K_e^l:12\times12,\quad
\mathbf d_e^l,\mathbf f_e^l:12\times1.
$$

정확한 산술에서

$$
(\mathbf K_e^l)^T=\mathbf K_e^l,
\qquad
(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l\ge0,
\qquad
\operatorname{rank}(\mathbf K_e^l)=6.
$$

나머지 6개 null modes는 자유 요소의 강체운동이다. 구조 경계조건을 적용한 뒤에만
전역 자유도 강성행렬이 positive definite일 수 있다.

Mass와 damping은 현재 선형 정적 범위에서 N/A이다.

## 12. External Force Vector

### 12.1 일반 일관 분포하중

이 절과 12.2절은 요소 formulation을 검증하는 constant local line-load kernel의
수학 계약이다. V0 parser와 CLI는 Abaqus `*DLOAD`를 지원하거나 Domain distributed-load
object를 만들지 않으므로 이 kernel을 해석 pipeline의 지원 하중으로 해석하지 않는다.

국부 단위길이당 하중 벡터를

$$
\mathbf p_l(x)=
\begin{bmatrix}p_x&p_y&p_z&m_x&0&0\end{bmatrix}^T
$$

로 두면 일관 요소하중은

$$
\boxed{
\mathbf f_{\mathrm{dist},e}^l
=\int_0^L\mathbf N_q^T(x)\mathbf p_l(x)\,dx}
$$

이다. 성분형은

$$
f_a^l=\int_0^L\sum_{i=1}^{6}N_{q,ia}(x)p_{l,i}(x)\,dx.
$$

분포 $m_y,m_z$까지 지원하려면 $\mathbf p_l$의 마지막 두 성분과 이에 대응하는
strong form을 별도 I/O 계약에서 명시해야 한다. 현재 baseline에서는 제외한다.

### 12.2 상수 분포하중의 닫힌형

$p_x,p_y,p_z,m_x$가 요소에서 상수이면

$$
\boxed{
\mathbf f_{\mathrm{dist},e}^l=
\begin{bmatrix}
\frac{p_xL}{2}\\
\frac{p_yL}{2}\\
\frac{p_zL}{2}\\
\frac{m_xL}{2}\\
-\frac{p_zL^2}{12}\\
\frac{p_yL^2}{12}\\
\frac{p_xL}{2}\\
\frac{p_yL}{2}\\
\frac{p_zL}{2}\\
\frac{m_xL}{2}\\
\frac{p_zL^2}{12}\\
-\frac{p_yL^2}{12}
\end{bmatrix}}
$$

이다. 벡터의 5, 11번째 $M_y$ 성분 부호는 $w$ 보간에 $-H_2,-H_4$가
사용되기 때문에 발생한다.

### 12.3 요소 내부 집중하중과 끝하중

$x_p$에 작용하는 국부 집중 일반화 하중 $\mathbf P_p$의 등가절점하중은

$$
\mathbf f_p^l=\mathbf N_q^T(x_p)\mathbf P_p.
$$

절점 끝하중은 해당 요소 자유도 성분에 직접 더한다. 동일한 물리 하중을 구조 절점
하중과 요소 끝하중으로 중복 조립하지 않아야 한다.

## 13. Numerical Integration

### 13.1 강성 적분

자연좌표에서

$$
\mathbf K_e^l
=\int_{-1}^{1}\mathbf B^T(\xi)\mathbf D\mathbf B(\xi)J_x\,d\xi.
$$

상수 물성에 대해 $\mathbf B$의 굽힘행은 $\xi$의 1차식이므로
$\mathbf B^T\mathbf D\mathbf B$는 최대 2차식이다. 2점 Gauss–Legendre 적분은
3차 다항식까지 정확하므로 전체 강성을 정확히 적분한다.

$$
\xi_1=-\frac1{\sqrt3},\qquad
\xi_2=+\frac1{\sqrt3},\qquad
w_1=w_2=1.
$$

행렬 및 성분형은

$$
\boxed{
\mathbf K_e^l
=\sum_{g=1}^{2}\mathbf B^T(\xi_g)\mathbf D\mathbf B(\xi_g)J_xw_g}
$$

$$
K_{ab}^l
=\sum_{g=1}^{2}\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}(\xi_g)D_{ij}B_{jb}(\xi_g)J_xw_g.
$$

1점 적분은 굽힘 부분에 spurious zero-energy mode를 만들 수 있으므로 사용하지
않는다. 구현 초기에는 2점 Gauss 결과를 11.3절 닫힌형 행렬과 성분별 비교한다.

### 13.2 하중 적분

상수 분포하중에 대한 cubic Hermite 함수 적분은 2점 Gauss로 정확하다.

$$
\mathbf f_{\mathrm{dist},e}^l
=\sum_{g=1}^{2}\mathbf N_q^T(\xi_g)
\mathbf p_l(\xi_g)J_xw_g.
$$

요소 내에서 물성이나 하중이 변하는 기능을 추가할 때는 다항 차수와 discontinuity를
기준으로 적분차수를 다시 정해야 하며, 현재 2점을 일반 규칙으로 오해하면 안 된다.

## 14. Local-to-Global Transformation

### 14.1 3×3 회전행렬

전역 벡터를 국부 성분으로 변환하는 행렬을

$$
\boxed{
\mathbf R=
\begin{bmatrix}
\mathbf e_x^T\\
\mathbf e_y^T\\
\mathbf e_z^T
\end{bmatrix}}
$$

로 정의한다.

$$
\mathbf v^l=\mathbf R\mathbf v^g,
\qquad
v_i^l=\sum_{I=1}^{3}R_{iI}v_I^g.
$$

검사 조건은

$$
\mathbf R\mathbf R^T=\mathbf I_3,\qquad
\det\mathbf R=+1.
$$

### 14.2 12×12 요소 변환행렬

병진과 small-rotation vector는 같은 회전행렬로 변환한다.

$$
\boxed{
\mathbf T=\operatorname{diag}(\mathbf R,\mathbf R,\mathbf R,\mathbf R)}
$$

첫 두 block은 node 1의 병진/회전, 다음 두 block은 node 2의 병진/회전이다.

$$
\mathbf d_e^l=\mathbf T\mathbf d_e^g,
\qquad
d_a^l=\sum_{A=1}^{12}T_{aA}d_A^g.
$$

### 14.3 전역 요소식

가상일 불변성에서

$$
(\delta\mathbf d_e^g)^T\mathbf f_e^g
=(\delta\mathbf d_e^l)^T\mathbf f_e^l
$$

이므로

$$
\boxed{
\mathbf K_e^g=\mathbf T^T\mathbf K_e^l\mathbf T,
\qquad
\mathbf f_e^g=\mathbf T^T\mathbf f_e^l.}
$$

성분형은

$$
K_{AB}^g
=\sum_{a=1}^{12}\sum_{b=1}^{12}
T_{aA}K_{ab}^lT_{bB},
\qquad
f_A^g=\sum_{a=1}^{12}T_{aA}f_a^l.
$$

## 15. Global Assembly and Linear Static Solution

요소 전역 자유도를 전체 자유도에서 추출하는 Boolean/scatter 행렬을
$\mathbf A_e$라 하면

$$
\mathbf d_e^g=\mathbf A_e\mathbf d.
$$

조립식은

$$
\boxed{
\mathbf K=\sum_e\mathbf A_e^T\mathbf K_e^g\mathbf A_e,
\qquad
\mathbf F=\sum_e\mathbf A_e^T\mathbf f_e^g+\mathbf F_{\mathrm{nodal}}.}
$$

성분형은

$$
K_{IJ}=\sum_e\sum_{A=1}^{12}\sum_{B=1}^{12}
A_{e,AI}K_{e,AB}^gA_{e,BJ},
$$

$$
F_I=\sum_e\sum_{A=1}^{12}A_{e,AI}f_{e,A}^g+F_{\mathrm{nodal},I}.
$$

자유/구속 자유도로 partition하면

$$
\begin{bmatrix}
\mathbf K_{ff}&\mathbf K_{fc}\\
\mathbf K_{cf}&\mathbf K_{cc}
\end{bmatrix}
\begin{bmatrix}\mathbf d_f\\\mathbf d_c\end{bmatrix}
=
\begin{bmatrix}\mathbf F_f\\\mathbf F_c\end{bmatrix}.
$$

주어진 $\mathbf d_c$에 대해

$$
\boxed{
\mathbf K_{ff}\mathbf d_f
=\mathbf F_f-\mathbf K_{fc}\mathbf d_c}
$$

를 풀고, 구속 자유도 반력은 원래 조립식에서

$$
\boxed{
\mathbf R_c
=\mathbf K_{cf}\mathbf d_f+\mathbf K_{cc}\mathbf d_c-\mathbf F_c}
$$

로 복구한다. 성분형은

$$
R_I=\sum_JK_{IJ}d_J-F_I,\qquad I\in\mathcal C.
$$

## 16. Output Recovery

### 16.1 변위

전역 요소변위에서

$$
\mathbf d_e^l=\mathbf T\mathbf d_e^g,
\qquad
\mathbf q(\xi)=\mathbf N_q(\xi)\mathbf d_e^l
$$

로 국부 nodal/section 변위를 복구한다.

### 16.2 일반화 변형률과 단면력

두 Gauss point $\xi_g=\pm1/\sqrt3$에서 일반화 변형률과 generalized section
resultant를 다음과 같이 복구한다.

$$
\boxed{
\boldsymbol\varepsilon(\xi_g)=\mathbf B(\xi_g)\mathbf d_e^l,
\qquad
\mathbf s(\xi_g)=\mathbf D\mathbf B(\xi_g)\mathbf d_e^l.}
$$

성분형은

$$
\varepsilon_i(\xi_g)=\sum_aB_{ia}(\xi_g)d_a^l,
\qquad
s_i(\xi_g)=\sum_jD_{ij}\varepsilon_j(\xi_g).
$$

endpoint section resultant는 별도로 $\xi=-1,+1$에서 $\mathbf D\mathbf B\mathbf d_e^l$을
직접 평가한 local $[N,T,M_y,M_z]$이다. 이는 양의 local $x$ 면을 기준으로 한
section-cut 부호를 가지며 endpoint 외향 법선 부호를 곱하지 않는다. 상수 단면
요소에서 $N,T$는 요소 내 상수이고 $M_y,M_z$는 선형이다. 따라서 Gauss generalized
values와 endpoint section resultant는 위치와 row identity가 서로 다른 결과다.

### 16.3 축응력

단면점 $(y,z)$의 axial `S11`은

$$
\boxed{
S11(\xi,y,z)=\sigma_{xx}(\xi,y,z)
=E\left[\varepsilon_0(\xi)
+z\kappa_y(\xi)-y\kappa_z(\xi)\right].}
$$

`S11`은 두 Gauss point와 입력 `*SECTION POINTS`의 각 단면점에서 평가한다. 입력
section point가 없으면 centroid $(y,z)=(0,0)$을 `source=fesa-default`로 사용한다.
Abaqus section-point 좌표와의 승인 mapping은 $x1=y$, $x2=z$이다.

Euler–Bernoulli 가정은 횡전단변형을 제거하므로 이 정식화만으로 transverse shear
stress 분포를 constitutive recovery하지 않는다. 또한 일반 단면의 비틀림 전단응력
분포는 $J$ 하나만으로 결정되지 않으므로 baseline output은 $T$까지로 제한한다.
원형 단면처럼 별도 승인된 section recovery가 있을 때만
$\tau=T\rho/J$ 형태를 사용할 수 있으나 이는 V0 범위 밖이다. `S11`은 analytical/unit
test와 HDF5 contract로 검증하고 Abaqus beam-stress reference comparison은 명시적
N/A이다.

### 16.4 요소 끝력

분포하중을 포함하는 요소의 equilibrium end action은

$$
\boxed{
\mathbf f_{\mathrm{end},e}^l
=\mathbf K_e^l\mathbf d_e^l-\mathbf f_{\mathrm{dist},e}^l}
$$

이고 전역 성분은

$$
\mathbf f_{\mathrm{end},e}^g
=\mathbf T^T\mathbf f_{\mathrm{end},e}^l.
$$

각 endpoint row의 local component 순서는 $[F_X,F_Y,F_Z,M_X,M_Y,M_Z]$이다.
V0 parser/CLI는 distributed load를 생성하지 않으므로 승인된 해석 pipeline에서는
$\mathbf f_{\mathrm{dist},e}^l=\mathbf0$이다. 이 경우 $\xi_1=-1$, $\xi_2=+1$에서
외향 법선 부호를 각각 $n_{\mathrm{out},1}=-1$, $n_{\mathrm{out},2}=+1$로 두면
5.3절의 자연경계 부호와

$$
F_{X,i}=n_{\mathrm{out},i}N(\xi_i),\qquad
M_{X,i}=n_{\mathrm{out},i}T(\xi_i),\qquad
M_{Y,i}=n_{\mathrm{out},i}M_y(\xi_i),\qquad
M_{Z,i}=n_{\mathrm{out},i}M_z(\xi_i),
$$

$$
F_{Y,i}=-n_{\mathrm{out},i}
\left.\frac{dM_z}{dx}\right|_{\xi_i},\qquad
F_{Z,i}=n_{\mathrm{out},i}
\left.\frac{dM_y}{dx}\right|_{\xi_i}
$$

가 일치해야 한다. 따라서 이 V0 조건의 equilibrium end action은 요소 끝에 작용하는
outward action이고, endpoint section resultant는 동일 section을 양의 local $x$ 면으로
자른 값이다. Abaqus node-station force row에는 equilibrium end action이 아니라 section
resultant를 동일한 section-cut 부호로 정규화해 비교하며, interior node의 두 endpoint가
승인 tolerance 안에서 먼저 일치해야 한다. 승인된 component mapping은
`SF1 -> N`, `SM1 -> My`, `SM2 -> Mz`, `SM3 -> T`이고, CSV에 없는 transverse
`SF2`, `SF3`는 reference comparison 대상이 아니다.

$\mathbf D\mathbf B\mathbf d_e^l$로 얻는 section resultant와 위 equilibrium end
action은 분포하중이 있을 때 서로 다른 recovery 표현이다. cubic displacement만으로
uniform transverse load의 quadratic exact moment diagram을 한 요소에서 표현할 수
없으므로 formulation-only line-load kernel에서는 위 V0 endpoint 관계를 적용하거나
두 결과를 같은 row로 동일시하지 않는다.

### 16.5 반력

반력은 요소 끝력을 임의로 합하는 별도 경로가 아니라 15절의 전체 조립 잔차
$\mathbf K\mathbf d-\mathbf F$에서 구한다. 이렇게 해야 nodal load, element load 및
비영 prescribed displacement가 모두 반영된다.

### 16.6 Nodal extrapolation

- $N,T$: 상수이므로 extrapolation 불필요
- $M_y,M_z,\kappa_y,\kappa_z$: 선형이므로 두 Gauss point 값에서 끝점으로 선형
  extrapolation할 수 있으나, 직접 $\mathbf B(\xi=\pm1)$를 평가하는 방법을 기준으로 한다.
- 여러 요소의 nodal stress averaging은 재료/단면 불연속을 숨길 수 있으므로 현재
  baseline에서 수행하지 않는다.

## 17. Algorithm Pseudocode

```text
for each element e:
    read X1, X2, Abaqus first section axis n1, E, nu, A, I11, I12, I22, J
    require I12 == 0
    map a = n1, Iy = I11, Iz = I22, G = E / (2 * (1 + nu))

    DeltaX = X2 - X1
    L = norm(DeltaX)
    coordinate_scale = max(1, norm(X1), norm(X2))
    require L > 1e-12 * coordinate_scale

    ex = DeltaX / L
    ey_trial = a - dot(a, ex) * ex
    require norm(ey_trial) > 1e-12 * max(1, norm(a))
    ey = ey_trial / norm(ey_trial)
    ez = cross(ex, ey)

    R = rows(ex, ey, ez)
    require R * transpose(R) approximately equals identity
    require determinant(R) approximately equals +1
    T = block_diag(R, R, R, R)

    validate E, G, A, Iy, Iz, J are positive
    D = diag(E*A, G*J, E*Iy, E*Iz)

    initialize Kl[12,12] = 0
    initialize fl[12] = 0

    for (xi, weight) in {(-1/sqrt(3), 1), (+1/sqrt(3), 1)}:
        Jx = L / 2
        r = (1 + xi) / 2

        evaluate N1L, N2L
        evaluate H1...H4 and first/second x-derivatives
        form Nq[6,12]
        form B[4,12]

        Kl += transpose(B) * D * B * Jx * weight
        if verifying the formulation-only constant local line-load kernel:
            fl += transpose(Nq) * local_line_load(xi) * Jx * weight

    in verification builds/checks:
        compare Kl component-wise with the closed-form matrix in Section 11.3
        check symmetry, rank 6, and rigid-body null modes

    Kg = transpose(T) * Kl * T
    fg = transpose(T) * fl
    assemble Kg and fg with the element-to-global DOF map

partition global free/constrained DOFs
solve Kff * df = Ff - Kfc * dc
recover Rc = Kcf * df + Kcc * dc - Fc

for each element output recovery:
    de_global = gather global nodal DOFs
    de_local = T * de_global

    for each mandatory Gauss point xi_g:
        epsilon = B(xi_g) * de_local
        section_resultant = D * epsilon
        recover S11 at each input section point or the fesa-default centroid

    for each endpoint xi in {-1, +1}:
        recover endpoint section_resultant = D * B(xi) * de_local
    recover local equilibrium end actions with the endpoint outward-action convention
recover reactions from the assembled residual, not an end-action sum
```

이 의사코드는 수학적 계산 순서만 정의하며 C++ class, signature 또는 storage API를
정의하지 않는다.

## 18. Verification and Benchmarks

[[Finite Element Program Implementation]]은 element-local 계산, 전역 조립, 해,
결과 복구를 분리하고 새 요소를 patch/benchmark test로 검증할 것을 요구한다.
[[Finite Element Modeling and Convergence Checks]]는 compatibility, equilibrium,
convergence 및 결과 해석을 함께 검사해야 한다는 근거를 제공한다.

### 18.1 형상함수 단위검사

1. 8.3절의 값/기울기 Kronecker 조건을 machine precision 범위에서 검사한다.
2. $N_1^L+N_2^L=1$, $H_1+H_3=1$을 검사한다.
3. 상수 $u,\theta_x$를 정확히 재현한다.
4. 선형 $v=c_0+c_1x$에 대해 $\theta_z=c_1$, $\kappa_z=0$인지 검사한다.
5. 선형 $w=c_0+c_1x$에 대해 $\theta_y=-c_1$, $\kappa_y=0$인지 검사한다.

### 18.2 행렬 불변량

승인된 normalized matrix 기준은

$$
\operatorname{normalized\_error}
(\mathbf K_e^l,(\mathbf K_e^l)^T)\le10^{-12},
$$

$$
\operatorname{normalized\_error}
(\mathbf K_{\mathrm{Gauss}}^l,\mathbf K_{\mathrm{closed}}^l)\le10^{-12}
$$

를 검사한다. 자유 요소의 고유치 또는 rank 검사에서는 6개 rigid-body mode와 6개
positive deformation mode가 있어야 한다.

대표 국부 강체모드는

$$
\mathbf r_{T_x}=
\begin{bmatrix}1&0&0&0&0&0&1&0&0&0&0&0\end{bmatrix}^T,
$$

$$
\mathbf r_{T_y}=
\begin{bmatrix}0&1&0&0&0&0&0&1&0&0&0&0\end{bmatrix}^T,
$$

$$
\mathbf r_{T_z}=
\begin{bmatrix}0&0&1&0&0&0&0&0&1&0&0&0\end{bmatrix}^T,
$$

$$
\mathbf r_{R_x}=
\begin{bmatrix}0&0&0&1&0&0&0&0&0&1&0&0\end{bmatrix}^T,
$$

$$
\mathbf r_{R_y}=
\begin{bmatrix}0&0&0&0&1&0&0&0&-L&0&1&0\end{bmatrix}^T,
$$

$$
\mathbf r_{R_z}=
\begin{bmatrix}0&0&0&0&0&1&0&L&0&0&0&1\end{bmatrix}^T.
$$

각 $k$에 대해

$$
\operatorname{normalized\_residual}
(\mathbf K_e^l\mathbf r_k)\le10^{-10}
$$

이어야 한다. Analytical axial, torsion 및 두 굽힘 평면 응답은 relative $10^{-9}$를
적용한다. Norm의 무차원화와 rank scale은 Numerical Review에서 이 승인 기준의 의미를
바꾸지 않는 계산 규약으로 명시한다.

### 18.3 기본 해석해 benchmark

길이 $L$의 cantilever 하나를 사용하여 다음을 검사한다.

| 하중 | 비교량 | 해석해 |
| --- | --- | --- |
| 끝 축력 $P_x$ | $u(L)$ | $P_xL/(EA)$ |
| 끝 torque $M_x$ | $\theta_x(L)$ | $M_xL/(GJ)$ |
| 끝 횡력 $P_y$ | $v(L)$ | $P_yL^3/(3EI_z)$ |
| 끝 횡력 $P_y$ | $\theta_z(L)$ | $P_yL^2/(2EI_z)$ |
| 끝 횡력 $P_z$ | $w(L)$ | $P_zL^3/(3EI_y)$ |
| 끝 횡력 $P_z$ | $\theta_y(L)$ | $-P_zL^2/(2EI_y)$ |
| 상수 분포하중 $p_y$ | $v(L)$ | $p_yL^4/(8EI_z)$ |
| 상수 분포하중 $p_y$ | $\theta_z(L)$ | $p_yL^3/(6EI_z)$ |

한 개의 cubic Hermite 요소는 끝 집중력/모멘트에 대한 cubic 이하 변위장을 정확히
재현한다. 상수 분포 횡하중의 exact displacement는 4차이므로 mesh refinement에
따른 수렴도 함께 확인한다.

### 18.4 좌표변환 검사

동일한 beam 문제를 전역 $x,y,z$축과 임의의 3차원 방향으로 rigid rotation하여
해석한다. 변위와 힘을 원래 좌표계로 되돌렸을 때 일치해야 한다.

$$
\mathbf T^T\mathbf T=\mathbf I_{12},
\qquad
(\mathbf d_e^g)^T\mathbf K_e^g\mathbf d_e^g
=(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l.
$$

### 18.5 평형 검사

$$
\left\|\sum\mathbf F_{\mathrm{applied}}+
\sum\mathbf R\right\|\le\epsilon_F,
$$

$$
\left\|\sum\mathbf M_{\mathrm{applied}}+
\sum\left(\mathbf X\times\mathbf F\right)+
\sum\mathbf M_{\mathrm{reaction}}\right\|\le\epsilon_M.
$$

승인된 B33 reference row는 같은 model, step/frame, quantity, component별로 Abaqus
reference 값만 사용해

$$
\operatorname{reference\_scale}=\max_i|r_i|,\qquad
\operatorname{row\_tolerance}=\operatorname{absolute\_floor}
+10^{-6}\operatorname{reference\_scale}
$$

를 적용한다. 승인된 SI bundle의 displacement/rotation floor는 $10^{-9}$,
force/moment floor는 $10^{-3}$이다. Reference 값을 zero-clamp하거나 row를 제거하지
않으며 missing, extra, duplicate, nonfinite, schema/identity mismatch는 수치 판정 전에
실패한다. 이 reference policy는 위 formulation/analytical tolerance를 대체하지 않는다.

## 19. Numerical Risks

| 위험 | 영향 | 필수 대응/검사 |
| --- | --- | --- |
| $L\approx0$ | Jacobian singularity, $L^{-3}$ 발산 | 조립 전 길이 검증 |
| guide vector와 $\mathbf e_x$ 평행 | 국부 $y,z$축 미정 | projection norm 검사 후 입력 오류 |
| $E,G,A,I_y,I_z,J\le0$ | 비물리/indefinite stiffness | 물성·단면 검증 |
| 짧거나 깊은 beam | 전단변형 누락으로 과도하게 stiff | Euler–Bernoulli 적용 한계로 기록; 승인되지 않은 자동 warning/error threshold를 두지 않고 B31은 지원 요소로 재해석하지 않음 |
| open section warping | $GJ/L$만으로 비틀림 부족 | 현재 범위 밖으로 명시 |
| 하중이 shear center를 벗어남 | 굽힘–비틀림 결합 누락 | eccentric load를 명시적 force+torque로 환산하거나 후속 기능화 |
| $I_y/I_z$ 또는 stiffness scale 극단값 | 전역 condition number 악화 | 단위/스케일 및 solver conditioning 검사 |
| 1점 굽힘 적분 | spurious mechanism | 2점 Gauss 또는 닫힌형 사용 |
| $\theta_y=-w'$ 부호 누락 | $w$–$M_y$ 결과 반전 | B, K, cantilever sign test |
| local/global 변환 방향 혼동 | 회전된 요소 결과 오류 | $\mathbf d^l=\mathbf T\mathbf d^g$ 계약 및 energy invariance test |
| free structure 또는 불충분 구속 | singular global matrix | 6 rigid modes와 구속조건 진단 |
| nodal averaging | 단면 불연속 stress 오염 | baseline에서 averaging 금지 |

Euler–Bernoulli 요소에는 shear strain 자유도가 없으므로 shear locking은 N/A이다.
Volumetric locking과 hourglass control도 N/A이지만, 1점 굽힘 적분에 의한 rank
deficiency는 별도의 under-integration 위험이다.

## 20. Evidence and Wiki Gap

### 20.1 사용한 FEM wiki 근거

- [[Beam and Frame Finite Elements]]: Euler–Bernoulli cubic displacement, nodal
  displacement/rotation, axial·bending·torsion 결합 및 공간좌표 변환
- [[Isoparametric Finite Elements]]: 자연좌표, 형상함수, Jacobian, derivative mapping,
  quadrature 및 distortion/under-integration 위험
- [[Finite Element Method]]: strong/variational statement에서 요소방정식, 조립, 해로
  이어지는 displacement FEM 흐름
- [[Finite Element Program Implementation]]: element-local 계산, 전역 조립, 결과
  복구 및 patch/benchmark verification
- [[Finite Element Modeling and Convergence Checks]]: compatibility, equilibrium,
  mesh convergence 및 결과 해석
- [[A-First-Course-in-the-Finite-Element-Method|A First Course in the Finite Element Method]]:
  beam, frame, grid, space-frame 및 isoparametric formulation의 교재 근거

### 20.2 근거의 적용 경계

Tier와 provenance가 정리된 research brief는
`docs/linear-static-3d-euler-beam/research.md`에 있다. Wiki와 원출처는 핵심
beam 이론과 Abaqus component 의미를 제공하고, 이 문서의 DOF 순서,
$\theta_y=-w'$ 부호, guide-vector 축, 12×12 행렬, line-load vector, 결과 위치와
tolerance는 승인된 project requirement/design과 결합한 FESA 계약이다. 따라서 exact
Abaqus B33 내부 알고리즘 동등성이나 full compatibility를 주장하지 않는다.

Research brief의 PARDISO manual 및 deterministic COO-to-CSR 근거 gap은 solver adapter와
assembly 구현 계획의 handoff이며 이 요소의 수학 정식화를 변경하지 않는다. Endpoint
sign-normalization gap은 5.3절과 16.4절의 outward end action/positive-face section-cut
관계로 명시했으며 Numerical Review가 독립 재유도해야 한다.

## 21. Resolved Decisions, Remaining Limitations, and Downstream Handoff

### 21.1 Resolved decisions and remaining limitations

1. Requirement와 research source는 Metadata의 실제 Step 0–1 산출물로 해결되었다.
2. Orientation과 section mapping은 `n1 -> local y`, `t x n1 -> local z`,
   $I_y=I11$, $I_z=I22$, $I12=0$으로 해결되었다.
3. Geometry validation은 3.1절의 scale-aware length/projection 기준으로 해결되었다.
4. Matrix, residual, analytical 및 B33 reference tolerance는 18절의 승인 기준으로
   해결되었다.
5. equilibrium end action, endpoint section resultant, Gauss generalized values 및
   section-point `S11`은 16절에서 위치·component·부호가 구분되었다. Stress reference
   comparison은 N/A이다.
6. Euler–Bernoulli slenderness는 적용 limitation으로 남는다. 공식 자료의 약 $1/15$
   지침을 FESA warning/error threshold로 바꾸지 않는다.
7. Torsional/transverse shear stress, warping, $I12\ne0$, B31/Timoshenko 및 `*DLOAD`
   parser support는 승인 범위 밖이다.

현재 Numerical Review를 막는 미해결 수학 입력은 없다. `ready-for-numerical-review`는
독립 수치 검토 준비 상태이며 정식화 최종 승인이나 release readiness를 뜻하지 않는다.

### 21.2 Numerical Review Agent

- $\theta_y=-w'$와 $M_y=-EI_yw''$ 부호를 독립적으로 재유도한다.
- 4×12 $\mathbf B$ 적분이 12×12 닫힌형과 일치하는지 symbolic/numerical 검사한다.
- 6 rigid-body modes, rank 6, symmetry 및 energy invariance를 검사한다.
- 자연경계력 $F_y,F_z,M_y,M_z$와 요소 끝력 부호를 검사한다.
- 16.4절의 outward equilibrium end action과 positive-face endpoint section resultant
  관계를 양 endpoint에서 독립 재유도한다.
- 2점 Gauss exactness와 1점 적분 rank deficiency를 확인한다.

### 21.3 I/O Definition Agent

- node coordinates, `n1`, $E,\nu,A,I11,I12,I22,J$ 입력과 승인 mapping을 정의한다.
- guide vector 평행, zero length, nonpositive property validation을 정의한다.
- local/global displacement, reaction, equilibrium end action, endpoint section resultant,
  Gauss generalized strain/resultant, section-point axial `S11`의 HDF5 schema와 위치를
  구분한다.
- Abaqus keyword subset은 별도 승인 범위만 정의하며 full compatibility를 주장하지 않는다.

### 21.4 Reference Model Agent

- axial cantilever, torsion cantilever, $y/z$ 양방향 bending cantilever를 분리한다.
- local 축과 일치하는 모델 및 임의 3차원 회전 모델을 모두 준비한다.
- formulation-only uniform line-load kernel은 unit/analytical coverage로 두고 `*DLOAD`
  reference support로 확대하지 않는다.
- 승인된 `reference/cantilever beam/`의 exact legacy path를 read-only로 유지하고
  reference artifact를 생성·수정·복원하지 않는다.

### 21.5 Implementation Planning Agent

TDD 순서는 다음이 적합하다.

1. Hermite 값/미분 및 mapping 단위검사
2. $\mathbf B$ 성분과 닫힌형 $\mathbf K_e^l$ 비교검사
3. rigid modes, symmetry, rank 및 positive-energy 검사
4. 좌표변환 orthogonality와 rotated-element invariance 검사
5. consistent uniform line-load vector 검사
6. axial/torsion/two-plane cantilever 해석해 검사
7. 조립, prescribed displacement 및 reaction 통합검사
8. 결과 복구와 end-action sign 검사

Production C++ 파일을 변경할 때는 관련 C++ test를 먼저 작성해 RED 실패를 확인한
뒤 GREEN 구현과 MSVC x64 Debug/CTest VERIFY를 수행해야 한다.