等参单元与坐标变换（Julia实现）

矩形单元比三角形单元有更高的精度，但它不适应不规则形状。而任意直线四边形单元可以适应不规则形状。但问题是，如果有一边的边界方程为 $y=bx+c$，代入位移插值函数 $u={\alpha }_1+{\alpha }_{2}x+{\alpha }_{3}y+{\alpha }_{4}xy={\alpha }_1+{\alpha }_{2}x+{\alpha }_3(bx+c)+{\alpha }_{4}x(bx+c)={\beta }_1+{\beta }_{2}x+{\beta }_3{x}^2$，显然，在边界上位移不再是线性分布，因而不能由节点位移惟一确定，即：即使相邻两单元节点位移相同，也不能保证边界位移连续。这就限制了四边形单元的应用范围。因此，为了实际应用，寻找适当方式，将规则单元转化为曲线单元，是非常必要的。普遍采用的方法是 等参变换，即单元 几何形状和单元内的 位移函数采用相同的数目的结点参数和插值函数进行变换。采用等参变换后，关键问题是将整体坐标系里的单元刚度、荷载等单元特性变换为局部坐标系中计算。

等参单元的概念

如图所示的两套坐标系，一套是 Oxy ，用于实际单元，称为子单元；一套是Oξη ，它是标准后化的正方形单元，称为母单元。从母单元项子单元变换，实际上就是建立两个坐标系之间的变换关系，即:
$$A^{{}^{\prime }}(x,y)=f(A(\xi ,\eta ))$$

将整体坐标用插值形式表示，即：
$$x=\sum _{i=1}{N}_i(\xi ,\eta )x_i,y=\sum _{i=1}{N}_i(\xi ,\eta )y_i$$
$(x_i,y_i)$为整体坐标系下各顶点的坐标；其中，Ni称为插值函数（形函数)。通过这一变换，对母单元上每一点(ξ,η)可以在实际单元中找到一对应点(x，y)，这样就在两个单元间建立了一一对应关系。
在实际单元中建立位移函数的插值形式：
$$u=\sum _{i=1}{N}_i(\xi ,\eta )u_i,v=\sum _{i=1}{N}_i(\xi ,\eta )v_i\text{\hspace{1em}(1)}$$
因为坐标变换和位移变换都用同样个数的相应的节点值做参数，并且具有完全相同的形函数作为变换系数，所以称这种单元为等参数单元。

坐标变换

$$x=\sum _{i=n}{N}_i(\xi ,\eta )x_i(n=i,j,k,m...),y=\sum _{i=n}{N}_i(\xi ,\eta )y_i(n=i,j,k,m...)\text{\hspace{1em}(2)}$$

几何方程及雅可比矩阵

$$\left\{\begin{array}{c}\epsilon \end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\end{array}\right\}\text{\hspace{1em}(3)}$$
将（1）式带入（3）式可得：
$$\left\{\begin{array}{c}\epsilon \end{array}\right\}=\left[\begin{array}{cc}\frac{\partial }{\partial x}& 0\\ 0& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]\ast \left[\begin{array}{ccccc}{N}_1& N_2& N_3& N_4& ...\end{array}\right]\ast \left[\begin{array}{cc}{u}_1& v_1\\ u_2& v_2\\ u_3& v_3\\ u_4& v_4\\ ...& ...\end{array}\right]$$
求上式$\frac{\partial N_i}{\partial x}，\frac{\partial N_i}{\partial y}$,由复合函数求导可得：
$$\begin{gathered} \frac{\partial N_i}{\partial \xi }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \xi }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \xi } \\ \frac{\partial N_i}{\partial \eta }=\frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \eta }+\frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \eta }\text{\hspace{1em}(4)} \end{gathered}$$
即:
$$\left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\end{array}\right\}=[J]\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}$$
其中矩阵[J]为雅可比矩阵
$$\left[\begin{array}{c}J\end{array}\right]=\left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]=\left[\begin{array}{cc}\sum \frac{\partial N_i}{\partial \xi }{x}_i& \sum \frac{\partial N_i}{\partial \xi }{y}_i\\ \sum \frac{\partial N_i}{\partial \eta }{x}_i& \sum \frac{\partial N_i}{\partial \eta }{y}_i\end{array}\right]=\left\{\begin{array}{cccc}\frac{\partial N_1}{\partial \xi }& \frac{\partial N_2}{\partial \xi }& \frac{\partial N_3}{\partial \xi }& \frac{\partial N_4}{\partial \xi }\\ \frac{\partial N_1}{\partial \eta }& \frac{\partial N_2}{\partial \eta }& \frac{\partial N_3}{\partial \eta }& \frac{\partial N_4}{\partial \eta }\end{array}\right\}\ast \left\{\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\\ x_4& y_4\end{array}\right\}$$
因此：
$$\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}=[J{]}^{-1}\ast \left\{\begin{array}{c}\frac{\partial N_i}{\partial \xi }\\ \frac{\partial N_i}{\partial \eta }\end{array}\right\}$$

