Work in progress on parallel domain decomposition.

AMDiS/src/DataCollector.cc

@@ -6,6 +6,7 @@
 #include "Projection.h"
 
 namespace AMDiS {
+
   DataCollector::DataCollector(const FiniteElemSpace *feSpace,
 			       DOFVector<double> *values,
 			       int level,
@@ -61,20 +62,17 @@ namespace AMDiS {
 
   void DataCollector::fillAllData()
   {
-    if (!elementDataCollected_) {
+    if (!elementDataCollected_)
       startCollectingElementData();
-    }
 
-    if (!periodicDataCollected_) {
+    if (!periodicDataCollected_)
       startCollectingPeriodicData();
-    }
 
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
   }
 
-  int DataCollector::startCollectingElementData()
+  void DataCollector::startCollectingElementData()
   {
     FUNCNAME("DataCollector::startCollectingElementData()");
 
@@ -99,26 +97,24 @@ namespace AMDiS {
 
     // Traverse elements to create element information
     elInfo = stack.traverseFirst(mesh_, level_, flag);
+
     while (elInfo) {
-      if (!writeElem_ || writeElem_(elInfo)) {
+      if (!writeElem_ || writeElem_(elInfo))
 	addElementData(elInfo);
-      }
+
       elInfo = stack.traverseNext(elInfo);
     }
 
     elementDataCollected_ = true;
-
-    return(0);
   }
 
-  int DataCollector::startCollectingValueData()
+  void DataCollector::startCollectingValueData()
   {
     FUNCNAME("DataCollector::startCollectingValueData()");
 
     DOFVector<int>::Iterator intPointIt(interpPointInd_, USED_DOFS);
-    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
+    for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
       (*intPointIt) = -1;
-    }
 
     interpPoints_.clear();
 
@@ -130,13 +126,11 @@ namespace AMDiS {
 
     // Traverse elements to add value information and to mark
     // interpolation points.
-    ElInfo *elInfo = stack.traverseFirst(mesh_, 
-					 level_, 
+    ElInfo *elInfo = stack.traverseFirst(mesh_, level_, 
 					 traverseFlag_ | Mesh::FILL_COORDS);
     while (elInfo) {
-      if (!writeElem_ || writeElem_(elInfo)) {
+      if (!writeElem_ || writeElem_(elInfo))
 	addValueData(elInfo);
-      }
       elInfo = stack.traverseNext(elInfo);
     }
 
@@ -146,11 +140,9 @@ namespace AMDiS {
       // Remove all interpolation marks and, instead, set to each
       // interpolation point its continous index starting from 0.
       int i = 0;
-      for (intPointIt.reset(); !intPointIt.end(); ++intPointIt) {
-	if (*intPointIt == -3) {
+      for (intPointIt.reset(); !intPointIt.end(); ++intPointIt)
+	if (*intPointIt == -3)
 	  *intPointIt = i++;
-	}
-      }
       
       // Traverse elements to create interpolation values.
       elInfo = stack.traverseFirst(mesh_, level_, traverseFlag_ | Mesh::FILL_COORDS);
@@ -163,11 +155,9 @@ namespace AMDiS {
    
     delete [] localDOFs_;
     valueDataCollected_ = true;
-    
-    return 0;
   }
 
-  int DataCollector::startCollectingPeriodicData()
+  void DataCollector::startCollectingPeriodicData()
   {    
     FUNCNAME("DataCollector::startCollectingPeriodicData()");
 
@@ -210,11 +200,9 @@ namespace AMDiS {
     }   
 
     periodicDataCollected_ = true;
-
-    return(0);
   }
 
-  int DataCollector::addElementData(ElInfo* elInfo)
+  void DataCollector::addElementData(ElInfo* elInfo)
   {
     FUNCNAME("DataCollector::addElementData()");
 
@@ -280,19 +268,15 @@ namespace AMDiS {
 	outputIndices_[elInfo->getNeighbour(i)->getIndex()] :
 	-1;
     }
-
     
-    if (dim_ == 3) {
+    if (dim_ == 3)
       elementInfo.type = (dynamic_cast<ElInfo3d*>(elInfo)->getType());
-    }
 
     // remember element info
     elements_.push_back(elementInfo);
-
-    return(0);
   }
 
-  int DataCollector::addValueData(ElInfo *elInfo)
+  void DataCollector::addValueData(ElInfo *elInfo)
   {
     FUNCNAME("DataCollector::addValueData()");
     
@@ -345,28 +329,24 @@ namespace AMDiS {
 	  nInterpPoints_++;
 	}
       }      
-    }
-   
-    return(0);
+    }  
   }
 
-  int DataCollector::addInterpData(ElInfo *elInfo)
+  void DataCollector::addInterpData(ElInfo *elInfo)
   {
     FUNCNAME("DataCollector::addInterpData()");
     
     basisFcts_->getLocalIndices(elInfo->getElement(), localAdmin_, localDOFs_);
 
     std::vector<int> elemInterpPoints(0);
-    for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++) {
+    for (int i = mesh_->getGeo(VERTEX); i < nBasisFcts_; i++)
       elemInterpPoints.push_back((*interpPointInd_)[localDOFs_[i]]);
-    }
     
-    interpPoints_.push_back(elemInterpPoints);
-  
-    return(0);
+    interpPoints_.push_back(elemInterpPoints); 
   }
 
-  int DataCollector::addPeriodicData(ElInfo *elInfo) {
+  void DataCollector::addPeriodicData(ElInfo *elInfo) 
+  {
     FUNCNAME("DataCollector::addPeriodicData");
 
     LeafDataPeriodic *ldp = dynamic_cast<LeafDataPeriodic*>
@@ -423,69 +403,60 @@ namespace AMDiS {
 	}
       }
     }
-      
-    return(0);
   }
 
   std::list<ElementInfo>* DataCollector::getElementInfos()
   {        
-    if (!elementDataCollected_) {
+    if (!elementDataCollected_)
       startCollectingElementData();
-    }
       
     return &elements_;
   }
 
   DOFVector< std::list<VertexInfo> >* DataCollector::getVertexInfos()
   {
-    if (!elementDataCollected_) {
+    if (!elementDataCollected_)
       startCollectingElementData();
-    }
 
     return vertexInfos_;
   }
 
   std::list<PeriodicInfo>* DataCollector::getPeriodicInfos()
   {
-    if (!periodicDataCollected_) {
+    if (!periodicDataCollected_)
       startCollectingPeriodicData();
-    }
 
     return &periodicInfos_;
   }
 
   int DataCollector::getNumberVertices()
   {
-    if (!elementDataCollected_) {
+    if (!elementDataCollected_)
       startCollectingElementData();
-    }
 
     return nVertices_;
   }
 
   int DataCollector::getNumberElements()
   {
-    if (!elementDataCollected_) {
+    if (!elementDataCollected_)
       startCollectingElementData();
-    }
 
     return nElements_;
   }
 
   int DataCollector::getNumberInterpPoints()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return nInterpPoints_;
   }
   
   int DataCollector::getNumberConnections()
   {
-    if (!periodicDataCollected_) {
+    if (!periodicDataCollected_)
       startCollectingPeriodicData();
-    }
 
     return nConnections_;
   }
@@ -502,45 +473,40 @@ namespace AMDiS {
 
   DOFVector<double>* DataCollector::getValues()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return values_;
   }
 
   DOFVector< std::list<WorldVector<double> > >* DataCollector::getDofCoords()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return dofCoords_;
   }
 
   DOFVector<int>* DataCollector::getInterpPointInd()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return interpPointInd_;
   }
 
   DOFVector< std::list<WorldVector<double> > >* DataCollector::getInterpPointCoords()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return interpPointCoords_;
   }
 
   std::vector< std::vector<int> >* DataCollector::getInterpPoints()
   {
-    if (!valueDataCollected_) {
+    if (!valueDataCollected_)
       startCollectingValueData();
-    }
 
     return &interpPoints_;
   }

