'''
-----------------------------------------------------------------------------
 Two dimensional model of a single edged notch specimen with crack
 along interface between two dissimilar elastic materials modeled
 using plane strain elements (CPE8).
-----------------------------------------------------------------------------
'''

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 = '2DDissimilarMaterialsCPE8'
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=(0.0, -40.0), point2=(20.0, 40.0))

myPlate = myModel.Part(name='Plate',
    dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
myPlate.BaseShell(sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a partition for the seam cracks

face = myPlate.faces.findAt((10,0,0),)
t = myPlate.MakeSketchTransform(sketchPlane=face,
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', 
    sheetSize=164.92, gridSpacing=4.12, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myPlate.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
mySketch.Line(point1=(0.0, 0.0), point2=(20.0, 0.0))
pickedFaces = myPlate.faces.findAt((10,0,0),)
myPlate.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Partition the centre edge

pickedEdge = myPlate.edges.findAt((10,0,0),)
myPlate.PartitionEdgeByPoint(edge=pickedEdge,
    point=myPlate.InterestingPoint(edge=pickedEdge, rule=MIDDLE))

# Create a set referring to the whole part

faces = myPlate.faces.findAt(((10,10,0),), ((10,-10,0),))
myPlate.Set(faces=faces, name='All')

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

# Assign material properties

# Create linear elastic material for the top section

myModel.Material(name='A1')
myModel.materials['A1'].Elastic(table=((1000000.0, 0.3), ))

# Create linear elastic material for the bottom section

myModel.Material(name='A2')
myModel.materials['A2'].Elastic(table=((2000000.0, 0.2), ))

# Create sets for the top and bottom faces

faces = myPlate.faces
pickedFace = myPlate.faces.findAt((10,10,0),)
face1 = faces[pickedFace.index:(pickedFace.index+1)]
myPlate.Set(faces=face1, name='Top')

pickedFace = myPlate.faces.findAt((10,-10,0),)
face1 = faces[pickedFace.index:(pickedFace.index+1)]
myPlate.Set(faces=face1, name='Bottom')

# Create the sections with the above material properties

myModel.HomogeneousSolidSection(name='SolidHomogeneousTop',
    material='A1', thickness=1.0)
myModel.HomogeneousSolidSection(name='SolidHomogeneousBottom',
    material='A2', thickness=1.0)

# Assign the above section to the faces

region = myPlate.sets['Top']
myPlate.SectionAssignment(region=region,
    sectionName='SolidHomogeneousTop')
region = myPlate.sets['Bottom']
myPlate.SectionAssignment(region=region,
    sectionName='SolidHomogeneousBottom')

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

# Create an assembly

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

# Create circular partitions for sweep mesh around the crack tip

faces1 = myPlateInstance.faces.findAt((10,10,0),)
t = myAssembly.MakeSketchTransform(sketchPlane=faces1, sketchPlaneSide=SIDE1,
    origin=(10.0, 20.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', sheetSize=89.44,
    gridSpacing=2.23, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
r1 = mySketch.referenceVertices
mySketch.sketchOptions.setValues(gridOrigin=(0.0, -20.0))
mySketch.CircleByCenterPerimeter(center=(0.0, -20.0), point1=(5.0, -20.0))
mySketch.constraintReferences(vertex1=r1.findAt((0.0, -20.0), 1))
faces = myPlateInstance.faces
f1 = myPlateInstance.faces.findAt((10,10,0),)
f2 = myPlateInstance.faces.findAt((10,-10,0),)
pickedFaces = faces[f1.index:(f1.index+1)] + faces[f2.index:(f2.index+1)]
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a set for the fixed point

vert = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((0,0,0),)
verts1 = vert[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='fixedPt')

# Create a set for the crack tip

v1 = myPlateInstance.vertices.findAt((10,0,0),)
verts1 = vert[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='crackTip')

# Create a set for the seam crack edge

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((12.5,0,0),)
e2 = myPlateInstance.edges.findAt((17.5,0,0),)
edges1 = edges[e1.index:(e1.index+1)]+edges[e2.index:(e2.index+1)]
myAssembly.Set(edges=edges1, name='seamCrack')

# Create a set for the top surface to apply the load

e1 = myPlateInstance.edges.findAt((10,40,0),)
side1Edges1 = edges[e1.index:(e1.index+1)]
myAssembly.Surface(side1Edges=side1Edges1, name='topSurf')

# Create a set for the bottom surface to apply the load

e1 = myPlateInstance.edges.findAt((10,-40,0),)
side1Edges1 = edges[e1.index:(e1.index+1)]
myAssembly.Surface(side1Edges=side1Edges1, name='bottomSurf')

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

# Create a step for applying a load

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

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

# Create interaction properties

# Assign seam crack property to the crack

pickedRegions = myAssembly.sets['seamCrack']
myAssembly.engineeringFeatures.assignSeam(regions=pickedRegions)

crackFront = crackTip = myAssembly.sets['crackTip']
vert = myPlateInstance.vertices
v21 = myPlateInstance.vertices.findAt((10,0,0))
v11 = myPlateInstance.vertices.findAt((20,0,0))
myAssembly.engineeringFeatures.ContourIntegral(name='Crack',
    symmetric=OFF, crackFront=crackFront, crackTip=crackTip, 
    extensionDirectionMethod=Q_VECTORS, qVectors=((v11, v21),), 
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

fPt = myAssembly.sets['fixedPt']
myModel.DisplacementBC(name='FixedXandY', createStepName='Initial',
    region=fPt, u1=0, u2=0, distributionType=UNIFORM,
    localCsys=None)

# Assign load conditions

region = myAssembly.surfaces['topSurf']
myModel.Pressure(name='topLoad', createStepName='LoadPlate',
    region=region, distributionType=UNIFORM, magnitude=-1.0)
region = myAssembly.surfaces['bottomSurf']
myModel.Pressure(name='bottomSurf', createStepName='LoadPlate',
    region=region, distributionType=UNIFORM, magnitude=-1.0)

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

# Create a mesh 

# Seed all the edges

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((12.5,0,0))
e2 = myPlateInstance.edges.findAt((7.5,0,0))
pickedEdges1 = edges[e1.index:(e1.index+1)]
pickedEdges2 = edges[e2.index:(e2.index+1)]

myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, end2Edges=pickedEdges2,
    ratio=2.0, number=6, constraint=FIXED)

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

pickedEdges = myPlateInstance.edges.findAt(((0,20,0),), ((0,-20,0),),
    ((20,20,0),), ((20,-20,0),))
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=15, constraint=FIXED)

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

pickedEdges = myPlateInstance.edges.findAt(((2.5,0,0),))
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FIXED)

pickedEdges = myPlateInstance.edges.findAt(((17.5,0,0),))
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=3, constraint=FIXED)

# Assign meshing controls to the respective regions
faces1 = myPlateInstance.faces
pickedRegions = faces1.getSequenceFromMask(mask=('[#f ]', ), )
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD) 
pickedFaces = myPlateInstance.faces.findAt(((10,2.5,0),), ((10,-2.5,0),))
myAssembly.setMeshControls(regions=pickedFaces, elemShape=QUAD_DOMINATED, 
    technique=SWEEP)
elemType1 = mesh.ElemType(elemCode=CPE8, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPE6M, elemLibrary=STANDARD)
pickedRegions =(faces1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2))
partInstances =(myPlateInstance, )
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=7)
myModel.HistoryOutputRequest(name='StressInt-Factors',
    createStepName='LoadPlate', contourIntegral='Crack',
    numberOfContours=7, contourType=K_FACTORS)

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

# Create the job 

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

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