'''
-----------------------------------------------------------------------------
 Symmetric model of a two dimensional single edged notch specimen
 for fracture analysis modeled using plane strain elements (CPE8).
-----------------------------------------------------------------------------
'''

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 = '2DSingleEdSymmGlCPE8'
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=(0.0, 0.0), point2=(20.0, 40.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 for the crack

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

# Create a set referring to the whole part

faces1 = myPlate.faces.findAt(((10,20,0),),)
myPlate.Set(faces=faces1, 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)
myPlateInstance = myAssembly.instances['myPlate-1']

# Create circular partitions for sweep mesh around the crack tip

face1 = myPlateInstance.faces.findAt((10,20,0))
t = myAssembly.MakeSketchTransform(sketchPlane=face1, sketchPlaneSide=SIDE1,
    origin=(10.0, 20.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', sheetSize=89.44,
    gridSpacing=2.23, transform=t)
g, v, d = mySketch.geometry, mySketch.vertices, mySketch.dimensions
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(0.0, -20.0))
mySketch.ArcByCenterEnds(center=(0.0, -20.0), point1=(5.0, -20.0),
    point2=(-5.0, -20.0), direction=COUNTERCLOCKWISE)
faces = myPlateInstance.faces
face1 = myPlateInstance.faces.findAt((10,20,0))
pickedFaces = faces[face1.index:(face1.index+1)]
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a set for the crack tip

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

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

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((12.5,0,0))
e2 = myPlateInstance.edges.findAt((17.5,0,0))
edges1 = edges[e1.index:(e1.index+1)]+edges[e2.index:(e2.index+1)]
myAssembly.Set(edges=edges1, name='fixedEdge')

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

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

# Create a set for the top surface

e1 = myPlateInstance.edges.findAt((10,40,0))
side1Edges1 = edges[e1.index:(e1.index+1)]
myAssembly.Surface(side1Edges=side1Edges1, name='loadSurf')

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

# Create a step for applying a load

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

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

# Create interaction properties

# Assign the crack properties

crackFront = crackTip = myAssembly.sets['crackTip']
verts = myPlateInstance.vertices
v11 = myPlateInstance.vertices.findAt((10,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

region = myAssembly.sets['fixedEdge']
myModel.DisplacementBC(name='FixedInY', 
    createStepName='Initial', region=region, u2=0.0,
    fixed=OFF, distributionType=UNIFORM, localCsys=None)

region = myAssembly.sets['fixedPt']
myModel.DisplacementBC(name='FixedInX', 
    createStepName='Initial', region=region, u1=0.0,
    fixed=OFF, distributionType=UNIFORM, localCsys=None)

# Assign load conditions

tSurf = myAssembly.surfaces['loadSurf']
myModel.Pressure(name='TopLoad', createStepName='LoadPlate',
    region=tSurf, distributionType=UNIFORM, magnitude=-1.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=2,
    constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt((10,5,0),)
pickedEdges = e1[pickedEdges.index:(pickedEdges.index+1)]
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=10,
    constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt((17.5,0,0),)
pickedEdges = e1[pickedEdges.index:(pickedEdges.index+1)]
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4,
    constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt((2.5,0,0),)
pickedEdges = e1[pickedEdges.index:(pickedEdges.index+1)]
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=3,
    constraint=FIXED)

pickedEdges1 = myPlateInstance.edges.findAt((20,20,0),)
pickedEdges2 = myPlateInstance.edges.findAt((0,20,0),)
pickedEdges = (e1[pickedEdges1.index:(pickedEdges1.index+1)]
    +  e1[pickedEdges2.index:(pickedEdges2.index+1)])
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=15,
    constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt((10,40,0),)
pickedEdges = e1[pickedEdges.index:(pickedEdges.index+1)]
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=5,
    constraint=FIXED)

# Assign meshing controls to the respective regions

faces1 = myPlateInstance.faces
f1 = myPlateInstance.faces.findAt((10,2.5,0))
pickedRegions = faces1[f1.index:(f1.index+1)]
myAssembly.setMeshControls(regions=pickedRegions,
    elemShape=QUAD_DOMINATED, technique=SWEEP)
elemType1 = mesh.ElemType(elemCode=CPE8, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPE6M, elemLibrary=STANDARD)

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=2)
myModel.HistoryOutputRequest(name='StressInt', createStepName='LoadPlate',
    contourIntegral='Crack', numberOfContours=2,
    contourType=K_FACTORS)
myModel.HistoryOutputRequest(name='Tstr', createStepName='LoadPlate',
    contourIntegral='Crack', numberOfContours=2,
    contourType=T_STRESS)

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

# Create the job 

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

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