'''
-----------------------------------------------------------------------------
 Semi elliptic crack in a rectangular plate modeled using continuum
 elements (C3D8).
-----------------------------------------------------------------------------
'''

from abaqus import *
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 = 'RectBlEllipticCrC3D8'
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.rectangle(point1=(0.0, 0.0), point2=(13.333, 26.667))
myBlock = myModel.Part(name='block', dimensionality=THREE_D, 
    type=DEFORMABLE_BODY)
myBlock.BaseSolidExtrude(sketch=mySketch, depth=1.6667)
mySketch.unsetPrimaryObject()
del myModel.sketches['barProfile']

# Create a partition for the crack line 

face1 = myBlock.faces.findAt((5,0,0.6),)
edge1 = myBlock.edges.findAt((13.333,0,0.83335),)
t = myBlock.MakeSketchTransform(sketchPlane=face1, sketchUpEdge=edge1,
    sketchPlaneSide=SIDE1, origin=(6.6665, 0.0, 0.83335))
mySketch = myModel.Sketch(name='barProfile', 
    sheetSize=26.87, gridSpacing=0.67, transform=t)
g = mySketch.geometry
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myBlock.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(-6.6665, 0.83335))
mySketch.EllipseByCenterPerimeter(center=(-6.6665, 0.83335),
    axisPoint1=(-2.4998, 0.83335), axisPoint2=(-6.6665, -0.16665))
f = myBlock.faces
pickedFaces = f[face1.index:(face1.index+1)]
myBlock.PartitionFaceBySketch(sketchUpEdge=edge1, faces=pickedFaces,
    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=50.0)
g1 = mySketch1.geometry
mySketch1.setPrimaryObject(option=STANDALONE)
mySketch1.EllipseByCenterPerimeter(center=(0.0, 0.0),
    axisPoint1=(4.1667, 0.0), axisPoint2=(0.0, -1.0))
mySketch1.Line(point1=(0.0, 0.0), point2=(0.0, -5.0))
mySketch1.Line(point1=(0.0, 0.0), point2=(6.0, 0.0))
mySketch1.autoTrimCurve(curve=g1[3], parameter1=0.655947506427765)
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=50.0, 
    transform=(-1.46951586845991e-17, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 
    1.46951586845991e-17, 0.0, 2.55128364723585e-16, -1.0, 0.0))
g2 = mySketch2.geometry
mySketch2.setPrimaryObject(option=SUPERIMPOSE)
mySketch2.ObliqueConstructionLine(point1=(-23.57, 0.0),
    point2=(23.57, 0.0))
mySketch2.ObliqueConstructionLine(point1=(0.0, -23.57),
    point2=(0.0, 23.57))
mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0),
    point1=(0.005, 0.0))
myInnerTube = myModel.Part(name='innerTube',
    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='outerCircleProfile', sheetSize=50.0, 
    transform=(-1.46951586845991e-17, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 
    1.46951586845991e-17, 0.0, -7.65385094170755e-16, -1.0, 0.0))
g3 = mySketch2.geometry
mySketch2.setPrimaryObject(option=SUPERIMPOSE)
mySketch2.ObliqueConstructionLine(point1=(-23.57, 0.0),
    point2=(23.57, 0.0))
mySketch2.ObliqueConstructionLine(point1=(0.0, -23.57),
    point2=(0.0, 23.57))
mySketch2.CircleByCenterPerimeter(center=(0.0, 0.0),
    point1=(0.2, 0.0))
myOuterTube = myModel.Part(name='outerTube', 
    dimensionality=THREE_D, type=DEFORMABLE_BODY)
myOuterTube.BaseShellSweep(sketch=mySketch2, path=mySketch1)
mySketch2.unsetPrimaryObject()
myViewport.setValues(displayedObject=myOuterTube)
del myModel.sketches['outerCircleProfile']
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='myBlock-1', part=myBlock)
myBlockInstance = myAssembly.instances['myBlock-1']

# Create an instance of the outer partition tube

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

# Position the outer tube around the crack tip

myOuterTubeInstance.rotateAboutAxis(axisPoint=(0.0, 0.0, 0.0),
    axisDirection=(13.333, 0.0, 0.0), angle=90.0)
myOuterTubeInstance.translate(vector=
    (8.88178419700125e-16, 0.0, 1.6667))

# Create an instance of the inner partition tube

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

# Position the inner tube around the crack tip

