*Heading
** Job name: InputFileTemplate
***************************************************************
** Abaqus Code for a Residual Control Staggered Solution Scheme for the Phase-Field Modeling of Brittle Fracture
** Version 3, DOI: 10.17632/p77tsyrbx2.3
** https://data.mendeley.com/datasets/p77tsyrbx2/3
**
** K. Seles, 2019
** karlo.seles@fsb.hr
***************************************************************
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
*Part, name=Beam
*Node
[YOUR INPUT HERE - node number, x coord, y coord]
***************************************************************
** UEL (U1) - the first layer of user phase-field elements
** 2D 4-node full integration linear element
** Properties - length scale parameter (xL) and Griffiths critical energy release rate (Gc)
**
*User element, nodes=4, type=U1, properties=2, coordinates=2
***************************************************************
*Element, type=U1
[YOUR INPUT HERE - Element number, element conectivity - 4 nodes]
***************************************************************
** UMAT (CPE4) - the 2nd layer of Abaqus standard displacement elements
**
** 2D 4-node full integration linear element for plane strain
** Properties - Young's modulus (EMOD) and Poisson's ratio (ENU)
**
*Element, TYPE=CPE4, elset=Displacement
[YOUR INPUT HERE - use same element conectivity as in first layer, but change the element numbers as "Element number + first layer element number count"]
***************************************************************
** UEL (U3) - the third layer of user elements
** 2D 4-node full integration linear element with identity stiffness matrix and dummy DOFs
** Properties - length scale parameter (xL) and Griffiths critical energy release rate (Gc)
**
** Variables - History parameter of maximum deformation energy (1x4 integration points)
**
*User element, nodes=4, type=U3, properties=2, coordinates=2, VARIABLES=4
***************************************************************
*Element, type=U3
[YOUR INPUT HERE - use same element conectivity as in first layer, but change the element numbers as "Element number + total first layer element number count + total second layer element number count"]
***************************************************************
*Nset, nset=Set-1, generate
[YOUR INPUT HERE - 1, total number of nodes, 1]
***************************************************************
** Creating layer element sets
*Elset, elset=Phase, generate
[YOUR INPUT HERE - 1, total first layer element number count, 1]
*Elset, elset=Displacement, generate
[YOUR INPUT HERE - total first layer element number count + 1, total first layer element number count + total second layer element number count, 1]
*Elset, elset=Phase_2nd, generate
[YOUR INPUT HERE - total first layer element number count + total second layer element number count + 1, total first layer element number count + total second layer element number count + total third layer element number count, 1]
***************************************************************
*Solid Section, elset=Displacement, material=Displacement
1.0
***************************************************************
** Material properties for the 1st and 3rd element layer
** Note: Must be equal
*Uel property, elset=Phase
[YOUR INPUT HERE - xL(in mm if geometry is in mm), gc(in kN/mm if geometry is in mm)]
*Uel property, elset=Phase_2nd
[YOUR INPUT HERE - xL(in mm if geometry is in mm), gc(in kN/mm if geometry is in mm)]
**xL, gc
***************************************************************
*End Part
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Beam, part=Beam
*End Instance
***************************************************************
** Creating Boundary Conditions Node Sets
*Nset, nset=name, instance=beam
[YOUR INPUT HERE]
** Node for the phase-field DOF monitor in status file (*.sta)
*Nset, nset=monitor, instance=beam
[YOUR INPUT HERE]
*End Assembly
***************************************************************
** Material properties for the 2nd element layer
** Elastic properties for Displacement UMAT
*Material, name=Displacement
*Depvar
*User Material, constants=2
[YOUR INPUT HERE]
**[E(in kN/mm^2 if geometry is in mm),   nu]
***************************************************************
** BOUNDARY CONDITIONS
*Boundary
[YOUR INPUT HERE - nodesetname, DOF, DOF]
***************************************************************
** STEP: Static
** Important to set "extrapolation=no"!
*Step, name=Static, nlgeom=NO, inc=1000000, extrapolation=no
*Static
[YOUR INPUT HERE - initial time step size, total time, minimal time step size, maximal time step size]
[YOUR INPUT HERE - *Boundary for displacement loading or *Cload for Force loading] *Boundary
[YOUR INPUT HERE - nodesetname, DOF, value]
************************************************************************
*Controls, reset
*Controls, parameters=time incrementation
 7001, 7001, 7001, 7000, 1500, 400
 , , 0.75, , , , 1.5, ,
************************************************************************
** Explanations:
** Second line:
** IMPORTANT - 1), and 2) have to be set to a higher number than 4) to allow the non-monotonic convergence typical for staggered phase-field fracture models!
** 4) is a maximum allowed number of iterations - if the system does not converge in 4) number of iterations, the increment size is reduced and the increment is attempted again!
** 5) and 6) are used for automatic loading incrementation together with 3) and 7) in the third line!
*****
** 1) - I_0 iteration after which a check is made whether the residuals increase in two consecutive iterations (default: 4). If so, the increment is reattempted with D_f times its size.
** 2) - I_R iteration after which a logarithmic convergence check is performed in each iteration (default: 8). If more than I_C iterations are needed, the increment is reattempted with D_C its size.
** 3) - I_P iteration after which the residual tolerance R^alpha_p is used instead of R^\alpha_n (default: 9). - Can be used to lower the tolerance after some number of iteration.
** 4) - I_C maximum number of iterations allowed (default: 16).
** 5) - I_L number of iterations after which the size of the subsequent increment will be reduced (default: 10).
** 6) - I_G maximum number of iterations allowed in two consecutive increments for the size of the next increment to be increased (default: 4).
*****
** Third line:
** 3) - D_B Cutback factor for the next increment if more than I_L iterations were needed in the current increment (default: 0.75).
** 7) - D_D Factor by which the next increment will be increased if less than I_G iterations are needed in two consecutive increments (default: 1.5).
******************************************************************
*Controls, parameters=field, field=DISPLACEMENT
,,,,,1e-14,1e-14
*Controls, parameters=field, field=ROTATION
,,,,,1e-14,1e-14
******************************************************************
** Explanation:
** IMPORTANT: 6) and 7) are set to very small numbers to avoid "zero force" and "zero moment" calculation state!
** In this version of the code we have opted for the Abaqus standard convergance tolerance - however, one can also use the
** fixed tolerance as used in previous version of the code, and in the paper. To use that, 4) has to be set as it is "user defined value of average moment norm" and then, 1) is multiplied with 4) and fixed tolerance is set!
*****
**1) - R^alpha_n Convergence criterion for the ratio of the largest residual to the average moment  (default: 0.005). - in this case, average moment is substituted with user-defined value of average moment norm.
**4) - q^alpha_u user-defined average force. If defined, the calculation of the average force is replaced by this value.
**5) - R^alpha_p Alternative residual convergence criterion to be used after I_P iterations instead of R^alpha_n (default: 0.02).
**6) - epsilon^alpha Criterion for zero flux relative to q^alpha (default: 10^-5).
**7) - C^alpha_epsilon Convergence criterion for the ratio of the largest solution correction to the largest incremental solution value in case of zero flux (default: 10^-3).
************************************************************************
** OUTPUT **
*MONITOR, NODE=Monitor, DOF=4
*Restart, write, frequency=0
*Output, field, frequency=1
*element output, elset=BEAM.DISPLACEMENT
S,E,SENER,SDV
*Node Output
RF, U
*End Step