be2die.for

两维弹性边界元程序

	ShpTTipY1N0(1,1,1,1)=-2.d0
	ShpTTipY1N0(2,1,1,1)=0.d0
	ShpTTipY1N0(3,1,1,1)=0.d0
	ShpTTipY1N0(0,2,1,1)=0.d0
	ShpTTipY1N0(1,2,1,1)=0.d0
	ShpTTipY1N0(2,2,1,1)=0.d0
	ShpTTipY1N0(3,2,1,1)=0.d0
	ShpTTipY1N0(0,3,1,1)=0.d0
	ShpTTipY1N0(1,3,1,1)=0.d0
	ShpTTipY1N0(2,3,1,1)=0.d0
	ShpTTipY1N0(3,3,1,1)=0.d0
!linear
!ordinary element
	T1=1.d0/4.d0
	T2=3.d0/4.d0
	ShpTTipY1N0(0,1,0,2)=T2/(T2-T1)
	ShpTTipY1N0(1,1,0,2)=1.d0/(T1-T2)
	ShpTTipY1N0(2,1,0,2)=0.d0
	ShpTTipY1N0(3,1,0,2)=0.d0
	ShpTTipY1N0(0,2,0,2)=T1/(T1-T2)
	ShpTTipY1N0(1,2,0,2)=1.d0/(T2-T1)
	ShpTTipY1N0(2,2,0,2)=0.d0
	ShpTTipY1N0(3,2,0,2)=0.d0
	ShpTTipY1N0(0,3,0,2)=0.d0
	ShpTTipY1N0(1,3,0,2)=0.d0
	ShpTTipY1N0(2,3,0,2)=0.d0
	ShpTTipY1N0(3,3,0,2)=0.d0
!crack tip element
	T1=1.d0/4.d0
	T2=3.d0/4.d0
	T3=1.d0
	ShpTTipY1N0(0,1,1,2)=T2*T3/(T1-T2)/(T1-T3)
	ShpTTipY1N0(1,1,1,2)=-(T2+T3)/(T1-T2)/(T1-T3)
	ShpTTipY1N0(2,1,1,2)=1.d0/(T1-T2)/(T1-T3)
	ShpTTipY1N0(3,1,1,2)=0.d0
	ShpTTipY1N0(0,2,1,2)=T1*T3/(T2-T1)/(T2-T3)
	ShpTTipY1N0(1,2,1,2)=-(T1+T3)/(T2-T1)/(T2-T3)
	ShpTTipY1N0(2,2,1,2)=1.d0/(T2-T1)/(T2-T3)
	ShpTTipY1N0(3,2,1,2)=0.d0
	ShpTTipY1N0(0,3,1,2)=0.d0
	ShpTTipY1N0(1,3,1,2)=0.d0
	ShpTTipY1N0(2,3,1,2)=0.d0
	ShpTTipY1N0(3,3,1,2)=0.d0
!quadrature
!ordinary element
	T1=1.d0/6.d0
	T2=0.5d0
	T3=5.d0/6.d0
	ShpTTipY1N0(0,1,0,3)=T2*T3/(T1-T2)/(T1-T3)
	ShpTTipY1N0(1,1,0,3)=-(T2+T3)/(T1-T2)/(T1-T3)
	ShpTTipY1N0(2,1,0,3)=1.d0/(T1-T2)/(T1-T3)
	ShpTTipY1N0(3,1,0,3)=0.d0
	ShpTTipY1N0(0,2,0,3)=T1*T3/(T2-T1)/(T2-T3)
	ShpTTipY1N0(1,2,0,3)=-(T1+T3)/(T2-T1)/(T2-T3)
	ShpTTipY1N0(2,2,0,3)=1.d0/(T2-T1)/(T2-T3)
	ShpTTipY1N0(3,2,0,3)=0.d0
	ShpTTipY1N0(0,3,0,3)=T1*T2/(T1-T3)/(T2-T3)
	ShpTTipY1N0(1,3,0,3)=-(T1+T2)/(T1-T3)/(T2-T3)
	ShpTTipY1N0(2,3,0,3)=1.d0/(T1-T3)/(T2-T3)
	ShpTTipY1N0(3,3,0,3)=0.d0
