anchor.diff:        label = lbE
   
   c anisotropy of material numbers 1,2,3
  -      if (label.EQ.1 .or. label.EQ.2 .or.
  -     & label.EQ.3) then
  -
  -      Call matrRotate(XY1,XY2,XY3,XY4,label,
  -     & DATAFEM,DATAFEM1)
  +      if (label .EQ. 1 .or. 
  +     &    label .EQ. 2 .or.
  +     &    label .EQ. 3) then
  +
  +        Call matrRotate(XY1,XY2,XY3,XY4,label,
  +     &                  DATAFEM,DATAFEM1)
   
   c     A(1:4,1:4) is elemental vector elliptic operator;
   c     in other words, for the bilinear form <grad(P1), grad(P1)>
   
  -      Call fem3Dtet(XY1, XY2, XY3, XY4,
  +        Call fem3Dtet(XY1, XY2, XY3, XY4,
        &              GRAD, FEM_P1, GRAD, FEM_P1,
        &              label, Ddiff, DATAFEM1, iSYS, 2,
        &              LDA, A, ir, ic)
   
   c ... compute the right hand side vector using external function Drhs
   c     in other words, the linear form <Drhs(x), P1>
  -      Call fem3Dtet(XY1, XY2, XY3, XY4,
  +        Call fem3Dtet(XY1, XY2, XY3, XY4,
        &              IDEN, FEM_P0, IDEN, FEM_P1,
        &              lbE, Drhs, DATAFEM1, iSYS, 5,
        &              LDA, F, ir, ic)
   
   
  -        Else
  -
  -      Call fem3Dtet(XY1, XY2, XY3, XY4,
  +      Else
  +
  +        Call fem3Dtet(XY1, XY2, XY3, XY4,
        &              GRAD, FEM_P1, GRAD, FEM_P1,
        &              label, Ddiff, DATAFEM, iSYS, 2,
        &              LDA, A, ir, ic)
   
  -      Call fem3Dtet(XY1, XY2, XY3, XY4,
  +        Call fem3Dtet(XY1, XY2, XY3, XY4,
        &              IDEN, FEM_P0, IDEN, FEM_P1,
        &              lbE, Drhs, DATAFEM, iSYS, 5,
        &              LDA, F, ir, ic)
   
  -        End if
  +      End if
   
   
   c ... impose the Neumann boundary conditions
         Do i = 1, 3
  -         XYP(i, 1) = XY1(i)
  -         XYP(i, 2) = XY2(i)
  -         XYP(i, 3) = XY3(i)
  -         XYP(i, 4) = XY4(i)
  +        XYP(i, 1) = XY1(i)
  +        XYP(i, 2) = XY2(i)
  +        XYP(i, 3) = XY3(i)
  +        XYP(i, 4) = XY4(i)
         End do
   
         Do k = 1, 4
  -         If(lbF(k).GT.0) Then
  +        If(lbF(k) .GT. 0) Then
               l = iref(k + 1)
               m = iref(l + 1)
   
  
positive.code:      & -1D0,XY3(1),XY3(2),-1D0,D3)
  
        if (D.NE.0) then
        A=D1/D
        B=D2/D
        C=D3/D
        in=1
        else
        stop 'determinant=0'
        end if
  
        Else if (t1.EQ.0) then
  
negative.code: c ... impose Dirichlet boundary conditions at triangle nodes
  c     this condition should go the last one
        Do k = 1, 4
           If(lbP(k) .NE. 0) Then
              x = XYP(1, k)
              y = XYP(2, k)
              z = XYP(3, k)!
              ibc = Dbc(x, y, z, lbP(k), DATAFEM, iSYS, eBC)
  
              If(ibc .EQ. BC_DIRICHLET) Then
                 Call applyDIR(LDA, nRow, A, F, k, eBC(1))
              End if
           End if
  