# build initial CAE model and input file

from abaqus import *
from abaqusConstants import *
session.Viewport(name='Viewport: 1')
session.viewports['Viewport: 1'].makeCurrent()
session.viewports['Viewport: 1'].maximize()
from caeModules import *
from driverUtils import executeOnCaeStartup
executeOnCaeStartup()

s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
    sheetSize=200.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.sketchOptions.setValues(viewStyle=AXISYM)
s.setPrimaryObject(option=STANDALONE)
s.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))
s.FixedConstraint(entity=g[2])
s.rectangle(point1=(42.0, 0.0), point2=(48.0, 21.0))
p = mdb.models['Model-1'].Part(name='topPipe_weld',
    dimensionality=AXISYMMETRIC, type=DEFORMABLE_BODY, twist=ON)
p = mdb.models['Model-1'].parts['topPipe_weld']
p.BaseShell(sketch=s)
s.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['topPipe_weld']
p.ReferencePoint(point=(0.0, 21.0, 0.0))
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
    sheetSize=200.0)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.sketchOptions.setValues(viewStyle=AXISYM)
s1.setPrimaryObject(option=STANDALONE)
s1.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))
s1.FixedConstraint(entity=g[2])
s1.rectangle(point1=(42.0, 0.0), point2=(48.0, -21.0))
p = mdb.models['Model-1'].Part(name='bottomPipe_weld',
    dimensionality=AXISYMMETRIC, type=DEFORMABLE_BODY, twist=ON)
p = mdb.models['Model-1'].parts['bottomPipe_weld']
p.BaseShell(sketch=s1)
s1.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['bottomPipe_weld']
p.ReferencePoint(point=(0.0, -21.0, 0.0))
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']

session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=ON,
    engineeringFeatures=ON)
mdb.models['Model-1'].Material(name='Material-1')
mdb.models['Model-1'].materials['Material-1'].Elastic(table=((1.0, 0.2), ))
mdb.models['Model-1'].Material('astroloy')
mdb.models['Model-1'].materials['astroloy'].Density(table=((7.8e-09, ), ))
mdb.models['Model-1'].materials['astroloy'].Conductivity(
    temperatureDependency=ON, table=((14.7, 20.0), (28.0, 1200.0)))
mdb.models['Model-1'].materials['astroloy'].InelasticHeatFraction()
mdb.models['Model-1'].materials['astroloy'].SpecificHeat(
    temperatureDependency=ON, table=((350000000.0, 20.0), (510000000.0,
    300.0), (980000000.0, 600.0), (1300000000.0, 1000.0), (680000000.0,
    1200.0)))
