fix bug

3 files changed, 51 insertions, 98 deletions

diff --git a/header.hpp b/header.hpp
index fffd4ae..213d210 100644
--- a/header.hpp
+++ b/header.hpp
@@ -16,14 +16,16 @@
 #include <vector>
 
 #define VERBOSE 1
+#define VERBOSE_HIGH 1
+
 #define N_EDGES 3
 #define SURFACE_N_EDGES 3
 #define SUBSURFACE_N_EDGES 9
 #define SUBSURFACE_N_NODES 6
-#define SUBSURFACE_N_NEIGHBORS 5
+#define SUBSURFACE_N_NEIGHBORS 5 // 5 von Neumann neighbors, 11 if Moore neighborhood
 #define DIM 3
 
-#define RADIUS 1000.0
+#define RADIUS 200.0
 
 #define N_AGENT 300
 #define N_ITERATION 5000
diff --git a/main.cpp b/main.cpp
index efd2912..e16886f 100644
--- a/main.cpp
+++ b/main.cpp
@@ -27,6 +27,8 @@ int main(void)
   Mesh subsurface_mesh = Mesh();
   subsurface_mesh.build_subsurface(subsurface_fp);
 
+  printf("[OK] # of cells: %d\n\n", subsurface_mesh.get_n_cells());
+
   subsurface_fp.close();
 #endif
   
diff --git a/mesh.cpp b/mesh.cpp
index c71a6a8..16c103a 100644
--- a/mesh.cpp
+++ b/mesh.cpp
@@ -479,71 +479,58 @@ int Mesh::build_subsurface(H5::H5File fp)
   hsize_t e71;
   hsize_t e80;
   hsize_t e81;
-  hsize_t e90;
-  hsize_t e91;
-
+  
   hsize_t _e0;
   hsize_t _e1;
 
