位置（全体座標）は、配列変数として実装される。
/* 位置変数 whV_x はベクトルで、 例えば以下のように定義される。 double whV_x[3]; whV_x[0] は x, whV_x[1] は y, whV_x[2] は zを表す。 */
体積座標から位置(全体座標)を求める 関数は以下のようになる。
/* 体積座標に対応する位置 */
void WH_Fem__Simplex3D__whV_x
(double whV_x0[3],
double whV_x1[3],
double whV_x2[3],
double whV_x3[3],
double L0, double L1, double L2, double L3,
double OUT__whV_x[3])
/*
入力引数:
whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
L0, L1, L2, L3 はスカラーで、体積座標
出力引数:
OUT__whV_x はベクトルで、位置
*/
{
OUT__whV_x[0]
= L0 * whV_x0[0] + L1 * whV_x1[0] + L2 * whV_x2[0];
OUT__whV_x[1]
= L0 * whV_x0[1] + L1 * whV_x1[1] + L2 * whV_x2[1];
OUT__whV_x[2]
= L0 * whV_x0[2] + L1 * whV_x1[2] + L2 * whV_x2[2];
}
位置(全体座標)から体積座標を求める 関数は以下のようになる。
/* 係数 a */
void WH_Fem__Simplex3D__a
(double whV_x0[3],
double whV_x1[3],
double whV_x2[3],
double whV_x3[3],
double V,
double *OUT__a0, double *OUT__a1,
double *OUT__a2, double *OUT__a3)
/*
入力引数:
whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
V はスカラーで、四面体の体積
出力引数:
OUT__a0, OUT__a1, OUT__a2, OUT__a3 はスカラーで、係数
*/
{
double x0 = whV_x0[0];
double x1 = whV_x1[0];
double x2 = whV_x2[0];
double x3 = whV_x3[0];
double y0 = whV_x0[1];
double y1 = whV_x1[1];
double y2 = whV_x2[1];
double y3 = whV_x3[1];
double z0 = whV_x0[2];
double z1 = whV_x1[2];
double z2 = whV_x2[2];
double z3 = whV_x3[2];
double factor = 1.0 / (6.0 * V);
*OUT__a0 = factor *
(+ x1 * y2 * z3 + x3 * y1 * z2 + x2 * y3 * z1
- x1 * y3 * z2 - x2 * y1 * z3 - x3 * y2 * z1);
*OUT__a1 = factor *
(+ x0 * y2 * z3 + x3 * y0 * z2 + x2 * y3 * z0
- x0 * y3 * z2 - x2 * y0 * z3 - x3 * y2 * z0);
*OUT__a2 = factor *
(+ x0 * y1 * z3 + x3 * y0 * z1 + x1 * y3 * z0
- x0 * y3 * z1 - x1 * y0 * z3 - x3 * y1 * z0);
*OUT__a3 = factor *
(+ x0 * y1 * z2 + x2 * y0 * z1 + x1 * y2 * z0
- x0 * y2 * z1 - x1 * y0 * z2 - x3 * y1 * z0);
}
/* 係数 b, c, d */
void WH_Fem__Simplex3D__bcd
(double whV_x0[3],
double whV_x1[3],
double whV_x2[3],
double whV_x3[3],
double V,
double *OUT__b0, double *OUT__b1,
double *OUT__b2, double *OUT__b3,
double *OUT__c0, double *OUT__c1,
double *OUT__c2, double *OUT__c3,
double *OUT__d0, double *OUT__d1,
double *OUT__d2, double *OUT__d3)
/*
入力引数:
whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
V はスカラーで、四面体の体積
出力引数:
OUT__b0, OUT__b1, OUT__b2, OUT__b3 はスカラーで、係数 b
OUT__c0, OUT__c1, OUT__c2, OUT__c3 はスカラーで、係数 c
OUT__d0, OUT__d1, OUT__d2, OUT__d3 はスカラーで、係数 d
*/
{
double x0 = whV_x0[0];
double x1 = whV_x1[0];
double x2 = whV_x2[0];
double x3 = whV_x3[0];
double y0 = whV_x0[1];
double y1 = whV_x1[1];
double y2 = whV_x2[1];
double y3 = whV_x3[1];
double z0 = whV_x0[2];
double z1 = whV_x1[2];
double z2 = whV_x2[2];
double z3 = whV_x3[2];
double factor = 1.0 / (6.0 * V);
*OUT__b0 = factor *
(- y1 * z2 + y1 * z3 + y2 * z1
- y2 * z3 - y3 * z1 + y3 * z2);
*OUT__b1 = factor *
(+ y0 * z2 - y0 * z3 - y2 * z0
+ y2 * z3 + y3 * z0 - y3 * z2);
*OUT__b2 = factor *
(- y0 * z1 + y0 * z3 + y1 * z0
- y1 * z3 - y3 * z0 + y3 * z1);
*OUT__b3 = factor *
(+ y0 * z1 - y0 * z2 - y1 * z0
+ y1 * z2 + y2 * z0 - y2 * z1);
*OUT__c0 = factor *
(+ z2 * x1 - z3 * x1 - z1 * x2
+ z3 * x2 + z1 * x3 - z2 * x3);
*OUT__c1 = factor *
(- z2 * x0 + z3 * x0 + z0 * x2
- z3 * x2 - z0 * x3 + z2 * x3);
*OUT__c2 = factor *
(+ z1 * x0 - z3 * x0 - z0 * x1
+ z3 * x1 + z0 * x3 - z1 * x3);
*OUT__c3 = factor *
(- z1 * x0 + z2 * x0 + z0 * x1
- z2 * x1 - z0 * x2 + z1 * x2);
*OUT__d0 = factor *
(- x2 * y1 + x3 * y1 + x1 * y2
- x3 * y2 - x1 * y3 + x2 * y3);
*OUT__d1 = factor *
(+ x2 * y1 - x3 * y1 - x1 * y2
+ x3 * y2 + x1 * y3 - x2 * y3);
*OUT__d2 = factor *
(- x2 * y1 + x3 * y1 + x1 * y2
- x3 * y2 - x1 * y3 + x2 * y3);
*OUT__d3 = factor *
(+ x2 * y1 - x3 * y1 - x1 * y2
+ x3 * y2 + x1 * y3 - x2 * y3);
}
/* 位置(全体座標)に対応する体積座標 */
void WH_Fem__Simplex3D__L0123
(double whV_x0[3],
double whV_x1[3],
double whV_x2[3],
double whV_x3[3],
double whV_x[3],
double *OUT__L0, double *OUT__L1,
double *OUT__L2, double *OUT__L3)
/*
入力引数:
whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
whV_x はベクトルで、位置(全体座標)
出力引数:
OUT__L0, OUT__L1, OUT__L2, OUT__L3 はスカラーで、体積座標
*/
{
double x = whV_x[0];
double y = whV_x[1];
double z = whV_x[2];
double V;
double a0, a1, a2, a3;
double b0, b1, b2, b3;
double c0, c1, c2, c3;
double d0, d1, d2, d3;
WH_Fem__Simplex3D__V
(whV_x0, whV_x1, whV_x2, whV_x3,
&V);
WH_Fem__Simplex3D__a
(whV_x0, whV_x1, whV_x2, whV_x3, V,
&a0, &a1, &a2, &a3);
WH_Fem__Simplex3D__bcd
(whV_x0, whV_x1, whV_x2, whV_x3, V,
&b0, &b1, &b2, &b3,
&c0, &c1, &c2, &c3,
&d0, &d1, &d2, &d3);
*OUT__L0 = a0 + b0 * x + c0 * y + d0 * z;
*OUT__L1 = a1 + b1 * x + c1 * y + d1 * z;
*OUT__L2 = a2 + b2 * x + c2 * y + d2 * z;
*OUT__L3 = a3 + b3 * x + c3 * y + d3 * z;
}
体積座標の勾配を求める 関数は以下のようになる。
/* 体積座標の勾配 */
void WH_Fem__Simplex3D__whP_L0123_whV_x
(double whV_x0[3],
double whV_x1[3],
double whV_x2[3],
double whV_x3[3],
double OUT__whP_L0_whV_x[3],
double OUT__whP_L1_whV_x[3],
double OUT__whP_L2_whV_x[3],
double OUT__whP_L3_whV_x[3])
/*
入力引数:
whV_x0, whV_x1, whV_x2, whV_x3 はベクトルで、各頂点の位置
出力引数:
OUT__whP_L0_whV_x,
OUT__whP_L1_whV_x,
OUT__whP_L2_whV_x,
OUT__whP_L3_whV_x はベクトルで、体積座標の勾配
*/
{
double V;
double b0, b1, b2, b3;
double c0, c1, c2, c3;
double d0, d1, d2, d3;
WH_Fem__Simplex3D__V
(whV_x0, whV_x1, whV_x2, whV_x3,
&V);
WH_Fem__Simplex3D__bcd
(whV_x0, whV_x1, whV_x2, whV_x3, V,
&b0, &b1, &b2, &b3,
&c0, &c1, &c2, &c3,
&d0, &d1, &d2, &d3);
OUT__whP_L0_whV_x[0] = b0;
OUT__whP_L1_whV_x[0] = b1;
OUT__whP_L2_whV_x[0] = b2;
OUT__whP_L3_whV_x[0] = b3;
OUT__whP_L0_whV_x[1] = c0;
OUT__whP_L1_whV_x[1] = c1;
OUT__whP_L2_whV_x[1] = c2;
OUT__whP_L3_whV_x[1] = c3;
OUT__whP_L0_whV_x[2] = d0;
OUT__whP_L1_whV_x[2] = d1;
OUT__whP_L2_whV_x[2] = d2;
OUT__whP_L3_whV_x[2] = d3;
}