mdb.models['Model-1'].materials['astroloy'].Expansion(table=((0.0, ), ))
mdb.models['Model-1'].materials['astroloy'].Elastic(table=((180000.0, 0.3), ))
mdb.models['Model-1'].materials['astroloy'].Plastic(temperatureDependency=ON,
    rate=ON, table=((4829.85, 0.0, 0.0, 20.0), (2360.95, 0.0, 0.0, 476.667),
    (1307.55, 0.0, 0.0, 665.824), (616.726, 0.0, 0.0, 810.97), (200.911, 0.0,
    0.0, 933.333), (46.6828, 0.0, 0.0, 1041.14), (13.1001, 0.0, 0.0, 1138.6), (
    8.17139, 0.0, 0.0, 1228.23), (0.332254, 0.0, 0.0, 1311.65), (1e-06, 0.0,
    0.0, 1390.0), (4829.85, 0.0, 1e-05, 20.0), (2360.95, 0.0, 1e-05, 476.667),
    (1307.55, 0.0, 1e-05, 665.824), (616.726, 0.0, 1e-05, 810.97), (200.911,
    0.0, 1e-05, 933.333), (46.6828, 0.0, 1e-05, 1041.14), (13.1001, 0.0,
    1e-05, 1138.6), (8.17139, 0.0, 1e-05, 1228.23), (0.332254, 0.0, 1e-05,
    1311.65), (1e-06, 0.0, 1e-05, 1390.0), (5072.73, 0.0, 50.0, 20.0), (
    3189.17, 0.0, 50.0, 476.667), (2472.78, 0.0, 50.0, 665.824), (1916.36,
    0.0, 50.0, 810.97), (1268.67, 0.0, 50.0, 933.333), (788.886, 0.0, 50.0,
    1041.14), (421.254, 0.0, 50.0, 1138.6), (262.763, 0.0, 50.0, 1228.23), (
    333.01, 0.0, 50.0, 1311.65), (5.0, 0.0, 50.0, 1390.0), (5083.93, 0.0,
    100.0, 20.0), (3232.56, 0.0, 100.0, 476.667), (2544.61, 0.0, 100.0,
    665.824), (2016.52, 0.0, 100.0, 810.97), (1378.2, 0.0, 100.0, 933.333), (
    895.757, 0.0, 100.0, 1041.14), (492.352, 0.0, 100.0, 1138.6), (307.111,
    0.0, 100.0, 1228.23), (454.267, 0.0, 100.0, 1311.65), (10.0, 0.0, 100.0,
    1390.0), (5235.08, 0.0, 1000000.0, 20.0), (3868.33, 0.0, 1000000.0,
    476.667), (3722.68, 0.0, 1000000.0, 665.824), (3968.33, 0.0, 1000000.0,
    810.97), (4141.94, 0.0, 1000000.0, 933.333), (4845.65, 0.0, 1000000.0,
    1041.14), (3910.89, 0.0, 1000000.0, 1138.6), (2439.47, 0.0, 1000000.0,
    1228.23), (28133.3, 0.0, 1000000.0, 1311.65), (100000.0, 0.0, 1000000.0,
    1390.0)))
mdb.models['Model-1'].HomogeneousSolidSection(name='pipe material',
    material='astroloy', thickness=1.0)


p = mdb.models['Model-1'].parts['bottomPipe_weld']
f = p.faces
faces = f.getSequenceFromMask(mask=('[#1 ]', ), )
region = regionToolset.Region(faces=faces)
p = mdb.models['Model-1'].parts['bottomPipe_weld']
p.SectionAssignment(region=region, sectionName='pipe material', offset=0.0)
p = mdb.models['Model-1'].parts['topPipe_weld']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
p = mdb.models['Model-1'].parts['topPipe_weld']
f = p.faces
faces = f.getSequenceFromMask(mask=('[#1 ]', ), )
region = regionToolset.Region(faces=faces)
p = mdb.models['Model-1'].parts['topPipe_weld']
p.SectionAssignment(region=region, sectionName='pipe material', offset=0.0)


a = mdb.models['Model-1'].rootAssembly
session.viewports['Viewport: 1'].setValues(displayedObject=a)
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByThreePoints(coordSysType=CYLINDRICAL, origin=(0.0, 0.0, 0.0),
    point1=(1.0, 0.0, 0.0), point2=(0.0, 0.0, -1.0))
p = mdb.models['Model-1'].parts['bottomPipe_weld']
a.Instance(name='bottomPipe_weld-1', part=p, dependent=OFF)
a = mdb.models['Model-1'].rootAssembly
p = mdb.models['Model-1'].parts['topPipe_weld']
a.Instance(name='topPipe_weld-1', part=p, dependent=OFF)

a = mdb.models['Model-1'].rootAssembly
r1 = a.instances['topPipe_weld-1'].referencePoints
refPoints1=(r1[2], )
a.Set(referencePoints=refPoints1, name='topFlywheelAttachment')


a = mdb.models['Model-1'].rootAssembly
r1 = a.instances['bottomPipe_weld-1'].referencePoints
refPoints1=(r1[2], )
a.Set(referencePoints=refPoints1, name='bottomFlywheelAttachment')


a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['topPipe_weld-1'].edges
edges1 = e1.getSequenceFromMask(mask=('[#4 ]', ), )
a.Set(edges=edges1, name='topAttachment')
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['bottomPipe_weld-1'].edges
edges1 = e1.getSequenceFromMask(mask=('[#4 ]', ), )
a.Set(edges=edges1, name='bottomAttachment')


