*heading
 Crack initiation and propagation in a plate with a hole (XFEM)
**
*preprint, model=yes, echo=yes, history=yes
**
*parameter
** Interface properties
** Ultimate strenth in tensile and mode II:
 ultI   =  22.0E6
 ultII  =  22.0E6
** Fracture toughness (N/m):
 GIc   = 2.87E3
 GIIc  = 2.87E3
** B-K parameter:
 eta = 1.0
**
** Dimensions
**
*parameter
 length = 100.0e-3
** Width in the plane strain direction
 width = 20.0e-3
** Hole radius
 radius = length/5
** Height
 height = 3*length + 2*radius
** A small number in comparison to the model dimensions
 small = 1.e-8*length
**
 hlen = length/2
 ht2 = height/2
**
*parameter
** Keep ex below 0.5*ey_a to minimize distortion
 ex = 20
** Number of elements in the vertical direction in the three model blocks
 ey_a = 40
 ey_b = 6
** ey_c should be odd for the crack to cut through the center of an element
 ey_c = 11
**
  ey = 2*(ey_a + ey_b) + ey_c
**
  hinc = 1000
**
  nl1 = 1
  nl2 = nl1 + ey_a
  nl3 = nl2 + ey_b
  nl4 = nl3 + ey_c
  nl5 = nl4 + ey_b
  nl6 = nl5 + ey_a
**
  nr1 = nl1 + ex*hinc
  nr2 = nl2 + ex*hinc
  nr3 = nl3 + ex*hinc
  nr4 = nl4 + ex*hinc
  nr5 = nl5 + ex*hinc
  nr6 = nl6 + ex*hinc
**
  db  = 9000001
  dt  = 9000002
**
** Coordinate system offset
  yoff = -ht2
  yl1 = 0 + yoff
  yl2 = yl1 + 1.5*length
  yl3 = yl2 + 0.5*radius
  yl4 = yl3 + radius
  yl5 = yl4 + 0.5*radius
  yl6 = yl5 + 1.5*length
**
  yr1 = yl1
  yr2 = yl2 - 0.5*length
  yr3 = yl3
  yr4 = yl4
  yr5 = yl5 + 0.5*length
  yr6 = yl6
** x-crd given by cos(30) * radius at locations 3 and 4
  xl3 = 0.866 * radius
  xl4 = xl3
**
  elab = 1001
  e1_1 = nl1
  e1_2 = nl1 + hinc
  e1_3 = nl1 + 1 + hinc
  e1_4 = nl1 + 1
**
*node
 <nl1>, 0., <yl1>
 <nr1>, <length>, <yr1>
 <nl2>, 0., <yl2>
 <nr2>, <length>, <yr2>
 <nl3>, <xl3>, <yl3>
 <nr3>, <length>, <yr3>
 <nl4>, <xl4>, <yl4>
 <nr4>, <length>, <yr4>
 <nl5>, 0., <yl5>
 <nr5>, <length>, <yr5>
 <nl6>, 0., <yl6>
 <nr6>, <length>, <yr6>
 <db>,    <hlen>,  <yl1>
 <dt>,    <hlen>,  <yl6>
**
*ngen, nset=left
 <nl1>, <nl2>
*ngen, nset=left, line=c
 <nl2>, <nl3>, , , 0., 0., 0., 0., 0., 1.0
*ngen, nset=left, line=c
 <nl3>, <nl4>, , , 0., 0., 0., 0., 0., 1.0
*ngen, nset=left, line=c
 <nl4>, <nl5>, , , 0., 0., 0., 0., 0., 1.0
*ngen, nset=left
 <nl5>, <nl6>
**
*ngen, nset=right
 <nr1>, <nr2>
*ngen, nset=right
 <nr2>, <nr3>
*ngen, nset=right
 <nr3>, <nr4>
*ngen, nset=right
 <nr4>, <nr5>
*ngen, nset=right
 <nr5>, <nr6>
**
*nfil, nset=nall
 left, right, <ex>, <hinc>
**
*nset,nset=tdisp
 <dt>
*nset,nset=bdisp
 <db>
*nset,nset=rdisp
 tdisp, bdisp
*nset, nset=bot, generate
<nl1>, <nr1>, <hinc>
*nset, nset=top, generate
<nl6>, <nr6>, <hinc>
*nset, nset=xsymm, generate
<nl1>, <nl2>
<nl5>, <nl6>
**
*equation
 bot, 2, 1.0, <db>, 2, -1.0
*equation
 top, 2, 1.0, <dt>, 2, -1.0
**
*element, type=CPE4
 <elab>, <e1_1>, <e1_2>, <e1_3>, <e1_4>
*elgen,elset=solid
 <elab>, <ex>, <hinc>, <hinc>, <ey>, 1, 1
*solid section, elset=solid, mat=elas, orient=global
 <width>
*orientation, name=global
*material, name=elas
*elastic
 3.24E9, 0.3
*damage initiation, criterion=maxps,tol=0.01
 <ultI>
*enrichment, type=propagation crack, name=enr1, elset=solid,interaction=lower
*surface interaction,name=lower
*surface behavior
*FRACTURE CRITERION,TYPE=VCCT,MIXED MODE BEHAVIOR=POWER,tol=0.1,viscosity=1.0E-5
<GIc>, <GIIc>,0,1.0,1.0,1.0
*boundary
 xsymm, 1
**
*step, nlgeom,inc=100000
*static
 0.01, 1.0, 1.e-9, 0.01
*controls, analysis=discontinuous
*controls,parameters=time incrementation
, , , , , , ,20
*boundary
 tdisp, 2, 2,  5.5e-4
 bdisp, 2, 2, -5.5e-4
*output,field, freq=50
*node output
 u,philsm
*element output
 s, statusxfem
*output,history,freq=50
*node output, nset=rdisp
 u2, rf2
*end step