Work on pdd and some other small code refactorings.

AMDiS/src/AMDiS_fwd.h

@@ -52,6 +52,7 @@ namespace AMDiS {
   class ITL_BasePreconditioner;
   class LeafDataPeriodic;
   class LevelAdmin;
+  class Line;
   class MacroElement;
   class MacroInfo;
   class Marker;
@@ -80,9 +81,13 @@ namespace AMDiS {
   class RCNeighbourList;
   class RefinementManager;
   class RobinBC;
+  class SubElInfo;
+  class SurfaceOperator;
   class SMIAdapter;
   class SystemVector;
+  class Tetrahedron;
   class TraverseStack;
+  class Triangle;
   class VertexInfo;
   class VertexVector;
 

AMDiS/src/DOFMatrix.cc

@@ -222,6 +222,17 @@ namespace AMDiS {
 	  DegreeOfFreedom col = colIndices[j];
 	  double entry = elMat[i][j];
 
+	  /*
+	  if (MPI::COMM_WORLD.Get_rank() == 0 && (row >= 156 || col >= 156))
+	    std::cout << "PROB 0: " << row << " " << col << std::endl;
+	  if (MPI::COMM_WORLD.Get_rank() == 1 && (row >= 151 || col >= 151))
+	    std::cout << "PROB 1: " << row << " " << col << std::endl;
+	  if (MPI::COMM_WORLD.Get_rank() == 2 && (row >= 146 || col >= 146))
+	    std::cout << "PROB 2: " << row << " " << col << std::endl;
+	  if (MPI::COMM_WORLD.Get_rank() == 3 && (row >= 153 || col >= 153))
+	    std::cout << "PROB 3: " << row << " " << col << std::endl;
+	  */
+
 	  if (add)
 	    ins[row][col] += sign * entry;
 	  else

AMDiS/src/ParallelDomainProblem.cc

@@ -88,11 +88,9 @@ namespace AMDiS {
 
     // === Global refinements. ===
 
-//     refinementManager->globalRefine(mesh, 1);
+    refinementManager->globalRefine(mesh, 12);
 
-//     updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
-
-//     exit(0);
+    updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
 
     // === Create petsc matrix. ===
 
@@ -155,6 +153,8 @@ namespace AMDiS {
 
   void ParallelDomainProblemBase::solvePetscMatrix(DOFVector<double> *vec)
   {
+    clock_t t = clock();
+
     KSP ksp;
     PC pc;
 
@@ -220,7 +220,9 @@ namespace AMDiS {
     }
     
     for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
-      delete [] sendBuffers[i];    
+      delete [] sendBuffers[i];
+
+    std::cout << "SOLUTION = " << TIME_USED(t,clock()) << std::endl;
   }
 
 
@@ -513,8 +515,8 @@ namespace AMDiS {
 	   sendIt = sendNewDofs.begin();
 	 sendIt != sendNewDofs.end();
 	 ++sendIt, i++) {
-      int nSendDOFs = sendIt->second.size() * 2;
-      sendBuffers[i] = new int[nSendDOFs];
+      int nSendDofs = sendIt->second.size() * 2;
+      sendBuffers[i] = new int[nSendDofs];
       int c = 0;
       for (std::map<const DegreeOfFreedom*, DegreeOfFreedom>::iterator 
 	     dofIt = sendIt->second.begin();
@@ -526,7 +528,7 @@ namespace AMDiS {
 	sendDofs[sendIt->first].push_back(dofIt->first);
       }
 
-      mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_INT, sendIt->first, 0);
+      mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
     }
 
     i = 0;
@@ -602,32 +604,6 @@ namespace AMDiS {
 
       delete [] recvBuffers[i];
     }
-
-#if 1
-    // === Create local information from sendDofs and recvDofs
-
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   it = sendDofs.begin();
-	 it != sendDofs.end();
-	 ++it)
-      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
-	   dofIt != it->second.end();
-	   dofIt++) {
-	//	std::cout << "AND SET S " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
-	//	const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
-      }
-
-    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
-	   it = recvDofs.begin();
-	 it != recvDofs.end();
-	 ++it)
-      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = it->second.begin();
-	   dofIt != it->second.end();
-	   dofIt++) {
-	//	std::cout << "AND SET R " << (*dofIt)[0] << " = " << mapGlobalLocalDOFs[(*dofIt)[0]] << std::endl;
-	//	const_cast<DegreeOfFreedom*>(*dofIt)[0] = mapGlobalLocalDOFs[(*dofIt)[0]];
-      }
-#endif
   }
 
 
