版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/walkandthink/article/details/73351565
一维形函数核心的数学知识就是Lagrange插值.在有限元计算中,我们常常提到的一阶单元,二阶单元就来源于形函数。
比如一维线性单元,因为是线性单元,所以每个单元只需要两个节点就可以了,两个点对应的Lagrange插值多项式,就可以写为如下形式:
将其写为自然坐标
可以验证:
同理,二阶的形函数(每个单元需要3个节点)可以写为:
可以验证:
可以画一下二次形函数的取值,代码如下:
import numpy as np
import matplotlib.pyplot as plt
xi=np.linspace(-1.0,1.0,100)
y0=xi*0.0
N1=xi*(xi-1)/2.0
N2=(1+xi)*(1-xi)
N3=xi*(xi+1)/2.0
plt.plot(xi,N1)
plt.plot(xi,N2)
plt.plot(xi,N3)
plt.plot(xi,y0,'k--')
plt.ylabel('$N_{i}$',fontsize=20)
plt.xlabel('$xi$',fontsize=20)
plt.title('$N_{i}(p=2)$',fontsize=21)
plt.savefig('shp1d.png',dpi=400)
结果跟我们预计的一样:
一般,我们常用的是一阶,二阶的形函数,更高阶的形函数可以根据Lagrange插值多项式来构造。而写成自然坐标是为了以后方便用高斯积分进行积分运算。
在进行数值模拟的过程中,往往会涉及到梯度项还有拉普拉斯项等高阶导数的运算,所以除了形函数本身之外,我们还需要考虑下
方法很简单,利用链式法则就可以了。
比如我们知道某个单元上所有节点的位移
那么位移的梯度(1D)就可以这么求:
为了得到位移的梯度,我们只需要计算
所以就可以很容易地得到:
而
所以
在AsFem中,我将
最后,我们可以得到形函数对
更高阶的偏导,比如拉普拉斯算符,可以利用同样的方法得到,这里不再重复。
在AsFem中对应的实现代码如下:
void Shp1D(const int &nDims,const double &xi,const MatrixXd &Coords,MatrixXd &shp, double &DetJac)
{
int nNodes;
double dxdxi, dydxi, dzdxi;
int i, j;
nNodes = int(Coords.rows());
if (nDims<1||nDims>3)
{
cout << "Wrong dimension case !!!" << endl;
abort();
}
if (nNodes<2 || nNodes>7)
{
cout << "Wrong number of nodes per element!!!" << endl;
cout<<"AsFem only support 1~6 order 1D element!!!"<<endl;
abort();
}
shp.setZero();
if (nNodes == 2)
{
// Linear element
shp(0, 0) = 0.5*(1.0 - xi);
shp(0, 1) =-0.5;
shp(1, 0) = 0.5*(xi + 1.0);
shp(1, 1) = 0.5;
}
else if (nNodes == 3)
{
// Quadratic line element(3 nodes)
shp(0, 0) = 0.5*xi*(xi - 1.0);
shp(0, 1) = 0.5*(2.0*xi - 1.0);
shp(1, 0) = -(xi + 1.0)*(xi - 1.0);
shp(1, 1) = -2.0*xi;
shp(2, 0) = 0.5*xi*(xi + 1.0);
shp(2, 1) = 0.5*(2.0*xi + 1.0);
}
// Compute the derivatives over x
// dN/dx
// actually, there should be positive and negative
// but I want user to decide the direction of norm
// so the DetJac is always positive here
if (nDims == 1)
{
dxdxi = 0.0;
for (i = 0; i < nNodes; i++)
{
dxdxi += Coords(i, 0)*shp(i, 1);
}
DetJac = fabs(dxdxi);
}
else if (nDims == 2)
{
// go to asfem
}
else if (nDims == 3)
{
// go to asfem
}
if (DetJac == 0.0)
{
cout << "Singular in element!!!" << endl;
abort();
}
for (i = 0; i < nNodes; i++)
{
shp(i,1)=shp(i,1)/DetJac;
}
}