# === 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 = 'shear_xfem_3d_dyn'

length = 3.0                              # plate half-length
height = 6.0                              # plate height
width  = .1                              # plate width

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

YM    = 3240                         # Young's modulus
NU    = 0.35                               # Poisson's ratio
MAXPS = 100                          # Damage initiation
DTOL  = 0.05                              # Damage tolerance
GI    = .7                           # 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=width  )

plateSketch.unsetPrimaryObject()
del plateSketch


# Part for crack geometry

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

crackSketch.Line(point1=(0.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=width  )

crackSketch.unsetPrimaryObject()
del crackSketch

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

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

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

plateMatl.Elastic(table=((YM, NU), ))
plateMatl.MaxpsDamageInitiation(table=((MAXPS, ), ), tolerance=DTOL)
plateMatl.maxpsDamageInitiation.DamageEvolution(
    mixedModeBehavior=BK, power=ETA, table=((GI, GI, GI), ), type=ENERGY)
plateMatl.maxpsDamageInitiation.DamageStabilizationCohesive(
    cohesiveCoeff=5E-7)
plateMatl.Density(table=((1.19e-09, ), 
    ))
plateModel.HomogeneousSolidSection(material='elas', name='solid', thickness=width)

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

# ===  Partition plate  ===

plate_part = mdb.models['shear_xfem_3d_dyn'].parts['plate']
faces=plate_part.faces
plate_part.DatumPlaneByOffset(plane=faces[1], flip=SIDE2, offset=height/2+.1)
plate_part.DatumPlaneByOffset(plane=faces[1], flip=SIDE2, offset=height/2-.1)
plate_part.DatumPlaneByOffset(plane=faces[2], flip=SIDE2, offset=length/2+.1 )
plate_part.DatumPlaneByOffset(plane=faces[2], flip=SIDE2, offset=length/2-.1 )


pickedCells= platePart.cells.findAt((0.0,-height/2,0.0))
datums1 = plate_part.datums
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[2], cells=pickedCells)
pickedCells= platePart.cells.findAt((0.0,height/2,0.0))
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[3], cells=pickedCells)
pickedCells= platePart.cells.findAt(((0.0,height/2,0.0),),((0.0,0.0,0.0),),((0.0,-height/2,0.0),))
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[4], cells=pickedCells)
pickedCells= platePart.cells.findAt(((length,height/2,0.0),),((length,0.0,0.0),),((length,-height/2,0.0),))
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[5], cells=pickedCells)

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

# === Assign section ===

platePart.Set(cells=platePart.cells[:], name='All')
platePart.SectionAssignment(region=platePart.sets['All'], sectionName='solid')

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

# === Create Assembly ===

plate_assem = mdb.models['shear_xfem_3d_dyn'].rootAssembly
part_add = mdb.models['shear_xfem_3d_dyn'].parts['plate']
plate_assem.Instance(name='plate-1', part=part_add, dependent=OFF)
part_add = mdb.models['shear_xfem_3d_dyn'].parts['crack']
plate_assem.Instance(name='crack-1', part=part_add, dependent=ON)

plate_assem.translate(instanceList=('crack-1', ), vector=(0.0, -crack_y, 0.0))

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

# === Create Seed & Mesh ===

pickedEdges=plate_assem.instances['plate-1'].edges.findAt(((0.0,height/4,0.0),),
    ((0.0,-height/4,0.0),),((length ,height/4,0.0),),((length ,-height/4,0.0),),
    ((length/2+.1 ,-height/4,0.0),),((length/2-.1 ,-height/4,0.0),), 
    ((length/2+.1 ,height/4,0.0),),((length/2-.1 ,height/4,0.0),),   
    ((0.0,height/4,width  ),),
    ((0.0,-height/4,width  ),),((length ,height/4,width  ),),((length ,-height/4,width  ),),
    ((length/2+.1 ,-height/4,width  ),),((length/2-.1 ,-height/4,width  ),),
    ((length/2+.1 ,height/4,width  ),),((length/2-.1 ,height/4,width  ),),)    
   

plate_assem.seedEdgeByNumber(edges=pickedEdges, number=40, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt( 
    ((length/4,height/2,0.0),), ((length/4,-height/2,0.0),),((length*3/4,height/2,0.0),), ((length*3/4,-height/2,0.0),), 
    ((length/4,.1,0.0),), ((length/4,-.1,0.0),),((length*3/4,.1,0.0),), ((length*3/4,-.1,0.0),),  
    ((length/4,height/2,0.1),), ((length/4,-height/2,0.1),),((length*3/4,height/2,0.1),), ((length*3/4,-height/2,0.1),), 
    ((length/4,.1,0.1),), ((length/4,-.1,0.1),),((length*3/4,.1,0.1),), ((length*3/4,-.1,0.1),), )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=20, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt( 
      ((0,0,0),), ((length/2+.1 ,0,0),),((length/2-.1 ,0,0),),((length ,0,0),),
      ((0,0,0.1),), ((length/2+.1 ,0,0.1),),((length/2-.1 ,0,0.1),),((length ,0,0.1),),
      ((length/2,height/2,0),), ((length/2,-height/2,0),), ((length/2,.1,0),), ((length/2,-.1,0),), 
      ((length/2,height/2,0.1),), ((length/2,-height/2,0.1),), ((length/2,.1,0.1),), ((length/2,-.1,0.1),),
      )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=15, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt( 
      ((0,0.1,0.05),), ((0 ,-.1,.05),),((0,height/2,.05),),((0 ,-height/2,.05),),  
      ((length/2-.1,0.1,0.05),), ((length/2-.1 ,-.1,.05),),((length/2-.1,height/2,.05),),((length/2-.1 ,-height/2,.05),),   
      ((length/2+.1,0.1,0.05),), ((length/2+.1 ,-.1,.05),),((length/2+.1,height/2,.05),),((length/2+.1 ,-height/2,.05),),   
      ((length,0.1,0.05),), ((length,-.1,.05),),((length,height/2,.05),),((length ,-height/2,.05),),   
      )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FIXED)

