implement surface environment

6 files changed, 373 insertions, 138 deletions

diff --git a/new/agent.cpp b/new/agent.cpp
index 8a68653..a6bb9e9 100644
--- a/new/agent.cpp
+++ b/new/agent.cpp
@@ -66,7 +66,7 @@ void Agent::advance(Mesh* mesh, size_t step)
 
 void Agent::advance_subsurface(Mesh* mesh, size_t step)
 {
-  Cell *location = mesh->get_cell(cell_index);
+  Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG);
   
   std::vector<double> velocity(DIM);
   std::vector<double> coordinates(DIM);
@@ -103,11 +103,11 @@ void Agent::advance_subsurface(Mesh* mesh, size_t step)
   double distance_norm = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]);
   
   // check if cell has changed
-  size_t n_cells = mesh->get_n_cells();
+  size_t n_cells = mesh->get_n_cells(SUBSURFACE_FLAG);
   for (size_t i = 0; i < n_cells; i ++) {
 
       for (size_t j = 0; j < DIM; j ++) {
-	coordinates_neighbor[j] = mesh->get_cell(i)->get_coordinate(j);
+	coordinates_neighbor[j] = mesh->get_cell(i, SUBSURFACE_FLAG)->get_coordinate(j);
 	distance[j] = coordinates_neighbor[j] - coordinates[j];
       }
 
diff --git a/new/cell.cpp b/new/cell.cpp
index 5d9c4b6..c872442 100644
--- a/new/cell.cpp
+++ b/new/cell.cpp
@@ -76,6 +76,10 @@ int Cell::get_cell_id() {
 void Cell::set_coupled_neighbor(hsize_t value) {
   coupled_neighbor = value;
 }
+
+hsize_t Cell::get_coupled_neighbor() {
+  return coupled_neighbor;
+}
   
 void Cell::set_water_depth(double value) {
   water_depth = value;
diff --git a/new/header.hpp b/new/header.hpp
index f4e4e70..b8ef032 100644
--- a/new/header.hpp
+++ b/new/header.hpp
@@ -5,6 +5,7 @@
 
 #include <iostream>
 #include <vector>
+#include <string>    // string class
 #include <omp.h>
 #include <algorithm> // min
 #include <float.h>   // DBL_MAX
@@ -27,15 +28,21 @@
 #define N_VERTICES 6
 #define N_EDGES 9
 
-#define SURFACE_MESH_FILE "../test/visdump_surface_mesh.h5"
+#define SURFACE_NCOL 4
+#define SUBSURFACE_NCOL 7
 
+#define SURFACE_MESH_FILE "../test/visdump_surface_mesh.h5"
+#define SURFACE_DATA_FILE "../test/visdump_surface_data.h5"
 #define SUBSURFACE_MESH_FILE "../test/visdump_mesh.h5"
 #define SUBSURFACE_DATA_FILE "../test/visdump_data.h5"
+#define SURFACE_DATA_GROUP_NAME_H5 "11841"
 #define SUBSURFACE_DATA_GROUP_NAME_H5 "11841"
 
 #define MESH_MIXED_ELEMENTS "/0/Mesh/MixedElements"
 #define MESH_NODES "/0/Mesh/Nodes"
 
+#define AREADIV 0.01
+#define FUZZ 0.1
 #define RADIUS 200.0
 
 #define SUBSURFACE_FLAG 1
diff --git a/new/main.cpp b/new/main.cpp
index 468b765..83aaf13 100644
--- a/new/main.cpp
+++ b/new/main.cpp
@@ -13,20 +13,26 @@ int main(int argc, char* argv[])
   fp_subsurface.close();
 
   H5::H5File fp_surface = H5::H5File(SURFACE_MESH_FILE, H5F_ACC_RDONLY);
+  mesh.build_surface_mesh(fp_surface);
+  fp_surface.close();
   
-  FILE *mesh_fp = fopen("mesh.csv", "w");
-  fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
-
-  Cell *c;
-  for (int n = 0; n < mesh.get_n_cells(); n ++) {
-    c = mesh.get_cell(n);
-    fprintf(mesh_fp, "%f, %f, %f, %d\n", c->get_coordinate(0),
-	    c->get_coordinate(1), c->get_coordinate(2), n);
-  }
-  fclose(mesh_fp);
+  //  FILE *mesh_fp = fopen("mesh.csv", "w");
+  //  fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
+
+//  Cell *c;
+//  for (int n = 0; n < mesh.get_n_cells(); n ++) {
+//    c = mesh.get_cell(n);
+//    fprintf(mesh_fp, "%f, %f, %f, %d\n", c->get_coordinate(0),
+//	    c->get_coordinate(1), c->get_coordinate(2), n);
+//  }
+//  fclose(mesh_fp);
 
   H5::H5File dp_subsurface = H5::H5File(SUBSURFACE_DATA_FILE, H5F_ACC_RDONLY);
-  mesh.build_subsurface_environment(dp_subsurface);
+  mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG);
+  dp_subsurface.close();
+
+  H5::H5File dp_surface = H5::H5File(SURFACE_DATA_FILE, H5F_ACC_RDONLY);
+  mesh.build_environment(dp_surface, SURFACE_FLAG);
   dp_subsurface.close();
 
   // Place agents in the domain and advance N_ITERATION steps
@@ -34,7 +40,7 @@ int main(int argc, char* argv[])
 
   // Bootstrap 
   srand(SEED);
-  int len = mesh.get_n_cells();
+  int len = mesh.get_n_cells(SUBSURFACE_FLAG);
   int index;
 
   std::cout << "[~>] Releasing " << N_AGENT << " agents" << std::endl;
@@ -43,16 +49,16 @@ int main(int argc, char* argv[])
 
     index = rand()%len;
     agents[i].set_cell_index(index);
-    agents[i].set_path(mesh.get_cell(index)->get_coordinate(0),
-		       mesh.get_cell(index)->get_coordinate(1),
-		       mesh.get_cell(index)->get_coordinate(2), 0);
+    agents[i].set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
+		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
+		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 0);
   
     agents[i].set_domain(SUBSURFACE_FLAG, 0);
     agents[i].set_domain(SUBSURFACE_FLAG, 1);
 
-    agents[i].set_path(mesh.get_cell(index)->get_coordinate(0),
-		       mesh.get_cell(index)->get_coordinate(1),
-		       mesh.get_cell(index)->get_coordinate(2), 1);
+    agents[i].set_path(mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(0),
+		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(1),
+		       mesh.get_cell(index, SUBSURFACE_FLAG)->get_coordinate(2), 1);
 
     size_t step;
 
diff --git a/new/mesh.cpp b/new/mesh.cpp
index cc69638..6180ea8 100644
--- a/new/mesh.cpp
+++ b/new/mesh.cpp
@@ -2,21 +2,146 @@
 
 #include "mesh.hpp"
 
+
+// ------------------------------------------------------------------------
+// -- PUBLIC METHODS
+// ------------------------------------------------------------------------
+
+void Mesh::build_surface_mesh(H5::H5File fp)
+{
+  std::cout << "[~>] Building the surface mesh" << std::endl;
+
+  get_mesh_metadata(fp, SURFACE_NCOL);
+
+  handle_surface_node_data(fp);
+
+  build_surface_cells();
+  
+}
+
+void Mesh::build_subsurface_mesh(H5::H5File fp)
+{
+  
+  std::cout << "[~>] Building the subsurface mesh" << std::endl;
+
+  get_mesh_metadata(fp, SUBSURFACE_NCOL);
+
+  handle_subsurface_node_data(fp);
+
+  build_subsurface_cells();
+  
+}
+
+void Mesh::build_environment(H5::H5File ep, int flag)
+{
+
+  std::string groupname;
+  std::string datasetname;
+  
+  // Read in velocities in x direction
+    
+  if (flag == SUBSURFACE_FLAG) {
+    std::cout << "[~>] Building the subsurface environment" << std::endl;
+    groupname = "darcy_velocity.cell.0";
+    datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
+  } else if (flag == SURFACE_FLAG) {
+    std::cout << "[~>] Building the surface environment" << std::endl;
+    groupname = "surface-velocity.cell.0";
+    datasetname = SURFACE_DATA_GROUP_NAME_H5;
+  }
+
+  read_velocity(ep, groupname, datasetname, flag, 0);
+
+  // Read in velocities in y direction
+
+  if (flag == SUBSURFACE_FLAG) {
+    groupname = "darcy_velocity.cell.1";
+    datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
+  } else if (flag == SURFACE_FLAG) {
+    groupname = "surface-velocity.cell.1";
+    datasetname = SURFACE_DATA_GROUP_NAME_H5;
+  }
+
+  read_velocity(ep, groupname, datasetname, flag, 1);
+
+  // Read velocities in z direction
+  
+  if (flag == SUBSURFACE_FLAG) {
+    groupname = "darcy_velocity.cell.2";
+    datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
+  } else if (flag == SURFACE_FLAG) {
+    groupname = "surface-surface_subsurface_flux.cell.0";
+    datasetname = SURFACE_DATA_GROUP_NAME_H5;
+  }
+
+  read_velocity(ep, groupname, datasetname, flag, 2);
+
+  
+  // Surface environment needs adjustments to the vertical velocity
+  // and needs to read in water depth
+  
+  if (flag == SURFACE_FLAG) {
+    
+    // Rescale the exchange flux from mole to m/s
+    for (int i = 0; i < n_cells_surface; i ++) {
+      surface_cells[i].set_velocity(2, surface_cells[i].get_velocity(2) * 1.8e-5 * AREADIV);
+      if (i % (n_cells_surface / MAXPRINT) == 0)
+  	printf("[~~~~>] Surface velocity in cell[%d] = { %e, %e, %e } m/s\n", i,
+	       surface_cells[i].get_velocity(0),
+	       surface_cells[i].get_velocity(1),
+	       surface_cells[i].get_velocity(2));
+    }
+
+    // Read in water depth (TODO)
+    read_water_depth(ep);
+    
+  } else if (flag == SUBSURFACE_FLAG) {
+    
+    for (int i = 0; i < n_cells_subsurface; i ++) {
+      if (i % (n_cells_subsurface / MAXPRINT) == 0)
+  	printf("[~~~~>] Darcy velocity in cell[%d] = { %e, %e, %e } m/s\n", i,
+	       subsurface_cells[i].get_velocity(0),
+	       subsurface_cells[i].get_velocity(1),
+	       subsurface_cells[i].get_velocity(2));
+    }
+    
+  }
+  
+}
+
+// ------------------------------------------------------------------------
+// GETTERS AND SETTERS
+// ------------------------------------------------------------------------
+
 Mesh::Mesh()
 {
-  n_cells = 0;
+  n_cells_surface = 0;
+  n_cells_subsurface = 0;
   high_coordinates.resize(DIM);
   low_coordinates.resize(DIM);
 }
 
-int Mesh::get_n_cells()
+int Mesh::get_n_cells(int flag)
 {
+  int n_cells;
+  
+  if (flag == SUBSURFACE_FLAG) {
+    n_cells = n_cells_subsurface;
+  } else {
+    n_cells = n_cells_subsurface;
+  }
+  
   return n_cells;
 }
 
-Cell* Mesh::get_cell(int index)
+Cell* Mesh::get_cell(int index, int flag)
 {
-  return &(cells[index]);
+
+  if (flag == SUBSURFACE_FLAG) {
+    return &(subsurface_cells[index]);
+  } else {
+    return &(surface_cells[index]);
+  } 
 }
 
 void Mesh::set_high_coordinate(int index, double value)
@@ -39,138 +164,122 @@ double Mesh::get_low_coordinate(int index)
   return low_coordinates[index];
 }
 
