Added option for generating stiffness matrix and modified example3 to…

… not generate a stiffness matrix (uses imported matrix instead)

src/core/cuda/cutil.cu

@@ -165,50 +165,35 @@ void __global__ compute_ele_indices_kernel(IndexType* tri0, IndexType* tri1, Ind
 //  int a = meiyongde + 1;
 //}
 
-template<>
-void tetmesh2ell<Matrix_ell_d_CG>(TetMesh* meshPtr, Matrix_ell_d_CG &A_d, bool verbose)
+template<typename MeshType>
+void populateStiffnessMatrix(MeshType mesh, Matrix_ell_d_CG &stiffnessMatrix, int numVerts, bool verbose)
 {
-  typedef typename Matrix_ell_d_CG::value_type ValueType;
-  typedef typename Matrix_ell_d_CG::index_type IndexType;
+  int maxsize = 0;
+  int num_entries = 0;
   const int X = Matrix_ell_d_CG::invalid_index;
 
-  int nv = meshPtr->vertices.size();
-  int ne = meshPtr->tets.size();
-
-
-
-  meshPtr->need_neighbors();
-  for(int i = 0; i < nv; i++)
+  for(int i = 0; i < numVerts; i++)
   {
-    std::sort(meshPtr->neighbors[i].begin(), meshPtr->neighbors[i].end());
+    num_entries += (int)mesh->neighbors[i].size();
+    maxsize = std::max(maxsize, (int)mesh->neighbors[i].size());
   }
-
-  int maxsize = 0;
-  int num_entries = 0;
-  for(int i = 0; i < nv; i++)
-  {
-    num_entries += (int)meshPtr->neighbors[i].size();
-    maxsize = std::max(maxsize, (int)meshPtr->neighbors[i].size());
-  }
-  num_entries += nv;
+  num_entries += numVerts;
   maxsize += 1; // should include itself
+
   if( verbose )
-    std::cout << "Constructing Matrix_ell_h_CG A... ";
-  Matrix_ell_h_CG A(nv, nv, num_entries, maxsize, 32);
-  if( verbose )
-    std::cout << "done." << std::endl;
-  //A.resize(nv, nv, num_entries, maxsize, 32);
+    std::cout << "Constructing Matrix_ell_h_CG A";
+  Matrix_ell_h_CG A(numVerts, numVerts, num_entries, maxsize, 32);
   if( verbose )
-    std::cout << "Adding values to matrix A... ";
-  for(int i = 0; i < nv; i++)
+    std::cout << "Adding values to matrix A";
+  for(int i = 0; i < numVerts; i++)
   {
     A.column_indices(i, 0) = i;
     for(int j = 1; j < maxsize; j++)
     {
       A.values(i, j) = 0.0;
-      if(j < meshPtr->neighbors[i].size() + 1)
+      if(j < mesh->neighbors[i].size() + 1)
       {
-        A.column_indices(i, j) = meshPtr->neighbors[i][j - 1];
+        A.column_indices(i, j) = mesh->neighbors[i][j - 1];
       }
       else
       {
@@ -217,16 +202,27 @@ void tetmesh2ell<Matrix_ell_d_CG>(TetMesh* meshPtr, Matrix_ell_d_CG &A_d, bool v
     }
   }
   if( verbose )
-    std::cout << "done." << std::endl;
+    std::cout << "Copying A to device";
+  //A_d = Matrix_ell_d_CG(A);
+  stiffnessMatrix = A;
+}
 
+template<>
+void tetmesh2ell<Matrix_ell_d_CG>(TetMesh* meshPtr, Matrix_ell_d_CG &A_d,
+		                          bool generateStiffnessMatrix, bool verbose)
+{
+  int nv = meshPtr->vertices.size();
 
-  if( verbose )
-    std::cout << "Copying A to device... ";
-//  A_d = Matrix_ell_d_CG(A);
-  A_d = A;
-  if( verbose )
-    std::cout << "done." << std::endl;
-//	Matrix_ell_d_CG A_tmp = A;
+  meshPtr->need_neighbors();
+  for(int i = 0; i < nv; i++)
+  {
+    std::sort(meshPtr->neighbors[i].begin(), meshPtr->neighbors[i].end());
+  }
+
+  if( generateStiffnessMatrix )
+  {
+    populateStiffnessMatrix<TetMesh*>(meshPtr, A_d, nv, verbose);
+  }
 }
 
 template<>
@@ -295,56 +291,21 @@ void trimesh2csr<Matrix_d_CG>(TriMesh* meshPtr, Matrix_d_CG &A_d)
 }
 
 template<>
-void trimesh2ell<Matrix_ell_d_CG>(TriMesh* meshPtr, Matrix_ell_d_CG &A_d)
+void trimesh2ell<Matrix_ell_d_CG>(TriMesh* meshPtr, Matrix_ell_d_CG &A_d,
+		                          bool generateStiffnessMatrix, bool verbose)
 {
-  typedef typename Matrix_ell_d_CG::value_type ValueType;
-  typedef typename Matrix_ell_d_CG::index_type IndexType;
-  const int X = Matrix_ell_d_CG::invalid_index;
-
   int nv = meshPtr->vertices.size();
-  int ne = meshPtr->faces.size();
-
-
 
   meshPtr->need_neighbors();
   for(int i = 0; i < nv; i++)
   {
     std::sort(meshPtr->neighbors[i].begin(), meshPtr->neighbors[i].end());
   }
 
-  int maxsize = 0;
-  int num_entries = 0;
-  for(int i = 0; i < nv; i++)
-  {
-    num_entries += (int)meshPtr->neighbors[i].size();
-    maxsize = std::max(maxsize, (int)meshPtr->neighbors[i].size());
-  }
-  num_entries += nv;
-  maxsize += 1; // should include itself
-
-  Matrix_ell_h_CG A;
-  A.resize(nv, nv, num_entries, maxsize, 32);
-  for(int i = 0; i < nv; i++)
+  if( generateStiffnessMatrix )
   {
-    A.column_indices(i, 0) = i;
-    for(int j = 1; j < maxsize; j++)
-    {
-      A.values(i, j) = 0.0;
-      if(j < meshPtr->neighbors[i].size() + 1)
-      {
-        A.column_indices(i, j) = meshPtr->neighbors[i][j - 1];
-      }
-      else
-      {
-        A.column_indices(i, j) = X;
-      }
-    }
+    populateStiffnessMatrix<TriMesh*>(meshPtr, A_d, nv, verbose);
   }
-
-  A_d = A;
-
-  A.resize(0, 0, 0, 0);
-//  meshPtr->neighbors.clear();
 }
 
 template<typename IndexType, typename ValueType>

src/core/cuda/setup_solver.cu

@@ -62,7 +62,9 @@ void checkMatrixForValidContents(Matrix_ell_h* A_h, const bool verbose)
   }
 }
 
