形函数

形函数也称插值函数，是构造单元内某点坐标、位移与单元顶点坐标或位移关系的函数。构造时需先假定单元的位移模式，一般用多项式表达，如三节点三角形单元采用 1+u+v，四节点四边形单元用1+u+v+u*v。插值函数是节点坐标与位移模式多项式的函数。

位移模式

我们以三角形单元为例说明形函数的构造方法：

1.三角形单元的位移模式为1+u+v；
2. 三个顶点的坐标分别为(xi，yi),(xj，yj),(xm，ym)(x_{i}，y_{i}),(x_{j}，y_{j}),(x_{m}，y_{m})(xi，yi),(xj，yj),(xm，ym)，对应局部坐标系下的坐标为(0，0),(1，0),(0，1)；
3. 单元内任意点的位移可通过下式表达：
(1)u=α1+α2x+α3yv=α4+α5x+α6yu=\alpha_{1}+ \alpha_{2} x+\alpha_{3}y\newline v=\alpha_{4}+ \alpha_{5} x+\alpha_{6}y \tag{1}u=α1+α2x+α3yv=α4+α5x+α6y(1)
4.将节点坐标带入（1）式，可得顶点处的位移为：
ui=α1+α2xi+α3yiuj=α1+α2xj+α3yjum=α1+α2xm+α3ymvi=α4+α5xi+α6yivj=α4+α5xj+α6yjvm=α4+α5xm+α6ymu_{i}=\alpha_{1}+ \alpha_{2} x_{i}+\alpha_{3}y_{i}\newline u_{j}=\alpha_{1}+ \alpha_{2} x_{j}+\alpha_{3}y_{j}\newline u_{m}=\alpha_{1}+ \alpha_{2} x_{m}+\alpha_{3}y_{m}\newline v_{i}=\alpha_{4}+ \alpha_{5} x_{i}+\alpha_{6}y_{i}\newline v_{j}=\alpha_{4}+ \alpha_{5} x_{j}+\alpha_{6}y_{j}\newline v_{m}=\alpha_{4}+ \alpha_{5} x_{m}+\alpha_{6}y_{m}\newline ui=α1+α2xi+α3yiuj=α1+α2xj+α3yjum=α1+α2xm+α3ymvi=α4+α5xi+α6yivj=α4+α5xj+α6yjvm=α4+α5xm+α6ym
即：
(2){uiujum}=[1xiyi1xjyj1xmym]∗{α1α2α3}\left\{\begin{matrix} u_{i}\\ u_{j} \\ u_{m} \\ \end{matrix}\right\} = \left[ \begin{matrix} 1&x_{i}&y_{i}\\ 1&x_{j}&y_{j}\\ 1&x_{m}&y_{m} \\ \end{matrix} \right] * \left\{\begin{matrix} \alpha_{1}\\ \alpha_{2} \\ \alpha_{3} \\ \end{matrix}\right\}\tag{2} ⎩⎨⎧uiujum⎭⎬⎫=⎣⎡111xixjxmyiyjym⎦⎤∗⎩⎨⎧α1α2α3⎭⎬⎫(2)
令：
[T]=[1xiyi1xjyj1xmym]\left[ \begin{matrix} T \end{matrix} \right] = \left[ \begin{matrix} 1&x_{i}&y_{i}\\ 1&x_{j}&y_{j}\\ 1&x_{m}&y_{m} \\ \end{matrix} \right] [T]=⎣⎡111xixjxmyiyjym⎦⎤

JuliaFEM中[T]矩阵的实现

[T]中“列”为位移模式中多项式在i,j,m处的数值。JuliaFEM中包FEMBasis中的vandermonde.jl主要用于实现[T]矩阵。