-void Mesh::build_surface_mesh(H5::H5File fp)
-{
-  std::cout << "[~>] Building the surface mesh" << std::endl;
 
-  H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS);
-  H5::DataType datatype_e = dataset_e.getDataType();
-  H5::DataSpace dataspace_e = dataset_e.getSpace();
-  
-}
 
-void Mesh::build_subsurface_mesh(H5::H5File fp)
-{
-  
-  std::cout << "[~>] Building the subsurface mesh" << std::endl;
-  
-  H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS);
-  H5::DataType datatype_e = dataset_e.getDataType();
-  H5::DataSpace dataspace_e = dataset_e.getSpace();
-  
-  int rank = dataspace_e.getSimpleExtentNdims();
+// ------------------------------------------------------------------------
+// -- PRIVATE METHODS
+// ------------------------------------------------------------------------
 
-  hsize_t *dims_e = new hsize_t[rank];
-  dataspace_e.getSimpleExtentDims(dims_e);
-
-  for (int i = 0; i < rank; i ++)
-    std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
+void Mesh::handle_surface_node_data(H5::H5File fp)
+{
 
-  int *data_e = new int[dims_e[0]*dims_e[1]];
-  dataset_e.read(data_e, H5::PredType::STD_I32LE);
+  H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
+  H5::DataType datatype_n = dataset_n.getDataType();
+  H5::DataSpace dataspace_n = dataset_n.getSpace();
 
-  size_t nx_e = 7;
-  size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
-  
-  element_node_indices.resize(ny_e);
+  int rank = dataspace_n.getSimpleExtentNdims();
 
-  for (size_t i = 0; i < ny_e; i ++) {
-    element_node_indices[i].resize(nx_e);
-    for (size_t j = 0; j < nx_e; j ++)
-      element_node_indices[i][j] = data_e[i * nx_e + j];
-  }
+  hsize_t *dims_n = new hsize_t[rank];
+  dataspace_n.getSimpleExtentDims(dims_n);
 
-  n_cells = ny_e;
-  cells.resize(ny_e);
+  for (int i = 0; i < rank; i ++)
+    std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
 
-  std::cout << "[~~~>] Number of cells: " << n_cells << std::endl;
+  double *data_n = new double[dims_n[0]*dims_n[1]];
+  dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
 
-  // Nodes
-  handle_subsurface_node_data(fp);
+  nodes.resize(dims_n[0]);
 
-  FILE *node_fp = fopen("nodes.csv", "w");
-  fprintf(node_fp, "x(m), y(m), z(m), index\n");
-  for (int n = 0; n < nodes.size(); n ++) {
-    fprintf(node_fp, "%f, %f, %f, %d\n", nodes[n][0], nodes[n][1], nodes[n][2], n);
+  for (size_t i = 0; i < dims_n[0]; i ++) {
+    nodes[i].resize(dims_n[1]);
+    for (size_t j = 0; j < dims_n[1]; j ++)
+      nodes[i][j] = data_n[i * dims_n[1] + j];
   }
-  fclose(node_fp);
-
-  // Cells
-  build_subsurface_cells();
   
 }
 
-void Mesh::build_surface_environment(H5::H5File ep)
+void Mesh::build_surface_cells()
 {
-  std::cout << "[~>] Building the surface environment" << std::endl;
-}
 
-void Mesh::build_subsurface_environment(H5::H5File ep)
-{
-  std::cout << "[~>] Building the subsurface environment" << std::endl;
+  std::vector<double> coords1(DIM);
+  std::vector<double> coords2(DIM);
+  std::vector<double> coords3(DIM);
+  std::vector<double> coords_avg(DIM);
 
-  // Read in Darcy velocity in x direction
-  
-  H5::Group group1         = ep.openGroup("darcy_velocity.cell.0");
-  H5::DataSet dataset1     = group1.openDataSet(SUBSURFACE_DATA_GROUP_NAME_H5);
-  H5::DataType datatype1   = dataset1.getDataType();
-  H5::DataSpace dataspace1 = dataset1.getSpace();
+  for (int i = 0; i < n_cells_surface; i ++) {
 
-  int rank1 = dataspace1.getSimpleExtentNdims();
+    int n1 = element_node_indices[i][1];
+    int n2 = element_node_indices[i][2];
+    int n3 = element_node_indices[i][3];
 
-  hsize_t *dims1 = new hsize_t[rank1];
-  dataspace1.getSimpleExtentDims(dims1);
+    std::vector<double> edge_length(DIM);
 
-  for (int i = 0; i < rank1; i ++)
-    std::cout << "[~~~>] Dimensions in " << i << " of Darcy velocity in x: " << dims1[i] << std::endl;
+    for (size_t k = 0; k < DIM; k ++)
+      edge_length[k] = DBL_MAX;
 
-  double *data1 = new double[dims1[0]*dims1[1]];
-  dataset1.read(data1, H5::PredType::IEEE_F64LE);
+    for (size_t j = 0; j < DIM; j ++) {
 
-  for (int i = 0; i < n_cells; i ++)
-    cells[i].set_velocity(0, data1[i]);
+      coords1[j] = nodes[n1][j];
+      coords2[j] = nodes[n2][j];
+      coords3[j] = nodes[n3][j];
 
-  H5::Group group2         = ep.openGroup("darcy_velocity.cell.1");
-  H5::DataSet dataset2     = group2.openDataSet(SUBSURFACE_DATA_GROUP_NAME_H5);
-  H5::DataType datatype2   = dataset2.getDataType();
-  H5::DataSpace dataspace2 = dataset2.getSpace();
+      coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
 
-  int rank2 = dataspace2.getSimpleExtentNdims();
-  hsize_t *dims2 = new hsize_t[rank2];
-  dataspace2.getSimpleExtentDims(dims2);
-  
-  for (int i = 0; i < rank2; i ++)
-    std::cout << "[~~~>] Dimensions in " << i << " of Darcy velocity in y: " << dims2[i] << std::endl;
+      edge_length[0] = fabs(coords1[j] - coords_avg[j]);
+      edge_length[1] = fabs(coords2[j] - coords_avg[j]);
+      edge_length[2] = fabs(coords3[j] - coords_avg[j]);
+      
+    }
+
+    double smallest_distance = edge_length[0];
+    double largest_distance  = edge_length[0];
+    
+    for (size_t k = 1; k < DIM; k ++) {
+      smallest_distance = std::min(smallest_distance, edge_length[k]);
+      largest_distance  = std::max(largest_distance, edge_length[k]);
+    }
 
-  double *data2 = new double[dims2[0]*dims2[1]];
-  dataset2.read(data2, H5::PredType::IEEE_F64LE);
+    surface_cells[i].set_characteristic_length(smallest_distance);
+    surface_cells[i].set_maximum_length(largest_distance);
 
-  for (int i = 0; i < n_cells; i ++)
-    cells[i].set_velocity(1, data2[i]);
+    for (size_t k = 0; k < DIM; k ++)
+      surface_cells[i].set_coordinate(k, coords_avg[k]);
 
-  H5::Group group3         = ep.openGroup("darcy_velocity.cell.2");
-  H5::DataSet dataset3     = group3.openDataSet(SUBSURFACE_DATA_GROUP_NAME_H5);
-  H5::DataType datatype3   = dataset3.getDataType();
-  H5::DataSpace dataspace3 = dataset3.getSpace();
+    surface_cells[i].set_cell_id((int) i);
 
-  int rank3 = dataspace3.getSimpleExtentNdims();
-  hsize_t *dims3 = new hsize_t[rank3];
-  dataspace3.getSimpleExtentDims(dims3);
-  
-  for (int i = 0; i < rank3; i ++)
-    std::cout << "[~~~>] Dimensions in " << i << " of Darcy velocity in z: " << dims2[i] << std::endl;
+    // find the coupled cell in the subsurface
+    
+    std::vector<double> distance(DIM);
+    hsize_t coupled_neighbor = -1; // store current coupled neighbor here
+    double dz = DBL_MAX;       // find the nearest cell in the vertical direction
+    for (int k = 0; k < n_cells_subsurface; k ++) {
+
+      for (size_t j = 0; j < DIM; j ++)
+	distance[j] = subsurface_cells[k].get_coordinate(j) - coords_avg[j];
+
+      double xy_distance = sqrt(distance[0]*distance[0] + distance[1]*distance[1]);
+
+      if (xy_distance <= FUZZ) {
+	dz = std::min(dz, fabs(distance[2]));
+	coupled_neighbor = (hsize_t) k;
+	if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length())
+	  break;
+      }
+      
+    }
 
-  double *data3 = new double[dims3[0]*dims3[1]];
-  dataset3.read(data3, H5::PredType::IEEE_F64LE);
+    surface_cells[i].set_coupled_neighbor(coupled_neighbor);
+    subsurface_cells[coupled_neighbor].set_coupled_neighbor(i);
 
-  for (int i = 0; i < n_cells; i ++) {
-    cells[i].set_velocity(2, data3[i]);
-    if (i % (n_cells / MAXPRINT) == 0)
-      printf("[~~~~>] Darcy velocity in cell[%d] = { %e, %e, %e }\n", i, cells[i].get_velocity(0), cells[i].get_velocity(1), cells[i].get_velocity(2));
+    if (i%(n_cells_surface/5)==0) {
+      std::cout << "[~~~>] Coupling surface and subsurface meshes" << std::endl;
+      std::cout << "[~~~~>] Subsurface neighbor of cell " << i << " is cell " << surface_cells[i].get_coupled_neighbor() << std::endl;
+      std::cout << "[~~~~>] Surface neighbor of cell " << coupled_neighbor << " is cell " << subsurface_cells[coupled_neighbor].get_coupled_neighbor() << std::endl;
+    }
+    
   }
-  
-}
 
+}
 
 void Mesh::handle_subsurface_node_data(H5::H5File fp)
 {
@@ -217,7 +326,7 @@ void Mesh::build_subsurface_cells()
     set_low_coordinate(i, DBL_MAX);
   }
 
