userc_geom.for

两维弹性边界元程序

!For User to prescribe
! GEOMETRY
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!It is assumed that the boundary consists of STRAIGHT SEGMENTS.
!For each straight segment, following parameters are used
! by MAIN PROGRAM to generate mesh:
!JBnd0Crk1WithTip2: =2 if a crack segment with a crack tip at the end;
!				  =1 if a crack segment without a crack tip;
!				  =0 if a regular external surface or an interface
!X1Stt!X2Stt: coordinates of start point||The start and end points should be 
!X1End!X2End: coordinates of end point  ||chosen such that the boundary outward 
!									  ||normal points to one's right-hand side 
!									  ||who walks from start to end. In the case
!									  ||of crack, this does not matter, unless
!									  ||it has a crack tip. In this case, the tip
!									  ||must be the end point.
!NDiv: number of elements in the segment
!NodEle: number of nodes (1, 2, or 3) in each element, i.e. order of interpretation
!Gradient: linear scaling factor of element size from start to end point
!		 =1 if uniform elements are defined through the segment
!NthGrp: Nth group to which all elements of the segment belong
!	Segments with crack tips should be assigned with different group numbers
!	It has been my own convention that noncrack groups are numbered from 1-100
!	while crack groups from 101-200.
!NthBd: Nth body to which all elements of the segment belong
!Mst2NthBd: =NthBd if no slave segment to the segment of NthBd;
!		  otherwise, matching nodes are assigned to the body of #Mst2NthBd
!	Therefore, that part of boundary for the body of #Mst2NthBd 
!	must not be repeated (see examples of multiple bodies).
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!NtotBndSegm: total number of straight segments
!NEAddCap: maximum additional number of elements for crack growth
!NthSegm: Nth segment
!Basically, through this subroutine, MAIN PROGRAM acquires
! NtotBndSegm, NEAddCap, and
! values of those parameters for each (NthSegm) segment.
! The format to implement this task can be flexible, depending on geometry.
***************************************************
	Subroutine UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd) !output
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
	Common/Example/NthExample
*================================================================
	If(NthExample.eq.1) then
	 Call Ex1_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 return
	elseif(NthExample.eq.2) then
	 Call Ex2_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 return
	elseif(NthExample.eq.3) then
	 Call Ex3_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 return
	elseif(NthExample.eq.4) then
	 Call Ex4_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd)
	 return
	endif
*====================================================
	NtotBndSegm=5
!===================================================
	NEAddCap=0
!============================================================
	If(NthSegm.eq.1) then
	 JBnd0Crk1Tip2=0
	 X1Stt=0.
	 X2Stt=0.
	 X1End=1.
	 X2End=0.
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=1
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.2) then
	elseif(NthSegm.eq.3) then
	elseif(NthSegm.eq.4) then
	elseif(NthSegm.eq.5) then
	endif
*============================================================
	return
	end
***************************************************
	Subroutine Ex1_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd) !output
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
*============================================================
	NtotBndSegm=8
!===================================================
	NEAddCap=0
!============================================================
	If(NthSegm.eq.1) then
	 JBnd0Crk1Tip2=0
	 X1Stt=0.
	 X2Stt=0.
	 X1End=1.
	 X2End=0.
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=1
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.2) then
	 JBnd0Crk1Tip2=0
	 X1Stt=1.
	 X2Stt=0.
	 X1End=1.
	 X2End=1.
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=2
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.3) then
	 JBnd0Crk1Tip2=0
	 X1Stt=1.
	 X2Stt=1.
	 X1End=0.
	 X2End=1.
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=3
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.4) then
	 JBnd0Crk1Tip2=0
	 X1Stt=0.
	 X2Stt=1.
	 X1End=0.
	 X2End=0.5
	 NDiv=2
	 NodE=3
	 Gradient=1.
	 NthGrp=4
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.5) then
	 JBnd0Crk1Tip2=0
	 X1Stt=0.
	 X2Stt=0.5
	 X1End=0.
	 X2End=0.
	 NDiv=2
	 NodE=3
	 Gradient=1.
	 NthGrp=4
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.6) then
	 JBnd0Crk1Tip2=2
	 X1Stt=0.
	 X2Stt=0.5
	 X1End=0.2
	 X2End=0.6
	 NDiv=3
	 NodE=3
	 Gradient=0.9
	 NthGrp=101
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.7) then
	 JBnd0Crk1Tip2=2
	 X1Stt=0.5
	 X2Stt=0.2
	 X1End=0.3
	 X2End=0.2
	 NDiv=3
	 NodE=3
	 Gradient=0.9
	 NthGrp=102
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.8) then
	 JBnd0Crk1Tip2=2
	 X1Stt=0.5
	 X2Stt=0.2
	 X1End=0.7
	 X2End=0.2
	 NDiv=3
	 NodE=3
	 Gradient=0.9
	 NthGrp=103
	 NthBd=1
	 Mst2NthBd=1
	endif
*============================================================
	return
	end
***************************************************
	Subroutine Ex2_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd) !output
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
*============================================================
	NtotBndSegm=2
