'''
-----------------------------------------------------------------------------
 Symmetric three dimensional model of a single edged notch specimen
 modeled using reduced integration continuum elements (C3D20R).
 
 Fracture analysis is done on this model subjected to thermal stresses
 obtained by running a Heat transfer analysis on this model.
 
-----------------------------------------------------------------------------
'''

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 = 'SingleEdThMesh2C3D20R'
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=(-10.0, 0.0), point2=(10.0, 40.0))

myPlate = myModel.Part(name='Plate', 
    dimensionality=THREE_D, type=DEFORMABLE_BODY)
myPlate.BaseSolidExtrude(sketch=mySketch, depth=1.0)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

myViewport.setValues(displayedObject=myPlate)

# Partition the edge for the crack

pickedEdge = myPlate.edges.findAt((1,0,0))
myPlate.PartitionEdgeByParam(edges=pickedEdge, parameter=0.5)

# Partition the plate

face1 = myPlate.faces.findAt((0,20,1))
edge1 = myPlate.edges.findAt((10.,20.,1.))
t = myPlate.MakeSketchTransform(sketchPlane=face1,
    sketchUpEdge=edge1, sketchPlaneSide=SIDE1,
    origin=(0.0,20.0,1.0))

mySketch = myModel.Sketch(name='partitionProfile',
    sheetSize=89.46, gridSpacing=2.23, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myPlate.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.sketchOptions.setValues(gridOrigin=(0.0, -20.0))
mySketch.ArcByCenterEnds(center=(0.0,-20.0), point1=(0.28,-20.0),
    point2=(-0.28,-20.0), direction=COUNTERCLOCKWISE)
mySketch.Line(point1=(10.0,-10.0), point2=(-10.0,-10.0))
mySketch.Line(point1=(0.0,-20.0), point2=(10.0,-10.0))
mySketch.Line(point1=(0.0,-20.0), point2=(-10.0,-10.0))

pickedFaces = myPlate.faces.findAt((0,20,1))
myPlate.PartitionFaceBySketch(sketchUpEdge=edge1,
    faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['partitionProfile']

# Extrude the partitions created above to partition
# the entire part

pickedCells = myPlate.cells
e1 = myPlate.edges.findAt((0.258686,0.107151,1.))
e2 = myPlate.edges.findAt((0.,0.280,1.))
e3 = myPlate.edges.findAt((-0.258686,0.107151,1.))
sPath = myPlate.edges.findAt((10,0,0.25))
pickedEdges =(e1,e2,e3)
myPlate.PartitionCellBySweepEdge(sweepPath=sPath, cells=pickedCells, 
    edges=pickedEdges)

pickedCells = myPlate.cells
e1 = myPlate.edges.findAt((0.098995,0.098995,1.))
e2 = myPlate.edges.findAt((5.098995,5.098995,1.))
e3 = myPlate.edges.findAt((0.,10.,1.))
e4 = myPlate.edges.findAt((-5.098995,5.098995,1.))
e5 = myPlate.edges.findAt((-0.098995,0.098995,1.))
pickedEdges = (e1,e2,e3,e4,e5)
sPath = myPlate.edges.findAt((10,0,0.25))
myPlate.PartitionCellBySweepEdge(sweepPath=sPath, cells=pickedCells, 
    edges=pickedEdges)

# Create a set for the entire part

pickedCells = myPlate.cells
myPlate.Set(cells=pickedCells, name='fullPart')

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

# Assign material properties

# Create linear elastic material

myModel.Material(name='LinearElastic')
myModel.materials['LinearElastic'].Elastic(table=((30000000.0,0.3),))
myModel.materials['LinearElastic'].Expansion(table=((7.5e-6,),))
myModel.HomogeneousSolidSection(name='SolidHomogeneous',
    material='LinearElastic', thickness=1.0)

region = myPlate.sets['fullPart']

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

# Create a set for the entire instance

pickedCells = myPlateInstance.cells
myAssembly.Set(cells=pickedCells, name='All')

# Create a set for the X face

faces1 = myPlateInstance.faces.findAt(((-0.140,0.,0.500),),
    ((-5.14,0.,0.500),),)
myAssembly.Set(faces=faces1, name='xFace')

# Create a set for the top face

faces1 = myPlateInstance.faces.findAt(((0.,40.,0.5),),)
myAssembly.Set(faces=faces1, name='topFace')

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

edge1 = myPlateInstance.edges.findAt(((-10.,40.,0.500),),)
myAssembly.Set(edges=edge1, name='topEdge')

# Create a set for the crack line

edge1 = myPlateInstance.edges.findAt(((0.,0.,0.500),),)
myAssembly.Set(edges=edge1, name='crackLine')

# Create a set for the front face

faces1 = myPlateInstance.faces.findAt(((0.,25.,1.),), ((-5.14,2.57,1.),),
    ((-0.020503,5.119497,1.),), ((5.098995,2.549497,1.),),
    ((-0.140,0.070,1.),), ((0.,0.140,1.),),
    ((0.140,0.070,1.),),)
myAssembly.Set(faces=faces1, name='frontFace')

# Create a set for the back face

faces1 = myPlateInstance.faces.findAt(((0.140,0.070,0.),),
    ((0.,0.140,0.),), ((-0.140,0.070,0.),),
    ((5.970671,2.985336,0.),), ((0.,5.970671,0.),),
    ((-5.970671,2.985336,0.),), ((0.,25.,0.),),)
myAssembly.Set(faces=faces1, name='backFace')


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

# Create a step for applying the thermal stresses
# obtained in the heat transfer analysis

myModel.StaticStep(name='ApplyThermalLoad', previous='Initial',
    description='Apply the thermal stresses')

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

# Create interaction properties

v1 = myPlateInstance.vertices.findAt((0,0,1),)
v2 = myPlateInstance.vertices.findAt((-10,0,1),)
crackFront = crackTip = myAssembly.sets['crackLine']
myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=ON, crackFront=crackFront, crackTip=crackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((v1, v2),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

region = myAssembly.sets['xFace']
myModel.DisplacementBC(name='yFixed', createStepName='Initial',
    region=region, u2=SET, distributionType=UNIFORM, localCsys=None)

region = myAssembly.sets['topFace']
myModel.DisplacementBC(name='topFaceYFixed', createStepName='Initial',
    region=region, u2=SET, distributionType=UNIFORM, localCsys=None)

region = myAssembly.sets['topEdge']
myModel.DisplacementBC(name='topEdgeXFixed', createStepName='Initial',
    region=region, u1=SET, distributionType=UNIFORM, localCsys=None)

region = myAssembly.sets['backFace']
myModel.DisplacementBC(name='backFaceZFixed', createStepName='Initial',
    region=region, u3=SET, distributionType=UNIFORM, localCsys=None)

region = myAssembly.sets['frontFace']
myModel.DisplacementBC(name='frontFaceZFixed', createStepName='Initial',
    region=region, u3=SET, distributionType=UNIFORM, localCsys=None)

# Create a temperature field.
# The *.odb file from the heat transfer analysis will be used
# to drive the analysis.

myModel.Temperature(name='TempField',
    createStepName='ApplyThermalLoad',
    distributionType=FROM_FILE,
    fileName='3DSingleEdgedThermMesh2.odb',
    beginStep=None, beginIncrement=None, endStep=None,
    endIncrement=None, interpolate=ON,
    absoluteExteriorTolerance=0.0, exteriorTolerance=0.05)

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

# Create a mesh

# Assign meshing controls to different regions

pickedRegions = myPlateInstance.cells.findAt(((0.,0.140,0.500),),
    ((0.140,0.,0.500),), ((-0.140,0.,0.500),),)
myAssembly.setMeshControls(regions=pickedRegions,
    elemShape=WEDGE, technique=SWEEP)

pickedRegions = myPlateInstance.cells.findAt(((0.,25.,0.500),),
    ((0.,5.14,0.500),), ((-10.,5.,0.500),), ((10.,5.,0.500),),)
myAssembly.setMeshControls(regions=pickedRegions,
    elemShape=HEX, technique=STRUCTURED)

# Seed all the edges

pickedEdges = myPlateInstance.edges.findAt(((0.098995,0.098995,1.),),
    ((0.098995,0.098995,0.),), ((-0.098995,0.098995,0.),),
    ((-0.098995,0.098995,1.),), ((0.140,0.,0.),),
    ((0.140,0.,1.),), ((-0.140,0.,0.),), ((-0.140,0.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=1, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((-0.19799,0.19799,0.500),),
    ((0.19799,0.19799,0.500),), ((0.,0.,0.500),),
    ((0.280,0.,0.500),), ((-0.280,0.,0.500),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=1, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.258686,0.107151,0.),),
    ((-0.258686,0.107151,0.),), ((0.258686,0.107151,1.),),
    ((-0.258686,0.107151,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=3, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.,0.280,0.),),
    ((0.,0.280,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=6, constraint=FIXED)

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

pickedEdges = myPlateInstance.edges.findAt(((0.,10.,0.),),
    ((0.,10.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=6, constraint=FIXED)

pickedEdges1 = myPlateInstance.edges.findAt(((-2.033918,0.,1.),),
    ((-1.070671,1.070671,1.),), ((1.070671,1.070671,1.),),)
pickedEdges2 = myPlateInstance.edges.findAt(((2.141342,0.,1.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=4, number=5, constraint=FIXED)

pickedEdges1 = myPlateInstance.edges.findAt(((2.141844,0.,0.),),
    ((1.070671,1.070671,0.),),)
pickedEdges2 = myPlateInstance.edges.findAt(((-2.141342,0.,0.),),
    ((-1.070671,1.070671,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=4, number=5, constraint=FIXED)

pickedEdges2 = myPlateInstance.edges.findAt(((10.,13.668598,1.),),)
pickedEdges1 = myPlateInstance.edges.findAt(((10.,14.075607,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3, number=3, constraint=FIXED)

pickedEdges1 = myPlateInstance.edges.findAt(((-10.,14.784095,1.),),)
pickedEdges2 = myPlateInstance.edges.findAt(((-10.,14.698612,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=3, number=3, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((0.,40.,1.),),
    ((0.,40.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges,
    number=6, constraint=FIXED)

elemType1 = mesh.ElemType(elemCode=C3D20R, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=C3D15, elemLibrary=STANDARD)
elemType3 = mesh.ElemType(elemCode=C3D10M, elemLibrary=STANDARD)
cells1 = myPlateInstance.cells
pickedRegions =(cells1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1,elemType2,elemType3))
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=7)

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

# Create the job and submit it for analysis

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

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