myInnerTubeInstance.rotateAboutAxis(axisPoint=(0.0, 0.0, 0.0),
    axisDirection=(13.333, 0.0, 0.0), angle=90.0)
myInnerTubeInstance.translate(vector=
    (8.88178419700125e-16, 0.0, 1.6667))

# Subtract the inner and outer tubes from the block

myAssembly.PartFromBooleanCut(name='partitionedBlock', 
    instanceToBeCut=myBlockInstance, 
    cuttingInstances=(myInnerTubeInstance, myOuterTubeInstance, ))

# Create an instance of the partitioned block

myPtBlock = myModel.parts['partitionedBlock']
myAssembly.Instance(name='myPtBlock-1', part=myPtBlock, dependent=OFF)
myPtBlInstance = myAssembly.instances['myPtBlock-1']

# Create a new assembly of the parts to be suppressed

myAssembly1 = myModel.rootAssembly
myAssembly1.suppressFeatures(('myInnerTube-1',
    'myOuterTube-1', 'myBlock-1', ))

# Create a set for the new bar instance

cells = myPtBlock.cells
myPtBlock.Set(cells=cells, name='All')

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

# Assign material properties

myViewport.setValues(displayedObject=myPtBlock)

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

# Assign the above section to the part

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

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

# Create a step for applying a load

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

# Create a partition to make the entire part meshable

c1 = myPtBlInstance.cells
e1 = myPtBlInstance.edges.findAt((2.217181,0.,0.82003))
pickedEdges = (e1,)
sPath = myPtBlInstance.edges.findAt((0.,13.3335,1.6667))

myAssembly.PartitionCellBySweepEdge(sweepPath=sPath, cells=c1,
    edges=pickedEdges)

myViewport.setValues(displayedObject=myAssembly)

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

# Create interaction properties

# Create a set for the top surface

surf = myPtBlInstance.faces
surf1 = myPtBlInstance.faces.findAt((6.6665,26.667,0.6))
surf2 = myPtBlInstance.faces.findAt((2.08335,26.667,0.9))
surfAll = (surf[surf1.index:(surf1.index+1)]
    + surf[surf2.index:(surf2.index+1)])
myAssembly.Surface(side1Faces=surfAll, name='topSurf')

# Create a set for the X face

faces1 = (myPtBlInstance.faces.findAt(((6.6665,0,0.8),), ((2.3,0,0.7),),
    ((2.179725,0.,0.81192),),))
myAssembly.Set(faces=faces1, name='xFace')

# Create a set for Y face

faces1 = (myPtBlInstance.faces.findAt(((0,13.3335,0.4),),
    ((0,13.3335,1.0),), ((0,0.1,0.541),), ((0,0.1,0.808),),
    ((0,0.02,0.641),), ((0,0.02,0.691),),))
myAssembly.Set(faces=faces1, name='yFace')

# Create a set for the vertex to be fixed in Z

v1 = myPtBlInstance.vertices.findAt(((0.,0.,1.6667),),)
myAssembly.Set(vertices=v1, name='pt')

# Create a set for the crack line

edges = myPtBlInstance.edges
edges1 = myPtBlInstance.edges.findAt((2.217181,0.,0.82003))
edges1 = edges[edges1.index:(edges1.index+1)]
myAssembly.Set(edges=edges1, name='crackLine')

# Create the contour integral definition for the crack

crackFront = crackTip = myAssembly.sets['crackLine']
v1 = myPtBlInstance.vertices.findAt((0.,0.,0.6667))
v2 = myPtBlInstance.vertices.findAt((0.,26.667,0.6667))
myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=ON, crackFront=crackFront, crackTip=crackTip, 
    extensionDirectionMethod=CRACK_NORMAL, crackNormal=(v1, v2), 
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions

# Fix the X face in the Y directions