!crack tip element
	T1=1.d0/6.d0
	T2=0.5d0
	T3=5.d0/6.d0
	T4=1.d0
	ConDown=(T1-T2)*(T1-T3)*(T1-T4)
	ShpTTipY1N0(0,1,1,3)=-T2*T3*T4/ConDown
	ShpTTipY1N0(1,1,1,3)=(T2*T3+T3*T4+T4*T2)/ConDown
	ShpTTipY1N0(2,1,1,3)=-(T2+T3+T4)/ConDown
	ShpTTipY1N0(3,1,1,3)=1.d0/ConDown
	ConDown=(T2-T1)*(T2-T3)*(T2-T4)
	ShpTTipY1N0(0,2,1,3)=-T1*T3*T4/ConDown
	ShpTTipY1N0(1,2,1,3)=(T1*T3+T3*T4+T4*T1)/ConDown
	ShpTTipY1N0(2,2,1,3)=-(T1+T3+T4)/ConDown
	ShpTTipY1N0(3,2,1,3)=1.d0/ConDown
	ConDown=(T3-T1)*(T3-T2)*(T3-T4)
	ShpTTipY1N0(0,3,1,3)=-T1*T2*T4/ConDown
	ShpTTipY1N0(1,3,1,3)=(T1*T2+T2*T4+T4*T1)/ConDown
	ShpTTipY1N0(2,3,1,3)=-(T1+T2+T4)/ConDown
	ShpTTipY1N0(3,3,1,3)=1.d0/ConDown
*=================================================================
	return
	end
*********************************************************************
*********************************************************************
*********************************************************************
	Subroutine BInflCoef(
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NtotCfK,NthCfk,CfUP)
	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)
	Dimension CfUP(NCmp,NCmp*2,0:NtotCfKCap)
*********************************************************************
	NthCfK=0
	CfUP=0
*==========================================================
	Do 60 NES=1,NtotEle
	If(NES/(NtotEle/10)*(NtotEle/10).eq.NES) then
	 write(*,1010) (NES-1.d0)/NtotEle
	endif
	Do 60 NNS=1,NodEle(NES)
	 Do 40 NEF=1,NtotEle
	 Do 40 NNF=1,NodEle(NEF)
	  If(NthBdEle(NES).ne.NthBdEle(NEF)) goto 40
				JDeriv=0
				If(JCrkN0C1Tip2Ele(NES).gt.0) JDeriv=1
	  Call CfIntInf(
	1	NCmp,NthBdEle(NES),XNod(1,NNS,NES),OEle(NES),
	3	XEle(1,1,NEF),OEle(NEF),ShpTNod(0,NNF,NEF),
	5	JDeriv,CfUP(1,1,NtotCfK+1))
	  NtotCfK=NtotCfK+1	  
	  NthCfK(NES,NNS,NEF,NNF,0)=NtotCfK
	  If(JDeriv.eq.1) NtotCfK=NtotCfK+1
	  If(JDeriv.eq.1) NthCfK(NES,NNS,NEF,NNF,1)=NtotCfK
40	 continue
60	continue
*=========================================================
1010	format(1x,f4.2,' -- portion of InflCoef')
********************************************************************
	return
	end
****************************************************************
****************************************************************
****************************************************************
	Subroutine CfIntInf(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,JDeriv,CfUP)
      Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XSP(2),XFE(2,2),ShpTFE(0:3)
	Dimension CfUP(NCmp,NCmp*2,0:JDeriv)
!LOCAL ARRAYS
	Allocatable Fint(:,:,:),Fint0(:,:,:)
	Allocatable TN(:),FN(:,:,:,:),N0_E(:),N1_E(:),N2_E(:),JActE(:)
	Allocatable VNormSP(:),VNormFE(:)
	NEleCap=10000
	Allocate(Fint(NCmp,NCmp*2,0:JDeriv),Fint0(NCmp,NCmp*2,0:JDeriv))
	Allocate(TN(0:NEleCap*2),FN(NCmp,NCmp*2,0:JDeriv,0:NEleCap*2))
	Allocate(N0_E(NEleCap),N1_E(NEleCap),N2_E(NEleCap),JActE(NEleCap))
	Allocate(VNormSP(NCmp),VNormFE(NCmp))
