from abaqus import *
from abaqusConstants import *

import __main__
import section
import regionToolset
import displayGroupMdbToolset as dgm
import part
import material
import assembly
import step
import interaction
import load
import mesh
import optimization
import job
import sketch
import visualization
import xyPlot
import displayGroupOdbToolset as dgo
import connectorBehavior

Beam_h = [20, 50,30] #Beam Height mm
Beam_w = [10, 10,10] #Beam width mm 

mesh_size = 2.5 #mesh size mm
Beam_d = 400.0 #Beam length mm
Z1 = [50,385,350]
Z2 = [50,5,10]
Pressure = 50.0 # Applied Pressure [MPa]
Young_m = 190000 #Young's Modulus [MPa]
Poisson_r = 0.32 #Poisson's ratio

session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)
for i in range(len(Beam_h)):
    mdb.Model(name='Model-%s'%(Beam_h[i]), modelType=STANDARD_EXPLICIT)
    session.viewports['Viewport: 1'].setValues(displayedObject=None)
    
    session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=OFF, 
        engineeringFeatures=OFF)
    session.viewports['Viewport: 1'].partDisplay.geometryOptions.setValues(
        referenceRepresentation=ON)
        
    s = mdb.models['Model-%s'%(Beam_h[i])].ConstrainedSketch(name='__profile__', 
        sheetSize=200.0)
        
    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)

    #Create a rectangle using two points
    s.setPrimaryObject(option=STANDALONE)
    s.Spot(point=(-Beam_w[i]*0.5, Beam_h[i]*0.5))
    s.Spot(point=(-Beam_w[i]*0.5, -Beam_h[i]*0.5))
    s.Spot(point=(Beam_w[i]*0.5, -Beam_h[i]*0.5))
    s.Spot(point=(Beam_w[i]*0.5, Beam_h[i]*0.5))
    s.Line(point1=(-Beam_w[i]*0.5, Beam_h[i]*0.5), point2=(-Beam_w[i]*0.5, -Beam_h[i]*0.5))
    s.VerticalConstraint(entity=g[2], addUndoState=False)
    s.Line(point1=(-Beam_w[i]*0.5, -Beam_h[i]*0.5), point2=(Beam_w[i]*0.5, -Beam_h[i]*0.5))
    s.HorizontalConstraint(entity=g[3], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[2], entity2=g[3], addUndoState=False)
    s.Line(point1=(Beam_w[i]*0.5, -Beam_h[i]*0.5), point2=(Beam_w[i]*0.5, Beam_h[i]*0.5))
    s.VerticalConstraint(entity=g[4], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[3], entity2=g[4], addUndoState=False)
    s.Line(point1=(Beam_w[i]*0.5, Beam_h[i]*0.5), point2=(-Beam_w[i]*0.5, Beam_h[i]*0.5))
    s.HorizontalConstraint(entity=g[5], addUndoState=False)
    s.PerpendicularConstraint(entity1=g[4], entity2=g[5], addUndoState=False)
        
    p = mdb.models['Model-%s'%(Beam_h[i])].Part(name='Beam', dimensionality=THREE_D, 
        type=DEFORMABLE_BODY)
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']

    #Extrude the sketch "s" over the depth "Beam_d"
    p.BaseSolidExtrude(sketch=s, depth=Beam_d)

    #Define elastic material    
    mdb.models['Model-%s'%(Beam_h[i])].Material(name='Steel')
    mdb.models['Model-%s'%(Beam_h[i])].materials['Steel'].Elastic(table=((Young_m, Poisson_r), ))

    #Define homogeneous solid section 
    mdb.models['Model-%s'%(Beam_h[i])].HomogeneousSolidSection(name='BeamSection', 
        material='Steel', thickness=None)

    #Create a set into part level, used for section assignment    
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    c = p.cells
    cells = c.getSequenceFromMask(mask=('[#1 ]', ), )
    region = p.Set(cells=cells, name='Beam')
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']

    #Assign the section to a defined "region"
    p.SectionAssignment(region=region, sectionName='BeamSection', offset=0.0, 
        offsetType=MIDDLE_SURFACE, offsetField='', 
        thicknessAssignment=FROM_SECTION)
    session.viewports['Viewport: 1'].partDisplay.setValues(sectionAssignments=OFF, 
        engineeringFeatures=OFF, mesh=ON)
    session.viewports['Viewport: 1'].partDisplay.meshOptions.setValues(
        meshTechnique=ON)
        
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    p.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=Z1[i])###first partition in xy plane
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    p.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=Beam_d-Z2[i])###2nd partition over xy plane which is total depth-z1
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    
    c = p.cells
    pickedCells = c.getSequenceFromMask(mask=('[#1 ]', ), )
    d1 = p.datums
    p.PartitionCellByDatumPlane(datumPlane=d1[3], cells=pickedCells)
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    c = p.cells
    pickedCells = c.getSequenceFromMask(mask=('[#1 ]', ), )
    d2 = p.datums
    p.PartitionCellByDatumPlane(datumPlane=d2[4], cells=pickedCells)
    
    #Apply general mesh seed    
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    p.seedPart(size=mesh_size, deviationFactor=0.1, minSizeFactor=0.1)

    #Generate Mesh
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    p.generateMesh()

    #Import part into the instance
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    session.viewports['Viewport: 1'].setValues(displayedObject=a)
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(
        optimizationTasks=OFF, geometricRestrictions=OFF, stopConditions=OFF)
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    a.DatumCsysByDefault(CARTESIAN)
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    a.Instance(name='Beam-1', part=p, dependent=ON)

    # a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    # f1 = a.instances['Beam-1'].faces
    # faces1 = f1.getSequenceFromMask(mask=('[#20 ]', ), )
    # a.Set(faces=faces1, name='BC_Set')

    # Create a surface set
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    s1 = a.instances['Beam-1'].faces
    #side1Faces1 = s1.findAt(((0.0,Beam_h[i]*0.5,(Beam_d-Z2[i]-Z1[i]), ))for pressure in the middle of the surface
    side1Faces1 = s1.findAt(((0.0,Beam_h[i]*0.5,(Z1[i]+0.1)), )) #pressure in the left surface and +0.1 is just a number to flip on right side
    
    a.Surface(side1Faces=side1Faces1, name='PressSurface')
    
 
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(
        adaptiveMeshConstraints=ON)

    #Create analysis    
    mdb.models['Model-%s'%(Beam_h[i])].StaticStep(name='Step-1', previous='Initial', 
        timePeriod=1.0, maxNumInc=100, initialInc=1.0, minInc=1e-05, 
        maxInc=1.0, nlgeom=OFF)
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(step='Step-1')
    session.viewports['Viewport: 1'].assemblyDisplay.setValues(loads=ON, bcs=ON, 
        predefinedFields=ON, connectors=ON, adaptiveMeshConstraints=OFF)

    #Define a region using a set    
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    region = a.surfaces['PressSurface']

    #Applied Load
    mdb.models['Model-%s'%(Beam_h[i])].Pressure(name='Load-1', createStepName='Step-1', 
        region=region, distributionType=UNIFORM, field='', magnitude=Pressure, 
        amplitude=UNSET)

    #Boundary conditions
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly  #access the root assembly
    e1 = a.instances['Beam-1'].edges               #access the edges,e1 contains all of the edges
    edges1 = e1.findAt(((0.0,-Beam_h[i]*0.5,0.0), ))#but we need only one edge so edge1. findAt needs exact coordinates x,y,z.of the edge
    
    region = a.Set(edges=edges1, name='Set-2')
    mdb.models['Model-%s'%(Beam_h[i])].DisplacementBC(name='BC_Set1', createStepName='Step-1', 
        region=region, u1=0.0, u2=0.0, u3=0.0, ur2=0.0, ur3=0.0, 
        amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
        localCsys=None)
        
    a = mdb.models['Model-%s'%(Beam_h[i])].rootAssembly
    e1 = a.instances['Beam-1'].edges
    edges2 = e1.findAt(((0.0,-Beam_h[i]*0.5,Beam_d), ))       

    region = a.Set(edges=edges2, name='Set-3')
    mdb.models['Model-%s'%(Beam_h[i])].DisplacementBC(name='BC_Set2', createStepName='Step-1', 
        region=region, u1=0.0, u2=0.0, ur2=0.0, ur3=0.0, 
        amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', 
        localCsys=None)
    #create node set
    
    p = mdb.models['Model-%s'%(Beam_h[i])].parts['Beam']
    de = p.nodes
    node_set_1 = de.getByBoundingBox(Beam_w[i]/2,Beam_h[i]/2,0.0,Beam_w[i]/2,Beam_h[i]/2,Beam_d/2)
    node_set_2 = de.getByBoundingBox(Beam_w[i]/2,Beam_h[i]/2,Beam_d/2,Beam_w[i]/2,Beam_h[i]/2,Beam_d)
    node_set=(node_set_1,node_set_2)
    p.Set(nodes=node_set, name='path-1')
	
    node_start = de.getByBoundingBox(Beam_w[i]/2,Beam_h[i]/2,0.0)
    p.Set(nodes=node_start, name='path-1-start')   
         
        
    mdb.Job(name='Job-%s'%(i+1), model='Model-%s'%(Beam_h[i]), description='', type=ANALYSIS, 
        atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, 
        memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, 
        explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, 
        modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', 
        scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, 
        numGPUs=0)
        
   # mdb.jobs['Job-%s'%(i+1)].writeInput()