anchor.diff:          i = getElasticType(gauss%pMaterial%mtype)
         endif
         
  -			if( i==0 ) then
  -          if( present(temperature) ) then
  -            call calElasticMatrix( matl, sectTYPE, matrix, temperature  )
  -          else
  -            call calElasticMatrix( matl, sectTYPE, matrix )
  -          endif
  -        elseif(  i==1 ) then
  -          if( present(temperature) ) then
  -            call calElasticMatrix_ortho( gauss%pMaterial, sectTYPE, cdsys, matrix, temperature )
  -          else
  -            call calElasticMatrix_ortho( gauss%pMaterial, sectTYPE, cdsys, matrix )
  -          endif
  +      if( i==0 ) then
  +        if( present(temperature) ) then
  +          call calElasticMatrix( matl, sectTYPE, matrix, temperature  )
           else
  -          print *, "Elasticity type", matl%mtype, "not supported"
  -          stop
  +          call calElasticMatrix( matl, sectTYPE, matrix )
           endif
  -
  -      elseif( matl%mtype==NEOHOOKE .or. matl%mtype==MOONEYRIVLIN ) then
  -        call calElasticMooneyRivlin( matl, sectType, cijkl, gauss%strain  )
  -        call mat_c2d( cijkl, matrix, sectType )
  -      elseif( matl%mtype==ARRUDABOYCE )  then
  -        call calElasticArrudaBoyce( matl, sectType, cijkl, gauss%strain )
  -        call mat_c2d( cijkl, matrix, sectType )
  -      elseif( matl%mtype==USERHYPERELASTIC )  then
  -        call uElasticMatrix( matl%variables(101:), gauss%strain, matrix )
  -      elseif( isElastoplastic(matl%mtype) )  then
  -        if( present( temperature ) ) then
  -          call calElastoPlasticMatrix( matl, sectType, gauss%stress,  &
  -	           gauss%istatus(1), gauss%fstatus, matrix, temperature  )
  +      elseif(  i==1 ) then
  +        if( present(temperature) ) then
  +          call calElasticMatrix_ortho( gauss%pMaterial, sectTYPE, cdsys, matrix, temperature )
           else
  -          call calElastoPlasticMatrix( matl, sectType, gauss%stress,  &
  -	         gauss%istatus(1), gauss%fstatus, matrix  )
  +          call calElasticMatrix_ortho( gauss%pMaterial, sectTYPE, cdsys, matrix )
           endif
  -      elseif( matl%mtype==USERMATERIAL ) then
  -        call uMatlMatrix( matl%name, matl%variables(101:), gauss%strain,  &
  -            gauss%stress, gauss%fstatus, matrix, dt, gauss%ttime )
  -      elseif( matl%mtype==NORTON ) then
  -        if( present( temperature ) ) then
  -          call iso_creep( matl, sectTYPE, gauss%stress, gauss%strain, gauss%fstatus,  &
  -		       gauss%plstrain, dt, gauss%ttime, matrix, temperature  )
  -        else
  -           call iso_creep( matl, sectTYPE, gauss%stress, gauss%strain, gauss%fstatus,  &
  -		       gauss%plstrain, dt, gauss%ttime, matrix  )
  -        endif
  -      else
  -        stop "Material type not supported!"
  -      endif
  +      else
  +        print *, "Elasticity type", matl%mtype, "not supported"
  +        stop
  +      endif
  +
  +    elseif( matl%mtype==NEOHOOKE .or. matl%mtype==MOONEYRIVLIN ) then
  +      call calElasticMooneyRivlin( matl, sectType, cijkl, gauss%strain  )
  +      call mat_c2d( cijkl, matrix, sectType )
  +    elseif( matl%mtype==ARRUDABOYCE )  then
  +      call calElasticArrudaBoyce( matl, sectType, cijkl, gauss%strain )
  +      call mat_c2d( cijkl, matrix, sectType )
  +    elseif( matl%mtype==USERHYPERELASTIC )  then
  +      call uElasticMatrix( matl%variables(101:), gauss%strain, matrix )
  +    elseif( isElastoplastic(matl%mtype) )  then
  +      if( present( temperature ) ) then
  +        call calElastoPlasticMatrix( matl, sectType, gauss%stress,  &
  +           gauss%istatus(1), gauss%fstatus, matrix, temperature  )
  +      else
  +        call calElastoPlasticMatrix( matl, sectType, gauss%stress,  &
  +         gauss%istatus(1), gauss%fstatus, matrix  )
  +      endif
  +    elseif( matl%mtype==USERMATERIAL ) then
  +      call uMatlMatrix( matl%name, matl%variables(101:), gauss%strain,  &
  +          gauss%stress, gauss%fstatus, matrix, dt, gauss%ttime )
  +    elseif( matl%mtype==NORTON ) then
  +      if( present( temperature ) ) then
  +        call iso_creep( matl, sectTYPE, gauss%stress, gauss%strain, gauss%fstatus,  &
  +         gauss%plstrain, dt, gauss%ttime, matrix, temperature  )
  +      else
  +         call iso_creep( matl, sectTYPE, gauss%stress, gauss%strain, gauss%fstatus,  &
  +         gauss%plstrain, dt, gauss%ttime, matrix  )
  +      endif
  +    else
  +      stop "Material type not supported!"
  +    endif
   
       end subroutine
   
  
