'''
-----------------------------------------------------------------------------
 Three dimensional model of an elliptic crack in a half space modeled
 using reduced integration elements (C3D20R).
-----------------------------------------------------------------------------
'''

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
import inpReader

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

# Create a model

Mdb()
modelName = '3DEllipticCrackC3D20R'
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='barProfile',sheetSize=200.0)
mySketch.sketchOptions.setValues(viewStyle=REGULAR)
mySketch.setPrimaryObject(option=STANDALONE)

mySketch.Line(point1=(0.0, 0.0), point2=(100.0, 0.0))
mySketch.Line(point1=(0.0, 0.0), point2=(0.0, -100.0))
mySketch.ArcByCenterEnds(center=(0.0, 0.0), point1=(100.0, 0.0),
    point2=(0.0, -100.0), direction=CLOCKWISE)
myBar = myModel.Part(name='Bar', dimensionality=THREE_D, 
    type=DEFORMABLE_BODY)
myBar.BaseSolidExtrude(sketch=mySketch, depth=50.0)
mySketch.unsetPrimaryObject()
myViewport.setValues(displayedObject=myBar)
del myModel.sketches['barProfile']

# Create a partition for the cracks 

face1 = myBar.faces.findAt((16,-10,50),)
edge1 = myBar.edges.findAt((0,-56.5,50),)
t = myBar.MakeSketchTransform(sketchPlane=face1, sketchUpEdge=edge1,
    sketchPlaneSide=SIDE1, sketchOrientation=LEFT, 
    origin=(42.441318, -42.441318, 50.0))
mySketch = myModel.Sketch(name='barProfile', 
    sheetSize=300.0, gridSpacing=7.5, transform=t)
g = mySketch.geometry
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myBar.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(-42.441318, 42.441318))
mySketch.EllipseByCenterPerimeter(center=(-42.441318, 42.441318),
    axisPoint1=(-17.441318, 42.441318), axisPoint2=(-42.441318, 32.441318))
mySketch.autoTrimCurve(curve=g[3], parameter1=0.737435400485992)
f = myBar.faces
pickedFaces = f[face1.index:(face1.index+1)]
myBar.PartitionFaceBySketch(sketchUpEdge=edge1, faces=pickedFaces,
    sketchOrientation=LEFT, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['barProfile']

# Create a shell tube for the inner partition (around the crack tip)

# First create the sweep path

mySketch1 = myModel.Sketch(name='sweepPath', sheetSize=200.0)
g1 = mySketch1.geometry
mySketch1.setPrimaryObject(option=STANDALONE)
mySketch1.EllipseByCenterPerimeter(center=(0.0, 0.0),
    axisPoint1=(25.0, 0.0), axisPoint2=(0.0, -10.0))
mySketch1.Line(point1=(0.0, 0.0), point2=(30.0, 0.0))
mySketch1.Line(point1=(0.0, 0.0), point2=(0.0, 23.4493026733398))
mySketch1.autoTrimCurve(curve=g1[3], parameter1=0.570547997951508)
mySketch1.delete(objectList=(g1[5], g1[7]))
mySketch1.unsetPrimaryObject()

# Create a circle and sweep it along the curve created above

mySketch2 = myModel.Sketch(name='innerCircleProfile', sheetSize=200.0, 
    transform=(-1.0, 6.12303176911188e-16, 0.0, 0.0, 0.0, 1.0, 
    6.12303176911188e-16, 1.0, 0.0, 25.0, -2.44921270764475e-15, 0.0))
g2 = mySketch2.geometry
mySketch2.setPrimaryObject(option=SUPERIMPOSE)
mySketch2.ObliqueConstructionLine(point1=(-94.28, 0.0),
    point2=(94.28, 0.0))
mySketch2.ObliqueConstructionLine(point1=(0.0, -94.28),
    point2=(0.0, 94.28))
mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0),
    point1=(1.0, 0.0))
myInnerTube = myModel.Part(name='InnerPartition',
    dimensionality=THREE_D, type=DEFORMABLE_BODY)
myInnerTube.BaseShellSweep(sketch=mySketch2, path=mySketch1)
mySketch2.unsetPrimaryObject()
myViewport.setValues(displayedObject=myInnerTube)
del myModel.sketches['innerCircleProfile']

# Create a shell tube for the outer partition (around the inner tube)

# For outer partition the sweep path is same as inner partition

mySketch2 = myModel.Sketch(name='outerCirclePartition', sheetSize=200.0, 
    transform=(-1.0, 6.12303176911188e-16, 0.0, 0.0, 0.0, 1.0, 
    6.12303176911188e-16, 1.0, 0.0, 25.0, -2.44921270764475e-15, 0.0))
g3 = mySketch2.geometry
mySketch2.setPrimaryObject(option=SUPERIMPOSE)
mySketch2.ObliqueConstructionLine(point1=(-94.28, 0.0),
    point2=(94.28, 0.0))
mySketch2.ObliqueConstructionLine(point1=(0.0, -94.28),
    point2=(0.0, 94.28))
mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0),
    point1=(3.0, 0.0))