'''
vandermonde_matrix(polynomial, coordinates)
Given some polynomial and coordinates points (1-3 dimensions), create a Vandermonde
matrix.
# Example
To genererate a Vandermonde matrix for a reference quadrangle `[-1.0, 1.0]^2` for
polynomial `p(u,v) = 1 + u + v + u*v`, one writes:polynomial = :(1 + u + v + u*v)
coordinates = ((-1.0,-1.0), (1.0,-1.0), (1.0,1.0), (-1.0,1.0))
V = vandermonde_matrix(polynomial, coordinates)
# output
[
1.0 -1.0 -1.0  1.0
1.0  1.0 -1.0 -1.0
1.0  1.0  1.0  1.0
1.0 -1.0  1.0 -1.0
]
'''
function vandermonde_matrix(polynomial::Expr, coordinates::NTuple{N, NTuple{D, T}}) where {N, D, T<:Number}A = zeros(N, N)first(polynomial.args) == :+ || error("Use only summation between terms of polynomial")args = polynomial.args[2:end]for i in 1:NX = coordinates[i]if D == 1data = (:u => X[1],)elseif D == 2data = (:u => X[1], :v => X[2])elseif D == 3data = (:u => X[1], :v => X[2], :w => X[3])endfor (j, term) in enumerate(args)A[i,j] = subs(term, data)endendreturn A
end

代码中subs（）主要实现用数值代替表达式中的变量进行计算；

'''
subs(expression, data)
Given expression and pair(s) of `symbol => value` data, substitute to expression.
# Examples
Let us have polynomial `1 + u + v + u*v^2`, and substitute u=1 and v=2:expression = :(1 + u + v + u*v^2)
data = (:u => 1.0, :v => 2.0)
subs(expression, data)
8.0
'''
function subs(p::Expr, data::NTuple{N,Pair{Symbol, T}}) where {N, T}for di in datap = subs(p, di)endreturn simplify(p)
end

上述三角形单元生产的[T]矩阵为：

