fix bug

3 files changed, 204 insertions, 35 deletions

diff --git a/header.hpp b/header.hpp
index 8119137..2a9aed5 100644
--- a/header.hpp
+++ b/header.hpp
@@ -19,6 +19,7 @@
 #define N_EDGES 3
 #define SURFACE_N_EDGES 3
 #define SUBSURFACE_N_EDGES 9
+#define SUBSURFACE_N_NODES 6
 #define DIM 3
 
 #define N_AGENT 300
diff --git a/makefile b/makefile
index 79666b8..1c13939 100644
--- a/makefile
+++ b/makefile
@@ -1,6 +1,6 @@
 PRG=agent
 
-CC=g++
+CC=g++-10
 CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
 LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
 SRC=mesh.cpp main.cpp tools.cpp agent.cpp
diff --git a/mesh.cpp b/mesh.cpp
index 3c3e1e2..125c4de 100644
--- a/mesh.cpp
+++ b/mesh.cpp
@@ -309,6 +309,8 @@ int Mesh::build_subsurface(H5::H5File fp)
       int n5 = md_e[i][5];
       int n6 = md_e[i][6];
 
+      //      std::cout << n1 << " " << n2 << " " << n3 << " " <<  n4 << " " << n5 << " " << n6 << std::endl;
+
       for (size_t j = 0; j < dims_n[1]; j++)
 	{
 	  coords1[j] = md_n[n1][j];
@@ -365,7 +367,7 @@ int Mesh::build_subsurface(H5::H5File fp)
 
   printf("[OK] initialize edges {\n");
   
-  int ns[SUBSURFACE_N_EDGES];
+  int ns[SUBSURFACE_N_NODES];
   
   //  hsize_t **szudzik_ids = new hsize_t*[ny_e];
   //  std::vector<std::vector<hsize_t>> edge_ids(ny_e, SUBSURFACE_N_EDGES, 2);
@@ -378,29 +380,105 @@ int Mesh::build_subsurface(H5::H5File fp)
   for (size_t i = 0; i < ny_e; i ++)
     {
       
-      for (size_t j = 0; j < SUBSURFACE_N_EDGES; j ++)
+      for (size_t j = 0; j < SUBSURFACE_N_NODES; j ++) {
 	ns[j] = md_e[i][j+1];
-	
-      for (size_t j = 0; j < SUBSURFACE_N_EDGES - 1; j ++) {
-
-	if (ns[j] < ns[j+1]) {
-	  edge_ids[i][j][0] = ns[j];
-	  edge_ids[i][j][1] = ns[j+1];
-	} else {
-	  edge_ids[i][j][0] = ns[j+1];
-	  edge_ids[i][j][1] = ns[j];
-	}
-	//	szudzik_ids[i][j] = szudzik(ns[j], ns[j+1]);
-	
       }
 
-      if (ns[SUBSURFACE_N_EDGES-1] < ns[0]) {
-	edge_ids[i][SUBSURFACE_N_EDGES-1][0] = ns[SUBSURFACE_N_EDGES-1];
-	edge_ids[i][SUBSURFACE_N_EDGES-1][1] = ns[0];
+      // top of the cell
+
+      if (ns[0] < ns[1]) {
+	edge_ids[i][0][0] = ns[0];
+	edge_ids[i][0][1] = ns[1];
+      } else {
+	edge_ids[i][0][1] = ns[0];
+	edge_ids[i][0][0] = ns[1];
+      }
+
+      if (ns[1] < ns[2]) {
+	edge_ids[i][1][0] = ns[1];
+	edge_ids[i][1][1] = ns[2];
+      } else {
+	edge_ids[i][1][1] = ns[2];
+	edge_ids[i][1][0] = ns[1];
+      }
+
+      if (ns[2] < ns[0]) {
+	edge_ids[i][2][0] = ns[2];
+	edge_ids[i][2][1] = ns[0];
+      } else {
+	edge_ids[i][2][1] = ns[0];
+	edge_ids[i][2][0] = ns[2];
+      }
+
+      // bottom of the cell
+
+      if (ns[3] < ns[4]) {
+	edge_ids[i][3][0] = ns[3];
+	edge_ids[i][3][1] = ns[4];
+      } else {
+	edge_ids[i][3][0] = ns[4];
+	edge_ids[i][3][1] = ns[3];
+      }
+
+      if (ns[4] < ns[5]) {
+	edge_ids[i][4][0] = ns[4];
+	edge_ids[i][4][1] = ns[5];
+      } else {
+	edge_ids[i][4][0] = ns[5];
+	edge_ids[i][4][1] = ns[4];
+      }
+
+      if (ns[5] < ns[3]) {
+	edge_ids[i][5][0] = ns[5];
+	edge_ids[i][5][1] = ns[3];
+      } else {
+	edge_ids[i][5][0] = ns[3];
+	edge_ids[i][5][1] = ns[5];
+      }
+
+      // sides of the cell
+
+      if (ns[0] < ns[3]) {
+	edge_ids[i][6][0] = ns[0];
+	edge_ids[i][6][1] = ns[3];
+      } else {
+	edge_ids[i][6][0] = ns[3];
+	edge_ids[i][6][1] = ns[0];
+      }
+
+      if (ns[1] < ns[4]) {
+	edge_ids[i][7][0] = ns[1];
+	edge_ids[i][7][1] = ns[4];
       } else {
-	edge_ids[i][SUBSURFACE_N_EDGES-1][0] = ns[0];
-	edge_ids[i][SUBSURFACE_N_EDGES-1][1] = ns[SUBSURFACE_N_EDGES-1];
+	edge_ids[i][7][0] = ns[4];
+	edge_ids[i][7][1] = ns[1];
       }
+
+      if (ns[2] < ns[5]) {
+	edge_ids[i][8][0] = ns[2];
+	edge_ids[i][8][1] = ns[5];
+      } else {
+	edge_ids[i][8][0] = ns[5];
+	edge_ids[i][8][1] = ns[2];
+      }
+	
+      // for (size_t j = 0; j < SUBSURFACE_N_EDGES - 1; j ++) {
+      // 	if (ns[j] < ns[j+1]) {
+      // 	  edge_ids[i][j][0] = ns[j];
+      // 	  edge_ids[i][j][1] = ns[j+1];
+      // 	} else {
+      // 	  edge_ids[i][j][0] = ns[j+1];
+      // 	  edge_ids[i][j][1] = ns[j];
+      // 	}
+      // 	//	szudzik_ids[i][j] = szudzik(ns[j], ns[j+1]);
+      // }
+      // if (ns[SUBSURFACE_N_EDGES-1] < ns[0]) {
+      // 	edge_ids[i][SUBSURFACE_N_EDGES-1][0] = ns[SUBSURFACE_N_EDGES-1];
+      // 	edge_ids[i][SUBSURFACE_N_EDGES-1][1] = ns[0];
+      // } else {
+      // 	edge_ids[i][SUBSURFACE_N_EDGES-1][0] = ns[0];
+      // 	edge_ids[i][SUBSURFACE_N_EDGES-1][1] = ns[SUBSURFACE_N_EDGES-1];
+      // }
       //     szudzik_ids[i][SUBSURFACE_N_EDGES-1] = szudzik(ns[SUBSURFACE_N_EDGES - 1], ns[0]);
       
     }
@@ -418,6 +496,16 @@ int Mesh::build_subsurface(H5::H5File fp)
 	}
       printf("\n");
     }
+
+//    for (size_t i = 0; i < maxprint; i ++)
+//    {
+//      printf("[ok]");
+//      for (size_t j = 0; j < SUBSURFACE_N_EDGES; j ++)
+//	{
+//	  printf(" {%lld %f %f %f} \n",  edge_ids[i][j][0], md_n[edge_ids[i][j][0]][0], md_n[edge_ids[i][j][0]][1], md_n[edge_ids[i][j][0]][2]);
+//	}
+//      printf("\n");
+//    }
   
   #endif
   
@@ -430,50 +518,130 @@ int Mesh::build_subsurface(H5::H5File fp)
   hsize_t e10;
   hsize_t e11;
   
+  hsize_t e20;
+  hsize_t e21;
+  hsize_t e30;
+  hsize_t e31;
+  
 #pragma omp parallel for
   for (size_t i = 0; i < ny_e; i ++)
     {
-
-      printf("[ok] searching neighbors of cell %ld/%ld\n", i+1, ny_e);
+#if VERBOSE
+      if (i % 100 == 0)
+	printf("[ok] searching neighbors of cell %ld/%ld\n", i+1, ny_e);
+#endif
       
       for (size_t j = 0; j < SUBSURFACE_N_EDGES; j ++)
 	{
 
+	  int flag = -1;
+
+#if DEBUG
 	  printf("[ok] searching for edge %ld/%d\n", j+1, SUBSURFACE_N_EDGES);
-	  
-	  //	  e0 = szudzik_ids[i][j];
+#endif
 	  
 	  e00 = edge_ids[i][j][0];
 	  e01 = edge_ids[i][j][1];
+
+#if DEBUG
+	  printf(">> %ld:%ld [%lld %lld]\n", i, j, e00, e01);
+#endif
 	  
 	  cells[i].set_neighbor(j, -1);
 
 	  for (size_t k = 0; k < ny_e; k ++)
 	    {
 
-	      int counter = 0;
+	      int valid_neighbor = 1;
 	      
-	      for (size_t m = 0; m < SUBSURFACE_N_EDGES; m ++)
-		{
+	      for (size_t jj = 0; jj < j; jj ++) {
+		if (cells[i].get_neighbor(jj) == k) {
+		  // already a neighbor?
+#if DEBUG
+		  printf(">> >> >> cell already stored as a neighbor: %ld \n", k);
+#endif
+		  valid_neighbor = -1;
+		}
+	      }
+
+	      if (valid_neighbor > 0) {
+	      
+		for (size_t m = 0; m < SUBSURFACE_N_EDGES; m ++) {
 		  
-		  //		  e1 = szudzik_ids[k][m];
 		  e10 = edge_ids[k][m][0];
 		  e11 = edge_ids[k][m][1];
 
-		  //		  printf("[%lld %lld] == [%lld %lld]?\n", e00, e01, e10, e11);
+#if DEBUG
+		  printf(">> >> %ld:%ld [%lld %lld]\n", k, m, e10, e11);
+#endif
+
+		  if (e00 == e10 && e01 == e11 && i != k) { // found one match
+
+#if DEBUG
+		    printf(">> >> >> [ok] found one match for edge %ld of cell %ld\n", m, k);
+#endif
+
+		    for (size_t r = 0; r < SUBSURFACE_N_EDGES; r ++) {
+
+		      if (r != j) {
+			
+			e20 = edge_ids[i][r][0];
+			e21 = edge_ids[i][r][1];
+
+#if DEBUG
+			printf(">> >> >> >> %ld:%ld [%lld %lld]\n", i, j, e20, e21);
+#endif
+
+			for (size_t p = 0; p < SUBSURFACE_N_EDGES; p ++) {
+			  
+			  e30 = edge_ids[k][p][0];
+			  e31 = edge_ids[k][p][1];
+
+#if DEBUG
+			  printf(">> >> >> >> >> %ld:%ld [%lld %lld]\n", i, j, e30, e31);
+#endif
+
+			  if (e20 == e30 && e21 == e31) {
+#if DEBUG
+			    printf(">> >> >> >> >> [ok] found second match for edge %ld of cell %ld\n", p, k);
+#endif
+#if VERBOSE
+			    if (i % 100 == 0)
+			      printf("[ok] found neighbor %ld\n", k);
+#endif
+			    flag = 1;
+			    cells[i].set_neighbor(j, (hsize_t) k);
+			    break;
+			  }
+			  
+			}
+			
+		      }
+
+		      if (flag > 0) {
+			break;
+		      }
+
+		    }
+
+		    if (flag > 0) {
+		      break;
+		    }
 
-		  if (e00 == e10 && e01 == e11 && i != k) {
-		    //		    printf("[ok] found one match for edge %ld of cell %ld\n", m, k);
-		    counter = counter + 1;
 		  }
 
-		  if (counter > 2) {
-		    //		    printf("[ok] found neighbor of cell %ld/%ld for edge %ld\n", i+1, ny_e, j+1);
-		    cells[i].set_neighbor(j, (hsize_t) k);
+		  if (flag > 0) {
 		    break;
 		  }
 		  
 		}
+
+	      }
+
+	      if (flag > 0) {
+		break;
+	      }
+			  
 	    }
 	}
     }