実装

形状関数の勾配を体積座標で評価する関数は それぞれ以下のようになる。 これは、全節点の形状関数について一度に計算される。

/* 形状関数の勾配(体積座標で評価) */
void WH_Fem__Simplex3D10N__whP_N_whV_x_arrayIn
(double whV_x0[3], 
 double whV_x1[3], 
 double whV_x2[3],
 double whV_x3[3],
 double L0, double L1, double L2, double L3,
 double OUT__whP_N_whV_x_arrayIn[10][3])
/*
  入力引数:
    whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
    L0, L1, L2, L3 はスカラーで、体積座標
  出力引数:
    OUT__whP_N_whV_x_arrayIn はベクトルの節点配列で、
      形状関数の勾配
*/
{
  double V;
  double b0, b1, b2, b3;
  double c0, c1, c2, c3;
  double d0, d1, d2, d3;
  
  WH_Fem__Simplex3D__V
    (whV_x0, whV_x1, whV_x2, whV_x3, 
     &V);
  WH_Fem__Simplex3D__bcd
    (whV_x0, whV_x1, whV_x2, whV_x3, V,
     &b0, &b1, &b2, &b3, 
     &c0, &c1, &c2, &c3,
     &d0, &d1, &d2, &d3);

  OUT__whP_N_whV_x_arrayIn[0][0] 
    = 4.0 * b0 * L0 - b0;
  OUT__whP_N_whV_x_arrayIn[1][0] 
    = 4.0 * b1 * L1 - b1;
  OUT__whP_N_whV_x_arrayIn[2][0] 
    = 4.0 * b2 * L2 - b2;
  OUT__whP_N_whV_x_arrayIn[3][0] 
    = 4.0 * b3 * L3 - b3;
  OUT__whP_N_whV_x_arrayIn[4][0] 
    = 4.0 * b1 * L0 + 4.0 * b0 * L1;
  OUT__whP_N_whV_x_arrayIn[5][0] 
    = 4.0 * b2 * L0 + 4.0 * b0 * L2;
  OUT__whP_N_whV_x_arrayIn[6][0] 
    = 4.0 * b3 * L0 + 4.0 * b0 * L3;
  OUT__whP_N_whV_x_arrayIn[7][0] 
    = 4.0 * b2 * L1 + 4.0 * b1 * L2;
  OUT__whP_N_whV_x_arrayIn[8][0] 
    = 4.0 * b3 * L2 + 4.0 * b2 * L3;
  OUT__whP_N_whV_x_arrayIn[9][0] 
    = 4.0 * b1 * L3 + 4.0 * b3 * L1;

  OUT__whP_N_whV_x_arrayIn[0][1] 
    = 4.0 * c0 * L0 - c0;
  OUT__whP_N_whV_x_arrayIn[1][1] 
    = 4.0 * c1 * L1 - c1;
  OUT__whP_N_whV_x_arrayIn[2][1] 
    = 4.0 * c2 * L2 - c2;
  OUT__whP_N_whV_x_arrayIn[3][1] 
    = 4.0 * c3 * L3 - c3;
  OUT__whP_N_whV_x_arrayIn[4][1] 
    = 4.0 * c1 * L0 + 4.0 * c0 * L1;
  OUT__whP_N_whV_x_arrayIn[5][1] 
    = 4.0 * c2 * L0 + 4.0 * c0 * L2;
  OUT__whP_N_whV_x_arrayIn[6][1] 
    = 4.0 * c3 * L0 + 4.0 * c0 * L3;
  OUT__whP_N_whV_x_arrayIn[7][1] 
    = 4.0 * c2 * L1 + 4.0 * c1 * L2;
  OUT__whP_N_whV_x_arrayIn[8][1] 
    = 4.0 * c3 * L2 + 4.0 * c2 * L3;
  OUT__whP_N_whV_x_arrayIn[9][1] 
    = 4.0 * c1 * L3 + 4.0 * c3 * L1;

  OUT__whP_N_whV_x_arrayIn[0][2] 
    = 4.0 * d0 * L0 - d0;
  OUT__whP_N_whV_x_arrayIn[1][2] 
    = 4.0 * d1 * L1 - d1;
  OUT__whP_N_whV_x_arrayIn[2][2] 
    = 4.0 * d2 * L2 - d2;
  OUT__whP_N_whV_x_arrayIn[3][2] 
    = 4.0 * d3 * L3 - d3;
  OUT__whP_N_whV_x_arrayIn[4][2] 
    = 4.0 * d1 * L0 + 4.0 * d0 * L1;
  OUT__whP_N_whV_x_arrayIn[5][2] 
    = 4.0 * d2 * L0 + 4.0 * d0 * L2;
  OUT__whP_N_whV_x_arrayIn[6][2] 
    = 4.0 * d3 * L0 + 4.0 * d0 * L3;
  OUT__whP_N_whV_x_arrayIn[7][2] 
    = 4.0 * d2 * L1 + 4.0 * d1 * L2;
  OUT__whP_N_whV_x_arrayIn[8][2] 
    = 4.0 * d3 * L2 + 4.0 * d2 * L3;
  OUT__whP_N_whV_x_arrayIn[9][2] 
    = 4.0 * d1 * L3 + 4.0 * d3 * L1;
}