Work on pdd

AMDiS/src/DOFMatrix.cc

@@ -14,7 +14,6 @@
 #include "BoundaryCondition.h"
 #include "BoundaryManager.h"
 #include "Assembler.h"
-#include "Utilities.h"
 
 namespace AMDiS {
 

AMDiS/src/ParallelDomainProblem.cc

@@ -276,6 +276,10 @@ namespace AMDiS {
   void ParallelDomainProblemBase::createInteriorBoundaryInfo(std::vector<const DegreeOfFreedom*>& rankDOFs,
 							     std::map<const DegreeOfFreedom*, int>& boundaryDOFs)
   {
+    FUNCNAME("ParallelDomainProblemBase::createInteriorBoundaryInfo()");
+
+    // === First, create all the information about the interior boundaries. ===
+
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
     while (elInfo) {
@@ -323,20 +327,86 @@ namespace AMDiS {
 	    /// === And add the part of the interior boundary. ===
 
 	    AtomicBoundary& bound = 
-	      interiorBoundary.getNewAtomicBoundary(ranksBoundary ? mpiRank :
-						    partitionVec[elInfo->getNeighbour(i)->getIndex()]);
+	      (ranksBoundary ?
+	       myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) :
+	       otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]));
+
 	    bound.rankObject.el = element;
 	    bound.rankObject.subObjAtBoundary = EDGE;
 	    bound.rankObject.ithObjAtBoundary = i;
 	    bound.neighbourObject.el = elInfo->getNeighbour(i);
 	    bound.neighbourObject.subObjAtBoundary = EDGE;
 	    bound.neighbourObject.ithObjAtBoundary = -1;
+	    std::cout << "ADD IN " << mpiRank << ": " << element->getIndex() << " " << elInfo->getNeighbour(i)->getIndex() << std::endl;
  	  }
 	}
       }
 
       elInfo = stack.traverseNext(elInfo);
     }
+
+
+    // === Once we have this information, we must care about their order. Eventually ===
+    // === all the boundaries have to be in the same order on both ranks (because    ===
+    // === each boundary is shared by exactly two ranks).                            ===
+
+    std::vector<int*> sendBuffers;
+    std::vector<int*> recvBuffers;
+
+    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = myIntBoundary.boundary.begin();
+	 rankIt != myIntBoundary.boundary.end();
+	 ++rankIt) {
+      int* buffer = new int[rankIt->second.size()];
+      for (int i = 0; i < static_cast<int>(rankIt->second.size()); i++)
+	buffer[i] = (rankIt->second)[i].rankObject.el->getIndex();
+      sendBuffers.push_back(buffer);
+      
+      mpiComm.Isend(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);
+    }
+
+    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
+	 rankIt != otherIntBoundary.boundary.end();
+	 ++rankIt) {
+      int *buffer = new int[rankIt->second.size()];
+      recvBuffers.push_back(buffer);
+      
+      mpiComm.Irecv(buffer, rankIt->second.size(), MPI_INT, rankIt->first, 0);      
+    }
+
+    mpiComm.Barrier();
+
+    int i = 0;
+    for (std::map<int, std::vector<AtomicBoundary> >::iterator rankIt = otherIntBoundary.boundary.begin();
+	 rankIt != otherIntBoundary.boundary.end();
+	 ++rankIt) {
+
+      // === We have received from rank "rankIt->first" the ordered list of element ===
+      // === indices. We now have to sort the corresponding list in this rank to    ===
+      // === get the same order.                                                    ===
+      
+      for (int j = 0; j < static_cast<int>(rankIt->second.size()); j++) {
+	// If the expected object is not at place, search for it.
+	if ((rankIt->second)[j].neighbourObject.el->getIndex() != recvBuffers[i][j]) {
+	  int k = j + 1;
+	  for (; k < static_cast<int>(rankIt->second.size()); k++)
+	    if ((rankIt->second)[k].neighbourObject.el->getIndex() == recvBuffers[i][j])
+	      break;
+
+	  // The element must always be found, because the list is just in another order.
+	  TEST_EXIT(k < rankIt->second.size())("Should never happen!\n");
+
+	  // Swap the current with the found element.
+	  AtomicBoundary tmpBound = (rankIt->second)[k];
+	  (rankIt->second)[k] = (rankIt->second)[j];
+	  (rankIt->second)[j] = tmpBound;	
+	}
+      }
+
+      delete [] recvBuffers[i++];
+    }
+
+    for (int i = 0; i < static_cast<int>(sendBuffers.size()); i++)
+      delete [] sendBuffers[i];
   }
 
 
