Added some catches for bad VTK files.

examples/ResampleWithDataset.cpp

         }
 
     } else {
-
         vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
             std::cout << "Reading" << argv[1] << std::endl;
             Reader->SetFileName(argv[1]);
         Resample->SetValidPointMaskArrayName("HomogeneousDirichlet");
         Resample->Update();
 
-    std::cout << *Resample << std::endl;
+    //std::cout << *Resample << std::endl;.
+
+   std::cout << "Resampled grid, ncells=" << Resample->GetOutput()->GetNumberOfCells()  
+             << " npoints=" << Resample->GetOutput()->GetNumberOfPoints() << std::endl; 
 
     vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
         Writer->SetInputData( Resample->GetOutput() );

examples/VTKEdgeGsphere.cpp

         uGrid->GetPointData()->SetScalars( G );
     }
 
+    std::cout << "Writing file with ncells=" << uGrid->GetNumberOfCells() << "\tnpoints=" << uGrid->GetNumberOfPoints() << std::endl;
+
     // Write out new file
     vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
         //Writer->SetInputConnection(Reader->GetOutputPort());
         Writer->SetInputData( uGrid );
         Writer->SetFileName( argv[4] );
+        Writer->Update();
         Writer->Write();
 
     //Reader->Delete();

src/FEM4EllipticPDE.cpp

     //      Method:  SetVTUGridFile
     //--------------------------------------------------------------------------------------
     void FEM4EllipticPDE::SetVTUGridFile ( std::string& fname  ) {
-        std::cout  << "Loading VTK .vtu file " << fname << std::endl;
+        std::cout  << "Loading VTK .vtu file " << fname;
         vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
         MeshReader->SetFileName(fname.c_str());
         MeshReader->Update();
         //vtkGrid->ShallowCopy( MeshReader->GetOutput() );
         vtkGrid = MeshReader->GetOutput();
         MeshReader->Delete();
+        std::cout  << " Finished! " << std::endl;
         return ;
     }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
 
     //      Method:  SetVTUGridFile
     //--------------------------------------------------------------------------------------
     void FEM4EllipticPDE::SetVTUGridFile ( char* fname  ) {
-        std::cout  << "Loading VTK .vtu file " << fname << std::endl;
+        std::cout  << "Loading VTK .vtu file " << fname;
         vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
         MeshReader->SetFileName(fname);
         MeshReader->Update();
         //vtkGrid->DeepCopy( MeshReader->GetOutput() );
         //vtkGrid->ShallowCopy( MeshReader->GetOutput() );
         vtkGrid = MeshReader->GetOutput();
-        MeshReader->Delete();
+        if (!vtkGrid->GetNumberOfCells()) {
+            throw std::runtime_error("In FEM4EllipticPDE::SetVTUGridFile NUMBER OF CELLS = 0");
+	}
+        if (!vtkGrid->GetNumberOfPoints()) {
+            throw std::runtime_error("In FEM4EllipticPDE::SetVTUGridFile NUMBER OF POINTS = 0");
+	}
+	MeshReader->Delete();
+        std::cout  << " Finished! " << std::endl;
         return ;
     }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
 

src/LinearMag.cpp

 //--------------------------------------------------------------------------------------
 void LinearMag::CalculateRHS ( const std::string& susName ) {
 
+    std::cout << "Calculating RHS...";
+    std::cout.flush();
+
     if ( !vtkGrid->GetCellData()->GetScalars(susName.c_str()) ) {
         std::string err("No cell data by name ");
         err.append(susName);
         throw std::runtime_error(err.c_str());
     }
 
+    if (!vtkGrid->GetNumberOfPoints()) {
+        throw std::runtime_error("Number of points zero in input grid!");
+    }
+
     vtkDoubleArray* G = vtkDoubleArray::New();
     G->SetNumberOfComponents(1);
     G->SetNumberOfTuples( vtkGrid->GetNumberOfPoints() );
     G->SetName("G");
     //g.resize(vtkGrid->GetNumberOfPoints());
     VectorXr GG = VectorXr::Zero( vtkGrid->GetNumberOfPoints() );
-
+    
     // Iterate over all the points or all of the cells?
     for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
         if ( vtkGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
             throw std::runtime_error("Non-tetrahedral mesh encountered!");
         }
-
         Real cellSus = vtkGrid->GetCellData()->GetScalars(susName.c_str())->GetTuple(ic)[0];
 
         Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
         */
         // Face 0, ID 0,1,2
         Eigen::Matrix<Real, 3, 2> CC = Eigen::Matrix<Real, 3, 2>::Ones() ;
-        {
-
             CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
             CC.col(2) = C.row(0).tail<3>() - C.row(2).tail<3>();
             Vector3r nhat = CC.col(1).cross(CC.col(2));
             GG(ID[0]) += flux;
             GG(ID[1]) += flux;
             GG(ID[2]) += flux;
-        }
-        // Face 1, ID 0,1,3
+	// Face 1, ID 0,1,3
         {
             CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
             CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
             GG(ID[1]) += flux;
             GG(ID[3]) += flux;
         }
-        // Face 1, ID 0,2,3
+        // Face 2, ID 0,2,3
         {
             CC.col(1) = C.row(0).tail<3>() - C.row(2).tail<3>();
             CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
             GG(ID[2]) += flux;
             GG(ID[3]) += flux;
         }
-        // Face 1, ID 1,2,3
+        // Face 3, ID 1,2,3
         {
             CC.col(1) = C.row(1).tail<3>() - C.row(2).tail<3>();
             CC.col(2) = C.row(1).tail<3>() - C.row(3).tail<3>();
 
     vtkGrid->GetPointData()->AddArray( G );
     vtkGrid->GetPointData()->SetScalars( G );
+    std::cout << "finished" << std::endl;
+
     return ;
 }		// -----  end of method LinearMag::CalculateRHS  -----