a3 = mdb.models['Model-1'].rootAssembly
s1 = a3.instances['topPipe_weld-1'].edges
side1Edges1 = s1.getSequenceFromMask(mask=('[#b ]', ), )
a3.Surface(side1Edges=side1Edges1, name='topWeldSurface')


a3 = mdb.models['Model-1'].rootAssembly
s1 = a3.instances['bottomPipe_weld-1'].edges
side1Edges1 = s1.getSequenceFromMask(mask=('[#b ]', ), )
a3.Surface(side1Edges=side1Edges1, name='bottomWeldSurface')


mdb.models['Model-1'].CoupledTempDisplacementStep(name='Weld step',
    previous='Initial', description="Weld step", response=TRANSIENT,
    creepIntegration=CREEP_OFF, timePeriod=timeRemaining, maxNumInc=10000, stabilization=None,
    timeIncrementationMethod=AUTOMATIC, initialInc=0.02, minInc=2e-05, maxInc=outputFrequency,
    deltmx=300, cetol=None, amplitude=STEP,
    extrapolation=PARABOLIC, matrixStorage=UNSYMMETRIC,nlgeom=ON)
mdb.models['Model-1'].steps['Weld step'].Restart(frequency=1,
    numberIntervals=0, overlay=ON, timeMarks=OFF)
mdb.models['Model-1'].FieldOutputRequest(name='F-Output-1',
    createStepName='Weld step', variables=('S', 'U', 'NT', 'PEEQ'))
mdb.models['Model-1'].HistoryOutputRequest(name='H-Output-1',
    createStepName='Weld step', variables=('ALLWK', 'ALLKE' ))
regionDef=mdb.models['Model-1'].rootAssembly.sets['topFlywheelAttachment']
mdb.models['Model-1'].HistoryOutputRequest(name='upset',
    createStepName='Weld step', variables=('U2', ), region=regionDef,
    sectionPoints=DEFAULT, rebar=EXCLUDE)

mdb.saveAs(primaryJobName)

mdb.models['Model-1'].ContactProperty('weld contact')
mdb.models['Model-1'].interactionProperties['weld contact'].NormalBehavior(
    pressureOverclosure=TABULAR, table=((0.0, 0.0), (weldContactPressure, weldContactDistance)),
    constraintEnforcementMethod=DEFAULT)
mdb.models['Model-1'].interactionProperties['weld contact'].TangentialBehavior(
    formulation=USER_DEFINED, nStateDependentVars=0, useProperties=ON, table=((
    360., ), ))

mdb.models['Model-1'].ContactProperty('self contact')
mdb.models['Model-1'].interactionProperties['self contact'].NormalBehavior(
    pressureOverclosure=TABULAR, table=((0.0, -selfContactDistance), (selfContactPressure, 0.0)),
    constraintEnforcementMethod=DEFAULT)

a4 = mdb.models['Model-1'].rootAssembly
region1=a4.surfaces['bottomWeldSurface']
a4 = mdb.models['Model-1'].rootAssembly
region2=a4.surfaces['topWeldSurface']
mdb.models['Model-1'].SurfaceToSurfaceContactStd(name='weld contact 1',
    createStepName='Weld step', master=region1, slave=region2,
    sliding=FINITE,
    enforcement=NODE_TO_SURFACE,
    interactionProperty='weld contact', adjustMethod=NONE, smooth=0.2,
    initialClearance=OMIT, datumAxis=None, clearanceRegion=None)
mdb.models['Model-1'].SurfaceToSurfaceContactStd(name='weld contact 2',
    createStepName='Weld step', master=region2, slave=region1,
    sliding=FINITE,
    enforcement=NODE_TO_SURFACE,
    interactionProperty='weld contact', adjustMethod=NONE, smooth=0.2,
    initialClearance=OMIT, datumAxis=None, clearanceRegion=None)

