userk_postprcs.for

两维弹性边界元程序

!For User to post-process 
! whatever is wanted, based on solution at each loading step
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!NLoad: Nth loading step
!NECap: maximum number of elements allowed
!NtotEle: total number of elements
!XEle(x1/2,ends,NthE): coordinates (first dimension) of two ends 
!	(second dimension) of Nth Element (third dimension).
!NodEle(NthE): number of nodes in Nth element, which determines
!	order of interpolation of element.
!OEle(NthE): orientation angle of Nth Element.
!NthBdEle(NthE): Nth body to which Nth Element belongs.
!NthGrpEle(NthE): Nth group to which Nth Element belongs.
!JCrkN0C1Tip2Ele(NthEle): 0 if non-crack, 1 if crack (non-tip), or
!	2 if crack with tip at second end, for Nth Element.
!NthSlvMstEle(NthEle): Nth element which is slave to NthEle.
!XNod(x1/2,NthN,NthE): coordinates of Nth Node in Nth Element.
!	NthN must not exceed NodEle(NthE).
!ShpNod(0:3,NthN,NthE): shape function in Tylor expansion 
!	for Nth Node in Nth Element.
!NthCfK,CfUP: recording influence coefficient.
!JBCU1P2(NthC,NthN,NthE): index of boundary condition in Nth Component
!	of Nth Node of Nth Element. 
!	If =1, displacement U is given;
!	If =2, traction P is given.
!UNod(NthC,NthN,NthE), PNod(): displacement and traction in Nth 
!	Component of Nth Node of Nth Element.
!UHomoNod(NthC,NthN,NthE), PHomoNod(): homogeneous field due to far-field
!	loading, initial strain field (eigenstrain), etc.
!UCrkPair1(NthC,NthN,NthE),UCrkPair2(): displacements at face 1 and 2 
!	of a crack element.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!SUBROUTINE QuadraTangentSs(
!	1	ES0,PR0,EigenSn0,		!elastic constants and eigenstrain
!	2	XKnot1,UKnot1,PKnot1,	!coordinates, displacement and traction at point 1
!	2	XKnot2,UKnot2,PKnot2,	!----at point 2
!	2	XKnot3,UKnot3,PKnot3,	!----at point 3
!	3	J123,Ssij,Snij,SnEng)	!J123=1,2,3: point 1,2 or 3 as calculation location
! To compute stress, strain and strain energy at a boundary point 
! given displacements and tractions at THREE points
! For QUADRATIC elements, it is convenient to carry out this task within each element.
!SUBROUTINE LinearTangentSs(
!	1	ES0,PR0,EigenSn0,
!	2	XKnot1,UKnot1,PKnot1,
!	2	XKnot2,UKnot2,PKnot2,
!	3	J12,Ssij,Snij,SnEng)
! To compute stress, strain and strain energy at a boundary point 
! given displacements and tractions at TWO points
! For LINEAR elements, it is convenient to carry out this task within each element.
!HOWEVER, for CONSTANT elements, neighbor nodes need to be found. Then, either
! quadra_ or linear_subroutine may be used.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Subroutine for DISPLACEMENT and STRESS at an internal point
! Subroutine InterDisplStress(
!	1	X1Pnt,X2Pnt,NthBd,UPnt,JAcqSsY1N0,SsPnt,
!	1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
!     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
!	1	UNod,PNod,UHomoNod,PHomoNod)
! X1Pnt,X2Pnt: coordinates of the point
! NthBd: Nth body to which this point belongs
! UPnt(:): displacement in the global coordinates
! JAcqSsY1N0: =1 if to acquire stress; =0 if not
! SsPnt(:,:): stress in the global coordinates
*********************************************************
	Subroutine UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
     	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Common/Example/NthExample
	Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
	Dimension NthBdEle(NECap),NthGrpEle(NECap)
	Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
	Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
	Dimension NthCfK(NECap,3,NECap,3,0:1),CfUP(2,4,0:NtotCfKCap)
	Dimension JBCU1P2(2,3,NECap),UNod(2,3,NECap),PNod(2,3,NECap)
	Dimension UHomoNod(2,3,NECap),PHomoNod(2,3,NECap)
	Dimension UCrkPair1(2,3,NECap),UCrkPair2(2,3,NECap)
*================================================================
	If(NthExample.eq.1) then
	 Call Ex1_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
	 return
	elseif(NthExample.eq.2) then
	 return !no post-processing
	elseif(NthExample.eq.3) then
	 Call Ex3_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
	 return
	elseif(NthExample.eq.4) then
	 Call Ex4_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
	 return
	endif
*================================================================
*================================================================
	return
	end
*********************************************************
	Subroutine Ex1_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
     	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
	Dimension NthBdEle(NECap),NthGrpEle(NECap)
	Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
	Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
	Dimension NthCfK(NECap,3,NECap,3,0:1),CfUP(2,4,0:NtotCfKCap)
	Dimension JBCU1P2(2,3,NECap),UNod(2,3,NECap),PNod(2,3,NECap)
	Dimension UHomoNod(2,3,NECap),PHomoNod(2,3,NECap)
	Dimension UCrkPair1(2,3,NECap),UCrkPair2(2,3,NECap)
