flag=1 で c_hat の結合（オフ対角）成分の副対称性が欠落
Bug #2 (closed)-2 — flag=1 で c_hat の結合（オフ対角）成分の副対称性が欠落
対象バージョン: FrontISTR 5.9 — commit 6d3a1684 (2026-05-19) 修正コミット: ce5ba6de（ブランチ debug/bug2-c-hat-symmetry、Bug #2 (closed)-1 fdf8efbe に続く）
1. 位置づけ（Bug #2 (closed)-1 との関係）
Bug #2 (closed)-1（fdf8efbe）は flag=1 ブランチで構成テンソル c_hat の副対称性が欠落する問題を、 変換ループ直後に副対称 9 成分を明示代入することで補完した。ただしその補完は 純せん断の対角成分（G12 / G23 / G31、すなわち c_hat(1,2,1,2) c_hat(2,3,2,3) c_hat(3,1,3,1) とその入れ替え）に限られていた。
D_hat には対角のせん断成分のほかに結合（オフ対角）成分がある。直交異方性でも 繊維角 angle ≠ 0 の場合、構成行列の回転によって法線–せん断の連成項（せん断–伸び結合）が 生じ、これらオフ対角成分が非ゼロになる。Bug #2 (closed)-1 の対角限定の補完ではこの結合成分の 副対称性が補完されず、依然として対称性が崩れる。
2. 原因
fistr1/src/lib/physics/ElasticLinear.f90、flag=1 ブランチの D_hat(5×5) → c_hat 変換ループ。 Bug #2 (closed)-1 はループ外で対角せん断成分のみ補完していた。結合成分を含む全成分に対して 副対称性を保証するには、ループ内で各 (is, js) 成分を代入する際に 4 つの副対称 置換すべてに同じ値を設定する必要がある。
c(i,j,k,l) = c(j,i,k,l) = c(i,j,l,k) = c(j,i,l,k)繊維角 0°（軸整列直交異方性）では結合成分が 0 のため Bug #2 (closed)-1 の補完で足りるが、 θ≠0 では結合成分が非ゼロになり不足する。
3. 修正
ElasticLinear.f90 の変換ループ本体で、各 D_hat 成分を 4 つの副対称インデックスへ 同時に代入する（Bug #2 (closed)-1 のループ外 9 行補完は本修正に置き換え）。
c_hat(ii, ij, ik, il) = D_hat(is, js)
+ c_hat(ij, ii, ik, il) = D_hat(is, js) ! i↔j
+ c_hat(ii, ij, il, ik) = D_hat(is, js) ! k↔l
+ c_hat(ij, ii, il, ik) = D_hat(is, js) ! both4. 検証の考え方（0 固有値の出現）
副対称性を欠いた c_hat は θ≠0 で剛性 K を欠陥化し、拘束された片持ち梁に本来現れない 零エネルギーモード（λ=0, 0 Hz）が生じる。正しい c_hat なら現れない。これを修正前後で比較する。
angle はラジアン入力：本バイナリ（
fdf8efbe/ce5ba6de）は Bug#3（度→ラジアン変換）未修正の ため angle はラジアン解釈。+45°=+0.785398 rad、−45°=−0.785398 rad を入力。
5. テストモデル・ケース構成
強異方性素材を ±45° で配置した片持ち梁（MITC3 デローニ不規則格子）。x=0 端全 DOF 拘束、 5 モードの固有値解析（!SOLUTION, TYPE=EIGEN）。
物性（flag=1, !ITEM=1, SUBITEM=9）: E1=140000, E2=10000, ν12=0.3, G12=G23=G31=5000 MPa, ρ=1.6E-9 t/mm³。要素 731 / 165 節点・268 要素 / 板厚 10mm。
| シナリオ | angle 入力 | 物理繊維角 |
|---|---|---|
case_p45（01_before / 03_after） |
+0.785398 rad | +45° |
case_m45（02_before / 04_after） |
−0.785398 rad | −45° |
フラット構成（before=01,02 / after=03,04）。before=fdf8efbe（②適用前）、after=ce5ba6de（②適用後）。 各ケースは入力一式（Job_1.msh / Job_1.cnt / hecmw_ctrl.dat）と run.log / FSTR.msg / 結果 Job_1.res.0.* を含む。
6. 検証結果
固有値 λ と振動数（5 モード）。before=fdf8efbe、after=ce5ba6de。TOTAL MASS=8.0000E-04。
−45°（case_m45）— 0 固有値の出現
| 実行 | mode1 | mode2 | mode3 | mode4 | mode5 |
|---|---|---|---|---|---|
02_case_m45_before |
12.216 | 0.000（λ=0, 残差15.7＝非収束の零モード） | 362.30 | 364.74 | 528.72 |
04_case_m45_after |
29.561 | 172.82 | 258.49 | 274.38 | 520.00 |
修正前は λ=0 の零エネルギーモードが出現（残差も収束しない）。修正後は全モード正値で零モード消失
+45°（case_p45）— 値は変わるが零モードは出ない
| 実行 | mode1 | mode2 | mode3 | mode4 | mode5 |
|---|---|---|---|---|---|
01_case_p45_before |
30.244 | 204.66 | 245.29 | 289.02 | 516.50 |
03_case_p45_after |
24.227 | 175.63 | 206.73 | 270.59 | 455.87 |
欠落していた結合項は sin2θ に odd（符号が方位で反転）。−45° 対角では零モードを生み、 +45° 対角では正定値を保ったまま誤値、と方位で症状が異なる。