!=============================================================
! 平面8节点等参元程序
! 清华大学土木工程系：陆新征
!============================================================

module Elem_Rect8 ! 八节点等参元
implicit none
integer (kind(1)),parameter ::ikind=(kind(1))
integer (kind(1)),parameter ::rkind=(kind(0.d0))

type :: typ_Kcol
real(rkind),pointer :: Row(:)
end type typ_Kcol

type :: typ_GValue !总体控制变量
integer(ikind) :: NNode, NElem, NLoad, NMat, NSupport
integer(ikind) :: NGlbDOF !整体自由度总数
integer(ikind) :: NGENS, NodeDOF,ElemNodeNo
integer(ikind) :: NInt
end type typ_GValue

type Typ_Node !定义节点类型
real(rkind) :: coord(2) !节点坐标
integer(ikind) :: GDOF(2) !整体自由度编码
real(rkind) :: DISP(2) !节点位移
real(rkind) :: dDISP(2) !节点位移增量
real(rkind) :: dForce(2) !节点不平衡力

end type typ_Node

!=============================================================================
Type typ_IntPoint !定义积分点参数
real(rkind) :: EPS(3) !应变
real(rkind) :: SIG(3) !应力
real(rkind) :: D(3,3) !本构矩阵
real(rkind) :: B(3,16) !几何矩阵
real(rkind) :: DETJ !雅克比行列式
end type Typ_IntPoint

type Typ_Rect8 !定义实体单元
integer(ikind) :: NodeNo(8) !节点编号
real(rkind) :: E !弹性模量
real(rkind) :: u !泊松比
real(rkind) :: t !单元厚度
real(rkind) :: EK(16,16) !单元刚度矩阵
type(typ_intpoint) :: IntP(9) !积分点
!........................
end type Typ_Rect8

type Typ_Load ! 定义荷载类型
integer(ikind) :: NodeNo !荷载作用节点编号
real(rkind) :: Value(2) !荷载大小
end type Typ_Load

type Typ_Support ! 定义支座类型
integer(ikind) :: NodeNo !支座节点编号
integer(ikind) :: DOF !支座约束自由度
real(rkind) :: Value !支座位移大小
end type Typ_Support

contains

subroutine Get_Elem_K(GValue,Elem,Node) ! 核心程序，得到单元的刚度矩阵
implicit none
type(typ_GValue) :: GValue
type(typ_Node) :: Node(:)
type(typ_Rect8) :: Elem(:)
real*8 :: InvJ(2,2)
real*8:: Coord(8,2),IntPCoord(9,2),IntPwt(9),dN(2,8)
integer :: I,ElemNo

call Get_IntP_Prop(IntPCoord,IntPWt) ! 计算积分点信息

do ElemNo=1,size(Elem)
Elem(ElemNo)%EK=0.0
do I=1, GValue%ElemNodeNo; !得到节点坐标
coord(I,:)=Node(Elem(ElemNo)%NodeNo(i))%coord;
end do
Elem(ElemNo)%EK=0d0
DO I=1,size(Elem(ElemNo)%IntP)
! 计算本构矩阵
call Get_D(Elem(ElemNo)%IntP(I)%D,Elem(ElemNo)%E,Elem(ElemNo)%u)
! 计算形函数对局部坐标导数
call Get_dN_dxi(dN,IntPCoord(i,1),IntPCoord(i,2))
InvJ=matmul(dN, Coord)
! 得到雅克比行列式
Elem(ElemNo)%IntP(I)%DETJ=InvJ(1,1)*InvJ(2,2)-InvJ(1,2)*InvJ(2,1)
InvJ=matinv(InvJ); ! 对雅克比行列式求逆
dN = matmul(InvJ,dN) ! 得到形函数对整体坐标的导数
call Get_B (Elem(ElemNo)%IntP(I)%B,dN) ! 得到几何矩阵
Elem(ElemNo)%EK=Elem(ElemNo)%EK+&
matmul(matmul(transpose(Elem(ElemNo)%IntP(I)%B),&
Elem(ElemNo)%IntP(I)%D),Elem(ElemNo)%IntP(I)%B)*&
Elem(ElemNo)%IntP(I)%DetJ*IntPWt(i)*Elem(ElemNo)%t

