begin sierra newtoncradlefinal

# DEFINE DIRECTION

  define direction Y with vector 0.0 1.0 0.0

# DEFINE MATERIALS:STEEL, E= 200E9 N/m2 = 200E9 Pa , Steel=7840 kg/m3

  begin material outerBalls
    density = 7.48e3
    begin parameters for model elastic
      #youngs modulus = 200.0e9
      youngs modulus = 200e5
      poissons ratio = 0.3
    end parameters for model elastic
  end material outerBalls

  begin material wireAndInnerSphere
    density = 7.48e3
    begin parameters for model elastic
      youngs modulus = 100.0e0
      poissons ratio = 0.3
    end parameters for model elastic
  end material wireAndInnerSphere

  begin rigid body RB1
    reference location = 0 15 0
    include nodes in nodeset_301
  end
  begin rigid body RB2
    reference location = 2 15 0
    include nodes in nodeset_302
  end
  begin rigid body RB3
    reference location = 4 15 0
    include nodes in nodeset_303
  end
  begin rigid body RB4
    reference location = 6 15 0
    include nodes in nodeset_304
  end
  begin rigid body RB5
    reference location = 8 15 0
    include nodes in nodeset_305
  end

  begin truss section wiredup1
    area = 0.1
  end truss section wiredup1

    begin function func_normal_contact_energyAbs_mask
    type = analytic
    expression variable: ts = global timestep
    expression variable: nce_prev = nodal normal_contact_energyAbs_mask_prev
    expression variable: vx = nodal velocity(1)
    expression variable: vy = nodal velocity(2)
    expression variable: vz = nodal velocity(3)
    expression variable: fcn = nodal contact_normal_force_magnitude_mask
    expression variable: ncx = nodal contact_normal_direction(1)
    expression variable: ncy = nodal contact_normal_direction(2)
    expression variable: ncz = nodal contact_normal_direction(3)
    evaluate expression = "nce_prev + (abs((vx*ts*fcn*ncx) + (vy*ts*fcn*ncy) + (vz*ts*fcn*ncz)))"
  end

  begin function func_normal_contact_energyAbs_mask_prev
    type = analytic
    expression variable: time = global time
    expression variable: lastUpdateTime = global lastUpdateTime
    expression variable: nce_old = nodal normal_contact_energyAbs_mask_prev
    expression variable: nce_current = nodal normal_contact_energyAbs_mask
    evaluate expression = "time > lastUpdateTime ? nce_current : nce_old"
  end

  begin function func_mid_to_upper_corset_2_normal_contact_energyAbs_density_mask
    type = analytic
    expression variable: nce = nodal normal_contact_energyAbs_mask
    expression variable: area = nodal surface_area
    evaluate expression = "nce/area"
  end

  begin function func_contact_energy_mask
    type = analytic
    expression variable: mask = nodal contact_status
    expression variable: ce = nodal contact_energy
    evaluate expression = " (mask == 1) ? ce : 0"
  end

  begin function func_contact_normal_force_magnitude_mask
    type = analytic
    expression variable: mask = nodal contact_status
    expression variable: cnfm = nodal contact_normal_force_magnitude
    evaluate expression = " (mask == 1) ? cnfm : 0"
  end

  begin function func_contact_area_mask
    type = analytic
    expression variable: mask = nodal contact_status
    expression variable: area = nodal contact_area
    evaluate expression = " (mask == 1) ? area : 0"
  end

# DEFINE FEM MODEL
#Note: Comment out gensesis files for desired test geometery configuration(1, 2, or 3 balls dropped)

  begin finite element model mesh
    Database Name =   newtoncradle_final_1balldrop.g
    Database Type = exodusII

# DEFINE BLOCKS(block set 1 series are inner volumes, series 100 are outer volumes, and 300 are set
# of curves/truss)

    begin parameters for block block_101 block_102 block_103 block_104 block_105
      material = outerBalls
      model = elastic
    end parameters for block block_101 block_102 block_103 block_104 block_105

    begin parameters for block block_1 block_2 block_3 block_4 block_5
      material = wireAndInnerSphere
      model= elastic
    end parameters for block block_1 block_2 block_3 block_4 block_5

    begin parameters for block block_300
      material = wireAndInnerSphere
      model = elastic
      section = wiredup1
    end parameters for block block_300

  end finite element model mesh

# DEFINE PROBLEM TIME AND TIME STEP PARAMETERS

  begin presto procedure myProcedure
    begin time control
      begin time stepping block p0
        start time = 0.0
        begin parameters for presto region myRegion
        end parameters for presto region myRegion
      end time stepping block p0
      #termination time = 9
      termination time = 4.5
    end time control

