#
#    Coupled Eulerian-Lagrangian Analysis with Abaqus/Explicit
#    Elastic Dam Break model
#
from abaqus import *
from abaqusConstants import *
session.viewports['Viewport: 1'].makeCurrent()
session.viewports['Viewport: 1'].maximize()
from caeModules import *
from driverUtils import executeOnCaeStartup
import uti, osutils

def getDriver():
    return os.environ.get('ABA_COMMAND', 'abaqus')

executeOnCaeStartup()
Mdb()

mdb.Model(name='DamBreech')
del mdb.models['Model-1']

# Import dam and reservoir parts
abaArgs = ["fetch", "-j"]
uti.spawnAndWait(getDriver(), abaArgs + ['dambreach_dam.sat',])
uti.spawnAndWait(getDriver(), abaArgs + ['dambreach_reservoir.sat',])

acis = mdb.openAcis('dambreach_dam.sat')
mdb.models['DamBreech'].PartFromGeometryFile(name='Dam', geometryFile=acis,  dimensionality=THREE_D, type=DEFORMABLE_BODY)
osutils.remove('dambreach_dam.sat')

acis = mdb.openAcis('dambreach_reservoir.sat')
mdb.models['DamBreech'].PartFromGeometryFile(name='Reservoir', geometryFile=acis,  dimensionality=THREE_D, type=EULERIAN)
osutils.remove('dambreach_reservoir.sat')

#    Define materials
#
mdb.models['DamBreech'].Material(name='Water')
mdb.models['DamBreech'].materials['Water'].Density(
    table=((1000.0, ), ))
mdb.models['DamBreech'].materials['Water'].Eos(
    table=((1500.0, 0.0, 0.0), ),
    type=USUP)
mdb.models['DamBreech'].materials['Water'].Viscosity(table=((0.001, ), ))

mdb.models['DamBreech'].Material(name='Elastic')
mdb.models['DamBreech'].materials['Elastic'].Density(
    table=((1100.0, ), ))
mdb.models['DamBreech'].materials['Elastic'].Elastic(
    table=((12.e6, 0.4), ))

#    Assign sections
#
mdb.models['DamBreech'].HomogeneousSolidSection(
    material='Elastic',
    name='Solid')

mdb.models['DamBreech'].EulerianSection(
    data={'water-1': 'Water'},
    name='Eulerian')

p = mdb.models['DamBreech'].parts['Reservoir']
c = p.cells
cells = c
region = (cells, )

p.SectionAssignment(
    region=region,
    sectionName='Eulerian')

p = mdb.models['DamBreech'].parts['Dam']
c = p.cells
cells = c
region = (cells, )

p.SectionAssignment(
    region=region,
    sectionName='Solid')

#    Assemble parts
#
a = mdb.models['DamBreech'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)

p = mdb.models['DamBreech'].parts['Dam']
a.Instance(dependent=OFF, name='Dam-1', part=p)

p = mdb.models['DamBreech'].parts['Reservoir']
a.Instance(dependent=OFF, name='Reservoir-1', part=p)

c1 = a.instances['Reservoir-1'].cells
pickedCells = c1
v1 = a.instances['Dam-1'].vertices
e1 = a.instances['Dam-1'].edges
a.PartitionCellByPlanePointNormal(
    point=v1.findAt(coordinates=(0.005, 0.15, 0.005)),
    normal=e1.findAt(coordinates=(0.00125, 0.15, 0.005)), 
    cells=pickedCells)

f1 = a.instances['Reservoir-1'].faces
a.DatumPlaneByOffset(
    plane=f1.findAt(coordinates=(0.038333, 0.15, 0.003333)), 
    flip=SIDE2, offset=0.01)

cells = c1.findAt(((-0.015, 0.0, 0.003333), ), ((0.105, 0.1, 0.003333), ))
pickedCells = cells
d1 = a.datums
a.PartitionCellByDatumPlane(datumPlane=d1[7], cells=pickedCells)

c1 = a.instances['Dam-1'].cells
cells1 = c1
c2 = a.instances['Reservoir-1'].cells
cells2 = c2
pickedCells = cells1+cells2
f1 = a.instances['Reservoir-1'].faces
a.PartitionCellByExtendFace(
    extendFace=f1.findAt(coordinates=(-0.036667, 0.079, 0.003333)),
    cells=pickedCells)

