実装

表面の法線方向ベクトルを求める 関数は以下のようになる。

/* 表面の法線方向ベクトル */
void WH_Fem__Simplex3D__whV_n
(double whV_x0[3], 
 double whV_x1[3], 
 double whV_x2[3],
 double whV_x3[3],
 int Iface,
 double OUT__whV_n[3])
/*
  入力引数:
    whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
    Iface は0から3までの整数で、表面ID
  出力引数:
    OUT__whV_n はベクトルで、法線方向
*/
{
  double whV_tmp0[3];
  double whV_tmp1[3];
  double whV_tmp2[3];

  switch (Iface) {
  case 0:
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x1, whV_x2,
       whV_tmp0);
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x2, whV_x3,
       whV_tmp1);    
    break;
  case 1:
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x0, whV_x3,
       whV_tmp0);
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x3, whV_x2,
       whV_tmp1);    
    break;
  case 2:
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x0, whV_x1,
       whV_tmp0);
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x1, whV_x3,
       whV_tmp1);    
    break;
  case 3:
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x0, whV_x2,
       whV_tmp0);
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_x2, whV_x1,
       whV_tmp1);    
    break;
  default:
    assert(0);
    break;
  }
  WH_Tensor3D__whV_times_whV_OUT_whV
    (whV_tmp0, whV_tmp1,
     whV_tmp2);
  WH_Tensor3D__whNormalize_whV_OUT_whV
    (whV_tmp2,
     OUT__whV_n);
}