!===================================================
	NEAddCap=80
!============================================================
	If(NthSegm.eq.1) then
	 JBnd0Crk1Tip2=2
	 X1Stt=0.
	 X2Stt=0.
	 X1End=SQRT(2.)/2.
	 X2End=SQRT(2.)/2.
	 NDiv=10
	 NodE=3
	 Gradient=0.9
	 NthGrp=101
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.2) then
	 JBnd0Crk1Tip2=2
	 X1Stt=0.
	 X2Stt=0.
	 X1End=-SQRT(2.)/2.
	 X2End=-SQRT(2.)/2.
	 NDiv=10
	 NodE=3
	 Gradient=0.9
	 NthGrp=102
	 NthBd=1
	 Mst2NthBd=1
	endif
*============================================================
	return
	end
***************************************************
	Subroutine Ex3_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd) !output
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
*============================================================
	NtotBndSegm=16
!===================================================
	NEAddCap=0
!============================================================
	HeightQW=0.1
	ThickSpacer=0.15
	ThickWetting=0.05
	X2Lvl1=ThickWetting
	X2Lvl2=0.d0
	X2Lvl3=-ThickSpacer
	X2Lvl4=X2Lvl3-ThickWetting
	X2Lvl5=-1.
	X1Width=1.
!----------------------------------------------------------
	If(NthSegm.eq.1) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-X1Width/2.
	 X2Stt=X2Lvl4
	 X1End=-X1Width/2.
	 X2End=X2Lvl5
	 NDiv=5
	 NodE=3
	 Gradient=1.25
	 NthGrp=1
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.2) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-X1Width/2.
	 X2Stt=X2Lvl5
	 X1End=X1Width/2.
	 X2End=X2Lvl5
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=1
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.3) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl5
	 X1End=X1Width/2.
	 X2End=X2Lvl4
	 NDiv=5
	 NodE=3
	 Gradient=0.8
	 NthGrp=3
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.4) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl4
	 X1End=-X1Width/2.
	 X2End=X2Lvl4
	 NDiv=10
	 NodE=3
	 Gradient=1.
	 NthGrp=4
	 NthBd=1
	 Mst2NthBd=2
	elseif(NthSegm.eq.5) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl4
	 X1End=X1Width/2.
	 X2End=X2Lvl3
	 NDiv=1
	 NodE=3
	 Gradient=1.
	 NthGrp=5
	 NthBd=2
	 Mst2NthBd=2
	elseif(NthSegm.eq.6) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl3
	 X1End=HeightQW
	 X2End=X2Lvl3
	 NDiv=5
	 NodE=3
	 Gradient=0.8
	 NthGrp=6
	 NthBd=2
	 Mst2NthBd=3
	elseif(NthSegm.eq.7) then
	 JBnd0Crk1Tip2=0
	 X1Stt=HeightQW
	 X2Stt=X2Lvl3
	 X1End=0.
	 X2End=X2Lvl3+HeightQW
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=7
	 NthBd=2
	 Mst2NthBd=3
	elseif(NthSegm.eq.8) then
	 JBnd0Crk1Tip2=0
	 X1Stt=0.
	 X2Stt=X2Lvl3+HeightQW
	 X1End=-HeightQW
	 X2End=X2Lvl3
	 NDiv=3
	 NodE=3
	 Gradient=1.
	 NthGrp=8
	 NthBd=2
	 Mst2NthBd=3
	elseif(NthSegm.eq.9) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-HeightQW
	 X2Stt=X2Lvl3
	 X1End=-X1Width/2.
	 X2End=X2Lvl3
	 NDiv=5
	 NodE=3
	 Gradient=1.25
	 NthGrp=9
	 NthBd=2
	 Mst2NthBd=3
	elseif(NthSegm.eq.10) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-X1Width/2.
	 X2Stt=X2Lvl3
	 X1End=-X1Width/2.
	 X2End=X2Lvl4
	 NDiv=1
	 NodE=3
	 Gradient=1.
	 NthGrp=10
	 NthBd=2
	 Mst2NthBd=2
	elseif(NthSegm.eq.11) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl3
	 X1End=X1Width/2.
	 X2End=X2Lvl2
	 NDiv=2
	 NodE=3
	 Gradient=1.
	 NthGrp=11
	 NthBd=3
	 Mst2NthBd=3
	elseif(NthSegm.eq.12) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl2
	 X1End=-X1Width/2.
	 X2End=X2Lvl2
	 NDiv=10
	 NodE=3
	 Gradient=1.
	 NthGrp=12
	 NthBd=3
	 Mst2NthBd=4
	elseif(NthSegm.eq.13) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-X1Width/2.
	 X2Stt=X2Lvl2
	 X1End=-X1Width/2.
	 X2End=X2Lvl3
	 NDiv=2
	 NodE=3
	 Gradient=1.
	 NthGrp=13
	 NthBd=3
	 Mst2NthBd=3
	elseif(NthSegm.eq.14) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl2
	 X1End=X1Width/2.
	 X2End=X2Lvl1
	 NDiv=1
	 NodE=3
	 Gradient=1.
	 NthGrp=14
	 NthBd=4
	 Mst2NthBd=4
	elseif(NthSegm.eq.15) then
	 JBnd0Crk1Tip2=0
	 X1Stt=X1Width/2.
	 X2Stt=X2Lvl1
	 X1End=-X1Width/2.
	 X2End=X2Lvl1
	 NDiv=20
	 NodE=3
	 Gradient=1.
	 NthGrp=15
	 NthBd=4
	 Mst2NthBd=4
	elseif(NthSegm.eq.16) then
	 JBnd0Crk1Tip2=0
	 X1Stt=-X1Width/2.
	 X2Stt=X2Lvl1
	 X1End=-X1Width/2.
	 X2End=X2Lvl2
	 NDiv=1
	 NodE=3
	 Gradient=1.
	 NthGrp=16
	 NthBd=4
	 Mst2NthBd=4
	endif
