/BATCH,LIST
FINISH
/CLEAR
/PREP7

*AFUN,DEG

!--Input variables--

A=2.5		!Crack length
H=0.05		!Crack width
LT_X=3.75	!Left tip x coordinate
LT_Y=10		!Left tip y coordinate
THLGL=180	!Initial angle of the local coordinate system relative to the global coordinate system
THLGR=0		!Initial angle of the local coordinate system relative to the global coordinate system
B=20		!Slab height
W=10		!Slab width

MESHSIZE=0.15	!Element size
RADIUS=0.00833	!Radius of first row of elements
NUM=8		!Number of elements in the circunference
NSTEP=200	!Number of steps

!--Geometry--

! Material parameters
MPTEMP,,,,,,,, 										
MPTEMP,1,0										
MPDATA,EX,1,,474000	!Elastic module									
MPDATA,PRXY,1,,0.3 	!Poisson ratio

! Element type
ET,1,PLANE183 										
! Element key options
KEYOPT,1,1,0 		!8-node quadrilateral								
KEYOPT,1,3,2		!Plane strain
KEYOPT,1,6,1		!Mixed u-P formulation

! KP Geometry
! - Boundary
K,,W	,0	,0,	!Generates the Slab's keypoints
K,,W	,B/2    ,0,
K,,W	,B	,0,
K,,W/2	,B	,0,
K,,0	,B	,0,
K,,0	,B/2    ,0,
K,,0	,0  	,0,
K,,W/2	,0	,0,

FLST,2,8,3
FITEM,2,1
FITEM,2,2
FITEM,2,3
FITEM,2,4
FITEM,2,5
FITEM,2,6
FITEM,2,7
FITEM,2,8
A,P51X			!Generates the area
	
! - Crack
K,,LT_X+H*cos(30)	,LT_Y-0.5*H	,0,	!LL
K,,LT_X			,LT_Y		,0,	!LT
K,,LT_X+H*cos(30)	,LT_Y+0.5*H	,0,	!LR

K,,LT_X+A-H*cos(30)	,LT_Y+0.5*H	,0,	!RL
K,,LT_X+A		,LT_Y		,0,	!RT
K,,LT_X+A-H*cos(30)	,LT_Y-0.5*H	,0,	!RR

FLST,2,6,3
FITEM,2,9
FITEM,2,10
FITEM,2,11
FITEM,2,12
FITEM,2,13
FITEM,2,14
A,P51X			!Generates the crack's area

ASBA,1,2		!Substract the crack's from the slab's area

!Crack tip circle
K,,LT_X+0.5*H*cos(30)	,LT_Y-0.25*H	,0,	!LLC
K,,LT_X-0.5*H		,LT_Y		,0,	!LTC
K,,LT_X+0.5*H*cos(30)	,LT_Y+0.25*H	,0,	!LRC
K,,LT_X+A-0.5*H*cos(30)	,LT_Y+0.25*H	,0,	!RLC
K,,LT_X+A+0.5*H		,LT_Y		,0,	!RTC	
K,,LT_X+A-0.5*H*cos(30)	,LT_Y-0.25*H	,0,	!RRC

*GET,KPc,KP,0,COUNT

LT=KPc-10		!Defines the left tip keypoint
RT=KPc-7		!Defines the right tip keypoint

LSTR,KPc-10,KPc-4				
LARC,KPc-5,KPc-4,KPc-10	,H/2
LARC,KPc-4,KPc-3,KPc-10	,H/2

LSTR,KPc-7,KPc-1
LARC,KPc-2,KPc-1,KPc-7	,H/2
LARC,KPc-1,KPc	,KPc-10	,H/2

*GET,Lc,LINE,0,COUNT

LSBL,Lc-11,	Lc-4,,	DELETE,KEEP
LSBL,Lc-10,	Lc-3,,	DELETE,KEEP
LSBL,Lc-8,	Lc-1,,	DELETE,KEEP
LSBL,Lc-7,	Lc,,	DELETE,KEEP

AL,Lc-4	,Lc-5,Lc+2
AL,Lc-11,Lc-5,Lc-3
AL,Lc-2	,Lc-1,Lc+4
AL,Lc-8	,Lc-2,Lc
AOVLAP,ALL

