📄 be2die.for

📁 两维弹性边界元程序
💻 FOR
📖 第 1 页 / 共 5 页
字号:
12 3 4 5 下一页
*******************************************************************
*******************************************************************
*******************************************************************
	USE numerical_libraries
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Integer*2 IYear1,IMonth1,IDay1,IHour1,IMin1,ISec1,IMilSec1
	Integer*2 IYear2,IMonth2,IDay2,IHour2,IMin2,ISec2,IMilSec2
!geometry
	Allocatable XEle(:,:,:),NodEle(:),OEle(:)
	Allocatable NthBdEle(:),NthGrpEle(:)
	Allocatable JCrkN0C1Tip2Ele(:),NthSlvMstEle(:)
	Allocatable XNod(:,:,:),ShpTNod(:,:,:),DNormGapMstNod(:,:)
!coefficient
	Allocatable NthCfK(:,:,:,:,:),CfUP(:,:,:)
!solution
	Allocatable JBCU1P2(:,:,:)
	Allocatable UNod(:,:,:),PNod(:,:,:)
	Allocatable UHomoNod(:,:,:),PHomoNod(:,:,:)
****************************************************************
	Call User_Example
	Call User_Open
***************************************************************
	CALL GETDAT(IYear1,IMonth1,IDay1)
	CALL GETTIM(IHour1,IMin1,ISec1,IMilSec1)
****************************************************************
	Call GetNECap(NCmp,NEAddCap,NECap)
	Allocate(XEle(2,2,NECap),NodEle(NECap),OEle(NECap))
	Allocate(NthBdEle(NECap),NthGrpEle(NECap))
	Allocate(JCrkN0C1Tip2Ele(NECap),NthSlvMstEle(NECap))
	Allocate(XNod(2,3,NECap),ShpTNod(0:3,3,NECap))
	Allocate(DNormGapMstNod(3,NECap))
*==============================================================
	Call AGeom(
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,DNormGapMstNod)
*==============================================================
	NtotNod=0
	Do NE=1,NtotEle
	NtotNod=NtotNod+NodEle(NE)
	enddo
	write(*,*) 'NtotEle=',NtotEle,', NtotNod=',NtotNod
*****************************************************************
	NtotCfKCap=0
	Do 110 NES=1,NtotEle
	Do 110 NNS=1,NodEle(NES)
	Do 110 NEF=1,NtotEle
	Do 110 NNF=1,NodEle(NEF)
	If(JCrkN0C1Tip2Ele(NES).eq.0) NtotCfKCap=NtotCfKCap+1
110	If(JCrkN0C1Tip2Ele(NES).ne.0) NtotCfKCap=NtotCfKCap+2
	Do 120 NES=1,NEAddCap
	Do 120 NNS=1,3
	Do 120 NEF=1,NtotEle
	Do 120 NNF=1,NodEle(NEF)
120	NtotCfKCap=NtotCfKCap+2
	Do 130 NES=1,NtotEle
	Do 130 NNS=1,NodEle(NES)
	Do 130 NEF=1,NEAddCap
	Do 130 NNF=1,3
	If(JCrkN0C1Tip2Ele(NES).eq.0) NtotCfKCap=NtotCfKCap+1
130	If(JCrkN0C1Tip2Ele(NES).ne.0) NtotCfKCap=NtotCfKCap+2
	Do 140 NES=1,NEAddCap
	Do 140 NNS=1,3
	Do 140 NEF=1,NEAddCap
	Do 140 NNF=1,3
140	NtotCfKCap=NtotCfKCap+2
*==============================================================
	Allocate(NthCfK(NECap,3,NECap,3,0:1))
	Allocate(CfUP(NCmp,NCmp*2,0:NtotCfKCap))
	Call BInflCoef(
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	NtotCfKCap,NtotCfK,NthCfk,CfUP)
*****************************************************************
	Allocate(JBCU1P2(NCmp,3,NECap))
	Allocate(UNod(NCmp,3,NECap),PNod(NCmp,3,NECap))
	Allocate(UHomoNod(NCmp,3,NECap),PHomoNod(NCmp,3,NECap))
