add subsurface domain - not working

6 files changed, 429 insertions, 34 deletions

diff --git a/agent.cpp b/agent.cpp
index 11fdcfe..fbb0f22 100644
--- a/agent.cpp
+++ b/agent.cpp
@@ -55,7 +55,7 @@ int HydroAgent::advance_surface(Mesh *mesh, size_t step)
 
   double vx = location->get_v(0);
   double vy = location->get_v(1);
-  
+
   double x_agent = x[step-1] + DT_SURFACE * vx;
   double y_agent = y[step-1] + DT_SURFACE * vy;
 
diff --git a/header.hpp b/header.hpp
index a4622ff..8119137 100644
--- a/header.hpp
+++ b/header.hpp
@@ -15,20 +15,27 @@
 #include <math.h>
 #include <vector>
 
-#define VERBOSE 0
+#define VERBOSE 1
 #define N_EDGES 3
+#define SURFACE_N_EDGES 3
+#define SUBSURFACE_N_EDGES 9
 #define DIM 3
 
-#define N_AGENT 30
-#define N_ITERATION 200
-#define DT_SURFACE 1.0
-#define DT_SUBSURFACE 3600.0
+#define N_AGENT 300
+#define N_ITERATION 5000
+#define CFL 0.5
+#define DT_SURFACE 1800.0
+#define DT_SUBSURFACE 1800.0
+
+#define SUBSURFACE 1
 
 #define SURFACE_FLAG 1
 #define SUBSURFACE_FLAG 2
 #define OUTOFBOUNDS_FLAG 3
 
-#define MESH_FILE_NAME "data/visdump_surface_mesh.h5"
-#define DATA_FILE_NAME "data/visdump_surface_data.h5"
+#define MESH_FILE_NAME "test/visdump_surface_mesh.h5"
+#define DATA_FILE_NAME "test/visdump_surface_data.h5"
+#define SUBSURFACE_MESH_FILE_NAME "test/visdump_mesh.h5"
+#define SUBSURFACE_DATA_FILE_NAME "test/visdump_data.h5"
 
 #endif
diff --git a/main.cpp b/main.cpp
index 7628088..efd2912 100644
--- a/main.cpp
+++ b/main.cpp
@@ -15,11 +15,20 @@ int main(void)
   H5::H5File fp = H5::H5File(MESH_FILE_NAME, H5F_ACC_RDONLY);
   
   Mesh mesh = Mesh();
-  mesh.build(fp);
+  mesh.build_surface(fp);
   
   printf("[OK] # of cells: %d\n\n", mesh.get_n_cells());
   
   fp.close();
+
+#if SUBSURFACE
+  H5::H5File subsurface_fp = H5::H5File(SUBSURFACE_MESH_FILE_NAME, H5F_ACC_RDONLY);
+
+  Mesh subsurface_mesh = Mesh();
+  subsurface_mesh.build_subsurface(subsurface_fp);
+
+  subsurface_fp.close();
+#endif
   
   H5::H5File dp = H5::H5File(DATA_FILE_NAME, H5F_ACC_RDONLY);
   
@@ -39,13 +48,21 @@ int main(void)
   printf("[OK] inject agents\n");
   
   for (size_t i = 0; i < N_AGENT; i ++) {
+    
     index = rand() % len;
-    agents[i].set_location(&(mesh.get_cells()[index]));
-    agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 0);
-    agents[i].set_domain(SURFACE_FLAG, 0);
-    // bootstrap
-    agents[i].set_domain(SURFACE_FLAG, 1);
-    agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 1);
+
+    if (fabs(mesh.get_cells()[index].get_v(0)) > 0.0 ||
+	fabs(mesh.get_cells()[index].get_v(1)) > 0.0) {
+        
+      agents[i].set_location(&(mesh.get_cells()[index]));
+      agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 0);
+      agents[i].set_domain(SURFACE_FLAG, 0);
+      // bootstrap
+      agents[i].set_domain(SURFACE_FLAG, 1);
+      agents[i].set_path(agents[i].get_location()->get_x(), agents[i].get_location()->get_y(), 0.0, 1);
+      
+    }
+    
   }
 
   for (size_t step = 1; step < N_ITERATION; step ++) {
diff --git a/makefile b/makefile
index 153a637..79666b8 100644
--- a/makefile
+++ b/makefile
@@ -1,7 +1,7 @@
 PRG=agent
 
 CC=g++
-CXXFLAGS=-Wall -Wextra -pedantic -O0
+CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
 LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
 SRC=mesh.cpp main.cpp tools.cpp agent.cpp
 LIB=-lm -lhdf5 -lhdf5_cpp
diff --git a/mesh.cpp b/mesh.cpp
index 65da4c5..3c3e1e2 100644
--- a/mesh.cpp
+++ b/mesh.cpp
@@ -23,6 +23,13 @@ Cell::Cell(const Cell &cell)
   z = cell.z;
 }
 
