Fixed a small bug for parallel computations + periodic boundaries + global refinements.

AMDiS/src/Mesh.cc

@@ -1429,6 +1429,7 @@ namespace AMDiS {
   {
     map<BoundaryType, VertexVector*>::iterator it;
     map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
+    
     for (it = periodicAssociations.begin(); it != end; ++it)
       if ((*(it->second))[dof1] == dof2)
 	return true;

AMDiS/src/RCNeighbourList.h

@@ -150,22 +150,16 @@ namespace AMDiS {
       coarseningManager = cm; 
     }
 
-    /** \brief
-     * Fills \ref rclist[i].neigh and \ref rclist[i].oppVertex  infos (0 <= i 
-     * < n_neigh)
-     */
+    /// Fills \ref rclist[i].neigh and \ref rclist[i].oppVertex  infos (0 <= i 
+    /// < n_neigh)
     void fillNeighbourRelations(int n_neigh, int bound);
-
-    /** \brief
-     * Adds those dof's on the parent that are handed on by the
-     * children and adds the dof in the midpoint of the coarsening edge (3d) 
-     */
+    
+    /// Adds those dof's on the parent that are handed on by the
+    /// children and adds the dof in the midpoint of the coarsening edge (3d) 
     void addDOFParent(int elIndex, DegreeOfFreedom* dof);
 
-    /** \brief
-     * If DOFs for higher order are been removed on parents during refinement
-     * they are now added again (2d)
-     */
+    /// If DOFs for higher order are been removed on parents during refinement
+    /// they are now added again (2d)
     void addDOFParents(int n_neigh);
 
     /// Removes DOFs during refinement (3d)

AMDiS/src/RefinementManager.cc

@@ -63,12 +63,14 @@ namespace AMDiS {
       ElInfo *elInfo = 
 	stack->traverseFirst(mesh, -1, 
 			     Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND);
+
       while (elInfo) {
 	if (elInfo->getElement()->getMark() > 0) {	  
 	  doMoreRecursiveRefine = 
 	    doMoreRecursiveRefine || (elInfo->getElement()->getMark() > 1);
 	  elInfo = refineFunction(elInfo);
 	}	
+
 	elInfo = stack->traverseNext(elInfo);
       }
     }

AMDiS/src/RefinementManager3d.cc

@@ -409,14 +409,14 @@ namespace AMDiS {
     int localNeighbour = 3 - direction;
     Tetrahedron *el = 
       dynamic_cast<Tetrahedron*>(const_cast<Element*>((*elInfo)->getElement()));
-    
+
     if ((*elInfo)->getNeighbour(localNeighbour) == NULL)
       return true;    
   
     int oppVertex = (*elInfo)->getOppVertex(localNeighbour);
-    #if DEBUG
+#if DEBUG
     int testIndex = (*elInfo)->getNeighbour(localNeighbour)->getIndex();
-    #endif
+#endif
     ElInfo *neighInfo = stack->traverseNeighbour3d((*elInfo), localNeighbour);
     int neighElType = neighInfo->getType();
 
@@ -562,9 +562,9 @@ namespace AMDiS {
 	if (neighInfo->getNeighbour(i)) {
 	  oppVertex = neighInfo->getOppVertex(i);
     
-    #if DEBUG
+#if DEBUG
 	  int testIndex = neighInfo->getNeighbour(i)->getIndex();
-    #endif
+#endif
 
 	  neighInfo = stack->traverseNeighbour3d(neighInfo, i);
 
@@ -596,9 +596,9 @@ namespace AMDiS {
 	("While looping back domains boundary was reached or i == 0\n");
       oppVertex = refineList.getOppVertex(i--, 0);
 
-      #if DEBUG
+#if DEBUG
       int testIndex = neighInfo->getNeighbour(oppVertex)->getIndex();
-      #endif
+#endif
       neighInfo = stack->traverseNeighbour3d(neighInfo, oppVertex);
 
       int edgeDof0, edgeDof1;
@@ -653,7 +653,7 @@ namespace AMDiS {
       edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
     }
 
-#if HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     Element *otherEl = NULL;
     int otherEdge = -1;
     FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
@@ -667,7 +667,7 @@ namespace AMDiS {
       bound = true;
     }
 
-#if HAVE_PARALLEL_DOMAIN_AMDIS
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
     // === If the refinement edge must be fixed, add also the other part of this ===
     // === edge to the refinement patch.                                         ===
 
@@ -739,11 +739,9 @@ namespace AMDiS {
     }
 #endif
 
-
     // fill neighbour information inside the patch in the refinement list
     refineList.fillNeighbourRelations(n_neigh, bound);
 
-
     // ============ Check for periodic boundary ============
 
     DegreeOfFreedom *next_edge[2];
@@ -803,7 +801,7 @@ namespace AMDiS {
 	}
       }
     }
-  
+
     stack->update();
 
     return elInfo;

AMDiS/src/Traverse.cc

@@ -525,9 +525,9 @@ namespace AMDiS {
   {
     FUNCNAME("TraverseStack::traverseNeighbour3d()");
 
-    Element *el2;
-    ElInfo *elinfo2;
-    int stack2_used;
+    Element *el2 = NULL;
+    ElInfo *elinfo2 = NULL;
+    int stack2_used = 0;
     int sav_neighbour = neighbour;
 
     // father.neigh[coarse_nb[i][j]] == child[i - 1].neigh[j]
@@ -611,10 +611,12 @@ namespace AMDiS {
       if (traverse_mel == NULL)  
 	return NULL;
     
-      if (nb < 2 && save_stack_used > 1)
-	stack2_used = 2;                /* go down one level in OLD hierarchy */
-      else
+      if (nb < 2 && save_stack_used > 1) {
+	// go down one level in OLD hierarchy
+	stack2_used = 2;                
+      } else {
 	stack2_used = 1;
+      }
       
       elinfo2 = save_elinfo_stack[stack2_used];
       el2 = elinfo2->getElement();
@@ -625,9 +627,11 @@ namespace AMDiS {
     } else {                                                
       // Goto other child.
 
-      stack2_used = stack_used + 1;
-      if (save_stack_used > stack2_used)
-	stack2_used++;                  /* go down one level in OLD hierarchy */
+      stack2_used = stack_used + 1;      
+      if (save_stack_used > stack2_used) {
+	// go down one level in OLD hierarchy
+	stack2_used++;
+      }
 
       elinfo2 = save_elinfo_stack[stack2_used];
       el2 = elinfo2->getElement();
@@ -680,7 +684,8 @@ namespace AMDiS {
       if (save_stack_used > stack2_used) { 
 	// `refine' both el and el2.
 
-	TEST_EXIT_DBG(el->getChild(0))("invalid new refinement?\n");
+	TEST_EXIT_DBG(el->getChild(0))
+	  ("Element %d has no children!\n", el->getIndex());
 
 	int i = 0;
 	if (el->getDof(0) == el2->getDof(0))
@@ -730,11 +735,12 @@ namespace AMDiS {
 
 	if (save_stack_used > stack2_used) {
 	  const DegreeOfFreedom *dof = el2->getDof(1);
+
 	  if (dof != el->getDof(1) && dof != el->getDof(2) &&
 	      !traverse_mesh->associated(dof[0], el->getDof(1, 0)) &&
 	      !traverse_mesh->associated(dof[0], el->getDof(2, 0))) {
-	  
-	    stack2_used++;               /* go down one level in OLD hierarchy */
+	    // go down one level in OLD hierarchy
+	    stack2_used++;               
 	    elinfo2 = save_elinfo_stack[stack2_used];
 	    el2 = elinfo2->getElement();
 	  } 
@@ -782,9 +788,9 @@ namespace AMDiS {
   {
     FUNCNAME("TraverseStack::traverseNeighbour2d()");
 
-    Triangle *el2;
-    ElInfo *elinfo2;
-    int stack2_used;
+    Triangle *el2 = NULL;
+    ElInfo *elinfo2 = NULL;
+    int stack2_used = 0;
     int sav_neighbour = neighbour;
 
     // father.neigh[coarse_nb[i][j]] == child[i-1].neigh[j]

AMDiS/src/parallel/MeshDistributor.cc

@@ -314,10 +314,12 @@ namespace AMDiS {
     MSG("Debug mode tests finished!\n");
 #endif
 
-    // === Create periodic DOF mapping, if there are periodic boundaries. ===
-
+    // Create periodic DOF mapping, if there are periodic boundaries.
     createPeriodicMap();
 
+    // Remove periodic boundary conditions in sequential problem definition. 
+    removePeriodicBoundaryConditions();
+
 #if (DEBUG != 0)
     ParallelDebug::testPeriodicBoundary(*this);
 #endif
@@ -331,11 +333,10 @@ namespace AMDiS {
       refineManager->globalRefine(mesh, globalRefinement);
 
       updateLocalGlobalNumbering();
-     
+
       // === Update periodic mapping, if there are periodic boundaries. ===     
 
       createPeriodicMap();
-
 #if (DEBUG != 0)
     ParallelDebug::testPeriodicBoundary(*this);
 #endif
@@ -344,9 +345,6 @@ namespace AMDiS {
     // Set DOF rank information to all matrices and vectors.
     setRankDofs();
 
-    // Remove periodic boundary conditions in sequential problem definition.
-    removePeriodicBoundaryConditions();
-
     initialized = true;
   }