myOuterTube = myModel.Part(name='OuterPartition', 
    dimensionality=THREE_D, type=DEFORMABLE_BODY)
myOuterTube.BaseShellSweep(sketch=mySketch2, path=mySketch1)
mySketch2.unsetPrimaryObject()
myViewport.setValues(displayedObject=myOuterTube)
del myModel.sketches['outerCirclePartition']
del myModel.sketches['sweepPath']

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

# Create an assembly

# Create an instance of the bar

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

# Create an instance of the inner partition tube

myAssembly.Instance(name='myInnerPartition-1',
    part=myInnerTube)
myInnerTubeInstance = myAssembly.instances['myInnerPartition-1']

# Position the inner tube around the crack tip

myInnerTubeInstance.translate(vector=
    (3.5527136788005e-15, 2.44921270764475e-15, 50.0))
myInnerTubeInstance.rotateAboutAxis(axisPoint=(100.0, 0.0, 50.0),
    axisDirection=(-100.0, 0.0, 0.0), angle=180.0)

# Create an instance of the outer partition tube

myAssembly.Instance(name='myOuterPartition-1', part=myOuterTube)
myOuterTubeInstance = myAssembly.instances['myOuterPartition-1']

# Position the outer tube around the crack tip

myOuterTubeInstance.translate(vector=
    (3.5527136788005e-15, 2.44921270764475e-15, 50.0))
myOuterTubeInstance.rotateAboutAxis(axisPoint=(100.0, 0.0, 50.0),
    axisDirection=(-100.0, 0.0, 0.0), angle=180.0)

# Subtract the inner and outer tubes from the bar

myAssembly.PartFromBooleanCut(name='PartitionedBar', 
    instanceToBeCut=myBarInstance, 
    cuttingInstances=(myInnerTubeInstance, myOuterTubeInstance, ))

# Create an instance of the partitioned bar

myPartitionedBar = myModel.parts['PartitionedBar']
myAssembly.Instance(name='myPartitionedBar-1', part=myPartitionedBar)
myNewBarInstance = myAssembly.instances['myPartitionedBar-1']

# Create a new assembly of the parts to be suppressed

myAssembly1 = myModel.rootAssembly
myAssembly1.suppressFeatures(('myInnerPartition-1',
    'myOuterPartition-1', 'myBar-1', ))

# Create a set for the new bar instance

cells = myNewBarInstance.cells
myAssembly.Set(cells=cells, name='FullPart')

# Create a partition to make the entire part meshable

pickedEdges1 = myNewBarInstance.edges.findAt((0,0,25.0))
pickedEdges2 = (myNewBarInstance.edges.findAt((14.230512,-8.221849,50.)),)

myAssembly.PartitionCellByExtrudeEdge(line=pickedEdges1, cells=cells,
    edges=pickedEdges2, sense=REVERSE)

myViewport.setValues(displayedObject=myAssembly)

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

# Assign material properties

myViewport.setValues(displayedObject=myBar)

# Create a set for the bar

c = myPartitionedBar.cells
myPartitionedBar.Set(cells=c, name='All')

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

# Assign the above section to the part

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

# Create a step for applying a load

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

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

# Create interaction properties

# Create a set for the crack front

edges = myNewBarInstance.edges
edges1 = myNewBarInstance.edges.findAt((14.230512,-8.221849,50.))
edges1 = edges[edges1.index:(edges1.index+1)]
myAssembly.Set(edges=edges1, name='CrackFront')

# Create the contour integral definition for the top crack

crackFront = crackTip = myAssembly.sets['CrackFront']
v1 = myNewBarInstance.vertices.findAt((0,0,50))
v2 = myNewBarInstance.vertices.findAt((0,0,0))
myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=ON, crackFront=crackFront, crackTip=crackTip, 
    extensionDirectionMethod=CRACK_NORMAL, crackNormal=(v1, v2), 
                                               midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE,crackTipSense=REVERSE)

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

# Create loads and boundary conditions 

# Create a set for the mid plane

faces = myNewBarInstance.faces
allFaces = (myNewBarInstance.faces.findAt(((60,-60,50),),
    ((15,-9,50),), ((15.5,-10,50),)))
myAssembly.Set(faces=allFaces, name='midPlane')

# Create a set for the side plane

faces = myNewBarInstance.faces
allFaces = (myNewBarInstance.faces.findAt(((0,-55,25),),
    ((0,-3.5,25),), ((0,-11.75,49),), ((0,-8.25,49),),
    ((0,-10.25,49.8),), ((0,-9.75,49.8),)))
myAssembly.Set(faces=allFaces, name='sidePlane')

# Create a set for the top plane

faces = myNewBarInstance.faces
allFaces = (myNewBarInstance.faces.findAt(((50,-55,0),),
    ((7,-5,0),)))
myAssembly.Set(faces=allFaces, name='topPlane')

# Create a set for a point to be fixed in Y

vert = myNewBarInstance.vertices
v1 = myNewBarInstance.vertices.findAt((0,-100,50))
v1 = vert[v1.index:(v1.index+1)]
myAssembly.Set(vertices=v1, name='endPt')

