実装：マトリクス標記

要素熱伝導行列を数値積分により求める 関数は、以下のようになる。

/* 要素熱伝導行列：三次元 */
void WH_ThermalFem3D__Linear__K_cnd_arrayInJn_AS_matrix
(WH_Fem__Shape3DType shapeType,
 WH_Fem__Element3DType elementType,
 int nNodes,
 double whV_x_arrayIn[/* nNodes */][3],
 double whM_lambda[3][3],
 double OUT__K_cnd_arrayInJn[/* nNodes * nNodes */])
/*
  入力引数:
    shapeType は要素形状のタイプ
    elementType は要素のタイプ
    nNodes は整数で、要素節点数
    whV_x_arrayIn はベクトルの節点配列で、位置（節点座標）
    whM_lambda は行列で、熱伝導率の行列
  出力引数:
    OUT__K_cnd_arrayInJn はスカラーの配列(節点、節点)で、
      要素熱伝導行列
*/
{
  int order;
  int nPoints;
  int Ip;
  int In, Jn;
 
  switch (elementType) {
  case WH_FEM__ELEMENT3D__TETRA4N:    
    order = 1; 
    break;
  case WH_FEM__ELEMENT3D__TETRA10N:
  case WH_FEM__ELEMENT3D__HEXA8N:
    order = 2;
    break;
  case WH_FEM__ELEMENT3D__HEXA20N:
  case WH_FEM__ELEMENT3D__HEXA27N:
    order = 3;
    break;
  default:
    assert(0);
    break;
  }

  for (In = 0; In < nNodes; In++) {
    for (Jn = 0; Jn < nNodes; Jn++) {
      OUT__K_cnd_arrayInJn[In * nNodes + Jn] = 0.0;
    }
  }

  WH_Fem__NumInt__Gauss3D__nPoints 
    (shapeType, order,
     &nPoints);
  for (Ip = 0; Ip < nPoints; Ip++) {
    double whV_xi_Ip[3];
    double w_Ip;
    double whP_N_whV_xi_arrayIn[WH_FEM__MAX_ELEMENT_NODES][3];
    double whPt_whV_x_whV_xi[3][3];
    double J_V;
    double whCv_N_x_arrayIn[WH_FEM__MAX_ELEMENT_NODES][3];
    
    WH_Fem__NumInt__Gauss3D__whV_xi_Ip
      (shapeType, order, Ip,
       whV_xi_Ip);
    WH_Fem__NumInt__Gauss3D__w_Ip
      (shapeType, order, Ip,
       &w_Ip);
    WH_Fem__Shape3D__whP_N_whV_xi_arrayIn
      (elementType, nNodes, whV_xi_Ip,
       whP_N_whV_xi_arrayIn);

    WH_Fem__Isoparam3D__Volume__whPt_whV_x_whV_xi
      (nNodes, whP_N_whV_xi_arrayIn, whV_x_arrayIn,
       whPt_whV_x_whV_xi);
    WH_Fem__Isoparam3D__Volume__J_V
      (shapeType, whPt_whV_x_whV_xi,
       &J_V);
    WH_Fem__Isoparam3D__whP_N_whV_x_arrayIn
      (whPt_whV_x_whV_xi, nNodes, whP_N_whV_xi_arrayIn,
       whCv_N_x_arrayIn);
    /* column vector (3x1) ... 3-D vector */

    for (In = 0; In < nNodes; In++) {
      for (Jn = 0; Jn < nNodes; Jn++) {
        double whCv_tmp[3];
        double K_cnd;

        WH_Matrix__whM_mul_whCv_OUT_whCv
          (3, 3, whM_lambda, 
           3, whCv_N_x_arrayIn[Jn],
           3, whCv_tmp);
        WH_Matrix__whTrans_whCv_mul_whCv_OUT_s
          (3, whCv_N_x_arrayIn[In], 
           3, whCv_tmp,
           &K_cnd);
        OUT__K_cnd_arrayInJn[In * nNodes + Jn] 
          += K_cnd * J_V * w_Ip;
      }
    }
  }
}