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

形函数

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

位移模式

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

1.三角形单元的位移模式为1+u+v；
2. 三个顶点的坐标分别为$(x_i，y_i),(x_j，y_j),(x_m，y_m)$，对应局部坐标系下的坐标为(0，0),(1，0),(0，1)；
3. 单元内任意点的位移可通过下式表达：
$$\begin{gathered} u={\alpha }_1+{\alpha }_{2}x+{\alpha }_{3}y \\ v={\alpha }_4+{\alpha }_{5}x+{\alpha }_{6}y\text{\hspace{1em}(1)} \end{gathered}$$
4.将节点坐标带入（1）式，可得顶点处的位移为：
$$\begin{gathered} u_i={\alpha }_1+{\alpha }_2{x}_i+{\alpha }_3{y}_i \\ {u}_j={\alpha }_1+{\alpha }_2{x}_j+{\alpha }_3{y}_j \\ {u}_m={\alpha }_1+{\alpha }_2{x}_m+{\alpha }_3{y}_m \\ {v}_i={\alpha }_4+{\alpha }_5{x}_i+{\alpha }_6{y}_i \\ {v}_j={\alpha }_4+{\alpha }_5{x}_j+{\alpha }_6{y}_j \\ {v}_m={\alpha }_4+{\alpha }_5{x}_m+{\alpha }_6{y}_m \end{gathered}$$
即：
$$\left\{\begin{array}{c}{u}_i\\ u_j\\ u_m\end{array}\right\}=\left[\begin{array}{ccc}1& x_i& y_i\\ 1& x_j& y_j\\ 1& x_m& y_m\end{array}\right]\ast \left\{\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\end{array}\right\}\text{\hspace{1em}(2)}$$
令：
$$\left[\begin{array}{c}T\end{array}\right]=\left[\begin{array}{ccc}1& x_i& y_i\\ 1& x_j& y_j\\ 1& x_m& y_m\end{array}\right]$$

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:N
        X = coordinates[i]
        if D == 1
            data = (:u => X[1],)
        elseif D == 2
            data = (:u => X[1], :v => X[2])
        elseif D == 3
            data = (:u => X[1], :v => X[2], :w => X[3])
        end
        for (j, term) in enumerate(args)
            A[i,j] = subs(term, data)
        end
    end
    return 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 data
        p = subs(p, di)
    end
    return 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）进行推导可得：
$$\left\{\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\end{array}\right\}={\left[\begin{array}{c}T\end{array}\right]}^{-1}\ast \left\{\begin{array}{c}{u}_i\\ u_j\\ u_m\end{array}\right\}\text{\hspace{1em}(3)}$$
将（3）式带入（1）式可得：
$$u=\left[\begin{array}{ccc}1& u& v\end{array}\right]\ast {\left[\begin{array}{c}T\end{array}\right]}^{-1}\ast \left\{\begin{array}{c}{u}_i\\ u_j\\ u_m\end{array}\right\}\text{\hspace{1em}(4)}$$
定义形函数：$$\left[\begin{array}{c}N\end{array}\right]=\left[\begin{array}{ccc}1& u& v\end{array}\right]\ast {\left[\begin{array}{c}T\end{array}\right]}^{-1}$$
这里的坐标为局部坐标，由$$\left[\begin{array}{c}T\end{array}\right]\ast {\left[\begin{array}{c}T\end{array}\right]}^{-1}=\left[\begin{array}{c}E\end{array}\right]$$
那么 $\left[\begin{array}{c}T\end{array}\right]\ast [X1]=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]$这里的X1为[T]^{-1}的第一列，同理 $\left[\begin{array}{c}T\end{array}\right]\ast [X2]=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]这里的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:N
        fill!(b, 0.0)
        b[i] = 1.0
        #计算逆矩阵的单列
        solution = V \ b
        N = Expr(:call, :+)
        for (ai, bi) in zip(solution, args)
            isapprox(ai, 0.0) && continue
            #位移模式多项式如[1 u v]与逆矩阵的单列相乘并放入数组中
            push!(N.args, simplify( :( $ai * $bi ) ))
        end
        push!(basis, N)
    end
    return 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 basis
        partial_derivatives = []
        for j in 1:D
        #分别对u v w进行求导
            push!(partial_derivatives, simplify(differentiate(N, vars[j])))
        end
        push!(dbasis, partial_derivatives)
    end
    return dbasis
end