Hi,
I'm new to the list and have done very little with symbolic math so
I'm hoping there is a simple way to do what I need to do. I have two
scalar functions f and g defined on the 6 dimensional space of
symmetric tensors. I need to calculate Df, Dg, and D^2 g and then
generate optimized C++ code to evaluate everything. Df and Dg can be
stored as 6x1 matrices and D^2g as 6x6. I wrote a simple script to
do that by building up the expressions for f and g and then calling
diff to build lists of the derivatives and then ccode to write
everything to an inline function. The problem is that the resulting
expressions are extremely large and my naive attempt at using the
simplify function crashed after using up all the memory on my machine
(8G).
The functions f and g are really just functions of the 3 tensor
invariants so if I were to do this by hand I'd build the gradient and
Hessian of each invariant and the gradient and Hessian of f and g and
then build Df, Dg, and D^2g by the chain rule. That would catch a lot
of repeated subexpressions and a few more simplifications could be
identified due to the fact that the tensor invariants contain rational
and trigonometric functions that tend to be repeated in the
derivatives. I doubt the hand generated code would be optimal though.
Does this seem like something the cse module or some of the other
specialized simplification support in sympy can handle? I'm pasting
the short python script in below.
Thanks,
Chris
mohrCoulomb.py
--------------------------
from sympy import *
from sympy import pprint
#
# Evaluate the yield surface, plastic flow rule, and derivatives, and
generate an inline C++ function
#
cfile = open("tmp.h",'w')
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) #t is
positive so perturb the denominator to prevent NaNs for zero stress
intial conditions
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")
#
#form the plastic potential and the plastic flow rule and the gradient
in stress space
#
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()
test.cpp
-------------
#include <cfloat>
#include <cmath>
#include <iostream>
#include <fstream>
#define nSymTen 6
#define sXX 0
#define sYY 1
#define sZZ 2
#define sYZ 3
#define sZY 3
#define sXZ 4
#define sZX 4
#define sXY 5
#define sYX 5
#include "tmp.h"
int main()
{
std::ofstream ffile;
ffile.open ("f.txt");
double
c(10.0),phi(2.0*M_PI*20.0/360.0),psi(0.0),stress[nSymTen],f,df[nSymTen],r[nSymTen],dr[nSymTen]
[nSymTen];
double sInit=1.0;
for (int i=0;i<101;i++)
for (int j=0;j<101;j++)
for (int k=0;k<101;k++)
{
stress[sXX]=sInit - k;
stress[sYY]=sInit - j;
stress[sZZ]=sInit - i;
stress[sYZ]=0.0;
stress[sXZ]=0.0;
stress[sXY]=0.0;
mohrCoulomb(phi,psi,c,stress,f,df,r,dr);
ffile<<stress[sXX]<<","<<stress[sYY]<<","<<stress[sZZ]<<","<<f<<std::endl;
}
ffile.close();
return 0;
}
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