*===============================================================
	NLoad=0
101	NLoad=NLoad+1
	Call CBCond(NLoad,NLoadMax,
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod)
*===============================================================
	Call DSolver(NLoad,
	1	NCmp,NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,DNormGapMstNod,
	1	NtotCfKCap,NtotCfK,NthCfk,CfUP,
	1	JBCU1P2,UNod,PNod,UHomoNod,PHomoNod)
*===============================================================
	If(NLoad.lt.NLoadMax) goto 101
*************************************************************
	Call GETDAT(IYear2,IMonth2,IDay2)
	Call GETTIM(IHour2,IMin2,ISec2,IMilSec2)
	Seconds21=(IDay2-IDay1)*24*3600+(IHour2-IHour1)*3600
     1	+(IMin2-IMin1)*60+(ISec2-ISec1)+(IMilSec2-IMilSec1)/100.0
	Min21=Int(Seconds21/60)
	MSecR21=Seconds21-Int(Seconds21/60)*60
	write(*,*) ' Total time to finish',Min21,MSecR21
******************************************************************
	DeAllocate(XEle,NodEle,OEle)
	DeAllocate(NthBdEle,NthGrpEle)
	DeAllocate(JCrkN0C1Tip2Ele,NthSlvMstEle)
	DeAllocate(XNod,ShpTNod)
	DeAllocate(DNormGapMstNod)
	DeAllocate(NthCfK,CfUP)
	DeAllocate(JBCU1P2)
	DeAllocate(UNod,PNod)
	DeAllocate(UHomoNod,PHomoNod)
*****************************************************************
	stop
	end
*********************************************************************
*********************************************************************
*********************************************************************
	Subroutine AGeom(
	1	NECap,NtotEle,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     1	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod,DNormGapMstNod)
	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 DNormGapMstNod(3,NECap)
**********************************************************************
	PI=4.d0*DATAN(1.d0)
