*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
 2,
 bot, 2, 1.0, <db>, 2, -1.0 
*equation
 2,
 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
 1,0,0, 0,1,0
*material, name=elas
*elastic
 3.24E9, 0.3
*damage initiation, criterion=maxps,tol=0.001
 <ultI>
*damage evolution,type=energy,mixed mode behavior=bk,power=<eta>
 <GIc>, <GIIc>
*damage stabilization
 1.e-5
*enrichment, type=propagation crack, name=enr1, elset=solid
*boundary
 xsymm, 1
*boundary
 tdisp, 1, 2
 bdisp, 1, 2
 right,1,1
**
*elset,elset=hole,generate
1041,1063
*surface,name=hole_inner,type=element
hole,s4
*surface,name=total,type=xfem
enr1
*step, nlgeom,inc=100000
*static
 0.01, 1.0, 1.e-9, 0.01
*controls, analysis=discontinuous
*controls,parameters=time incrementation
 , , , , , , ,20
*dsload
hole_inner,PNU
total,PNU
*output,field, freq=25
*node output
 u,philsm
*element output
 s, statusxfem
*output,history
*node output, nset=rdisp
 u2, rf2
*end step