a.rotate(
    angle=90.0,
    axisDirection=(10.0, 0.0, 0.0),
    axisPoint=(0.0, 0.0, 0.0),
    instanceList=('Dam-1', 'Reservoir-1'))

f1 = a.instances['Reservoir-1'].faces

faces1 = f1.findAt(
    ((-0.015, -0.003333, 0.0), ),
    ((0.038333, -0.001667, 0.0), ))
a.Set(faces=faces1, name='Reservoir-Bottom')

faces1 = f1.findAt(
    ((0.105, -0.003333, 0.052667), ),
    ((0.105, -0.001667, 0.099333), ),
    ((0.105, -0.001667, 0.143333), ))
a.Set(faces=faces1, name='Reservoir-Right')

faces1 = f1.findAt(
    ((0.003333, 0.0, 0.052667), ),
    ((0.001667, -0.005, 0.099333), ),
    ((0.038333, 0.0, 0.052667), ),
    ((0.038333, -0.005, 0.099333), ),
    ((0.071667, 0.0, 0.143333), ),
    ((0.001667, -0.005, 0.143333), ),
    ((0.071667, -0.005, 0.052667), ),
    ((0.003333, 0.0, 0.099333), ),
    ((0.003333, 0.0, 0.143333), ),
    ((0.038333, -0.005, 0.143333), ),
    ((0.003333, -0.005, 0.052667), ),
    ((0.071667, 0.0, 0.099333), ))
a.Set(faces=faces1, name='Reservoir-v2')

c1 = a.instances['Reservoir-1'].cells

cells1 = c1.findAt(
    ((0.105, -0.003333, 0.052667), ),
    ((0.105, -0.001667, 0.099333), ))
a.Set(cells=cells1, name='Water-initial')

cells1 = c1
a.Set(cells=cells1, name='AllReservoir')


f1 = a.instances['Dam-1'].faces

faces1 = f1.findAt(
    ((0.0, -0.003333, 0.102667), ),
    ((0.001667, -0.003333, 0.15), ))
a.Set(faces=faces1, name='Dam-Fix')

faces1 = f1.findAt(
    ((0.001667, 0.0, 0.052667), ),
    ((0.001667, -0.005, 0.102667), ),
    ((0.003333, -0.005, 0.052667), ),
    ((0.003333, 0.0, 0.102667), ))
a.Set(faces=faces1, name='Dam-v2')

v1 = a.instances['Dam-1'].vertices

verts1 = v1.findAt(((0.0, -0.005, 0.0), ))
a.Set(vertices=verts1, name='Dam-Tip')

c1 = a.instances['Dam-1'].cells

cells1 = c1
a.Set(cells=cells1, name='AllDam')

#    Define step and output
#
mdb.models['DamBreech'].ExplicitDynamicsStep(
    description='gravity flow and deformation of plate',
    name='Flow',
    previous='Initial', 
    timePeriod=0.4)

mdb.models['DamBreech'].fieldOutputRequests['F-Output-1'].setValues(
    region=a.sets['AllDam'],
    timeInterval=0.01,
    timeMarks=OFF)

mdb.models['DamBreech'].FieldOutputRequest(
    createStepName='Flow',
    name='F-Output-2',
    region=a.sets['AllReservoir'],
    timeInterval=0.01,
    variables=('EVF', 'DENSITYVAVG'))

mdb.models['DamBreech'].HistoryOutputRequest(
    createStepName='Flow',
    name='H-Output-2',
    region=a.sets['Dam-Tip'],
    sectionPoints=DEFAULT, 
    timeInterval=0.01, variables=('U1', 'U3'))

mdb.models['DamBreech'].steps['Flow'].Restart(
    numberIntervals=1,
    overlay=ON, 
    timeMarks=OFF)   

#    Define general contact
#
mdb.models['DamBreech'].ContactProperty('nofric')

mdb.models['DamBreech'].ContactExp(
    createStepName='Flow',
    name='cel-contact')

mdb.models['DamBreech'].interactions['cel-contact'].includedPairs.setValuesInStep(
    stepName='Flow',useAllstar=ON)

mdb.models['DamBreech'].interactions['cel-contact'].contactPropertyAssignments.appendInStep(
    assignments=((GLOBAL, SELF, 'nofric'), ), stepName='Flow')