julia> p= :(1+u+v)
:(1 + u + v)
julia> X=(
(0.0, 0.0), # i
(1.0, 0.0), # j
(0.0, 1.0), # m
)
((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
julia> V=vandermonde_matrix(p,X)
3×3 Array{Float64,2}:
1.0  0.0  0.0
1.0  1.0  0.0
1.0  0.0  1.0

形函数定义

将（2）进行推导可得：
(3){α1α2α3}=[T]−1∗{uiujum}\left\{ \begin{matrix} \alpha_{1}\\ \alpha_{2} \\ \alpha_{3} \\ \end{matrix} \right\}= \left[ \begin{matrix} T \end{matrix} \right] ^{-1} * \left\{\begin{matrix} u_{i}\\ u_{j} \\ u_{m} \\ \end{matrix}\right\} \tag{3} ⎩⎨⎧α1α2α3⎭⎬⎫=[T]−1∗⎩⎨⎧uiujum⎭⎬⎫(3)
将（3）式带入（1）式可得：
(4)u=[1uv]∗[T]−1∗{uiujum}u= \left[ \begin{matrix} 1&u&v \end{matrix} \right] * \left[ \begin{matrix} T \end{matrix} \right] ^{-1} * \left\{\begin{matrix} u_{i}\\ u_{j} \\ u_{m} \\ \end{matrix}\right\} \tag{4}u=[1uv]∗[T]−1∗⎩⎨⎧uiujum⎭⎬⎫(4)
定义形函数：[N]=[1uv]∗[T]−1\left[ \begin{matrix} N \end{matrix} \right]=\left[ \begin{matrix} 1&u&v \end{matrix} \right] * \left[ \begin{matrix} T \end{matrix} \right] ^{-1} [N]=[1uv]∗[T]−1
这里的坐标为局部坐标，由[T]∗[T]−1=[E]\left[ \begin{matrix} T \end{matrix} \right]* \left[ \begin{matrix} T \end{matrix} \right] ^{-1} = \left[ \begin{matrix} E \end{matrix} \right][T]∗[T]−1=[E]
那么 [T]∗[X1]=[100]\left[ \begin{matrix} T \end{matrix} \right]*[X1]= \left[ \begin{matrix} 1\\ 0\\ 0 \end{matrix} \right][T]∗[X1]=⎣⎡100⎦⎤这里的X1为[T]^{-1}的第一列，同理 [T]∗[X2]=[010]这里的X2为[T]−1的第二列。\left[ \begin{matrix} T \end{matrix} \right]*[X2]= \left[ \begin{matrix} 0\\ 1\\ 0 \end{matrix} \right]这里的X2为[T]^{-1}的第二列。[T]∗[X2]=⎣⎡010⎦⎤这里的X2为[T]−1的第二列。

Julia中形函数的实现

create_basis.jl中的calculate_interpolation_polynomials(p,v)函数用于计算形函数，参数p为位移模式多项式(p::Expr)如:(1+u+v)，参数v即为上述的[T]矩阵。

function calculate_interpolation_polynomials(p, V)basis = []first(p.args) == :+ || error("Use only summation between terms of polynomial")args = p.args[2:end]N = size(V, 1)b = zeros(N)for i in 1:Nfill!(b, 0.0)b[i] = 1.0#计算逆矩阵的单列solution = V \ bN = Expr(:call, :+)for (ai, bi) in zip(solution, args)isapprox(ai, 0.0) && continue#位移模式多项式如[1 u v]与逆矩阵的单列相乘并放入数组中push!(N.args, simplify( :( $ai * $bi ) ))endpush!(basis, N)endreturn basis
end

Julia中形函数导数的实现

主要通过Calculus模块中differentiate(N,var)函数，N为形函数，var为[:u, :v, :w]

function calculate_interpolation_polynomial_derivatives(basis, D)vars = [:u, :v, :w]dbasis = []for N in basispartial_derivatives = []for j in 1:D#分别对u v w进行求导push!(partial_derivatives, simplify(differentiate(N, vars[j])))endpush!(dbasis, partial_derivatives)endreturn dbasis
end

有限元形函数及JuliaFEM中的实现方式相关推荐

形函数的构造原理-有限元形函数的几个种类
在有限元法中,形函数是一个十分重要的概念.它不仅可以用做单元的内插函数,把单元内任一点的位移用节点位移表示,而且可作为加权余量法中的加权函数,可以处理外载荷,将分布力等效为节点上的集中力和力矩,此外, ...
C++成员函数在内存中的存储方式
用类去定义对象时,系统会为每一个对象分配存储空间.如果一个类包括了数据和函数,要分别为数据和函数的代码分配存储空间.按理说,如果用同一个类定义了10个对象,那么就需要分别为10个对象的数据和函数代码分 ...
c++全局类对象_C++ 类在内存中的存储方式(一)
说了这么久的 C++ 终于说到类了,还是从内存出发来讨论一下 C++ 的类在内存中的存储方式(之前写过一篇内存对齐的文章,类同样在一定程度上遵循内存对齐原则,不过比结构体复杂一下) 如有侵权,请联系 ...
JuliaFEM中的数据格式——fields.jl
文章目录 Fields(场) 四种离散场四种连续场字典格式定义 new_field函数 Fields(场) 在JuliaFEM中,从项目一开始,我们脑海中就有一个清晰的概念:"一切都是场 ...
c2064 项不会计算为接受0个参数的函数_无网格法理论与Matlab程序设计（6）——传统径向基点插值（RPIM）形函数...
参考资料 G.R.Liu Y.T.GU著王建明周学军译 <无网格法理论及程序设计> 数值实现 Matlab 2019a 前情回顾形式主义的居士:无网格法理论与Matlab程序设计(1 ...
python里删除range里的数字_python中range函数与列表中删除元素
一.range函数使用 range(1,5) 代表从1到4(不包含5),结果为:1,2,3,4 ,默认步长为1 range(1,5,2) 结果为:1, 3 (同样不包含5) ,步长为2 ...
matlab绘制心形函数
matlab 7.0 绘制二维.三维心形函数又到周六,下周就要迎来春节小长假了,想想都有些激动.在外漂了一整年,总于可以回家和父母团聚了,还有吃好吃的...,哎呀~想想都流口水呢.不过先不要激动,假 ...
OpenGL绘制心形函数
OpenGL绘制心形函数用的最后一个 r =(float) (r_beishu*(Math.sin(Math.PI*i/180f)*Math.sqrt(Math.abs(Math.cos(Math. ...
matlab实例——动态心形函数及其涉及的知识点
心形特别感谢心形函数1 心形函数2 知识点 figure(...)函数的一些用法第一种用法最简单第二种用法第三种用法最后一种用法 set函数 num2str(n)解释 ezplot一维绘图 ...
笛卡尔心形函数表达式_几何画板制作笛卡尔心形函数的详细操作方法
朋友们或许不知道几何画板怎样制作笛卡尔心形函数的详细操作,那么今天绿软吧就讲解几何画板制作笛卡尔心形函数的详细操作方法哦,希望能够帮助到大家呢. 1.新建参数.右键绘图区空白处,"新建参数& ...

有限元形函数及JuliaFEM中的实现方式