実装

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

/* 剛性行列：２次元 */
void WH_SolidFemLinear2D__K_InIdJnJd
(double whT_N_epsilon_InId[2][2], 
 double whT_N_epsilon_JnId[2][2],
 double whT4_C_e[2][2][2][2],
 double *OUT__K_InIdJnJd)
/*
  入力引数:
    whT_N_epsilon_InId, whT_N_epsilon_JnJd はテンソルで、
      それぞれ節点 In 自由度 Id 、および、
      節点 Jn 自由度 Jd の歪みの補間関数勾配テンソル
    whT4_C_e は４階テンソルで、弾性定数テンソル
  出力引数:
    OUT__K_InIdJnJd はスカラーで、
      剛性行列の節点 In 自由度 Id 節点 Jn 自由度 Jd 成分
*/
{
  double whT_tmp[2][2];

  WH_Tensor2D__whT4_colon_whT_OUT_whT
    (whT4_C_e, whT_N_epsilon_JnJd,
     whT_tmp);
  WH_Tensor2D__whT_colon_whT_OUT_s
    (whT_tmp, whT_N_epsilon_InId,
     OUT__K_InIdJnJd);
}

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

#define NNODES 6
#define NDOFS 2

/* 剛性行列：２次元 */
void WH_SolidFemLinear2D__K_arrayInIdJnJd
(double whV_X_arrayIn[NNODES][2],
 double whT4_C_e[2][2][2][2],
 double a,
 double OUT__K_arrayInIdJnJd[NNODES * NDOFS][NNODES * NDOFS])
/*
  入力引数:
    whV_X_arrayIn はベクトルの節点配列で、位置（節点座標）
    whT4_C_e は４階テンソルで、弾性定数テンソル
    a はスカラーで、厚み
  出力引数:
    OUT__K_arrayInIdJnJd はスカラーの
      配列(節点、自由度、節点、自由度)で、剛性行列
*/
{
  int nPoints;
  int Ip;
  int In, Jn;
  int Id, Jd;
  const int order = 2;
 
  for (In = 0; In < NNODES; In++) {
    for (Id = 0; Id < NDOFS; Id++) {
      for (Jn = 0; Jn < NNODES; Jn++) {
        for (Jd = 0; Jd < NDOFS; Jd++) {
          OUT__K_arrayInIdJnJd[In][Id][Jn][Jd] = 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_arrayIn[NNODES][2];
    double whT_N_X_arrayInId[NNODES][NDOFS][2][2];
    double whT_N_epsilon_arrayInId[NNODES][NDOFS][2][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);
    WH_Element2D__whT_N_x_arrayInId
      (NNODES, whP_N_whV_X_arrayIn, 
       whT_N_X_arrayInId);
    WH_SolidFemLinear2D__whT_N_epsilon_arrayInId
      (NNODES, whT_N_X_arrayInId,
       whT_N_epsilon_arrayInId);

    for (In = 0; In < NNODES; In++) {
      for (Id = 0; Id < NDOFS; Id++) {
        for (Jn = 0; Jn < NNODES; Jn++) {
          for (Jd = 0; Jd < NDOFS; Jd++) {
            double K;
      
            WH_SolidFemLinear2D__K_InIdJnJd
              (whT_N_epsilon_arrayIn[In][Id], 
               whT_N_epsilon_arrayIn[Jn][Jd],
               whT4_C_e,
               &K);
            OUT__K_arrayInIdJnJd[In][Id][Jn][Jd]
              += K * J_S * w_Ip;
          }
        }
      }
    }
  }
}