end do
end do
return
end subroutine

subroutine Get_D(D,e,v) ! 得到本构矩阵，平面应力
implicit none
real*8 :: e,v
real*8 :: D (:,:)
D=0.d0
D(1,1)=1D0; D(2,2)=1d0
D(1,2)=v; D(2,1)=v
D(3,1:2)=0d0; D(1:2,3)=0d0;
D(3,3)=(1.d0-v)/2.d0;
D=E/(1.d0-v*v)*D;
return
end subroutine Get_D

subroutine Get_IntP_Prop(IntPCoord,IntPWt) ! 得到高斯积分点的参数
implicit none
real*8 ::IntPCoord(:,:),IntPWt(:) ! 高斯积分点坐标，高斯积分点权重
real*8 :: z,x(3),y(3),w(3)
integer :: i,j
z=.2d0*sqrt(15.d0)
x=(/-1.d0,0.d0,1.d0/); y=(/1.d0,0.d0,-1.d0/)
w=(/5.d0/9.d0,8.d0/9.d0,5.d0/9.d0/)
do i=1,3
do j=1,3
IntPCoord((i-1)*3+j,1)=x(j)*z
IntPCoord((i-1)*3+j,2)=y(i)*z
IntPWt((i-1)*3+j)=w(i)*w(j)
end do
end do

return
end subroutine Get_IntP_Prop

subroutine Get_dN_dxi(dN_dxi,xi,eta) ! 得到形函数对局部坐标系的导数
implicit none
real*8 :: dN_dxi(:,:), xi,eta
real*8:: x1, x2, x3, x4, x5 ,x6 ,x7 ,x8
x1=.25*(1.-eta);
x2=.25*(1.+eta);
x3=.25*(1.-xi);
x4=.25*(1.+xi)

dN_dxi(1,1)=x1*(2.*xi+eta)
dN_dxi(1,2)=-8.*x1*x2
dN_dxi(1,3)=x2*(2.*xi-eta)
dN_dxi(1,4)=-4.*x2*xi
dN_dxi(1,5)=x2*(2.*xi+eta)
dN_dxi(1,6)=8.*x2*x1
dN_dxi(1,7)=x1*(2.*xi-eta)
dN_dxi(1,8)=-4.*x1*xi
dN_dxi(2,1)=x3*(xi+2.*eta)
dN_dxi(2,2)=-4.*x3*eta
dN_dxi(2,3)=x3*(2.*eta-xi)
dN_dxi(2,4)=8.*x3*x4
dN_dxi(2,5)=x4*(xi+2.*eta)
dN_dxi(2,6)=-4.*x4*eta
dN_dxi(2,7)=x4*(2.*eta-xi)
dN_dxi(2,8)=-8.*x3*x4
return
end subroutine Get_dN_dxi

subroutine Get_B(B,dN_dx) ! 得到几何矩阵
implicit none
real*8 :: B(:,:), dN_dx(:,:)
integer::k,l,m,n , ih,nod;
real*8 :: x,y,z
B=0.
do m=1,8
k=2*m; l=k-1; x=dN_dx(1,m); y=dN_dx(2,m)
B(1,l)=x; B(3,k)=x; B(2,k)=y; B(3,l)=y
end do
return
end subroutine Get_B

function matinv(A) result (B) ! 计算2x2逆矩阵
real(rkind) ,intent (in)::A(:,:)
!real(rkind) , allocatable::B(:,:)
real(rkind) , pointer::B(:,:)
real(rkind) :: x
integer :: N
N=size(A,dim=2)
allocate(B(N,N))
B(1,1)=A(2,2)
B(2,2)=A(1,1)
B(1,2)=-A(1,2)
B(2,1)=-A(2,1)
x=A(1,1)*A(2,2)-A(1,2)*A(2,1)
B=B/x
return
end function matinv

end module