!--Initial meshing--

FLST,5,6,4,ORDE,6
FITEM,5,Lc+5
FITEM,5,Lc+3
FITEM,5,Lc+1
FITEM,5,Lc-2
FITEM,5,Lc-5
FITEM,5,Lc-10
CM,_Y,LINE
LSEL, , , ,P51X
CM,_Y1,LINE
CMSEL,,_Y
LESIZE,_Y1, , ,4,3, , , ,1
FLST,5,4,4,ORDE,4
FITEM,5,Lc-4
FITEM,5,Lc-3
FITEM,5,Lc-1
FITEM,5,Lc
CM,_Y,LINE
LSEL, , , ,P51X
CM,_Y1,LINE
CMSEL,,_Y
LESIZE,_Y1, , ,8, , , , ,1 

KSCON,LT,RADIUS,1,NUM,,
KSCON,RT,RADIUS,1,NUM,,

ESIZE,MESHSIZE,0,
AMESH,ALL

!--------------------------DO-LOOP-------------------------

*DO,i,1,NSTEP

!--Boundary conditions--

/SOL

FLST,2,2,4,ORDE,2   
FITEM,2,7   
FITEM,2,8  
/GO 
DL,P51X, ,UY,0   	!Fixed Y

FLST,2,2,4,ORDE,2   
FITEM,2,7   
FITEM,2,8  
/GO 
DL,P51X, ,UX,0   	!Fixed X

FLST,2,2,4,ORDE,2   
FITEM,2,3   
FITEM,2,4  
/GO 
SFL,P51X,PRES,-18000,   !Load

!--Solver operation

!--Define crack tip--
*GET,LT_X,KP,LT,LOC,X
*GET,LT_Y,KP,LT,LOC,Y
NSEL,S,LOC,X,LT_X
NSEL,R,LOC,Y,LT_Y
CM,LCTA,NODE			!Left crack tip area
*GET,LTN,NODE,LCTA,NUM,MIN	!Left tip node
KSEL,ALL

*GET,RT_X,KP,RT,LOC,X
*GET,RT_Y,KP,RT,LOC,Y
NSEL,S,LOC,X,RT_X
NSEL,R,LOC,Y,RT_Y
CM,RCTA,NODE			!Right crack crack tip area
*GET,RTN,NODE,RCTA,NUM,MIN	!Right crack tip node
KSEL, ALL

LOCAL,1001,0,LT_X,LT_Y,0,THLGL,0,0
ALLSEL,ALL

CINT,DELE,1

!--Contour integral
/SOL
CINT,NEW,1
CINT,TYPE,SIFS
CINT,CTNC,LCTA
CINT,NORM,1001,2
CINT,NCON,6

SOLVE

*GET,K1L_1,CINT,1,CTIP,LTN,CONTOUR,1,DTYPE,K1
*GET,K1L_2,CINT,1,CTIP,LTN,CONTOUR,2,DTYPE,K1
*GET,K1L_3,CINT,1,CTIP,LTN,CONTOUR,3,DTYPE,K1
*GET,K1L_4,CINT,1,CTIP,LTN,CONTOUR,4,DTYPE,K1
*GET,K1L_5,CINT,1,CTIP,LTN,CONTOUR,5,DTYPE,K1
*GET,K1L_6,CINT,1,CTIP,LTN,CONTOUR,6,DTYPE,K1
*GET,K2L_1,CINT,1,CTIP,LTN,CONTOUR,1,DTYPE,K2
*GET,K2L_2,CINT,1,CTIP,LTN,CONTOUR,2,DTYPE,K2
*GET,K2L_3,CINT,1,CTIP,LTN,CONTOUR,3,DTYPE,K2
*GET,K2L_4,CINT,1,CTIP,LTN,CONTOUR,4,DTYPE,K2
*GET,K2L_5,CINT,1,CTIP,LTN,CONTOUR,5,DTYPE,K2
*GET,K2L_6,CINT,1,CTIP,LTN,CONTOUR,6,DTYPE,K2