# DEFINE BOUNDARY CONDITIONS
# ONLY DEFINE RIGID BODY AS FIXED DISPLACEMENT DUE TO REFERENCE LOCATION ON RIGID BODY (do not want
# entire nodeset rigid)

    begin presto region myRegion

      extensive rigid body vars output = off      

      begin node based time step parameters set1
        step interval = 100
      end node based time step parameters set1

      use finite element model mesh

      begin fixed displacement
        rigid body = RB1
        components = x y z
      end
      begin fixed displacement
        rigid body = RB2
        components = x y z
      end
      begin fixed displacement
        rigid body = RB3
        components= x y z
      end
      begin fixed displacement
        rigid body = RB4
        components= x y z
      end
      begin fixed displacement
        rigid body = RB5
        components = x y z
      end

      begin fixed rotation
        rigid body = RB1 
        components= x y
      end 
      begin fixed rotation
        rigid body= RB2
        components= x y
      end
      begin fixed rotation
        rigid body = RB3
        components= x y
      end 
      begin fixed rotation
        rigid body = RB4
        components= x y
      end 
      begin fixed rotation
        rigid body = RB5
        components= x y
      end 

      begin gravity
        direction = y
        gravitational constant = -9.81
      end gravity

# DEFINE PROBLEM REGIONS  

      begin contact definition contacts
        skin all blocks = on exclude block_300
        search = dash
        begin interaction defaults
          friction model = sticky
          general contact = on
          self contact  = off
          constraint formulation = node_face
        end
        BEGIN CONSTANT FRICTION MODEL sticky
          FRICTION COEFFICIENT = 0.3 
        END 
      end contact definition contacts


      begin user variable normal_contact_energyAbs_mask
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end
      begin user variable normal_contact_energyAbs_mask_prev
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end
      begin user variable normal_contact_energyAbs_density_mask
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end
      begin user variable contact_normal_force_magnitude_mask
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end
      begin user variable contact_area_mask
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end
      begin user variable lastUpdateTime
        type = global real length = 1
        global operator = max
        initial value = 0.0
        use with restart
      end
      begin user variable contact_energy_mask
        type = node real length = 1
        global operator = sum
        initial value = 0.0
        use with restart
      end

      begin user output
        include all blocks
        compute nodal normal_contact_energyAbs_mask_prev as function func_normal_contact_energyAbs_mask_prev
        compute global lastUpdateTime from expression "time"
        compute at every step
      end

      begin user output
        include all blocks
        compute nodal normal_contact_energyAbs_mask as function func_normal_contact_energyAbs_mask
        compute nodal contact_normal_force_magnitude_mask as function func_contact_normal_force_magnitude_mask
        compute nodal contact_area_mask as function func_contact_area_mask
        compute nodal contact_energy_mask as function func_contact_energy_mask
        compute at every step
      end

      begin user output
        include all blocks
        compute global nce as sum of nodal normal_contact_energyAbs_mask
        compute global cf as sum of nodal contact_normal_force_magnitude_mask
        compute global my_area as sum of nodal contact_area_mask
        compute global ce as sum of nodal contact_energy_mask
        compute at every step
      end

      begin history output
        Database Name = newton_cradle.h
        Database Type = exodusII
        At time 0, Increment = .1
        variable = global timestep
        variable = global lastUpdateTime
        variable = global nce
        variable = global cf
        variable = global my_area
        variable = global ce
      end

    Begin Heartbeat Output normalized
        Stream Name = NewtonCradle.csv
        Format = SpyHis
        Start Time = 0.0
        At Time 0 Increment = 0.005
        Termination Time = 9.0
        Global kinetic_energy as KineticEnergy
        Global External_energy as ExternalEnergy
        Global Internal_energy as InternalEnergy
        Global SE as StrainEnergy
        global normIESE as Normalized_Dissipation
        global PEabs_negEEplusMaxEE as PEabs_negEEplusMaxEE
        global AbsTotalE as AbsTotalE
        global PE as PE
        global contact_energy as contact_energy
        global TotalE as TotalE
        global hourglass_energy as hourglass_energy

        global Kinetic_energyKJ as Kinetic_energyKJ
        global External_energyKJ as External_energyKJ
        global Internal_energyKJ as Internal_energyKJ
        global Potential_energyKJ as Potential_energyKJ
        global TotalE_KJ as TotalE_KJ

        global Sum_SE_EE_KE as Sum_SE_EE_KE
        global Sum_IE_EE_KE as Sum_IE_EE_KE
        global Et_Minus_E as Et_Minus_E
        global normEnergy as normEnergy
        global energyNorm as energyNorm
        global DissEnergyOVERKEo as DissEnergyOVERKEo
        global Dissipation_energy as Dissipation_Energy
      End Heartbeat Output normalized

     Begin User Output
       compute global SE as sum of element strain_energy
       compute global Sum_SE_EE_KE from expression "kinetic_energy+external_energy+SE"
       compute global Sum_IE_EE_KE from expression "kinetic_energy+external_energy+internal_energy"
       compute global IE_minus_SE from expression "internal_energy-SE"
       compute global normIESE from expression "IE_minus_SE/606990" #606990  is max EE

       # shift and flip EE to show PE
       compute global PEabs_negEEplusMaxEE from expression "abs(-1*external_energy+606990)"

       compute global AbsTotalE from expression "kinetic_energy+internal_energy+PEabs_negEEplusMaxEE"
       compute global PE from expression "(-1*external_energy)+606990"
       compute global TotalE from expression "kinetic_energy+internal_energy+PE"

       compute global Kinetic_energyKJ from expression "kinetic_energy/1000"
       compute global External_energyKJ from expression "external_energy/1000"
       compute global Internal_energyKJ from expression "internal_energy/1000"
       compute global Potential_energyKJ from expression "PE/1000"
       compute global TotalE_KJ from expression "TotalE/1000"

       compute global Et_Minus_E from expression "Sum_IE_EE_KE-Sum_SE_EE_KE"
       compute global normEnergy from expression "Et_Minus_E/1.22634E6"

       compute global energyNorm from expression "max(abs(Kinetic_Energy), abs(Internal_Energy))"
       compute global max_KineticEnergy from expression "max(abs(Kinetic_Energy), abs(0))"
       compute global max_InternalEnergy from expression "max(abs(Internal_Energy), abs(0))"
       compute global max_ExternalEnergy from expression "max(abs(External_Energy), abs(0))"

       compute global max_Sum_KE_IE from expression "max_KineticEnergy+max_InternalEnergy"
       compute global max_Sum_EE_KE_PE from expression \#
                                   "max_KineticEnergy+max_InternalEnergy+max_ExternalEnergy"
       compute global Dissipation_Energy from expression "max_Sum_KE_IE-max_ExternalEnergy"
       compute global DissEnergyOVERKEo from expression "max_Sum_KE_IE/602660" #pulled 602660 as max KE
     End User Output

      begin Results Output output_presto
        Database Name =   newtoncradle_final_1balldrop.e 
        Database Type = exodusII
        At time 0, Increment = .1

        nodal variables   = force_external        as f_ext
        nodal variables   = velocity              as vel
        nodal variables   = displacement          as displ
        nodal variables   = reaction              as reactions
        nodal variables   = force_internal        as f_int
        nodal variables   = force_contact         as f_cont
        nodal variables   = contact_status        as contact_stat

        element variables = stress                as stress
        element variables = von_mises             as vonmises
        element variables = effective_log_strain  as EffLogStrain
        element variables = truss_force           as TrussForce
        element variables = strain_energy         as StrainEnergy

        global variables = timestep               as timestep
        global variables = external_energy        as ExternalEnergy
        global variables = internal_energy        as InternalEnergy
        global variables = kinetic_energy         as KineticEnergy
        global variables = hourglass_energy       as HourglassEnergy #should be 0
        global variables = contact_energy         as Contactenergy
        global variables = momentum               as Momentum

        global variables = energyNorm as energyNorm
        global variables = max_KineticEnergy as max_KineticEnergy
        global variables = max_InternalEnergy as max_InternalEnergy
        global variables = max_ExternalEnergy as max_ExternalEnergy
        global variables = max_Sum_EE_KE_PE as max_Sum_EE_KE_PE

        global variables = max_Sum_KE_IE as max_Sum_KE_IE
        global variables = Dissipation_energy as Dissipation_Energy
        global variables = DissEnergyOVERKEo as DissEnergyOVERKEo

        global variables = normEnergy   as normEnergy
        global variables = Et_Minus_E   as Et_Minus_E
        global variables = Sum_SE_EE_KE as Sum_SE_EE_KE
        global variables = Sum_IE_EE_KE as Sum_IE_EE_KE
        global variables = SE as strain_energy_summedElements
        global variables = PEabs_negEEplusMaxEE as PEabs_negEEplusMaxEE
        global variables = AbsTotalE as AbsTotalE
        global variables = PE as PE
        global variables = TotalE as TotalE
        global variables = normIESE as normIESE
        global variables = IE_minus_SE as IE_minus_SE

      end results output output_presto

    end presto region myRegion

  end presto procedure myProcedure

end sierra newtoncradlefinal