Fixed fixing 3d mesh refinement for very complicated macro mesh partitionings.

AMDiS/src/RefinementManager3d.cc

@@ -684,7 +684,6 @@ namespace AMDiS {
 
     for (int edgeIndex = 0; 
  	 edgeIndex < static_cast<unsigned int>(refineEdges.size()); edgeIndex++) {
-
       MSG("  MAKE -> %d\n", edgeIndex);
 
       Element *otherEl = refineEdges[edgeIndex].first;   

AMDiS/src/parallel/MeshDistributor.cc

@@ -291,11 +291,10 @@ namespace AMDiS {
 
     updateLocalGlobalNumbering();
 
-    // === In 3D we have to make some test, if the resulting mesh is valid.  ===
-    // === If it is not valid, there is no possiblity yet to fix this        ===
-    // === problem, just exit with an error message.                         ===
 
-    check3dValidMesh();
+    // === In 3D we have to fix the mesh to allow local refinements. ===
+
+    fix3dMeshRefinement();
 
 
     // === If in debug mode, make some tests. ===
@@ -658,13 +657,17 @@ namespace AMDiS {
   }
 
 
-  void MeshDistributor::check3dValidMesh()
+  void MeshDistributor::fix3dMeshRefinement()
   {
-    FUNCNAME("MeshDistributor::check3dValidMesh()");
+    FUNCNAME("MeshDistributor::fix3dMeshRefinement()");
 
     if (mesh->getDim() != 3)
       return;
 
+    MPI::COMM_WORLD.Barrier();
+    double first = MPI::Wtime();
+
+
     // === Create set of all edges and all macro element indices in rank's ===
     // === subdomain.                                                      ===
 
@@ -674,7 +677,7 @@ namespace AMDiS {
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
     while (elInfo) {
-      for (int i = 0; i < 6; i++) {
+      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
 	ElementObjectData elData(elInfo->getElement()->getIndex(), i);
 	DofEdge e = elObjDb.getEdgeLocalMap(elData);
 	allEdges.insert(e);
@@ -686,11 +689,10 @@ namespace AMDiS {
     }
 
 
-    // === Reset fixed refinement edges and assume the mesh is valid. Search ===
-    // === for the opposite.                                                 ===
+    // === Reset fixed edges, and start new search for edges that must  ===
+    // === be fixed.                                                    ===
 
     FixRefinementPatch::connectedEdges.clear();
-    bool valid3dMesh = true;
 
     // Traverse all edges
     for (std::set<DofEdge>::iterator it = allEdges.begin(); it != allEdges.end(); ++it) {
@@ -702,62 +704,48 @@ namespace AMDiS {
 
 
       // We have now to check, if the current edge connects two macro elements,
-      // which are otherwise not connected. The basic idea to check this is very
-      // simple: We take the first macro element in rank's subdomain that contain
-      // this edge and add it to the set variable "el0". All other macro elements
-      // which share this edge are added to "el1". 
+      // which are otherwise not connected. 
 
-      std::set<int> el0, el1;
+      // All elements on this edge.
+      std::set<int> allEls;
+      // Maps from each element index to the local edge index of the common edge.
       map<int, int> edgeNoInEl;
 
       for (unsigned int i = 0; i < edgeEls.size(); i++) {
 	if (rankMacroEls.count(edgeEls[i].elIndex)) {
-	  if (el0.empty())
-	    el0.insert(edgeEls[i].elIndex);
-	  else
-	    el1.insert(edgeEls[i].elIndex);
-
+	  allEls.insert(edgeEls[i].elIndex);
 	  edgeNoInEl[edgeEls[i].elIndex] = edgeEls[i].ithObject;
 	}
       }
            
-      TEST_EXIT_DBG(!el0.empty())("Should not happen!\n");
-
-      // If el1 is empty, there is only on macro element in the mesh which
-      // contains this edge. Hence, we can continue with checking another edge.
-      if (el1.empty())
-	continue;
+      TEST_EXIT_DBG(!allEls.empty())("Should not happen!\n");
 
-      bool found = false;      
-      do {
-	found = false;
-	for (std::set<int>::iterator elIt = el0.begin(); elIt != el0.end(); ++elIt) {
-	  for (int i = 0; i < 4; i++) {
-	    int neighEl = macroElementNeighbours[*elIt][i];
-	    if (neighEl != -1 && el1.count(neighEl)) {
-	      el0.insert(neighEl);
-	      el1.erase(neighEl);
-	      found = true;
-	    }
-	  }
 
-	  if (found)
-	    break;
-	}
-      } while (found);
+      // If there is only one element in rank on this edge, there is nothing 
+      // to do.
+      if (allEls.size() == 1)
+	continue;
       
-      if (!el1.empty()) {
-	valid3dMesh = false;
+      // Create a set of element index sets. All element within one set are 
+      // neighbours of each other, but can not reach any other element of the
+      // other index sets by neighbourhood relations.
+      vector<std::set<int> > disconnectedEls;
+      helpToFix(allEls, edgeNoInEl, disconnectedEls);   
       
-	MSG_DBG("Unvalid 3D mesh with %d (%d - %d) elements on this edge!\n", 
-		edgeEls.size(), el0.size(), el1.size());
+      // Fix the edges
+      for (int dc0 = 0; dc0 < static_cast<int>(disconnectedEls.size()); dc0++) {
+	for (int dc1 = 0; dc1 < static_cast<int>(disconnectedEls.size()); dc1++) {
+	  if (dc0 == dc1)
+	    continue;   
 	  
-	for (std::set<int>::iterator elIt0 = el0.begin(); elIt0 != el0.end(); ++elIt0) {
-	  for (std::set<int>::iterator elIt1 = el1.begin(); elIt1 != el1.end(); ++elIt1) {
+	  int elIdx = *(disconnectedEls[dc1].begin());
 	  pair<Element*, int> edge0 = 
-	      make_pair(elObjDb.getElementPtr(*elIt0), edgeNoInEl[*elIt0]);
+	    make_pair(elObjDb.getElementPtr(elIdx), edgeNoInEl[elIdx]);
+	  
+	  for (std::set<int>::iterator elIt = disconnectedEls[dc0].begin();
+	       elIt != disconnectedEls[dc0].end(); ++elIt) {
 	    pair<Element*, int> edge1 = 
-	      make_pair(elObjDb.getElementPtr(*elIt1), edgeNoInEl[*elIt1]);
+	      make_pair(elObjDb.getElementPtr(*elIt), edgeNoInEl[*elIt]);
 
 #if (DEBUG != 0)
 	    DofEdge dofEdge0 = edge0.first->getEdge(edge0.second);
@@ -767,26 +755,60 @@ namespace AMDiS {
 	    mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0);
 	    mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1);
 	    
-	    MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
-		edge0.first->getIndex(), edge0.second,
-		edge1.first->getIndex(), edge1.second);
+	    // 	      MSG("FOUND EDGE %d/%d <-> %d/%d\n", 
+	    // 		  edge0.first->getIndex(), edge0.second,
+	    // 		  edge1.first->getIndex(), edge1.second);
 	    
-	    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]);
-#endif
+	    // 	      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]);
 	    
 	    TEST_EXIT_DBG(dofEdge0 == dofEdge1)("Should noth happen!\n");
+#endif
 
-	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge0, edge1));
 	    FixRefinementPatch::connectedEdges.push_back(make_pair(edge1, edge0));
 	  }
 	}	  
       }
     }
 
-    MSG_DBG("Mesh is 3d valid mesh: %d\n", valid3dMesh);
+
+    MPI::COMM_WORLD.Barrier();
+    MSG("Fix 3D mesh refinement needed %.5f seconds\n", MPI::Wtime() - first);
+  }
+
+
+  void MeshDistributor::helpToFix(std::set<int> &elems, 
+				  map<int, int> &edgeNoInEl,
+				  vector<std::set<int> > &disconnectedEls)
+  {
+    std::set<int> firstElem;
+    firstElem.insert(*(elems.begin()));
+
+    std::set<int> otherElems(++(elems.begin()), elems.end());
+    
+    bool found = false;      
+    do {
+      found = false;
+      for (std::set<int>::iterator elIt = firstElem.begin(); elIt != firstElem.end(); ++elIt) {
+	for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
+	  int neighEl = macroElementNeighbours[*elIt][i];
+	  if (neighEl != -1 && otherElems.count(neighEl)) {
+	    firstElem.insert(neighEl);
+	    otherElems.erase(neighEl);
+	    found = true;
+	  }
+	}
+	
+	if (found)
+	  break;
+      }
+    } while (found);
+
+    disconnectedEls.push_back(firstElem);
+    if (otherElems.size())
+      helpToFix(otherElems, edgeNoInEl, disconnectedEls);
   }
 
 
@@ -1614,7 +1636,7 @@ namespace AMDiS {
 #endif
 
     // In 3D we have to make some test, if the resulting mesh is valid.
-    check3dValidMesh();
+    fix3dMeshRefinement();
 
 
     MPI::COMM_WORLD.Barrier();

AMDiS/src/parallel/MeshDistributor.h

@@ -282,7 +282,23 @@ namespace AMDiS {
     /// only on one edge are not neighours by definition. This functions checks
     /// for this situation and fix the problem. For this, the mesh is search for
     /// all edges connecting two elements that are otherwise not connected.
-    void check3dValidMesh();
+    void fix3dMeshRefinement();
+
+    /** \brief Is used only within \ref fix3dMeshRefinement.
+     * 
+     * \param[in]  elems            Set of macro element indices.
+     * \param[in]  edgeNoInEl       Maps from each element index in elems to a
+     *                              local edge number.
+     * \param[out] disconnectedEls  On output, this vector contains sets of
+     *                              element indices. The union is equal to elems.
+     *                              Each set contains all element indices, which
+     *                              are reachable among each other by neighbour
+     *                              relations. Elements within two different sets
+     *                              cannot be reached via neigbourhood relation.
+     */
+    void helpToFix(std::set<int> &elems, 
+		   map<int, int> &edgeNoInEl,
+		   vector<std::set<int> > &disconnectedEls);
 
     void setBoundaryDofRequirement(Flag flag)
     {