Work on adaptivity in parallelization code and some code refactoring in this part too.

AMDiS/libtool

@@ -82,13 +82,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.3.2-gcc//bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"
 
 # A language-specific compiler.
-CC="gcc"
+CC="/usr/lib/openmpi/1.3.2-gcc//bin/mpicc"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -174,7 +174,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6801,13 +6801,13 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.3.2-gcc//bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"
 
 # A language-specific compiler.
-CC="g++"
+CC="/usr/lib/openmpi/1.3.2-gcc//bin/mpiCC"
 
 # Is the compiler the GNU C compiler?
 with_gcc=yes
@@ -6893,7 +6893,7 @@ dlopen_self=unknown
 dlopen_self_static=unknown
 
 # Compiler flag to prevent dynamic linking.
-link_static_flag="-static"
+link_static_flag=""
 
 # Compiler flag to turn off builtin functions.
 no_builtin_flag=" -fno-builtin"
@@ -6960,11 +6960,11 @@ predeps=""
 
 # Dependencies to place after the objects being linked to create a
 # shared library.
-postdeps="-lstdc++ -lm -lgcc_s -lc -lgcc_s"
+postdeps="-lmpi_cxx -lmpi -lopen-rte -lopen-pal -ldl -lnsl -lutil -ldl -lstdc++ -lm -lgcc_s -lpthread -lc -lgcc_s"
 
 # The library search path used internally by the compiler when linking
 # a shared library.
