Work on pdd.

AMDiS/src/DOFVector.h

@@ -40,10 +40,6 @@
 #include "DOFMatrix.h" 
 #include "BasisFunction.h"
 
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-#include "petscao.h"
-#endif
- 
 namespace AMDiS {
 
   template<typename T> 
@@ -547,14 +543,6 @@ namespace AMDiS {
 
     DOFVector<WorldVector<T> >* getRecoveryGradient(DOFVector<WorldVector<T> >*) const;
 
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    /// Sets the petsc application ordering object to map dof indices.
-    void useApplicationOrdering(AO *ao) 
-    {
-      applicationOrdering = ao;
-    }
-#endif
-
   protected: 
     /// Used by Int while mesh traversal
     static int Int_fct(ElInfo* elinfo);
@@ -590,11 +578,6 @@ namespace AMDiS {
     /// Used for mesh traversal
     static DOFVector<T> *traverseVector;
 
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    /// Petsc application ordering to map dof indices.
-    AO *applicationOrdering;
-#endif
-
   protected:
     /// Used while calculating vector norms
     static FastQuadrature *quad_fast;

AMDiS/src/DOFVector.hh

@@ -110,10 +110,6 @@ namespace AMDiS {
       coarsenOperation(NO_OPERATION)
   {
     init(f, n);
-
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    applicationOrdering = NULL;
-#endif
   } 
 
   template<typename T>

AMDiS/src/DataCollector.cc

@@ -215,11 +215,10 @@ namespace AMDiS {
     // read element region
     ElementData *ed = elInfo->getElement()->getElementData(ELEMENT_REGION);
     
-    if (ed) {
+    if (ed)
       elementInfo.elementRegion = dynamic_cast<ElementRegion_ED*>(ed)->getRegion();
-    } else {
-      elementInfo.elementRegion = -1;
-    }
+    else 
+      elementInfo.elementRegion = -1;    
    
     // read surface regions to element info
     ed = elInfo->getElement()->getElementData(SURFACE_REGION);

AMDiS/src/FixVec.h

@@ -25,23 +25,18 @@
 #include <iostream>
 #include "Global.h"
 #include "MatrixVector.h"
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
-  class Mesh;
-  template<typename T> class DimMat;
-  template<typename T> class WorldVector;
-  template<typename T> class WorldMatrix;
-
-  /** determines how to initialize a FixVec when constructed */
-  enum InitType 
-    {
-      NO_INIT = 0,       /**< no initialisation */
-      VALUE_LIST = 1,    /**< a complete value list is given */
-      DEFAULT_VALUE = 2, /**< all values ar set to a given default value */
-      UNIT_VECTOR = 3,   /**< the i-th value is 1, all other values are 0 */
-      UNIT_MATRIX = 4    /**< values at the diagonal of a matrix are set to one */
-    };
+  /// determines how to initialize a FixVec when constructed.
+  enum InitType {
+    NO_INIT = 0,       /**< no initialisation */
+    VALUE_LIST = 1,    /**< a complete value list is given */
+    DEFAULT_VALUE = 2, /**< all values ar set to a given default value */
+    UNIT_VECTOR = 3,   /**< the i-th value is 1, all other values are 0 */
+    UNIT_MATRIX = 4    /**< values at the diagonal of a matrix are set to one */
+  };
 
   /** \ingroup Common
    * \brief
@@ -121,10 +116,6 @@ namespace AMDiS {
   };
 
 
-  // ===========================================================================
-  // ===== class VectorOfFixVecs ===============================================
-  // ===========================================================================
-
   /** \ingroup Common
    * \brief
    * Contains an vector of FixVecs of same type. To simply allocate an array of 
@@ -267,9 +258,6 @@ namespace AMDiS {
     std::vector<FixVecType*> vec;
   };
 
-  // ===========================================================================
-  // ===== class MatrixOfFixVecs ===============================================
-  // ===========================================================================
 
   /** \ingroup Common
    * \brief
@@ -401,9 +389,6 @@ namespace AMDiS {
     }
   };
 
-  // ===========================================================================
-  // ===== class WorldVector ===================================================
-  // ===========================================================================
 
   /** \ingroup Common
    * \brief
@@ -450,10 +435,6 @@ namespace AMDiS {
   };
 
 
-  // ===========================================================================
-  // ===== class WorldMatrix ===================================================
-  // ===========================================================================
-
   /** \ingroup Common
    * \brief
    * A WorldMatrix is a FixVec with DIM_OF_WORLD WorldVector entries.
@@ -513,10 +494,6 @@ namespace AMDiS {
   };
 
 
-  // ===========================================================================
-  // ===== global functions ====================================================
-  // ===========================================================================
-
 
   /// returns the euclidian distance of a and b
   template<typename T, GeoIndex d>
@@ -526,9 +503,8 @@ namespace AMDiS {
   template<typename T, GeoIndex d>
   std::ostream& operator <<(std::ostream& out, const FixVec<T,d>& fixvec)
   {
-    for (int i = 0; i < fixvec.getSize()-1; i++) {
+    for (int i = 0; i < fixvec.getSize() - 1; i++)
       out << fixvec[i] << " ";
-    }
     out << fixvec[fixvec.getSize()-1];
     return out;
   }
@@ -585,9 +561,10 @@ namespace AMDiS {
 			 const WorldVector<double>& v2) 
   {
     int dow = Global::getGeo(WORLD);
-    for (int i = 0; i < dow; i++) {
-      if (abs(v1[i] - v2[i]) > DBL_TOL) return false;
-    }
+    for (int i = 0; i < dow; i++)
+      if (abs(v1[i] - v2[i]) > DBL_TOL) 
+	return false;
+
     return true;
   }
 

AMDiS/src/ParallelDomainProblem.cc

@@ -90,11 +90,17 @@ namespace AMDiS {
 
     int globalRefinement = 0;
     GET_PARAMETER(0, "testgr", "%d", &globalRefinement);
-    std::cout << "GR = " << globalRefinement << std::endl;
 
-    refinementManager->globalRefine(mesh, globalRefinement);
+    if (globalRefinement > 0) {
+      refinementManager->globalRefine(mesh, globalRefinement);
 
-    updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
+      updateLocalGlobalNumbering(nRankDOFs, nOverallDOFs);
+    }
+
+
+#if (DEBUG != 0)
+    testInteriorBoundary();
+#endif
 
     // === Create petsc matrix. ===
 
@@ -257,7 +263,7 @@ namespace AMDiS {
 	if (partitionData->getLevel() == 0) {
 	  elNum = element->getIndex();
 	}
-	TEST_EXIT(elNum != -1)("invalid element number\n");
+	TEST_EXIT_DBG(elNum != -1)("invalid element number\n");
 	if (element->isLeaf()) {
 	  elemWeights[elNum] += 1.0;
 	  localWeightSum += 1.0;
@@ -418,7 +424,7 @@ namespace AMDiS {
 	      break;
 
 	  // The element must always be found, because the list is just in another order.
-	  TEST_EXIT(k < static_cast<int>(rankIt->second.size()))
+	  TEST_EXIT_DBG(k < static_cast<int>(rankIt->second.size()))
 	    ("Should never happen!\n");
 
 	  // Swap the current with the found element.
@@ -627,7 +633,7 @@ namespace AMDiS {
 	  }
 	}
 
-	TEST_EXIT(found)("Should not happen!\n");
+	TEST_EXIT_DBG(found)("Should not happen!\n");
       }
 
       delete [] recvBuffers[i];
@@ -690,9 +696,9 @@ namespace AMDiS {
 	  ERROR_EXIT("Should never happen!\n");
 	}
 
-	TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
+	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
 	  ("Should never happen!\n");
-	TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
+	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
 	  ("Should never happen!\n");
 
 	newBoundaryDOFs[dof1] = boundaryDOFs[dof1];
@@ -736,9 +742,9 @@ namespace AMDiS {
 	  ERROR_EXIT("Should never happen!\n");
 	}
 
-	TEST_EXIT(boundaryDOFs.find(dof1) != boundaryDOFs.end())
+	TEST_EXIT_DBG(boundaryDOFs.find(dof1) != boundaryDOFs.end())
 	  ("Should never happen!\n");
-	TEST_EXIT(boundaryDOFs.find(dof2) != boundaryDOFs.end())
+	TEST_EXIT_DBG(boundaryDOFs.find(dof2) != boundaryDOFs.end())
 	  ("Should never happen!\n");
 
 	rankDOFs.erase(dof1);
@@ -749,7 +755,8 @@ namespace AMDiS {
 	std::vector<const DegreeOfFreedom*> boundDOFs;
 	addAllDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary,
 		   boundDOFs);
-	for (int i = 0; i < static_cast<int>(boundDOFs.size()); i++) {
+	for (int i = static_cast<int>(boundDOFs.size()) - 1; i >= 0; i--) {
+	//	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]);
@@ -847,9 +854,6 @@ namespace AMDiS {
 	   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;
@@ -964,6 +968,22 @@ namespace AMDiS {
   }
 
 
+  void ParallelDomainProblemBase::testInteriorBoundary()
+  {
+    // Maps to each neighbour rank an array of WorldVectors. This array contain the 
+    // coordinates of all dofs this rank shares on the interior boundary with the 
+    // neighbour rank.
+    std::map<int, std::vector<WorldVector<double> > > sendCoords;
+    std::map<int, int> nRecvCoords;
+    
+    for (std::map<const DegreeOfFreedom*, int>::iterator it = boundaryDOFs.begin();
+	 it != boundaryDOFs.end();
+	 ++it) {
+      
+    }
+  }
+
+
   ParallelDomainProblemScal::ParallelDomainProblemScal(const std::string& name,
 						       ProblemScal *problem,
 						       ProblemInstatScal *problemInstat)
@@ -989,23 +1009,22 @@ namespace AMDiS {
 
   Flag ParallelDomainProblemScal::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
   {
-    //    return iterationIF->oneIteration(adaptInfo, toDo);
+    FUNCNAME("ParallelDomainProblemScal::oneIteration()");
 
-    Flag flag = dynamic_cast<StandardProblemIteration*>(iterationIF)->buildAndAdapt(adaptInfo, toDo);
+    Flag flag =  dynamic_cast<StandardProblemIteration*>(iterationIF)->
+      buildAndAdapt(adaptInfo, toDo);
 
-    fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());
-
-    solvePetscMatrix(probScal->getSolution());
-
-//     if (toDo.isSet(SOLVE))
-//       iterationIF->getProblem()->solve(adaptInfo, false);
-
-//     if (toDo.isSet(SOLVE_RHS))
-//       iterationIF->getProblem()->solve(adaptInfo, true);
+    if (toDo.isSet(SOLVE)) {      
+      fillPetscMatrix(probScal->getSystemMatrix(), probScal->getRHS());      
+      solvePetscMatrix(probScal->getSolution());
+    }
 
-//     if (toDo.isSet(ESTIMATE)) 
-//       iterationIF->getProblem()->estimate(adaptInfo);
+    if (toDo.isSet(SOLVE_RHS)) {
+      ERROR_EXIT("Not yet implemented!\n");
+    }
 
+    if (toDo.isSet(ESTIMATE))
+      iterationIF->getProblem()->estimate(adaptInfo);
 
     return flag;
   }

AMDiS/src/ParallelDomainProblem.h

@@ -167,7 +167,15 @@ namespace AMDiS {
     void createDOFMemberInfo(std::map<const DegreeOfFreedom*, std::set<int> >& partitionDOFs,
 			     std::vector<const DegreeOfFreedom*>& rankDOFs,
 			     std::map<const DegreeOfFreedom*, int>& boundaryDOFs);
+
      
+    /** \brief
+     * This function is used for debugging only. It traverses all interior boundaries
+     * and compares the dof indices on them with the dof indices of the boundarys
+     * neighbours. The function fails, when dof indices on an interior boundary does
+     * not fit together.
+     */
+    void testInteriorBoundary();
 
   protected:
     ///

AMDiS/src/ProblemScal.cc

@@ -141,15 +141,6 @@ namespace AMDiS {
       rhs->getBoundaryManager()->addBoundaryCondition(periodic);
   }
 
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-  void ProblemScal::useApplicationOrdering(AO *ao) 
-  {
-    applicationOrdering = ao;
-    rhs->useApplicationOrdering(ao);
-  }
-#endif
-
-
   void ProblemScal::createMesh()
   {
     TEST_EXIT(Parameters::initialized())("parameters not initialized\n");

AMDiS/src/ProblemScal.h

@@ -31,10 +31,6 @@
 #include "Boundary.h"
 #include "StandardProblemIteration.h"
 
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-#include "petscao.h"
-#endif
-
 namespace AMDiS {
 
   class ProblemScal : public ProblemStatBase,
@@ -56,9 +52,6 @@ namespace AMDiS {
 	useGetBound(true),
 	refinementManager(NULL),
 	coarseningManager(NULL),
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-        applicationOrdering(NULL),
-#endif
 	info(10)
     {}
 
@@ -204,11 +197,6 @@ namespace AMDiS {
     /// Adds periodic boundary conditions.
     virtual void addPeriodicBC(BoundaryType type);
     
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    /// Sets the petsc application ordering object to map dof indices.
-    void useApplicationOrdering(AO *ao);
-#endif
-
     /** \name getting methods
      * \{ 
      */
@@ -381,11 +369,6 @@ namespace AMDiS {
      */
     CoarseningManager *coarseningManager;
   
-#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    /// Petsc application ordering to map dof indices.
-    AO *applicationOrdering;
-#endif
-
     /// Info level.
     int info;
 

AMDiS/src/StandardProblemIteration.cc

@@ -12,9 +12,9 @@ namespace AMDiS {
   {
     FUNCNAME("StandardProblemIteration::beginIteration()");
 
-    INFO(info,4)("\n");
-    INFO(info,4)("begin of iteration number: %d\n", adaptInfo->getSpaceIteration() + 1);
-    INFO(info,4)("=============================\n");
+    INFO(info, 4)("\n");
+    INFO(info, 4)("begin of iteration number: %d\n", adaptInfo->getSpaceIteration() + 1);
+    INFO(info, 4)("=============================\n");
   }
 
   Flag StandardProblemIteration::oneIteration(AdaptInfo *adaptInfo, Flag toDo) 
@@ -39,9 +39,9 @@ namespace AMDiS {
   {
     FUNCNAME("StandardProblemIteration::endIteration()");
 
-    INFO(info,4)("\n");
-    INFO(info,4)("end of iteration number: %d\n", adaptInfo->getSpaceIteration() + 1);
-    INFO(info,4)("=============================\n");
+    INFO(info, 4)("\n");
+    INFO(info, 4)("end of iteration number: %d\n", adaptInfo->getSpaceIteration() + 1);
+    INFO(info, 4)("=============================\n");
   }
 
   Flag StandardProblemIteration::buildAndAdapt(AdaptInfo *adaptInfo, Flag toDo)
@@ -50,11 +50,10 @@ namespace AMDiS {
 
     Flag flag = 0, markFlag = 0;
 
-    if (toDo.isSet(MARK)) {
+    if (toDo.isSet(MARK))
       markFlag = problem->markElements(adaptInfo);
-    }  else {
+    else
       markFlag = 3;
-    }
 
     if (toDo.isSet(BUILD))
       problem->buildBeforeRefine(adaptInfo, markFlag);