# === Modules ===

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


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

# === Parameters ===

modelName = 'Crack_propagation'

length = 100.0e-03                        # plate half-length
radius = length/5.0                       # hole radius
height = 3.0*length + 2.0*radius          # plate height
width  = 20.0e-03                         # plate width

YM    = 3.24e+09                          # Young's modulus
NU    = 0.3                               # Poisson's ratio
MAXPS = 22.0e+06                          # Damage initiation
DTOL  = 0.01                              # Damage tolerance
GI    = 2870.0                            # Fracture energy
ETA   = 1.0                               # B-K exponent


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

# === Create model ===

Mdb()

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

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


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

# === Create part ===

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

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

plateSketch.ConstructionLine(angle=90.0, point1=(0.0, 0.0))
plateSketch.ArcByCenterEnds(center=(0.0, 0.0), direction=COUNTERCLOCKWISE, 
    point1=(0.0, -radius), point2=(0.0, +radius))
plateSketch.Line(point1=(0.0, radius), 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))
plateSketch.Line(point1=(0.0, -height/2.0), point2=(0.0, -radius))

platePart = plateModel.Part(dimensionality=TWO_D_PLANAR, name='plate', 
                type=DEFORMABLE_BODY)
platePart.BaseShell(sketch=plateModel.sketches['plateProfile'])

plateSketch.unsetPrimaryObject()
del plateSketch


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

# === Partition part for meshing ===

f1 = platePart.faces.findAt((length/2.0, 0.0, 0.0))
t1 = platePart.MakeSketchTransform(sketchPlane=f1, origin=(0.0, 0.0, 0.0))

plateSketch = plateModel.ConstrainedSketch(gridSpacing=0.01, name='partitionProfile', 
                 sheetSize=height, transform=t1)
plateSketch.setPrimaryObject(option=SUPERIMPOSE)
plateSketch.sketchOptions.setValues(decimalPlaces=3)

platePart.projectReferencesOntoSketch(filter=COPLANAR_EDGES, sketch=plateSketch)

plateSketch.Line(point1=(cos(pi/6)*radius, -radius/2.0), point2=(length, -radius/2.0))
plateSketch.Line(point1=(cos(pi/6)*radius, +radius/2.0), point2=(length, +radius/2.0))
plateSketch.Line(point1=(0.0, radius), point2=(length, length/2.0+radius))
plateSketch.Line(point1=(0.0, -radius), point2=(length, -length/2.0-radius))

platePart.PartitionFaceBySketch(faces=f1, sketch=plateSketch)

plateSketch.unsetPrimaryObject()
del plateSketch


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

# === Create geometry sets ===

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

e1  = platePart.edges.findAt(( (0.0,  radius+length, 0.0), ))
e1 += platePart.edges.findAt(( (0.0, -radius-length, 0.0), ))
platePart.Set(edges=e1, name='xsymm')

e1 = platePart.edges.findAt(( (length/2.0, -height/2.0, 0.0), ))
platePart.Set(edges=e1, name='bottom')

e1 = platePart.edges.findAt(( (length/2.0, +height/2.0, 0.0), ))
platePart.Set(edges=e1, name='top')


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

# === 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=1e-05)

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


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

# === 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 part ===

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

# Mesh technique
platePart.setMeshControls(elemShape=QUAD, regions=platePart.faces[:], technique=STRUCTURED)


# Seed mesh
ex = 20; ey_a = 40;  ey_b = 6; ey_c = 11;

e1 = platePart.edges.findAt(( (length/2.0, -height/2.0, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ex)

e1 = platePart.edges.findAt(( (length, -length-radius, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey_a)

e1 = platePart.edges.findAt(( (length, -length/2.0+radius, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey_b)

e1 = platePart.edges.findAt(( (length, 0.0, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey_c)

e1 = platePart.edges.findAt(( (length, +length/2.0-radius, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey_b)

e1 = platePart.edges.findAt(( (length, +length+radius, 0.0), ))
platePart.seedEdgeByNumber(edges=e1, number=ey_a)

platePart.generateMesh()

# === End part ===


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

# === Assemble ===

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

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


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

# === Create assembly sets ===

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

v1 = (plateModel.rootAssembly.referencePoints[rp_dt.id], )
plateModel.rootAssembly.Set(name='tdisp', referencePoints=v1)

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

# === End assembly ===


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

# === Create constraint equations ===

plateModel.Equation(name='ce_bot', terms=((1.0, 'plate_1.bottom', 2), (-1.0, 'bdisp', 2)))
plateModel.Equation(name='ce_top', terms=((1.0, 'plate_1.top', 2), (-1.0, 'tdisp', 2)))


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

# === Create step ===

plateModel.StaticStep(initialInc=0.01, maxInc=0.01, maxNumInc=100000,
    minInc=1e-09, name='Static', nlgeom=ON, previous='Initial')

plateModel.steps['Static'].control.setValues(
    allowPropagation=OFF, discontinuous=ON, resetDefaultValues=OFF, 
    timeIncrementation=(8.0, 10.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['rdisp'], sectionPoints=
    DEFAULT, variables=('U2', 'RF2'))


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

# === Apply boundary conditions ===

plateModel.DisplacementBC(amplitude=UNSET, createStepName=
    'Initial', distributionType=UNIFORM, fieldName='', localCsys=None, name=
    'symm', region=plateModel.rootAssembly.instances['plate_1'].sets['xsymm']
    , u1=SET, u2=UNSET, ur3=UNSET)

plateModel.DisplacementBC(amplitude=UNSET, createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='bot', region=plateModel.rootAssembly.sets['bdisp'], u1=UNSET, u2=
    -0.00055, ur3=UNSET)

plateModel.DisplacementBC(amplitude=UNSET, createStepName=
    'Static', distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None
    , name='top', region=plateModel.rootAssembly.sets['tdisp'], u1=UNSET, u2=
    0.00055, ur3=UNSET)


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

# === Define enrichment ===

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


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

# === Create Job === 

mdb.Job(model=modelName, name='crack_xfem',
    description='Crack initiation and propagation in a plate with a hole (XFEM)')


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