'''
-----------------------------------------------------------------------------
 Two dimensional symmetric model of a double edged notch specimen (linear
 elastic material) modeled using reduced integration plane strain
 elements (CPE8R).
-----------------------------------------------------------------------------
'''

from abaqus import *
import testUtils
testUtils.setBackwardCompatibility()
from abaqusConstants import *

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

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

# Create a model

Mdb()
modelName = '2DDoubleEdSymmCPE8R'
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.sketchOptions.setValues(viewStyle=REGULAR)
mySketch.setPrimaryObject(option=STANDALONE)

mySketch.rectangle(point1=(0.0, 0.0), point2=(20.0, 50.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)

# Partition the edge on the symmetry plane

e1 = myPlate.edges.findAt((10,0,0))
(myPlate.PartitionEdgeByPoint(edge=e1,
    point=myPlate.InterestingPoint(edge=e1, rule=MIDDLE)))

# Create a set referring to the whole part

faces = myPlate.faces.findAt(((10,25,0),))
myPlate.Set(faces=faces, name='All')

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

# Assign material properties

# Create linear elastic material

myModel.Material(name='LinearElastic')
myModel.materials['LinearElastic'].Elastic(table=((30000000.0, 0.3), ))
myModel.HomogeneousSolidSection(name='SolidHomogeneous',
    material='LinearElastic', 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, dependent=OFF)
myPlateInstance = myAssembly.instances['myPlate-1']

# Create a set for the top surface of the plate

edges = myPlateInstance.edges
side1Edge1 = myPlateInstance.edges.findAt((10,50,0))
side1Edge1 = edges[side1Edge1.index:(side1Edge1.index+1)]
myAssembly.Surface(side1Edges=side1Edge1, name='topSurf')

# Create circular partitions for sweep mesh around the crack tip

face1 = myPlateInstance.faces.findAt((10,25,0),)
t = myAssembly.MakeSketchTransform(sketchPlane=face1,
    sketchPlaneSide=SIDE1, origin=(10.0, 25.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', sheetSize=107.7,
    gridSpacing=2.69, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
r, r1 = mySketch.referenceGeometry, mySketch.referenceVertices
mySketch.sketchOptions.setValues(gridOrigin=(0.0,-25.0))
mySketch.ArcByCenterEnds(center=(0.0, -25.0), point1=(5.0,-25.0),
    point2=(-5.0,-25.0), direction=COUNTERCLOCKWISE)
f1 = myPlateInstance.faces
pickedFaces = f1[face1.index:(face1.index+1)]
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a set for the X edge of the plate

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((2.5,0,0))
e2 = myPlateInstance.edges.findAt((7.5,0,0))
edge = edges[e1.index:(e1.index+1)] + edges[e2.index:(e2.index+1)]
myAssembly.Set(edges=edge, name='xEdge')

# Create a set for the Y edge of the plate

edges = myPlateInstance.edges
edge1 = myPlateInstance.edges.findAt((0,25,0))
edge1 = edges[edge1.index:(edge1.index+1)]
myAssembly.Set(edges=edge1, name='yEdge')

# Create a set for the crack tip

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

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

# Create a step for applying a load

myModel.StaticStep(name='ApplyLoad', previous='Initial',
    description='Apply a pressure load')

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

# Create interaction properties

# Create the contour integral definition for the crack

crackFront = crackTip = myAssembly.sets['crackTip']
verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((10,0,0))
v2 = myPlateInstance.vertices.findAt((0,0,0))
myAssembly.engineeringFeatures.ContourIntegral(name='Crack', symmetric=ON,
    crackFront=crackFront, crackTip=crackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((v1,v2),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

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

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

# Assign the loads

region = myAssembly.surfaces['topSurf']
myModel.Pressure(name='topLoad', createStepName='ApplyLoad',
    region=region, distributionType=UNIFORM, magnitude=-100.0)

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

# Create a mesh 

# Seed all the edges

e1 = myPlateInstance.edges
pickedEdges1 = myPlateInstance.edges.findAt((11,0,0),)
pickedEdges2 = myPlateInstance.edges.findAt((9,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=6, constraint=FIXED)

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

pickedEdges = myPlateInstance.edges.findAt(((0,25,0),), ((20,25,0),))
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=15, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((10,50,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=5, constraint=FIXED)

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

pickedEdges = myPlateInstance.edges.findAt(((2.5,0,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=3, constraint=FIXED)

# Assign meshing controls to the respective regions

faces = myPlateInstance.faces
f1 = myPlateInstance.faces.findAt((10,2.5,0))
pickedRegions = faces[f1.index:(f1.index+1)]
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED,
    technique=SWEEP)
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

myModel.historyOutputRequests.changeKey(fromName='H-Output-1',
    toName='JInt')
myModel.historyOutputRequests['JInt'].setValues(contourIntegral='Crack',
    numberOfContours=7)

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

# Create the job 

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

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