LOCAL,1002,0,RT_X,RT_Y,0,THLGR,0,0
ALLSEL,ALL

CINT,DELE,2

/SOL
CINT,NEW,2
CINT,TYPE,SIFS
CINT,CTNC,RCTA
CINT,NORM,1002,2
CINT,NCON,6

SOLVE

*GET,K1R_1,CINT,2,CTIP,RTN,CONTOUR,1,DTYPE,K1
*GET,K1R_2,CINT,2,CTIP,RTN,CONTOUR,2,DTYPE,K1
*GET,K1R_3,CINT,2,CTIP,RTN,CONTOUR,3,DTYPE,K1
*GET,K1R_4,CINT,2,CTIP,RTN,CONTOUR,4,DTYPE,K1
*GET,K1R_5,CINT,2,CTIP,RTN,CONTOUR,5,DTYPE,K1
*GET,K1R_6,CINT,2,CTIP,RTN,CONTOUR,6,DTYPE,K1
*GET,K2R_1,CINT,2,CTIP,RTN,CONTOUR,1,DTYPE,K2
*GET,K2R_2,CINT,2,CTIP,RTN,CONTOUR,2,DTYPE,K2
*GET,K2R_3,CINT,2,CTIP,RTN,CONTOUR,3,DTYPE,K2
*GET,K2R_4,CINT,2,CTIP,RTN,CONTOUR,4,DTYPE,K2
*GET,K2R_5,CINT,2,CTIP,RTN,CONTOUR,5,DTYPE,K2
*GET,K2R_6,CINT,2,CTIP,RTN,CONTOUR,6,DTYPE,K2

/POST1

!--Crack increment and direction

!--Effective SIF
K1L=(K1L_3+K1L_4+K1L_5)/3
K2L=(K2L_3+K2L_4+K2L_5)/3
K1R=(K1R_3+K1R_4+K1R_5)/3
K2R=(K2R_3+K2R_4+K2R_5)/3

DTHL=-2*atan((2*K2L)/(K1L +sqrt(K1L**2+8*K2L**2)))
DTHR=-2*atan((2*K2R)/(K1R +sqrt(K1R**2+8*K2R**2)))

THLGL=THLGL+DTHL
THLGR=THLGR+DTHR

LDA=0.1
RDA=0.1

/PREP7

!-- Define new geometry

KSEL,ALL
CSDEL,ALL

*GET,LT_X,KP,LT,LOC,X !Noetig?
*GET,LT_Y,KP,LT,LOC,Y
*GET,RT_X,KP,RT,LOC,X
*GET,RT_Y,KP,RT,LOC,Y

ACLEAR,ALL

A,LT-1,LT,LT+1
A,RT-1,RT,RT+1

NUMCMP,AREA

FLST,2,7,5,ORDE,7
FITEM,2,1
FITEM,2,2
FITEM,2,3
FITEM,2,4
FITEM,2,5
FITEM,2,6
FITEM,2,7
AADD,P51X

LDX=LDA*cos(THLGL)
LDY=LDA*sin(THLGL)
RDX=RDA*cos(THLGR)
RDY=RDA*sin(THLGR)

!*GET,KPc,KP,0,COUNT
KPnum=30+i*20
NUMSTR,KP,KPnum

! New crack tip
K,,LT_X+LDX+H*(cos(150)*cos(THLGL)-sin(150)*sin(THLGL))	,LT_Y+LDY+H*(cos(150)*sin(THLGL)+sin(150)*cos(THLGL))	,0
K,,LT_X+LDX						,LT_Y+LDY						,0
K,,LT_X+LDX+H*(cos(210)*cos(THLGL)-sin(210)*sin(THLGL))	,LT_Y+LDY+H*(cos(210)*sin(THLGL)+sin(210)*cos(THLGL))	,0

K,,RT_X+RDX+H*(cos(150)*cos(THLGR)-sin(150)*sin(THLGR))	,RT_Y+RDY+H*(cos(150)*sin(THLGR)+sin(150)*cos(THLGR))	,0
K,,RT_X+RDX						,RT_Y+RDY						,0
K,,RT_X+RDX+H*(cos(210)*cos(THLGR)-sin(210)*sin(THLGR))	,RT_Y+RDY+H*(cos(210)*sin(THLGR)+sin(210)*cos(THLGR))	,0

