Fixed 3D refinement problem in parallel computations with mesh repartitioning.

AMDiS/src/CoarseningManager3d.cc

@@ -19,36 +19,40 @@
 #include "RCNeighbourList.h"
 #include "FixVec.h"
 #include "DOFIndexed.h"
+#include "Debug.h"
 
 namespace AMDiS {
 
-  void CoarseningManager3d::coarsenFunction(ElInfo *el_info)
+  void CoarseningManager3d::coarsenFunction(ElInfo *elInfo)
   {
+    FUNCNAME("CoarseningManager3d::coarsenFunction()");
+
     Tetrahedron *el = 
-      dynamic_cast<Tetrahedron*>(const_cast<Element*>(el_info->getElement()));
+      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
 
+    // If element must not be coarsend, return.
     if (el->getMark() >= 0) 
-      return; // el must not be coarsend, return :-(
+      return; 
 
+    // Single leaves don't get coarsened.
     if (el->isLeaf()) 
-      return;       // single leaves don't get coarsened
+      return;
 
     if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) {
-      // one of the children must not be coarsend; return :-(
+      // One of the children must not be coarsend, so return.
 
       el->setMark(0);
       return;
     }
 
     if (!(el->getChild(0)->isLeaf()) || !(el->getChild(1)->isLeaf())) {
-      // one of the children is not a leaf element; try again later on
+      // One of the children is not a leaf element. Try again later on.
       doMore = true;
       return;
     }
 
     DegreeOfFreedom *edge[2];
     int n_neigh, bound = 0;
-    ElInfo *elinfo = el_info;
 
     /****************************************************************************/
     /*  get a list for storing all elements at the coarsening edge and fill it  */
@@ -60,24 +64,95 @@ namespace AMDiS {
     /*  give the refinement edge the right orientation                          */
     /****************************************************************************/
 
-    if (el->getDof(0,0) < el->getDof(1,0)) {
-      edge[0] = const_cast<int*>(el->getDof(0));
-      edge[1] = const_cast<int*>(el->getDof(1));
+    if (el->getDof(0, 0) < el->getDof(1, 0)) {
+      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(0));
+      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(1));
     } else {
-      edge[1] = const_cast<int*>(el->getDof(0));
-      edge[0] = const_cast<int*>(el->getDof(1));
+      edge[1] = const_cast<DegreeOfFreedom*>(el->getDof(0));
+      edge[0] = const_cast<DegreeOfFreedom*>(el->getDof(1));
     }
 
     coarsenList.setElement(0, el, true);
     n_neigh = 1;
 
     coarsenList.setOppVertex(0, 0, 0);
-    coarsenList.setElType(0, el_info->getType());
+    coarsenList.setElType(0, elInfo->getType());
     bound = false;
-    if (getCoarsenPatch(elinfo, edge, 0, coarsenList, &n_neigh)) {
-      getCoarsenPatch(elinfo, edge, 1, coarsenList, &n_neigh);
+    if (getCoarsenPatch(elInfo, edge, 0, coarsenList, &n_neigh)) {
+      getCoarsenPatch(elInfo, edge, 1, coarsenList, &n_neigh);
       bound = true;
     }
+
+
+#if HAVE_PARALLEL_DOMAIN_AMDIS
+    Element *otherEl = NULL;
+    int otherEdge = -1;      
+    FixRefinementPatch::getOtherEl(stack, &otherEl, otherEdge);
+    
+    // === If the refinement edge must be fixed, add also the other part of this ===
+    // === edge to the refinement patch.                                         ===
+    
+    if (otherEl) {
+      // TODO: Remove these two lines and make something more meaningful!!
+      el->setMark(0);
+      return;
+      
+      TraverseStack stack2;
+      ElInfo *elInfo2 = 
+	stack2.traverseFirstOneMacro(mesh, otherEl->getIndex(), -1, 
+				     Mesh::CALL_EVERY_EL_PREORDER | 
+				     Mesh::FILL_NEIGH |
+				     Mesh::FILL_BOUND);
+      
+      bool foundEdge = false;
+      
+      while (elInfo2) {
+	Element *el2 = elInfo2->getElement();
+	for (int i = 0; i < 6; i++) {
+	  DofEdge edge2 = elInfo2->getElement()->getEdge(i);
+	  if (edge2.first == *(edge[0]) && edge2.second == *(edge[1]) && !el2->isLeaf()) {
+	    
+	    if (!el2->getChild(0)->isLeaf() || !el2->getChild(1)->isLeaf()) {
+	      int edgeNo0 = el->getEdgeOfChild(0, i, elInfo2->getType());
+	      int edgeNo1 = el->getEdgeOfChild(1, i, elInfo2->getType());
+	      
+	      bool refineChildFirst = 
+		!(i > 0 &&
+		  (edgeNo0 >= 0 && !el2->getChild(0)->isLeaf()) || 
+		  (edgeNo1 >= 0 && !el2->getChild(1)->isLeaf()));
+	      
+	      if (refineChildFirst) {
+		// TODO: WAS SOLL ICH NUR MACHEN???
+		el->setMark(0);
+		return;
+	      }
+	    } else {
+	      coarsenList.setElType(n_neigh, elInfo2->getType());
+	      coarsenList.setElement(n_neigh, elInfo2->getElement(), true);
+	      n_neigh++;
+	      
+	      foundEdge = true;
+	      
+	      TraverseStack *tmpStack = stack;
+	      stack = &stack2;
+	      
+	      if (getCoarsenPatch(elInfo2, edge, 0, coarsenList, &n_neigh)) {
+		getCoarsenPatch(elInfo2, edge, 1, coarsenList, &n_neigh);
+		bound = true;
+	      }
+	      stack = tmpStack;
+	      break;
+	    }
+	  }
+	}
+	
+	elInfo2 = stack2.traverseNext(elInfo2);
+      }
+      
+      TEST_EXIT_DBG(foundEdge)("Should not happen!\n");
+    }  
+#endif
+    
     coarsenList.fillNeighbourRelations(n_neigh, bound);
 
     /****************************************************************************/
@@ -95,8 +170,8 @@ namespace AMDiS {
 
       while (edge[0] != NULL) {
 	coarsenList.periodicSplit(edge, next_edge,
-				   &n_neigh, &n_neigh_periodic,
-				   periodicList);
+				  &n_neigh, &n_neigh_periodic,
+				  periodicList);
 
 	coarsenPatch(periodicList, n_neigh_periodic, bound);
 
@@ -119,22 +194,19 @@ namespace AMDiS {
     Tetrahedron *el = 
       dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getElement(index)));
     Tetrahedron *child[2];
-    Tetrahedron *neigh;
-    int dir, el_type, i, node, opp_v;
 
     child[0] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(0)));
     child[1] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el->getChild(1)));
