anchor.diff:          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(:)
   
           rot = 0.0D0
           rot(1, 2)= 0.5D0*( Fbar(1, 2)-Fbar(2, 1) );  rot(2, 1) = -rot(1, 2)
  
positive.code:       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)
        end do
  
        TEMPC = dot_product( H(1:nn),TT(1:nn) )
  
negative.code:     !---------------------------------------------------------------------
  
      real(kind=kreal) :: ALP, alp0, D(6, 6), B(6, ndof*nn)
      real(kind=kreal) :: Z1(3), det, ecoord(3, nn)
      integer(kind=kint) :: j, LX, serr
      real(kind=kreal) :: estrain(6), SGM(6), H(nn)
      real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
      real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
      real(kind=kreal) :: TEMPC, TEMP0, V0, jacob_ave, THERMAL_EPS, tm(6,6)
      logical :: ierr, matlaniso
  
      !---------------------------------------------------------------------
  