如何计算出有限元函数在某一点的值
令为的局部坐标，
为参考坐标，从参考单元到局部单元的坐标变换可表示为：
$$
\begin{gathered}
\begin{pmatrix} x\ y \ z\end{pmatrix} =
\begin{pmatrix}
x_2-x_1 & x_3-x_1 & x_4-x_1 \
y_2-y_1 & y_3-y_1 & y_4-y_1 \
z_2-z_1 & z_3-z_1 & z_4-z_1
\end{pmatrix}
\begin{pmatrix}
\hat x\ \hat y \\hat z
\end{pmatrix}+
\begin{pmatrix}
x_1\x_2\x_3
\end{pmatrix}
\end{gathered}
$$
可以简写成，反过来可以写成
。
局部坐标系下的的任意一函数，我们都可以写成关于
的多项式，这个多项式也可以用
表示。我们知道，对应参考坐标
和局部坐标系
是相对应的，即
。根据
可以计算出多项式的系数。上式左边的矩阵G可以提前算出，记为
，右端项记为
。
回到原来的问题，如何计算出有限元函数在某一点的值？
- 获取坐标
，函数值
- 计算矩阵
，再根据给定坐标
，计算
。
- 用
计算出基函数(1 x y z xy xz yz x^2 y^2 z^2)
- 根据给定的
和矩阵
计算出系数
- 基函数和系数做点积。
代码
首先要有一个计算基函数的函数
std::vector<double>& evaluate_basis_at_point(std::vecor<double> &point){
std::vector<double> basis(10);
basis[0] = 1.0;
basis[1] = point[0];
basis[2] = point[1];
basis[3] = point[2];
basis[4] = point[0]*point[1];
basis[5] = point[0]*point[2];
basis[6] = point[1]*point[2];
basis[7] = point[0]*point[0];
basis[8] = point[1]*point[1];
basis[9] = point[2]*point[2];
return basis;
}
1. 计算
假设我们的局部坐标是4个3维的向量，用matlab进行符号计算：
syms x0 x1 x2 x3;
syms y0 y1 y2 y3;
syms z0 z1 z2 z3;
F = [x0-x1 x0-x2 x0-x3;
y0-y1 y0-y2 y0-y3;
z0-z1 z0-z2 z0-z3];
simplify(inv(F))
化简可以得到：
det = (x0*y1*z2 - x0*y2*z1 - x1*y0*z2 + x1*y2*z0 + x2*y0*z1 - x2*y1*z0 - x0*y1*z3 + x0*y3*z1 + x1*y0*z3 - x1*y3*z0 - x3*y0*z1 + x3*y1*z0 + x0*y2*z3 - x0*y3*z2 - x2*y0*z3 + x2*y3*z0 + x3*y0*z2 - x3*y2*z0 - x1*y2*z3 + x1*y3*z2 + x2*y1*z3 - x2*y3*z1 - x3*y1*z2 + x3*y2*z1)
[ (y0*z2 - y2*z0 - y0*z3 + y3*z0 + y2*z3 - y3*z2)/det,
-(x0*z2 - x2*z0 - x0*z3 + x3*z0 + x2*z3 - x3*z2)/det,
(x0*y2 - x2*y0 - x0*y3 + x3*y0 + x2*y3 - x3*y2)/det]
[ -(y0*z1 - y1*z0 - y0*z3 + y3*z0 + y1*z3 - y3*z1)/det,
(x0*z1 - x1*z0 - x0*z3 + x3*z0 + x1*z3 - x3*z1)/det,
-(x0*y1 - x1*y0 - x0*y3 + x3*y0 + x1*y3 - x3*y1)/det]
[ (y0*z1 - y1*z0 - y0*z2 + y2*z0 + y1*z2 - y2*z1)/det,
-(x0*z1 - x1*z0 - x0*z2 + x2*z0 + x1*z2 - x2*z1)/det,
(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1)/det]
写成C++代码是
double det = (x[0]*y[1]*z[2] - x[0]*y[2]*z[1] - x[1]*y[0]*z[2] + x[1]*y[2]*z[0] + x[2]*y[0]*z[1] - x[2]*y[1]*z[0] - x[0]*y[1]*z[3] + x[0]*y[3]*z[1] + x[1]*y[0]*z[3] - x[1]*y[3]*z[0] - x[3]*y[0]*z[1] + x[3]*y[1]*z[0] + x[0]*y[2]*z[3] - x[0]*y[3]*z[2] - x[2]*y[0]*z[3] + x[2]*y[3]*z[0] + x[3]*y[0]*z[2] - x[3]*y[2]*z[0] - x[1]*y[2]*z[3] + x[1]*y[3]*z[2] + x[2]*y[1]*z[3] - x[2]*y[3]*z[1] - x[3]*y[1]*z[2] + x[3]*y[2]*z[1])
double F[3][3] ={
{ (y[0]*z[2] - y[2]*z[0] - y[0]*z[3] + y[3]*z[0] + y[2]*z[3] - y[3]*z[2])/det,
-(x[0]*z[2] - x[2]*z[0] - x[0]*z[3] + x[3]*z[0] + x[2]*z[3] - x[3]*z[2])/det,
(x[0]*y[2] - x[2]*y[0] - x[0]*y[3] + x[3]*y[0] + x[2]*y[3] - x[3]*y[2])/det},
{-(y[0]*z[1] - y[1]*z[0] - y[0]*z[3] + y[3]*z[0] + y[1]*z[3] - y[3]*z[1])/det,
(x[0]*z[1] - x[1]*z[0] - x[0]*z[3] + x[3]*z[0] + x[1]*z[3] - x[3]*z[1])/det,
-(x[0]*y[1] - x[1]*y[0] - x[0]*y[3] + x[3]*y[0] + x[1]*y[3] - x[3]*y[1])/det},
{ (y[0]*z[1] - y[1]*z[0] - y[0]*z[2] + y[2]*z[0] + y[1]*z[2] - y[2]*z[1])/det,
-(x[0]*z[1] - x[1]*z[0] - x[0]*z[2] + x[2]*z[0] + x[1]*z[2] - x[2]*z[1])/det,
(x[0]*y[1] - x[1]*y[0] - x[0]*y[2] + x[2]*y[0] + x[1]*y[2] - x[2]*y[1])/det}
}
2. 计算
std::vector<double>& evaluate_basis_at_point(std::vecor<double> &point, double **F){
std::vector<double> hat_point(3);
for(size_t i=0;i<3;i++){
hat_point[0] += F[0][i]*point[i];
hat_point[1] += F[1][i]*point[i];
hat_point[2] += F[2][i]*point[i];
}
for(size_t i=0;i<3;i++){
hat_point[0] -= F[0][i]*point1[i];
hat_point[1] -= F[1][i]*point1[i];
hat_point[2] -= F[2][i]*point1[i];
}
return hat_point;
}
3. 
std::vector<double> p0{0.0,0.0,0.0};
std::vector<double> p1{0.5,0.0,0.0};
std::vector<double> p2{1.0,0.0,0.0};
std::vector<double> p3{0.0,0.5,0.0};
std::vector<double> p4{0.0,1.0,0.0};
std::vector<double> p5{0.5,0.5,0.0};
std::vector<double> p6{0.0,0.0,0.5};
std::vector<double> p7{0.0,0.5,0.5};
std::vector<double> p8{0.5,0.0,0.5};
std::vector<double> p9{0.0,0.0,1.0};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
auto b1 = evaluate_basis_at_point{p1};
// 利用eigen求逆得到一个10*10的矩阵G_inv,以备用。
4.基函数
// hat_point 是计算出来的参考点
// B 是各个自由度上的值
b = evaluate_basis_at_point(hat_point);
parameters = multiply(B, G_inv);
result = multiply(b, parameters);
// 如果B是多维的，那么就有多个B。
网友评论