@@ -527,8 +597,10 @@ namespace AMDiS {
   void ParallelDomainProblemBase::updateLocalGlobalNumbering(int& nRankDOFs, 
 							     int& nOverallDOFs)
   {
+    FUNCNAME("ParallelDomainProblemBase::updateLocalGlobalNumbering()");
+
     std::set<const DegreeOfFreedom*> rankDOFs;
-    std::set<const DegreeOfFreedom*> boundaryDOFs;
+    std::map<const DegreeOfFreedom*, int> boundaryDOFs;
 
     
     /// === Get all DOFs in ranks partition. ===
@@ -545,20 +617,105 @@ namespace AMDiS {
     }
 
 
-    // === Traverse on interior boundaries and move all not ranked owned DOFs from
-    //     rankDOFs to boundaryDOFs === //
+    // === Traverse on interior boundaries and move all not ranked owned DOFs from ===
+    // === rankDOFs to boundaryDOFs                                                ===
 
     for (std::map<int, std::vector<AtomicBoundary> >::iterator it = 
-	   interiorBoundary.boundary.begin();
-	 it != interiorBoundary.boundary.end();
+	   myIntBoundary.boundary.begin();
+	 it != myIntBoundary.boundary.end();
 	 ++it) {
-      for (int i = 0; i < it->second.size(); i++) {
-	
+      for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin();
+	   boundIt != it->second.end();
+	   ++boundIt) {
+	const DegreeOfFreedom *dof1 = NULL;
+	const DegreeOfFreedom *dof2 = NULL;
+
+	switch (boundIt->rankObject.ithObjAtBoundary) {
+	case 0:
+	  dof1 = boundIt->rankObject.el->getDOF(1);
+	  dof2 = boundIt->rankObject.el->getDOF(2);
+	  break;
+	case 1:
+	  dof1 = boundIt->rankObject.el->getDOF(0);
+	  dof2 = boundIt->rankObject.el->getDOF(2);
+	  break;
+	case 2:
+	  dof1 = boundIt->rankObject.el->getDOF(0);
+	  dof2 = boundIt->rankObject.el->getDOF(1);
+	  break;
+	default:
+	  ERROR_EXIT("Should never happen!\n");
+	}
+
+	std::vector<const DegreeOfFreedom*> boundDOFs;
+	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
+		   boundDOFs);
       }
     }    
   }
 
 
