'''
-----------------------------------------------------------------------------
 Nonuniform crack-face loading on an edge crack in a half space modeled
 using reduced integration plane strain elements (CPE8R).
-----------------------------------------------------------------------------
'''

from abaqus import *
from abaqusConstants import *

import part, material, section, assembly, step, interaction
import regionToolset, displayGroupMdbToolset as dgm, mesh, load, job 
import os

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

# Create a model

Mdb()
modelName = '2DEdgeCrackCPE8R'
myModel = mdb.Model(name=modelName)
    
# Create a new viewport in which to display the model
# and the results of the analysis.

myViewport = session.Viewport(name=modelName)
myViewport.makeCurrent()
myViewport.maximize()
    
#---------------------------------------------------------------------------

# Create a part

# Create a sketch for the base feature

mySketch = myModel.Sketch(name='plateProfile',sheetSize=200.0)
mySketch.setPrimaryObject(option=STANDALONE)
mySketch.rectangle(point1=(0.0, 0.0), point2=(20.0, 20.0))

myPlate = myModel.Part(name='plate',
    dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
myPlate.BaseShell(sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

myViewport.setValues(displayedObject=myPlate)

# Create a partition for the sweep mesh

f = myPlate.faces.findAt((10,10,0))
t = myPlate.MakeSketchTransform(sketchPlane=f,
    sketchPlaneSide=SIDE1, origin=(10.0, 10.0, 0.0))
mySketch = myModel.Sketch(name='partitionProfile', 
    sheetSize=40.0, gridSpacing=1.0, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myPlate.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(-10.0,-10.0))
mySketch.ArcByCenterEnds(center=(-9.0,-10.0), point1=(-8.5,-10.0),
    point2=(-9.5,-10.0), direction=COUNTERCLOCKWISE)
pickedFaces = myPlate.faces.findAt((10,10,0))
myPlate.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['partitionProfile']

e1 = myPlate.edges.findAt((1.25,0,0))
e2 = myPlate.edges.findAt((1,0.5,0))
myPlate.PartitionEdgeByPoint(edge=e1,
    point=myPlate.InterestingPoint(edge=e2, rule=CENTER))

# Create a set for the whole part

faces1 = myPlate.faces.findAt(((1,0.25,0),), ((10,10,0),),)
myPlate.Set(faces=faces1, name='all')

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

# Assign material properties

# Create linear elastic material

myModel.Material(name='linearElas')
myModel.materials['linearElas'].Elastic(table=((200000000000.0, 0.3),))

myModel.HomogeneousSolidSection(name='solidHomogeneous',
    material='linearElas', thickness=1.0)

region = myPlate.sets['all']

# Assign the above section to the part

myPlate.SectionAssignment(region=region, sectionName='solidHomogeneous')

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

# Create an assembly

myAssembly = myModel.rootAssembly
myViewport.setValues(displayedObject=myAssembly)
myAssembly.DatumCsysByDefault(CARTESIAN)
myAssembly.Instance(name='myPlate-1', part=myPlate)
myPlateInstance = myAssembly.instances['myPlate-1']

# Create a set for the crack tip

verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((1,0,0))
verts1 = verts[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='crackTip')

# Create a set for the X edge

e1 = myPlateInstance.edges.findAt(((1.25,0,0),), ((12,0,0),),)
myAssembly.Set(edges=e1, name='xEdge')

# Create a set for the vertex to be fixed in X

verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((20,0,0))
verts1 = verts[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='pt')

# Create a set for the load surface

side1Edges1 = myPlateInstance.edges.findAt(((0.25,0,0),),
    ((0.75,0,0),),)
myAssembly.Surface(side1Edges=side1Edges1, name='loadSurf')

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

# Create steps for applying a load
# 1st step: constant load
# 2nd step: linearly varying load
# 3rd step: quadratically varying load
# 4th step: cubically varying load
# 5th step: quartically varying load

myModel.StaticStep(name='Step-1', previous='Initial')
myModel.StaticStep(name='Step-2', previous='Step-1')
myModel.StaticStep(name='Step-3', previous='Step-2')
myModel.StaticStep(name='Step-4', previous='Step-3')
myModel.StaticStep(name='Step-5', previous='Step-4')

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

# Create interaction properties

crackFront = crackTip = myAssembly.sets['crackTip']
v11 = myPlateInstance.vertices.findAt((1,0,0))
v21 = myPlateInstance.vertices.findAt((20,0,0))
myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=ON, crackFront=crackFront, crackTip=crackTip, 
    extensionDirectionMethod=Q_VECTORS, qVectors=((v11, v21),), 
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

# Fix the X edge in the Y direction

region = myAssembly.sets['xEdge']
myModel.DisplacementBC(name='yFixed', createStepName='Initial',
    region=region, u2=SET, distribution=UNIFORM,
    localCsys=None)

# Fix a vertex in the X direction

region = myAssembly.sets['pt']
myModel.DisplacementBC(name='xFixed', createStepName='Initial',
    region=region, u1=SET, distribution=UNIFORM,
    localCsys=None)

# Assign load conditions

region = myAssembly.surfaces['loadSurf']
myModel.Pressure(name='prLoad1', createStepName='Step-1',
    region=region, distribution=USER_DEFINED, 
    magnitude=0.0, amplitude=UNSET)
myModel.loads['prLoad1'].deactivate('Step-2')

myModel.Pressure(name='prLoad2', createStepName='Step-2',
    region=region, distribution=USER_DEFINED, 
    magnitude=0.0, amplitude=UNSET)
myModel.loads['prLoad2'].deactivate('Step-3')

myModel.Pressure(name='prLoad3', createStepName='Step-3',
    region=region, distribution=USER_DEFINED, 
    magnitude=0.0, amplitude=UNSET)
myModel.loads['prLoad3'].deactivate('Step-4')

myModel.Pressure(name='prLoad', createStepName='Step-4',
    region=region, distribution=USER_DEFINED, 
    magnitude=0.0, amplitude=UNSET)
myModel.loads['prLoad'].deactivate('Step-5')

myModel.Pressure(name='prLoad5', createStepName='Step-5',
    region=region, distribution=USER_DEFINED, 
    magnitude=0.0, amplitude=UNSET)

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

# Create a mesh 

# Create a partition 

f1 = myPlateInstance.faces.findAt((10,10,0))
t = myAssembly.MakeSketchTransform(sketchPlane=f1,
    sketchPlaneSide=SIDE1, origin=(10.008844,10.009619,0.0))
mySketch = myModel.Sketch(name='partitionProfile', 
    sheetSize=40.0, gridSpacing=1.0, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(-10.008844,-10.009619))
mySketch.ArcByCenterEnds(center=(-9.008844,-10.009619),
    point1=(-6.008844,-10.009619),
    point2=(-10.5508360682907,-7.21368909124756), 
    direction=COUNTERCLOCKWISE)
pickedFaces = myPlateInstance.faces.findAt(((10,10,0),),
    ((1,0.25,0),),)
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['partitionProfile']

# Seed all the edges

e1 = myPlateInstance.edges
pickedEdges1 = myPlateInstance.edges.findAt((1.031788,0,0),)
pickedEdges2 = myPlateInstance.edges.findAt((0.967369,0,0),)
pickedEdges1 = e1[pickedEdges1.index:(pickedEdges1.index+1)]
pickedEdges2 = e1[pickedEdges2.index:(pickedEdges2.index+1)]
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2, 
    ratio=2.0, number=5, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((1,0.5,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=12, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.25,0,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.,1.414214,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=8, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((2.75,0,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=8, constraint=FIXED)

pickedEdges2 = myPlateInstance.edges.findAt(((0,4,0),),)
myAssembly.seedEdgeByBias(end2Edges=pickedEdges2, ratio=3, number=8,
    constraint=FIXED)

pickedEdges1 = myPlateInstance.edges.findAt(((5,0,0),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=1, number=8,
    constraint=FIXED)

# Assign meshing controls

pickedFaces = myPlateInstance.faces.findAt(((1,0.25,0),),)
myAssembly.setMeshControls(regions=pickedFaces,
    elemShape=QUAD_DOMINATED, technique=SWEEP)

pickedFaces = myPlateInstance.faces.findAt(((2.04152,1.274165,0.),),
    ((6.29152,11.274165,0.),),)
myAssembly.setMeshControls(regions=pickedFaces,
    elemShape=QUAD, technique=STRUCTURED)
elemType1 = mesh.ElemType(elemCode=CPE8R, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPE6M, elemLibrary=STANDARD)
faces1 = myPlateInstance.faces
pickedRegions =(faces1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2))
partInstances =(myPlateInstance, )
myAssembly.generateMesh(regions=partInstances)

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

# Request history output for the crack

del myModel.historyOutputRequests['H-Output-1']
myModel.HistoryOutputRequest(name='JInt1', createStepName='Step-1',
    contourIntegral='Crack', numberOfContours=4)
myModel.historyOutputRequests['JInt1'].deactivate('Step-2')
myModel.HistoryOutputRequest(name='JInt2', createStepName='Step-2',
    contourIntegral='Crack', numberOfContours=4)
myModel.historyOutputRequests['JInt2'].deactivate('Step-3')
myModel.HistoryOutputRequest(name='JInt3', createStepName='Step-3',
    contourIntegral='Crack', numberOfContours=4)
myModel.historyOutputRequests['JInt3'].deactivate('Step-4')
myModel.HistoryOutputRequest(name='JInt4', createStepName='Step-4',
    contourIntegral='Crack', numberOfContours=4)
myModel.historyOutputRequests['JInt4'].deactivate('Step-5')
myModel.HistoryOutputRequest(name='JInt5', createStepName='Step-5',
    contourIntegral='Crack', numberOfContours=4)

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

# Create the job 

# Change the name of the user subroutine file based on platform
# If a NTI platform is used the fortran file with *.for extension
# is used, whereas on all other platforms a *.f extension is
# necessary.

if os.environ['PLATFORM'] == 'nti':
    userFile = '2DEdgeCrackCPE8R.for'
else:
    userFile = '2DEdgeCrackCPE8R.f'

myJob = mdb.Job(name=modelName, model=modelName,
    description='Contour integral analysis',
    userSubroutine=userFile)
mdb.saveAs(pathName=modelName)

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