実装

/* 領域での積分点数 */
void WH_Fem__NumInt__Gauss3D__Hexa__nPoints 
(int order,
 int *OUT__nPoints)
/*
  入力引数:
    order は整数で、ガウス積分の次数
  出力引数:
    OUT__nPoints は整数で、積分点数
*/
{
  int nLinePoints;

  WH_Fem__NumInt__Gauss1D__nPoints 
    (order,
     &nLinePoints);
  *OUT__nPoints = nLinePoints * nLinePoints * nLinePoints;
}

/* 領域での積分点の自然座標 */
void WH_Fem__NumInt__Gauss3D__Hexa__whV_xi_Ip
(int order, 
 int Ip,
 double OUT__whV_xi_Ip[3])
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Ip は整数で、積分点ID
  出力引数:
    OUT__whV_xi_Ip はベクトルで、自然座標
*/
{
  int nLinePoints;
  int Ip_xi, Ip_eta, Ip_zeta;
  double xi, eta, zeta;

  WH_Fem__NumInt__Gauss1D__nPoints 
    (order,
     &nLinePoints);
  Ip_xi = (Ip % (nLinePoints * nLinePoints)) % nLinePoints;
  Ip_eta = (Ip % (nLinePoints * nLinePoints)) / nLinePoints;
  Ip_zeta = Ip / (nLinePoints * nLinePoints);
  WH_Fem__NumInt__Gauss1D__xi_Ip
    (order, Ip_xi,
     &xi);
  WH_Fem__NumInt__Gauss1D__xi_Ip
    (order, Ip_eta,
     &eta);
  WH_Fem__NumInt__Gauss1D__xi_Ip
    (order, Ip_zeta,
     &zeta);
  OUT__whV_xi_Ip[0] = xi;
  OUT__whV_xi_Ip[1] = eta;
  OUT__whV_xi_Ip[2] = zeta;
}

/* 領域での積分点の重み係数 */
void WH_Fem__NumInt__Gauss3D__Hexa__w_Ip
(int order, 
 int Ip,
 double *OUT__w_Ip)
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Ip は整数で、積分点ID
  出力引数:
    OUT__w_Ip はスカラーで、重み係数
*/
{
  int nLinePoints;
  int Ip_xi, Ip_eta, Ip_zeta;
  double w_xi, w_eta, w_zeta;

  WH_Fem__NumInt__Gauss1D__nPoints 
    (order,
     &nLinePoints);
  Ip_xi = (Ip % (nLinePoints * nLinePoints)) % nLinePoints;
  Ip_eta = (Ip % (nLinePoints * nLinePoints)) / nLinePoints;
  Ip_zeta = Ip / (nLinePoints * nLinePoints);
  WH_Fem__NumInt__Gauss1D__w_Ip
    (order, Ip_xi,
     &w_xi);
  WH_Fem__NumInt__Gauss1D__w_Ip
    (order, Ip_eta,
     &w_eta);
  WH_Fem__NumInt__Gauss1D__w_Ip
    (order, Ip_zeta,
     &w_zeta);
  *OUT__w_Ip = w_xi * w_eta * w_zeta;
}

/* 境界表面上での積分点数 */
void WH_Fem__NumInt__Gauss3D__Hexa__nPoints_Iface
(int order,
 int Iface,
 int *OUT__nPoints)
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Iface は整数で、要素面ID
  出力引数:
    OUT__nPoints は整数で、積分点数
*/
{
  WH_Fem__NumInt__Gauss2D__Quad__nPoints 
    (order,
     OUT__nPoints);
}

/* 境界表面上での積分点の自然座標 */
void WH_Fem__NumInt__Gauss3D__Hexa__whV_xi_IfaceIp
(int order, 
 int Iface,
 int Ip,
 double OUT__whV_xi_Ip[3])
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Iface は整数で、要素面ID
    Ip は整数で、積分点ID
  出力引数:
    OUT__whV_xi_Ip はベクトルで、自然座標
*/
{
  double whV_xi_quad[2];

  WH_Fem__NumInt__Gauss2D__Quad__whV_xi_Ip
    (order, Ip,
     whV_xi_quad);
  switch (Iface) {
  case 5:   /* -xi */
    OUT__whV_xi_Ip[0] = -1.0;
    OUT__whV_xi_Ip[1] = whV_xi_quad[0];
    OUT__whV_xi_Ip[2] = whV_xi_quad[1];
    break;
  case 3:   /* +xi */
    OUT__whV_xi_Ip[0] = 1.0;
    OUT__whV_xi_Ip[1] = whV_xi_quad[0];
    OUT__whV_xi_Ip[2] = whV_xi_quad[1];
    break;
  case 2:   /* -eta */
    OUT__whV_xi_Ip[0] = whV_xi_quad[1];
    OUT__whV_xi_Ip[1] = -1.0;
    OUT__whV_xi_Ip[2] = whV_xi_quad[0];
    break;
  case 4:   /* +eta */
    OUT__whV_xi_Ip[0] = whV_xi_quad[1];
    OUT__whV_xi_Ip[1] = 1.0;
    OUT__whV_xi_Ip[2] = whV_xi_quad[0];
    break;
  case 0:   /* -zeta */
    OUT__whV_xi_Ip[0] = whV_xi_quad[0];
    OUT__whV_xi_Ip[1] = whV_xi_quad[1];
    OUT__whV_xi_Ip[2] = -1.0;
    break;
  case 1:   /* +zeta */
    OUT__whV_xi_Ip[0] = whV_xi_quad[0];
    OUT__whV_xi_Ip[1] = whV_xi_quad[1];
    OUT__whV_xi_Ip[2] = 1.0;
    break;
  default:
    assert(0);
    break;
  }
}

/* 境界表面上での積分点の重み係数 */
void WH_Fem__NumInt__Gauss3D__Hexa__w_IfaceIp
(int order, 
 int Iface,
 int Ip,
 double *OUT__w_Ip)
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Iface は整数で、要素面ID
    Ip は整数で、積分点ID
  出力引数:
    OUT__w_Ip はスカラーで、重み係数
*/
{
  WH_Fem__NumInt__Gauss2D__Quad__w_Ip
    (order, Ip, 
     OUT__w_Ip);
}