*=================================================
	PI=4.d0*DATAN(1.d0)
	EscSplit=1.d-6
*=================================================
	DStoF12=(XSP(1)-XFE(1,1))*(XSP(2)-XFE(2,2))-
	1	(XSP(1)-XFE(1,2))*(XSP(2)-XFE(2,1))
	If(DABS(DStoF12).gt.1.d-9) goto 5
	Call CfAnaInf(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,JDeriv,Fint)
	goto 9999
*=================================================
5	NSplitStt=1
!=================================================
	TStt=0.d0
	TEnd=1.d0
!=================INTEGRATION=================
!INITIATION STEP	 
	NEle=2**NSplitStt
! NODES
	Do 40 NN=0,NEle*2
	TN(NN)=TStt+NN*(TEnd-TStt)/(NEle*2)
40	Call CfGnPnt(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,TN(NN),
	1	JDeriv,FN(1,1,0,NN))
! Eleents
	Do 45 NE=1,NEle
	N0_E(NE)=NE*2-2
	N1_E(NE)=NE*2-1
45	N2_E(NE)=NE*2
	Do 50 I=1,NCmp
	Do 50 J=1,NCmp*2
	Do 50 K=0,JDeriv
	Fint(I,J,K)=0.0
	Do 50 NE=1,NEle
	N0=N0_E(NE)
	N1=N1_E(NE)
	N2=N2_E(NE)
	FETemp=(FN(I,J,K,N0)+4.d0*FN(I,J,K,N1)
	1	+FN(I,J,K,N2))*(TN(N2)-TN(N0))/6.d0
50	Fint(I,J,K)=Fint(I,J,K)+FETemp
	
	Do 60 NE=1,NEle
60	JActE(NE)=1
!
	NSplit=NSplitStt
!================== 
101	NSplit=NSplit+1

	Fint0=Fint

	NEleTemp=NEle
 	Do 170 NE=1,NEleTemp
	 If(JActE(NE).EQ.0) GOTO 170
	 N0=N0_E(NE)
	 N1=N1_E(NE)
	 N2=N2_E(NE)
! Two new nodes
	 If(NEle.ge.NEleCap) then
	  WRITE(*,*)
	  WRITE(*,*) 'CfInf: NEle.gt.Allowed'
	  goto 899
	 endIf
	 NEle=NEle+1
	 TN(NEle*2-1)=(TN(N0)+TN(N1))/2.d0
	 TN(NEle*2)=(TN(N1)+TN(N2))/2.d0
	 Do 140 NN=NEle*2-1,NEle*2
140	 Call CfGnPnt(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,TN(NN),
	1	JDeriv,FN(1,1,0,NN))
! A NEW Element
	 N0_E(NE)=N0
	 N1_E(NE)=NEle*2-1
	 N2_E(NE)=N1
	 N0_E(NEle)=N1
	 N1_E(NEle)=NEle*2
	 N2_E(NEle)=N2
! UPDATE Fint()
	 FintImprMax=0.d0
	 Do 155 I=1,NCmp
	 Do 155 J=1,NCmp*2
	 Do 155 K=0,JDeriv
	 FED1=(FN(I,J,K,N0)+4.d0*FN(I,J,K,NEle*2-1)
	1	+FN(I,J,K,N1))*(TN(N1)-TN(N0))/6.d0						   
	 FED2=(FN(I,J,K,N1)+4.d0*FN(I,J,K,NEle*2)
	1	+FN(I,J,K,N2))*(TN(N2)-TN(N1))/6.d0
	 FED0=(FN(I,J,K,N0)+4.d0*FN(I,J,K,N1)
	1	+FN(I,J,K,N2))*(TN(N2)-TN(N0))/6.d0
	 FintImpr=FED1+FED2-FED0
	 Fint(I,J,K)=Fint(I,J,K)+FintImpr
