有限元基础(一) Jacobian 矩阵和高斯积分

           好久没有更新博客了~前一段又有点对人生迷茫，玩了一段三国志13，fm，很是懈怠，五一前本来想在把读书时的代码找出来，转换成python代码的任务迟迟没有完成。今天，迎来了三年股市最大跌幅，三月清仓，月定投华能国际，竟然在四月底跑赢了指数，定投的指数基金也在上月成功减仓！准备月定投中国中铁和中国石油，现在这里mark一下！

&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

          回归正题，今天要用Python实现了任意四边形单元的面积计算，需要理解的几个基本概念：形函数，等参单元、jacobian矩阵和高斯积分。

         形函数：这个概念的提出前提是单元内所感兴趣的函数是连续的，因此可以建立一组插值函数表达式，方便获得单元内某个坐标所对应的函数值。

 形函数和等餐单元如图所示：

         

          等参单元：它是作为唯一一个标准的四边形单元存在的，其他形状的四边形都是由它乘以Jacobian矩阵变换而来的，打个不恰当的比喻，等参单元就好比是父亲，他所有的孩子都是父亲基因依照某种变换法则得到的。

         Jacobian矩阵：它表征的就一个任意形状的四边形单元和等参单元之间的转换关系（也就是上面所说的变换法则），包括各个方向的拉伸和旋转，其行列式，则是该单元与等参单元面积的比值。 

         高斯积分：前提要求，四边形内的函数必须是一个连续函数，而采用高斯积分的目的，是将连续积分转换成离散积分，便于计算机编程实现，因此高斯积分点越多，其精度越高，当然，精度高会带来计算时间的增加。如果一个单元内要积分的函数非线性程度不是很强的话，可以考虑采用较少的高斯积分点。进一步可以得出，如果单元尺寸足够小，其内部的非线性将会进一步降低，因此可以考虑选择较少的高斯积分点进行后续的积分计算。参见下图。

# -*- coding: utf-8 -*-
"""
Created on Sat Apr 27 07:05:44 2019

@author: Yujin.Wang

"""
import numpy as np
import pandas as pd


def funcQuadShapeDiff(eps,eta):
'''形函数导数'''
	NDiff = []

	NDiff = [[-0.25*(1-eta),-0.25*(1+eta),0.25*(1+eta),0.25*(1-eta)],
		   [-0.25*(1-eps),0.25*(1-eps),0.25*(1-eps),-0.25*(1-eps)]]
	return np.matrix(NDiff)


def funcQuadJacobian(QuadNode):
''''Jacobian矩阵计算和高斯积分'''
	a = pow(1./3,1./2)
	GaussPoint = [[-a,a],[a,a],[a,-a],[-a,-a]] #等参单元的二阶高斯积分点坐标
	GaussSum = 0 
	for i in range(4):
		J = np.matrix([[0.,0.],[0.,0.]])
		NDiff = funcQuadShapeDiff(GaussPoint[i][0],GaussPoint[i][1])

		J[0,0] = np.matmul(NDiff[0],QuadNode[:,1])
		J[0,1] = np.matmul(NDiff[0],QuadNode[:,2])
		J[1,0] = np.matmul(NDiff[1],QuadNode[:,1])
		J[1,1] = np.matmul(NDiff[1],QuadNode[:,2])
		GaussSum +=np.linalg.det(J)
	return GaussSum



if __name__ == "__main__":
	QuadNode = [[1,1.,-0.],[2,1.,3.],[3,-1.,3.],[4,-1.,0.]]#节点编号，X坐标，Y坐标
	QuadNode = np.matrix(QuadNode)#转变成矩阵形式
	a = funcQuadJacobian(QuadNode)#计算单元面积
	print "Area:",a

	

对应的文献可参加I.M.Smith的有限元方法编程。

参考文献：

       I.M,Smith 著，张新春等 译. 有限单元方法编程（第五版）