# Fix the symmetry plane in the Z direction

region = myAssembly.sets['midPlane']
myModel.DisplacementBC(name='ZFixed', createStepName='ApplyLoad',
    region=region, u3=0.0, fixed=OFF, distributionType=UNIFORM,
    localCsys=None)

# Fix the side plane in the X direction

region = myAssembly.sets['sidePlane']
myModel.DisplacementBC(name='XFixed', createStepName='ApplyLoad',
    region=region, u1=0.0, fixed=OFF, 
    distributionType=UNIFORM, localCsys=None)

# Fix a point in the Y direction

region = myAssembly.sets['endPt']
myModel.DisplacementBC(name='YFixed', 
    createStepName='ApplyLoad', region=region, u1=UNSET, u2=0.0, u3=UNSET, 
    ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude=UNSET, fixed=OFF, 
    distributionType=UNIFORM, localCsys=None)

# Create a set for the top surface

faces = myNewBarInstance.faces
allFaces = (myNewBarInstance.faces.findAt(((50,-55,0),),
    ((7,-5,0),)))
myAssembly.Surface(side1Faces=allFaces, name='topSurface')

# Apply a pressure load on the top face

region = myAssembly.surfaces['topSurface']
myModel.Pressure(name='PressureLoad', createStepName='ApplyLoad',
    region=region, distributionType=UNIFORM, 
    magnitude=-100.0, amplitude=UNSET)

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

# Create a mesh 

# Seed all the edges

pickedEdges1 = (myNewBarInstance.edges.findAt(((64,0,50),),((62.5,0,0),),
    ((0,-55,0),),((0,-56.5,50),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=18, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((0,-10,23.5),),((25,0,23.5),),
    ((100,0,25),), ((0,-100,25),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=10)

pickedEdges1 = myNewBarInstance.edges.findAt(((0,0,25),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=12, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((70.710678,-70.710678,50),),
    ((70.710678,-70.710678,0),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=10)

pickedEdges1 = myNewBarInstance.edges.findAt(((0,-5,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=4)

pickedEdges1 = myNewBarInstance.edges.findAt(((12.5,0,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=8)

pickedEdges1 = myNewBarInstance.edges.findAt(((14.230512,-8.221849,0),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

pickedEdges1 = myNewBarInstance.edges.findAt(((0,-3.5,50),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=3)

pickedEdges1 = myNewBarInstance.edges.findAt(((11,0,50),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5)

pickedEdges1 = (myNewBarInstance.edges.findAt(((27.121325,0,47.878684),),
    ((0,-12.12132,47.87868),), ((0,-7.87868,47.87868),),
    ((22.878679,0,47.87868),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=3, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((0,-11.75,50),),
    ((23.25,0,50),), ((0,-8.25,50),), ((26.75,0,50),), ((0,-10,48.25),),
    ((25,0,48.25),),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((0,-10,49.75),),
    ((25,0,49.75),), ((0,-10.25,50),), ((24.75,0,50),), ((0,-9.75,50),),
    ((25.25,0,50),),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=1, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((0,-10.707107,49.292893),),
    ((24.292917,0,49.29287),), ((25.707149,0,49.292935),),
    ((0,-9.292893,49.292893),),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=3, constraint=FIXED)

pickedEdges1 = (myNewBarInstance.edges.findAt(((14.992905,-9.044854,50),),
    ((13.466124,-7.3923,50),), ((14.230512,-8.221849,50),),
    ((16.512692,-10.672038,50),), ((11.932046,-5.712304,50),),
    ((14.230498,-8.221853,47),), ((14.230512,-8.221849,0),),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

# Assign meshing controls to the respective cells

cells = myNewBarInstance.cells.findAt(((14,-8,50),), ((14.4,-8.5,50),))
myAssembly.setMeshControls(regions=cells, elemShape=WEDGE, technique=SWEEP)

# All the remaining cells will be meshed with hex elements using
# structured meshing. This is also the default
 
elemType1 = mesh.ElemType(elemCode=C3D20R, elemLibrary=STANDARD, 
    kinematicSplit=AVERAGE_STRAIN, secondOrderAccuracy=OFF, 
    hourglassControl=STIFFNESS, distortionControl=OFF)
elemType2 = mesh.ElemType(elemCode=C3D15, elemLibrary=STANDARD)
elemType3 = mesh.ElemType(elemCode=C3D10M, elemLibrary=STANDARD)
cells1 = myNewBarInstance.cells
pickedRegions =(cells1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2, elemType3))
partInstances = (myNewBarInstance, )
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=3)
myModel.HistoryOutputRequest(name='StressInt', createStepName='ApplyLoad',
    contourIntegral='Crack', numberOfContours=3, contourType=K_FACTORS)
myModel.HistoryOutputRequest(name='TStr', createStepName='ApplyLoad',
    contourIntegral='Crack', numberOfContours=3, 
    contourType=T_STRESS)

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

# Create a job and write an input file to get the orphan mesh

mdb.Job(name=modelName, model=modelName, type=ANALYSIS, 
    description='Submit the job for analysis')
mdb.saveAs(pathName=modelName)

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