Principal Strain Direction

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

Principal Strain Direction

sugan-2
Hi 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 in the graphics window but I would like to get the numbers with respect to the global axis.
   
  Any help is greatly appreciated!
   
  Thank you

       

[Non-text portions of this message have been removed]

Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

Benz-2
The answer is easy : it's not an available
output........................................ And that's
incredible.................

benz


--- In [hidden email], sugan <sugan_1978@...> wrote:
>
> Hi 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 in the graphics window but I would like to get the numbers
with respect to the global axis.
>    
>   Any help is greatly appreciated!
>    
>   Thank you
>
>        
>
> [Non-text portions of this message have been removed]
>


Reply | Threaded
Open this post in threaded view
|

Re: Re: Principal Strain Direction

sugan-2
Hi everyone,
   
  Would it be possible to calculate the directions of the principal strains using an ABAQUS script? If so, would anyone be having a script to do the same? I am not familiar with writing scripts and would greatly appreciate if there is a template available .
   
  Thanks!

BenZ <[hidden email]> wrote:
          The answer is easy : it's not an available
output........................................ And that's
incredible.................

benz

--- In [hidden email], sugan <sugan_1978@...> wrote:
>
> Hi 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 in the graphics window but I would like to get the numbers
with respect to the global axis.
>
> Any help is greatly appreciated!
>
> Thank you
>
>
>
> [Non-text portions of this message have been removed]
>



                           

       

[Non-text portions of this message have been removed]

Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

Dave Lindeman
In reply to this post by sugan-2
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


sugan wrote:

>
>
> Hi 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
> in the graphics window but I would like to get the numbers with respect
> to the global axis.
>
> Any help is greatly appreciated!
>
> Thank you
>
> [Non-text portions of this message have been removed]
>
>
Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

Frank Richter-2

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
> Von: Dave Lindeman <[hidden email]>
> An: [hidden email]
> 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
>
>
> sugan wrote:
> >
> >
> > Hi 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
> > in the graphics window but I would like to get the numbers with respect
> > to 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
Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

Dave Lindeman
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.

A typical screenshot for the use of frin.py is shown below:

 >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


Frank Richter wrote:

>
>
>
> 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
>  > Von: Dave Lindeman <[hidden email] <mailto:ddlindeman%40mmm.com>>
>  > An: [hidden email] <mailto:Abaqus%40yahoogroups.com>
>  > 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
>  >
>  >
>  > sugan wrote:
>  > >
>  > >
>  > > Hi 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
>  > > in the graphics window but I would like to get the numbers with
> respect
>  > > to 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>
>
>
Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

sugan-2
In reply to this post by Dave Lindeman
Thanks for the info Dave! That was of great help.
Regards,
Sai

Dave Lindeman <[hidden email]> wrote:                             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
 
 sugan wrote:
 >
 >
 > Hi 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
 > in the graphics window but I would like to get the numbers with respect
 > to the global axis.
 >
 > Any help is greatly appreciated!
 >
 > Thank you
 >
 > [Non-text portions of this message have been removed]
 >
 >
 
     
                                       

       

[Non-text portions of this message have been removed]

Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

cpolindara
In reply to this post by Dave Lindeman
Hi Dave, thank you for sharing this information. It has been really useful. I do have one last question nonetheless... I've been looking into simula.com>Support>MySupport using the AnswerID as you suggested, but no results are shown for 1628. When looking into Abaqus documentation for sprind.f no  results are shown either, and when searching in google all I get is the link to this forum ¿could you tell me please where can I find that fortran subroutine?

Thanks a lot,

César.

Dave Lindeman wrote
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.

A typical screenshot for the use of frin.py is shown below:

 >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


Frank Richter wrote:
>
>
>
> 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
>  > Von: Dave Lindeman <ddlindeman@mmm.com <mailto:ddlindeman%40mmm.com>>
>  > An: Abaqus@yahoogroups.com <mailto:Abaqus%40yahoogroups.com>
>  > 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
>  >
>  >
>  > sugan wrote:
>  > >
>  > >
>  > > Hi 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
>  > > in the graphics window but I would like to get the numbers with
> respect
>  > > to 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>
>
>
Reply | Threaded
Open this post in threaded view
|

Re: Principal Strain Direction

AGGOUN
This post has NOT been accepted by the mailing list yet.
In reply to this post by Dave Lindeman
Hello sir. David,
In fact, I seek how to determine the main directions of the constraint, so I used the script you applied for long time, the problem  Is that I have an error message when I execute it, telling me that "No module named Sprind", do you have any idea where this message is, knowing that I tried to download the Sprind module, but it does not exist. I thank you in advance for your response.