Fixed bug in 3d refinement in parallel computations.

AMDiS/src/Debug.h

@@ -60,7 +60,7 @@ namespace AMDiS {
 
     /** \brief
      * Create a vtu file with name 'dofindex.vtu'. All nodes in the mesh are colored
-     * by the global dof index.
+     * by the global DOF index.
      *
      * \param[in]  feSpace  The FE space to be used.
      */

AMDiS/src/Global.h

@@ -304,6 +304,12 @@ namespace AMDiS {
   /// prints a message
 #define MSG Msg::print_funcname(funcName), Msg::print
 
+#if (DEBUG == 0) 
+  #define MSG_DBG
+#else
+  #define MSG_DBG Msg::print_funcname(funcName), Msg::print
+#endif
+
   /// prints a message, if min(Msg::msgInfo, info) >= noinfo
 #define INFO(info,noinfo)						\
   if (Msg::getMsgInfo() && (std::min(Msg::getMsgInfo(), (info)) >= (noinfo))) MSG

AMDiS/src/RefinementManager3d.cc

@@ -453,19 +453,19 @@ namespace AMDiS {
 	      break;
 	    
 	  TEST_EXIT_DBG(edgeDof0 < vertices)
- 	    ("dof %d not found on element %d with nodes (%d %d %d %d)\n", 
+ 	    ("DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
  	     edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
  	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
 
 	  TEST_EXIT_DBG(edgeDof1 < vertices)
- 	    ("dof %d not found on element %d with nodes (%d %d %d %d)\n", 
+ 	    ("DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
  	     edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
  	     neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
 	}
       }
 
       TEST_EXIT_DBG(edgeDof0 < vertices && edgeDof1 < vertices)
-	("dof %d or dof %d not found on element %d with nodes (%d %d %d %d)\n", 
+	("DOF %d or DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
 	 edge[0][0], edge[1][0], neigh->getIndex(), neigh->getDof(0, 0),
 	 neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
 
@@ -621,8 +621,6 @@ namespace AMDiS {
     if (el->getMark() <= 0)  
       return elInfo;
 
-    //    MSG("Refine: %d\n", elInfo->getElement()->getIndex());
-
     int bound = false;
     DegreeOfFreedom *edge[2];
 
@@ -681,24 +679,52 @@ namespace AMDiS {
 	  DofEdge edge2 = elInfo2->getElement()->getEdge(i);
 	  if (edge2.first == *(edge[0]) &&
 	      edge2.second == *(edge[1])) {
-	    refineList.setElType(n_neigh, elInfo2->getType());
-	    refineList.setElement(n_neigh, elInfo2->getElement());
-	    n_neigh++;
 
-	    foundEdge = true;
-	    TraverseStack *tmpStack = stack;
-	    stack = &stack2;
-	    
-	    if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
-	      getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
-	      bound = true;
-	    }
+	    // We found the edge on the other element side on leaf level. Two
+	    // cases may occure: either it is already a refinement edge, i.e,
+	    // it has local edge number 0, or it is not a refinement edge. In
+	    // the first case, we are fine and this leaf element may be set
+	    // to all elements on the refinement edge. Otherwise, the element
+	    // must be refined at least once to get a refinement edge.
+
+	    if (i == 0) {
+	      // Edge is refinement edge, so add it to refine list.
+
+	      refineList.setElType(n_neigh, elInfo2->getType());
+	      refineList.setElement(n_neigh, elInfo2->getElement());
+	      n_neigh++;
+	      
+	      foundEdge = true;
+	      TraverseStack *tmpStack = stack;
+	      stack = &stack2;
+	      
+	      if (getRefinePatch(&elInfo2, edge, 0, refineList, &n_neigh)) {
+		getRefinePatch(&elInfo2, edge, 1, refineList, &n_neigh);
+		bound = true;
+	      }
+	      
+	      stack = tmpStack;
+	      break;
+	    } else {
+	      // Edge i not refinement edge, so refine the element further.
 
-	    stack = tmpStack;
-	    break;
+	      Element *el2 = elInfo->getElement();
+	      el2->setMark(std::max(el2->getMark(), 1));
+
+	      TraverseStack *tmpStack = stack;
+	      stack = &stack2;
+
+	      elInfo2 = refineFunction(elInfo2);
+
+	      stack = tmpStack;
+	      break;
+	    }
 	  }
 	}
 
+	if (foundEdge)
+	  break;
+
 	elInfo2 = stack2.traverseNext(elInfo2);
       }
 
@@ -778,8 +804,8 @@ namespace AMDiS {
 
 
   void FixRefinementPatch::getOtherEl(TraverseStack *stack, 
-				     Element **otherEl, 
-				     int &otherEdge)
+				      Element **otherEl, 
+				      int &otherEdge)
   {
     FUNCNAME("FixRefinementPatch::getOtherEl()");
 

AMDiS/src/parallel/MeshDistributor.cc

@@ -170,7 +170,7 @@ namespace AMDiS {
       }
 
 #if (DEBUG != 0)
-    ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
+      ParallelDebug::writeDebugFile(*this, debugOutputDir + "mpi-dbg", "dat");
 #endif    
 
       initialized = true;
@@ -558,6 +558,7 @@ namespace AMDiS {
 
 
     FixRefinementPatch::connectedEdges.clear();
+    bool valid3dMesh = true;
 
     for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
       vector<ElementObjectData>& edgeEls = elObjects.getElements(*it);
@@ -603,8 +604,10 @@ namespace AMDiS {
       } while (found);
       
       if (!el1.empty()) {
-	MSG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", 
-	    edgeEls.size(), el0.size(), el1.size());
+	valid3dMesh = false;
+
+	MSG_DBG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", 
+		edgeEls.size(), el0.size(), el1.size());
 
 	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
 	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
@@ -626,10 +629,10 @@ namespace AMDiS {
 	    mesh->getDofIndexCoords(dofEdge0.first, feSpace, c0);
 	    mesh->getDofIndexCoords(dofEdge0.second, feSpace, c1);
 
-	    MSG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
-		dofEdge0.first, dofEdge0.second, 
-		c0[0], c0[1], c0[2],
-		c1[0], c1[1], c1[2]);
+	    MSG_DBG("FOUND EDGE: %d %d with coords %f %f %f   %f %f %f\n",
+		    dofEdge0.first, dofEdge0.second, 
+		    c0[0], c0[1], c0[2],
+		    c1[0], c1[1], c1[2]);
 
 	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
 
@@ -639,6 +642,8 @@ namespace AMDiS {
 	}
       }
     }
+
+    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
   }
 
 
@@ -768,9 +773,7 @@ namespace AMDiS {
 
 
       // === Check the boundaries and adapt mesh if necessary. ===
-#if (DEBUG != 0)
-      MSG("Run checkAndAdaptBoundary ...\n");
-#endif
+      MSG_DBG("Run checkAndAdaptBoundary ...\n");
 
       meshChanged |= checkAndAdaptBoundary(allBound);
 

AMDiS/src/parallel/MeshPartitioner.cc

@@ -47,7 +47,7 @@ namespace AMDiS {
 	vertexElements[el->getEdge(0)].insert(elIndex);
       } else {
 	int elInRank = std::min(elIndex / elPerRank, mpiSize - 1);
-	
+
 	elementInRank[elIndex] = (elInRank == mpiRank);
 	partitionMap[elIndex] = elInRank;	
       }

AMDiS/src/parallel/ParallelDebug.cc

@@ -805,6 +805,9 @@ namespace AMDiS {
 
     // === Write to all elements in ranks mesh the included dofs. ===
 
+    file << "\n\n";
+    file << "# First line containes number of elements in mesh, second line contain the number of DOFs per element.\n";
+    file << "# Than, each entry contains of two lines. The first is the element index, the second line is a list with the local DOF indices of this element.\n";
     file << elDofMap.size() << "\n";
     file << basisFcts->getNumber() << "\n";
     for (ElDofMap::iterator it = elDofMap.begin(); it != elDofMap.end(); ++it) {