-  for (int i = 0; i < n_cells; i ++) {
+  for (int i = 0; i < n_cells_subsurface; i ++) {
 
     int n1 = element_node_indices[i][1];
     int n2 = element_node_indices[i][2];
@@ -261,8 +370,8 @@ void Mesh::build_subsurface_cells()
       largest_distance  = std::max(largest_distance, edge_length[k]);
     }
 
-    cells[i].set_characteristic_length(smallest_distance);
-    cells[i].set_maximum_length(largest_distance);
+    subsurface_cells[i].set_characteristic_length(smallest_distance);
+    subsurface_cells[i].set_maximum_length(largest_distance);
     
     for (size_t k = 0; k < DIM; k ++) {
 
@@ -270,12 +379,111 @@ void Mesh::build_subsurface_cells()
       high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]);
       low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]);
 
-      cells[i].set_coordinate(k, coords_avg[k]);
+      subsurface_cells[i].set_coordinate(k, coords_avg[k]);
       
     }
     
-    cells[i].set_cell_id((int) i);
+    subsurface_cells[i].set_cell_id((int) i);
     
   }
   
 }
+
+void Mesh::get_mesh_metadata(H5::H5File fp, int flag)
+{
+
+  H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS);
+  H5::DataType datatype_e = dataset_e.getDataType();
+  H5::DataSpace dataspace_e = dataset_e.getSpace();
+
+  int rank = dataspace_e.getSimpleExtentNdims();
+
+  hsize_t *dims_e = new hsize_t[rank];
+  dataspace_e.getSimpleExtentDims(dims_e);
+
+  for (int i = 0; i < rank; i ++)
+      std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
+
+  int *data_e = new int[dims_e[0]*dims_e[1]];
+  dataset_e.read(data_e, H5::PredType::STD_I32LE);
+
+  size_t nx_e = flag;
+  size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
+  
+  element_node_indices.resize(ny_e);
+
+  for (size_t i = 0; i < ny_e; i ++) {
+    element_node_indices[i].resize(nx_e);
+    for (size_t j = 0; j < nx_e; j ++)
+      element_node_indices[i][j] = data_e[i * nx_e + j];
+  }
+  
+  if (flag == SUBSURFACE_NCOL) {
+    n_cells_subsurface = ny_e;
+    subsurface_cells.resize(ny_e);
+  } else if (flag == SURFACE_NCOL) {
+    n_cells_surface = ny_e;
+    surface_cells.resize(ny_e);
+  }
+
+  std::cout << "[~~~>] Number of cells: " << ny_e << std::endl;
+  
+}
+
+void Mesh::read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index)
+{
+  
+  H5::Group group		= fp.openGroup(groupname.c_str());
+  H5::DataSet dataset		= group.openDataSet(datasetname.c_str());
+  H5::DataType datatype		= dataset.getDataType();
+  H5::DataSpace dataspace	= dataset.getSpace();
+  
+  int rank = dataspace.getSimpleExtentNdims();
+
+  hsize_t *dims = new hsize_t[rank];
+  dataspace.getSimpleExtentDims(dims);
+
+  for (int i = 0; i < rank; i ++)
+    std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+  
+  double *data = new double[dims[0]*dims[1]];
+  dataset.read(data, H5::PredType::IEEE_F64LE);
+
+  if (flag == SUBSURFACE_FLAG) {
+    for (int i = 0; i < n_cells_subsurface; i ++)
+      subsurface_cells[i].set_velocity(index, data[i]);
+  } else if (flag == SURFACE_FLAG) {
+    for (int i = 0; i < n_cells_surface; i ++)
+      surface_cells[i].set_velocity(index, data[i]);
+  }
+
+}
+
+void Mesh::read_water_depth(H5::H5File fp)
+{
+
+  H5::Group group		= fp.openGroup("surface-ponded_depth.cell.0");
+  H5::DataSet dataset		= group.openDataSet(SURFACE_DATA_GROUP_NAME_H5);
+  H5::DataType datatype		= dataset.getDataType();
+  H5::DataSpace dataspace	= dataset.getSpace();
+  
+  int rank = dataspace.getSimpleExtentNdims();
+
+  hsize_t *dims = new hsize_t[rank];
+  dataspace.getSimpleExtentDims(dims);
+
+  for (int i = 0; i < rank; i ++)
+    std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+  
+  double *data = new double[dims[0]*dims[1]];
+  dataset.read(data, H5::PredType::IEEE_F64LE);
+
+  for (int i = 0; i < n_cells_surface; i ++) {
+    surface_cells[i].set_water_depth(data[i]);
+    if (i%(n_cells_surface/5)==0) {
+      printf("[~~~~>] Surface ponded depth in cell[%d] =  %e m\n", i,
+	     surface_cells[i].get_water_depth());
+    }
+  }
+  
+}
diff --git a/new/mesh.hpp b/new/mesh.hpp
index 2ee2680..dc26ac1 100644
--- a/new/mesh.hpp
+++ b/new/mesh.hpp
@@ -16,33 +16,43 @@ public:
   void build_surface_mesh(H5::H5File fp);
   void build_subsurface_mesh(H5::H5File fp);
 
