Element Errors

[ad10]

Thanks for replying.

I am facing one issue here. I have formulated one 4-node simple linear elastic element. I am getting entirely different results for the following loop,

BATCh
  DT,,1
  LOOP,,10      ! Number of Time steps
    TIME
    !loop,,10    ! Number of NR iterations
    Tang,,1
    !NEXT
  NEXT
ENd

If I don't use a time loop means simply Tang,,1, both results (FEAP and user element) are the same results. I can't be able to find the answer. Your suggestion is beneficial. Another question is related to the formulation of residual p,

I have seen two distinct cases for this,
1.  by using stiffness
         do i =1,8
              p(i) = 0.d0
              do j = 1,8
                  p(i) = p(i) - s(i,j)*ul(j)
              end do
          end do
For this I able to get results only for Tang,,1 Not for others (time loop)
2. by calculating, tr(B)*stress, but for this case, both look the same!

[Prof. R.L. Taylor]

very hard to tell what the problem is. 

All elements need to be formed as a residual in the p(*) array and a tangent in the s(*,*) array.  A linear problem is just a special case.

If you include transients you have to use the ctan(3) values to pick up the time coefficient.

Post your element if you want a more specific answer.

[ad10]

Dear all,

I have formulated the S and p array in isw 3 as follows

Code: [Select]

                  do i =1,ndim*nnode
                      do j = 1,ndim*nnode
                          s(i,j) = s(i,j) + kdd(i,j)*factor
                      end do   
                  end do
                 
                  do i =1,nst
                      do j = 1,nst
                          p(i) = p(i) - s(i,j)*utmp(j)
                      end do
                  end do
where, utmp are the current step displacements.

It works gives when only one NR iteration (simply Tang,,1) is used. But, it diverges for more NR loops used! What should I do so I will get a convergence?

[Prof. S. Govindjee]

Helpful commands:

Code: [Select]

SHOW,S -- tangent matrix in last element
SHOW,P -- RHS contribution from last element
SHOW,UL -- displacement vectors for last elementyou can use these command to check that you have the correct values in your arrays (check that they are consistent in matlab or python).

There is also the command

Code: [Select]

NTAN,ELEM,<element #>which shows you the difference between a numerical tangent and your element tangent (though I am not certain that this feature is working correctly in the current version of the code).

[Prof. R.L. Taylor]

Post your elements we can look to see what you are doing that may cause the error

[Prof. S. Govindjee]

His element is attached to his last message: elmt03.f

[ad10]

Prof.,

I have checked for S and P.

In both cases, S is the same so I can say stiffness is calculations are right. But, p is different for both. In the case of displacement loading through edisp on edge, both results are the same! but results change for eforce/ force loading condition.

I have attached the S and P for both cases as follows.

[Prof. R.L. Taylor]

Look at utmp values.  They must be wrong?

[ad10]

Yeah, utmp values are wrong!

can you please tell me, what's I am missing here?

Thanks

[Prof. R.L. Taylor]

How did you set utmp?  Should come from ul array

[ad10]

Prof,

I have assigned the ul values to the utmp

Code: [Select]

!     Initializing the displacement vector
          utmp(1) = ul(1,1)
          utmp(2) = ul(2,1)
! ---------------------------------
          utmp(3) = ul(1,2)
          utmp(4) = ul(2,2)
! ---------------------------------
          utmp(5) = ul(1,3)
          utmp(6) = ul(2,3)
! ---------------------------------
          utmp(7) = ul(1,4)
          utmp(8) = ul(2,4)


Like this and directly called these to formulate the load vector.

Is this the right way?

[ad10]

My isw 3 in short is this,

Code: [Select]

if(isw.eq.3 .or.isw.eq.6) then   ! Compute residual/tangent
!     Material Properties

          emod = d(1)
          enu  = d(2)
          thk  = d(3)
          nnode = nen
          ndim  = ndm
          ngp = 2
!     Initializing the displacement vector
          utmp(1) = ul(1,1)
          utmp(2) = ul(2,1)
! ---------------------------------
          utmp(3) = ul(1,2)
          utmp(4) = ul(2,2)
! ---------------------------------
          utmp(5) = ul(1,3)
          utmp(6) = ul(2,3)
! ---------------------------------
          utmp(7) = ul(1,4)
          utmp(8) = ul(2,4)
         
          do i=1,nnode
              do j=1,ndim
                  xref(j,i) = xl(j,i)  !Reference coordinaates
              end do
          end do
!     Initializing s(nst,nst) and p(nst)           
          do i =1,ndim*nnode
              p(i)    = zero
          end do
                 
          do i =1,ndim*nnode
              do j =1,ndim*nnode
                  s(i,j) = zero
              end do
          end do
!     Setting the integration points         
          call int2d(ngp,lint,sg)
!     Integration Loop
          do l =1,lint
              ss(1)=sg(1,l)
              !write(*,*)'xi:  ',ss(1)
              ss(2)=sg(2,l)
              !write(*,*)'eta: ',ss(2)
              !write(*,*)'Step Integration:    ',l
              call kshape2d(shapef,dshape,ss(1),ss(2),ngp,nnode,ndim)
              !write(*,*)'shapefunction:      ',shapef

              call kjac2d(xjac,xjacinv,detxjac,ndim,nnode,shapef
     *            ,dshape,xref)
             
              call kstrain(stran,utmp,Bmat,ndim,nnode)



         
!     Taking material model into consideration
              call pzero(CmatVoigt2d,3*3)
                  CmatVoigt2D(1,1)=(emod/((one+enu)*(one-two*enu)))*
     *                              (one-enu)
                  CmatVoigt2D(1,2)=(emod/((one+enu)*(one-two*enu)))*enu
                  CmatVoigt2D(2,1)=CmatVoigt2D(1,2)
                  CmatVoigt2D(2,2)=CmatVoigt2D(1,1)
                  CmatVoigt2D(3,3)=(emod/((one+enu)*(one-two*enu)))
     *                             *((one-two*enu)/two)
                 
                  call pzero(stress,3)
                  do ii = 1,3
                      do jj =1,3
                          stress(ii) = stress(ii) + CmatVoigt2d(ii,jj)*
     *                                 stran(jj)
                      end do
                  end do
                 
                  !write(*,*)'CmatVoigt2D',CmatVoigt2D

                 
                 
                  factor = detxjac*sg(3,l)*thk
                 
                  call kBmat2d(Bmat,ndim,nnode,ngp,shapef
     *            ,dshape,xjac,xjacinv)
                  !write(*,*)'Bmat:        ',Bmat
                 
                  call kCompTrans(Bmat,ndim*nnode,3,BmatT)

                  call kmatrixtripl(kdd,ndim*nnode,ndim*nnode,BmatT,
     &             nnode*ndim,3,CmatVoigt2d,3,3,Bmat,3,ndim*nnode)
                  !write(*,*)'kdd',kdd
                  !write(*,*)'kdd:     ',kdd
                  do i =1,ndim*nnode
                      do j = 1,ndim*nnode
                          s(i,j) = s(i,j) + kdd(i,j)*factor
                      end do   
                  end do
                 
                  do i =1,nst
                      do j = 1,nst
                          p(i) = p(i) - s(i,j)*utmp(j)
                      end do
                  end do
                 
          end do
         
      endif


I have nothing added to other ISWs. apart from 1.

[Prof. R.L. Taylor]

Post the entire element for more help

[Prof. R.L. Taylor]

I see it sorry

[Prof. R.L. Taylor]

You have the computation of the p inside the quadrature loop, move it outside.
You do not need to initialize the p and s arrays, it is done before the call.