@@ -656,8 +632,8 @@ namespace AMDiS {
     // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
     // === rankDOFs to boundaryDOFs                                                ===
 
-    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDOFs;
-    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDOFs;
+    std::map<int, std::vector<const DegreeOfFreedom*> > sendNewDofs = sendDofs;
+    std::map<int, std::vector<const DegreeOfFreedom*> > recvNewDofs = recvDofs;
 
     for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
 	   myIntBoundary.boundary.begin();
@@ -699,7 +675,7 @@ namespace AMDiS {
 		   boundDOFs);
 	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
 	  newBoundaryDOFs[boundDOFs[i]] = mpiRank;
-	  sendNewDOFs[it->first].push_back(boundDOFs[i]);
+	  sendNewDofs[it->first].push_back(boundDOFs[i]);
 	}
       }
     }    
@@ -748,7 +724,7 @@ namespace AMDiS {
 	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
 	  rankDOFs.erase(boundDOFs[i]);
 	  newBoundaryDOFs[boundDOFs[i]] = it->first;
-	  recvNewDOFs[it->first].push_back(boundDOFs[i]);
+	  recvNewDofs[it->first].push_back(boundDOFs[i]);
 	}
       }
 
@@ -787,13 +763,72 @@ namespace AMDiS {
 
     // === Send new DOF indices. ===
 
-    for (int rank = 0; rank < mpiSize; rank++) {
-      if (rank == mpiRank) 
-	continue;
+    std::vector<int*> sendBuffers(sendNewDofs.size());
+    std::vector<int*> recvBuffers(recvNewDofs.size());
+
+    i = 0;
+    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
+	   sendIt = sendNewDofs.begin();
+	 sendIt != sendNewDofs.end(); 
+	 ++sendIt, i++) {
+      int nSendDofs = sendIt->second.size();
+      sendBuffers[i] = new int[nSendDofs];
+      int c = 0;
+      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = sendIt->second.begin();
+	   dofIt != sendIt->second.end();
+	   ++dofIt)
+	sendBuffers[i][c++] = (*dofIt)[0] + rstart;
 
+      mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0);
+    }
+
+    i = 0;
+    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
+	   recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end(); ++recvIt, i++) {
+      int nRecvDofs = recvIt->second.size();
+      recvBuffers[i] = new int[nRecvDofs];
+	
+      mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0);
+    }
+
+    mpiComm.Barrier();
+
+    i = 0;
+    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
+	   sendIt = sendNewDofs.begin();
+	 sendIt != sendNewDofs.end(); 
+	 ++sendIt)
+      delete [] sendBuffers[i++];
+
+    i = 0;
+    int newDofIndex = nRankDOFs;
+    for (std::map<int, std::vector<const DegreeOfFreedom*> >::iterator 
+	   recvIt = recvNewDofs.begin();
+	 recvIt != recvNewDofs.end(); ++recvIt) {
       
+      int j = 0;
+      for (std::vector<const DegreeOfFreedom*>::iterator dofIt = recvIt->second.begin();
+	   dofIt != recvIt->second.end();
+	   ++dofIt) {
+	const_cast<DegreeOfFreedom*>(*dofIt)[0] = newDofIndex;
+// 	if (mpiRank == 1) 
+// 	  std::cout << "SET TO " << newDofIndex << std::endl;
+
+	mapLocalGlobalDOFs[newDofIndex] = recvBuffers[i][j];
+	mapGlobalLocalDOFs[recvBuffers[i][j]] = newDofIndex;
+	isRankDOF[newDofIndex] = false;
+	newDofIndex++;
+	j++;
+      }
+
+      delete [] recvBuffers[i++];
     }
 
+    // === Update list of dofs that must be communicated for solution exchange. ===
+
+    sendDofs = sendNewDofs;
+    recvDofs = recvNewDofs;
   }
 
 

AMDiS/src/Recovery.h

@@ -31,10 +31,6 @@
 using namespace std;
 using namespace AMDiS;
 
-// ============================================================================
-// ===== class Monomial =======================================================
-// ============================================================================
-
 class Monomial : public
 BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> >
 {
@@ -46,9 +42,9 @@ public:
       exponent(expon)
   {
     degree_ = exponent[0];
-    for (int i=1; i<exponent.size(); i++)
+    for (int i = 1; i<exponent.size(); i++)
       degree_ += exponent[i];
-  };
+  }
 
   virtual ~Monomial() {};
 
@@ -60,39 +56,40 @@ public:
       result *= pow(y[i] - z[i], exponent[i]);
 
     return result;
