実装

形状関数の微分を自然座標で評価する関数は 以下のようになる。 これは、全節点の形状関数について一度に計算される。

/* 形状関数の自然座標についての微分 */
void WH_Fem__Shape3D__Hexa8N__whP_N_whV_xi_arrayIn
(double whV_xi[3],
 double OUT__whP_N_whV_xi_arrayIn[8][3])
/*
  入力引数:
    whV_xi はベクトルで、自然座標
  出力引数:
    OUT__whP_N_whV_xi_arrayIn はベクトルの節点配列で、
      形状関数の自然座標に関する微分
*/
{
  double xi = whV_xi[0];
  double eta = whV_xi[1];
  double zeta = whV_xi[2];

  OUT__whP_N_whV_xi_arrayIn[0][0] 
    = -0.125 * (1.0 - eta) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[1][0] 
    =  0.125 * (1.0 - eta) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[2][0] 
    =  0.125 * (1.0 + eta) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[3][0] 
    = -0.125 * (1.0 + eta) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[4][0] 
    = -0.125 * (1.0 - eta) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[5][0] 
    =  0.125 * (1.0 - eta) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[6][0] 
    =  0.125 * (1.0 + eta) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[7][0] 
    = -0.125 * (1.0 + eta) * (1.0 + zeta);

  OUT__whP_N_whV_xi_arrayIn[0][1] 
    = -0.125 * (1.0 - xi) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[1][1] 
    = -0.125 * (1.0 + xi) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[2][1] 
    =  0.125 * (1.0 + xi) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[3][1] 
    =  0.125 * (1.0 - xi) * (1.0 - zeta);
  OUT__whP_N_whV_xi_arrayIn[4][1] 
    = -0.125 * (1.0 - xi) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[5][1] 
    = -0.125 * (1.0 + xi) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[6][1] 
    =  0.125 * (1.0 + xi) * (1.0 + zeta);
  OUT__whP_N_whV_xi_arrayIn[7][1] 
    =  0.125 * (1.0 - xi) * (1.0 + zeta);

  OUT__whP_N_whV_xi_arrayIn[0][2] 
    = -0.125 * (1.0 - xi) * (1.0 - eta);
  OUT__whP_N_whV_xi_arrayIn[1][2] 
    = -0.125 * (1.0 + xi) * (1.0 - eta);
  OUT__whP_N_whV_xi_arrayIn[2][2] 
    = -0.125 * (1.0 + xi) * (1.0 + eta);
  OUT__whP_N_whV_xi_arrayIn[3][2] 
    = -0.125 * (1.0 - xi) * (1.0 + eta);
  OUT__whP_N_whV_xi_arrayIn[4][2] 
    =  0.125 * (1.0 - xi) * (1.0 - eta);
  OUT__whP_N_whV_xi_arrayIn[5][2] 
    =  0.125 * (1.0 + xi) * (1.0 - eta);
  OUT__whP_N_whV_xi_arrayIn[6][2] 
    =  0.125 * (1.0 + xi) * (1.0 + eta);
  OUT__whP_N_whV_xi_arrayIn[7][2] 
    =  0.125 * (1.0 - xi) * (1.0 + eta);
}