# === Modules ===

from part import *
from material import *
from section import *
from assembly import *
from step import *
from regionToolset import *
from interaction import *
from load import *
from mesh import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *


#---------------------------------------------------------------------------------

# === Parameters ===

modelName = 'modeII_3d_lefm'

length = 3.0                              # plate half-length
height = 6.0                              # plate height
width  = 1.0                              # plate width

crack_len = 1.5                           # crack length
crack_y   = 0.083333                      # crack offset in y-direction

YM    = 2.10e+11                          # Young's modulus
NU    = 0.3                               # Poisson's ratio
MAXPS = 8.44E7                            # Damage initiation
DTOL  = 0.05                              # Damage tolerance
GI    = 42200.0                           # Fracture energy
ETA   = 1.0                               # Power-law exponent


#---------------------------------------------------------------------------------

# === Create model ===

Mdb()

viewportName = session.Viewport(name=modelName)
viewportName.makeCurrent()
viewportName.maximize()

plateModel = mdb.Model(name=modelName)
del mdb.models['Model-1']


#---------------------------------------------------------------------------------

# === Create parts ===

plateSketch = plateModel.ConstrainedSketch(name='plateProfile',sheetSize=height)

plateSketch.setPrimaryObject(option=STANDALONE)
plateSketch.sketchOptions.setValues(decimalPlaces=3)

plateSketch.Line(point1=(0.0, -height/2.0), point2=(0.0, height/2.0))
plateSketch.Line(point1=(0.0, height/2.0), point2=(length, height/2.0))
plateSketch.Line(point1=(length, height/2.0), point2=(length, -height/2.0))
plateSketch.Line(point1=(length, -height/2.0), point2=(0.0, -height/2.0))

platePart = plateModel.Part(dimensionality=THREE_D, name='plate', 
                type=DEFORMABLE_BODY)
platePart.BaseSolidExtrude(sketch=plateModel.sketches['plateProfile'], depth=1.0)

plateSketch.unsetPrimaryObject()
del plateSketch


# Part for crack geometry

crackSketch =  plateModel.ConstrainedSketch(name='crackProfile',sheetSize=height)

crackSketch.Line(point1=(0, crack_y), point2=(crack_len, crack_y))

crackPart = plateModel.Part(dimensionality=THREE_D, name='crack',
                type=DEFORMABLE_BODY)
crackPart.BaseShellExtrude(sketch=plateModel.sketches['crackProfile'], depth=1.0)

crackSketch.unsetPrimaryObject()
del crackSketch

#---------------------------------------------------------------------------------

# === Create Partition ===

plate_part = mdb.models[modelName].parts['plate']
faces=plate_part.faces
plate_part.DatumPlaneByOffset(plane=faces[2], flip=SIDE2, offset=length/2.0 )

plate_part = mdb.models[modelName].parts['plate']
faces=plate_part.faces
plate_part.DatumPlaneByOffset(plane=faces[1], flip=SIDE2, offset=height/2.0 )

pickedCells= platePart.cells.findAt((0.0,0,0.0))
datums1 = plate_part.datums
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[2], cells=pickedCells)

pickedCells= platePart.cells.findAt(((0.0,0,0.0),),  ((length,0,0.0),), ) 
datums1 = plate_part.datums
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[3], cells=pickedCells)
#---------------------------------------------------------------------------------


# === Create geometry sets ===

platePart.Set(cells=platePart.cells[:], name='All')

pickface = platePart.faces.findAt(( (length/2.0,  -height/2.0, width/2.0), ),( (length*3/4,  -height/2.0, width/2.0), ))
platePart.Set(faces=pickface , name='bottom')

pickface  = platePart.faces.findAt( ((length/4.0,  +height/2.0, width/2.0), ),((length/2.0,  +height/2.0, width/2.0), )  )
platePart.Set(faces=pickface , name='top')

v1 = platePart.vertices.findAt(( (0.0,  -height/2.0, 0.0), ),
         ((0.0, +height/2.0, 0.0), ), ((length,  -height/2.0, 0.0), ),
         ((length,  +height/2.0, 0.0), ))
