'''
-----------------------------------------------------------------------------
 Axisymmetric sub model of a conical crack in a block modeled using
 axisymmetric reduced integration elements (CAX8R).

 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:
 AxisymmConeCrackGl_model.py and AxisymmConeCrackGl_job.py
-----------------------------------------------------------------------------
'''

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 

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

# Copy the global model into a new model

Mdb()
globalModelName = 'AxisymmConeCrackGl'
openMdb(globalModelName)
subModelName = 'AxisymmConeCrackSub'
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='Submodel of an axisymmetric cone crack',
    globalJob='AxisymmConeCrackGl.odb')

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

# Create the sub model section of the global model

# Create a sketch for the base feature

myBlock = mySubModel.parts['Block']
myViewport.setValues(displayedObject=myBlock)

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

mySubSketch.ObliqueConstructionLine(point1=(0.0, -100.0),
    point2=(0.0, 100.0))
mySubSketch.rectangle(point1=(0.0, 0.0),
    point2=(42.0, -24.0))
mySubBlock = mySubModel.Part(name='subBlock', 
    dimensionality=AXISYMMETRIC, type=DEFORMABLE_BODY)

mySubBlock.BaseShell(sketch=mySubSketch)
mySubSketch.unsetPrimaryObject()
del mySubModel.sketches['blockProfile']

myViewport.setValues(displayedObject=mySubBlock)

# Partition the block to create the crack

face = mySubBlock.faces
t = mySubBlock.MakeSketchTransform(sketchPlane=face[0],
    sketchPlaneSide=SIDE1, origin=(21.0,-12.0,0.0))
mySubSketch = mySubModel.Sketch(name='crackProfile', 
    sheetSize=96.74, gridSpacing=2.41, transform=t)
mySubSketch.setPrimaryObject(option=SUPERIMPOSE)
mySubBlock.projectReferencesOntoSketch(sketch=mySubSketch,
    filter=COPLANAR_EDGES)
r = mySubSketch.referenceGeometry
r1 = mySubSketch.referenceVertices
mySubSketch.sketchOptions.setValues(gridOrigin=(-21.0,12.0))
mySubSketch.AngularConstructionLine(point=(-11.0,12.0), angle=135.0)
mySubSketch.Line(point1=(-11.0,12.0),
    point2=(-0.393398282201787,1.39339828220179))
mySubSketch.Line(point1=(-0.393398282201787,1.39339828220179),
    point2=(-21.0,1.39339828220179))
pickedFaces = mySubBlock.faces
mySubBlock.PartitionFaceBySketch(faces=pickedFaces,
    sketch=mySubSketch)
mySubSketch.unsetPrimaryObject()
del mySubModel.sketches['crackProfile']

# Create the remaining partitions

f = mySubBlock.faces
t = mySubBlock.MakeSketchTransform(sketchPlane=f[0],
    sketchPlaneSide=SIDE1, origin=(23.503219,-13.167747,0.0))
mySubSketch = mySubModel.Sketch(name='partitionProfile', 
    sheetSize=96.74, gridSpacing=2.41, transform=t)
mySubSketch.setPrimaryObject(option=SUPERIMPOSE)
mySubBlock.projectReferencesOntoSketch(sketch=mySubSketch,
    filter=COPLANAR_EDGES)
r = mySubSketch.referenceGeometry
r1 = mySubSketch.referenceVertices
mySubSketch.sketchOptions.setValues(gridOrigin=(-23.503219, 13.167747))
mySubSketch.CircleByCenterPerimeter(
    center=(-2.89661728220179, 2.56114528220179), 
    point1=(-7.89661728220179, 2.56114528220179))
mySubSketch.Line(point1=(-2.89661728220179,2.56114528220179),
    point2=(18.4967810000001,2.56114528220179))
mySubSketch.AngularConstructionLine(
    point=(2.10338271779821,2.56114528220179), 
    angle=45.0)
mySubSketch.AngularConstructionLine(
    point=(-7.89661728220179, 2.56114528220179), 
    angle=60.0)
mySubSketch.AngularConstructionLine(
    point=(2.10338271779823, 2.5611452822018), 
    angle=135.0)
mySubSketch.Line(point1=(-23.503219,-5.442253),
    point2=(18.4967810000001,-5.442253))
mySubSketch.Line(point1=(-7.89661728220179,2.56114528220179),
    point2=(-12.5173814348628,-5.442253))
mySubSketch.Line(point1=(2.10338271779822,2.5611452822018),
    point2=(10.1067810000001,-5.44225300000005))
