実装

#define MAX_TETRA_ORDER 3
#define MAX_TETRA_POINTS 5

static int Tetra__nPoints[MAX_TETRA_ORDER] = {
  1, 4, 5
};

#define TETRA_A2  0.13819660

static double Tetra__L0123_Ip
[MAX_TETRA_ORDER][MAX_TETRA_POINTS][4] = {
  { { 1.0/4, 1.0/4, 1.0/4, 1.0/4 } },
  { { 1.0 - 2 * TETRA_A2, TETRA_A2, TETRA_A2, TETRA_A2 },
    { TETRA_A2, 1.0 - 2 * TETRA_A2, TETRA_A2, TETRA_A2 },
    { TETRA_A2, TETRA_A2, 1.0 - 2 * TETRA_A2, TETRA_A2 },
    { TETRA_A2, TETRA_A2, TETRA_A2, 1.0 - 2 * TETRA_A2 } },
  { { 1.0/4, 1.0/4, 1.0/4, 1.0/4 },
    { 1.0/2, 1.0/6, 1.0/6, 1.0/6 },
    { 1.0/6, 1.0/2, 1.0/6, 1.0/6 },
    { 1.0/6, 1.0/6, 1.0/2, 1.0/6 },
    { 1.0/6, 1.0/6, 1.0/6, 1.0/2 } }
};

static double Tetra__w_Ip
[MAX_TETRA_ORDER][MAX_TETRA_POINTS] = {
  {  1.0 },
  {  1.0/4, 1.0/4, 1.0/4, 1.0/4 },
  {  -16.0/20, 
       9.0/20, 9.0/20, 9.0/20, 9.0/20 }
};

/* 領域での積分点数 */
void WH_Fem__NumInt__Gauss3D__Tetra__nPoints 
(int order,
 int *OUT__nPoints)
/*
  入力引数:
    order は整数で、ガウス積分の次数
  出力引数:
    OUT__nPoints は整数で、積分点数
*/
{
  assert(0 < order);
  assert(order <= MAX_TETRA_ORDER);

  *OUT__nPoints = Tetra__nPoints[order - 1];
}

/* 領域での積分点の自然座標 */
void WH_Fem__NumInt__Gauss3D__Tetra__whV_xi_Ip
(int order, 
 int Ip,
 double OUT__whV_xi_Ip[3])
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Ip は整数で、積分点ID
  出力引数:
    OUT__whV_xi_Ip はベクトルで、自然座標
*/
{
  double L0, L1, L2, L3;

  assert(0 < order);
  assert(order <= MAX_TETRA_ORDER);
  assert(0 <= Ip);
  assert(Ip < Tetra__nPoints[order - 1]);

  L0 = Tetra__L0123_Ip[order - 1][Ip][0];
  L1 = Tetra__L0123_Ip[order - 1][Ip][1];
  L2 = Tetra__L0123_Ip[order - 1][Ip][2];
  L3 = Tetra__L0123_Ip[order - 1][Ip][3];

  OUT__whV_xi_Ip[0] = L1;
  OUT__whV_xi_Ip[1] = L2;
  OUT__whV_xi_Ip[2] = L3;
}

/* 領域での積分点の重み係数 */
void WH_Fem__NumInt__Gauss3D__Tetra__w_Ip
(int order, 
 int Ip,
 double *OUT__w_Ip)
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Ip は整数で、積分点ID
  出力引数:
    OUT__w_Ip はスカラーで、重み係数
*/
{
  assert(0 < order);
  assert(order <= MAX_TETRA_ORDER);
  assert(0 <= Ip);
  assert(Ip < Tetra__nPoints[order - 1]);

  *OUT__w_Ip = Tetra__w_Ip[order - 1][Ip];
}

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

/* 境界表面上での積分点の自然座標 */
void WH_Fem__NumInt__Gauss3D__Tetra__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_tri[2];
  double L0_tri, L1_tri, L2_tri;
  double L0, L1, L2, L3;

  WH_Fem__NumInt__Gauss2D__Tri__whV_xi_Ip
    (order, Ip,
     whV_xi_tri);
  L0_tri = 1.0 - whV_xi_tri[0] - whV_xi_tri[1];
  L1_tri = whV_xi_tri[0];
  L2_tri = whV_xi_tri[1];

  switch (Iface) {
  case 0:   /* L0 */
    L0 = 0.0;
    L1 = L0_tri;
    L2 = L1_tri;
    L3 = L2_tri;
    break;
  case 1:   /* L1 */
    L0 = L2_tri;
    L1 = 0.0;
    L2 = L0_tri;
    L3 = L1_tri;
    break;
  case 2:   /* L2 */
    L0 = L1_tri;
    L1 = L2_tri;
    L2 = 0.0;
    L3 = L0_tri;
    break;
  case 3:   /* L3 */
    L0 = L0_tri;
    L1 = L1_tri;
    L2 = L2_tri;
    L3 = 0.0;
    break;
  default:
    assert(0);
    break;
  }

  OUT__whV_xi_Ip[0] = L1;
  OUT__whV_xi_Ip[1] = L2;
  OUT__whV_xi_Ip[2] = L3;
}

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