-    el_type = coarsenList.getType(index);
+    int el_type = coarsenList.getType(index);
 
     /****************************************************************************/
     /*  Information about patch neighbours is still valid! But edge and face    */
     /*  dof's in a common face of patch neighbours have to be removed           */
     /****************************************************************************/
 
-    for (dir = 0; dir < 2; dir++) {
-      neigh = 
+    for (int dir = 0; dir < 2; dir++) {
+      Tetrahedron *neigh = 
 	dynamic_cast<Tetrahedron*>(const_cast<Element*>(coarsenList.getNeighbourElement(index, dir)));
-      opp_v = coarsenList.getOppVertex(index, dir);
 
       if (!neigh || neigh->isLeaf()) {
 	/****************************************************************************/
@@ -142,11 +214,11 @@ namespace AMDiS {
 	/****************************************************************************/
 
 	if (mesh->getNumberOfDofs(EDGE)) {
-	  node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir];
+	  int node = mesh->getNode(EDGE) + Tetrahedron::nChildEdge[el_type][0][dir];
 	  mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), EDGE);
 	}
 	if (mesh->getNumberOfDofs(FACE)) {
-	  node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
+	  int node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][0][dir];
 	  mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
 	  node = mesh->getNode(FACE) + Tetrahedron::nChildFace[el_type][1][dir];
 	  mesh->freeDof(const_cast<int*>(child[1]->getDof(node)), FACE);
@@ -160,14 +232,14 @@ namespace AMDiS {
     /****************************************************************************/
 
     if (mesh->getNumberOfDofs(FACE)) {
-      node = mesh->getNode(FACE);
+      int node = mesh->getNode(FACE);
       mesh->freeDof(const_cast<int*>(child[0]->getDof(node)), FACE);
     }
 
 
     if (mesh->getNumberOfDofs(CENTER)) {
-      node = mesh->getNode(CENTER);
-      for (i = 0; i < 2; i++)
+      int node = mesh->getNode(CENTER);
+      for (int i = 0; i < 2; i++)
 	mesh->freeDof(const_cast<int*>(child[i]->getDof(node)), CENTER);
     }
 
@@ -191,7 +263,7 @@ namespace AMDiS {
 
   /****************************************************************************/
   /*  get_coarse_patch:  gets the patch for coarsening starting on element    */
-  /*  el_info->el in direction of neighbour [3-dir]; returns 1 if a boundary  */
+  /*  elInfo->el in direction of neighbour [3-dir]; returns 1 if a boundary  */
   /*  reached and 0 if we come back to the starting element.                  */
   /*                                                                          */
   /*  if NEIGH_IN_EL we only can find the complete coarsening patch if the    */
@@ -205,16 +277,13 @@ namespace AMDiS {
   /*  the starting element before we return                                   */
   /****************************************************************************/
  
-  bool CoarseningManager3d::getCoarsenPatch(ElInfo *el_info, 
+  bool CoarseningManager3d::getCoarsenPatch(ElInfo *elInfo, 
 					    DegreeOfFreedom *edge[2],
 					    int dir, 
 					    RCNeighbourList &coarsenList, 
 					    int *n_neigh)
   {
-    FUNCNAME("CoarseningManager3d::getCoarsenPatch");
-    ElInfo *neigh_info;
-    Tetrahedron *el, *neigh;
-    int i, j, k, opp_v, edge_no;
+    FUNCNAME("CoarseningManager3d::getCoarsenPatch()");
 
     static unsigned char next_el[6][2] = {{3,2},
 					  {1,3},
@@ -223,15 +292,16 @@ namespace AMDiS {
 					  {0,2},
 					  {0,1}};
 
-    el = dynamic_cast<Tetrahedron*>(const_cast<Element*>(el_info->getElement()));
-    neigh = 
-      dynamic_cast<Tetrahedron*>(const_cast<Element*>(el_info->getNeighbour(3 - dir)));
+    Tetrahedron *el = 
+      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getElement()));
+    Tetrahedron *neigh = 
+      dynamic_cast<Tetrahedron*>(const_cast<Element*>(elInfo->getNeighbour(3 - dir)));
     if (neigh == NULL)
       return true;
 
-    opp_v = el_info->getOppVertex(3 - dir);
+    int opp_v = elInfo->getOppVertex(3 - dir);
+    ElInfo *neigh_info = stack->traverseNeighbour3d(elInfo, 3 - dir);
 
-    neigh_info = stack->traverseNeighbour3d(el_info, 3 - dir);
     TEST_EXIT_DBG(neigh == neigh_info->getElement())
       ("neigh %d and neigh_info->el %d are not identical\n", 
        neigh->getIndex(), neigh_info->getElement()->getIndex());
@@ -244,6 +314,7 @@ namespace AMDiS {
     coarsenList.setElType(*n_neigh, neigh_info->getType()); 
   
     int n_vertices = mesh->getGeo(VERTEX);
+    int i, j, k, edge_no;
 			 
     while (neigh != el) {
       for (j = 0; j < n_vertices; j++)
@@ -261,10 +332,15 @@ namespace AMDiS {
 	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
 	    break;
      
-	TEST_EXIT_DBG(j < n_vertices  &&  k < n_vertices)
-	  ("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));
+	TEST_EXIT_DBG(j < n_vertices)
+	  ("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(k < n_vertices)
+	  ("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));
       }
       edge_no = Tetrahedron::edgeOfDofs[j][k];
       coarsenList.setCoarsePatch(*n_neigh, edge_no == 0);
@@ -366,6 +442,12 @@ namespace AMDiS {
 
     // === And now start to coarsen the patch: remove dof's of the coarsening edge. ===
 
+    WorldVector<double> c;
+    FiniteElemSpace *feSpace = FiniteElemSpace::provideFeSpace(mesh);
+    mesh->getDofIndexCoords(el->getChild(0)->getDof(3), feSpace, c);
+
+    //    MSG("Free dof %p   %f %f %f\n", el->getChild(0)->getDof(3), c[0], c[1], c[2]);
+
     mesh->freeDof(const_cast<int*>(el->getChild(0)->getDof(3)), VERTEX);
     mesh->incrementNumberOfVertices(-1);
 

AMDiS/src/Debug.cc

@@ -74,6 +74,19 @@ namespace AMDiS {
     }
 
 
+    void colorEdgeInMesh(FiniteElemSpace *feSpace,
+			 Element *el, 
+			 int localEdgeNo, 
+			 std::string filename)
+    {
+      DOFVector<double> tmp(feSpace, "tmp");
+      tmp.set(0.0);
+      tmp[el->getEdge(localEdgeNo).first] = 1.0;
+      tmp[el->getEdge(localEdgeNo).second] = 1.0;
+      VtkWriter::writeFile(tmp, filename);
+    }
+
+
     void writeElementIndexMesh(Mesh *mesh, std::string filename, int level)
     {
       FUNCNAME("debug::writeElementIndexMesh()");

AMDiS/src/Debug.h

@@ -66,6 +66,11 @@ namespace AMDiS {
      */
     void writeDofIndexMesh(FiniteElemSpace *feSpace);
 
+    void colorEdgeInMesh(FiniteElemSpace *feSpace,
+			 Element *el, 
+			 int localEdgeNo, 
+			 std::string filename);
+
     /** \brief
      * Creates a vtu file where all elements in the mesh are colored by the global
      * element indices.

AMDiS/src/ElInfo3d.cc

@@ -483,7 +483,7 @@ namespace AMDiS {
       normal[1] /= det;
       normal[2] /= det;
     } else {
-      MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld);
+      MSG("Not implemented for DIM_OF_WORLD = %d in 3D!\n", dimOfWorld);
     }
 
     return det;
@@ -593,37 +593,11 @@ namespace AMDiS {
 	if ((nb = const_cast<Element*>((*neigh_old)[cv[i]]))) {
 	  TEST_EXIT_DBG(nb->getChild(0))("nonconforming triangulation\n");
 
-	  if (elOld->getIndex() == 15376 && ichild == 0) {
-	    MSG("NEIGH %d OF EL is %d\n", i, nb->getIndex());
-	    MSG("EL HAS DOFS: %d %d %d %d\n",
-		elOld->getDof(0, 0),
-		elOld->getDof(1, 0),
-		elOld->getDof(2, 0),
-		elOld->getDof(3, 0));
-	    MSG("NEIGH-EL HAS DOFS: %d %d %d %d\n",
-		nb->getDof(0, 0),
-		nb->getDof(1, 0),
-		nb->getDof(2, 0),
-		nb->getDof(3, 0));
-	  }
-
-	  
 	  int k;
 	  for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */
 	    nbk = const_cast<Element*>(nb->getChild(k));
 
-	    if (elOld->getIndex() == 15376 && ichild == 0) {
-	      MSG("TEST k = %d, nbk = %d, el = %d\n", k, 
-		  nbk->getDof(0,0), elOld->getDof(ichild,0));
-	      MSG("TEST PTR nbk = %p, el = %p\n",
-		  nbk->getDof(0), elOld->getDof(ichild));
-	    }
-
 	    if (nbk->getDof(0) == elOld->getDof(ichild)) {
-	      if (elOld->getIndex() == 15376 && ichild == 0) {
-		MSG("DRIN!\n");
-	      }
-
 	      /* opp. vertex */
 	      dof = const_cast<int*>(nb->getDof(elInfoOld->oppVertex[cv[i]])); 
 	      
@@ -641,17 +615,13 @@ namespace AMDiS {
 		  }
 		  (*neigh_local)[i] = nbk->getChild(0);
 		  oppVertex[i] = 3;
-		  if (elOld->getIndex() == 15376 && ichild == 0) {
-		    MSG("BREAK 1\n");
-		  }
+
 		  break;
 		}
 	      } else {
 		if (dof != nbk->getDof(2)) { 
 		  ov = -1; 
-		  if (elOld->getIndex() == 15376 && ichild == 0) {
-		    MSG("BREAK 2\n");
-		  }
+
 		  break; 
 		}
 		ov = 2;
@@ -663,18 +633,12 @@ namespace AMDiS {
 
 	      (*neigh_local)[i] = nbk;
 	      oppVertex[i] = ov;
-	      if (elOld->getIndex() == 15376 && ichild == 0) {
-		MSG("BREAK 3\n");
-	      }
+
 	      break;
 	    }
 	    
 	  } /* end for k */
 
-	  if (elOld->getIndex() == 15376 && ichild == 0) {
-	    MSG("k = %d  ov = %d\n", k, ov);
-	  }
-
 
 	  // === Test if periodic. ===
 

AMDiS/src/FiniteElemSpace.cc

@@ -89,7 +89,7 @@ namespace AMDiS {
   }
 
 
-  FiniteElemSpace *FiniteElemSpace::provideFeSpace(DOFAdmin *admin,
+  FiniteElemSpace* FiniteElemSpace::provideFeSpace(DOFAdmin *admin,
 						   const BasisFunction *basFcts,
 						   Mesh *mesh,
 						   std::string name_)
@@ -104,6 +104,22 @@ namespace AMDiS {
   }
 
 
+#if DEBUG
+  FiniteElemSpace* FiniteElemSpace::provideFeSpace(Mesh *mesh)
+  {
+    FUNCNAME("FiniteElemSpace::provideFeSpace()");
+
+    for (unsigned int i = 0; i < feSpaces.size(); i++)
+      if (feSpaces[i]->mesh == mesh)
+	return feSpaces[i];
+
+    ERROR_EXIT("FE space not found!\n");
+    
+    return NULL;
+  }