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

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.05                          # crack offset in y-direction

YM    = 2.10e+11                          # Young's modulus
NU    = 0.3                               # Poisson's ratio
MAXPS = 2.2e+08                          # 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=TWO_D_PLANAR, name='plate', 
                type=DEFORMABLE_BODY)
platePart.BaseShell(sketch=plateModel.sketches['plateProfile'])

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=TWO_D_PLANAR, name='crack',
                type=DEFORMABLE_BODY)
crackPart.BaseWire(sketch=plateModel.sketches['crackProfile'])

crackSketch.unsetPrimaryObject()
del crackSketch


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


# === Create geometry sets ===

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

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=POWER_LAW, power=ETA, table=((GI, GI, GI), ), type=ENERGY)
plateMatl.maxpsDamageInitiation.DamageStabilizationCohesive(
    cohesiveCoeff=0.001)

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 plate ===

# 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 = 30;  ey = 60;

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

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

platePart.generateMesh()

# === End part ===


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

# === 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, -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(initialInc=0.005, maxInc=0.01, maxNumInc=10000,
    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['bdisp'], sectionPoints=
    DEFAULT, variables=('U1', 'RF1'))


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

# === Apply boundary conditions ===

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

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

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


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

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

edges = plateModel.rootAssembly.instances['crack_1'].edges
e1 = edges.findAt(( (crack_len/2.0, crack_y, 0.0), ))
crackLocation = Region(edges=e1)

plateModel.ContactProperty('contact')

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


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

# === Create Job === 

mdb.Job(model=modelName, name='mixed_mode_xfem_cpe4',
    description='Crack propagation in a plate under mixed-mode loading (XFEM)')


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