実装

#define MAX_TRI_ORDER 5
#define MAX_TRI_POINTS 7

static int Tri__nPoints[MAX_TRI_ORDER] = {
  1, 3, 4, 6, 7
};

#define TRI_A4  0.445948490915965
#define TRI_B4  (0.111690794839005 * 2)
#define TRI_C4  0.091576213509771
#define TRI_D4  (0.054975871827661 * 2)

#define TRI_A5  0.470142064105115
#define TRI_B5  (0.0661970763942530 * 2)
#define TRI_C5  0.101286507323456
#define TRI_D5  (0.0629695902724135 * 2)

static double Tri__L012_Ip
[MAX_TRI_ORDER][MAX_TRI_POINTS][3] = {
  { { 1.0/3, 1.0/3, 1.0/3 } },
  { { 0.0/2, 1.0/2, 1.0/2 },
    { 1.0/2, 0.0/2, 1.0/2 },
    { 1.0/2, 1.0/2, 0.0/2 } },
  { { 1.0/3, 1.0/3, 1.0/3 },
    { 11.0/15,  2.0/15,  2.0/15 },
    {  2.0/15, 11.0/15,  2.0/15 },
    {  2.0/15,  2.0/15, 11.0/15 } },
  { { 1.0 - 2 * TRI_A4, TRI_A4, TRI_A4 },
    { TRI_A4, 1.0 - 2 * TRI_A4, TRI_A4 },
    { TRI_A4, TRI_A4, 1.0 - 2 * TRI_A4 },
    { 1.0 - 2 * TRI_C4, TRI_C4, TRI_C4 },
    { TRI_C4, 1.0 - 2 * TRI_C4, TRI_C4 },
    { TRI_C4, TRI_C4, 1.0 - 2 * TRI_C4 } },
  { { 1.0/3, 1.0/3, 1.0/3 },
    { 1.0 - 2 * TRI_A5, TRI_A5, TRI_A5 },
    { TRI_A5, 1.0 - 2 * TRI_A5, TRI_A5 },
    { TRI_A5, TRI_A5, 1.0 - 2 * TRI_A5 },
    { 1.0 - 2 * TRI_C5, TRI_C5, TRI_C5 },
    { TRI_C5, 1.0 - 2 * TRI_C5, TRI_C5 },
    { TRI_C5, TRI_C5, 1.0 - 2 * TRI_C5 } }
};

static double Tri__w_Ip
[MAX_TRI_ORDER][MAX_TRI_POINTS] = {
  {  1.0 },
  {  1.0/3, 1.0/3, 1.0/3 },
  {  -27.0/48, 
      25.0/48, 25.0/48, 25.0/48  },
  {  TRI_B4, TRI_B4, TRI_B4,
     TRI_D4, TRI_D4, TRI_D4  },
  {  9.0/40, 
     TRI_B5, TRI_B5, TRI_B5,
     TRI_D5, TRI_D5, TRI_D5  }
};

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

  *OUT__nPoints = Tri__nPoints[order - 1];
}

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

  assert(0 < order);
  assert(order <= MAX_TRI_ORDER);
  assert(0 <= Ip);
  assert(Ip < Tri__nPoints[order - 1]);

  L0 = Tri__L012_Ip[order - 1][Ip][0];
  L1 = Tri__L012_Ip[order - 1][Ip][1];
  L2 = Tri__L012_Ip[order - 1][Ip][2];

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

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

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

/* 境界稜線上での積分点数 */
void WH_Fem__NumInt__Gauss2D__Tri__nPoints_Iedge
(int order,
 int Iedge,
 int *OUT__nPoints)
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Iedge は整数で、要素辺ID
  出力引数:
    OUT__nPoints は整数で、積分点数
*/
{
  WH_Fem__NumInt__Gauss1D__nPoints 
    (order,
     OUT__nPoints);
}

/* 境界稜線上での積分点の自然座標 */
void WH_Fem__NumInt__Gauss2D__Tri__whV_xi_IedgeIp
(int order, 
 int Iedge,
 int Ip,
 double OUT__whV_xi_Ip[2])
/*
  入力引数:
    order は整数で、ガウス積分の次数
    Iedge は整数で、要素辺ID
    Ip は整数で、積分点ID
  出力引数:
    OUT__whV_xi_Ip はベクトルで、自然座標
*/
{
  double xi_line;
  double L0_line, L1_line;
  double L0, L1, L2;

  WH_Fem__NumInt__Gauss1D__xi_Ip
    (order, Ip,
     &xi_line);
  L0_line = 1.0 - xi_line;
  L1_line = xi_line;

  switch (Iedge) {
  case 0:   /* L0 */
    L0 = 0.0;
    L1 = L0_line;
    L2 = L1_line;
    break;
  case 1:   /* L1 */
    L0 = L1_line;
    L1 = 0.0;
    L2 = L0_line;
    break;
  case 2:   /* L2 */
    L0 = L0_line;
    L1 = L1_line;
    L2 = 0.0;
    break;
  default:
    assert(0);
    break;
  }

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

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