155	 If(DABS(FintImpr).GT.FintImprMax) FintImprMax=DABS(FintImpr)
!Split CONTROL
	 If(FintImprMax.LE.EscSplit/DSQRT(2.d0)**NSplit) then
	  JActE(NE)=0
	  JActE(NEle)=0
	 else
	  JActE(NE)=1
	  JActE(NEle)=1
	 endIf
170	continue
!=============================================================
      NActEle=0
	Do 195 NE=1,NEle
195	If(JActE(NE).NE.0) NActEle=NActEle+1
	If(NActEle.NE.0) GOTO 101
!==================================================
899	CfEr=0.d0	
	Do 900 I=1,NCmp
	Do 900 J=1,NCmp*2
	Do 900 K=0,JDeriv
	CfErTemp=DABS(Fint(I,J,K)-Fint0(I,J,K))
900	If(CfEr.lt.CfErTemp) CfEr=CfErTemp
*=============================================================
!transform to local coordinates on both source and field points
9999	VNormSP(1)=DCOS(OSP-Pi/2.d0)
	VNormSP(2)=DSIN(OSP-Pi/2.d0)
	If(NCmp.eq.3) VNormSP(3)=0.d0
	VNormFE(1)=DCOS(OFE-Pi/2.d0)
	VNormFE(2)=DSIN(OFE-Pi/2.d0)
	If(NCmp.eq.3) VNormFE(3)=0.d0
	Do 1100 K=0,JDeriv
	Call TransCftoLocal(
	1	NCmp,Fint(1,1,K),VNormSP,VNormFE,CfUP(1,1,K))
1100	Call TransCftoLocal(
	1	NCmp,Fint(1,NCmp+1,K),VNormSP,VNormFE,CfUP(1,NCmp+1,K))
*==========================================================
	DEAllocate(Fint,Fint0,TN,FN,N0_E,N1_E,N2_E,JActE)
	DeAllocate(VNormSP,VNormFE)
	RETURN
	END
****************************************************************
****************************************************************
	Subroutine CfGnPnt(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,T,JDeriv,FN)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XSP(2),XFE(2,2),ShpTFE(0:3),FN(NCmp,NCmp*2,0:JDeriv)
	Allocatable GnU(:,:),GnP(:,:),GnUD(:,:,:),GnPD(:,:,:)
	Allocatable VNormSP(:),VNormFE(:)
	Allocate(GnU(NCmp,NCmp),GnP(NCmp,NCmp))
	Allocate(GnUD(NCmp,NCmp,NCmp),GnPD(NCmp,NCmp,NCmp))
	Allocate(VNormSP(NCmp),VNormFE(NCmp))
*=================================================================
	PI=4.d0*DATAN(1.d0)
*============================================
	RX12=DSQRT((XFE(1,2)-XFE(1,1))**2+(XFE(2,2)-XFE(2,1))**2)
	Wt=RX12*(ShpTFE(0)+ShpTFE(1)*T+ShpTFE(2)*T*T+ShpTFE(3)*T*T*T)
*=================================================================
	X1SP=XSP(1)
	X2SP=XSP(2)
	X1FP=XFE(1,1)+T*(XFE(1,2)-XFE(1,1))
	X2FP=XFE(2,1)+T*(XFE(2,2)-XFE(2,1))
	VNormSP(1)=DCOS(OSP-Pi/2.d0)
	VNormSP(2)=DSIN(OSP-Pi/2.d0)
	If(NCmp.eq.3) VNormSP(3)=0.d0
	VNormFE(1)=DCOS(OFE-Pi/2.d0)
	VNormFE(2)=DSIN(OFE-Pi/2.d0)
	If(NCmp.eq.3) VNormFE(3)=0.d0
*=================================================================
	Call UserGn(
	1	NthBd,X1SP,X2SP,X1FP,X2FP,VNormFE,
     1	GnU,GnP,JDeriv,GnUD,GnPD)
