'''
-----------------------------------------------------------------------------
 Two dimensional model of a slant crack in a plate (lamina elasticity)
 modeled using plane stress elements (CPS8).
-----------------------------------------------------------------------------
'''

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 = 'SlantCrackLamCPS8'
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=50.0)
mySketch.sketchOptions.setValues(viewStyle=REGULAR)
mySketch.setPrimaryObject(option=STANDALONE)

mySketch.rectangle(point1=(-2.5, -5.0), point2=(2.5, 5.0))

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

myViewport.setValues(displayedObject=myPlate)

# Create a partition for the seam cracks and the boundary conditions

face = myPlate.faces.findAt((0,0,0),)
t = myPlate.MakeSketchTransform(sketchPlane=face,
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', 
    sheetSize=22.36, gridSpacing=0.55, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myPlate.projectReferencesOntoSketch(sketch=mySketch, filter=COPLANAR_EDGES)
mySketch.AngularConstructionLine(point=(0.0, 0.0), angle=45.0)
mySketch.Line(point1=(0.0, 0.0), point2=(1.1, 1.1))
mySketch.Line(point1=(0.0, 0.0), point2=(-1.1, -1.1))
v = mySketch.vertices
d = mySketch.dimensions
mySketch.ObliqueDimension(vertex1=v[4], vertex2=v[5],
    textPoint=(-0.512194991111755, 1.4282363653183))
mySketch.ObliqueDimension(vertex1=v[4], vertex2=v[7],
    textPoint=(0.0590995773673058, -1.84193241596222))
mySketch.changeDimension(dimension=d[9], value=1.0,
    vertexList=(v[5], ))
mySketch.changeDimension(dimension=d[10], value=1.0,
    vertexList=(v[7], ))
faces = myPlate.faces
pickedFaces = myPlate.faces.findAt((0,0,0))
pickedFaces = faces[pickedFaces.index:(pickedFaces.index+1)]
myPlate.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

# Create a set referring to the whole part

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

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

# Assign material properties

# Create a local coordinate system to define the material orientation

myPlate.DatumCsysByThreePoints(name='Rect', coordSysType=CARTESIAN,
    origin=(0.0, 0.0, 0.0), point1=(1.0, 0.0, 0.0),
    point2=(0.0, 1.0, 0.0))

# Create laminate elasticity

myModel.Material(name='Lam')
myModel.materials['Lam'].Elastic(type=LAMINA,
    table=((200000000000.0, 150000000000.0, 0.35, 20000000000.0, 
    9000000000.0, 8500000000.0), ))

myModel.HomogeneousSolidSection(name='SolidHomogeneous',
    material='Lam', thickness=1.0)

region = myPlate.sets['All']

# Assign the material orientation to the part

datums = myPlate.datums[4]
myPlate.assignMaterialOrientation(region=region, localCsys=datums)

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

# Create a set for the top edge of the plate

edges = myPlateInstance.edges
edge1 = myPlateInstance.edges.findAt((0,5,0))
edge1 = edges[edge1.index:(edge1.index+1)]
myAssembly.Set(edges=edge1, name='topEdge')

# Create a set for the bottom edge of the plate

edges = myPlateInstance.edges
edge1 = myPlateInstance.edges.findAt((0,-5,0))
edge1 = edges[edge1.index:(edge1.index+1)]
myAssembly.Set(edges=edge1, name='bottomEdge')

# Create a set for the top crack tip

verts1 = myPlateInstance.vertices
ver1 = myPlateInstance.vertices.findAt((0.707107,0.707107,0))
ver1 = verts1[ver1.index:(ver1.index+1)]
myAssembly.Set(vertices=ver1, name='topCrackTip')

# Create a set for the bottom crack tip

verts1 = myPlateInstance.vertices
ver1 = myPlateInstance.vertices.findAt((-0.707107,-0.707107,0))
ver1 = verts1[ver1.index:(ver1.index+1)]
myAssembly.Set(vertices=ver1, name='bottomCrackTip')

# Partition the top and bottom edges at their midpoints

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((0,5,0))
mPt1 = myPlateInstance.InterestingPoint(edge=e1, rule=MIDDLE)
myAssembly.PartitionEdgeByPoint(edge=e1, point=mPt1)
e1 = myPlateInstance.edges.findAt((0,-5,0))
mPt2 = myPlateInstance.InterestingPoint(edge=e1, rule=MIDDLE)
myAssembly.PartitionEdgeByPoint(edge=e1, point=mPt2)

# Create a set for the top and bottom points to fix in X

verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((0,5,0))
v2 = myPlateInstance.vertices.findAt((0,-5,0))
verts1 = (verts[v1.index:(v1.index+1)] +
    verts[v2.index:(v2.index+1)])
myAssembly.Set(vertices=verts1, name='fixedInX')

# Create circular partitions for sweep mesh around the two crack tips

faces1 = myPlateInstance.faces.findAt((0,0,0),)
t = myAssembly.MakeSketchTransform(sketchPlane=faces1,
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
mySketch = myModel.Sketch(name='plateProfile', sheetSize=22.36,
    gridSpacing=0.55, transform=t)
mySketch.setPrimaryObject(option=SUPERIMPOSE)
myAssembly.projectReferencesOntoSketch(sketch=mySketch,
    filter=COPLANAR_EDGES)
mySketch.CircleByCenterPerimeter(center=(0.707107, 0.707107),
    point1=(0.353553, 0.353553))
mySketch.CircleByCenterPerimeter(center=(-0.707107, -0.707107), 
    point1=(-0.353553, -0.353553))
g = mySketch.geometry
mySketch.RadialDimension(curve=g[2],
    textPoint=(0.16, 1.92))
mySketch.RadialDimension(curve=g[4],
    textPoint=(0.08, -2.33))
r = mySketch.referenceVertices
mySketch.constraintReferences(vertex1=r.findAt((-0.707107, -0.707107), 1))
faces = myPlateInstance.faces
pickedFaces = myPlateInstance.faces.findAt((0,2.5,0))
pickedFaces = faces[pickedFaces.index:(pickedFaces.index+1)]
myAssembly.PartitionFaceBySketch(faces=pickedFaces, sketch=mySketch)
mySketch.unsetPrimaryObject()
del myModel.sketches['plateProfile']

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

# Create a step for applying a load

myModel.StaticStep(name='Step-1', previous='Initial',
    description='Apply a displacement of 0.01 on the top and bottom edges')
myModel.StaticStep(name='Step-2', previous='Step-1', 
    description='Apply a displacement of 0.01 on the top and bottom edges')
myModel.StaticStep(name='Step-3', previous='Step-2', 
    description='Apply a displacement of 0.01 on the top and bottom edges')

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

# Create interaction properties

# Create a set for the seam crack edge

edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((-0.53033, -0.53033, 0))
e2 = myPlateInstance.edges.findAt((-0.176776, -0.176776, 0))
e3 = myPlateInstance.edges.findAt((0.176776, 0.176776, 0))
e4 = myPlateInstance.edges.findAt((0.53033, 0.53033, 0))

edges1 = (edges[e1.index:(e1.index+1)]
    + edges[e2.index:(e2.index+1)] + edges[e3.index:(e3.index+1)]
    + edges[e4.index:(e4.index+1)])
myAssembly.Set(edges=edges1, name='seamCrack')

# Assign seam crack property to the entire crack

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

# Create the contour integral definition for the top crack

topCrackFront = topCrackTip = myAssembly.sets['topCrackTip']
verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((0.707107,0.707107,0),)
edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((1.06066,1.06066,0))
v2 = myPlateInstance.InterestingPoint(edge=e1, rule=MIDDLE)
myAssembly.engineeringFeatures.ContourIntegral(name='TopCrack', symmetric=OFF,
    crackFront=topCrackFront, crackTip=topCrackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((v1,v2),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)

# Create the contour integral definition for the bottom crack

bottomCrackFront = bottomCrackTip = myAssembly.sets['bottomCrackTip']
verts = myPlateInstance.vertices
v1 = myPlateInstance.vertices.findAt((-0.707107,-0.707107,0),)
edges = myPlateInstance.edges
e1 = myPlateInstance.edges.findAt((-1.06066,-1.06066,0))
v2 = myPlateInstance.InterestingPoint(edge=e1, rule=MIDDLE)
myAssembly.engineeringFeatures.ContourIntegral(name='BottomCrack',
    symmetric=OFF, crackFront=bottomCrackFront,
    crackTip=bottomCrackTip, extensionDirectionMethod=Q_VECTORS,
    qVectors=((v1,v2),), midNodePosition=0.25,
    collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions 

# Assign boundary conditions

region = myAssembly.sets['fixedInX']
myModel.DisplacementBC(name='XFixed', createStepName='Initial',
    region=region, u1=SET, distribution=UNIFORM, localCsys=None)

region = myAssembly.sets['topEdge']
myModel.DisplacementBC(name='YDispTopEdge', createStepName='Step-1',
    region=region, u2=0.01, fixed=OFF, distribution=UNIFORM,
    localCsys=None)

region = myAssembly.sets['bottomEdge']
myModel.DisplacementBC(name='YDispBottomEdge', createStepName='Step-1',
    region=region, u2=-0.01, amplitude=UNSET, fixed=OFF,
    distribution=UNIFORM, localCsys=None)

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

# Create a mesh 

# Seed all the edges

e1 = myPlateInstance.edges
pickedEdges1 = myPlateInstance.edges.findAt((0.707107,0.707107,0),)
pickedEdges2 = myPlateInstance.edges.findAt((-0.707107,-0.707107,0),)
pickedEdges1 = e1[pickedEdges1.index:(pickedEdges1.index+1)]
pickedEdges2 = e1[pickedEdges2.index:(pickedEdges2.index+1)]

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

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

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

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

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

# Assign meshing controls to the respective regions

faces = myPlateInstance.faces
f1 = myPlateInstance.faces.findAt((-0.707107,-0.707107,0))
f2 = myPlateInstance.faces.findAt((0.707107,0.707107,0))
pickedRegions = (faces[f1.index:(f1.index+1)] +
    faces[f2.index:(f2.index+1)])
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED,
    technique=SWEEP)
elemType1 = mesh.ElemType(elemCode=CPS8, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPS6M, elemLibrary=STANDARD)
faces1 = myPlateInstance.faces
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='JIntBottom')
myModel.historyOutputRequests['JIntBottom'].setValues(contourIntegral='BottomCrack', numberOfContours=7)
myModel.HistoryOutputRequest(name='JIntTop', createStepName='Step-1',
    contourIntegral='TopCrack', numberOfContours=7)
myModel.HistoryOutputRequest(name='StressIntBottom', createStepName='Step-1',
    contourIntegral='BottomCrack', numberOfContours=7, 
    contourType=K_FACTORS)
myModel.HistoryOutputRequest(name='StressIntTop', createStepName='Step-1',
    contourIntegral='TopCrack', numberOfContours=7, contourType=K_FACTORS)

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

# Create the job 

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

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