elemType1 = ElemType(elemCode=C3D8R, elemLibrary=STANDARD, 
    kinematicSplit=AVERAGE_STRAIN, secondOrderAccuracy=OFF, 
    hourglassControl=DEFAULT, distortionControl=DEFAULT)

elemType2 = ElemType(elemCode=C3D6, elemLibrary=STANDARD, 
    distortionControl=DEFAULT)
elemType3 = ElemType(elemCode=C3D4, elemLibrary=STANDARD, 
    distortionControl=DEFAULT)
 
myplateInstance=plate_assem.instances['plate-1']
cells1 = myplateInstance.cells
pickedRegions =(cells1, )  
plate_assem.setElementType(regions=pickedRegions , elemTypes=(elemType1, elemType2, elemType3))
 
partInstances =(plate_assem.instances['plate-1'], )
plate_assem.generateMesh(regions=partInstances)

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

#=== Create Sets ===

pickedFaces=plate_assem.instances['plate-1'].faces.findAt(
    ((length/2,0,0.1),), ((length/2,height/4,0.1),),  ((length/2,-height/4,0.1),),
    ((length/4,0,0.1),), ((length/4,height/4,0.1),),  ((length/4,-height/4,0.1),),
    ((length*3/4,0,0.1),), ((length*3/4,height/4,0.1),),  ((length*3/4,-height/4,0.1),),
    )
plate_assem.Set(faces=pickedFaces, name='zsym')

pickedFaces=plate_assem.instances['plate-1'].faces.findAt(
    ((length/4,height/2,0.05),), ((length/2,height/2,0.01),),  ((length*3/4,height/2,0.01),),
     )
plate_assem.Set(faces=pickedFaces, name='top')

pickedFaces=plate_assem.instances['plate-1'].faces.findAt(
    ((0,-height/4,0.05),),  
     )
plate_assem.Set(faces=pickedFaces, name='velocity-BC')

pickNodesfrom = plate_assem.instances['plate-1'].nodes 
pickedNodes=pickNodesfrom.getByBoundingBox(-0.01,-.099,-.1,.01,-.001,.2)
#pickedNodes=pickNodesfrom.getByBoundingBox(-.01,-.09,-.1,.01,-,01,.2)
plate_assem.Set(nodes=pickedNodes, name='velocity-BC-nodes')

pickedFaces=plate_assem.instances['crack-1'].faces
plate_assem.Set(faces=pickedFaces,name='XFEM-Crack')

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

#=== Create Enrichment Region and XFEM Crack ===

plateModel.ContactProperty('XFEM-Crack')

crackDomain = plate_assem.instances['plate-1'].sets['All']

crackLocation = plate_assem.sets['XFEM-Crack']

plate_assem.engineeringFeatures.XFEMCrack(name='Crack-1', crackDomain=crackDomain, 
    interactionProperty='XFEM-Crack', crackLocation=crackLocation)

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

# === Create step ===

plateModel.ImplicitDynamicsStep(name='Step-1', 
    previous='Initial', description='Dynamic Shear', timePeriod=6e-06, 
    maxNumInc=10000, initialInc=1e-09, minInc=6e-11, maxInc=5e-09, 
    hafTolMethod=VALUE, haftol=500.0, nlgeom=ON)

del plateModel.fieldOutputRequests['F-Output-1']

plateModel.FieldOutputRequest(name='F-Output-2', 
    createStepName='Step-1', variables=('S', 'E', 'U', 'PHILSM', 'STATUSXFEM'))

# === Apply Boundary Conditions === 

plateModel.TabularAmplitude(name='Amp-1', timeSpan=STEP, 
    smooth=SOLVER_DEFAULT, data=((0.0, 0.0), (1e-07, 1.0), (6e-06, 1.0)))

 
plateModel.ZsymmBC(name='Symmetric', createStepName='Step-1', 
    region=plate_assem.sets['zsym']  )

plateModel.EncastreBC(name='Fixed Top', 
    createStepName='Step-1', region=plate_assem.sets['top']  )

 
plateModel.VelocityBC(name='Velocity', 
    createStepName='Step-1', region=plate_assem.sets['velocity-BC']  , v1=25000.0, v2=UNSET, v3=UNSET, 
    vr1=UNSET, vr2=UNSET, vr3=UNSET, amplitude='Amp-1', localCsys=None, 
    distributionType=UNIFORM, fieldName='')

plateModel.VelocityBC(name='Velocity_at_node', 
    createStepName='Step-1', region=plate_assem.sets['velocity-BC-nodes']  , v1=25000.0, v2=UNSET, v3=UNSET, 
    vr1=UNSET, vr2=UNSET, vr3=UNSET, amplitude='Amp-1', localCsys=None, 
    distributionType=UNIFORM, fieldName='')

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

# === Create Job === 

mdb.Job(name='Crackprop_shear_xfem_3D_dyn', model='shear_xfem_3d_dyn', 
    description='Dynamic crack initiation and propagation analysis with shear')

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