mySubSketch.Line(point1=(2.10338271779822,2.5611452822018),
    point2=(12.7099844355964,13.167747))
mySubSketch.Line(point1=(-12.5173814348628,-5.442253),
    point2=(-12.5173814348628,-10.832253))
mySubSketch.Line(point1=(10.1067810000001,-5.44225300000005),
    point2=(10.1067810000001,-10.832253))
pickedFaces = mySubBlock.faces
mySubBlock.PartitionFaceBySketch(faces=pickedFaces, sketch=mySubSketch)
mySubSketch.unsetPrimaryObject()
del mySubModel.sketches['partitionProfile']

# Create a set for the whole part

f = mySubBlock.faces
mySubBlock.Set(faces=f, name='subAll')

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

# Assign material properties

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

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

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

myViewport.setValues(displayedObject=myAssembly)
mySubBlockIns = myAssembly.Instance(name='subBlock-1',
    part=mySubBlock, dependent=OFF)
mySubBlockInstance = myAssembly.instances['subBlock-1']
del myAssembly.features['Block-1']

# Delete all the sets/step/BCs/loads/features
# for the full block
del mySubModel.steps['ApplyLoad']
del myAssembly.sets['axisymmEdge']
del myAssembly.sets['base']
del myAssembly.sets['crackTip']
del myAssembly.sets['seamCrackEdge']
del myAssembly.surfaces['topCrackSurf']
del myAssembly.engineeringFeatures.cracks['Crack']
del mySubModel.boundaryConditions['axisymmEdge']
del mySubModel.boundaryConditions['baseFixed']

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

# Create a step

mySubModel.StaticStep(name='subApplyLoad', previous='Initial',
    description='Apply the load on the submodel')

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

# Create the necessary sets

# Create a set for the crack tip

verts = mySubBlockIns.vertices
v1 = mySubBlockIns.vertices.findAt((20.606602,-10.606602,0.))
verts1 = verts[v1.index:(v1.index+1)]
myAssembly.Set(vertices=verts1, name='subCrackTip')

# Create a set for the axisymmetric edge

edges1 = mySubBlockIns.edges.findAt(((0.,-21.305,0.),),
    ((0.,-14.608301,0.),), ((0.,-5.303301,0.),),)
myAssembly.Set(edges=edges1, name='subSymmEdge')

# Create a set for the driven boundary

edges1 = mySubBlockIns.edges.findAt(((42.,-5.303301,0.),),
    ((37.805,-24.,0.),), ((42.,-21.305,0.),),
    ((5.492919,-24.,0.),), ((42.,-14.608301,0.),),
    ((22.297919,-24.,0.),),)
myAssembly.Set(edges=edges1, name='subDrivenBoundary')

# Create a set for the seam edges

edges1 = mySubBlockIns.edges.findAt(((13.535534,-3.535534,0.),),
    ((18.838835,-8.838835,0.),),)
myAssembly.Set(edges=edges1, name='subSeamEdge')

# Create a surface set for the loading surface

s1 = mySubBlockIns.edges.findAt((5.,0.,0.))
edges = mySubBlockIns.edges
side1Edges1 = edges[s1.index:(s1.index+1)]
myAssembly.Surface(side1Edges=side1Edges1, name='subLoadSurf')

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

# Create interaction properties

# Assign seam properties

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

# Create the contour integral definition for the crack

crackFront = crackTip = myAssembly.sets['subCrackTip']
verts = mySubBlockIns.vertices
v1 = mySubBlockIns.vertices.findAt((10,0,0))
v2 = mySubBlockIns.vertices.findAt((20.606602,-10.606602,0.))
myAssembly.engineeringFeatures.ContourIntegral(name='subCrack',
    crackFront=crackFront, crackTip=crackTip,
    extensionDirectionMethod=Q_VECTORS, qVectors=((v1,v2),),
    midNodePosition=0.27, collapsedElementAtTip=SINGLE_NODE)

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

# Create loads and boundary conditions

# Create the BC on the axisymmetric edge. Fix it in X.

region = myAssembly.sets['subSymmEdge']
mySubModel.DisplacementBC(name='subAxisymm', createStepName='Initial',
    region=region, u1=0.0, fixed=OFF,
    distributionType=UNIFORM, localCsys=None)

# Create the sub-modeling BC on the driven boundaries

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

# Apply the pressure load

