'''
-----------------------------------------------------------------------------
 Axisymmetric model of penny shaped crack in a round bar,
 modeled using axisymmetric quadrilateral reduced integration
 elements (CAX8).
-----------------------------------------------------------------------------
'''   

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 = '2DAxPennyCrackCAX8R'
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='roundBarProfile',sheetSize=200.0)
mySketch.sketchOptions.setValues(viewStyle=AXISYM)
mySketch.setPrimaryObject(option=STANDALONE)

mySketch.ObliqueConstructionLine(point1=(0.0, -100.0),
    point2=(0.0, 100.0))
mySketch.rectangle(point1=(0.0, 0.0), point2=(20.0, 30.0))

myRoundBar = myModel.Part(name='AxiRoundBar', 
    dimensionality=AXISYMMETRIC, type=DEFORMABLE_BODY)
myRoundBar.BaseShell(sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['roundBarProfile']

myViewport.setValues(displayedObject=myRoundBar)

faces = myRoundBar.faces.findAt(((10,10,0),))

# Create a set referring to the whole part

myRoundBar.Set(faces=faces, name='All')

e = myRoundBar.edges
pickedEdges = myRoundBar.edges.findAt(((2.75,0,0),))
myRoundBar.PartitionEdgeByParam(edges=pickedEdges, parameter=0.5)

f = myRoundBar.faces
t = myRoundBar.MakeSketchTransform(sketchPlane=f[0],
    sketchPlaneSide=SIDE1, origin=(10.0, 15.0, 0.0))
mySketch = myModel.Sketch(name='roundBarProfile', 
    sheetSize=72.11, gridSpacing=1.8, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)

myRoundBar.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
r, r1 = mySketch.referenceGeometry, mySketch.referenceVertices

mySketch.VerticalConstructionLine(point=(0.0, -15.0))
mySketch.sketchOptions.setValues(gridOrigin=(0.0, -15.0))
mySketch.ArcByCenterEnds(center=(0.0, -15.0), point1=(0.5, -15.0),
    point2=(-0.5, -15.0), direction=COUNTERCLOCKWISE)
mySketch.ArcByCenterEnds(center=(0.0, -15.0), point1=(4.5, -15.0),
    point2=(-4.5, -15.0), direction=COUNTERCLOCKWISE)

mySketch.Line(point1=(10.0, -5.0), point2=(-10.0, -5.0))
mySketch.Line(point1=(10.0, 5.0), point2=(-10.0, 5.0))
mySketch.Line(point1=(0.0, -15.0), point2=(0.0, 15.0))
mySketch.Line(point1=(0.0, -15.0), point2=(10.0, -5.0))
mySketch.Line(point1=(0.0, -15.0), point2=(-10.0, -5.0))
mySketch.constraintReferences(vertex1=r1.findAt((0.0, -15.0), 1))

faces = myRoundBar.faces.findAt(((10,10,0),))
myRoundBar.PartitionFaceBySketch(faces=faces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del mySketch

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

# Assign material properties

import material
import section

# 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 = myRoundBar.sets['All']

# Assign the above section to the part

myRoundBar.SectionAssignment(region=region,
    sectionName='SolidHomogeneous')

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

# Create an assembly

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

myAssembly.DatumCsysByThreePoints(coordSysType=CYLINDRICAL,
    origin=(0.0, 0.0, 0.0), point1=(1.0, 0.0, 0.0),
    point2=(0.0, 0.0, -1.0))

pickedEdges = myRoundBarInstance.edges.findAt(((10.25,0,0),), ((12.5,0,0),),
    ((17.25,0,0),))
myAssembly.Set(edges=pickedEdges, name='XAxis')

pickedEdges = myRoundBarInstance.edges.findAt(((0,25,0),),((0,15,0),),
    ((0,5,0),))
myAssembly.Set(edges=pickedEdges, name='YAxis')

topEdge = myRoundBarInstance.edges.findAt(((5,30,0),),((15,30,0),))
myAssembly.Set(edges=topEdge, name='topEdge')

topSurf = myAssembly.Surface(side1Edges=topEdge, name='topSurf')

f1 = myRoundBarInstance.faces

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

# Create a step for applying a load

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

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

# Create interaction properties

verts1 = myRoundBarInstance.vertices
vxx = myRoundBarInstance.vertices.findAt((10,0,0))
vyy = myRoundBarInstance.vertices.findAt((20,0,0))
crackTip = verts1[vxx.index:(vxx.index+1)]
crackFront = myAssembly.Set(vertices=crackTip, name='crackFront')

myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=ON, crackFront=crackFront, crackTip=crackFront, 
    extensionDirectionMethod=Q_VECTORS, qVectors=((vxx, vyy),), 
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

XAxis = myAssembly.sets['XAxis']
myModel.DisplacementBC(name='FixedXAxis', createStepName='LoadRoundBar',
    region=XAxis, u2=0.0, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

YAxis = myAssembly.sets['YAxis']
myModel.DisplacementBC(name='FixedYAxis', createStepName='LoadRoundBar',
    region=YAxis, u1=0.0, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

# Assign load conditions

myModel.Pressure(name='Load', createStepName='LoadRoundBar', region=topSurf,
    distribution=UNIFORM, magnitude=-100.0)

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

# Create a mesh 

# Seed all the edges

pickedEdges = myRoundBarInstance.edges.findAt(((0,5,0),), ((20,5,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

pickedEdges = myRoundBarInstance.edges.findAt(((5,10,0),), ((15,10,0),),
    ((5,20,0),), ((15,20,0),), ((5,30,0),), ((15,30,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

pickedEdges = myRoundBarInstance.edges.findAt(((5,10,0),), ((15,10,0),),
    ((5,20,0),), ((15,20,0),), ((5,30,0),), ((15,30,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

pickedEdges = myRoundBarInstance.edges.findAt(((2.75,0,0),),
    ((3.40901,6.59099,0),), ((10,7.25,0),), ((16.59099,6.59099,0),),
    ((17.25,0,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FIXED)

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

pickedEdges = myRoundBarInstance.edges.findAt(((9.75,0,0),),
    ((9.823223,0.176777,0),), ((10,0.25,0),), ((10.176777,0.176777,0),),
    ((10.25,0,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FIXED)

pickedEdges = myRoundBarInstance.edges.findAt(((9.53806,0.191342,0),),
    ((9.808658,0.46194,0),), ((10.191342,0.46194,0),),
    ((10.46194,0.191342,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

pickedEdges1 = myRoundBarInstance.edges.findAt(((10.,1.061012,0.),),
    ((10.779868,0.779868,0.),), ((11.100268,0.,0.),),)
pickedEdges2 = myRoundBarInstance.edges.findAt(((8.905941,0.,0.),),
    ((9.170525,0.829475,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3.0, number=4, constraint=FIXED)

pickedEdges = myRoundBarInstance.edges.findAt(((5.842542,1.722075,0),), 
    ((8.277925,4.157458,0),), ((11.722075,4.157458,0),),
    ((14.157458,1.722075,0),)) 
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=2, constraint=FIXED)

# Assign meshing controls to the respective regions

pickedRegions = myRoundBarInstance.faces.findAt(((10.25,0.1,0),),
    ((10.10,0.25,0),), ((10,0.25,0),), ((9.75,0.1,0),))
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED, 
    technique=SWEEP)

pickedRegions = myRoundBarInstance.faces.findAt(((12.5,1.25,0),),
    ((11,2.5,0),), ((10,2,0),), ((7.5,1.25,0),))
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED, 
    technique=SWEEP)

pickedRegions = myRoundBarInstance.faces.findAt(((5,15,0),), ((15,15,0),),
    ((5,25,0),), ((15,25,0),))
myAssembly.setMeshControls(regions=pickedRegions, technique=SWEEP)

elemType1 = mesh.ElemType(elemCode=CAX8R, elemLibrary=STANDARD, 
    secondOrderAccuracy=OFF, hourglassControl=STIFFNESS,
    distortionControl=OFF)
elemType2 = mesh.ElemType(elemCode=CAX6M, elemLibrary=STANDARD)

faces1 = myRoundBarInstance.faces
pickedRegions = (faces1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2))

partInstances = (myRoundBarInstance, )
myAssembly.generateMesh(regions=partInstances)

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

# Request history output for the crack

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

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

# Create the job 

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

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