'''
-----------------------------------------------------------------------------
 Two dimensional sub model of a double edged notch specimen (linear
 elastic material) modeled using plane strain elements (CPE8).

 First the global model job is completed. The *.odb file from the global
 model is used to drive this submodel.

 Global model scripts to be run:
 2DDoubleEdSymmGlCPE8_model.py and 2DDoubleEdSymmGlCPE8_job.py
-----------------------------------------------------------------------------
'''

from abaqus import *
from abaqusConstants import *

import part, material, section, assembly, step, interaction
import regionToolset, displayGroupMdbToolset as dgm, mesh, load, job 

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

# Copy the global model into a new model

Mdb()
globalModelName = '2DDoubleEdSymmGlCPE8'
openMdb(globalModelName)
subModelName = '2DDoubleEdSymmSubCPE8'
mySubModel = mdb.Model(name=subModelName)

mySubModel = (mdb.Model(name=subModelName,
    objectToCopy=mdb.models[globalModelName]))
    
# Create a new viewport in which to display the model
# and the results of the analysis.

myAssembly = mySubModel.rootAssembly
myViewport = session.viewports['Viewport: 1']
myViewport.setValues(displayedObject=myAssembly)
myViewport.makeCurrent()
myViewport.maximize()

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

# Edit the model attributes to refer to the global model ODB file
# to drive the sub model.

mySubModel.setValues(description='Double edged notch specimen submodel', 
    globalJob='2DDoubleEdSymmGlCPE8.odb')

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

# Create the sub model section of the global model

# Create a sketch for the base feature

myPlate = mySubModel.parts['Plate']
myViewport.setValues(displayedObject=myPlate)

mySubSketch = mySubModel.Sketch(name='plateProfile',sheetSize=200.0)
mySubSketch.sketchOptions.setValues(viewStyle=REGULAR)
mySubSketch.setPrimaryObject(option=STANDALONE)

mySubSketch.ArcByCenterEnds(center=(0.0, 0.0), point1=(5.0, 0.0),
    point2=(-5.0, 0.0), direction=COUNTERCLOCKWISE)
mySubSketch.Line(point1=(-5.0, 0.0), point2=(5.0, 0.0))

mySubPlate = mySubModel.Part(name='subPlate', 
    dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
mySubPlate.BaseShell(sketch=mySubSketch)
mySubSketch.unsetPrimaryObject()
del mySubModel.sketches['plateProfile']

myViewport.setValues(displayedObject=mySubPlate)

# Partition the edge to create the crack tip

e1 = mySubPlate.edges.findAt((0,0,0))
mySubPlate.PartitionEdgeByPoint(edge=e1,
    point=mySubPlate.InterestingPoint(edge=e1, rule=MIDDLE))

# Create a set referring to the whole part

faces = mySubPlate.faces.findAt(((0,2.5,0),))
mySubPlate.Set(faces=faces, name='subAll')

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

# Assign material properties

region = mySubPlate.sets['subAll']
mySubPlate.SectionAssignment(region=region,
    sectionName='SolidHomogeneous')

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

# Create an assembly
# Place the plate created above at the same position as in the
# global model. Then the instance of the full plate is deleted.

myViewport.setValues(displayedObject=myAssembly)
myAssembly.Instance(name='subPlate-1', part=mySubPlate, dependent=OFF)
mySubPlateInstance = myAssembly.instances['subPlate-1']
mySubPlateInstance.translate(vector=(10.0, -3.06151588455594e-16, 0.0))
del myAssembly.features['myPlate-1']

# Create a set for the crack tip

verts = mySubPlateInstance.vertices
v1 = mySubPlateInstance.vertices.findAt((10,0,0))
verts1 = verts[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='crackTipSub')

# Create a set for the X edge of the plate

edges = mySubPlateInstance.edges
e1 = mySubPlateInstance.edges.findAt((7.5,0,0))
e1 = edges[e1.index:(e1.index+1)]
myAssembly.Set(edges=e1, name='xEdgeSub')

# Create a set for the sub modeling boundary

e1 = mySubPlateInstance.edges.findAt((10,5,0))
e1 = edges[e1.index:(e1.index+1)]
myAssembly.Set(edges=e1, name='subBoundary')

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

# Create interaction properties

# Create the contour integral definition for the crack

crackFront = crackTip = myAssembly.sets['crackTipSub']
verts = mySubPlateInstance.vertices
v1 = mySubPlateInstance.vertices.findAt((10,0,0))
v2 = mySubPlateInstance.vertices.findAt((5,0,0))
myAssembly.engineeringFeatures.ContourIntegral(name='CrackSub',
    symmetric=ON, crackFront=crackFront, crackTip=crackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((v1,v2),),
    midNodePosition=0.25, collapsedElementAtTip=SINGLE_NODE)
del myAssembly.engineeringFeatures.cracks['Crack']

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

# Create loads and boundary conditions 

# Delete the existing boundary conditions

del mySubModel.boundaryConditions['xFixed']
del mySubModel.boundaryConditions['yFixed']
del mySubModel.loads['topLoad']
del myAssembly.sets['xEdge']
del myAssembly.sets['yEdge']
del myAssembly.sets['crackTip']

# Create the sub modeling boundary condition

region = myAssembly.sets['subBoundary']
mySubModel.SubmodelBC(name='subBC', createStepName='ApplyLoad',
    region=region, globalStep='1', globalIncrement=0,
    timeScale=OFF, dof=(1, 2), globalDrivingRegion='', 
    absoluteExteriorTolerance=0.0, exteriorTolerance=0.025)

# Create a boundary condition for X edge to be fixed in Y

region = myAssembly.sets['xEdgeSub']
mySubModel.DisplacementBC(name='yFixed', createStepName='Initial',
    region=region, u2=SET, distribution=UNIFORM, localCsys=None)

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

# Create a mesh 

# Seed all the edges

edges = mySubPlateInstance.edges
pickedEdges1 = mySubPlateInstance.edges.findAt((11,0,0))
pickedEdges2 = mySubPlateInstance.edges.findAt((9,0,0))
pickedEdges1 = edges[pickedEdges1.index:(pickedEdges1.index+1)]
pickedEdges2 = edges[pickedEdges2.index:(pickedEdges2.index+1)]
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1,
    end2Edges=pickedEdges2, ratio=3.0, number=6, constraint=FIXED)

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

# Assign meshing controls to the respective regions

faces = mySubPlateInstance.faces
f1 = mySubPlateInstance.faces.findAt((10,2.5,0))
pickedRegions = faces[f1.index:(f1.index+1)]
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED,
    technique=SWEEP)
elemType1 = mesh.ElemType(elemCode=CPE8, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CPE6M, elemLibrary=STANDARD)
faces1 = mySubPlateInstance.faces
pickedRegions =(faces1, )
myAssembly.setElementType(regions=pickedRegions,
    elemTypes=(elemType1, elemType2))
partInstances =(mySubPlateInstance, )
myAssembly.generateMesh(regions=partInstances)

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

# Request history output for the crack

mySubModel.historyOutputRequests['JInt'].setValues(
    contourIntegral='CrackSub', numberOfContours=6)

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

# Create the job 

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

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