#    Apply gravity load
#
mdb.models['DamBreech'].Gravity(
    comp3=-9.81,
    createStepName='Flow', 
    distributionType=UNIFORM,
    name='Gravity')

#    Apply boundary conditions
#
mdb.models['DamBreech'].VelocityBC(
    createStepName='Initial', 
    distributionType=UNIFORM,
    name='FixDam', 
    region=a.sets['Dam-Fix'],
    v1=0.0, 
    v3=0.0)

mdb.models['DamBreech'].VelocityBC(
    createStepName='Initial', 
    distributionType=UNIFORM,
    name='Dam-v2', 
    region=a.sets['Dam-v2'],
    v2=0.0)

mdb.models['DamBreech'].VelocityBC(
    createStepName='Initial', 
    distributionType=UNIFORM,
    name='Reservoir-v2',
    region=a.sets['Reservoir-v2'],
    v2=0.0)

mdb.models['DamBreech'].VelocityBC(
    createStepName='Initial', 
    distributionType=UNIFORM,
    name='Reservoir-bottom',
    region=a.sets['Reservoir-Bottom'], 
    v3=0.0)

mdb.models['DamBreech'].VelocityBC(
    createStepName='Initial', 
    distributionType=UNIFORM,
    name='Reservoir-right',
    region=a.sets['Reservoir-Right'],
    v1=0.0)

#    Assign Eulerian material
#
instList = (a.instances['Reservoir-1'], )
rgn1 = a.sets['Water-initial']
fract1 = (1, )
mdb.models['DamBreech'].MaterialAssignment(
    name='Water-initial', 
    instanceList=instList,
    useFields=False,
    assignmentList=((rgn1, fract1), ))


#    Add hydrostatic pressure definition via keywords editor
#
mdb.models['DamBreech'].keywordBlock.synchVersions(storeNodesAndElements=False)

mdb.models['DamBreech'].keywordBlock.insert(36, """*Initial Conditions, type=STRESS, GEOSTATIC
Water-initial, 0.0, 0.14, -1373.4, 0.0, 1.0, 1.0""")


#    Mesh instances
#
partInstances =(a.instances['Reservoir-1'], )
a.seedPartInstance(regions=partInstances, size=0.005, deviationFactor=0.1)
a.generateMesh(regions=partInstances)

elemType1 = mesh.ElemType(elemCode=EC3D8R, elemLibrary=EXPLICIT)
elemType2 = mesh.ElemType(elemCode=UNKNOWN_WEDGE, elemLibrary=EXPLICIT)
elemType3 = mesh.ElemType(elemCode=UNKNOWN_TET, elemLibrary=EXPLICIT)

c1 = a.instances['Reservoir-1'].cells
cells1 = c1
pickedRegions =(cells1, )
a.setElementType(
    regions=pickedRegions,
    elemTypes=(elemType1, elemType2, elemType3))


partInstances =(a.instances['Dam-1'], )
a.seedPartInstance(regions=partInstances, size=0.002, deviationFactor=0.1)
e1 = a.instances['Dam-1'].edges
pickedEdges = e1.findAt(
    ((0.00125, -0.005, 0.079), ),
    ((0.00375, 0.0, 0.079), ),
    ((0.00375, 0.0, 0.0), ),
    ((0.00125, -0.005, 0.15), ),
    ((0.00125, 0.0, 0.15), ),
    ((0.00375, -0.005, 0.0), ))
a.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FIXED)
a.generateMesh(regions=partInstances)


elemType1 = mesh.ElemType(elemCode=C3D8R, elemLibrary=EXPLICIT)
elemType2 = mesh.ElemType(elemCode=C3D6, elemLibrary=EXPLICIT)
elemType3 = mesh.ElemType(elemCode=C3D4, elemLibrary=EXPLICIT)

c1 = a.instances['Dam-1'].cells
cells1 = c1
pickedRegions =(cells1, )
a.setElementType(
    regions=pickedRegions,
    elemTypes=(elemType1, elemType2, elemType3))

# Create job and write input file
#
mdb.Job(name='dam_deflection_cel', model='DamBreech', type=ANALYSIS, 
    description='CEL Analysis-Elastic dam deflection')

mdb.jobs['dam_deflection_cel'].writeInput(consistencyChecking=OFF)

# Save model
#
mdb.saveAs('dam_deflection_example')