diff options
Diffstat (limited to 'new/mesh.cpp')
| -rw-r--r-- | new/mesh.cpp | 109 |
1 files changed, 76 insertions, 33 deletions
diff --git a/new/mesh.cpp b/new/mesh.cpp index 6180ea8..6ae6c22 100644 --- a/new/mesh.cpp +++ b/new/mesh.cpp @@ -7,7 +7,7 @@ // -- PUBLIC METHODS // ------------------------------------------------------------------------ -void Mesh::build_surface_mesh(H5::H5File fp) +void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep) { std::cout << "[~>] Building the surface mesh" << std::endl; @@ -15,7 +15,7 @@ void Mesh::build_surface_mesh(H5::H5File fp) handle_surface_node_data(fp); - build_surface_cells(); + build_surface_cells(ep); } @@ -92,7 +92,7 @@ void Mesh::build_environment(H5::H5File ep, int flag) surface_cells[i].get_velocity(2)); } - // Read in water depth (TODO) + // Read in water depth read_water_depth(ep); } else if (flag == SUBSURFACE_FLAG) { @@ -128,7 +128,7 @@ int Mesh::get_n_cells(int flag) if (flag == SUBSURFACE_FLAG) { n_cells = n_cells_subsurface; } else { - n_cells = n_cells_subsurface; + n_cells = n_cells_surface; } return n_cells; @@ -164,8 +164,6 @@ double Mesh::get_low_coordinate(int index) return low_coordinates[index]; } - - // ------------------------------------------------------------------------ // -- PRIVATE METHODS // ------------------------------------------------------------------------ @@ -198,16 +196,20 @@ void Mesh::handle_surface_node_data(H5::H5File fp) } -void Mesh::build_surface_cells() +void Mesh::build_surface_cells(H5::H5File ep) { - std::vector<double> coords1(DIM); - std::vector<double> coords2(DIM); - std::vector<double> coords3(DIM); - std::vector<double> coords_avg(DIM); + std::vector<double> coords1(DIM-1); + std::vector<double> coords2(DIM-1); + std::vector<double> coords3(DIM-1); + std::vector<double> coords_avg(DIM-1); + + read_elevation(ep); for (int i = 0; i < n_cells_surface; i ++) { + //std::cout << "build cell " << i << std::endl; + int n1 = element_node_indices[i][1]; int n2 = element_node_indices[i][2]; int n3 = element_node_indices[i][3]; @@ -217,7 +219,11 @@ void Mesh::build_surface_cells() for (size_t k = 0; k < DIM; k ++) edge_length[k] = DBL_MAX; - for (size_t j = 0; j < DIM; j ++) { + + double smallest_distance; + double largest_distance; + + for (size_t j = 0; j < DIM-1; j ++) { coords1[j] = nodes[n1][j]; coords2[j] = nodes[n2][j]; @@ -228,40 +234,45 @@ void Mesh::build_surface_cells() 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]); + smallest_distance = edge_length[0]; + 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]); + } + } surface_cells[i].set_characteristic_length(smallest_distance); surface_cells[i].set_maximum_length(largest_distance); - for (size_t k = 0; k < DIM; k ++) + for (size_t k = 0; k < DIM-1; k ++) surface_cells[i].set_coordinate(k, coords_avg[k]); - + surface_cells[i].set_cell_id((int) i); // 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 + int coupled_neighbor = -1; // store current coupled neighbor here + double dz = DBL_MAX; + 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]; + distance[j] = subsurface_cells[k].get_coordinate(j) - surface_cells[i].get_coordinate(j); + // Because the mesh is a stacked mesh, the only candidates for the subsurface + // neighbor are the cells located in the same column. These cells will have + // the same centroid in the xy-plane as the surface cell. The problem reduces + // to finding the nearest cell in the vertical direction in the column. + 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; + coupled_neighbor = k; if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length()) break; } @@ -326,6 +337,9 @@ void Mesh::build_subsurface_cells() set_low_coordinate(i, DBL_MAX); } + double smallest_distance; + double largest_distance; + for (int i = 0; i < n_cells_subsurface; i ++) { int n1 = element_node_indices[i][1]; @@ -359,15 +373,15 @@ void Mesh::build_subsurface_cells() edge_length[3] = fabs(coords4[j] - coords_avg[j]); edge_length[4] = fabs(coords5[j] - coords_avg[j]); edge_length[5] = fabs(coords6[j] - coords_avg[j]); - - } - double smallest_distance = edge_length[0]; - double largest_distance = edge_length[0]; - - for (size_t k = 1; k < N_VERTICES; k ++) { + smallest_distance = edge_length[0]; + largest_distance = edge_length[0]; + + for (size_t k = 1; k < N_VERTICES; k ++) { smallest_distance = std::min(smallest_distance, edge_length[k]); largest_distance = std::max(largest_distance, edge_length[k]); + } + } subsurface_cells[i].set_characteristic_length(smallest_distance); @@ -487,3 +501,32 @@ void Mesh::read_water_depth(H5::H5File fp) } } + +void Mesh::read_elevation(H5::H5File fp) +{ + + H5::Group group = fp.openGroup("surface-elevation.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_coordinate(2, data[i]); + if (i%(n_cells_surface/5)==0) { + printf("[~~~~>] Surface elevation in cell[%d] = %e m\n", i, + surface_cells[i].get_coordinate(2)); + } + } + +} |