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:         dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
          dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
          dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
  
          rot = 0.0D0
          rot(1, 2)= 0.5D0*( Fbar(1, 2)-Fbar(2, 1) );  rot(2, 1) = -rot(1, 2)
  
negative.code:       !**
        if( matlaniso ) then
          do j = 1, 3
            estrain(j) = ALPO(j)*(TEMPC-ref_temp)-alpo0(j)*(TEMP0-ref_temp)
          end do
          estrain(4:6) = 0.0D0
          call transformation(coordsys, tm)
          estrain(:) = matmul( estrain(:), tm  )      ! to global coord
          estrain(4:6) = estrain(4:6)*2.0D0
        else
          THERMAL_EPS = ALP*(TEMPC-ref_temp)-alp0*(TEMP0-ref_temp)
          estrain(1:3) = THERMAL_EPS
          estrain(4:6) = 0.0D0
        end if
  