+int Cell::init_neighbors(size_t n_neighbors)
+{
+  neighbors.resize(n_neighbors);
+  
+  return 0;
+}
+
 int Cell::set_coordinates(double* coordinates)
 {
   x = coordinates[0];
@@ -34,7 +41,7 @@ int Cell::set_coordinates(double* coordinates)
 
 int Cell::set_neighbors(hsize_t *value_array)
 {
-  for (int i = 0; i < N_EDGES; i ++)
+  for (size_t i = 0; i < neighbors.size(); i ++)
     neighbors[i] = value_array[i];
     
   return 0;
@@ -101,6 +108,7 @@ int Cell::get_cell_id()
 
 Mesh::Mesh()
 {
+  Mesh::n_cells = 0;
 }
 
 int Mesh::get_n_cells()
@@ -161,14 +169,14 @@ double Mesh::get_zmin()
   return zmin;
 }
 
-int Mesh::build(H5::H5File fp)
+int Mesh::build_subsurface(H5::H5File fp)
 {
 
-  printf("[OK] building the mesh\n");
+  printf("[OK] building the subsurface mesh\n");
   printf("[OK] reading element node map { \n");
     
-  H5::DataSet dataset_e     = fp.openDataSet("/0/Mesh/MixedElements");
-  H5::DataType datatype_e   = dataset_e.getDataType();
+  H5::DataSet dataset_e = fp.openDataSet("/0/Mesh/MixedElements");
+  H5::DataType datatype_e = dataset_e.getDataType();
   H5::DataSpace dataspace_e = dataset_e.getSpace();
   
   int rank = dataspace_e.getSimpleExtentNdims();
@@ -183,8 +191,9 @@ int Mesh::build(H5::H5File fp)
   
   int *data_e = new int[dims_e[0]*dims_e[1]];
 
-  size_t nx_e = 4;
-  size_t ny_e = dims_e[0] * dims_e[1] / 4;
+  size_t nx_e = 7;
+  size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
+
   int **md_e = new int*[ny_e];
   
   for (size_t i = 0; i < ny_e; i ++)
@@ -219,13 +228,14 @@ int Mesh::build(H5::H5File fp)
 
   printf("[OK] reading nodes { \n");
 
-  H5::DataSet dataset_n     = fp.openDataSet("/0/Mesh/Nodes");
-  H5::DataType datatype_n   = dataset_n.getDataType();
+  H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes");
+  H5::DataType datatype_n = dataset_n.getDataType();
   H5::DataSpace dataspace_n = dataset_n.getSpace();
 
   rank = dataspace_n.getSimpleExtentNdims();
 
   printf("[ok] rank of nodes: %d\n", rank);
+  
   hsize_t *dims_n = new hsize_t[rank];
   dataspace_n.getSimpleExtentDims(dims_n);
 
@@ -272,6 +282,9 @@ int Mesh::build(H5::H5File fp)
   double *coords1 = new double[dims_n[1]];
   double *coords2 = new double[dims_n[1]];
   double *coords3 = new double[dims_n[1]];
+  double *coords4 = new double[dims_n[1]];
+  double *coords5 = new double[dims_n[1]];
+  double *coords6 = new double[dims_n[1]];
 
   double *coords_avg = new double[dims_n[1]];
 
@@ -286,16 +299,26 @@ int Mesh::build(H5::H5File fp)
 
   for (size_t i = 0; i < ny_e; i ++)
     {
+
+      cells[i].init_neighbors(SUBSURFACE_N_EDGES);
+      
       int n1 = md_e[i][1];
       int n2 = md_e[i][2];
       int n3 = md_e[i][3];
+      int n4 = md_e[i][4];
+      int n5 = md_e[i][5];
+      int n6 = md_e[i][6];
 
       for (size_t j = 0; j < dims_n[1]; j++)
 	{
 	  coords1[j] = md_n[n1][j];
 	  coords2[j] = md_n[n2][j];
 	  coords3[j] = md_n[n3][j];
-	  coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
+	  coords4[j] = md_n[n4][j];
+	  coords5[j] = md_n[n5][j];
+	  coords6[j] = md_n[n6][j];
+
+	  coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0;
 	}
 
       for (size_t k = 0; k < DIM; k ++)
@@ -334,11 +357,345 @@ int Mesh::build(H5::H5File fp)
   delete[] coords1;
   delete[] coords2;
   delete[] coords3;
+  delete[] coords4;
+  delete[] coords5;
+  delete[] coords6;
   
   printf("[OK] } done\n");
 
   printf("[OK] initialize edges {\n");
+  
+  int ns[SUBSURFACE_N_EDGES];
+  
+  //  hsize_t **szudzik_ids = new hsize_t*[ny_e];
+  //  std::vector<std::vector<hsize_t>> edge_ids(ny_e, SUBSURFACE_N_EDGES, 2);
+
+  std::vector<std::vector<std::vector<hsize_t>>> edge_ids(ny_e, std::vector<std::vector<hsize_t>>(SUBSURFACE_N_EDGES, std::vector<hsize_t>(2)));
+  
+  //for (size_t i = 0; i < ny_e; i ++)
+  //  szudzik_ids[i] = new hsize_t[SUBSURFACE_N_EDGES];
+  
+  for (size_t i = 0; i < ny_e; i ++)
+    {
+      
+      for (size_t j = 0; j < SUBSURFACE_N_EDGES; 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];
+      } 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]);
+      
+    }
+
+  #if VERBOSE
+
+  maxprint = (ny_e < 10) ? ny_e : 10;
+    
+  for (size_t i = 0; i < maxprint; i ++)
+    {
+      printf("[ok]");
+      for (size_t j = 0; j < SUBSURFACE_N_EDGES; j ++)
+	{
+	  printf(" {%6lld %6lld} ", edge_ids[i][j][0], edge_ids[i][j][1]);
+	}
+      printf("\n");
+    }
+  
+  #endif
+  
+  printf("[OK] } done\n");
+
+  printf("[OK] initialize neighbors {\n");
+
+  hsize_t e00;
+  hsize_t e01;
+  hsize_t e10;
+  hsize_t e11;
+  
+#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);
+      
+      for (size_t j = 0; j < SUBSURFACE_N_EDGES; j ++)
+	{
+
+	  printf("[ok] searching for edge %ld/%d\n", j+1, SUBSURFACE_N_EDGES);
+	  
+	  //	  e0 = szudzik_ids[i][j];
+	  
+	  e00 = edge_ids[i][j][0];
+	  e01 = edge_ids[i][j][1];
+	  
+	  cells[i].set_neighbor(j, -1);
+
+	  for (size_t k = 0; k < ny_e; k ++)
+	    {
+
+	      int counter = 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 (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);
+		    break;
+		  }
+		  
+		}
+	    }
+	}
+    }
+
+  #if VERBOSE
+
+  maxprint = (ny_e < 10) ? ny_e : 10;
+  
+  for (size_t i = 0; i < maxprint; i ++)
+    {
+      printf("[ok]");
+      for (size_t j = 0; j < N_EDGES; j ++)
+	{
+	  printf(" %6lld ", cells[i].get_neighbor(j));
+	}
+      printf("\n");
+    }
+  
+  #endif
+  
+  printf("[OK] } done\n");
+
+  //  delete[] szudzik_ids;
+  delete[] md_e;
+  delete[] data_e;
+  delete[] dims_e;
+
+  delete[] md_n;
+  delete[] data_n;
+  delete[] dims_n;
+  
+  printf("[OK] mesh built\n\n");
+
+  return 0;
+  
+}
+
+int Mesh::build_surface(H5::H5File fp)
+{
+
+  printf("[OK] building the surface mesh\n");
+  printf("[OK] reading element node map { \n");
+    
+  H5::DataSet dataset_e = fp.openDataSet("/0/Mesh/MixedElements");
+  H5::DataType datatype_e = dataset_e.getDataType();
+  H5::DataSpace dataspace_e = dataset_e.getSpace();
+  
+  int rank = dataspace_e.getSimpleExtentNdims();
+
+  printf("[ok] rank of element node map: %d\n", rank);
+  
+  hsize_t *dims_e = new hsize_t[rank];
+  dataspace_e.getSimpleExtentDims(dims_e);
+
+  for (int i = 0; i < (int) rank; i ++)
+      printf("[ok] dims[%d] of element node map: %lld\n", i, dims_e[i]);
+  
+  int *data_e = new int[dims_e[0]*dims_e[1]];
+
+  // structure of hdf5 is <element type> <node 1> <node 2> <node 3>
+  size_t nx_e = 4;
+  size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
+
+  int **md_e = new int*[ny_e];
+  
+  for (size_t i = 0; i < ny_e; i ++)
+    md_e[i] = data_e + i * nx_e;
+  
+  printf("[ok] read the element map... ");
+
+  dataset_e.read(data_e, H5::PredType::STD_I32LE);
+
+  printf("done\n");
+
+  #if VERBOSE
+  printf("[ok] element map:\n");
+
+  size_t maxprint = 10;
+
+  if (ny_e < maxprint)
+    maxprint = ny_e;
+  
+  for (size_t i = 0; i < maxprint; i ++)
+    {
+      printf("[ok]");
+      for (size_t j = 0; j < nx_e; j ++)
+	  printf(" %d ", md_e[i][j]);
+      printf("\n");
+    }
+
+  printf("[ok] ...\n");
+  #endif
+  
+  printf("[OK] } done\n");
+
+  printf("[OK] reading nodes { \n");
+
+  H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes");
+  H5::DataType datatype_n = dataset_n.getDataType();
+  H5::DataSpace dataspace_n = dataset_n.getSpace();
+
+  rank = dataspace_n.getSimpleExtentNdims();
+
+  printf("[ok] rank of nodes: %d\n", rank);
+  
+  hsize_t *dims_n = new hsize_t[rank];
+  dataspace_n.getSimpleExtentDims(dims_n);
+
+  for (int i = 0; i < (int) rank; i ++)
+    printf("[ok] dims[%d] of nodes: %lld\n", i, dims_n[i]);
+
+  double *data_n = new double[dims_n[0]*dims_n[1]];
+  double **md_n  = new double*[dims_n[0]];
+
+  for (size_t i = 0; i < dims_n[0]; i ++)
+      md_n[i] = data_n + i * dims_n[1];
+
+  printf("[ok] read the nodes... ");
+
+  dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
+
+  printf("done\n");
+
+  #if VERBOSE
+  printf("[ok] nodes:\n");
+
+  maxprint = 10;
+  if (dims_n[0] < 10)
+    maxprint = dims_n[0];
+
+  for (size_t i = 0; i < maxprint; i ++)
+    {
+      printf("[ok]");
+      for (size_t j = 0; j < dims_n[1]; j ++)
+	  printf(" %f ", md_n[i][j]);
+      printf("\n");
+    }
+
+  printf("[ok] ...\n");
+  #endif
+  
+  printf("[OK] } done\n");
+
+  printf("[OK] building cells {\n");
+
+  Mesh::n_cells = ny_e;
+  Mesh::cells = new Cell[Mesh::n_cells];
+
+  double *coords1 = new double[dims_n[1]];
+  double *coords2 = new double[dims_n[1]];
+  double *coords3 = new double[dims_n[1]];
+
+  double *coords_avg = new double[dims_n[1]];
+
+  double *coords_max = new double[DIM];
+  double *coords_min = new double[DIM];
 
+  for (size_t i = 0; i < DIM; i ++)
+    {
+      coords_max[i] = DBL_MIN;
+      coords_min[i] = DBL_MAX;
+    }
+
+  for (size_t i = 0; i < ny_e; i ++)
+    {
+
+      cells[i].init_neighbors(SURFACE_N_EDGES);
+      
+      int n1 = md_e[i][1];
+      int n2 = md_e[i][2];
+      int n3 = md_e[i][3];
+
+      for (size_t j = 0; j < dims_n[1]; j++)
+	{
+	  coords1[j] = md_n[n1][j];
+	  coords2[j] = md_n[n2][j];
+	  coords3[j] = md_n[n3][j];
+
+	  coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
+	}
+
+      for (size_t k = 0; k < DIM; k ++)
+	{
+	  coords_max[k] = (coords_avg[k] > coords_max[k]) ? coords_avg[k] : coords_max[k];
+	  coords_min[k] = (coords_avg[k] < coords_min[k]) ? coords_avg[k] : coords_min[k];
+	}
+
+      cells[i].set_cell_id((int) i);
+      cells[i].set_coordinates(coords_avg);
+
+    }
+
+  set_min_coordinates(coords_min);
+  set_max_coordinates(coords_max);
+
+  printf("[ok] min: %f %f %f, max: %f %f %f\n", get_xmin(), get_ymin(), get_zmin(), get_xmax(), get_ymax(), get_zmax());
+
+  #if VERBOSE
+
+  maxprint = 10;
+
+  if (ny_e < maxprint)
+    maxprint = ny_e;
+
+  for (size_t i = 0; i < maxprint; i ++)
+    {
+      printf("[ok] cell[%ld] centroid: %f %f %f\n", i, cells[i].get_x(), cells[i].get_y(), cells[i].get_z());
+    }
+
+  printf("[ok] ...\n");
+  
+  #endif
+
+  delete[] coords_avg;
+  delete[] coords1;
+  delete[] coords2;
+  delete[] coords3;
+  
+  printf("[OK] } done\n");
+
+  printf("[OK] initialize edges {\n");
+  
   int ns[N_EDGES];
   hsize_t **szudzik_ids = new hsize_t*[ny_e];
   
@@ -383,15 +740,22 @@ int Mesh::build(H5::H5File fp)
     {
       for (size_t j = 0; j < N_EDGES; j ++)
 	{
+	  
 	  e0 = szudzik_ids[i][j];
 	  cells[i].set_neighbor(j, -1);
+
 	  for (size_t k = 0; k < ny_e; k ++)
 	    {
+	      
 	      for (size_t m = 0; m < N_EDGES; m ++)
 		{
+		  
 		  e1 = szudzik_ids[k][m];
-		  if (e0 == e1 && i != k)
+
+		  if (e0 == e1 && i != k) {
 		    cells[i].set_neighbor(j, (hsize_t) k);
+		  }
+		  
 		}
 	    }
 	}
@@ -442,7 +806,7 @@ int Mesh::read_data(H5::H5File dp)
   printf("[OK] reading velocity x { \n");
   printf("[ok] # of objects: %lld\n", group.getNumObjs());
 
-  H5::DataSet dataset     = group.openDataSet("11977");
+  H5::DataSet dataset     = group.openDataSet("70");
   H5::DataType datatype   = dataset.getDataType();
   H5::DataSpace dataspace = dataset.getSpace();
 
@@ -478,8 +842,11 @@ int Mesh::read_data(H5::H5File dp)
 
   #if VERBOSE
   maxprint = (dims[0] > 10) ? 10 : dims[0];
-  for (size_t i = 0; i < maxprint; i ++)
-    printf("[ok] %f\n", cells[i].get_v(0));
+  for (size_t i = 0; i < dims[0]; i ++) {
+    if (cells[i].get_v(0) > 0.0) {
+      printf("[ok] %f\n", cells[i].get_v(0));
+    }
+  }
   #endif
   
   printf("[OK] } done\n");
@@ -491,7 +858,7 @@ int Mesh::read_data(H5::H5File dp)
   printf("[OK] reading velocity y { \n");
   printf("[ok] # of objects: %lld\n", group_vy.getNumObjs());
 
-  H5::DataSet dataset_vy     = group_vy.openDataSet("11977");
+  H5::DataSet dataset_vy     = group_vy.openDataSet("70");
   H5::DataType datatype_vy   = dataset_vy.getDataType();
   H5::DataSpace dataspace_vy = dataset_vy.getSpace();
 
@@ -539,7 +906,7 @@ int Mesh::read_data(H5::H5File dp)
   printf("[OK] reading ponded water depth { \n");
   printf("[ok] # of objects: %lld\n", group_vy.getNumObjs());
 
-  H5::DataSet dataset_h     = group_h.openDataSet("11977");
+  H5::DataSet dataset_h     = group_h.openDataSet("70");
   H5::DataType datatype_h   = dataset_h.getDataType();
   H5::DataSpace dataspace_h = dataset_h.getSpace();
 
@@ -594,3 +961,4 @@ int Mesh::read_data(H5::H5File dp)
   return 0;
   
 }
+
diff --git a/mesh.hpp b/mesh.hpp
index a01fab7..3d34c13 100644
--- a/mesh.hpp
+++ b/mesh.hpp
@@ -44,6 +44,8 @@ public:
 
   int set_cell_id(int value);
 
+  int init_neighbors(size_t n_neighbors);
+
 protected:
 
 private:
@@ -58,7 +60,7 @@ private:
   double h;
   double v[3];
 
-  int neighbors[N_EDGES];
+  std::vector<hsize_t> neighbors;
   
 };
 
@@ -69,7 +71,8 @@ public:
   
   Mesh();
 
-  int build(H5::H5File fp);
+  int build_surface(H5::H5File fp);
+  int build_subsurface(H5::H5File fp);
   int get_n_cells();
   int read_data(H5::H5File dp);