region = myAssembly.sets['xFace']
myModel.DisplacementBC(name='yFixed', createStepName='Initial',
    region=region, u2=0.0, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

# Fix the Y face in the X direction

region = myAssembly.sets['yFace']
myModel.DisplacementBC(name='xFixed', createStepName='Initial',
    region=region, u1=0.0, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

# Fix a vertex in the Z direction

region = myAssembly.sets['pt']
myModel.DisplacementBC(name='zFixed', createStepName='Initial',
    region=region, u3=0.0, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

# Apply a pressure load on the top face

region = myAssembly.surfaces['topSurf']
myModel.Pressure(name='prLoad', createStepName='ApplyLoad',
    region=region, distribution=UNIFORM, 
    magnitude=-1.0, amplitude=UNSET)

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

# Create a mesh 

# Seed all the edges

pickedEdges1 = (myPtBlInstance.edges.findAt(((4.308121,0.141421,1.6667),),
    ((4.025279,0.141421,1.6667),), ((0.,0.141396,0.808146),),
    ((0.,0.141419,0.525276),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=6, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((4.1667,0.025,1.6667),),
    ((4.1417,0,1.6667),), ((4.1917,0,1.6667),), ((0,0.025,0.6667),),
    ((0.,0.,0.6417),), ((0.,0.,0.6917),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=1, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((0.,0.125,0.6667),),
    ((0.,0.,0.5417),), ((0.,0.,0.7917),), ((4.1667,0.125,1.6667),),
    ((4.0417,0.,1.6667),), ((4.2917,0.,1.6667),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((2.217181,0.,0.82003),),
    ((1.090774,0.,0.696563),), ((2.372976,0.,0.642024),),
    ((1.089261,0.,0.701475),), ((2.06162,0.,0.999576),),
    ((1.087582,0.,0.706376),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((1.98335,0.,1.6667),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,0.68455),))
pickedEdges2 = (myPtBlInstance.edges.findAt(((0.,0.014928,0.6667),),
    ((0.,0.,0.651772),)))
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3.0, number=10)
'''
pickedEdges1 = (myPtBlInstance.edges.findAt(((4.186988,0.,1.6667),),
    ((4.1667,0.030668,1.6667),)))
pickedEdges2 = myPtBlInstance.edges.findAt(((4.141774,0.,1.6667),))
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3.0, number=10)
'''

pickedEdges1 = myPtBlInstance.edges.findAt(((4.141774,0.,1.6667),))
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=3.0,
    number=10, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((4.1667,0.024933,1.6667),))
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=3.0,
    number=10, constraint=FIXED)

pickedEdges2 = myPtBlInstance.edges.findAt(((4.181774,0.,1.6667),))
myAssembly.seedEdgeByBias(end2Edges=pickedEdges2, ratio=3.0,
    number=10, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,1.2667),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5)

pickedEdges1 = myPtBlInstance.edges.findAt(((0.,0.,0.23335),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((13.333,0.,0.83335),),
    ((13.333,26.667,0.83335),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5)

pickedEdges1 = (myPtBlInstance.edges.findAt(((8.84985,0.,1.6667),),
    ((8.74985,26.667,1.6667),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=15, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((6.6665,0.,0),),
    ((6.6665,26.667,0),)))                                            
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=25, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((2.217181,26.667,0.82003),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((2.08335,26.667,1.6667),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=20, constraint=FIXED)

pickedEdges1 = myPtBlInstance.edges.findAt(((0.,26.667,1.1667),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5)

pickedEdges1 = myPtBlInstance.edges.findAt(((0.,26.667,0.33335),))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=5, constraint=FIXED)

pickedEdges1 = (myPtBlInstance.edges.findAt(((4.1667,13.4335,1.6667),),
    ((0.,13.4335,0.6667),), ((13.333,13.3335,1.6667),),
    ((0.,13.3335,1.6667),), ((0.,13.3335,0.),),
    ((13.333,13.3335,0.),)))
myAssembly.seedEdgeByNumber(edges=pickedEdges1, number=30, constraint=FIXED)

# Assign meshing controls to the respective cells

cells = myPtBlInstance.cells.findAt(((2.167847,0.,0.916656),),
    ((2.378009,0.,0.742006),),)
myAssembly.setMeshControls(regions=cells, elemShape=HEX, technique=STRUCTURED)

cells = myPtBlInstance.cells.findAt(((2.271117,0.,0.830836),),
    ((2.276242,0.,0.826576),),)
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=C3D8, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=C3D6, elemLibrary=STANDARD)
elemType3 = mesh.ElemType(elemCode=C3D4, elemLibrary=STANDARD)
cells1 = myPtBlInstance.cells
pickedRegions =(cells1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2, elemType3))
partInstances = (myPtBlInstance, )
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=6)
myModel.HistoryOutputRequest(name='StrInt',
    createStepName='ApplyLoad', contourIntegral='Crack',
    numberOfContours=6, contourType=K_FACTORS)
myModel.HistoryOutputRequest(name='TStr',
    createStepName='ApplyLoad', contourIntegral='Crack',
    numberOfContours=6, contourType=T_STRESS)

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

# Create the job 

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

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