topSurf = myAssembly.surfaces['subLoadSurf']
mySubModel.Pressure(name='subPrLoad', createStepName='subApplyLoad',
    region=topSurf, distributionType=UNIFORM, magnitude=10.0)

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

# Create a mesh 

# Assign mesh controls

# First assign QUAD mesh type with STRUCTURED meshing technique
# to the whole part and then assign QUAD_DOMINATED SWEEP meshing
# technique to the region around the crack tip

pickedRegions = mySubBlockIns.faces
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD,
    technique=STRUCTURED)

pickedRegions = mySubBlockIns.faces.findAt(((21.56331,-8.296903,0.),),
    ((20.606602,-13.106602,0.),), ((18.296903,-9.649893,0.),))
myAssembly.setMeshControls(regions=pickedRegions, elemShape=QUAD_DOMINATED, 
    technique=SWEEP)

# Seed all the edges

pickedEdges1 = mySubBlockIns.edges.findAt(((20.983011,-10.606602,0.),),
    ((19.891129,-10.606602,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1, ratio=7.0, number=12,
    constraint=FIXED)

pickedEdges2 = mySubBlockIns.edges.findAt(((20.13615,-10.13615,0.),),)
myAssembly.seedEdgeByBias(end2Edges=pickedEdges2, ratio=7.0, number=12,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((22.520019,-5.987204,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=11,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((20.606602,-15.606602,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=16,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((15.987204,-8.693185,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=5,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((23.106602,0.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=11,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((22.297919,-24.,0.),),
    ((22.297919,-18.61,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=16,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((5.,0.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=6,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((7.803301,-10.606602,0.),),
    ((5.492919,-18.61,0.),), ((5.492919,-24.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=10,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((5.,0.,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=6,
    constraint=FIXED)

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

pickedEdges = mySubBlockIns.edges.findAt(((0.,-14.608301,0.),),
    ((29.608301,-14.608301,0.),), ((42.,-14.608301,0.),),
    ((13.29622,-14.608301,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4,
    constraint=FIXED)

pickedEdges = mySubBlockIns.edges.findAt(((13.535534,-3.535534,0.),),
    ((30.909903,-5.303301,0.),), ((42.,-5.303301,0.),),)
myAssembly.seedEdgeByNumber(edges=pickedEdges, number=4,
    constraint=FIXED)

pickedEdges1 = mySubBlockIns.edges.findAt(((0.,-19.693965,0.),),
    ((10.985838,-19.693965,0.),), ((33.61,-19.693965,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1,
    ratio=2.0, number=2, constraint=FIXED)

pickedEdges1 = mySubBlockIns.edges.findAt(((34.515,-18.61,0.),),
    ((34.515,-24.,0.),),)
myAssembly.seedEdgeByBias(end1Edges=pickedEdges1,
    ratio=2.0, number=4, 
    constraint=FIXED)

pickedEdges2 = mySubBlockIns.edges.findAt(((42.,-19.151982,0.),),)
myAssembly.seedEdgeByBias(end2Edges=pickedEdges2,
    ratio=2.0, number=2, constraint=FIXED)

# Assign the element type and generate the mesh

elemType1 = mesh.ElemType(elemCode=CAX8R, elemLibrary=STANDARD, 
    secondOrderAccuracy=OFF, hourglassControl=STIFFNESS, distortionControl=OFF)
elemType2 = mesh.ElemType(elemCode=CAX6M, elemLibrary=STANDARD)
faces1 = mySubBlockIns.faces
pickedRegions =(faces1, )
myAssembly.setElementType(regions=pickedRegions, elemTypes=(elemType1,
    elemType2))
partInstances =(mySubBlockIns, )
myAssembly.generateMesh(regions=partInstances)

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

# Request history output for the crack

mySubModel.historyOutputRequests.changeKey(
    fromName='H-Output-1', toName='JInt')
mySubModel.historyOutputRequests['JInt'].setValues(
    contourIntegral='subCrack', sectionPoints=DEFAULT, rebar=EXCLUDE, 
    numberOfContours=5)
mySubModel.HistoryOutputRequest(name='StrInt', 
    createStepName='subApplyLoad', contourIntegral='subCrack', 
    sectionPoints=DEFAULT, rebar=EXCLUDE, numberOfContours=5, 
    contourType=K_FACTORS)
mySubModel.HistoryOutputRequest(name='TStr', 
    createStepName='subApplyLoad', contourIntegral='subCrack', 
    sectionPoints=DEFAULT, rebar=EXCLUDE, numberOfContours=5, 
    contourType=T_STRESS)

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

# Create the job 

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

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