Fixed coarsening bug with periodic BCs in sequentiell code.

AMDiS/src/CoarseningManager3d.cc

@@ -313,46 +313,70 @@ namespace AMDiS {
     coarsenList.setElement(*n_neigh, neigh);
     coarsenList.setElType(*n_neigh, neigh_info->getType()); 
   
-    int n_vertices = mesh->getGeo(VERTEX);
-    int i, j, k, edge_no;
-			 
+    int nVertices = mesh->getGeo(VERTEX);			 
+
     while (neigh != el) {
-      for (j = 0; j < n_vertices; j++)
+      // === Find the coarsening edge on the current neighbour element. ===
+
+      int i, j, k;
+
+      // === First, try to identify the edge DOFs directly. ===
+
+      for (j = 0; j < nVertices; j++)
 	if (neigh->getDof(j) == edge[0])
 	  break;
-      for (k = 0; k < n_vertices; k++)
+
+      for (k = 0; k < nVertices; k++)
 	if (neigh->getDof(k) == edge[1])
 	  break;
 
-      if (j > 3 || k > 3) {
-	for (j = 0; j < n_vertices; j++)
+      // === If one of the edge DOFs was not found, try to make use of periodic ===
+      // === DOF associations. First, use the direct associations between DOFs. ===
+      // === If this will not work, continue with testing on indirect           ===
+      // === associations.                                                      ===
+
+      if (j >= nVertices) {
+	for (j = 0; j < nVertices; j++)
 	  if (mesh->associated(neigh->getDof(j, 0), edge[0][0]))
 	    break;
-	for (k = 0; k < n_vertices; k++)
+	
+	if (j >= nVertices)
+	  for (j = 0; j < nVertices; j++)
+	    if (mesh->indirectlyAssociated(neigh->getDof(j, 0), edge[0][0]))
+	      break;
+
+	TEST_EXIT_DBG(j < nVertices)
+	  ("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
+	   el->getIndex(), edge[0][0], neigh->getIndex(), neigh->getDof(0, 0),
+	   neigh->getDof(1, 0), neigh->getDof(2, 0), neigh->getDof(3, 0));
+      }
+      
+      if (k >= nVertices) {
+	for (k = 0; k < nVertices; k++)
 	  if (mesh->associated(neigh->getDof(k, 0), edge[1][0]))
 	    break;
+
+	if (k >= nVertices)
+	  for (k = 0; k < nVertices; k++)
+	    if (mesh->indirectlyAssociated(neigh->getDof(k, 0), edge[1][0]))
+	      break;
      
-	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),
+	TEST_EXIT_DBG(k < nVertices)
+	  ("Process element %d: DOF %d not found on element %d with nodes (%d %d %d %d)\n", 
+	   el->getIndex(), 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);
 
-      /****************************************************************************/
-      /*  get the direction of the next neighbour				    */
-      /****************************************************************************/
+      int edgeNo = Tetrahedron::edgeOfDofs[j][k];
+      coarsenList.setCoarsePatch(*n_neigh, edgeNo == 0);
+
+
+      // === Get the direction of the next neighbour. ===
 
-      if (next_el[edge_no][0] != opp_v)
-	i = next_el[edge_no][0];
+      if (next_el[edgeNo][0] != opp_v)
+	i = next_el[edgeNo][0];
       else
-	i = next_el[edge_no][1];
+	i = next_el[edgeNo][1];
 
       ++*n_neigh;
 
@@ -383,7 +407,7 @@ namespace AMDiS {
     /*  the domain's boundary is reached; loop back to the starting el          */
     /****************************************************************************/
 
-    i = *n_neigh - 1;
+    int i = *n_neigh - 1;
     opp_v = coarsenList.getOppVertex(i, 0);
     do {
       TEST_EXIT_DBG(neigh_info->getNeighbour(opp_v)  &&  i > 0)