+  void ParallelDomainProblemBase::addAllDOFs(Element *el, int ithEdge, 
+					     std::vector<const DegreeOfFreedom*>& dofs,
+					     bool addVertices)
+  {
+    FUNCNAME("ParallelDomainProblemBase::addAllDOFs()");
+
+    const DegreeOfFreedom* boundDOF1 = NULL;
+    const DegreeOfFreedom* boundDOF2 = NULL;
+
+    if (addVertices) {
+      switch (ithEdge) {
+      case 0:
+	boundDOF1 = el->getDOF(1);
+	boundDOF2 = el->getDOF(2);
+	break;
+      case 1:
+	boundDOF1 = el->getDOF(0);
+	boundDOF2 = el->getDOF(2);
+	break;
+      case 2:
+	boundDOF1 = el->getDOF(0);
+	boundDOF2 = el->getDOF(1);
+	break;
+      default:
+	ERROR_EXIT("Should never happen!\n");
+      }
+
+      dofs.push_back(boundDOF1);
+    }
+
+    switch (ithEdge) {
+    case 0:
+      if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) {
+	addAllDOFs(el->getSecondChild()->getFirstChild(), 0, dofs, false);
+	dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2));
+	addAllDOFs(el->getSecondChild()->getSecondChild(), 1, dofs, false);
+      }
+      break;
+    case 1:
+      if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) {
+	addAllDOFs(el->getFirstChild()->getFirstChild(), 0, dofs, false);
+	dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2));
+	addAllDOFs(el->getFirstChild()->getSecondChild(), 1, dofs, false);
+      }
+      break;
+    case 2:      
+      if (el->getFirstChild()) {
+	addAllDOFs(el->getFirstChild(), 0, dofs, false);
+	dofs.push_back(el->getFirstChild()->getDOF(2));
+	addAllDOFs(el->getSecondChild(), 1, dofs, false);
+      }
+      break;      
+    default:
+      ERROR_EXIT("Should never happen!\n");
+    }
+
+    if (addVertices)
+      dofs.push_back(boundDOF2);		   
+  }
+
+
   void ParallelDomainProblemBase::createDOFMemberInfo(
 		       std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
 		       std::vector<const DegreeOfFreedom*>& rankDOFs,

AMDiS/src/ParallelDomainProblem.h

@@ -137,9 +137,15 @@ namespace AMDiS {
 
     void updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs);
 
+    void addAllDOFs(Element *el, int ithEdge, 
+		    std::vector<const DegreeOfFreedom*>& dofs,
+		    bool addVertices = true);
+
     /** \brief
      * This function traverses the whole mesh, i.e. before it is really partitioned,
-     * and collects information about which DOF corresponds to which rank.
+     * and collects information about which DOF corresponds to which rank. Can only
+     * be used, if \ref partitionVec is set correctly. This is only the case, when
+     * the macro mesh is partitioned.
      *
      * @param[out] partionDOFs   Stores to each DOF pointer the set of ranks the DOF is
      *                           part of.
@@ -215,10 +221,20 @@ namespace AMDiS {
     int nRankDOFs;
 
     /** \brief 
-     * Defines the interioir boundaries of the domain that result from partitioning
-     * the whole mesh.
+     * Defines the interior boundaries of the domain that result from partitioning
+     * the whole mesh. Contains only the boundaries, which are owned by the rank, i.e.,
+     * the object gives for every neighbour rank i the boundaries this rank owns and 
+     * shares with rank i.
+     */
+    InteriorBoundary myIntBoundary;
+    
+    /** \brief
+     * Defines the interior boundaries of the domain that result from partitioning
+     * the whole mesh. Contains only the boundaries, which are not owned by the rank,
+     * i.e., the object gives for every neighbour rank i the boundaries that are
+     * owned by rank i and are shared with this rank.
      */
-    InteriorBoundary interiorBoundary;
+    InteriorBoundary otherIntBoundary;
 
     /** \brief
      * This map contains for each rank the list of dofs the current rank must send

AMDiS/src/UmfPackSolver.h

@@ -51,7 +51,8 @@ namespace AMDiS {
       /** \brief
        * Returns a new UmfPackSolver object.
        */
-      OEMSolver* create() { 
+      OEMSolver* create() 
+      { 
 	return new UmfPackSolver(this->name); 
       }
     };
@@ -67,7 +68,11 @@ namespace AMDiS {
     }
     
     /// Destructor
-    ~UmfPackSolver() { if (solver) delete solver;}
+    ~UmfPackSolver() 
+    { 
+      if (solver) 
+	delete solver;
+    }
 
     /// Solves the system directly
     int solveSystem(const DOFMatrix::base_matrix_type&      A,

AMDiS/src/Utilities.h

-// ============================================================================
-// ==                                                                        ==
-// == AMDiS - Adaptive multidimensional simulations                          ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  crystal growth group                                                  ==
-// ==                                                                        ==
-// ==  Stiftung caesar                                                       ==
-// ==  Ludwig-Erhard-Allee 2                                                 ==
-// ==  53175 Bonn                                                            ==
-// ==  germany                                                               ==
-// ==                                                                        ==
-// ============================================================================
-// ==                                                                        ==
-// ==  http://www.caesar.de/cg/AMDiS                                         ==
-// ==                                                                        ==
-// ============================================================================
-
-/** \file Utilities.h */
-
-#ifndef AMDIS_UTILITIES_H
-#define AMDIS_UTILITIES_H
-
-#include <vector>
-#include <algorithm>
-
-
-
-template <typename T>
-std::vector<T> inline invert_reorder(const std::vector<T>& v)
-{
-    T my_max= *std::max_element(v.begin(), v.end()) + 1;
-    std::vector<T> u(my_max, T(-1));
-    for (int i= 0; i < v.size(); i++) 
-	if (v[i] >= 0) {
-	    assert(u[v[i]] == -1); // check double insertion
-	    u[v[i]]= T(i);
-    }
-    assert(find(u.begin(), u.end(), T(-1)) == u.end()); // check if all entries are set
-    return u;
-}
-
-#endif  // AMDIS_UTILITIES_H

AMDiS/src/ValueWriter.h

@@ -29,18 +29,10 @@
 #include "Flag.h"
 #include "Mesh.h"
 #include "DataCollector.h"
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
-  class DOFAdmin;
-  class ElInfo;
-
-  template<typename T> class DOFVector;
-
-  // ============================================================================
-  // ===== class ValueWriter ====================================================
-  // ============================================================================
-
   /** \ingroup Output
    * \brief
    * ValueWriter is a static class which writes the values of a DOFVector
@@ -53,9 +45,7 @@ namespace AMDiS {
   class ValueWriter
   {
   public:
-    /** \brief
-     * Writes DOFVector values to values->feSpace->mesh.
-     */
+    /// Writes DOFVector values to values->feSpace->mesh.
     static void writeValues(DataCollector *dc,
 			    const std::string filename,
 			    double time = 0.0,
@@ -64,34 +54,22 @@ namespace AMDiS {
 			    bool (*writeElem)(ElInfo*) = NULL);
 			   
   protected:
-    /** \brief
-     * File to which the values should be written
-     */
+    /// File to which the values should be written
     static FILE *valueFile;
 
-    /** \brief
-     * Global interpolation point indexing
-     */
+    /// Global interpolation point indexing
     static DOFVector<int> *interpPointInd;
 
-    /** \brief
-     * list of coords for each dof
-     */
+    /// list of coords for each dof
     static DOFVector< std::list<WorldVector<double> > > *dofCoords;
 
-    /** \brief
-     * DOFAdmin of values
-     */
+    /// DOFAdmin of values
     static const DOFAdmin *admin;
 
-    /** \brief
-     * Pointer to have a global access to values
-     */
+    /// Pointer to have a global access to values
     static DOFVector<double> *valueVec;
 
-    /** \brief
-     * Total number of interpolation points.
-     */
+    /// Total number of interpolation points.
     static int ni;
   };
 

AMDiS/src/VertexInfo.h

@@ -26,35 +26,27 @@
 
 namespace AMDiS {
   
-  /** \brief
-   * Stores coordinates and output index for one vertex.
-   */
+  /// Stores coordinates and output index for one vertex.
   class VertexInfo 
   {
   public:
-    /** \brief
-     * Coordinates for this vertex.
-     */
+    /// Coordinates for this vertex.
     WorldVector<double> coords;
     
-    /** \brief
-     * Index for the output file.
-     */
+    /// Index for the output file.
     int outputIndex;
     
-    /** \brief
-     * Used to check, whether coords are already stored for a given dof.
-     */ 
-    bool operator==(const WorldVector<double>& c) {
+    /// Used to check, whether coords are already stored for a given dof.
+    bool operator==(const WorldVector<double>& c) 
+    {
       return (c == coords);
-    };
+    }
     
-    /** \brief
-     * Used to check, whether coords are already stored for a given dof.
-     */ 
-    bool operator!=(const WorldVector<double>& c) {
+    /// Used to check, whether coords are already stored for a given dof.
+    bool operator!=(const WorldVector<double>& c) 
+    {
       return (c != coords);
-    };
+    }
   };
 }
 

AMDiS/src/VertexVector.h

@@ -21,24 +21,27 @@ namespace AMDiS {
 
     ~VertexVector();
 
-    const DOFAdmin *getAdmin() { return admin; };
+    const DOFAdmin *getAdmin() 
+    { 
+      return admin; 
+    }
 
-    void resize(int size) {
-      int i, oldSize = static_cast<int>(vec.size());
+    void resize(int size) 
+    {
+      int oldSize = static_cast<int>(vec.size());
       vec.resize(size);
-      for(i = oldSize; i < size; i++) {
+      for (int i = oldSize; i < size; i++)
 	vec[i] = i;
-      }
     }
 
     void set(DegreeOfFreedom val);
 
-    void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF) {
+    void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF) 
+    {
       DOFContainer::compressDOFContainer(size, newDOF);
       int totalSize = getAdmin()->getSize();
-      for(int i = size; i < totalSize; i++) {
+      for(int i = size; i < totalSize; i++)
 	vec[i] = i;
-      }
     }
 
   protected: