gistfile1.txt

from StringIO import StringIO

from sympy import *
from sympy import pprint

cfile = StringIO()
cfile.write("""
inline void mohrCoulomb(const double& phi, const double& psi, const
double& c, const double stress[nSymTen], //inputs
                              double& f, double df[nSymTen], double
r[nSymTen], double dr[nSymTen][nSymTen]) //outputs
{
 double theta,dtheta[nSymTen],ddtheta[nSymTen][nSymTen],
   stress_m,dstress_m[nSymTen],ddstress_m[nSymTen][nSymTen],
   stress_bar,dstress_bar[nSymTen],ddstress_bar[nSymTen][nSymTen],
   df_theta,df_stress_m,df_stress_bar,
   dg_theta,dg_stress_m,dg_stress_bar,
   ddg_theta,ddg_stress_m,ddg_stress_bar;
 const double& stress_xx(stress[sXX]),
   stress_yy(stress[sYY]),
   stress_zz(stress[sZZ]),
   stress_yz(stress[sYZ]),
   stress_zx(stress[sZX]),
   stress_xy(stress[sXY]);
""")
#set up symbols for stress tensor
stress = (stress_xx,stress_yy,stress_zz,stress_yz,stress_zx,stress_xy) = symbols(("stress_xx","stress_yy","stress_zz","stress_yz","stress_zx","stress_xy"))
stress_keys=['sXX','sYY','sZZ','sYZ','sZX','sXY']
#set up symbols for numerical constants
(phi,psi,c) = symbols(("phi","psi","c"))
#
#form expressions for invariants of the stress tensor
#
s = (stress_xx+stress_yy+stress_zz)/sqrt(3)
t = sqrt((stress_xx - stress_yy)**2 + (stress_yy - stress_zz)**2 +
(stress_zz - stress_xx)**2 + Rational(6)*stress_xy**2 +
Rational(6)*stress_yz**2 + Rational(6)*stress_zx**2)/sqrt(3)
s_x = (Rational(2)*stress_xx - stress_yy - stress_zz)/Rational(3)
s_y = (Rational(2)*stress_yy - stress_xx - stress_zz)/Rational(3)
s_z = (Rational(2)*stress_zz - stress_xx - stress_yy)/Rational(3)
J3 = s_x*s_y*s_z - s_x*stress_yz**2 - s_y*stress_zx**2 - s_z*stress_xy**2 + Rational(2)*stress_xy*stress_yz*stress_zx
theta = asin(-Rational(3)*sqrt(6)*J3/(t**3+1.0e-8))/Rational(3)
stress_m = s/sqrt(3)
stress_bar = sqrt(Rational(3,2))*t
#
#form the yield surface and gradient in stress space
#
f = stress_m*sin(phi)+stress_bar*(cos(theta)/sqrt(3) -
sin(theta)*sin(phi)/Rational(3)) - c*cos(phi)
df = diff(f,stress)
#
#write the C code for f and gradient in stress space
#
cfile.write("  f = "+ccode(f)+";\n\n")
for key,expr in zip(stress_keys,df):
   cfile.write("  df["+key+"] = "+ccode(expr)+";\n\n")

g = stress_m*sin(phi)+stress_bar*(cos(theta)/sqrt(3) -
sin(theta)*sin(phi)/Rational(3)) - c*cos(phi)
#dg = diff(g,stress)
#ddg = diff(dg,stress)
#for key,expr in zip(stress_keys,dg):
#   cfile.write("  r["+key+"] = "+ccode(expr)+";\n\n")
#for keyr,expr_list in zip(stress_keys,ddg):
#   for keyc,expr in zip(stress_keys,expr_list):
#       cfile.write("  dr["+keyr+"]["+keyc+"] = "+ccode(expr)+";\n\n")
cfile.write("}\n")
#cfile.close()

cfile.seek(0)
print cfile.read()