'''
-----------------------------------------------------------------------------
 Three dimensional symmetric model of a double edged notch specimen (linear
 elastic material) modeled using reduced integration continuum
 solid elements (C3D20R).
-----------------------------------------------------------------------------
'''

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

mySketch.rectangle(point1=(0.0, 0.0), point2=(20.0, 50.0))

myBlock = myModel.Part(name='Block', 
    dimensionality=THREE_D, type=DEFORMABLE_BODY)
myBlock.BaseShell(sketch=mySketch)
myBlock.BaseSolidExtrude(sketch=mySketch, depth=1.0)
mySketch.unsetPrimaryObject()
del myModel.sketches['blockProfile']

myViewport.setValues(displayedObject=myBlock)

# Create a set for the entire part

pickedCells = myBlock.cells
myBlock.Set(cells=pickedCells, 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 = myBlock.sets['All']

# Assign the above section to the part

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

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

# Create an assembly

myAssembly = myModel.rootAssembly
myViewport.setValues(displayedObject=myAssembly)
myAssembly.DatumCsysByDefault(CARTESIAN)
myAssembly.Instance(name='myBlock-1', part=myBlock, dependent=OFF)
myBlockInstance = myAssembly.instances['myBlock-1']

# Partition the face on the X-Z plane

face = myBlockInstance.faces.findAt((10.,0.,0.500))
edge = myBlockInstance.edges.findAt((20.,0.,0.500))
t = myAssembly.MakeSketchTransform(sketchPlane=face, sketchUpEdge=edge, 
    sketchPlaneSide=SIDE1, origin=(10.0, 0.0, 0.5))
mySketch = myModel.Sketch(name='partitionProfile', 
    sheetSize=40.04, gridSpacing=1.0, transform=t)
g, v, d = mySketch.geometry, mySketch.vertices, mySketch.dimensions
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
r, r1 = mySketch.referenceGeometry, mySketch.referenceVertices
mySketch.sketchOptions.setValues(gridOrigin=(-10.0, -0.5))
mySketch.Line(point1=(0.0, -0.5), point2=(0.0, 0.5))
f1 = myBlockInstance.faces.findAt((10,0,0.5))
e1 = myBlockInstance.edges.findAt((20,0,0.5))
myAssembly.PartitionFaceBySketch(sketchUpEdge=e1, faces=f1, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['partitionProfile']

# Partition the faces to create the inner and outer
# partitions for sweep meshing

f1 = myBlockInstance.faces.findAt((10,25,1))
e1 = myBlockInstance.edges.findAt((20,25,1))
t = myAssembly.MakeSketchTransform(sketchPlane=f1,
    sketchUpEdge=e1, sketchPlaneSide=SIDE1, origin=(10.0, 25.0, 1.0))
mySketch = myModel.Sketch(name='sweepPartitionProfile', 
    sheetSize=107.72, gridSpacing=2.69, transform=t)
g, v, d = mySketch.geometry, mySketch.vertices, mySketch.dimensions
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
r, r1 = mySketch.referenceGeometry, mySketch.referenceVertices
mySketch.sketchOptions.setValues(gridOrigin=(-10.0,-25.0))
mySketch.ArcByCenterEnds(center=(0.0,-25.0), point1=(0.5,-25.0),
    point2=(-0.5,-25.0), direction=COUNTERCLOCKWISE)
mySketch.ArcByCenterEnds(center=(0.0,-25.0), point1=(5.0,-25.0),
    point2=(-5.0,-25.0), direction=COUNTERCLOCKWISE)
pickedFaces = myBlockInstance.faces.findAt((10,25,1))
e1 = myBlockInstance.edges.findAt((20,25,1))
myAssembly.PartitionFaceBySketch(sketchUpEdge=e1, faces=pickedFaces,
    sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['sweepPartitionProfile']

# Extrude the partitions created above to partition the whole part

pickedCells = myBlockInstance.cells
e1 = myBlockInstance.edges.findAt((10.,0.500,1.))
pickedEdges =(e1, )
e2 = myBlockInstance.edges.findAt((10.,0.,0.25))
myAssembly.PartitionCellByExtrudeEdge(line=e2,
    cells=pickedCells, edges=pickedEdges, sense=REVERSE)

pickedCells = myBlockInstance.cells
e1 = myBlockInstance.edges.findAt((10,5,1))
pickedEdges =(e1, )
e2 = myBlockInstance.edges.findAt((10.,0.,0.25))
myAssembly.PartitionCellByExtrudeEdge(line=e2, cells=pickedCells,
    edges=pickedEdges, sense=REVERSE)

# Create a set for the X face

faces1 = myBlockInstance.faces.findAt(((9.75,0.,0.500),),
    ((7.25,0.,0.500),), ((2.5,0.,0.500),),)
myAssembly.Set(faces=faces1, name='xFace')

# Create a set for the Y face

faces1 = myBlockInstance.faces.findAt(((0.,25.,0.500),),)
myAssembly.Set(faces=faces1, name='yFace')

# Create a set for the crack tip

edges1 = myBlockInstance.edges.findAt(((10.,0.,0.500),),)
myAssembly.Set(edges=edges1, name='crackTip')

# Create a set for the front Z face

faces1 = myBlockInstance.faces.findAt(((10.,0.250,1.),), ((10.,2.75,1.),),
    ((10.,27.5,1.),),)
myAssembly.Set(faces=faces1, name='zFrontFace')

# Create a set for the rear Z face

faces1 = myBlockInstance.faces.findAt(((10.,27.5,0.),), ((10.,2.75,0.),),
    ((10.,0.250,0.),),)
myAssembly.Set(faces=faces1, name='zRearFace')

# Create a surface set for the loading surface

side1Faces1 = myBlockInstance.faces.findAt(((10,50,0.5),),)
myAssembly.Surface(side1Faces=side1Faces1, name='loadSurf')

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

# Create a step for applying a load

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

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

# Create interaction properties

# Create the contour integral definition for the crack

crackFront = crackTip = myAssembly.sets['crackTip']
verts = myBlockInstance.vertices
v1 = myBlockInstance.vertices.findAt((10,0,1))
v2 = myBlockInstance.vertices.findAt((0,0,1))
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, distribution=UNIFORM, localCsys=None)

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

region = myAssembly.sets['zFrontFace']
myModel.DisplacementBC(name='zFrontFixed', createStepName='Initial',
    region=region, u3=SET, distribution=UNIFORM, 
    localCsys=None)

region = myAssembly.sets['zRearFace']
myModel.DisplacementBC(name='zRearFixed', createStepName='Initial',
    region=region, u3=SET, distribution=UNIFORM, 
    localCsys=None)

# Assign the loads

region = myAssembly.surfaces['loadSurf']
myModel.Pressure(name='prLoad', createStepName='ApplyLoad',
    region=region, distribution=UNIFORM, magnitude=-100.0)

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

# Create a mesh 

# Assign meshing controls

pickedRegions = myBlockInstance.cells.findAt(((10,0.25,0.5),),)
myAssembly.setMeshControls(regions=pickedRegions, elemShape=WEDGE)

pickedRegions = myBlockInstance.cells.findAt(((10.,2.75,1.),),
    ((10.,27.5,1.),),)
myAssembly.setMeshControls(regions=pickedRegions, technique=SWEEP)

# Seed all the edges

pickedEdges = myBlockInstance.edges.findAt(((10.25,0.,0.),),
    ((9.75,0.,0.),), ((10.25,0.,1.),), ((9.75,0.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1,
    constraint=FIXED)

pickedEdges = myBlockInstance.edges.findAt(((10.5,0.,0.500),),
    ((9.5,0.,0.500),), ((10.,0.,0.500),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1,
    constraint=FIXED)

pickedEdges = myBlockInstance.edges.findAt(((10.,0.500,0.),),
    ((10.,0.500,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=12,
    constraint=FIXED)

pickedEdges1 = myBlockInstance.edges.findAt(((8.5,0,0),),)
pickedEdges2 = myBlockInstance.edges.findAt(((8.5,0,1),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1,
    end2Edges=pickedEdges2, ratio=2.0, number=5,
    constraint=FIXED)

pickedEdges1 = myBlockInstance.edges.findAt(((10.75,0,0),),)
pickedEdges2 = myBlockInstance.edges.findAt(((10.75,0,1),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1,
    end2Edges=pickedEdges2, ratio=2.0, number=5, constraint=FIXED)

pickedEdges = myBlockInstance.edges.findAt(((0,0,0.5),),
    ((20,0,0.5),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FIXED)

pickedEdges = myBlockInstance.edges.findAt(((10.,5.,0.),),
    ((10.,5.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=12)

pickedEdges = myBlockInstance.edges.findAt(((2.5,0.,0.),),
    ((2.5,0.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=3)

pickedEdges = myBlockInstance.edges.findAt(((17.5,0.,0.),),
    ((17.5,0.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4)

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

pickedEdges = myBlockInstance.edges.findAt(((10.,50.,0.),),
    ((10.,50.,1.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=5,
    constraint=FIXED)

pickedEdges = myBlockInstance.edges.findAt(((0.,50.,0.500),),
    ((20.,50.,0.500),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=1,
    constraint=FIXED)

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 = myBlockInstance.cells
pickedRegions =(cells1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2, elemType3))
partInstances =(myBlockInstance, )
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)

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

# Create the job 

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

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