*=================================================================
	Do 110 I=1,NCmp
	Do 110 J=1,NCmp
	FN(I,J,0)=GnU(I,J)
110	FN(I,J+NCmp,0)=GnP(I,J)
*=================================================================
	If(JDeriv.eq.0) goto 999
*=================================================================
	Do 210 I=1,NCmp
	Do 210 J=1,NCmp
	FN(I,J,1)=0.d0
	FN(I,J+NCmp,1)=0.d0
	Do 210 K=1,NCmp
	FN(I,J,1)=FN(I,J,1)+GnUD(I,K,J)*VNormSP(K)
210	FN(I,J+NCmp,1)=FN(I,J+NCmp,1)+GnPD(I,K,J)*VNormSP(K)
*=================================================================
999	Do 900 K=0,JDeriv
	Do 900 I=1,NCmp
	Do 900 J=1,NCmp*2
900	FN(I,J,K)=FN(I,J,K)*Wt
*================================================================
	DeAllocate(GnU,GnP,GnUD,GnPD)
	DeAllocate(VNormSP,VNormFE)
	return
	end
*******************************************************************
*******************************************************************
*******************************************************************
	Subroutine CfAnaInf(NCmp,NthBd,XSP,OSP,XFE,OFE,ShpTFE,JDeriv,Cf)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension XSP(2),XFE(2,2),ShpTFE(0:3),Cf(NCmp,NCmp*2,0:JDeriv)
	Allocatable ShpSFE(:)
	Allocatable VNormSP(:),VNormFE(:)
	Allocatable GnUP1(:,:),GnPP1(:,:),GnUDP1(:,:,:),GnPDP1(:,:,:)
	Allocatable GnUP2(:,:),GnPP2(:,:),GnUDP2(:,:,:),GnPDP2(:,:,:)
	Allocatable CPlnGnU(:,:),CLogGnU(:,:),CR_1GnP(:,:),
	1	CR_1GnUD(:,:,:),CR_2GnPD(:,:,:)
	Allocate(ShpSFE(0:3))
	Allocate(VNormSP(NCmp),VNormFE(NCmp))
	Allocate(GnUP1(NCmp,NCmp),GnPP1(NCmp,NCmp),
	1	GnUDP1(NCmp,NCmp,NCmp),GnPDP1(NCmp,NCmp,NCmp))
	Allocate(GnUP2(NCmp,NCmp),GnPP2(NCmp,NCmp),
	1	GnUDP2(NCmp,NCmp,NCmp),GnPDP2(NCmp,NCmp,NCmp))
	Allocate(CPlnGnU(NCmp,NCmp),CLogGnU(NCmp,NCmp),CR_1GnP(NCmp,NCmp),
	1	CR_1GnUD(NCmp,NCmp,NCmp),CR_2GnPD(NCmp,NCmp,NCmp))
*==============================================================
	PI=4.d0*DATAN(1.d0)
*==============================================================
	X1SP=XSP(1)
	X2SP=XSP(2)
*==============================================================
	VNormSP(1)=DCOS(OSP-Pi/2.d0)
	VNormSP(2)=DSIN(OSP-Pi/2.d0)
	If(NCmp.eq.3) VNormSP(3)=0.d0
	VNormFE(1)=DCOS(OFE-Pi/2.d0)
	VNormFE(2)=DSIN(OFE-Pi/2.d0)
	If(NCmp.eq.3) VNormFE(3)=0.d0
*==============================================================
	Call UserGn(
	1	NthBd,X1SP,X2SP,
     1	X1SP+(XFE(1,2)-XFE(1,1)),X2SP+(XFE(2,2)-XFE(2,1)),VNormFE,
     1	GnUP1,GnPP1,JDeriv,GnUDP1,GnPDP1)
	Call UserGn(
	1	NthBd,X1SP,X2SP,
     1	X1SP+(XFE(1,2)-XFE(1,1))*2,X2SP+(XFE(2,2)-XFE(2,1))*2,VNormFE,
     1	GnUP2,GnPP2,JDeriv,GnUDP2,GnPDP2)