User element stress projection correct in feap84, wrong in feap85

[crepes]

Dear Feap team, Dear FEAP‌ users

I wrote a user element for linear elasticity (plane strain). It works fine in feap8.4: the stress plots are the same as the stress plots of the build-in linear elastic model. In feap8.5 only the displacements are identical. The stresses differ strongly. I wrote a subroutine inside the elmt*.f. I am calling it like

Code: [Select]

      elseif(isw.eq.8) then

      iste = 3

        call stcn09(ix,d,xl,ul,shp,hr(nph),hr(nph+numnp),
     &              ndf,ndm,nel,numnp)


      end if
and

Code: [Select]

      if( isw.eq.1 ) then

(...)

      pstyp = 2

(...)

The projection routine is

Code: [Select]

      subroutine stcn09(ix,d,xl,ul,shp,dt,st,ndf,ndm,nel,numnp)

      implicit none

      include 'qudshp.h'

      integer  ndf,ndm,numnp,nel
      integer  j,l,ll,ix(*)
      real*8   xsj,xg
      real*8   cc(3)
      real*8   dt(numnp), st(numnp,*), xl(ndm,*),shp(3,4),d(*)
      real*8   sigma(4), epsi(4),ul(ndf,*),sg(3,9)
      real*8   sigv(6),dd(6,6)

c Lumped projection routine

      call quadr2d(d,.true.)

      do l = 1,lint !quadrature points

c       Shape functions
        call interp2d(l, xl,ix, ndm,nel, .false.)

c       calculate stresses sigv and moduli dd
        call strej09(d,ul,ndf,nel,l,dd,sigv)

        do j = 1,nel
          ll=ix(j)
          xg   = shp2(3,j,l)*jac(l)
          dt(ll) = dt(ll) + xg
c         Stress projections
          st(ll,1) = st(ll,1) + sigv(1)*xg
          st(ll,2) = st(ll,2) + sigv(2)*xg
          st(ll,3) = st(ll,3) + sigv(4)*xg
        end do


      end do !l

      end
The routine which is calculating the stresses (strej09) should be fine, since I am using the same one for calculating the (correct) divergence term of the residual.

I assume hr(nph) and hr(nph+numnp) is different in feap8.5. Am I right?

Kind regards,

crepes

[Prof. R.L. Taylor]

The stress projection changed some time ago.  Version 8.5 separates plot of stresses by material types and allows different values at the interface between elements.  I suggest you look at the projection routines for the built in elements and follow how they are done.  The projections are done into arrays that use the element vector 'p' and array 's' locations, although for projection 's' has a different pointer eventually.  Look at 'slcn2d.f (called from the element module, say sld2d1.f.

This should correct your problem.  Note if you have multi-material models the stresses you get will not be the same as in ver8.4 -- usually only near the material interface.

[crepes]

Dear Professor Taylor,

Thank you very much for your help. I did some small changes based on your suggestion. Now it is working very well.

Kind regards,

crepes