Stress Projection form GPs to nodes

[A.Harandi]

I have a user elemnt in feap for which I tried to extrapolate strees from Gauss points to the nodes the plot is nonsense ,however i compare it with another user elemnets I have known it and the s and p are calculated right and also the stresses and strains comes to this sub routine rightly I wonder where is the problem.


      subroutine nSE_ProjStre(ix,xl,sig,eps,nes,nel,nqp,p,s,se)
     
      implicit none
     
     
      include 'strnum.h'
      include 'cdata.h'
     
      integer :: ix(*)
      integer :: nes,nel
      integer :: nqp
      double precision :: sig(nes,*)
      double precision :: eps(3,*)
      double precision :: p(*),det,detx
      double precision :: s(nen,*)
      double precision :: se(nen)
      double precision :: xsj,xg
      double precision :: qdp(nqp,2+1)
      double precision :: XIvec(2)
      integer, parameter :: ndm = 2
      double precision :: shp(nel,nen),shpx(3,nel),shpxi(3,nel)
      double precision :: xl(ndm,nen)
      integer :: i,iqp
     

     
     
      call nSetQuadPt(qdp)
      p(1:4)=0.0d0
      write(*,*)'p',p(1:4)
      do iqp = 1,nqp
          XIvec = qdp(iqp,1:ndm)
          call SHP2D(XIvec ,xl, shpx, detx, ndm, nen, ix, .false.) ! N(x,y), N derv x, N derv y
          do i = 1,nel
            xg = shpx(3,i)*detx
            p(i) = p(i) + xg
           
c!           Stress projections

            s(i,1) = s(i,1) + sig(1,iqp)*xg
            s(i,2) = s(i,2) + sig(2,iqp)*xg
            s(i,3) = s(i,3) + sig(3,iqp)*xg
           
c!           Strain projections
            s(i, 4) = s(i, 4) + eps(1,iqp)*xg
            s(i, 5) = s(i, 5) + eps(2,iqp)*xg
            s(i, 6) = s(i, 6) + eps(3,iqp)*xg

          end do
         
      end do
      write(*,*)'s',s(1,1:6)
      write(*,*)'p',p(1:4)
      iste=6   
     
      end subroutine
Its is worthy to mention that after calling the strees in the case of Isw==8 I call this subroutine with
         call nSE_ProjStre(ix,xl,spkv,glsv,nes,nel,nqp,p,s,se)

[Prof. S. Govindjee]

What version of FEAP are you using.

[Prof. R.L. Taylor]

I would suggest you print the values of the stress you pass to the routine.  I suspect you have a dimension error or an undefined values on the call.  Also check that you need all the variables you have defined.  If you want feap to compute and plot principal stress then sig(1:4) should be sig_11, sig_22, sig_33, sig_12 in s(i,1), s(i,2), s(i,3), and s(i,4).  I would do the same with strains so that iste = 8.  Also not you have to set istv in isw = 1 part of your element.  One is for dimensioning (istv) and the other to tell how many to process from this element type.

This should work in ver8.3 and later I think

[A.Harandi]

Dear Prof Taylor and Prof Govindjee
Thank you very much for your help.
I am working with version 8,4.
I have set the values and about my stresses which call them in this subroutine I am sure since I compare them with the other user element I have.
Also I put istv = max(istv,6) in my isw=1 case.
I have also tried to remove the variables which are not used in this subroutine but the result is the same.
I have also my s and p value after I calling this subroutine in the element main code which is fine and same as the expected value.
I am pretty sure that s and p are calculated correctly, however the plot is sth else.

[Prof. R.L. Taylor]

I would suggest you do a 1- element problem and carefully compare your element results to those from what feap element and the module slcn2d produce to see where the problem lies.  One thing in your module is setting p=0 -- you should not have to do this, it does not hurt but should not be necessary.  Make sure the value of 'nes' matches what you have in your element.

Ii assume also you are comparing for a 4-node quadrilateral or a 3-node triangle -- any other element may not work with diagonal projection.

[A.Harandi]

Dear Prof Taylor,
I have checked every thing
in the main element subroutine the s(8,8) sinc it is used for stiffness when I call my subroutine in isw==8 i mentioned the new dimension in the subroutine. so the size of s in the main subroutine is changed again how can I overcome this problem since s is also used for residual and tangent also.

[Prof. S. Govindjee]

The leading dimension of s( , ) should be nen for version 8.4:

Code: [Select]

real*8  s(nen,*)


[Prof. R.L. Taylor]

When we call an element routine with isw = 8, the item passed as  "s" to the element hass dimension s(nen,*); when we call with most other 'isw' it has dimension s(nst,*).  Thus, you need to dimension in your projection routine as Professor Govindjee remarked.