-void getMatrixFromMesh(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A_h, const bool verbose) {
+void getMatrixFromMesh(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A_h,
+		               const bool generateStiffnessMatrix, const bool verbose)
+{
   srand48(0);
   meshPtr->set_verbose(verbose);
   //type is triangle mesh
@@ -79,7 +81,7 @@ void getMatrixFromMesh(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A_h, con
   Matrix_ell_d_CG Aell_d;
   Vector_d_CG RHS(meshPtr->vertices.size(), 0.0);
   //generate the unit constant mesh stiffness matrix
-  trimesh2ell<Matrix_ell_d_CG >(meshPtr, Aell_d);
+  trimesh2ell<Matrix_ell_d_CG >(meshPtr, Aell_d, generateStiffnessMatrix, verbose);
   cudaThreadSynchronize();
   //assembly step
   fem2d->initializeWithTriMesh(meshPtr);
@@ -111,7 +113,9 @@ int setup_solver(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A_d,
   return 0;
 }
 
-void getMatrixFromMesh(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_h* A_h, const bool verbose) {
+void getMatrixFromMesh(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_h* A_h,
+		               const bool generateStiffnessMatrix, const bool verbose)
+{
   srand48(0);
   meshPtr->set_verbose(verbose);
   //type is tet mesh
@@ -128,7 +132,7 @@ void getMatrixFromMesh(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_h* A_h, con
   Matrix_ell_d_CG Aell_d;
   Vector_d_CG RHS(meshPtr->vertices.size(), 0.0);
   //generate the unit constant mesh stiffness matrix
-  tetmesh2ell<Matrix_ell_d_CG >(meshPtr, Aell_d, verbose);
+  tetmesh2ell<Matrix_ell_d_CG >(meshPtr, Aell_d, generateStiffnessMatrix, verbose);
   cudaThreadSynchronize();
   //assembly step
   fem3d->initializeWithTetMesh(meshPtr);

src/core/include/cutil.h

@@ -143,13 +143,13 @@ struct cudaCSRGraph
 };
 
 template<class Matrix>
-void trimesh2ell(TriMesh* meshPtr, Matrix &A);
+void trimesh2ell(TriMesh* meshPtr, Matrix &A, bool generateStiffnessMatrix, bool verbose);
 
 template<class Matrix>
 void trimesh2csr(TriMesh* meshPtr, Matrix &A);
 
 template<class Matrix>
-void tetmesh2ell(TetMesh* meshPtr, Matrix &A, bool verbose);
+void tetmesh2ell(TetMesh* meshPtr, Matrix &A, bool generateStiffnessMatrix, bool verbose);
 
 
 

src/core/include/setup_solver.h

@@ -31,8 +31,8 @@ int setup_solver(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_d* A,
 int setup_solver(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_d* A,
     Vector_h_CG* x_h, Vector_h_CG* b_h, const bool verbose);
 
-void getMatrixFromMesh(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_h* A, const bool verbose);
-void getMatrixFromMesh(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A, const bool verbose);
+void getMatrixFromMesh(AMG_Config& cfg, TetMesh* meshPtr, Matrix_ell_h* A, const bool generateStiffnessMatrix, const bool verbose);
+void getMatrixFromMesh(AMG_Config& cfg, TriMesh* meshPtr, Matrix_ell_h* A, const bool getMatrixFromMesh, const bool verbose);
 
 class SparseEntry_t {
   public:

src/examples/example1.cu

@@ -79,7 +79,7 @@ int main(int argc, char** argv)
   //The stiffness matrix A 
   Matrix_ell_h A_h;
   //get the basic stiffness matrix (constant) by creating the mesh matrix
-  getMatrixFromMesh(cfg, tetmeshPtr, &A_h, verbose);
+  getMatrixFromMesh(cfg, tetmeshPtr, &A_h, true, verbose);
   //intialize the b matrix to ones for now. TODO @DEBUG
   Vector_h_CG b_h(A_h.num_rows, 1.0);
   //The answer vector.

src/examples/example2.cu

@@ -35,7 +35,7 @@ int main(int argc, char** argv)
   //The stiffness matrix A 
   Matrix_ell_h A;
   //get the basic stiffness matrix (constant) by creating the mesh matrix
-  getMatrixFromMesh(cfg, tetmeshPtr, &A, verbose);
+  getMatrixFromMesh(cfg, tetmeshPtr, &A, true, verbose);
   //intialize the b matrix to ones for now. TODO @DEBUG
   Vector_h_CG b_h(A.num_rows, 1.0);
   //The answer vector.

src/examples/example3.cu

@@ -150,27 +150,28 @@ int main(int argc, char** argv)
       (filename + ".node").c_str(),
       (filename + ".ele").c_str(), zero_based, verbose);
   //The stiffness matrix A 
-  Matrix_ell_h A_h;
+//  Matrix_ell_h A_h;
   //get the basic stiffness matrix (constant) by creating the mesh matrix
-  getMatrixFromMesh(cfg, tetmeshPtr, &A_h, verbose);
+//  getMatrixFromMesh(cfg, tetmeshPtr, &A_h, false, verbose);
   //intialize the b matrix to ones for now. TODO @DEBUG
-  Vector_h_CG b_h(A_h.num_rows, 1.0);
+//  Vector_h_CG b_h(A_h.num_rows, 1.0);
   //The answer vector.
-  Vector_h_CG x_h(A_h.num_rows, 0.0); //intial X vector
+//  Vector_h_CG x_h(A_h.num_rows, 0.0); //intial X vector
   //************************ DEBUG : creating identity matrix for stiffness properties for now.
-  Matrix_ell_h identity(A_h.num_rows, A_h.num_rows, A_h.num_rows, 1);
-  for (int i = 0; i < A_h.num_rows; i++) {
-    identity.column_indices(i, 0) = i;
-    identity.values(i, 0) = 1;
-  }
+//  Matrix_ell_h identity(A_h.num_rows, A_h.num_rows, A_h.num_rows, 1);
+//  for (int i = 0; i < A_h.num_rows; i++) {
+//    identity.column_indices(i, 0) = i;
+//    identity.values(i, 0) = 1;
+//  }
   //multiply the mesh matrix by the stiffness properties matrix.
   vector<double> test;
-  Matrix_ell_h out, A_fromFile;
-  Matrix_ell_h my_A(A_h);
-  cusp::multiply(identity, my_A, out);
-  A_h = Matrix_ell_h(out);
+//  Matrix_ell_h out, A_fromFile;
+//  Matrix_ell_h my_A(A_h);
+//  cusp::multiply(identity, my_A, out);
+//  A_h = Matrix_ell_h(out);
 
   //Import right-hand-side single-column array (b)
+  Vector_h_CG b_h;
   if( importRhsVectorFromFile(bFilename, &test, b_h, verbose) < 0 )
 	  return 0;
 
@@ -179,6 +180,9 @@ int main(int argc, char** argv)
   if( importStiffnessMatrixFromFile(aFilename, &A_h_imported, verbose) < 0 )
 	  return 0;
 
+  Vector_h_CG x_h(A_h_imported.num_rows, 0.0); //intial X vector
+  getMatrixFromMesh(cfg, tetmeshPtr, &A_h_imported, false, verbose);
+
   if( verbose )
     std::cout << "Calling setup_solver." << std::endl;
   //The final call to the solver