実装

表面の面積を求める 関数は以下のようになる。

/* 表面の面積 */
void WH_Fem__Simplex3D__S
(double whV_x0[3], 
 double whV_x1[3], 
 double whV_x2[3],
 double whV_x3[3],
 int Iface,
 double *OUT__S)
/*
  入力引数:
    whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
    Iface は0から3までの整数で、稜線ID
  出力引数:
    OUT__S はスカラーで、表面の面積
*/
{
  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__whLen_whV_OUT_s
    (whV_tmp2,
     OUT__S);
  *OUT__S /= 2;
}