AMDiS/src/DataCollector.h

@@ -101,25 +101,25 @@ namespace AMDiS {
 
   protected:  
     /// Start collecting element and vertex data of the problem.
-    int startCollectingElementData();
+    void startCollectingElementData();
 
     /// Start collecting value data of the problem.
-    int startCollectingValueData();
+    void startCollectingValueData();
 
     /// Start collecting periodic data of the problem.
-    int startCollectingPeriodicData();
+    void startCollectingPeriodicData();
 
     /// Adds information about one element and its vertices.
-    int addElementData(ElInfo* elInfo);
+    void addElementData(ElInfo* elInfo);
 
     /// Adds value information of one element.
-    int addValueData(ElInfo *elInfo);
+    void addValueData(ElInfo *elInfo);
 
     /// Adds information about interpolation points of vertices.
-    int addInterpData(ElInfo *elInfo);
+    void addInterpData(ElInfo *elInfo);
 
     /// Adds value information of one element.
-    int addPeriodicData(ElInfo *elInfo);
+    void addPeriodicData(ElInfo *elInfo);
 
     /// Vector with vertex values
     DOFVector<double> *values_;

AMDiS/src/ParallelDomainProblem.cc

@@ -10,6 +10,9 @@
 #include "PartitionElementData.h"
 #include "DOFMatrix.h"
 #include "DOFVector.h"
+#include "VtkWriter.h"
+
+#include "petscksp.h"
 
 namespace AMDiS {
 
@@ -138,9 +141,8 @@ namespace AMDiS {
       PartitionElementData *partitionData = 
 	dynamic_cast<PartitionElementData*>
 	((*it)->getElement()->getElementData(PARTITION_ED));
-      if (partitionData->getPartitionStatus() != IN) {
+      if (partitionData->getPartitionStatus() != IN)
 	macrosToRemove.push_back(*it);
-      }
     }
 
     mesh->removeMacroElements(macrosToRemove);
@@ -151,17 +153,15 @@ namespace AMDiS {
     int *lOrder = (int*)(malloc(sizeof(int) * rankDofs.size()));
 
     for (std::vector<const DegreeOfFreedom*>::iterator it = rankDofs.begin();
-	 it != rankDofs.end(); ++it) {
+	 it != rankDofs.end(); ++it)
       gOrder[nRankDOFs++] = (*it)[0];
-    }
 
     int rstart = 0;
     MPI_Scan(&nRankDOFs, &rstart, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
     rstart -= nRankDOFs;
    
-    for (int i = 0; i < nRankDOFs; i++) {
+    for (int i = 0; i < nRankDOFs; i++)
       lOrder[i] = rstart + i;
-    }
 
     AOCreateBasic(PETSC_COMM_WORLD, nRankDOFs, gOrder, lOrder, &applicationOrdering);
     
@@ -275,6 +275,10 @@ namespace AMDiS {
     ierr = VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
     ierr = VecSetSizes(petscRhsVec, rankDofs.size(), partitionDofs.size());    
     ierr = VecSetType(petscRhsVec, VECMPI);
+
+    ierr = VecCreate(PETSC_COMM_WORLD, &petscSolVec);
+    ierr = VecSetSizes(petscSolVec, rankDofs.size(), partitionDofs.size());    
+    ierr = VecSetType(petscSolVec, VECMPI);
   }
 
   void ParallelDomainProblemBase::exitParallelization(AdaptInfo *adaptInfo)
@@ -285,14 +289,27 @@ namespace AMDiS {
   void ParallelDomainProblemBase::fillPetscMatrix(DOFMatrix *mat, 
 						  DOFVector<double> *vec)
   {
-    /*   DOFMatrix::Iterator rowIt(mat, USED_DOFS);
-    for (rowIt.reset(); !rowIt.end(); ++rowIt) {
-      for (int i = 0; i < static_cast<int>((*rowIt).size()); i++) {
-	if ((*rowIt)[i].col >= 0) {
-	  MatSetValues(petscMatrix, 1, &i, 1, &((*rowIt)[i].col), &((*rowIt)[i].entry), ADD_VALUES);
+    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
+    namespace traits= mtl::traits;
+    typedef DOFMatrix::base_matrix_type Matrix;
+
+    traits::row<Matrix>::type row(mat->getBaseMatrix());
+    traits::col<Matrix>::type col(mat->getBaseMatrix());
+    traits::const_value<Matrix>::type value(mat->getBaseMatrix());
+
+    typedef traits::range_generator<major, Matrix>::type cursor_type;
+    typedef traits::range_generator<nz, cursor_type>::type icursor_type;
+
+    std::cout.precision(10);
+    for (cursor_type cursor = begin<major>(mat->getBaseMatrix()), cend = end<major>(mat->getBaseMatrix()); cursor != cend; ++cursor)
+      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
+	if (value(*icursor) != 0.0) {
+	  int r = row(*icursor);
+	  int c = col(*icursor);
+	  double v = value(*icursor);
+	  MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES);
 	}
-      }
-    }
+
 
     MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
@@ -303,7 +320,38 @@ namespace AMDiS {
       double value = *dofIt;
 
       VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES);
-      }*/
+    }
+  }
+
+  void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
+  {
+    KSP ksp;
+    PC pc;