Fixed small bug for mesh partitioning with periodic boundary conditions.

AMDiS/src/ElInfo3d.cc

@@ -501,18 +501,18 @@ namespace AMDiS {
     const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
     int *ce;                    /* ce = child_edge[el_type][ichild] */
     Element *nb, *nbk;
-    Element *el_old = elInfoOld->element;
+    Element *elOld = elInfoOld->element;
     Flag fillFlag_local = elInfoOld->fillFlag;
     DegreeOfFreedom *dof;
     int ov = -1;
     Mesh *mesh = elInfoOld->getMesh();
 
-    TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); 
+    TEST_EXIT_DBG(elOld->getChild(0))("missing child?\n"); 
 
-    element = const_cast<Element*>(el_old->getChild(ichild));
+    element = const_cast<Element*>(elOld->getChild(ichild));
     macroElement = elInfoOld->macroElement;
     fillFlag = fillFlag_local;
-    parent  = el_old;
+    parent = elOld;
     level = elInfoOld->level + 1;
     iChild = ichild;
     int el_type_local = 0;
@@ -539,8 +539,8 @@ namespace AMDiS {
 	for (int j = 0; j < dimOfWorld; j++)
 	  coord[i][j] = elInfoOld->coord[cv[i]][j];
       }
-      if (el_old->getNewCoord()) {
-	coord[3] = *(el_old->getNewCoord());
+      if (elOld->getNewCoord()) {
+	coord[3] = *(elOld->getNewCoord());
       } else {
 	for (int j = 0; j < dimOfWorld; j++)
 	  coord[3][j] = (elInfoOld->coord[0][j] + elInfoOld->coord[1][j]) / 2;
@@ -556,11 +556,10 @@ namespace AMDiS {
       Flag fill_opp_coords;
       fill_opp_coords.setFlags(fillFlag_local & Mesh::FILL_OPP_COORDS);
       
-      /*----- nb[0] is other child --------------------------------------------*/
+      // === nb[0] is other child ===
       
-      /*    if (nb = el_old->child[ochild])   old version  */
-      if (el_old->getChild(0)  &&  
-	  (nb = const_cast<Element*>(el_old->getChild(ochild)))) {
+      if (elOld->getChild(0) &&  
+	  (nb = const_cast<Element*>(elOld->getChild(ochild)))) {
 	
 	if (nb->getChild(0)) {         /* go down one level for direct neighbour */
 	  if (fill_opp_coords.isAnySet()) {
@@ -588,16 +587,43 @@ namespace AMDiS {
       }
 
 
-      /*----- nb[1],nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
+      /*----- nb[1], nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
       
       for (int i = 1; i < 3; i++) {
 	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 (nbk->getDof(0) == el_old->getDof(ichild)) {
+
+	    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]])); 
 	      
@@ -615,11 +641,17 @@ 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;
@@ -631,18 +663,30 @@ namespace AMDiS {
 
 	      (*neigh_local)[i] = nbk;
 	      oppVertex[i] = ov;
+	      if (elOld->getIndex() == 15376 && ichild == 0) {
+		MSG("BREAK 3\n");
+	      }
 	      break;
 	    }
 	    
 	  } /* end for k */
 
-	  // periodic ?
+	  if (elOld->getIndex() == 15376 && ichild == 0) {
+	    MSG("k = %d  ov = %d\n", k, ov);
+	  }
+
+
+	  // === Test if periodic. ===
+
 	  if (k == 2 || ov == -1) {
-	    for (k = 0; k < 2; k++) {  /* look at both childs of old neighbour */
+
+	    // Look at both childs of old neighbour.
+	    for (k = 0; k < 2; k++) {
 	      nbk = const_cast<Element*>(nb->getChild(k));
-	      if (nbk->getDof(0) == el_old->getDof(ichild) ||
-		  mesh->associated(nbk->getDof(0, 0), el_old->getDof(ichild, 0))) {
-		/* opp. vertex */
+	      if (nbk->getDof(0) == elOld->getDof(ichild) ||
+		  mesh->associated(nbk->getDof(0, 0), elOld->getDof(ichild, 0))) {
+
+		// opp. vertex 
 		dof = const_cast<int*>(nb->getDof(elInfoOld->oppVertex[cv[i]])); 
 		
 		if (dof == nbk->getDof(1) || 
@@ -680,7 +724,9 @@ namespace AMDiS {
 	      
 	    } /* end for k */
 
-	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
+	    TEST_EXIT_DBG(k < 2)
+	      ("Fill %d child of element %d (0-Level: %d): Child not found with vertex!\n",
+	       ichild, elOld->getIndex(), elInfoOld->getMacroElement()->getIndex());
 	  }
 	} else {
 	  (*neigh_local)[i] = NULL;

AMDiS/src/parallel/MeshDistributor.cc

@@ -204,11 +204,15 @@ namespace AMDiS {
       for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
 	if ((*it)->getNeighbour(i) && 
 	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
+
+	  int neighIndex = (*it)->getNeighbour(i)->getIndex();
+
 	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
 	  (*it)->setNeighbour(i, NULL);
 	  (*it)->setBoundary(i, 0);
 
 	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
+	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
 	}
       }
     }
@@ -218,11 +222,15 @@ namespace AMDiS {
       for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
 	if ((*it)->getNeighbour(i) && 
 	    mesh->isPeriodicAssociation((*it)->getBoundary(i))) {
+
+	  int neighIndex = (*it)->getNeighbour(i)->getIndex();
+
 	  (*it)->getNeighbour(i)->setNeighbour((*it)->getOppVertex(i), NULL);
 	  (*it)->setNeighbour(i, NULL);
 	  (*it)->setBoundary(i, 0);
 
 	  macroElementNeighbours[(*it)->getIndex()][i] = -1;
+	  macroElementNeighbours[neighIndex][(*it)->getOppVertex(i)] = -1;
 	}
       }
     }
@@ -974,7 +982,7 @@ namespace AMDiS {
     }
 
 
-    // === Run ParMETiS to calculate a new mesh partitioning.  ===
+    // === Run mesh partitioner to calculate a new mesh partitioning.  ===
 
     partitioner->setLocalGlobalDofMap(&mapLocalGlobalDofs);
     bool partitioningSucceed = partitioner->partition(elemWeights, ADAPTIVE_REPART);
@@ -983,6 +991,7 @@ namespace AMDiS {
       return;
     }
 
+
     // === Create map that maps macro element indices to pointers to the ===
     // === macro elements.                                               ===
 
@@ -1060,6 +1069,7 @@ namespace AMDiS {
       }
     }
 
+
     StdMpi<MeshCodeVec> stdMpi(mpiComm, true);
     stdMpi.send(sendCodes);
     for (std::map<int, std::vector<int> >::iterator it = partitioner->getRecvElements().begin();
@@ -1067,6 +1077,7 @@ namespace AMDiS {
       stdMpi.recv(it->first);
     stdMpi.startCommunication();
 
+
     if (interchangeVectors.size() == 0)
       WARNING("There are no interchange vectors defined!\n");
 

AMDiS/src/parallel/MeshManipulation.cc

@@ -141,7 +141,7 @@ namespace AMDiS {
 	    el1->getNonVertexDofs(feSpace, b1, dofs1);
 
 #if (DEBUG != 0)
-	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
+	    //	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
 #endif
 
  	    for (unsigned int i = 0; i < dofs0.size(); i++)
@@ -183,7 +183,7 @@ namespace AMDiS {
 	    el1->getNonVertexDofs(feSpace, b1, dofs1);
 
 #if (DEBUG != 0)
-	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
+	    //	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
 #endif
 
  	    for (unsigned int i = 0; i < dofs0.size(); i++)