anchor.diff:          write (*,*) "LSolver has finished!"
   
   c ===   write results to VTK file
  -        call GMVscalarTet(sle_solution, 
  -     &                    "mesh.vtk", 
  +        call GMVscalarVRT(sle_solution, 
  +     &                    "mesh.gmv", 
        &                    10,
        &                    total_mesh_nodes,
        &                    node_coordinates, 
  
positive.code: 
                    ibc = Dbc(x, y, z, lbF(k), DATAFEM, iSYS, eBC)
                     
                    label = lbF(k)
  
                    Call fem3Dtri(XYP(1, k), XYP(1, l), XYP(1, m),
  
negative.code:               end if
            End do
  
  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
        End do
  
            Return
        End
  
  