anchor.diff:  
           do i = 1, total_mesh_nodes
   c ===   again, no dirichlet boundary condition
  -            node_material_labels(i) = 1
  +            node_material_labels(i) = 0
  +        end do            
  +        do i = 1, 1
  +c ===   again, no dirichlet boundary condition
  +            node_material_labels(face_nodes(1,i)) = 2
  +            node_material_labels(face_nodes(2,i)) = 2
  +            node_material_labels(face_nodes(3,i)) = 2
           end do            
   
   c === mesh is ready; calculate dirac function
  
positive.code:         write (*,*) "LSolver has finished!"
  
  c ===   write results to VTK file
          call GMVscalarTet(sle_solution, 
       &                    "mesh.vtk", 
       &                    10,
       &                    total_mesh_nodes,
       &                    node_coordinates, 
  
negative.code:           iSYS(1) = 1
            iSYS(2) = 1
  
            coef(1, 1) = thermo_coefficient
  
            Ddiff = TENSOR_SCALAR
  
  