anchor.diff:        wg = getWeight(etype, LX)*det
         B(1:6, 1:nn*ndof) = 0.0D0
         do j = 1, nn
  -        B4 = ( Bbar(j, 1)-gderiv(j, 1) )/3.0D0
  -        B6 = ( Bbar(j, 2)-gderiv(j, 2) )/3.0D0
  -        B8 = ( Bbar(j, 3)-gderiv(j, 3) )/3.0D0
  -        B(1, 3*j-2) = gderiv(j, 1)+B4
  -        B(1, 3*j-1) = B6
  -        B(1, 3*j  ) = B8
  -
  -        B(2, 3*j-2) = B4
  -        B(2, 3*j-1) = gderiv(j, 2)+B6
  -        B(2, 3*j  ) = B8
  -
  -        B(3, 3*j-2) = B4
  -        B(3, 3*j-1) = B6
  -        B(3, 3*j  ) = gderiv(j, 3)+B8
  -
  -        B(4, 3*j-2) = gderiv(j, 2)
  -        B(4, 3*j-1) = gderiv(j, 1)
  -        B(5, 3*j-1) = gderiv(j, 3)
  -        B(5, 3*j  ) = gderiv(j, 2)
  -        B(6, 3*j-2) = gderiv(j, 3)
  -        B(6, 3*j  ) = gderiv(j, 1)
  +        B(1,3*j-2) = gderiv(j, 1)
  +        B(2,3*j-1) = gderiv(j, 2)
  +        B(3,3*j  ) = gderiv(j, 3)
  +        B(4,3*j-2) = gderiv(j, 2)
  +        B(4,3*j-1) = gderiv(j, 1)
  +        B(5,3*j-1) = gderiv(j, 3)
  +        B(5,3*j  ) = gderiv(j, 2)
  +        B(6,3*j-2) = gderiv(j, 3)
  +        B(6,3*j  ) = gderiv(j, 1)
  +      end do
  +
  +      do j = 1, nn
  +        Z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
  +        B(1,3*j-2:3*j) = B(1,3*j-2:3*j)+Z1(1:3)
  +        B(2,3*j-2:3*j) = B(2,3*j-2:3*j)+Z1(1:3)
  +        B(3,3*j-2:3*j) = B(3,3*j-2:3*j)+Z1(1:3)
         end do
   
         TEMPC = dot_product( H(1:nn),TT(1:nn) )
  
positive.code:     ecoord(2, :) = YY(:)
      ecoord(3, :) = ZZ(:)
  
      naturalCoord = 0.0D0
      call getGlobalDeriv(etype, nn, naturalcoord, ecoord, det, Bbar)
      call getShapeFunc( etype, naturalcoord, H(1:nn) )
      TEMPC0 = dot_product( H(1:nn), TT(1:nn) )
      TEMP00 = dot_product( H(1:nn), T0(1:nn) )
  
      ! LOOP FOR INTEGRATION POINTS
      do LX = 1, NumOfQuadPoints(etype)
  
negative.code:         dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
          dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
          dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
          dstrain(:) = dstrain(:)-EPSTH(:)
  
          gausses(LX)%strain(1:6) = dstrain(1:6)+EPSTH(:)
          gausses(LX)%stress(1:6) = matmul( D(1:6, 1:6), dstrain(1:6) )
  