platePart.Set(vertices=v1, name='fixZ')

#---------------------------------------------------------------------------------

# === Define material and section properties ===

plateMatl = plateModel.Material(name='elas')

plateMatl.Elastic(table=((YM, NU), ))
plateMatl.MaxpsDamageInitiation(table=((MAXPS, ), ), tolerance=DTOL)


plateModel.HomogeneousSolidSection(material='elas', name='solid')

#---------------------------------------------------------------------------------

# === Assign section and orientation ===

platePart.MaterialOrientation(fieldName='', localCsys=None, orientationType=GLOBAL, region=
    platePart.sets['All'], stackDirection=STACK_3)

platePart.SectionAssignment(region=platePart.sets['All'], sectionName='solid')

#---------------------------------------------------------------------------------

# === Assign mesh controls and mesh plate ===

# Element type
platePart.setElementType(elemTypes=(ElemType(elemCode=C3D8R, elemLibrary=STANDARD), 
    ), regions=platePart.sets['All'])

pickedEdges=platePart.edges.findAt(
  ((0.0,+height/2.0,width/2  ),), ((0.0,-height/2.0,width/2  ),), ((0.0,0,width/2  ),), 
   )
platePart.seedEdgeByNumber(edges=pickedEdges, number=8, constraint=FIXED) 

pickedEdges=platePart.edges.findAt(
  ((length/2.0,+height/2.0,width/2.0  ),), ((length/2.0,-height/2.0,width/2.0  ),), ((length/2.0,0,width/2.0  ),), 
  ((length,+height/2.0,width/2.0  ),), ((length,-height/2.0,width/2.0  ),), ((length,0,width/2.0  ),), 
   )
platePart.seedEdgeByNumber(edges=pickedEdges, number=8, constraint=FIXED) 

pickedEdges=platePart.edges.findAt(
    ((0,-height/4.0,0  ),),  ((length/2.0,-height/4.0,0  ),), ((length,-height/4.0,0  ),),
    ((0,-height/4.0,width  ),), ((length/2.0,-height/4.0,width  ),), ((length,-height/4.0,width    ),),
   )
platePart.seedEdgeByNumber(edges=pickedEdges, number=12, constraint=FIXED) 

pickedEdges=platePart.edges.findAt(
    ((0,height/4.0,0  ),),  ((length/2.0,height/4.0,0  ),), ((length,height/4.0,0  ),),
    ((0,height/4.0,width  ),), ((length/2.0,height/4.0,width  ),), ((length,height/4.0,width    ),),
   )
platePart.seedEdgeByNumber(edges=pickedEdges, number=18, constraint=FIXED) 

pickedEdges=platePart.edges.findAt(
    ((length*3/4,0,0  ),),  ((length*3/4.0,height/2.0,0  ),), ((length*3/4,-height/2,0  ),),
    ((length*3/4,0,width    ),),  ((length*3/4.0,height/2.0,width    ),), ((length*3/4,-height/2,width    ),),
   )
platePart.seedEdgeByNumber(edges=pickedEdges, number=13, constraint=FIXED) 

pickedEdges=platePart.edges.findAt(
    ((length/8,0,0  ),),  ((length/8.0,height/2.0,0  ),),  
      ((length/8.0,height/2.0,width    ),),  
   )

platePart.seedEdgeByBias(biasMethod=SINGLE, end2Edges=pickedEdges, ratio=2.0, 
    number=5, constraint=FINER)

pickedEdges=platePart.edges.findAt(
      ((length/8,-height/2,0  ),),
    ((length/8,0,width    ),),   ((length/8,-height/2,width    ),),
   )

platePart.seedEdgeByBias(biasMethod=SINGLE, end1Edges=pickedEdges, ratio=2.0, 
    number=5, constraint=FINER)
platePart.generateMesh()

#---------------------------------------------------------------------------------

# === Assemble ===