!local
	Allocatable UPnt(:),SsPnt(:,:)
	Allocate(UPnt(2),SsPnt(2,2))
*================================================================
!Acquire U&Ss along (x1,x2=0.9)
	Do NP=1,19
	X1Pnt=NP*1.0/20
	X2Pnt=0.9
	NthBd=1
	JAcqSsY1N0=1
	Call InterDisplStress(
	1	X1Pnt,X2Pnt,NthBd,UPnt,JAcqSsY1N0,SsPnt,
	1	NLoad,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
	1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	UNod,PNod,UHomoNod,PHomoNod)
	write(12,7010) 
	1	X1Pnt,UPnt(1),UPnt(2),SsPnt(1,1),SsPnt(2,2),SsPnt(1,2)
	enddo
7010	format(10(1x,f18.9))
*================================================================
	return
	end
*********************************************************
	Subroutine Ex3_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
     	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
	Dimension NthBdEle(NECap),NthGrpEle(NECap)
	Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
	Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
	Dimension NthCfK(NECap,3,NECap,3,0:1),CfUP(2,4,0:NtotCfKCap)
	Dimension JBCU1P2(2,3,NECap),UNod(2,3,NECap),PNod(2,3,NECap)
	Dimension UHomoNod(2,3,NECap),PHomoNod(2,3,NECap)
	Dimension UCrkPair1(2,3,NECap),UCrkPair2(2,3,NECap)
!local
	Allocatable SsijNod(:,:,:,:),SnijNod(:,:,:,:)
	Allocatable SnEngNod(:,:)
	Allocate(SsijNod(2,2,3,NECap),SnijNod(2,2,3,NECap))
	Allocate(SnEngNod(3,NECap))
*================================================================
	Do 6600 NE=1,NtotEle
	NthBd=NthBdEle(NE)
	Call UserElasConst(ES0,PR0,NthBd)
	EigenSn0=UserIsoEigenSn(NLoad,NthBd)
	Do 6550 NN=1,NodEle(NE)
	J123=NN
6550	Call QuadraTangentSs(
	1	ES0,PR0,EigenSn0,
	2	XNod(1,1,NE),UNod(1,1,NE),PNod(1,1,NE),
	2	XNod(1,2,NE),UNod(1,2,NE),PNod(1,2,NE),
	2	XNod(1,3,NE),UNod(1,3,NE),PNod(1,3,NE),
	3	J123,SsijNod(1,1,NN,NE),SnijNod(1,1,NN,NE),SnEngNod(NN,NE))
6600	continue
!solutions along the top surface
	Do 6700 NE=1,NtotEle
	Do 6700 NN=1,NodEle(NE)
	If(NthGrpEle(NE).eq.15) then
	 write(11,6701) XNod(1,NN,NE),UNod(1,NN,NE),UNod(2,NN,NE),
	1	SsijNod(2,2,NN,NE),SnijNod(1,1,NN,NE),SnijNod(2,2,NN,NE),
	1	SnijNod(1,1,NN,NE)+SnijNod(2,2,NN,NE),SnEngNod(NN,NE)
	endif
6700	continue
6701	format(10(1x,f18.9))
*================================================================
	return
	end
*********************************************************
	Subroutine Ex4_UserPostPrcs(NLoad,
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod,
     1	UCrkPair1,UCrkPair2)
     	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XEle(2,2,NECap),NodEle(NECap),OEle(NECap)
	Dimension NthBdEle(NECap),NthGrpEle(NECap)
	Dimension JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap)
	Dimension XNod(2,3,NECap),ShpTNod(0:3,3,NECap)
	Dimension NthCfK(NECap,3,NECap,3,0:1),CfUP(2,4,0:NtotCfKCap)
	Dimension JBCU1P2(2,3,NECap),UNod(2,3,NECap),PNod(2,3,NECap)
	Dimension UHomoNod(2,3,NECap),PHomoNod(2,3,NECap)
	Dimension UCrkPair1(2,3,NECap),UCrkPair2(2,3,NECap)
!local
	Allocatable SsijNod(:,:,:,:),SnijNod(:,:,:,:)
	Allocatable SnEngNod(:,:)
	Allocate(SsijNod(2,2,3,NECap),SnijNod(2,2,3,NECap))
	Allocate(SnEngNod(3,NECap))
*==============================================================
	Do 6800 NE=1,NtotEle
	Do 6800 NN=1,NodEle(NE)
	If(NthSlvMstEle(NE).gt.0) then
	 NESlv=NthSlvMstEle(NE)
	 NNSlv=NodEle(NE)-NN+1
	 write(12,6801) XNod(1,NN,NE),
	1	UNod(1,NN,NE)+UNod(1,NNSlv,NESlv),
	1	UNod(2,NN,NE)+UNod(2,NNSlv,NESlv),
	1	PNod(1,NN,NE),PNod(2,NN,NE)
	endif
6800	continue
6801	format(10(1x,f18.9))
*================================================================
	return
	end