**********************************************************************
	NB=1
	Call UserGeom(NtotBndSegm,NEAddCap,
	1	NB,JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
*===========================================================
	NE=0
!=============================================================
!straight segments of boundary and crack
	Do 100 NB=1,NtotBndSegm
	 Call UserGeom(NtotBndSegm,NEAddCap,
	1	NB,JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 If(NthBd.eq.Mst2NthBd) then
	  Call SegElemMst(NE,NthBd,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
     1	Gradient,NthGrp,
	2	NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     2	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
	 else
	  Call SegElemMstSlv(NE,NthBd,Mst2NthBd,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
     1	Gradient,NthGrp,
	2	NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     2	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
	 endif
100	continue
!===========================================================
	NtotEle=NE
*=============================================================
!Gap-distance of master-slave pair nodes
	Do 310 NE=1,NtotEle
	Do 310 NN=1,NodEle(NE)
	DNormGapMstNod(NN,NE)=0.d0
	If(NthSlvMstEle(NE).gt.0) then
	 NESlv=NthSlvMstEle(NE)
	 NNSlv=NodEle(NE)+1-NN
	 DNormGapMstNod(NN,NE)= !0.d0
	1	(XNod(1,NNSlv,NESlv)-XNod(1,NN,NE))*DSIN(OEle(NE))-
	1	(XNod(2,NNSlv,NESlv)-XNod(2,NN,NE))*DCOS(OEle(NE))
	endif
310	continue
*********************************************************************
	Open(51,file='ZAEleNods.txt',status='unknown')
	Do 410 NE=1,NtotEle
	write(51,*)
	write(51,9110) NE,XEle(1,1,NE),XEle(2,1,NE)
	write(51,9110) NE,XEle(1,2,NE),XEle(2,2,NE)
	Do 410 NN=1,NodEle(NE)	
	write(51,*)
410	write(51,9110) NN,XNod(1,NN,NE),XNod(2,NN,NE)
	close(51)
9110	format(1x,I4,2(1x,F10.6))
*********************************************************************
	return
	end
*********************************************************************
*********************************************************************
*********************************************************************
	Subroutine GetNECap(NCmp,NEAddCap,NECap)
*===========================================================
	NCmp=2
*===========================================================
	NB=1
	Call UserGeom(NtotBndSegm,NEAddCap,
	1	NB,JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
*===========================================================
	NE=0
!=============================================================
!straight segments of boundary and crack
	Do 100 NB=1,NtotBndSegm
	 Call UserGeom(NtotBndSegm,NEAddCap,
	1	NB,JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 If(NthBd.eq.Mst2NthBd) then
	  NE=NE+NDiv
	 else
	  NE=NE+NDiv*2
	 endif
100	continue
!===========================================================
	NECap=NE+NEAddCap
*===========================================================
	return
	end
*********************************************************************
*********************************************************************
*********************************************************************
	Subroutine SegElemMst(NE,NthBd,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
     1	Gradient,NthGrp,
	2	NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     2	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
	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)
!local
	Allocatable ShpTTipY1N0(:,:,:,:)
	Allocate(ShpTTipY1N0(0:3,3,0:1,3))
*=============================================================
	PI=4.d0*DATAN(1.d0)
*=========================================================
	Call ShpTTemplate(ShpTTipY1N0)
*=============================================================
	OAll=DATAN2(X2End-X2Stt,X1End-X1Stt)
	If(OAll.lt.0.d0) OAll=OAll+Pi*2.d0
!============================================================
	Do 100 ND=1,NDiv
	 NE=NE+1
	 If(DABS(Gradient-1.d0).gt.1.d-6) then
	  T1=(1.d0-Gradient**(ND-1))/(1.d0-Gradient**NDiv)
	  T2=(1.d0-Gradient**ND)/(1.d0-Gradient**NDiv)
	 else
	  T1=(ND-1.d0)/NDiv
	  T2=ND*1.d0/NDiv
	 endif
	 XEle(1,1,NE)=X1Stt+T1*(X1End-X1Stt)
	 XEle(2,1,NE)=X2Stt+T1*(X2End-X2Stt)
	 XEle(1,2,NE)=X1Stt+T2*(X1End-X1Stt)
	 XEle(2,2,NE)=X2Stt+T2*(X2End-X2Stt)
	 NodEle(NE)=NodE
	 OEle(NE)=OAll
	 NthBdEle(NE)=NthBd
	 NthGrpEle(NE)=NthGrp
	 JCrkN0C1Tip2Ele(NE)=0
	 If(JBnd0Crk1Tip2.gt.0) JCrkN0C1Tip2Ele(NE)=1
	 If(JBnd0Crk1Tip2.eq.2.and.ND.eq.NDiv) JCrkN0C1Tip2Ele(NE)=2
	 NthSlvMstEle(NE)=0
	 Do 30 NN=1,NodE
	 T=(NN-0.5d0)/NodE
	 XNod(1,NN,NE)=XEle(1,1,NE)+T*(XEle(1,2,NE)-XEle(1,1,NE))
	 XNod(2,NN,NE)=XEle(2,1,NE)+T*(XEle(2,2,NE)-XEle(2,1,NE))
	 JTip=0
	 If(JCrkN0C1Tip2Ele(NE).eq.2) JTip=1
	 Do 30 L=0,3
30	 ShpTNod(L,NN,NE)=ShpTTipY1N0(L,NN,JTip,NodE)
100	continue
*=========================================================
	DeAllocate(ShpTTipY1N0)
	return
	end
*********************************************************************
	Subroutine SegElemMstSlv(NE,NthBd,Mst2NthBd,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE,
     1	Gradient,NthGrp,
	2	NECap,XEle,NodEle,OEle,NthBdEle,NthGrpEle,
     2	JCrkN0C1Tip2Ele,NthSlvMstEle,XNod,ShpTNod)
	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)
!local
	Allocatable ShpTTipY1N0(:,:,:,:)
	Allocate(ShpTTipY1N0(0:3,3,0:1,3))
*=============================================================
	PI=4.d0*DATAN(1.d0)
*=========================================================
	Call ShpTTemplate(ShpTTipY1N0)
*=============================================================
	OAll=DATAN2(X2End-X2Stt,X1End-X1Stt)
	If(OAll.lt.0.d0) OAll=OAll+Pi*2.d0
!============================================================
	Do 100 ND=1,NDiv
	 NE=NE+1
	 If(DABS(Gradient-1.d0).gt.1.d-6) then
	  T1=(1.d0-Gradient**(ND-1))/(1.d0-Gradient**NDiv)
	  T2=(1.d0-Gradient**ND)/(1.d0-Gradient**NDiv)
	 else
	  T1=(ND-1.d0)/NDiv
	  T2=ND*1.d0/NDiv
	 endif
	 XEle(1,1,NE)=X1Stt+T1*(X1End-X1Stt)
	 XEle(2,1,NE)=X2Stt+T1*(X2End-X2Stt)
	 XEle(1,2,NE)=X1Stt+T2*(X1End-X1Stt)
	 XEle(2,2,NE)=X2Stt+T2*(X2End-X2Stt)
	 NodEle(NE)=NodE
	 OEle(NE)=OAll
	 NthBdEle(NE)=NthBd
	 NthGrpEle(NE)=NthGrp
	 JCrkN0C1Tip2Ele(NE)=0
	 Do 30 NN=1,NodE
	 T=(NN*2.d0-1.d0)/6.d0
	 XNod(1,NN,NE)=XEle(1,1,NE)+T*(XEle(1,2,NE)-XEle(1,1,NE))
	 XNod(2,NN,NE)=XEle(2,1,NE)+T*(XEle(2,2,NE)-XEle(2,1,NE))
	 JTip=0
	 Do 30 L=0,3
30	 ShpTNod(L,NN,NE)=ShpTTipY1N0(L,NN,JTip,NodE)
	 NthSlvMstEle(NE)=NE+1
!slave
	 NE=NE+1
	 XEle(1,1,NE)=XEle(1,2,NE-1)
	 XEle(2,1,NE)=XEle(2,2,NE-1)
	 XEle(1,2,NE)=XEle(1,1,NE-1)
	 XEle(2,2,NE)=XEle(2,1,NE-1)
	 OEle(NE)=OAll+Pi
	 NodEle(NE)=NodE
	 NthBdEle(NE)=Mst2NthBd
	 NthGrpEle(NE)=-NthGrp
	 JCrkN0C1Tip2Ele(NE)=0
	 Do 60 NN=1,NodE
	 T=(NN-0.5d0)/NodE
	 XNod(1,NN,NE)=XEle(1,1,NE)+T*(XEle(1,2,NE)-XEle(1,1,NE))
	 XNod(2,NN,NE)=XEle(2,1,NE)+T*(XEle(2,2,NE)-XEle(2,1,NE))
	 JTip=0
	 Do 60 L=0,3
60	 ShpTNod(L,NN,NE)=ShpTTipY1N0(L,NN,JTip,NodE)
	 NthSlvMstEle(NE)=-(NE-1)
100	continue
*=========================================================
	DeAllocate(ShpTTipY1N0)
	return
	end
**************************************************************
	Subroutine ShpTTemplate(ShpTTipY1N0)
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Dimension ShpTTipY1N0(0:3,3,0:1,3)
*=============================================================
!constant
!ordinary element
	ShpTTipY1N0(0,1,0,1)=1.d0
	ShpTTipY1N0(1,1,0,1)=0.d0
	ShpTTipY1N0(2,1,0,1)=0.d0
	ShpTTipY1N0(3,1,0,1)=0.d0
	ShpTTipY1N0(0,2,0,1)=0.d0
	ShpTTipY1N0(1,2,0,1)=0.d0
	ShpTTipY1N0(2,2,0,1)=0.d0
	ShpTTipY1N0(3,2,0,1)=0.d0
	ShpTTipY1N0(0,3,0,1)=0.d0
	ShpTTipY1N0(1,3,0,1)=0.d0
	ShpTTipY1N0(2,3,0,1)=0.d0
	ShpTTipY1N0(3,3,0,1)=0.d0
!crack tip element
	ShpTTipY1N0(0,1,1,1)=2.d0
12 3 4 5 下一页
💿 文件大小 35 K
👤 上传用户 5201314
📂 所属分类数学计算
🏷️ 相关标签

#弹性 #程序
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -