如何计算出有限元函数在某一点的值

作者: 马鹏飞_47c5 | 来源:发表于2020-10-21 21:31 被阅读0次

如何计算出有限元函数在某一点的值

令 $A_i$ 为的局部坐标， $\hat{A_i}$ 为参考坐标，从参考单元到局部单元的坐标变换可表示为：
$$
\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}
$$
可以简写成 $A=F\hat{A}+A_1$ ，反过来可以写成 $\hat{A}=F^{-1}A-F^{-1}A_1$ 。

局部坐标系下的的任意一函数 $f(A)$ ，我们都可以写成关于 $x,y,z$ 的多项式，这个多项式也可以用 $\hat x, \hat y ,\hat z$ 表示。我们知道，对应参考坐标 $A_i$ 和局部坐标系 $\hat{A_i}$ 是相对应的，即 $f(A_i)=\hat{f}(\hat{A_i})$ 。根据
$\begin{pmatrix} 1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\ 1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\1&\hat x_1&\hat y_1&\hat z_1&\hat x_1\hat y_1&\hat x_1\hat z_1&\hat y_1\hat z_1&\hat x_1^2&\hat y_1^2&\hat z_1^2\\ \end{pmatrix} \begin{pmatrix} a\\b\\c\\d\\e\\f\\g\\h\\i\\j \end{pmatrix} = \begin{pmatrix} f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\ f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\f(\hat{A_i})\\ \end{pmatrix}$
可以计算出多项式的系数 $a,b,c,d,e,f,g,h,i,j$ 。上式左边的矩阵G可以提前算出，记为 $G^{-1}$ ，右端项记为 $B$ 。

回到原来的问题，如何计算出有限元函数在某一点的值？

获取坐标 $A_i$ ，函数值 $f(A_i)$
计算矩阵 $F^{-1}$ ，再根据给定坐标 $A$ ，计算 $\hat{A}=F^{-1}A-F^{-1}A_1$ 。
用 $\hat A$ 计算出基函数(1 x y z xy xz yz x^2 y^2 z^2)
根据给定的 $B$ 和矩阵 $G^{-1}$ 计算出系数
基函数和系数做点积。

代码

首先要有一个计算基函数的函数

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. 计算 $F^{-1}$

假设我们的局部坐标 $A_i$ 是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. 计算 $\hat{A}$

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. $\hat A_i$

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。

如何计算出有限元函数在某一点的值

代码

1. 计算 $F^{-1}$

2. 计算 $\hat{A}$

3. $\hat A_i$

4.基函数

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读