'''
-----------------------------------------------------------------------------
 Two dimensional symmetric global model of a double edged notch specimen
 (linear elastic material) modeled using hybrid reduced integration plane
 strain elements (CPE8RH). Coarse mesh is used.
-----------------------------------------------------------------------------
'''

from abaqus import *
from abaqusConstants import *

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

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

# Create a model

Mdb()
modelName = '2DDoubleEdPlasGlCPE8RH'
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, 100.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='ElasticMaterial')
myModel.materials['ElasticMaterial'].Elastic(table=((30000000.0, 0.3), ))
myModel.HomogeneousSolidSection(name='SolidHomogeneous',
    material='ElasticMaterial', 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,100,0))
side1Edge1 = edges[side1Edge1.index:(side1Edge1.index+1)]
myAssembly.Surface(side1Edges=side1Edge1, name='topSurf')

# Create the circular partition for sweep meshing around the
# crack tip and the partition to be used while monitoring the
# plastic zone.

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)
mySketch.Line(point1=(10.0,-15.0), point2=(-10.0,-15.0))
mySketch.Line(point1=(0.0, -25.0), point2=(10.0, -15.0))
mySketch.Line(point1=(0.0, -25.0), point2=(-10.0, -15.0))
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.,5.,0.),),
    ((0.,23.52,0.),), ((0.,15.,0.),),)
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, distribution=UNIFORM, localCsys=None)

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

# Assign the loads

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

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

# Create a mesh 

# Seed all the edges

pickedEdges1 = myPlateInstance.edges.findAt(((10.491145,0.491145,0.),),
    ((10.745097,0.,0.),),)
pickedEdges2 = myPlateInstance.edges.findAt(((9.315948,0.,0.),),
    ((9.508855,0.491145,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2, 
    ratio=3.0, number=3, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.,5.,0.),),
    ((14.619398,1.913417,0.),), ((5.380602,1.913417,0.),),
    ((20.,5.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((2.5,0.,0.),),
    ((3.232233,6.767767,0.),), ((16.767767,6.767767,0.),),
    ((17.5,0.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FIXED)

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

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

pickedEdges1 = myPlateInstance.edges.findAt(((20.,28.2,0.),),)
pickedEdges2 = myPlateInstance.edges.findAt(((0,28.2,0),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3.0, number=4, constraint=FIXED)

# Assign meshing controls to the respective regions

pickedRegions = myPlateInstance.faces.findAt(((9.412402,0.245573,0.),),
    ((10.,0.491145,0.),), ((10.618121,0.245573,0.),),)
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED,
    technique=SWEEP)
pickedRegions = myPlateInstance.faces.findAt(((10.,15.,0.),),
    ((2.866117,3.383883,0.),), ((10.,6.767767,0.),),
    ((17.133883,3.383883,0.),), ((10.,27.04,0.),),)
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD,
    technique=STRUCTURED)
elemType1 = mesh.ElemType(elemCode=CPE8RH, 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=4)

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

# Create the job 

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

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