#
#   Python script to run parametric study 
#
# Run this from an abaqus command window with abaqus cae -noGUI impact_parameter_study.py
# Make sure that the cae model database file is in the ABAQUS working directory
#
from abaqus import *
from abaqusConstants import *
from caeModules import *
from odbAccess import *
import sys
openMdb('Sphere-film-impact.cae')

# Sphere radius, film shear modulus, and impact velocity
# We could read these from the mdb as well, of course, but simpler to hard code
R = 0.05
mu = 1
V0 = 0.05
minmass = 0.01
maxmass = 0.15
nmasses = 6

paramvals = []
tcvals = []
fmaxvals = []

# Loop over different values of sphere mass
for i in range (0,nmasses):
     m = minmass + (maxmass-minmass)*i/(nmasses-1)
     dimGroup = mu*R**3/(m*V0**2)
     paramvals.append(dimGroup)
     mdb.models['Model-1'].parts['Sphere'].engineeringFeatures.inertias['Inertia-1'].setValues( \
     mass=m)

     thisJob = mdb.jobs['Impact-Simulation']
     thisJob.submit(consistencyChecking=OFF)
     thisJob.waitForCompletion()

     odb = openOdb('Impact-Simulation.odb')

#  Extract the acceleration-time data for the sphere.  You may need to change the name of the node

     step1 = odb.steps['Impact']
     region = step1.historyRegions['Node SPHERE-1.65']
     acceldata = region.historyOutputs['A2'].data

     nn = len(acceldata)
#    Find the contact time by searching the accels from the end backwards to find the first nonzero value
     for n in reversed(xrange(nn)):
           if (acceldata[n][1]>0.000001):
              tcontact = acceldata[n][0]
              break

#    Find the max contact force as the max value of mass*acceleration
     maxforce = m*max([aval[1] for aval in acceldata])

#    Store data for printing later
     tcvals.append(tcontact*V0/R)
     fmaxvals.append(maxforce/(mu*R**2))
   
     odb.close()

resultsFile = open('contactData.csv','w')
for i in range(0, nmasses):
     resultsFile.write('%10.4E,%10.4E,%10.4E\n'% (paramvals[i],tcvals[i],fmaxvals[i]))
resultsFile.close()