*GET,Lc,LINE,0,count 
Lnum = Lc+i*40 
NUMSTR,LINE,Lnum

A,LT-1,KPnum,KPnum+1,KPnum+2,LT+1
A,RT-1,KPnum+3,KPnum+4,KPnum+5,RT+1
NUMCMP,AREA
ASBA,3,1
ASBA,4,2

NUMMRG,KP
NUMCMP,KP

*GET,LT,LINE,Lnum+1,KP,2
*GET,RT,LINE,Lnum+5,KP,2
*GET,KPc,KP,0,COUNT

K,,LT_X+LDX+0.5*H*(cos(150)*cos(THLGL)-sin(150)*sin(THLGL))	,LT_Y+LDY+0.5*H*(cos(150)*sin(THLGL)+sin(150)*cos(THLGL))	,0
K,,LT_X+LDX+0.5*H*(cos(0)*cos(THLGL)-sin(0)*sin(THLGL))		,LT_Y+LDY+0.5*H*(cos(0)*sin(THLGL)+sin(0)*cos(THLGL))		,0
K,,LT_X+LDX+0.5*H*(cos(210)*cos(THLGL)-sin(210)*sin(THLGL))	,LT_Y+LDY+0.5*H*(cos(210)*sin(THLGL)+sin(210)*cos(THLGL))	,0

K,,RT_X+RDX+0.5*H*(cos(150)*cos(THLGR)-sin(150)*sin(THLGR))	,RT_Y+RDY+0.5*H*(cos(150)*sin(THLGR)+sin(150)*cos(THLGR))	,0
K,,RT_X+RDX+0.5*H*(cos(0)*cos(THLGR)-sin(0)*sin(THLGR))		,RT_Y+RDY+0.5*H*(cos(0)*sin(THLGR)+sin(0)*cos(THLGR))		,0
K,,RT_X+RDX+0.5*H*(cos(210)*cos(THLGR)-sin(210)*sin(THLGR))	,RT_Y+RDY+0.5*H*(cos(210)*sin(THLGR)+sin(210)*cos(THLGR))	,0

LSTR,LT,KPc+2
LARC,KPc+1,KPc+2,LT,0.025,
LARC,KPc+2,KPc+3,LT,0.025,
LSTR,RT,KPc+5
LARC,KPc+4,KPc+5,RT,0.025,
LARC,KPc+5,KPc+6,RT,0.025,

LSBL,LNUM+1,LNUM+9,,DELETE,KEEP
LSBL,LNUM+2,LNUM+10,,DELETE,KEEP
LSBL,LNUM+5,LNUM+12,,DELETE,KEEP
LSBL,LNUM+6,LNUM+13,,DELETE,KEEP

AL,LNUM+15,LNUM+9,LNUM+8
AL,LNUM+1,LNUM+10,LNUM+8
AL,LNUM+17,LNUM+12,LNUM+11
AL,LNUM+5,LNUM+13,LNUM+11
AOVLAP,ALL

FLST,5,6,4,ORDE,6
FITEM,5,LNUM+2
FITEM,5,LNUM+8
FITEM,5,LNUM+11
FITEM,5,LNUM+14
FITEM,5,LNUM+16
FITEM,5,LNUM+18
CM,_Y,LINE
LSEL, , , ,P51X
CM,_Y1,LINE
CMSEL,,_Y
LESIZE,_Y1, , ,4,3, , , ,1
FLST,5,4,4,ORDE,4
FITEM,5,LNUM+9
FITEM,5,LNUM+10
FITEM,5,LNUM+12
FITEM,5,LNUM+13
CM,_Y,LINE
LSEL, , , ,P51X
CM,_Y1,LINE
CMSEL,,_Y
LESIZE,_Y1, , ,8, , , , ,1 

KSCON,LT,RADIUS,1,NUM,,
KSCON,RT,RADIUS,1,NUM,,

ESIZE,MESHSIZE,0,
AMESH,ALL

*ENDDO