実装

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

/* 体積 */
void WH_Fem__Simplex3D__V
(double whV_x0[3], 
 double whV_x1[3], 
 double whV_x2[3],
 double whV_x3[3],
 double *OUT__V)
/*
  入力引数:
    whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
  出力引数:
    OUT__V はスカラーで、体積
*/
{
  double x0 = whV_x0[0];
  double x1 = whV_x1[0];
  double x2 = whV_x2[0];
  double x3 = whV_x3[0];
  double y0 = whV_x0[1];
  double y1 = whV_x1[1];
  double y2 = whV_x2[1];
  double y3 = whV_x3[1];
  double z0 = whV_x0[2];
  double z1 = whV_x1[2];
  double z2 = whV_x2[2];
  double z3 = whV_x3[2];

  *OUT__V = 1.0 / 6 * 
    (- x0 * y1 * z2 + x0 * y1 * z3 + x0 * y2 * z1
     - x0 * y2 * z3 - x0 * y3 * z1 + x0 * y3 * z2
     + x1 * y0 * z2 - x1 * y0 * z3 - x1 * y2 * z0
     + x1 * y2 * z3 + x1 * y3 * z0 - x1 * y3 * z2
     - x2 * y0 * z1 + x2 * y0 * z3 + x2 * y1 * z0
     - x2 * y1 * z3 - x2 * y3 * z0 + x2 * y3 * z1
     + x3 * y0 * z1 - x3 * y0 * z2 - x3 * y1 * z0
     + x3 * y1 * z2 + x3 * y2 * z0 - x3 * y2 * z1);
}