-compiler_lib_search_path=`echo "-L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
+compiler_lib_search_path=`echo "-L/usr/lib/openmpi/1.3.2-gcc/lib -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2 -L/usr/lib/gcc/i386-redhat-linux/4.1.2/../../.." | $SED -e "s@${gcc_dir}@\${gcc_dir}@g;s@${gcc_ver}@\${gcc_ver}@g"`
 
 # Method to check whether dependent libraries are shared objects.
 deplibs_check_method="pass_all"
@@ -7109,7 +7109,7 @@ AR="ar"
 AR_FLAGS="cru"
 
 # A C compiler.
-LTCC="gcc"
+LTCC="/usr/lib/openmpi/1.3.2-gcc//bin/mpicc"
 
 # LTCC compiler flags.
 LTCFLAGS="-g -O2"

AMDiS/src/AMDiS_fwd.h

@@ -58,6 +58,7 @@ namespace AMDiS {
   class MacroInfo;
   class Marker;
   class Mesh; 
+  class MeshStructure;
   class OEMSolver;
   class Operator;
   class OperatorTerm;

AMDiS/src/DOFMatrix.cc

@@ -195,7 +195,7 @@ namespace AMDiS {
       if (condition && condition->isDirichlet()) {
 	if (condition->applyBoundaryCondition()) {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
- 	  if (rankDofs[rowIndices[i]]) 
+ 	  if ((*rankDofs)[rowIndices[i]]) 
  	    applyDBCs.insert(static_cast<int>(row));
 #else
 	  applyDBCs.insert(static_cast<int>(row));

AMDiS/src/DOFMatrix.h

@@ -406,7 +406,7 @@ namespace AMDiS {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
     {
-      rankDofs = dofmap;
+      rankDofs = &dofmap;
     }
 #endif
 
@@ -488,7 +488,7 @@ namespace AMDiS {
      int nnzPerRow;
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    std::map<DegreeOfFreedom, bool> rankDofs;
+    std::map<DegreeOfFreedom, bool> *rankDofs;
 #endif
 
     /// Inserter object: implemented as pointer, allocated and deallocated as needed

AMDiS/src/DOFVector.h

@@ -211,16 +211,15 @@ namespace AMDiS {
     }
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-    inline void setRankDofs(std::map<DegreeOfFreedom, bool> dofmap)
+    inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap)
     {
-      rankDofs.clear();
-
-      //      rankDofs = dofmap;
+      //      rankDofs = &dofmap;
     }
 
     inline bool isRankDof(DegreeOfFreedom dof)
     {
-      return rankDofs[dof];
+      TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n");
+      return (*rankDofs)[dof];
     }
 #endif
 
@@ -269,7 +268,7 @@ namespace AMDiS {
 
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
   public:
-    std::map<DegreeOfFreedom, bool> rankDofs;
+    std::map<DegreeOfFreedom, bool> *rankDofs;
 #endif
   };
 

AMDiS/src/DirichletBC.cc

@@ -54,7 +54,7 @@ namespace AMDiS {
 
     for (int i = 0; i < nBasFcts; i++) {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-      if (vector->isRankDof(dofIndices[i]))
+      //      if (vector->isRankDof(dofIndices[i]))
 #endif
       if (localBound[i] == boundaryType) {
 	if (f) {

AMDiS/src/ElInfo.h

@@ -344,10 +344,7 @@ namespace AMDiS {
      */
     double calcDet(const FixVec<WorldVector<double>, VERTEX> &coords) const;
 
-    /** \brief
-     * Checks whether flag is set in ElInfo's \ref fillFlag. If not, the program
-     * exits.
-     */
+    /// Checks whether flag is set in ElInfo's \ref fillFlag. If not, the program exits.
     void testFlag(const Flag& flag) const;
 
     /** \brief
@@ -411,7 +408,7 @@ namespace AMDiS {
 
       ERROR("virtual function not implemented in this sub-class ");
     
-      return(0.0);
+      return 0.0;
     }
 
     /// Get ElInfo's \ref elType.

AMDiS/src/Element.cc

@@ -5,6 +5,7 @@
 #include "FixVec.h"
 #include "ElementRegion_ED.h"
 #include "Serializer.h"
+#include "MeshStructure.h"
 
 namespace AMDiS {
 
@@ -565,4 +566,31 @@ namespace AMDiS {
     return result;
   }
 
+  void fitElementToMeshCode(RefinementManager *refineManager, MeshStructure &code, 
+			    Element *el, int ithSide, int elType)
+  {
+    if (code.isLeafElement())
+      return;
+
+    if (el->isLeaf()) {
+      el->setMark(1);
+      refineManager->refineMesh(el->getMesh());
+    }
+
+    int s1 = el->getSideOfChild(0, ithSide, elType);
+    int s2 = el->getSideOfChild(1, ithSide, elType);
+
+    code.nextElement();
+    
+    if (s1 != -1) 
+      fitElementToMeshCode(refineManager, code, 
+			   el->getFirstChild(), s1, el->getChildType(elType));
+
+    code.nextElement();
+    
+    if (s2 != -1)
+      fitElementToMeshCode(refineManager, code, 
+			   el->getSecondChild(), s2, el->getChildType(elType));
+  }
+
 }

AMDiS/src/Element.h

@@ -574,6 +574,11 @@ namespace AMDiS {
     friend class Mesh;
   };
 
+  void fitElementToMeshCode(RefinementManager *refineManager, 
+			    MeshStructure &code, 
+			    Element *el, 
+			    int ithSide, 
+			    int elType);
 }
 
 #endif  // AMDIS_ELEMENT_H

AMDiS/src/FiniteElemSpace.h

@@ -46,7 +46,7 @@ namespace AMDiS {
     static FiniteElemSpace *provideFESpace(DOFAdmin *admin,
 					   const BasisFunction *basFcts,
 					   Mesh *mesh,
-					   std::string name_ = "");
+					   std::string name = "");
 
     /// Destructor.
     ~FiniteElemSpace();
@@ -92,7 +92,7 @@ namespace AMDiS {
      * Constructs a FiniteElemSpace with name name_ and the given DOFAdmin,
      * BasisFunction and Mesh.
      */
-    FiniteElemSpace(DOFAdmin* admin_,
+    FiniteElemSpace(DOFAdmin* admin,
 		    const BasisFunction* basisFcts, 
 		    Mesh* mesh,
 		    std::string name = "");

AMDiS/src/Mesh.h

@@ -99,7 +99,7 @@ namespace AMDiS {
     }
 
     /// Returns \ref dim of the mesh
-    inline int getDim() const \
+    inline int getDim() const
     { 
       return dim; 
     }

AMDiS/src/MeshStructure.h

@@ -53,6 +53,10 @@ namespace AMDiS {
       reset();
     }
 
+    /** \brief
+     * Sets all position counters, that are used to traverse the code, to the starting
+     * position. The code itself is not changed.
+     */
     void reset();
 
     inline void commit() 

AMDiS/src/ParallelDomainBase.cc

@@ -40,13 +40,13 @@ namespace AMDiS {
   ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF,
 					 ProblemTimeInterface *tIF,
 					 FiniteElemSpace *fe,
-					 RefinementManager *refineManager)
+					 RefinementManager *refinementManager)
     : iterationIF(iIF),
       timeIF(tIF),
       name(iIF->getName()),
       feSpace(fe),
       mesh(fe->getMesh()),
-      refinementManager(refineManager),
+      refineManager(refinementManager),
       initialPartitionMesh(true),
       nRankDofs(0),
       rstart(0),
@@ -100,14 +100,7 @@ namespace AMDiS {
 
     // === Create new global and local DOF numbering. ===
 
-    // Set of all DOFs of the rank.
-    std::vector<const DegreeOfFreedom*> rankDofs;
-    // Number of DOFs in ranks partition that are owned by the rank.
-    nRankDofs = 0;
-    // Number of all DOFs in the macro mesh.
-    int nOverallDOFs = 0;
-
-    createLocalGlobalNumbering(rankDofs, nRankDofs, nOverallDOFs);
+    createLocalGlobalNumbering();
 
     // === Create interior boundary information ===
 
@@ -120,6 +113,7 @@ namespace AMDiS {
 #if (DEBUG != 0)
     dbgTestElementMap(elMap);
     dbgTestInteriorBoundary();
+    dbgTestCommonDofs(true);
 #endif
 
     // === Reset all DOFAdmins of the mesh. ===
@@ -136,33 +130,14 @@ namespace AMDiS {
     GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement);
 
     if (globalRefinement > 0) {
-      refinementManager->globalRefine(mesh, globalRefinement);
-
-#if (DEBUG != 0)
-      elMap.clear();
-      dbgCreateElementMap(elMap);
-#endif
-
-      updateLocalGlobalNumbering(nRankDofs, nOverallDOFs);
+      refineManager->globalRefine(mesh, globalRefinement);
 
-      updateDofAdmins();
-
-#if (DEBUG != 0)
-      dbgTestElementMap(elMap);
-#endif
+      updateLocalGlobalNumbering();
 
       // === Update periodic mapping, if there are periodic boundaries. ===
+
       createPeriodicMap();
     }
-
-    lastMeshChangeIndex = mesh->getChangeIndex();
-
-#if (DEBUG != 0)
-    dbgTestCommonDofs(true);
-#endif
-
-    nRankRows = nRankDofs * nComponents;
-    nOverallRows = nOverallDOFs * nComponents;
   }
 
 
@@ -759,6 +734,8 @@ namespace AMDiS {
 
   void ParallelDomainBase::checkMeshChange()
   {
+    FUNCNAME("ParallelDomainBase::checkMeshChange()");
+
     // === If mesh has not been changed, return. ===
 
     if (mesh->getChangeIndex() == lastMeshChangeIndex)
@@ -811,8 +788,14 @@ namespace AMDiS {
 	elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, 
 		    boundIt->rankObj.elType);
 
-	if (elCode.getCode() != recvCodes[i].getCode())
+	if (elCode.getCode() != recvCodes[i].getCode()) {
+	  TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n");
+	  //	  recvCodes[i].reset();
+// 	  fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el,
+// 			       boundIt->rankObj.ithObj, boundIt->rankObj.elType);
+
 	  meshFitTogether = false;
+	}
 
 	i++;
       }
@@ -823,7 +806,7 @@ namespace AMDiS {
       exit(0);    
     }
 
-    lastMeshChangeIndex = mesh->getChangeIndex();
+    updateLocalGlobalNumbering();
   }
 
   
@@ -1207,23 +1190,21 @@ namespace AMDiS {
   }
 
 
-  void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDofs,
-						      int& nRankDofs, 
-						      int& nOverallDOFs)
+  void ParallelDomainBase::createLocalGlobalNumbering()
   {
     FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()");
 
     // === Get rank information about DOFs. ===
 
     // Stores to each DOF pointer the set of ranks the DOF is part of.
-    std::map<const DegreeOfFreedom*, std::set<int> > partitionDOFs;
-    DofContainer rankAllDofs;
+    DofToPartitions partitionDofs;
+    DofContainer rankDofs, rankAllDofs;
     DofToRank boundaryDofs;
 
-    createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
+    createDofMemberInfo(partitionDofs, rankDofs, rankAllDofs, boundaryDofs, vertexDof);
 
     nRankDofs = rankDofs.size();
-    nOverallDOFs = partitionDOFs.size();
+    int nOverallDofs = partitionDofs.size();
 
 
     // === Get starting position for global rank dof ordering. ====
@@ -1284,8 +1265,8 @@ namespace AMDiS {
 	// If the boundary dof is a rank dof, it must be send to other ranks.
 
 	// Search for all ranks that have this dof too.
-	for (std::set<int>::iterator itRanks = partitionDOFs[it->first].begin();
-	     itRanks != partitionDOFs[it->first].end();
+	for (std::set<int>::iterator itRanks = partitionDofs[it->first].begin();
+	     itRanks != partitionDofs[it->first].end();
 	     ++itRanks) {
 	  if (*itRanks != mpiRank) {
 	    TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1)
@@ -1376,12 +1357,22 @@ namespace AMDiS {
 
     createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
 			rankDofsNewGlobalIndex);
+
+    nRankRows = nRankDofs * nComponents;
+    nOverallRows = nOverallDofs * nComponents;
+
+    lastMeshChangeIndex = mesh->getChangeIndex();
   }
 
-  void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDofs, int& nOverallDOFs)
+  void ParallelDomainBase::updateLocalGlobalNumbering()
   {
     FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()");
 
+#if (DEBUG != 0)
+    ElementIdxToDofs elMap;
+    dbgCreateElementMap(elMap);
+#endif
+
     typedef std::set<const DegreeOfFreedom*> DofSet;
 
     // === Get all DOFs in ranks partition. ===
@@ -1525,10 +1516,11 @@ namespace AMDiS {
 
     // === Calculate number of overall DOFs of all partitions. ===
 
-    mpiComm.Allreduce(&nRankDofs, &nOverallDOFs, 1, MPI_INT, MPI_SUM);
+    int nOverallDofs = 0;
+    mpiComm.Allreduce(&nRankDofs, &nOverallDofs, 1, MPI_INT, MPI_SUM);
 
 
-    // Do not change the indices now, but create a new indexing a store it here.
+    // Do not change the indices now, but create a new indexing and store it here.
     DofIndexMap rankDofsNewLocalIndex;
     isRankDof.clear();
     int i = 0;
@@ -1611,6 +1603,20 @@ namespace AMDiS {
 
     createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex,
 			rankDofsNewGlobalIndex);
+
+    nRankRows = nRankDofs * nComponents;
+    nOverallRows = nOverallDofs * nComponents;
+
+    // === Update dof admins due to new number of dofs. ===
+  
+    updateDofAdmins();
+
+    lastMeshChangeIndex = mesh->getChangeIndex();
+
+#if (DEBUG != 0)
+    dbgTestElementMap(elMap);
+    dbgTestCommonDofs(true);
+#endif
   }
 
   void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex,
@@ -1722,7 +1728,6 @@ namespace AMDiS {
     sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue);
   }
 
-
   void ParallelDomainBase::createPeriodicMap()
   {
     FUNCNAME("ParallelDomainBase::createPeriodicMap()");

AMDiS/src/ParallelDomainBase.h

@@ -201,20 +201,11 @@ namespace AMDiS {
     /// Removes all macro elements from the mesh that are not part of ranks partition.
     void removeMacroElements();
 
+    /// Creates from a macro mesh a correct local and global DOF index numbering.
+    void createLocalGlobalNumbering();
 
-    /** \brief
-     * Creates from a macro mesh a correct local and global DOF index numbering.
-     *
-     * \param[out] rankDOFs      Returns all DOFs from the macro mesh, which are owned
-     *                           by the rank after partitioning the macro mesh.
-     * \param[out] nRankDOFs     Number of DOFs owned by rank.
-     * \param[out] nOverallDOFs  Number of all DOFs in macro mesh.
-     */
-    void createLocalGlobalNumbering(DofContainer& rankDOFs,
-				    int& nRankDOFs, 
-				    int& nOverallDOFs);
-
-    void updateLocalGlobalNumbering(int& nRankDOFs, int& nOverallDOFs);
+    /// Updates the local and global DOF numbering after the mesh has been changed.
+    void updateLocalGlobalNumbering();
 
     /** \brief
      * Creates to all dofs in rank's partition that are on a periodic boundary the
@@ -452,12 +443,16 @@ namespace AMDiS {
     /// Mesh of the problem.
     Mesh *mesh;
 
+    /** \brief
+     * A refinement manager that should be used on the mesh. It is used to refine
+     * elements at interior boundaries in order to fit together with elements on the
+     * other side of the interior boundary.
+     */    
+    RefinementManager *refineManager;
+
     /// Info level.
     int info;
 
-    /// Refinement manager for the mesh.
-    RefinementManager *refinementManager;
-
     /// Pointer to the paritioner which is used to devide a mesh into partitions.
     ParMetisPartitioner *partitioner;
 

AMDiS/src/ProblemVec.cc

@@ -621,7 +621,6 @@ namespace AMDiS {
     FUNCNAME("ProblemVec::buildAfterCoarsen()");
 
     //    printOpenmpTraverseInfo(this, true);
-    //    buildAfterCoarsen_sebastianMode(adaptInfo, flag);
 
     clock_t first = clock();
 #ifdef _OPENMP
@@ -791,200 +790,6 @@ namespace AMDiS {
 #endif    
   }
 
-  void ProblemVec::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
-  {
-    FUNCNAME("ProblemVec::buildAfterCoarsen()");
-
-    clock_t first = clock();
-#ifdef _OPENMP
-    double wtime = omp_get_wtime();
-#endif
-
-    for (int i = 0; i < static_cast<int>(meshes.size()); i++)
-      meshes[i]->dofCompress();
-
-    Flag assembleFlag = 
-      flag | 
-      (*systemMatrix)[0][0]->getAssembleFlag() | 
-      rhs->getDOFVector(0)->getAssembleFlag()  |
-      Mesh::CALL_LEAF_EL                        | 
-      Mesh::FILL_COORDS                         |
-      Mesh::FILL_DET                            |
-      Mesh::FILL_GRD_LAMBDA |
-      Mesh::FILL_NEIGH;
-
-    if (useGetBound)
-      assembleFlag |= Mesh::FILL_BOUND;
-
-    traverseInfo.updateStatus();
-
-    // Used to calculate the overall number of non zero entries.
-    int nnz = 0;  
-
-
-    /// === INITIALIZE ===
-
-    for (int i = 0; i < nComponents; i++) {
-      MSG("%d DOFs for %s\n", 
-	  componentSpaces[i]->getAdmin()->getUsedSize(), 
-	  componentSpaces[i]->getName().c_str());
-
-      rhs->getDOFVector(i)->set(0.0);
-
-      for (int j = 0; j < nComponents; j++) {
-	// Only if this variable is true, the current matrix will be assembled.	
-	bool assembleMatrix = true;
-	// The DOFMatrix which should be assembled (or not, if assembleMatrix
-	// will be set to false).
-	DOFMatrix *matrix = (*systemMatrix)[i][j];
-
-	if (matrix) 
-	  matrix->calculateNnz();
-	
-	// If the matrix was assembled before and it is marked to be assembled
-	// only once, it will not be assembled.
-	if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
-	  assembleMatrix = false;
-	} else if (matrix) {
-	  matrix->getBaseMatrix().
-	    change_dim(componentSpaces[i]->getAdmin()->getUsedSize(), 
-		       componentSpaces[j]->getAdmin()->getUsedSize());
-
-	  set_to_zero(matrix->getBaseMatrix());	  
-	}
-
-	// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
-	if (!matrix || !assembleMatrix)
-	  assembleMatrix = false;
-
-	// If the matrix should not be assembled, the rhs vector has to be considered.
-	// This will be only done, if i == j. So, if both is not true, we can jump
-	// to the next matrix.
-	if (!assembleMatrix && i != j) {
-	  if (matrix)
-	    nnz += matrix->getBaseMatrix().nnz();
-
-	  continue;
-	}
-
-	if (assembleMatrix && matrix->getBoundaryManager())
-	  matrix->getBoundaryManager()->initMatrix(matrix);
-
-	if (matrix && assembleMatrix) 
-	  matrix->startInsertion(matrix->getNnz());
-      }
-    }
-
-
-    // === TRAVERSE ===
-
-    Mesh *mesh = componentMeshes[0];
-    const FiniteElemSpace *feSpace = componentSpaces[0];
-    const BasisFunction *basisFcts = feSpace->getBasisFcts();
-    ElementMatrix elMat(basisFcts->getNumber(), basisFcts->getNumber());
-    ElementMatrix tmpElMat(elMat);
-    ElementVector elVec(basisFcts->getNumber());
-    ElementVector tmpElVec(elVec);