実装

熱伝導行列の 節点 $In$ 節点 $Jn$ における成分を計算する 関数は、２次元と３次元問題についてそれぞれ以下のようになる。

/* 熱伝導行列：２次元 */
void WH_ThermalFemLinear2D__K_cnd_InJn
(double whV_N_X_In[2], double whV_N_X_Jn[2],
 double whT_lambda[2][2],
 double *OUT__K_cnd_InJn)
/*
  入力引数:
    whV_N_X_In, whV_N_X_Jn はベクトルで、
      それぞれ節点 In, Jn の補間関数勾配ベクトル
    whT_lambda はテンソルで、熱伝導率テンソル
  出力引数:
    OUT__K_cnd_InJn はスカラーで、
      熱伝導行列の節点 In 節点 Jn 成分
*/
{
  double whV_tmp[2];

  WH_Tensor2D__whT_cdot_whV_OUT_whV
    (whT_lambda, whV_N_X_Jn,
     whV_tmp);
  WH_Tensor2D__whV_cdot_whV_OUT_s
    (whV_N_X_In, whV_tmp,
     OUT__K_cnd_InJn);
}

/* 熱伝導行列：３次元 */
void WH_ThermalFemLinear3D__K_cnd_InJn
(double whV_N_X_In[3], double whV_N_X_Jn[3],
 double whT_lambda[3][3],
 double *OUT__K_cnd_InJn)
/*
  入力引数:
    whV_N_X_In, whV_N_X_Jn はベクトルで、
      それぞれ節点 In, Jn の補間関数勾配ベクトル
    whT_lambda はテンソルで、熱伝導率テンソル
  出力引数:
    OUT__K_cnd_InJn はスカラーで、
      熱伝導行列の節点 In 節点 Jn 成分
*/
{
  double whV_tmp[3];

  WH_Tensor3D__whT_cdot_whV_OUT_whV
    (whT_lambda, whV_N_X_Jn,
     whV_tmp);
  WH_Tensor3D__whV_cdot_whV_OUT_s
    (whV_N_X_In, whV_tmp,
     OUT__K_cnd_InJn);
}

さらに、 熱伝導行列全体を数値積分により求める 関数は、２次元と３次元問題についてそれぞれ以下のようになる。 なお、ここでは、 熱伝導率は要素内において一定であると仮定する。

#define NNODES 6

/* 熱伝導行列：２次元 */
void WH_ThermalFemLinear2D__K_cnd_arrayInJn
(double whV_X_arrayIn[NNODES][2],
 double whT_lambda[2][2],
 double a,
 double OUT__K_cnd_arrayInJn[NNODES][NNODES])
/*
  入力引数:
    whV_X_arrayIn はベクトルの節点配列で、位置（節点座標）
    whT_lambda はテンソルで、熱伝導率テンソル
    a はスカラーで、厚み
  出力引数:
    OUT__K_cnd_arrayInJn はスカラーの配列(節点、節点)で、
      熱伝導行列
*/
{
  int nPoints;
  int Ip;
  int In, Jn;
  const int order = 2;
 
  for (In = 0; In < NNODES; In++) {
    for (Jn = 0; Jn < NNODES; Jn++) {
      OUT__K_cnd_arrayInJn[In][Jn] = 0;
    }
  }

  nPoints = WH_Gauss2D__Triangle__nIntegrationPoints (order);
  for (Ip = 0; Ip < nPoints; Ip++) {
    double whV_xi_Ip[2];
    double w_Ip;
    double whP_N_whV_xi_arrayIn[NNODES][2];
    double whPt_whV_X_whV_xi[2][2];
    double J_S;
    double whV_N_X_arrayIn[NNODES][2];
    
    WH_Gauss2D__Triangle__whV_xi_Ip
     (order, Ip,
      whV_xi_Ip);
    WH_Gauss2D__Triangle__w_Ip
      (order, Ip,
       &w_Ip);
    WH_Shape__Triangle6Nodes__whP_N_whV_xi_arrayIn
     (whV_xi_Ip,
      whP_N_whV_xi_arrayIn);

    WH_Element2D__Volume__whPt_whV_x_whV_xi
      (NNODES, whP_N_whV_xi_arrayIn, whV_X_arrayIn,
       whPt_whV_X_whV_xi);
    WH_Element2D__Volume__J_S
      (whPt_whV_X_whV_xi, a,
        &J_S);
    WH_Element2D__whP_N_whV_x_arrayIn
      (whPt_whV_X_whV_xi, NNODES, whP_N_whV_xi_arrayIn,
       whV_N_X_arrayIn);

    for (In = 0; In < NNODES; In++) {
      for (Jn = 0; Jn < NNODES; Jn++) {
        double K_cnd;
      
        WH_ThermalFemLinear2D__K_cnd_InJn
          (whV_N_X_arrayIn[In], whV_N_X_arrayIn[Jn],
           whT_lambda,
           &K_cnd);
        OUT__K_cnd_arrayInJn[In][Jn] 
          += K_cnd * J_S * w_Ip;
      }
    }
  }
}