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

length = 0.2286                           # plate half-length
height = .0726                            # plate height
width  = .0254                            # plate width

crack_len = .019                          # crack length


YM    = 3.137e+10                         # Young's modulus
NU    = 0.2                               # Poisson's ratio
MAXPS = 1.045e+07                         # Damage initiation
DTOL  = 0.05                              # Damage tolerance
GI    = 19.58                             # Fracture energy
ETA   = 1.0                               # Power-law exponent
density=2400                              # Density

l=length/2-.0127                          # Distance from center to support 
rate=0.635                                # l1/l, l1 is the distance from center to crack
Crack_X=length/2-l*rate                   # X coordinate of the crack

delta=.0001

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

# === 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, 0), point2=(0.0, height))
plateSketch.Line(point1=(0.0, height), point2=(length,height ))
plateSketch.Line(point1=(length, height), point2=(length, 0))
plateSketch.Line(point1=(length, 0), point2=(0.0, 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=(Crack_X,0), point2=(Crack_X,crack_len))

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.Density(table=((density, ), 
    ))
plateModel.HomogeneousSolidSection(material='elas', name='solid')

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

# === Partition plate ===

plate_part = mdb.models[modelName].parts['plate']
faces=plate_part.faces
plate_part.DatumPlaneByOffset(plane=faces[1], flip=SIDE2, offset=height-crack_len)
pickedCells= platePart.cells.findAt((0.0,0,0.0))
datums1 = plate_part.datums
plate_part.PartitionCellByDatumPlane(datumPlane=datums1[2], 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[modelName].rootAssembly
 
plate_assem.Instance(name='plate-1', part=plate_part, dependent=OFF) 
part_add = mdb.models[modelName].parts['crack']
plate_assem.Instance(name='crack-1', part=part_add, dependent=ON)

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

# === Create Seed & Mesh ===

pickedEdges=plate_assem.instances['plate-1'].edges.getByBoundingBox(0,0,0,length ,height ,width  )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=90, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt(
  ((0.0,0,width/2  ),), ((0.0,crack_len,width/2  ),), ((0.0,height,width/2  ),), 
  ((length ,0,width/2  ),), ((length ,crack_len,width/2  ),), ((length ,height,width/2  ),), 
   )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=10, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt(
  ((0.0,crack_len/2,0),), ((length ,crack_len/2,0),),  
  ((0.0,crack_len/2,width  ),), ((length ,crack_len/2,width  ),),
   )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=7, constraint=FIXED)

pickedEdges=plate_assem.instances['plate-1'].edges.findAt(
  ((0.0,height/2,0),), ((length ,height/2,0),),  
  ((0.0,height/2,width  ),), ((length ,height/2,width  ),),
   )
plate_assem.seedEdgeByNumber(edges=pickedEdges, number=23, 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 ===

pickedNodes=plate_assem.instances['plate-1'].nodes.getByBoundingBox(1.27e-02-delta,  0,  0, 1.27e-02+delta, 0 ,  width)
plate_assem.Set(nodes=pickedNodes, name='bottom_left')

pickedNodes=plate_assem.instances['plate-1'].nodes.getByBoundingBox(.2159-delta,  0,   0, .2159+delta,  0,  width)
plate_assem.Set(nodes=pickedNodes, name='bottom_right')

pickedNodes=plate_assem.instances['plate-1'].nodes.getByBoundingBox( .2159-delta,  -delta,  width-delta,.2159+delta,  delta,  width+delta)

plate_assem.Set(nodes=pickedNodes, name='fix_x_left')
pickedNodes=plate_assem.instances['plate-1'].nodes.getByBoundingBox(  1.27e-02-delta,  0-delta,  width-delta,1.27e-02+delta,  delta,  width+delta)

plate_assem.Set(nodes=pickedNodes, name='fix_x_right')
pickedNodes=plate_assem.instances['plate-1'].nodes.getByBoundingBox(length/2-delta,  height-delta, -delta,length/2+delta,  height+delta,  width+delta)
plate_assem.Set(nodes=pickedNodes, name='top')

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=1.3e-3, 
    maxNumInc=10000, initialInc=5e-07, minInc=1.5e-15, maxInc=5e-06, 
    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'))
plateModel .steps['Step-1'].control.setValues(
    allowPropagation=OFF, resetDefaultValues=OFF, timeIncrementation=(4.0, 8.0, 
    9.0, 16.0, 10.0, 4.0, 12.0, 10.0, 6.0, 3.0, 50.0))
plateModel.steps['Step-1'].control.setValues(discontinuous=ON)

regionDef=plateModel.rootAssembly.sets['top']
plateModel.HistoryOutputRequest(name='H-Output-2', 
    createStepName='Step-1', variables=('RF2', ), region=regionDef, 
    sectionPoints=DEFAULT, rebar=EXCLUDE)
#---------------------------------------------------------------------------------

# === Apply Boundary Conditions === 

plateModel.TabularAmplitude(name='Amp-1', timeSpan=STEP, 
    smooth=SOLVER_DEFAULT, data=((0.0, 0.0), (1.96e-4, 1.0), (1.5e-3, 1.0)))
region = plate_assem.sets['top']
plateModel.VelocityBC(name='velocity', createStepName='Step-1', 
    region=region, v1=UNSET, v2=-0.06, v3=UNSET, vr1=UNSET, vr2=UNSET, 
    vr3=UNSET, amplitude='Amp-1', localCsys=None, distributionType=UNIFORM, 
    fieldName='')
region = plate_assem.sets['bottom_left'] 
plateModel.DisplacementBC(name='BC-1', createStepName='Step-1', 
    region=region, u1=UNSET, u2=0.0, u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, 
    amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
    localCsys=None)
region = plate_assem.sets['bottom_right'] 
plateModel.DisplacementBC(name='BC-2', createStepName='Step-1', 
    region=region, u1=0.0, u2=0.0, u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, 
    amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
    localCsys=None)

region = plate_assem.sets['fix_x_left']
plateModel.DisplacementBC(name='BC-3', createStepName='Step-1', 
    region=region, u1=UNSET, u2=UNSET, u3=0.0, ur1=UNSET, ur2=UNSET, ur3=UNSET, 
    amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
    localCsys=None)

region = plate_assem.sets['fix_x_right'] 
plateModel.DisplacementBC(name='BC-4', createStepName='Step-1', 
    region=region, u1=UNSET, u2=UNSET, u3=0.0, ur1=UNSET, ur2=UNSET, ur3=UNSET, 
    amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
    localCsys=None)

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

# === Create Job === 

mdb.Job(name='crackprop_tpb_xfem_c3d8r_dyn', model='tpb_xfem_dyn', 
    description='Crack propagation in a beam under impact loading simulated using XFEM' )