雅可比矩阵实现（JuliaFEM)

在JuliaFEM的FEMBasis包中math.jl中实现了雅可比矩阵的计算，函数需要三个参数分别为：

1. 单元类型
2. 整体坐标系下节点的坐标值，即上述$(x_i,y_i)$
3. 待求点的局部坐标系下的坐标值

首先需要计算$\frac{\partial N_i}{\partial \xi }，\frac{\partial N_i}{\partial \eta }$，需要用到eval_dbasis(B, xi)函数，B为单元类型，xi为待求点的局部坐标系下的坐标值。详见有限元形函数及JuliaFEM中的实现方式。

'''
jacobian(B, X, xi)
Given basis B, calculate jacobian at xi.
# Example
#单元类型
B = Quad4()
#整体坐标系下的节点坐标
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
#(0.0, 0.0)为待求点的局部坐标系下的坐标
jacobian(B, X, (0.0, 0.0))
# output
2×2 Array{Float64,2}:
 0.5  0.0
 0.0  0.5 
'''
function jacobian(B, X, xi)
	#dB求dN/dξ,dN/dη
    dB = eval_dbasis(B, xi)
    dim1, nbasis = size(B)
    dim2 = length(first(X))
    J = zeros(dim1, dim2)
    for i=1:dim1
        for j=1:dim2
            for k=1:nbasis
            	  #分项相乘后累加
                J[i,j] += dB[i,k]*X[k][j]
            end
        end
    end
    return J
end    

grad函数是对$\left\{\begin{array}{c}\frac{\partial N_i}{\partial x}\\ \frac{\partial N_i}{\partial y}\end{array}\right\}$的求解，

'''
grad(B, X, xi)
Given basis B, calculate gradient dN/dX at xi.
# Example
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
grad(B, X, (0.0, 0.0))
# output
2×4 Array{Float64,2}:
 -0.5   0.5  0.5  -0.5
 -0.5  -0.5  0.5   0.5 
'''
function grad(B, X, xi)
    J = jacobian(B, X, xi)
    dB = eval_dbasis(B, xi)
    G = inv(J) * dB
    return G
end

由（4）式可知，$$\left\{\begin{array}{c}\epsilon \end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u}{\partial x}\\ \frac{\partial v}{\partial y}\\ \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\end{array}\right\}$$ u，v为各节点的整体坐标系的位移，在math.jl中另一个grad(B, T, X, xi)函数即用于求解上式的应变，此处的T为单元各节点的整体坐标系下的位移向量。求得的结果为：$\left[\begin{array}{cc}\sum \frac{\partial N_i}{\partial x}{u}_i& \sum \frac{\partial N_i}{\partial x}{v}_i\\ \sum \frac{\partial N_i}{\partial y}{u}_i& \sum \frac{\partial N_i}{\partial y}{v}_i\end{array}\right]=\left[\begin{array}{cc}\frac{\partial u}{\partial x}& \frac{\partial v}{\partial x}\\ \frac{\partial u}{\partial y}& \frac{\partial v}{\partial y}\end{array}\right]$

grad(B, T, X, xi)
Calculate gradient of `T` with respect to `X` in point `xi` using basis `B`.
# Example
B = Quad4()
X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
grad(B, u, X, (0.0, 0.0))
# output
2×2 Array{Float64,2}:
1.5  0.5
1.0  2.0
"""
function grad(B, T, X, xi)
#B的行函数的梯度
G = grad(B, X, xi)
nbasis = length(B)
#与[T]的内积
dTdX = sum(kron(G[:,i], T[i]') for i=1:nbasis)
return dTdX'
end