plateModel.rootAssembly.DatumCsysByDefault(CARTESIAN)
plateModel.rootAssembly.Instance(dependent=ON, name='plate_1', part=platePart)
plateModel.rootAssembly.Instance(dependent=ON, name='crack_1', part=crackPart)

# Reference points for displacement application
rp_db = plateModel.rootAssembly.ReferencePoint(point=(length/2.0, -height/2.0, 
    0.0))
plateModel.rootAssembly.features.changeKey(fromName='RP-1', toName='db')

#---------------------------------------------------------------------------------

# === Create assembly sets ===

v1 = (plateModel.rootAssembly.referencePoints[rp_db.id], )
plateModel.rootAssembly.Set(name='bdisp', referencePoints=v1)

# === End assembly ===

#---------------------------------------------------------------------------------

# === Create constraint equation ===

plateModel.Equation(name='ce_bot', terms=((1.0, 'plate_1.bottom', 1), (-1.0, 'bdisp', 1)))

#---------------------------------------------------------------------------------

# === Create step ===

plateModel.StaticStep(timePeriod=0.5,initialInc=0.05, maxInc=0.05, maxNumInc=1000,
    minInc=1e-09, name='Static', nlgeom=ON, previous='Initial')

plateModel.steps['Static'].control.setValues(
    allowPropagation=OFF,  resetDefaultValues=OFF, 
    timeIncrementation=(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 20.0, 6.0, 3.0, 
    50.0))

plateModel.fieldOutputRequests['F-Output-1'].setValues(
    variables=('S', 'LE', 'U', 'RF', 'PHILSM', 'STATUSXFEM'))

plateModel.HistoryOutputRequest(createStepName='Static', 
    name='H-Output-2', rebar=EXCLUDE, region=
    mdb.models[modelName].rootAssembly.sets['bdisp'], sectionPoints=
    DEFAULT, variables=('U1', 'RF1'))


#---------------------------------------------------------------------------------

# === Apply boundary conditions ===

plateModel.TabularAmplitude(name='Amp-1', timeSpan=STEP, 
    smooth=SOLVER_DEFAULT, data=((0.0, 0.0), (1, 1.0)))
plateModel.DisplacementBC(amplitude=UNSET, createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='fixz', region=plateModel.rootAssembly.instances['plate_1'].sets['fixZ']
    , u3=0.0)

plateModel.DisplacementBC(  createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='rp',amplitude='Amp-1',  region=plateModel.rootAssembly.sets['bdisp'], u1=0.01)

plateModel.DisplacementBC(amplitude=UNSET, createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='bot', region=plateModel.rootAssembly.instances['plate_1'].sets['bottom']
    , u2=0)

plateModel.DisplacementBC( createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='top',amplitude='Amp-1',  region=plateModel.rootAssembly.instances['plate_1'].sets['top']
    , u1=-0.01, u2=0)


#---------------------------------------------------------------------------------

# === Define enrichment and initial crack ===

faces = plateModel.rootAssembly.instances['crack_1'].faces
f1 = faces.findAt(( (crack_len/2.0, crack_y, width/2.0), ))
crackLocation = Region(faces=f1)

plateModel.ContactProperty('contact')

plateModel.rootAssembly.engineeringFeatures.XFEMCrack( name='enr1',
    crackDomain=plateModel.rootAssembly.instances['plate_1'].sets['All']
    , interactionProperty='contact', crackLocation=crackLocation)

plateModel.interactionProperties['contact'].FractureCriterion(
    initTable=((GI, GI, GI, 1.0, 1.0, 1.0), ), 
    mixedModeBehavior=POWER, tolerance=0.2, viscosity=1e-04)

plateModel.rootAssembly.engineeringFeatures.cracks['enr1'].setValues(interactionProperty='contact')

#---------------------------------------------------------------------------------

# === Create Job === 

mdb.Job(model=modelName, name='modeII_lefm_xfem_c3d8r',
    description='Crack propagation in a plate under modeII loading (XFEM)')


#---------------------------------------------------------------------------------