実装

与えられた位置 $\{ x^{src} \}$ から 対応する自然座標 $\{ \xi^{dest} \}$ を求める 関数は以下のようになる。

/* 位置から自然座標への変換 */
void WH_Fem__Isoparam3D__Volume__whV_xi
(WH_Fem__Element3DType elementType,
 int nNodes,
 double whV_x_arrayIn[/* nNodes */][3],
 double whV_x_src[3],
 double whV_xi_ini[3],
 double epsilon,
 double OUT__whV_xi_dest[3])
/*
  入力引数:
    elementType は要素のタイプ
    nNodes は整数で、要素節点数
    whV_x_arrayIn はベクトルの節点配列で、位置（節点座標）
    whV_x_src はベクトルで、位置
    whV_xi_ini はベクトルで、自然座標の初期値
    epsilon はスカラーで、トレランス
  出力引数:
    OUT__whV_xi_dest はベクトルで、自然座標
*/
{
  WH_Tensor3D__assign_whV_OUT_whV
    (whV_xi_ini,
     OUT__whV_xi_dest);

  for (;;) {
    double N_arrayIn[WH_FEM__MAX_ELEMENT_NODES];
    double whV_f[3];
    double whP_N_whV_xi_arrayIn[WH_FEM__MAX_ELEMENT_NODES][3];
    double whT_J[3][3];
    double whInc_whV_xi[3];
    double whV_tmp0[3];
    double whT_tmp1[3][3];
    double whT_tmp2[3][3];
    double whV_tmp3[3];
    double tmp4;

    WH_Fem__Shape3D__N_arrayIn
      (elementType, nNodes, OUT__whV_xi_dest, 
       N_arrayIn);
    WH_Fem__Isoparam3D__Volume__whV_x
      (nNodes, N_arrayIn, whV_x_arrayIn,
       whV_tmp0);
    WH_Tensor3D__whV_sub_whV_OUT_whV
      (whV_tmp0, whV_x_src,
       whV_f);

    WH_Fem__Shape3D__whP_N_whV_xi_arrayIn
      (elementType, nNodes, OUT__whV_xi_dest,
       whP_N_whV_xi_arrayIn);
    WH_Fem__Isoparam3D__Volume__whPt_whV_x_whV_xi
      (nNodes, whP_N_whV_xi_arrayIn, whV_x_arrayIn,
       whT_J);

    WH_Tensor3D__whTrans_whT_OUT_whT
      (whT_J,
       whT_tmp1);
    WH_Tensor3D__whInv_whT_OUT_whT
      (whT_tmp1,
       whT_tmp2);
    WH_Tensor3D__whT_cdot_whV_OUT_whV
      (whT_tmp2, whV_f,
       whV_tmp3);
    WH_Tensor3D__s_mul_whV_OUT_whV
      (-1.0, whV_tmp3,
       whInc_whV_xi);    

    WH_Tensor3D__whV_add_whV_OUT_whV
      (OUT__whV_xi_dest, whInc_whV_xi,
       OUT__whV_xi_dest);

    WH_Tensor3D__whLen_whV_OUT_s
      (whInc_whV_xi,
       &tmp4);
    if (tmp4 < epsilon) break;
  }
}