a5 = mdb.models['Model-1'].rootAssembly
region=a5.surfaces['bottomWeldSurface']
mdb.models['Model-1'].SelfContactStd(name='Bottom self contact', createStepName='Initial',
    enforcement=NODE_TO_SURFACE,
    surface=region, interactionProperty='self contact', smooth=0.2)
region=a5.surfaces['topWeldSurface']
mdb.models['Model-1'].SelfContactStd(name='Top self contact', createStepName='Initial',
    enforcement=NODE_TO_SURFACE,
    surface=region, interactionProperty='self contact', smooth=0.2)

mdb.models['Model-1'].Equation(name='topFlywheelAttachment', terms=((1.0,
    'topAttachment', 5), (-1.0, 'topFlywheelAttachment', 5)))
mdb.models['Model-1'].Equation(name='topFlywheelAttachment2', terms=((1.0,
    'topAttachment', 2), (-1.0, 'topFlywheelAttachment', 2)))
mdb.models['Model-1'].Equation(name='bottomFlywheelAttachment5', terms=((1.0,
    'bottomAttachment', 5), (-1.0, 'bottomFlywheelAttachment', 5)))
mdb.models['Model-1'].Equation(name='bottomFlywheelAttachment2', terms=((1.0,
    'bottomAttachment', 2), (-1.0, 'bottomFlywheelAttachment', 2)))


mdb.models['Model-1'].ActuatorSensorProp(name='flywheel inertia',
                                         realProperties=(flywheelInertia, maxUpset
    ), integerProperties=())
#: The interaction property "IntProp-3" has been created.
a5 = mdb.models['Model-1'].rootAssembly
region=a5.sets['topFlywheelAttachment']
mdb.models['Model-1'].ActuatorSensor(name='top flywheel', createStepName='Initial',
    point=region, interactionProperty='flywheel inertia', noCoordComponents=3,
    unsymm=OFF, noSolutionDepVar=2, userSubUel='U1', dof='5,2', solutionDepVars=(
    flywheelStartVelocity, 0.0))
region=a5.sets['bottomFlywheelAttachment']
mdb.models['Model-1'].ActuatorSensor(name='bottom flywheel', createStepName='Initial',
    point=region, interactionProperty='flywheel inertia', noCoordComponents=3,
    unsymm=OFF, noSolutionDepVar=2, userSubUel='U2', dof='5', solutionDepVars=(
    0.0, 0.0))
mdb.save()


a6 = mdb.models['Model-1'].rootAssembly
region = a6.sets['bottomFlywheelAttachment']
mdb.models['Model-1'].EncastreBC(name='BC-1', createStepName='Initial',
    region=region)

a6 = mdb.models['Model-1'].rootAssembly
s1 = a6.instances['topPipe_weld-1'].edges
side1Edges1 = s1.getSequenceFromMask(mask=('[#4 ]', ), )
a6.Surface(side1Edges=side1Edges1, name='pressureSurface')

region = a6.surfaces['pressureSurface']
mdb.models['Model-1'].Pressure(name='Load-1', createStepName='Weld step',
    region=region, distributionType=UNIFORM, field='', magnitude=appliedPressure,
    amplitude=UNSET)

f1 = a6.instances['bottomPipe_weld-1'].faces
faces1 = f1.getSequenceFromMask(mask=('[#1 ]', ), )
f2 = a6.instances['topPipe_weld-1'].faces
faces2 = f2.getSequenceFromMask(mask=('[#1 ]', ), )
region = regionToolset.Region(faces=faces1+faces2)
mdb.models['Model-1'].Temperature(name='Predefined Field-1',
    createStepName='Initial', region=region, distributionType=UNIFORM,
    crossSectionDistribution=CONSTANT_THROUGH_THICKNESS, magnitudes=(20.0, ))



a6 = mdb.models['Model-1'].rootAssembly
partInstances =(a6.instances['topPipe_weld-1'],
    a6.instances['bottomPipe_weld-1'], )
a6.seedPartInstance(regions=partInstances, size=meshSize, deviationFactor=0.1)

