Good point Frank -- Thanks.
General instructions:
=============================================================================
The attached Python script sprind.py is the equivalent of the Fortran
file sprind.f. It should be imported as a module using
import sprind
or
from sprind import sprind
The module implements a function called sprind(). The function prototype is
sprind(s, lstr, ndi, nshr)
where:
s is the stress or strain tensor
lstr is a variable to distinguish between stress or strain. Valid values
are 1 for stress or 2 for strains
ndi is the number of direct tensor components
nshr is the number of shear components.
The function returns a tuple (pValues, pDirections) where pValues is a
sequence of principle values and pDirections is a sequence of sequences
of principal directions.
For example, the principal direction corresponding to pValues[i] is
pDirections[i] with i ={0,1,2}.
sprind.py just implements the sprind() function. The function has to be
called with valid arguments in order to compute the principal values and
principal directions.
The attached script fprin.py is the Python equivalent of the fprin.f
postprocessing program which accesses a specified output database and
reads the stress (or strain values) and calls the above sprind()
function. The principal values and principal directions obtained are
then printed to a text file (by default to pValues.dat).
Usage:
abaqus python fprin.py odbName.odb
If odbName is not specified on the command line then the user will be
prompted to enter the name of the output database. If the file
pValues.dat already exists then the user will be prompted to enter the
name of the text file to which the principal values and directions
should be written. Note that any existing file with the name specified
will be overwritten.
Once the program successfully runs to completion, all the principal
values and principal directions will be available from the text file.
This script can be easily changed such that the principal values and
principal directions are written to the output database instead of a
text file. For more information on how to write data to the .odb file,
please refer to Abaqus Scripting User's Manual.
abaqus python fprin.py
Enter the name of the output database: viewer_tutorial
Enter the name of the text file: pValues.dat
Successfully computed the principal values and directions.
These principal values and directions are now available in pValues.dat
Limitations:
The script assumes that the stress or strain values are written to the
.odb at the integration points. Since tensor data is written to the .odb
in single precision the results obtained using sprind.py may be less
accurate than the ones that are obtained using sprind.f which operates
on the double precision results file.
For more information see:
'Obtaining stress invariants principal stress/strain values and
directions and rotating tensors,'
Section 2.1.9 of the Version 6.7 Abaqus User Subroutines Reference Manual
Section 2.1.10 of the Version 6.8 Abaqus User Subroutines Reference Manual
=============================================================================
Listing frpin.py:
=============================================================================
# python script to test sprind.py
from sprind import sprind
from odbAccess import *
import sys
import math
import os
# Get odb Name
if len(sys.argv)>1:
odbName = sys.argv[1]
else:
odbName = raw_input('Enter the name of the output database:').strip()
if odbName.find('.odb')==-1:
odbName += '.odb'
odb = openOdb(odbName)
# Create the text file
fileName = 'pValues.dat'
if os.path.isfile(fileName):
fileName = raw_input('Enter the name of the text file: ').strip()
outFile = open(fileName,'w')
# Process all the frames of all the steps
steps = odb.steps.values()
ndiDict =
{TENSOR_3D_FULL:3,TENSOR_3D_PLANAR:3,TENSOR_3D_SURFACE:2,TENSOR_2D_PLANAR:3,TENSOR_2D_SURFACE:2}
varDict = {'S':'STRESS','E':'STRAINS'}
#Loop through each step
for step in steps:
stepName = step.name
#Loop through each frame
for frame in step.frames:
frameId = frame.frameId
outFile.write('**Step Name:\t'+stepName+'\tFrame
Number:\t'+str(frameId)+'\n')
fo = frame.fieldOutputs
#Loop through stress and strain
for var in varDict.keys():
if fo.has_key(var):
varFO = fo[var]
vals = varFO.values
val0 = vals[0]
ndi = ndiDict[varFO.type]
nshr = len(val0.data)-ndi
for val in vals:
#Loop through each element and integration point
line = '%s\tElement: %d\tIntegration Point: %d
\n'%(varDict[var],val.elementLabel,val.integrationPoint)
outFile.write(line)
s = val.data
ps,an = sprind(s,1,ndi,nshr)
for i in range(3):
line1 = 'PS%d =%+1.8E'%(i+1,ps[i])
line2 = '\nPD%d1 =%+1.8E PD%d2 =%+1.8E
PD%d3 =%+1.8E\n\n'%(i+1,an[i][0],i+1,an[i][1],i+1,an[i][2])
outFile.write(line1+line2)
outFile.close()
print 'Successfully computed the principal values and directions.'
print 'These principal values and directions are now available in %s' %
(fileName)
=============================================================================
Listing of sprind.py:
=============================================================================
#Python function to compute the principal directions and values
def sprind(s,lstr,ndi,nshr):
import math
pi23=2.094395102393195
cons1=10000.0
precis=2.22e-16
preciz=precis*cons1
zero=0.0
half=0.5
one=1.e0
two=2.e0
third=one/3.e0
asmall=1.e0/1.e36
ps=[zero,zero,zero]
an=[[zero,zero,zero],[zero,zero,zero],[zero,zero,zero]]
# LSTR=1 - STRESS
# LSTR=2 - STRAIN
#
# THE EIGENVALUES OF S(*) ARE PUT IN PS(K1),K1=1,3
#
# DIRECTION COSINES OF PS(K1) ARE PUT IN AN(K1,K2),K2=1,2,3
#
ndip1=ndi
ndip2=ndi+1
ndip3=ndi+2
#
# NO SHEAR COMPONENTS
#
if nshr==0:
an[0][0]=one
an[1][1]=one
an[2][2]=one
for i in range(3):
ps[i]=s[i]
#
# ONE SHEAR COMPONENT: FIRST NORMALIZE WITHOUT DEVIATORIC PART
#
elif nshr==1:
if ndi==0:
sh=zero
s11=zero
s22=zero
elif ndi==1:
sh=half*s[0]
s11=sh
s22=-sh
else:
sh=half*(s[0]+s[1])
s11=half*(s[0]-s[1])
s22=-s11
if lstr==1:
s12=s[ndip1]
else:
s12=half*s[ndip1]
facd= math.fabs(s11)
facs= math.fabs(s12)
fact= max(facd,facs)
#
# ESSENTIALLY UNIT MATRIX
#
if fact<=math.fabs(preciz*sh) or facs<preciz*facd:
ps[0]=s[0]
ps[1]=s[1]
an[0][0]=one
an[1][1]=one
elif fact<asmall:
# -- We've been given very small [zero] numbers
ps[0]=zero
ps[1]=zero
an[0][0]=one
an[1][1]=one
else:
#
# SCALE THE DEVIATORIC STRESS COMPONENTS
#
ofac=one/fact
s11=ofac*s11
s22=ofac*s22
s12=ofac*s12
#
# GET THE EIGENVALUES AND EIGENVECTORS
#
temp=s11**2+s12**2
d=math.sqrt(temp)
ps[0]=sh-fact*d
ps[1]=sh+fact*d
s11=s11+d
s22=s22+d
if math.fabs(s11)>=math.fabs(s22):
ofac=one/math.sqrt(s11**2+s12**2)
an[0][0]=ofac*s12
an[1][0]=-ofac*s11
else:
ofac=one/math.sqrt(s12**2+s22**2)
an[0][0]=ofac*s22
an[1][0]=-ofac*s12
an[0][1]=-an[1][0]
an[1][1]=an[0][0]
if ndi==3:
ps[2]=s[2]
an[2][2]=one
else:
#
# THREE SHEAR COMPONENTS, ALL DIRECTIONS UNKNOWN
#
# GET DEVIATORIC STRESSES
#
if ndi==0:
sh=zero
s11=zero
s22=zero
s33=zero
elif ndi==1:
sh=third*s[0]
s11=s[0]-sh
s22=-sh
s33=-sh
elif ndi==2:
sh=third*(s[0]+s[1])
s11=s[0]-sh
s22=s[1]-sh
s33=-sh
else:
sh=third*(s[0]+s[1]+s[2])
s11=s[0]-sh
s22=s[1]-sh
s33=s[2]-sh
if lstr==1:
s12=s[ndip1]
s13=s[ndip2]
if nshr>2:
s23=s[ndip3]
else:
s23=zero
else:
s12=half*s[ndip1]
s13=half*s[ndip2]
if nshr>2:
s23=half*s(ndip3)
else:
s23=zero
facd=max(math.fabs(s11),math.fabs(s22),math.fabs(s33))
facs=max(math.fabs(s12),math.fabs(s13),math.fabs(s23))
fact=max(facd,facs)
if fact<=math.fabs(preciz*sh) or facs<preciz*facd:
#
# ESSENTIALLY UNIT MATRIX
#
for k1 in range(ndi):
ps[k1]=s[k1]
an[0][0]=one
an[1][1]=one
an[2][2]=one
elif fact<asmall :
# -- We've been given very small [zero] numbers
for k1 in range(ndi):
ps[k1]=zero
an[0][0]=one
an[1][1]=one
an[2][2]=one
else:
#
# SCALE THE DEVIATORIC STRESS COMPONENTS
#
ofac=one/fact
s11=ofac*s11
s22=ofac*s22
s33=ofac*s33
s12=ofac*s12
s13=ofac*s13
s23=ofac*s23
#
# DO THE EIGENVALUE CALCULATION FOR THE SCALED STRESSES
#
q=math.sqrt(third*(s12**2+s13**2+s23**2+
half*(s11**2+s22**2+s33**2)))
r=(half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2)
+s12*s13*s23)/q**3
if(r>=one-precis) :
cos1=-half
cos2=-half
cos3= one
elif(r<=precis-one) :
cos1=-one
cos2= half
cos3= half
else:
ang = third*math.acos(r)
cos1= math.cos(ang)
cos2= math.cos(ang+pi23)
cos3=-cos1-cos2
ps[0]=two*q*cos1
ps[1]=two*q*cos2
ps[2]=two*q*cos3
#
# SPECIAL CASES: DOUBLE EIGENVALUES. SELECT THE UNIQUE ONE (K1)
#
if(ps[0]==ps[1]) :
k1=2
elif(ps[1]==ps[2]) :
k1=0
else:
k1=1
#
# SUBTRACT SELECTED EIGENVALUE
#
s11=s11-ps[k1]
s22=s22-ps[k1]
s33=s33-ps[k1]
#
# FIND THE FIRST PRINCIPAL DIRECTION BY CROSS PRODUCT OF
TWO ROWS
# TRY WHICH ROWS GIVE A NON-ZERO RESULT
# K0 INDICATES WHICH ROWS WERE USED
#
k0=1
an[0][k1]=s22*s33-s23*s23
an[1][k1]=s23*s13-s12*s33
an[2][k1]=s12*s23-s22*s13
anorm=an[0][k1]**2+an[1][k1]**2+an[2][k1]**2
a1=s12*s23-s13*s22
a2=s13*s12-s11*s23
a3=s11*s22-s12*s12
anormt=a1**2+a2**2+a3**2
if(math.fabs(anormt)>math.fabs(anorm)) :
k0=3
an[0][k1]=a1
an[1][k1]=a2
an[2][k1]=a3
anorm=anormt
a1=s12*s33-s23*s13
a2=s13*s13-s33*s11
a3=s11*s23-s13*s12
anormt=a1**2+a2**2+a3**2
if(math.fabs(anormt)>math.fabs(anorm)) :
k0=2
an[0][k1]=a1
an[1][k1]=a2
an[2][k1]=a3
anorm=anormt
if(anorm==zero) :
k0=0
an[0][k1]=one
an[1][k1]=zero
an[2][k1]=zero
anorm=one
anorm=one/math.sqrt(anorm)
an[0][k1]=an[0][k1]*anorm
an[1][k1]=an[1][k1]*anorm
an[2][k1]=an[2][k1]*anorm
if(k1!=1 or ps[0]==ps[2]) :
#
# CASE OF DOUBLE EIGENVALUES: ONLY REQUIREMENT IS THEY
ARE NORMAL TO THE
# FIRST DIRECTION. FIRST SELECT SECOND AND THIRD
EIGENVALUE (K2 AND K3)
#
if(k1!=1) :
k2=2-k1
k3=1
else:
k2=0
k3=2
# PICK UP ROW WHICH IS GUARANTEED NON-ZERO OR GENERATE
UNIT X DIRECTION
#
if(k0==0) :
an[0][k2]=zero
an[1][k2]=one
an[2][k2]=zero
anorm=one
elif(k0!=2) :
an[0][k2]=s12
an[1][k2]=s22
an[2][k2]=s23
anorm=s12**2+s22**2+s23**2
else:
an[0][k2]=s11
an[1][k2]=s12
an[2][k2]=s13
anorm=s11**2+s12**2+s13**2
else:
#
# THREE SEPARATE EIGENVALUES: REPEAT THE PROCESS
FOR THE FIRST ONE
#
k2=0
k3=2
s11=s11+ps[k1]-ps[k2]
s22=s22+ps[k1]-ps[k2]
s33=s33+ps[k1]-ps[k2]
an[0][k2]=s22*s33-s23*s23
an[1][k2]=s23*s13-s12*s33
an[2][k2]=s12*s23-s22*s13
anorm=an[0][k2]**2+an[1][k2]**2+an[2][k2]**2
a1=s12*s23-s13*s22
a2=s13*s12-s11*s23
a3=s11*s22-s12*s12
anormt=a1**2+a2**2+a3**2
if(math.fabs(anormt)>math.fabs(anorm)) :
an[0][k2]=a1
an[1][k2]=a2
an[2][k2]=a3
anorm=anormt
a1=s12*s33-s23*s13
a2=s13*s13-s33*s11
a3=s11*s23-s13*s12
anormt=a1**2+a2**2+a3**2
if(math.fabs(anormt)>math.fabs(anorm)) :
an[0][k2]=a1
an[1][k2]=a2
an[2][k2]=a3
anorm=anormt
if(anorm==zero) :
an[0][k2]=zero
an[1][k2]=one
an[2][k2]=zero
anorm=one
anorm=one/math.sqrt(anorm)
an[0][k2]=an[0][k2]*anorm
an[1][k2]=an[1][k2]*anorm
an[2][k2]=an[2][k2]*anorm
#
# GET THE LAST ONE BY #ROSSING THE FIRST TWO
#
an[0][k3]=an[1][k1]*an[2][k2]-an[2][k1]*an[1][k2]
an[1][k3]=an[2][k1]*an[0][k2]-an[0][k1]*an[2][k2]
an[2][k3]=an[0][k1]*an[1][k2]-an[1][k1]*an[0][k2]
#
# SCALE UP EIGENVALUES AND ADD HYDROSTATIC PART
#
ps[0]=fact*ps[0]+sh
ps[1]=fact*ps[1]+sh
ps[2]=fact*ps[2]+sh
return ps,an
=============================================================================
-------------------------
Dave Lindeman
Lead Research Specialist
3M Company
3M Center 235-3F-08
St. Paul, MN 55144
651-733-6383
Dear Dave,
this answer is not accessible to users without support. If you have the code,
please share it for everyone's benefit.
Thank you in advance,
Frank
-------- Original-Nachricht --------
Datum: Tue, 09 Dec 2008 11:02:48 -0600
Betreff: Re: [Abaqus] Principal Strain Direction
Login to the "My Support" section of simulia.com, and look for Answer ID
1628 -- there are already Python scripts available for this.
Regards,
Dave
-------------------------
Dave Lindeman
Lead Research Specialist
3M Company
3M Center 235-3F-08
St. Paul, MN 55144
651-733-6383
Post by suganHi All,
Is it possible to obtain the angle of the principal strains? In the
visualization module, the direction of the principal strains can be
seen
Post by suganin the graphics window but I would like to get the numbers with
respect
Post by suganto the global axis.
Any help is greatly appreciated!
Thank you
[Non-text portions of this message have been removed]
--
----------------------------------------------------------
Frank Richter
Institute of Materials Science
Ruhr-Universitaet Bochum
Bochum
Germany
Sensationsangebot verlängert: GMX FreeDSL - Telefonanschluss + DSL
für nur 16,37 Euro/mtl.!* http://dsl.gmx.de/?ac=OM.AD.PD003K1308T4569a
<http://dsl.gmx.de/?ac=OM.AD.PD003K1308T4569a>