'''
-----------------------------------------------------------------------------
 Full model of a two dimensional double edged notch specimen
 for fracture analysis modeled using plane strain elements with reduced
 integration (CPE8).
-----------------------------------------------------------------------------
'''

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

# use the old default mesh controls.
from mesh import *
session.defaultMesherOptions.setValues(elemShape2D=QUAD)
session.defaultMesherOptions.setValues(quadAlgorithm=MEDIAL_AXIS)

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

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

# Create a model

Mdb()
modelName = '2DDoubleEdgedNotchCPE8R'
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=AXISYM)
mySketch.setPrimaryObject(option=STANDALONE)

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

# Create a partition for the seam cracks and the boundary conditions

face = myPlate.faces.findAt((0,40,0),)
t = myPlate.MakeSketchTransform(sketchPlane=face,
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', 
    sheetSize=215.4, gridSpacing=5.38, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myPlate.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
mySketch.Line(point1=(-20.0, 0.0), point2=(20.0, 0.0))
pickedFaces = myPlate.faces.findAt((0,20,0),)
myPlate.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Partition the middle edge at the midpoint

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

# Create a set referring to the whole part

faces = myPlate.faces.findAt(((0,20,0),), ((0,-20,0),))
myPlate.Set(faces=faces, name='All')

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

# Assign material properties

import material
import section

# 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)
myPlateInstance = myAssembly.instances['myPlate-1']

# Create a set for the mid point of the plate

verts1 = myPlateInstance.vertices
vmp = myPlateInstance.vertices.findAt((0,0,0))
midPt = verts1[vmp.index:(vmp.index+1)]
myAssembly.Set(vertices=midPt, name='midPoint')

# Create circular partitions for sweep mesh around the crack tip

faces1 = myPlateInstance.faces.findAt((0,0,0),)
t = myAssembly.MakeSketchTransform(sketchPlane=faces1,
    sketchPlaneSide=SIDE1, origin=(0.0, -25.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', sheetSize=128.06,
    gridSpacing=3.2, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(0.0, 25.0))
mySketch.CircleByCenterPerimeter(center=(-10.0, 25.0),
    point1=(-15.0, 25.0))
mySketch.CircleByCenterPerimeter(center=(10.0, 25.0),
    point1=(15.0, 25.0))
faces1 = myPlateInstance.faces
f1 = myPlateInstance.faces.findAt((0,30,0),)
f2 = myPlateInstance.faces.findAt((0,-30,0),)
pickedFaces = (faces1[f1.index:(f1.index+1)] +
    faces1[f2.index:(f2.index+1)])
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a set for the top surface

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

# Create a set for the bottom surface

s1 = myPlateInstance.edges.findAt((0,-50,0),)
side1Edges1 = edges1[s1.index:(s1.index+1)]
myAssembly.Surface(side1Edges=side1Edges1, name='bottomSurf')

# Create a set for the edge to be fixed in the Y direction

e1 = myPlateInstance.edges.findAt((-2.5,0,0),)
e2 = myPlateInstance.edges.findAt((2.5,0,0),)
fixedEdges = edges1[e1.index:(e1.index+1)] + edges1[e2.index:(e2.index+1)]
myAssembly.Set(edges=fixedEdges, name='YfixedEdge')

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

# Create a step for applying a load

myModel.StaticStep(name='LoadPlate', previous='Initial',
    description='Apply the load')

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

# Create interaction properties

# Partition the edges to assign the create the left and right seam cracks

leftEdge = myPlateInstance.edges.findAt((-10,0,0),)
myAssembly.PartitionEdgeByPoint(edge=leftEdge,
    point=myPlateInstance.InterestingPoint(edge=leftEdge, rule=MIDDLE))
rightEdge = myPlateInstance.edges.findAt((10,0,0),)
myAssembly.PartitionEdgeByPoint(edge=rightEdge,
    point=myPlateInstance.InterestingPoint(edge=rightEdge, rule=MIDDLE))

# Assign seam crack property to the left crack

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((-17.5,0,0),)
e2 = myPlateInstance.edges.findAt((-12.5,0,0),)
pickedEdges = edges[e1.index:(e1.index+1)] + edges[e2.index:(e2.index+1)]
pickedRegions = regionToolset.Region(edges=pickedEdges)
myAssembly.engineeringFeatures.assignSeam(regions=pickedRegions)

# Assign seam crack property to the right crack

e1 = myPlateInstance.edges.findAt((17.5,0,0),)
e2 = myPlateInstance.edges.findAt((12.5,0,0),)
pickedEdges = edges[e1.index:(e1.index+1)] + edges[e2.index:(e2.index+1)]
pickedRegions = regionToolset.Region(edges=pickedEdges)
myAssembly.engineeringFeatures.assignSeam(regions=pickedRegions)

# Create a set for the left crack tip

verts1 = myPlateInstance.vertices
vxx = myPlateInstance.vertices.findAt((-10,0,0),)
vyy = myPlateInstance.vertices.findAt((20,0,0),)
crackTip = verts1[vxx.index:(vxx.index+1)]
leftCrackTip = leftCrackFront = myAssembly.Set(vertices=crackTip,
    name='leftPt')
myAssembly.engineeringFeatures.ContourIntegral(name='Crack-left',
    symmetric=OFF, crackFront=leftCrackFront, crackTip=leftCrackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((vxx, vyy),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

# Create a set for the right crack tip

verts1 = myPlateInstance.vertices
vxx = myPlateInstance.vertices.findAt((10,0,0),)
vyy = myPlateInstance.vertices.findAt((-20,0,0),)
crackTip = verts1[vxx.index:(vxx.index+1)]
rightCrackTip = rightCrackFront = myAssembly.Set(vertices=crackTip,
    name='rightPt')
myAssembly.engineeringFeatures.ContourIntegral(name='Crack-right',
    symmetric=OFF, crackFront=rightCrackFront, crackTip=rightCrackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((vxx, vyy),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

YEdge = myAssembly.sets['YfixedEdge']
myModel.DisplacementBC(name='YFixed', createStepName='LoadPlate',
    region=YEdge, u2=0.0, fixed=OFF, distributionType=UNIFORM,
    localCsys=None)

mPoint = myAssembly.sets['midPoint']
myModel.DisplacementBC(name='XFixed', createStepName='LoadPlate',
    region=mPoint, u1=0.0, amplitude=UNSET, fixed=OFF,
    distributionType=UNIFORM, localCsys=None)

# Assign load conditions

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

bSurf = myAssembly.surfaces['bottomSurf']
myModel.Pressure(name='BottomLoad', createStepName='LoadPlate',
    region=bSurf, distributionType=UNIFORM, magnitude=-100.0)

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

# Create a mesh 

# Seed all the edges

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

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(((20,25,0),), ((-20,25,0),),
    ((-20,-25,0),), ((20,-25,0),))
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=24, constraint=FIXED)

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

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

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

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

# Assign meshing controls to the respective regions

pickedFaces = myPlateInstance.faces.findAt(((-10,2.5,0),), ((-10,-2.5,0),),
    ((10,2.5,0),), ((10,-2.5,0),))
myAssembly.setMeshControls(regions=pickedFaces, 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='Left')
myModel.historyOutputRequests['Left'].setValues(contourIntegral='Crack-left',
    numberOfContours=8)
myModel.HistoryOutputRequest(name='Right', createStepName='LoadPlate',
    contourIntegral='Crack-right', numberOfContours=8)

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

# Create the job 

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

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