-  };
+  }
 
 private:
   WorldVector<int> exponent;
 };
 
 
-// ============================================================================
-// ===== class RecoveryStructure ==============================================
-// ============================================================================
-
 class RecoveryStructure
 {
 public:     // construtor, destructor
   // constructor
   RecoveryStructure() :
     coords(NULL),
-    A(NULL), rec_uh(NULL), rec_grdUh(NULL),
+    A(NULL), 
+    rec_uh(NULL), 
+    rec_grdUh(NULL),
     neighbors(NULL)
-  {};
+  {}
 
   // Destructor?
   ~RecoveryStructure()
   {
     clear();
-  };
-
-  RecoveryStructure(const RecoveryStructure& rhs) :
-    coords(NULL),
-    A(NULL), rec_uh(NULL), rec_grdUh(NULL),
-    neighbors(NULL) {
+  }
+
+  RecoveryStructure(const RecoveryStructure& rhs) 
+    : coords(NULL),
+      A(NULL), 
+      rec_uh(NULL), 
+      rec_grdUh(NULL),
+      neighbors(NULL) 
+  {
     *this=rhs;
-  };
+  }
 
   RecoveryStructure& operator=(const RecoveryStructure& rhs);
 

AMDiS/src/RecoveryEstimator.h

@@ -27,161 +27,120 @@
 
 namespace AMDiS {
 
-  // ============================================================================
-  // ===== class RecoveryEstimator ============================================
-  // ============================================================================
-
-  /** \brief
-   * Error estimator using the recovery gradient of the finite element solution.
-   */
+  /// Error estimator using the recovery gradient of the finite element solution.
   class RecoveryEstimator : public Estimator
   {
   public:
-    /** \brief
-     * Creator class.
-     */
+    /// Creator class.
     class Creator : public EstimatorCreator
     {
     public:
-      Creator() : EstimatorCreator(), uh(NULL) {};
+      Creator() 
+	: EstimatorCreator(), 
+	  uh(NULL) 
+      {}
 
-      virtual ~Creator() {};
+      virtual ~Creator() {}
 
       inline void setSolution(DOFVector<double> *uh_)
       {
 	uh = uh_;
-      };
+      }
 
-      /** \brief
-       * Returns a new Estimator object.
-       */
+      /// Returns a new Estimator object.
       Estimator* create()
       {
 	return new RecoveryEstimator(name, uh, row);
-      };
+      }
 
     protected:
       DOFVector<double> *uh;
     };
 
-    /** \brief
-     * constructor
-     */
+    /// constructor
     RecoveryEstimator(std::string name, DOFVector<double> *uh, int r = -1);
 
-    /** \brief
-     * destructor.
-     */
-    virtual ~RecoveryEstimator() {};
+    /// destructor.
+    virtual ~RecoveryEstimator() {}
 
-    /** \brief
-     * implements \ref Estimator::estimate().
-     */
+    /// implements \ref Estimator::estimate().
     virtual double estimate(double timestep = 0.0);
 
 
-    /*----- Setting functions -----*/
-
-    // Sets uh.
+    /// Sets uh.
     inline void setUh(DOFVector<double> *uh)
     {
       uh_ = uh;
-    };
+    }
 
-    // Sets f.
+    /// Sets f.
     inline void setFct(AbstractFunction<double, WorldVector<double> > *fct)
     {
       f_vec = fct;
-    };
+    }
 
+    ///
     inline void setFct(AbstractFunction<double, double> *fct)
     {
       f_scal = fct;
-    };
+    }
 
-    // Sets auxiliar vector.
+    /// Sets auxiliar vector.
     inline void setAuxVec(DOFVector<double> *uh)
     {
       aux_vec = uh;
-    };
-
+    }
 
-    /*----- Getting functions -----*/
 
-    // Gets recovery gradient.
+    /// Gets recovery gradient.
     inline DOFVector<WorldVector<double> >* getRecGrd()
     {
       return rec_grd;
-    };
+    }
 
-    // Gets higher-order recovery solution.
+    /// Gets higher-order recovery solution.
     inline DOFVector<double>* getRecUh()
     {
       return rec_uh;
-    };
-
+    }
 
-  protected:     // protected members
 
-    /** \brief
-     * finite element solution
-     */
+  protected:   
+    /// finite element solution
     DOFVector<double> *uh_;
 
-    /** \brief
-     * absolute or relative error?
-     */
+    /// absolute or relative error?
     int relative_;
 
-    /** \brief
-     * constant for scaling the estimator
-     */
+    /// constant for scaling the estimator
     double C;
 