*============================================================
	return
	end
***************************************************
	Subroutine Ex4_UserGeom(NtotBndSegm,NEAddCap,NthSegm,
	1	JBnd0Crk1Tip2,X1Stt,X2Stt,X1End,X2End,NDiv,NodE, !output
	1	Gradient,NthGrp,NthBd,Mst2NthBd) !output
	Implicit Real*8 (A-H,O-Z)
	Implicit Integer (I-N)
*============================================================
	PI=4.d0*DATAN(1.d0)
*============================================================
	RIncl=0.1d0
	X1CtrIncl1=0.3d0
	X2CtrIncl1=0.2d0
	X1CtrIncl2=0.35d0
	X2CtrIncl2=-0.2d0
	X1CtrIncl3=0.6d0
	X2CtrIncl3=0.1d0
!
	NDivIncl=36
	NtotBndSegm=2+NDivIncl*3
!===================================================
	NEAddCap=100
!============================================================
	If(NthSegm.eq.1) then
	 JBnd0Crk1Tip2=2
	 X1Stt=-0.25d0
	 X2Stt=0.d0
	 X1End=0.d0
	 X2End=0.d0
	 NDiv=5
	 NodE=3
	 Gradient=0.9d0
	 NthGrp=101
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.eq.2) then
	 JBnd0Crk1Tip2=2
	 X1Stt=-0.25d0
	 X2Stt=0.d0
	 X1End=-0.5d0
	 X2End=0.d0
	 NDiv=3
	 NodE=3
	 Gradient=1.d0
	 NthGrp=102
	 NthBd=1
	 Mst2NthBd=1
	elseif(NthSegm.ge.3.and.NthSegm.le.2+NDivIncl) then
	 JBnd0Crk1Tip2=0
	 AngleStt=(NthSegm-3.)/NDivIncl*PI*2.
	 AngleEnd=(NthSegm-2.)/NDivIncl*PI*2.
	 X1Stt=X1CtrIncl1+RIncl*DCOS(AngleStt)
	 X2Stt=X2CtrIncl1+RIncl*DSIN(AngleStt)
	 X1End=X1CtrIncl1+RIncl*DCOS(AngleEnd)
	 X2End=X2CtrIncl1+RIncl*DSIN(AngleEnd)
	 NDiv=1
	 NodE=1
	 Gradient=1.d0
	 NthGrp=1
	 NthBd=2
	 Mst2NthBd=1
	elseif(NthSegm.ge.3+NDivIncl.and.NthSegm.le.2+NDivIncl*2) then
	 JBnd0Crk1Tip2=0
	 AngleStt=(NthSegm-3.d0-NDivIncl)/NDivIncl*PI*2.d0
	 AngleEnd=(NthSegm-2.d0-NDivIncl)/NDivIncl*PI*2.d0
	 X1Stt=X1CtrIncl2+RIncl*DCOS(AngleStt)
	 X2Stt=X2CtrIncl2+RIncl*DSIN(AngleStt)
	 X1End=X1CtrIncl2+RIncl*DCOS(AngleEnd)
	 X2End=X2CtrIncl2+RIncl*DSIN(AngleEnd)
	 NDiv=1
	 NodE=1
	 Gradient=1.d0
	 NthGrp=2
	 NthBd=3
	 Mst2NthBd=1
	elseif(NthSegm.ge.3+NDivIncl*2.and.NthSegm.le.2+NDivIncl*3) then
	 JBnd0Crk1Tip2=0
	 AngleStt=(NthSegm-3.d0-NDivIncl*2)/NDivIncl*PI*2.d0
	 AngleEnd=(NthSegm-2.d0-NDivIncl*2)/NDivIncl*PI*2.d0
	 X1Stt=X1CtrIncl3+RIncl*DCOS(AngleStt)
	 X2Stt=X2CtrIncl3+RIncl*DSIN(AngleStt)
	 X1End=X1CtrIncl3+RIncl*DCOS(AngleEnd)
	 X2End=X2CtrIncl3+RIncl*DSIN(AngleEnd)
	 NDiv=1
	 NodE=1
	 Gradient=1.d0
	 NthGrp=3
	 NthBd=4
	 Mst2NthBd=1
	endif
*============================================================
	return
	end