positive.code:              strain, stress, gauss%fstatus, dt, gauss%ttime )
        end if
  
      end subroutine StressUpdate
  
  !> Transfer rank 4 constituive matrix to rank 2 form
  subroutine mat_c2d( cijkl, dij, itype )
    real(kind=kreal),   intent(in)  :: cijkl(3,3,3,3)
    real(kind=kreal),   intent(out) :: dij(6,6)
    integer,            intent(in)  :: itype
  
negative.code: 
  module m_MatMatrix
    use hecmw_util
    use mMaterial
    use mMechGauss
    use m_ElasticLinear
    use mHyperElastic
    use m_ElastoPlastic
    use mViscoElastic
    use mCreep
    use mUElastic
    use mUmat
  
    implicit none
  
    contains
  
  !> Fetch the nlgeom flag of the material
    integer function getNlgeomFlag( gauss )
    type( tGaussStatus ), intent(in) :: gauss      !> status of qudrature point
    getNlgeomFlag = gauss%pMaterial%nlgeom_flag
    end function
  
  !> Calculate constituive matrix
    subroutine MatlMatrix( gauss, sectType, matrix, dt, cdsys, temperature, isEp )
      type( tGaussStatus ), intent(in) :: gauss          !> status of qudrature point
      INTEGER, INTENT(IN)              :: sectType       !> plane strain/stress or 3D
      REAL(KIND=kreal), INTENT(OUT)    :: matrix(:,:)    !> constitutive matrix
      REAL(KIND=kreal), INTENT(IN)     :: dt             !> time increment
      REAL(kind=kreal), INTENT(IN)     :: cdsys(3,3)     !> material coordinate system
      REAL(KIND=kreal), INTENT(IN), optional  :: temperature   !> temperature
      INTEGER(KIND=kint), INTENT(IN), optional :: isEp
  
      integer :: i
      integer :: flag = 0
      real(kind=kreal)            :: cijkl(3,3,3,3)
      TYPE( tMaterial ), pointer  :: matl
      matl=>gauss%pMaterial
  
      if( present(isEp) )then
        if( isEp == 1 )flag = 1
      endif
  
      if( matl%mtype==USERELASTIC ) then
        call uElasticMatrix( matl%variables(101:), gauss%strain, matrix )
      elseif( isViscoelastic(matl%mtype) ) then
        if( present(temperature) ) then
          call calViscoelasticMatrix( matl, sectTYPE, dt, matrix, temperature )
        else
          call calViscoelasticMatrix( matl, sectTYPE, dt, matrix )
        endif
      elseif( isElastic(matl%mtype) .or. flag==1 ) then
        if(flag==1)then
          i = getElasticType(ELASTIC)
  