-  void build_surface_environment(H5::H5File ep);
+  void build_environment(H5::H5File ep, int flag);
   void build_subsurface_environment(H5::H5File ep);
 
-  Cell* get_cell(int index);
+  Cell* get_cell(int index, int flag);
 
   void set_high_coordinate(int index, double value);
   double get_high_coordinate(int index);
   void set_low_coordinate(int index, double value);
   double get_low_coordinate(int index);
 
-  int get_n_cells();
+  int get_n_cells(int flag);
 
 private:
 
-  int n_cells;
-  std::vector<Cell> cells;
+  int n_cells_subsurface;
+  int n_cells_surface;
+  
+  std::vector<Cell> subsurface_cells;
+  std::vector<Cell> surface_cells;
 
   std::vector<double> high_coordinates;
   std::vector<double> low_coordinates;
 
   std::vector<std::vector<double>> nodes;
   std::vector<std::vector<int>> element_node_indices;
-  std::vector<std::vector<std::vector<hsize_t>>> edge_ids;
+  //std::vector<std::vector<std::vector<hsize_t>>> edge_ids;
   
   void handle_subsurface_node_data(H5::H5File fp);
   void build_subsurface_cells();
-  
+
+  void handle_surface_node_data(H5::H5File fp);
+  void build_surface_cells();
+
+  void get_mesh_metadata(H5::H5File fp, int ncol);
+  void read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index);
+  void read_water_depth(H5::H5File fp);
+    
 };
 
 #endif