-  
 #pragma omp parallel for
   for (size_t i = 0; i < ny_e; i ++) {
-
     
     double xi = cells[i].get_x();
     double yi = cells[i].get_y();
     double zi = cells[i].get_z();
+
+    for (size_t j = 0; j < SUBSURFACE_N_NEIGHBORS; j ++)
+      cells[i].set_neighbor(j, -1);
       
-#if VERBOSE
+#if VERBOSE_HIGH
       if (i % 100 == 0)
-	printf("[ok] searching neighbors of cell %ld/%ld\n", i+1, ny_e);
+	printf("[ok] searching neighbors of cell: %ld/%ld\n", i+1, ny_e);
 #endif
 
-      // neighbor #1
-      
       e00 = edge_ids[i][0][0];
       e01 = edge_ids[i][0][1];
 
       e10 = edge_ids[i][1][0];
       e11 = edge_ids[i][1][1];
 
-      // neighbor #2
-
-      e20 = edge_ids[i][3][0];
-      e21 = edge_ids[i][3][1];
-
-      e30 = edge_ids[i][4][0];
-      e31 = edge_ids[i][4][1];
+      e20 = edge_ids[i][2][0];
+      e21 = edge_ids[i][2][1];
 
-      // neighbor #3
+      e30 = edge_ids[i][3][0];
+      e31 = edge_ids[i][3][1];
       
-      e40 = edge_ids[i][6][0];
-      e41 = edge_ids[i][6][1];
+      e40 = edge_ids[i][4][0];
+      e41 = edge_ids[i][4][1];
 
-      e50 = edge_ids[i][7][0];
-      e51 = edge_ids[i][7][1];
+      e50 = edge_ids[i][5][0];
+      e51 = edge_ids[i][5][1];
 
-      // neighbor #4
+      e60 = edge_ids[i][6][0];
+      e61 = edge_ids[i][6][1];
 
-      e60 = edge_ids[i][7][0];
-      e61 = edge_ids[i][7][1];
-
-      e70 = edge_ids[i][8][0];
-      e71 = edge_ids[i][8][1];
-
-      // neighbor #5
+      e70 = edge_ids[i][7][0];
+      e71 = edge_ids[i][7][1];
 
       e80 = edge_ids[i][8][0];
       e81 = edge_ids[i][8][1];
 
-      e90 = edge_ids[i][6][0];
-      e91 = edge_ids[i][6][1];
-
-      printf("%ld [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld]\n", i, e00, e01, e10, e11, e20, e21, e30, e31, e40, e41, e50, e51, e60, e61, e70, e71, e80, e81, e90, e91);
-
-      // printf("{ %lld, %lld, %lld, %lld, %lld, %lld, %lld, %lld, %lld }\n", e00, e10, e20, e30, e40, e50, e60, e70, e80, e90);
-      // printf("{ %lld, %lld, %lld, %lld, %lld, %lld, %lld, %lld, %lld }\n", e01, e11, e21, e31, e41, e51, e61, e71, e81, e91);
-
+#if VERBOSE_HIGH
+      if (i % 100 == 0)
+	printf("[ok] indices: [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld], [%lld, %lld]\n", e00, e01, e10, e11, e20, e21, e30, e31, e40, e41, e50, e51, e60, e61, e70, e71, e80, e81);
+#endif
+      
+      int counter = 0;
       
       for (size_t k = 0; k < ny_e; k ++) {
 
@@ -554,82 +541,44 @@ int Mesh::build_subsurface(H5::H5File fp)
 
 	double distance = sqrt((xk - xi) * (xk - xi) + (yk - yi) * (yk - yi) + (zk - zi) * (zk - zi));
 	
-	if (i != k) {
+	if (i != k && distance <= RADIUS) {
 
-	  // a neighbor is found if two edges of the same face match
-	  int counter1 = 0;
-	  int counter2 = 0;
-	  int counter3 = 0;
-	  int counter4 = 0;
-	  int counter5 = 0;
-
-#define VERBOSE_HIGH 1
+	  int _tmp_count = 0;
 	  
 	  for (size_t m = 0; m < SUBSURFACE_N_EDGES; m ++) {
 
 	    _e0 = edge_ids[k][m][0];
 	    _e1 = edge_ids[k][m][1];
 
-	    printf("%ld-%ld { %lld, %lld }\n", k, m, _e0, _e1);
+	    if ((e00 == _e0 && e01 == _e1) || (e10 == _e0 && e11 == _e1) ||
+		(e20 == _e0 && e21 == _e1) || (e30 == _e0 && e31 == _e1) ||
+		(e40 == _e0 && e41 == _e1) || (e50 == _e0 && e51 == _e1) ||
+		(e60 == _e0 && e61 == _e1) || (e70 == _e0 && e71 == _e1) ||
+		(e80 == _e0 && e81 == _e1)) {
 
-	    if ((e00 == _e0 && e01 == _e1) || (e10 == _e0 && e11 == _e1)) {
-	      counter1 ++;
-	      if (counter1 == 2) {
-		cells[i].set_neighbor(0, (hsize_t) k);
-#if VERBOSE_HIGH
-		if (i % 100 == 0)
-		  printf("[ok] found neighbor %ld\n", k+1);
-#endif
-		break;
-	      }
-	    } else if ((e20 == _e0 && e21 == _e1) || (e30 == _e0 && e31 == _e1)) {
-	      counter2 ++;
-	      if (counter2 == 2) {
-		cells[i].set_neighbor(1, (hsize_t) k);
-#if VERBOSE_HIGH
-		if (i % 100 == 0)
-		  printf("[ok] found neighbor %ld\n", k+1);
-#endif
-		break;
-	      }
-	    } else if ((e40 == _e0 && e41 == _e1) || (e50 == _e0 && e51 == _e1)) {
-	      counter3 ++;
-	      if (counter3 == 2) {
-		cells[i].set_neighbor(2, (hsize_t) k);
-#if VERBOSE_HIGH
-		if (i % 100 == 0)
-		  printf("[ok] found neighbor %ld\n", k+1);
-#endif
-		break;
-	      }
-	    } else if ((e60 == _e0 && e61 == _e1) || (e70 == _e0 && e71 == _e1)) {
-	      counter4 ++;
-	      if (counter4 == 2) {
-		cells[i].set_neighbor(3, (hsize_t) k);
-#if VERBOSE_HIGH
-		if (i % 100 == 0)
-		  printf("[ok] found neighbor %ld\n", k+1);
-#endif
-		break;
-	      }
-	    } else if ((e80 == _e0 && e81 == _e1) || (e90 == _e0 && e91 == _e1)) {
-	      counter5 ++;
-	      if (counter5 == 2) {
-		cells[i].set_neighbor(4, (hsize_t) k);
+	      _tmp_count ++; 
+
+	      if (_tmp_count == 2) { // a cell k is a neighbor of cell
+				     // i if 2 of their edges match,
+				     // which means that a surface is
+				     // shared between k and i
+		cells[i].set_neighbor(counter, (hsize_t) k);
+		counter ++;
 #if VERBOSE_HIGH
 		if (i % 100 == 0)
-		  printf("[ok] found neighbor %ld\n", k+1);
+		  printf("[ok] found neighbor %d/%d: %ld [%lld, %lld]\n", counter, SUBSURFACE_N_NEIGHBORS, k, _e0, _e1);
 #endif
 		break;
 	      }
+	      
 	    }
-	    
 	  }
-	  
 	}
+
+	if (counter == SUBSURFACE_N_NEIGHBORS)
+	  break;
 	
       }
-
   }
 
   #if VERBOSE
@@ -639,7 +588,7 @@ int Mesh::build_subsurface(H5::H5File fp)
   for (size_t i = 0; i < maxprint; i ++)
     {
       printf("[ok]");
-      for (size_t j = 0; j < N_EDGES; j ++)
+      for (size_t j = 0; j < SUBSURFACE_N_NEIGHBORS; j ++)
 	{
 	  printf(" %6lld ", cells[i].get_neighbor(j));
 	}