-    /** \brief
-     * recovery method
-     */
+    /// recovery method
     int method;
 
-    /** \brief
-     * Working finite element space
-     */
+    /// Working finite element space
     const FiniteElemSpace *feSpace;
 
-    /** \brief
-     * Degree of corresponding basic functions
-     */
+    /// Degree of corresponding basic functions
     int degree;
 
-    /** \brief
-     * Basis functions for recovery vector.
-     */
+    /// Basis functions for recovery vector.
     const BasisFunction *rec_basFcts;
 
-    /** \brief
-     * Recovery gradient
-     */
+    /// Recovery gradient
     DOFVector<WorldVector<double> > *rec_grd;
 
-    /** \brief
-     * Higher-order recovery solution
-     */
+    /// Higher-order recovery solution
     DOFVector<double> *rec_uh;
 
-    /** \brief
-     * Diffusion coefficient (for flux recovery)
-     */
+    /// Diffusion coefficient (for flux recovery)
     AbstractFunction<double, WorldVector<double> > *f_vec;
     AbstractFunction<double, double> *f_scal;
 
-    /** \brief
-     * auxiliar vector
-     */
+    /// auxiliar vector
     DOFVector<double> *aux_vec;
 
-    /** \brief
-     * Recovery structure.
-     */
+    /// Recovery structure.
     Recovery *rec_struct;
   };
 

AMDiS/src/RefinementManager.h

@@ -24,19 +24,10 @@
 
 #include "Global.h"
 #include "Flag.h"
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
-  class Element;
-  class Line;
-  class Triangle;
-  class Tetrahedron;
-  class Mesh;
-  class ElInfo;
-  class TraverseStack;
-  class RCNeighbourList;
-
-
   /** \ingroup Adaption 
    * \brief
    * Base class of RefinementManager1d, RefinementManager2d, RefinementManager3d.
@@ -60,7 +51,8 @@ namespace AMDiS {
      * Generates new coordinates on curved boundaries. Can be overriden by
      * sub classes if used.
      */
-    virtual void setNewCoords() {
+    virtual void setNewCoords() 
+    {
       FUNCNAME("RefinementManager::setNewCoords");
       ERROR_EXIT("called for base class!\n");
     }
@@ -78,11 +70,13 @@ namespace AMDiS {
     Flag globalRefine(Mesh *aMesh, int mark);
 
     /// Set \ref newCoords
-    inline void newCoord(bool nc) { 
+    inline void newCoord(bool nc) 
+    { 
       newCoords = nc; 
     }
 
-    inline bool newCoord() { 
+    inline bool newCoord() 
+    { 
       return newCoords; 
     }
 
@@ -90,21 +84,25 @@ namespace AMDiS {
      * Implements the refinement of el_info->el. Can be overriden by sub
      * classes if used.
      */
-    virtual ElInfo* refineFunction(ElInfo* ) {
+    virtual ElInfo* refineFunction(ElInfo*) 
+    {
       FUNCNAME("RefinementManager::refineFunction()");
       ERROR_EXIT("called for base class!\n");
       return NULL;
     }
 
-    inline void setMesh(Mesh *m) {
+    inline void setMesh(Mesh *m) 
+    {
       mesh = m;
     }
 
-    inline void setStack(TraverseStack *s) {
+    inline void setStack(TraverseStack *s) 
+    {
       stack = s; 
     }
 
-    inline TraverseStack *getStack() { 
+    inline TraverseStack *getStack() 
+    { 
       return stack; 
     }
 

AMDiS/src/RefinementManager1d.h

@@ -31,37 +31,25 @@ namespace AMDiS {
   class RefinementManager1d : public RefinementManager
   {
   public:
-    /** \brief
-     * Calls base class constructor.
-     */
+    /// Calls base class constructor.
     RefinementManager1d()
       : RefinementManager()
-    {};
+    {}
 
-    /** \brief
-     * destructor 
-     */
-    virtual ~RefinementManager1d() {};
+    /// destructor 
+    virtual ~RefinementManager1d() {}
 
-    /** \brief
-     * Implements RefinementManager::refineMesh.
-     */
+    /// Implements RefinementManager::refineMesh.
     Flag refineMesh(Mesh *aMesh);
 
-    /** \brief
-     * Implements RefinementManager::setNewCoords
-     */
+    /// Implements RefinementManager::setNewCoords
     void setNewCoords();
 
   protected:
-    /** \brief
-     * Used by refineMesh while mesh traversal
-     */  
+    /// Used by refineMesh while mesh traversal
     static int recursiveRefineFunction(ElInfo* el_info);