elemType1 = mesh.ElemType(elemCode=CGAX4HT, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=CGAX3T, elemLibrary=STANDARD)
a6 = mdb.models['Model-1'].rootAssembly
f1 = a6.instances['topPipe_weld-1'].faces
faces1 = f1.getSequenceFromMask(mask=('[#1 ]', ), )
f2 = a6.instances['bottomPipe_weld-1'].faces
faces2 = f2.getSequenceFromMask(mask=('[#1 ]', ), )
pickedRegions =((faces1+faces2), )
a6.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2))
a6 = mdb.models['Model-1'].rootAssembly
partInstances =(a6.instances['bottomPipe_weld-1'],
    a6.instances['topPipe_weld-1'], )
a6 = mdb.models['Model-1'].rootAssembly
a1 = mdb.models['Model-1'].rootAssembly
f1 = a1.instances['topPipe_weld-1'].faces
faces1 = f1.getSequenceFromMask(mask=('[#1 ]', ), )
f2 = a1.instances['bottomPipe_weld-1'].faces
faces2 = f2.getSequenceFromMask(mask=('[#1 ]', ), )
pickedRegions = faces1+faces2
a1.setMeshControls(regions=pickedRegions, technique=FREE, allowMapped=False)
a6 = mdb.models['Model-1'].rootAssembly
# Apply partitions
#
# Find the lowest point on the top pipe
e = a.instances['topPipe_weld-1'].edges
upperCutPosition = 0
lowestLocation = 99999
for k in range(len(e)):
    ((x,y,z),) = e[k].pointOn
    if y < lowestLocation:
	lowestLocation = y
upperCutPosition = nearWeldZone+lowestLocation
upperCutPosition = sliceTopInstance('topPipe_weld-1',upperCutPosition)
surfaceNearBottom('topPipe_weld-1','topWeldSurface',upperCutPosition - 0.001)

# Find the highest point on the lower pipe
e = a.instances['bottomPipe_weld-1'].edges
highestLocation = -99999
for k in range(len(e)):
    ((x,y,z),) = e[k].pointOn
    if y > highestLocation:
	highestLocation = y
lowerCutPosition = highestLocation-nearWeldZone
lowerCutPosition = sliceBottomInstance('bottomPipe_weld-1',lowerCutPosition)
surfaceNearTop('bottomPipe_weld-1','bottomWeldSurface',lowerCutPosition + 0.001)

removeEdgeSeeds('topPipe_weld-1')
removeEdgeSeeds('bottomPipe_weld-1')
seedNearZero('topPipe_weld-1',nearWeldMeshSize,(0.25*nearWeldZone) + upperCutPosition)
seedNearZero('bottomPipe_weld-1',nearWeldMeshSize,abs(lowerCutPosition - (0.25*nearWeldZone)))


partInstances =(a6.instances['bottomPipe_weld-1'],
    a6.instances['topPipe_weld-1'], )
a6.generateMesh(regions=partInstances)
mdb.save()
session.viewports['Viewport: 1'].view.fitView()
#: The model database has been saved to "D:\users\yhn\Projects\rezone\M_M.cae".

import job
# mdb.models['Model-1'].keywordBlock.synchVersions(storeNodesAndElements=False)
# mdb.models['Model-1'].keywordBlock.setValues(edited = 0)
# mdb.models['Model-1'].keywordBlock.synchVersions(storeNodesAndElements=False)
# historyBlock = whereIsLastBlock("*End Step") - 1
# remeshCard = "*remesh trigger,distortion = " + "%g" % (distortionLimit, )
# mdb.models['Model-1'].keywordBlock.insert(historyBlock, remeshCard)

mdb.Job(name=primaryJobName, model='Model-1', type=ANALYSIS, explicitPrecision=SINGLE,
        nodalOutputPrecision=SINGLE, description='',
        parallelizationMethodExplicit=DOMAIN, multiprocessingMode=DEFAULT,
        numDomains=1, userSubroutine=userSubFileName, numCpus=1, scratch='',
        echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF)

mdb.models['Model-1'].setValues(noPartsInputFile=OFF)
mdb.jobs[primaryJobName].writeInput(consistencyChecking=OFF)
mdb.save()