summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIlhan Özgen Xian <iozgen@lbl.gov>2021-05-25 11:17:08 -0700
committerIlhan Özgen Xian <iozgen@lbl.gov>2021-05-25 11:17:08 -0700
commitdf6adc77e0cbe568832f57212929a9acfb568a4c (patch)
tree0ab3aa5ed3468e011c344ccf655358b806cf0bd2
parent2833f2e67ca27ae6c0b6d1062a1de1be36701e75 (diff)
clean up
-rw-r--r--agent.cpp323
-rw-r--r--agent.hpp59
-rw-r--r--cell.cpp (renamed from new/cell.cpp)0
-rw-r--r--cell.hpp (renamed from new/cell.hpp)0
-rw-r--r--header.hpp69
-rw-r--r--main.cpp203
-rw-r--r--makefile10
-rwxr-xr-xmerge-results.sh (renamed from new/merge-results.sh)0
-rw-r--r--mesh.cpp1334
-rw-r--r--mesh.hpp120
-rw-r--r--new/agent.cpp245
-rw-r--r--new/agent.hpp41
-rw-r--r--new/header.hpp56
-rw-r--r--new/main.cpp133
-rw-r--r--new/makefile33
-rw-r--r--new/mesh.cpp546
-rw-r--r--new/mesh.hpp59
-rwxr-xr-xtest.cpp140
-rw-r--r--tools.cpp11
-rw-r--r--tools.hpp8
20 files changed, 758 insertions, 2632 deletions
diff --git a/agent.cpp b/agent.cpp
index edf5b8d..97d2e69 100644
--- a/agent.cpp
+++ b/agent.cpp
@@ -1,196 +1,245 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode: c++ -*-
-#include "header.hpp"
+#include "agent.hpp"
-HydroAgent::HydroAgent()
+// ------------------------------------------------------------------------
+// -- PUBLIC METHODS
+// ------------------------------------------------------------------------
+
+void Agent::advance(Mesh* mesh, size_t step)
{
- Cell *location;
+
+ switch(domain[step-1]) {
+ case SURFACE_FLAG:
+ advance_surface(mesh,step);
+ break;
+ case SUBSURFACE_FLAG:
+ advance_subsurface(mesh,step);
+ break;
+ default:
+ set_path(path[step-1][0], path[step-1][1], path[step-1][2], step);
+ time[step] = time[step-1];
+ domain[step] = OUTSIDE_FLAG;
+ break;
+ }
- probability_i = 0.0;
- probability_a = 0.0;
}
-Cell* HydroAgent::get_location()
+// ------------------------------------------------------------------------
+// GETTERS AND SETTERS
+// ------------------------------------------------------------------------
+
+Agent::Agent()
{
- return location;
+ time.resize(N_ITERATION);
+ domain.resize(N_ITERATION);
+ path.resize(N_ITERATION);
+ for (size_t i = 0; i < N_ITERATION; i ++)
+ path[i].resize(DIM);
}
-int HydroAgent::set_location(Cell *cell)
+void Agent::set_cell_index(int value)
{
- location = cell;
- return 0;
+ cell_index = value;
}
-int HydroAgent::set_path(double xcrd, double ycrd, double zcrd, size_t step)
+void Agent::set_domain(int value, size_t step)
{
- x[step] = xcrd;
- y[step] = ycrd;
- z[step] = zcrd;
- return 0;
+ domain[step] = value;
}
-int HydroAgent::advance(Mesh *mesh, size_t step)
+double Agent::get_time(size_t step)
{
- int flag = 0;
-
- if (domain[step-1] == SURFACE_FLAG) {
- flag = advance_surface(mesh, step);
- } else if (domain[step-1] == SUBSURFACE_FLAG) {
- flag = advance_subsurface(mesh, step);
- } else if (domain[step-1] == OUTOFBOUNDS_FLAG) {
- flag = 1;
- x[step] = x[step-1];
- y[step] = y[step-1];
- set_domain(OUTOFBOUNDS_FLAG, step);
- #if VERBOSE
- printf("[ok] agent is out of bounds. do not advect.\n");
- #endif
- }
+ return time[step];
+}
- return flag;
+void Agent::set_path(double x, double y, double z, size_t step)
+{
+ path[step][0] = x;
+ path[step][1] = y;
+ path[step][2] = z;
}
-int HydroAgent::advance_surface(Mesh *mesh, size_t step)
+double Agent::get_timestepsize()
{
+ return time_step_size;
+}
+
+double Agent::get_path(size_t step, size_t index)
+{
+ return path[step][index];
+}
+
+int Agent::get_domain(size_t step)
+{
+ return domain[step];
+}
+
+// ------------------------------------------------------------------------
+// -- PRIVATE METHODS
+// ------------------------------------------------------------------------
- double vx = location->get_v(0);
- double vy = location->get_v(1);
+void Agent::advance_surface(Mesh* mesh, size_t step)
+{
+ Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG);
- double x_agent = x[step-1] + DT_SURFACE * vx;
- double y_agent = y[step-1] + DT_SURFACE * vy;
+ std::vector<double> velocity(DIM);
+ std::vector<double> coordinates(DIM-1);
+ std::vector<double> coordinates_neighbor(DIM-1);
+ std::vector<double> distance(DIM-1);
+
+ domain[step] = SURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
- x[step] = x_agent;
- y[step] = y_agent;
+ for (size_t i = 0; i < DIM; i ++)
+ velocity[i] = location->get_velocity(i);
- double x_centroid = location->get_x();
- double y_centroid = location->get_y();
+ // CFL criteria
+ double length = location->get_characteristic_length();
+ double velocity_norm = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]);
- double dx = x_centroid - x_agent;
- double dy = y_centroid - y_agent;
+ if (velocity_norm > EPS) {
+ time_step_size = CFL * length / velocity_norm;
+ } else {
+ time_step_size = MAXTIMESTEP;
+ }
- double distance = sqrt(dx * dx + dy * dy);
+ time[step] = time[step-1] + time_step_size;
- int cell_index = -1;
+ // check if agent infiltrates
+ double h = location->get_water_depth();
+ int flag = INFILTRATE_FALSE;
- for (int n = 0; n < N_EDGES; n ++) {
+ if (h > EPS) { // only wet cell can infiltrate
+ double qe = std::min(velocity[2], 0.0);
+ double pi = fabs(qe) * time_step_size / h;
+ double trial = (double) (rand()%100) / 100.0;
+ if (trial < pi)
+ flag = INFILTRATE_TRUE;
+ }
- hsize_t neighbor_index = location->get_neighbor(n);
+ if (flag == INFILTRATE_FALSE) {
- double x_neighbor = mesh->get_cells()[neighbor_index].get_x();
- double y_neighbor = mesh->get_cells()[neighbor_index].get_y();
+ // advect the agent horizontally, which means the last dimension
+ // (z-direction) is omitted
+ for (size_t i = 0; i < DIM-1; i ++) {
+ coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
+ distance[i] = location->get_coordinate(i) - coordinates[i];
+ }
- dx = x_neighbor - x_agent;
- dy = y_neighbor - y_agent;
+ // horizontal distance to the centroid of current cell
+ double distance_norm = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
- double _distance = sqrt(dx * dx + dy * dy);
+ // check if cell has changed
+ size_t n_cells = mesh->get_n_cells(SURFACE_FLAG);
+ for (size_t i = 0; i < n_cells; i ++) {
- if (distance > _distance) {
- cell_index = n;
- distance = _distance;
+ for (size_t j = 0; j < DIM-1; j ++) {
+ coordinates_neighbor[j] = mesh->get_cell(i, SURFACE_FLAG)->get_coordinate(j);
+ distance[j] = coordinates_neighbor[j] - coordinates[j];
+ }
+
+ double distance_norm_ = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
+
+ if (distance_norm > distance_norm_) {
+ cell_index = i;
+ distance_norm = distance_norm_;
}
}
- if (cell_index >= 0) {
- location = &(mesh->get_cells()[location->get_neighbor(cell_index)]);
+ // check if out of the domain: if the agent travelled thrice the
+ // maximum length of the cell and still can't find a neighboring
+ // cell, it must be advected out of the domain
+ if (distance_norm > 10.0 * location->get_maximum_length())
+ domain[step] = OUTSIDE_FLAG;
+
+ } else {
+ // infiltration detected
+ domain[step] = SUBSURFACE_FLAG;
+ coordinates[0] = path[step-1][0];
+ coordinates[1] = path[step-1][1];
+ cell_index = mesh->get_cell(cell_index, SURFACE_FLAG)->get_coupled_neighbor();
}
- if ((x_agent > mesh->get_xmax()) || (x_agent < mesh->get_xmin()) ||
- (y_agent > mesh->get_ymax()) || (y_agent < mesh->get_ymin())) {
- domain[step] = OUTOFBOUNDS_FLAG;
- } else {
- domain[step] = SURFACE_FLAG;
- }
+ set_path(coordinates[0], coordinates[1], location->get_coordinate(2), step);
- return 1;
}
-int HydroAgent::advance_subsurface(Mesh *mesh, size_t step)
-{
- double vx = location->get_v(0);
- double vy = location->get_v(1);
- double vz = location->get_v(2);
+void Agent::advance_subsurface(Mesh* mesh, size_t step)
+{
+ Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG);
- double x_agent = x[step-1] + DT_SUBSURFACE * vx;
- double y_agent = y[step-1] + DT_SUBSURFACE * vy;
- double z_agent = z[step-1] + DT_SUBSURFACE * vz;
+ std::vector<double> velocity(DIM);
+ std::vector<double> coordinates(DIM);
+ std::vector<double> coordinates_neighbor(DIM);
+ std::vector<double> distance(DIM);
- x[step] = x_agent;
- y[step] = y_agent;
- z[step] = z_agent;
+ domain[step] = SUBSURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
+
+ for (size_t i = 0; i < DIM; i ++)
+ velocity[i] = location->get_velocity(i);
- double x_centroid = location->get_x();
- double y_centroid = location->get_y();
- double z_centroid = location->get_z();
+ // CFL criteria
+ double length = location->get_characteristic_length();
+ double velocity_norm = sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]);
- double dx = x_centroid - x_agent;
- double dy = y_centroid - y_agent;
- double dz = z_centroid - z_agent;
+ if (velocity_norm > EPS) {
+ time_step_size = CFL * length / velocity_norm;
+ } else {
+ time_step_size = MAXTIMESTEP;
+ }
- double distance = sqrt(dx * dx + dy * dy + dz * dz);
+ time[step] = time[step-1] + time_step_size;
- int cell_index = -1;
+ // advect the agent
+ for (size_t i = 0; i < DIM; i ++) {
+ coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
+ distance[i] = location->get_coordinate(i) - coordinates[i];
+ }
+
+ set_path(coordinates[0], coordinates[1], coordinates[2], step);
- for (int n = 0; n < SUBSURFACE_N_NEIGHBORS; n ++) {
+ // distance to the centroid of current cell
+ 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(SUBSURFACE_FLAG);
+ for (size_t i = 0; i < n_cells; i ++) {
- hsize_t neighbor_index = location->get_neighbor(n);
+ for (size_t j = 0; j < DIM; j ++) {
+ coordinates_neighbor[j] = mesh->get_cell(i, SUBSURFACE_FLAG)->get_coordinate(j);
+ distance[j] = coordinates_neighbor[j] - coordinates[j];
+ }
- double x_neighbor = mesh->get_cells()[neighbor_index].get_x();
- double y_neighbor = mesh->get_cells()[neighbor_index].get_y();
- double z_neighbor = mesh->get_cells()[neighbor_index].get_z();
+ double distance_norm_ = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]);
- dx = x_neighbor - x_agent;
- dy = y_neighbor - y_agent;
- dz = z_neighbor - z_agent;
+ if (distance_norm > distance_norm_) {
+ cell_index = i;
+ distance_norm = distance_norm_;
+ }
+
+ }
- double _distance = sqrt(dx * dx + dy * dy + dz * dz);
+ // check if out of the domain: if the agent travelled thrice the
+ // maximum length of the cell and still can't find a neighboring
+ // cell, it must be advected out of the domain
+ if (distance_norm > 10.0 * location->get_maximum_length()) {
+ domain[step] = OUTSIDE_FLAG;
+ }
- if (distance > _distance) {
- cell_index = n;
- distance = _distance;
+ // check if exfiltrates
+ int neighbor_index = mesh->get_cell(cell_index, SUBSURFACE_FLAG)->get_coupled_neighbor();
+ double z_surface;
+ if (neighbor_index > -1) {
+ // border to surface detected
+ z_surface = mesh->get_cell(neighbor_index, SURFACE_FLAG)->get_coordinate(2);
+ if (coordinates[2] >= z_surface) {
+ // exfiltration detected
+ cell_index = (int) neighbor_index;
+ domain[step] = SURFACE_FLAG;
+ set_path(coordinates[0], coordinates[1], z_surface, step);
}
-
}
- if (cell_index >= 0) {
- location = &(mesh->get_cells()[location->get_neighbor(cell_index)]);
- }
-
- if ((x_agent > mesh->get_xmax()) || (x_agent < mesh->get_xmin()) ||
- (y_agent > mesh->get_ymax()) || (y_agent < mesh->get_ymin()) ||
- (z_agent > mesh->get_zmax()) || (z_agent < mesh->get_zmin())) {
- domain[step] = OUTOFBOUNDS_FLAG;
- } else {
- domain[step] = SURFACE_FLAG;
- }
-
- return 2;
-}
-
-int HydroAgent::set_domain(int domainflag, size_t step)
-{
- domain[step] = domainflag;
-
- return 0;
-}
-
-int HydroAgent::get_domain(size_t step)
-{
- return domain[step];
-}
-
-double* HydroAgent::get_x()
-{
- return &x[0];
-}
-
-double* HydroAgent::get_y()
-{
- return &y[0];
-}
-
-double* HydroAgent::get_z()
-{
- return &z[0];
}
diff --git a/agent.hpp b/agent.hpp
index 65f778f..f99c2ab 100644
--- a/agent.hpp
+++ b/agent.hpp
@@ -1,55 +1,40 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode: c++ -*-
#ifndef AGENT_HPP
#define AGENT_HPP
#include "header.hpp"
+#include "cell.hpp"
+#include "mesh.hpp"
-// agent.hpp
-// =========
-// HydroAgent to be used in agent-based modeling to track flow paths
-// on a given flow field.
-// At each time step, the HydroAgent needs to decide whether to
-// infiltrate or to get advected following the flow. This is based
-// on computed probabilities, see (Reaney, 2008).
-// First, check if the agent infiltrates/exfiltrates. If not, check
-// if it advects. If yes, advect according to flow velocity.
-
-class HydroAgent
+class Agent
{
public:
-
- HydroAgent();
-
- Cell* get_location();
- int set_location(Cell *cell);
- int set_path(double x, double y, double z, size_t step);
- int set_domain(int domainflag, size_t step);
- int get_domain(size_t step);
-
- int advance(Mesh *mesh, size_t step);
- int advance_surface(Mesh *mesh, size_t step);
- int advance_subsurface(Mesh *mesh, size_t step);
- double *get_x();
- double *get_y();
- double *get_z();
+ Agent();
-protected:
+ void set_cell_index(int value);
+ void set_domain(int value, size_t step);
+ void set_path(double x, double y, double z, size_t step);
+ double get_path(size_t step, size_t index);
+ int get_domain(size_t step);
+ double get_timestepsize();
+ double get_time(size_t step);
+
+ void advance(Mesh* mesh, size_t step);
private:
- Cell *location;
+ int cell_index;
+ double time_step_size;
- int domain[N_ITERATION];
-
- double x[N_ITERATION];
- double y[N_ITERATION];
- double z[N_ITERATION];
-
- double probability_i; // infiltration probability
- double probability_a; // advection probability
+ std::vector<int> domain;
+ std::vector<double> time;
+ std::vector<std::vector<double>> path;
+
+ void advance_surface(Mesh *mesh, size_t step);
+ void advance_subsurface(Mesh *mesh, size_t step);
};
diff --git a/new/cell.cpp b/cell.cpp
index c872442..c872442 100644
--- a/new/cell.cpp
+++ b/cell.cpp
diff --git a/new/cell.hpp b/cell.hpp
index 7c2cd83..7c2cd83 100644
--- a/new/cell.hpp
+++ b/cell.hpp
diff --git a/header.hpp b/header.hpp
index 24922d9..cff75d2 100644
--- a/header.hpp
+++ b/header.hpp
@@ -1,49 +1,56 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode:c++ -*-
-#ifndef HEADER_HPP
+#ifndef CONSTANTS_HPP
+#define CONSTANTS_HPP
#include <iostream>
+#include <vector>
+#include <string> // string class
+//#include <omp.h>
+#include <algorithm> // min
+#include <float.h> // DBL_MAX
+#include <math.h> // fabs, sqrt
+#include <mpi.h>
#include <H5Cpp.h>
#include <H5File.h>
-#include "agent.hpp"
-#include "mesh.hpp"
-#include "tools.hpp"
-
-#include <float.h>
-#include <math.h>
-#include <vector>
+#define MAXPRINT 5
+#define MAXTIMESTEP 1.0e5
+#define EPS 1.0e-10
-#define VERBOSE 1
-#define VERBOSE_HIGH 0
+#define N_AGENT 100
+#define N_ITERATION 2000
+#define SEED 108
+#define CFL 0.5
-#define N_EDGES 3
-#define SURFACE_N_EDGES 3
-#define SUBSURFACE_N_EDGES 9
-#define SUBSURFACE_N_NODES 6
-#define SUBSURFACE_N_NEIGHBORS 5 // 5 von Neumann neighbors, 11 if Moore neighborhood
#define DIM 3
+#define N_NEIGHBORS 5
+#define N_VERTICES 6
+#define N_EDGES 9
-#define RADIUS 200.0
+#define SURFACE_NCOL 4
+#define SUBSURFACE_NCOL 7
-#define N_AGENT 300
-#define N_ITERATION 5000
-#define CFL 0.5
-#define DT_SURFACE 1800.0
-#define DT_SUBSURFACE 1800.0
+#define INFILTRATE_TRUE 1
+#define INFILTRATE_FALSE 2
-#define SUBSURFACE 1
+#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 SURFACE_FLAG 1
-#define SUBSURFACE_FLAG 2
-#define OUTOFBOUNDS_FLAG 3
+#define MESH_MIXED_ELEMENTS "/0/Mesh/MixedElements"
+#define MESH_NODES "/0/Mesh/Nodes"
-#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"
+#define AREADIV 0.01
+#define FUZZ 0.1
+#define RADIUS 200.0
-#define DATA_GROUP_NAME_H5 "11841"
+#define SUBSURFACE_FLAG 1
+#define SURFACE_FLAG 2
+#define OUTSIDE_FLAG 3
#endif
diff --git a/main.cpp b/main.cpp
index e5b4301..1e45aac 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,138 +1,133 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode: c++ -*-
-#include "header.hpp"
+#include "mesh.hpp"
+#include "agent.hpp"
-int main(void)
+int main(int argc, char* argv[])
{
- printf("\n\n");
- printf(" / /^ agent provocateur\n");
- printf(" o--o =================\n");
- printf("\n\n");
-
- printf("[OK] # of agents: %d, # of iterations: %d\n[OK] size of surface time step: %f s\n[OK] size of subsurface time step: %f s\n\n", N_AGENT, N_ITERATION, DT_SURFACE, DT_SUBSURFACE);
-
- H5::H5File fp = H5::H5File(MESH_FILE_NAME, H5F_ACC_RDONLY);
-
+ int rank;
+ int size;
+
+ MPI_Init(&argc, &argv);
+ MPI_Comm_size(MPI_COMM_WORLD, &size);
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
Mesh mesh = Mesh();
- mesh.build_surface(fp);
-
- printf("[OK] # of cells: %d\n\n", mesh.get_n_cells());
- fp.close();
+ H5::H5File fp_subsurface = H5::H5File(SUBSURFACE_MESH_FILE, H5F_ACC_RDONLY);
+ H5::H5File dp_subsurface = H5::H5File(SUBSURFACE_DATA_FILE, H5F_ACC_RDONLY);
+ H5::H5File fp_surface = H5::H5File(SURFACE_MESH_FILE, H5F_ACC_RDONLY);
+ H5::H5File dp_surface = H5::H5File(SURFACE_DATA_FILE, H5F_ACC_RDONLY);
-#if SUBSURFACE
- H5::H5File subsurface_fp = H5::H5File(SUBSURFACE_MESH_FILE_NAME, H5F_ACC_RDONLY);
+ if (rank == 0)
+ std::cout << "[~>] Building the subsurface mesh" << std::endl;
- Mesh subsurface_mesh = Mesh();
- subsurface_mesh.build_subsurface(subsurface_fp);
+ mesh.build_subsurface_mesh(fp_subsurface);
+ fp_subsurface.close();
- printf("[OK] # of cells: %d\n\n", subsurface_mesh.get_n_cells());
-
- subsurface_fp.close();
-#endif
-
- H5::H5File dp = H5::H5File(DATA_FILE_NAME, H5F_ACC_RDONLY);
- mesh.read_surface_data(dp);
- dp.close();
-
-#if SUBSURFACE
- H5::H5File sdp = H5::H5File(SUBSURFACE_DATA_FILE_NAME, H5F_ACC_RDONLY);
- subsurface_mesh.read_subsurface_data(sdp);
- sdp.close();
-#endif
+ if (rank == 0) {
+ std::cout << "[~~>] Number of cells: " << mesh.get_n_cells(SUBSURFACE_FLAG) << std::endl;
+ std::cout << "[~>] Building the surface mesh (includes coupling with subsurface and may take a while)" << std::endl;
+ }
- // Place agents in the domain and advance N_ITERATION steps
- std::vector<HydroAgent> agents(N_AGENT);
+ mesh.build_surface_mesh(fp_surface, dp_surface);
+ fp_surface.close();
- srand(108);
+ if (rank == 0) {
- int len = mesh.get_n_cells();
+ int len_surf = mesh.get_n_cells(SURFACE_FLAG);
- int index;
+ std::cout << "[~~>] Number of cells: " << len_surf << std::endl;
+
+ FILE *mesh_fp = fopen("mesh.csv", "w");
+ fprintf(mesh_fp, "x (m), y (m), z (m), index\n");
+ Cell *c;
- printf("[OK] inject agents\n");
-
- for (size_t i = 0; i < N_AGENT; i ++) {
- index = rand() % len;
-
- 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 (int n = 0; n < len_surf; n ++) {
+ c = mesh.get_cell(n, SURFACE_FLAG);
+ fprintf(mesh_fp, "%f, %f, %f, %d\n",
+ c->get_coordinate(0),
+ c->get_coordinate(1),
+ c->get_coordinate(2), n);
}
-
+
+ fclose(mesh_fp);
+
+ std::cout << "[~>] Reading the subsurface environment" << std::endl;
+
}
+
+ mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG);
+ dp_subsurface.close();
- for (size_t step = 1; step < N_ITERATION; step ++) {
+ mesh.build_environment(dp_surface, SURFACE_FLAG);
+ dp_subsurface.close();
- if (step % 10 == 0) {
- printf("[OK] time step %08lu {\n", step);
- }
-
- for (size_t i = 0; i < N_AGENT; i ++) {
-
- #if VERBOSE
- if (step % 100 == 0) {
- if (i < 10) {
- printf("[ok] advance agent #%lu\n", i+1);
- } else if (i == 10) {
- printf("[ok] ...\n");
- }
- }
- #endif
-
- agents[i].advance(&mesh, step);
-
- }
-
- if (step % 10 == 0) {
- printf("[OK] } done\n");
- }
-
- }
+ // Place agents in the domain and advance N_ITERATION steps
- printf("[OK] simulation completed\n\n");
+ // Bootstrap
+ srand(SEED + rank);
+ int len = mesh.get_n_cells(SUBSURFACE_FLAG);
+ int index;
- printf("[OK] write out agent data {\n");
+ std::cout << "[~>] Releasing " << N_AGENT << " agents" << std::endl;
- char output_fname[1024];
+ size_t partition = N_AGENT / size;
+ size_t low_bound = rank * partition;
+ size_t high_bound = (rank + 1) * partition;
+ if (rank == size - 1) // top processor picks up extras
+ high_bound = N_AGENT;
- // write out the results
- for (size_t i = 0; i < N_AGENT; i ++) {
+ for (size_t i = low_bound; i < high_bound; i ++) {
- sprintf(output_fname, "agent-%06lu.dat", i);
- FILE *result_fp = fopen(output_fname, "w");
+ Agent agent;
+
+ index = rand()%len;
+ agent.set_cell_index(index);
+ agent.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);
+
+ agent.set_domain(SUBSURFACE_FLAG, 0);
+ agent.set_domain(SUBSURFACE_FLAG, 1);
- #if VERBOSE
- if (i < 10) {
- printf("[ok] write out %s\n", output_fname);
- } else if (i == 10) {
- printf("[ok] ...\n");
- }
- #endif
+ agent.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;
- double* x = agents[i].get_x();
- double* y = agents[i].get_y();
+ char output_fname[1024];
- for (size_t step = 0; step < N_ITERATION; step ++) {
- fprintf(result_fp, "%f %f %f\n", step * DT_SURFACE, x[step], y[step]);
+ // move agent
+ for (step = 1; step < N_ITERATION; step ++) {
+ if (step % (N_ITERATION / 5) == 0 && i % (N_AGENT / 5) == 0) {
+ std::cout << "[~~~>] agent[" << i << "] step = " << step << " dt = " << agent.get_timestepsize() << " s" << std::endl;
+ }
+ agent.advance(&mesh, step);
}
+ // write out results
+ sprintf(output_fname, "results/agent-%06lu.csv", i);
+ FILE *result_fp = fopen(output_fname, "w");
+
+ fprintf(result_fp, "time (s), x (m), y (m), z (m), domain\n");
+ for (step = 0; step < N_ITERATION; step ++) {
+ fprintf(result_fp, "%f, %f, %f, %f, %d\n",
+ agent.get_time(step),
+ agent.get_path(step, 0),
+ agent.get_path(step, 1),
+ agent.get_path(step, 2),
+ agent.get_domain(step));
+ }
fclose(result_fp);
-
+
}
- printf("[OK] } done\n");
- printf("[OK] finish\n\n");
-
- return 0;
+ MPI_Finalize();
+ return MPI_SUCCESS;
+
}
diff --git a/makefile b/makefile
index ffdfe46..4a91e28 100644
--- a/makefile
+++ b/makefile
@@ -1,10 +1,10 @@
-PRG=agent
+PRG=agent-provocateur
-CC=g++-11
-#CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
-CXXFLAGS=-Wall -Wextra -pedantic -O0
+#CC=g++-11
+CC=mpic++
+CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
-SRC=mesh.cpp main.cpp tools.cpp agent.cpp
+SRC=cell.cpp mesh.cpp main.cpp agent.cpp
LIB=-lm -lhdf5 -lhdf5_cpp
OBJ=$(SRC:.cpp=.o)
diff --git a/new/merge-results.sh b/merge-results.sh
index 5459897..5459897 100755
--- a/new/merge-results.sh
+++ b/merge-results.sh
diff --git a/mesh.cpp b/mesh.cpp
index f3eb381..8112ede 100644
--- a/mesh.cpp
+++ b/mesh.cpp
@@ -1,1128 +1,546 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode: c++ -*-
-#include "header.hpp"
+#include "mesh.hpp"
-Cell::Cell()
-{
- x = 0.0;
- y = 0.0;
- z = 0.0;
-}
-Cell::Cell(double *coordinates)
-{
- x = coordinates[0];
- y = coordinates[1];
- z = coordinates[2];
-}
+// ------------------------------------------------------------------------
+// -- PUBLIC METHODS
+// ------------------------------------------------------------------------
-Cell::Cell(const Cell &cell)
+void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep)
{
- x = cell.x;
- y = cell.y;
- z = cell.z;
-}
-int Cell::init_neighbors(size_t n_neighbors)
-{
- neighbors.resize(n_neighbors);
-
- return 0;
-}
+ get_mesh_metadata(fp, SURFACE_NCOL);
-int Cell::set_coordinates(double* coordinates)
-{
- x = coordinates[0];
- y = coordinates[1];
- z = coordinates[2];
+ handle_surface_node_data(fp);
- return 0;
+ build_surface_cells(ep);
+
}
-int Cell::set_neighbors(hsize_t *value_array)
+void Mesh::build_subsurface_mesh(H5::H5File fp)
{
- for (size_t i = 0; i < neighbors.size(); i ++)
- neighbors[i] = value_array[i];
-
- return 0;
-}
-int Cell::set_neighbor(int index, hsize_t value)
-{
- neighbors[index] = value;
- return 0;
-}
+ get_mesh_metadata(fp, SUBSURFACE_NCOL);
-int Cell::set_h(double value)
-{
- h = value;
- return 0;
-}
+ handle_subsurface_node_data(fp);
-int Cell::set_v(int index, double value)
-{
- v[index] = value;
- return 0;
+ build_subsurface_cells();
+
}
-int Cell::set_cell_id(int value)
+void Mesh::build_environment(H5::H5File ep, int flag)
{
- cell_id = value;
- return 0;
-}
-hsize_t Cell::get_neighbor(int index)
-{
- return neighbors[index];
-}
+ std::string groupname;
+ std::string datasetname;
+
+ // Read in velocities in x direction
+
+ if (flag == SUBSURFACE_FLAG) {
+ groupname = "darcy_velocity.cell.0";
+ datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
+ } else if (flag == SURFACE_FLAG) {
+ groupname = "surface-velocity.cell.0";
+ datasetname = SURFACE_DATA_GROUP_NAME_H5;
+ }
-double Cell::get_x()
-{
- return x;
-}
+ read_velocity(ep, groupname, datasetname, flag, 0);
-double Cell::get_y()
-{
- return y;
-}
+ // Read in velocities in y direction
-double Cell::get_z()
-{
- return z;
-}
+ 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;
+ }
-double Cell::get_h()
-{
- return h;
-}
+ read_velocity(ep, groupname, datasetname, flag, 1);
-double Cell::get_v(int index)
-{
- return v[index];
-}
+ // 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;
+ }
-int Cell::get_cell_id()
-{
- return cell_id;
-}
+ read_velocity(ep, groupname, datasetname, flag, 2);
-Mesh::Mesh()
-{
- Mesh::n_cells = 0;
-}
+
+ // 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 VERBOSE
+ 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));
+#endif
+ }
-int Mesh::get_n_cells()
-{
- return Mesh::n_cells;
+ // Read in water depth
+ read_water_depth(ep);
+
+ }
+#if VERBOSE
+ 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));
+ }
+
+ }
+#endif
+
}
-Cell* Mesh::get_cells()
-{
- return Mesh::cells;
-}
+// ------------------------------------------------------------------------
+// GETTERS AND SETTERS
+// ------------------------------------------------------------------------
-int Mesh::set_min_coordinates(double *coordinates_min)
+Mesh::Mesh()
{
- xmin = coordinates_min[0];
- ymin = coordinates_min[1];
- zmin = coordinates_min[2];
-
- return 0;
+ n_cells_surface = 0;
+ n_cells_subsurface = 0;
+ high_coordinates.resize(DIM);
+ low_coordinates.resize(DIM);
}
-int Mesh::set_max_coordinates(double *coordinates_max)
+int Mesh::get_n_cells(int flag)
{
- xmax = coordinates_max[0];
- ymax = coordinates_max[1];
- zmax = coordinates_max[2];
-
- return 0;
+ int n_cells;
+
+ if (flag == SUBSURFACE_FLAG) {
+ n_cells = n_cells_subsurface;
+ } else {
+ n_cells = n_cells_surface;
+ }
+
+ return n_cells;
}
-double Mesh::get_xmax()
+Cell* Mesh::get_cell(int index, int flag)
{
- return xmax;
-}
-double Mesh::get_ymax()
-{
- return ymax;
+ if (flag == SUBSURFACE_FLAG) {
+ return &(subsurface_cells[index]);
+ } else {
+ return &(surface_cells[index]);
+ }
}
-double Mesh::get_zmax()
+void Mesh::set_high_coordinate(int index, double value)
{
- return zmax;
+ high_coordinates[index] = value;
}
-double Mesh::get_xmin()
+double Mesh::get_high_coordinate(int index)
{
- return xmin;
+ return high_coordinates[index];
}
-double Mesh::get_ymin()
+void Mesh::set_low_coordinate(int index, double value)
{
- return ymin;
+ low_coordinates[index] = value;
}
-double Mesh::get_zmin()
+double Mesh::get_low_coordinate(int index)
{
- return zmin;
+ return low_coordinates[index];
}
-#if SUBSURFACE
-int Mesh::build_subsurface(H5::H5File fp)
-{
-
- 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::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]];
-
- 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 ++)
- md_e[i] = data_e + i * nx_e;
-
- printf("[ok] read the element map... ");
-
- dataset_e.read(data_e, H5::PredType::STD_I32LE);
+// ------------------------------------------------------------------------
+// -- PRIVATE METHODS
+// ------------------------------------------------------------------------
- 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");
+void Mesh::handle_surface_node_data(H5::H5File fp)
+{
- H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes");
+ H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
H5::DataType datatype_n = dataset_n.getDataType();
H5::DataSpace dataspace_n = dataset_n.getSpace();
- rank = dataspace_n.getSimpleExtentNdims();
+ int 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]);
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
+#endif
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];
+ nodes.resize(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
+ 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];
+ }
- printf("[OK] } done\n");
-
- printf("[OK] building cells {\n");
-
- Mesh::n_cells = ny_e;
- Mesh::cells = new Cell[ny_e];
-
- 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]];
-
- 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(SUBSURFACE_N_NEIGHBORS);
-
- 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];
- 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 ++)
- {
- 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;
- delete[] coords4;
- delete[] coords5;
- delete[] coords6;
-
- printf("[OK] } done\n");
+void Mesh::build_surface_cells(H5::H5File ep)
+{
- printf("[OK] initialize edges {\n");
-
- int ns[SUBSURFACE_N_NODES];
-
- // Szudzik IDs for very large number of edges is not safe
- 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 ++) {
-
- for (size_t j = 0; j < SUBSURFACE_N_NODES; j ++) {
- ns[j] = md_e[i][j+1];
- }
+ 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);
- // top of the cell
+ read_elevation(ep);
- 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];
- }
+ for (int i = 0; i < n_cells_surface; i ++) {
- 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];
- }
+ int n1 = element_node_indices[i][1];
+ int n2 = element_node_indices[i][2];
+ int n3 = element_node_indices[i][3];
- 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];
- }
+ std::vector<double> edge_length(DIM);
- // bottom of the cell
+ for (size_t k = 0; k < DIM; k ++)
+ edge_length[k] = DBL_MAX;
- 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];
- }
+ double smallest_distance;
+ double largest_distance;
+
+ for (size_t j = 0; j < DIM-1; j ++) {
- 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];
- }
+ coords1[j] = nodes[n1][j];
+ coords2[j] = nodes[n2][j];
+ coords3[j] = nodes[n3][j];
- // sides of the cell
+ coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
- 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];
- }
+ 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]);
- if (ns[1] < ns[4]) {
- edge_ids[i][7][0] = ns[1];
- edge_ids[i][7][1] = ns[4];
- } else {
- edge_ids[i][7][0] = ns[4];
- edge_ids[i][7][1] = ns[1];
- }
+ smallest_distance = edge_length[0];
+ largest_distance = edge_length[0];
- 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 k = 1; k < DIM; k ++) {
+ smallest_distance = std::min(smallest_distance, edge_length[k]);
+ largest_distance = std::max(largest_distance, edge_length[k]);
}
}
-
- printf("[OK] } done\n");
-
- printf("[OK] initialize neighbors {\n");
-
- hsize_t e00;
- hsize_t e01;
- hsize_t e10;
- hsize_t e11;
- hsize_t e20;
- hsize_t e21;
- hsize_t e30;
- hsize_t e31;
- hsize_t e40;
- hsize_t e41;
- hsize_t e50;
- hsize_t e51;
- hsize_t e60;
- hsize_t e61;
- hsize_t e70;
- hsize_t e71;
- hsize_t e80;
- hsize_t e81;
-
- hsize_t _e0;
- hsize_t _e1;
- 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();
+ surface_cells[i].set_characteristic_length(smallest_distance);
+ surface_cells[i].set_maximum_length(largest_distance);
+
+ for (size_t k = 0; k < DIM-1; k ++)
+ surface_cells[i].set_coordinate(k, coords_avg[k]);
- for (size_t j = 0; j < SUBSURFACE_N_NEIGHBORS; j ++)
- cells[i].set_neighbor(j, -1);
-
-#if VERBOSE_HIGH
- if (i % 100 == 0)
- printf("[ok] searching neighbors of cell: %ld/%ld\n", i+1, ny_e);
-#endif
+ surface_cells[i].set_cell_id((int) i);
- e00 = edge_ids[i][0][0];
- e01 = edge_ids[i][0][1];
+ // find the coupled cell in the subsurface
+ std::vector<double> distance(DIM);
+ int coupled_neighbor = -1; // store current coupled neighbor here
+ double dz = DBL_MAX;
- e10 = edge_ids[i][1][0];
- e11 = edge_ids[i][1][1];
+ for (int k = 0; k < n_cells_subsurface; k ++) {
- e20 = edge_ids[i][2][0];
- e21 = edge_ids[i][2][1];
+ for (size_t j = 0; j < DIM; j ++)
+ distance[j] = subsurface_cells[k].get_coordinate(j) - surface_cells[i].get_coordinate(j);
- e30 = edge_ids[i][3][0];
- e31 = edge_ids[i][3][1];
+ // 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.
- e40 = edge_ids[i][4][0];
- e41 = edge_ids[i][4][1];
-
- e50 = edge_ids[i][5][0];
- e51 = edge_ids[i][5][1];
+ double xy_distance = sqrt(distance[0]*distance[0] + distance[1]*distance[1]);
- e60 = edge_ids[i][6][0];
- e61 = edge_ids[i][6][1];
-
- e70 = edge_ids[i][7][0];
- e71 = edge_ids[i][7][1];
-
- e80 = edge_ids[i][8][0];
- e81 = edge_ids[i][8][1];
-
-#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 ++) {
-
- // we will limit our search to a 3D sphere with radius RADIUS around the cell
- double xk = cells[k].get_x();
- double yk = cells[k].get_y();
- double zk = cells[k].get_z();
-
- double distance = sqrt((xk - xi) * (xk - xi) + (yk - yi) * (yk - yi) + (zk - zi) * (zk - zi));
-
- if (i != k && distance <= RADIUS) {
-
- 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];
-
- 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)) {
-
- _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 %d/%d: %ld [%lld, %lld]\n", counter, SUBSURFACE_N_NEIGHBORS, k, _e0, _e1);
-#endif
- break;
- }
-
- }
- }
- }
-
- if (counter == SUBSURFACE_N_NEIGHBORS)
+ if (xy_distance <= FUZZ) {
+ dz = std::min(dz, fabs(distance[2]));
+ coupled_neighbor = k;
+ if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length())
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 < SUBSURFACE_N_NEIGHBORS; 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::read_subsurface_data(H5::H5File dp)
-{
- printf("[OK] build subsurface environment\n");
+ surface_cells[i].set_coupled_neighbor(coupled_neighbor);
+ subsurface_cells[coupled_neighbor].set_coupled_neighbor(i);
- // read darcy velocity x
-
- H5::Group group = dp.openGroup("darcy_velocity.cell.0");
-
- printf("[OK] reading velocity x{ \n");
- printf("[ok] # of objects: %lld\n", group.getNumObjs());
-
- H5::DataSet dataset = group.openDataSet(DATA_GROUP_NAME_H5);
- H5::DataType datatype = dataset.getDataType();
- H5::DataSpace dataspace = dataset.getSpace();
-
- int rank = dataspace.getSimpleExtentNdims();
-
- printf("[ok] rank of data: %d\n", rank);
-
- hsize_t *dims = new hsize_t[rank];
- dataspace.getSimpleExtentDims(dims);
-
- for (int i = 0; i < (int) rank; i ++)
- printf("[ok] dims[%d] of data: %lld\n", i, dims[i]);
-
- double *data = new double[dims[0]*dims[1]];
- double **md = new double*[dims[0]];
-
- for (size_t i = 0; i < dims[0]; i ++)
- md[i] = data + i * dims[1];
-
- dataset.read(data, H5::PredType::IEEE_F64LE);
-
#if VERBOSE
- size_t maxprint = (dims[0] > 10) ? 10 : dims[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", data[i]);
-#endif
-
- printf("[OK] } done\n");
- printf("[OK] set velocity x { \n");
-
- for (int i = 0; i < n_cells; i ++)
- cells[i].set_v(0, data[i]);
-
-#if VERBOSE
- maxprint = (dims[0] > 10) ? 10 : dims[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));
+ 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;
}
- }
#endif
-
- printf("[OK] } done\n");
- return 0;
-}
-
-#endif
-
-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
+void Mesh::handle_subsurface_node_data(H5::H5File fp)
+{
- printf("[OK] } done\n");
-
- printf("[OK] reading nodes { \n");
-
- H5::DataSet dataset_n = fp.openDataSet("/0/Mesh/Nodes");
+ H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
H5::DataType datatype_n = dataset_n.getDataType();
H5::DataSpace dataspace_n = dataset_n.getSpace();
- rank = dataspace_n.getSimpleExtentNdims();
+ int 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]);
-
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
+#endif
+
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];
+ nodes.resize(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
+ 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];
+ }
- printf("[OK] } done\n");
+}
- printf("[OK] building cells {\n");
+void Mesh::build_subsurface_cells()
+{
- Mesh::n_cells = ny_e;
- Mesh::cells = new Cell[Mesh::n_cells];
+ std::vector<double> coords1(DIM);
+ std::vector<double> coords2(DIM);
+ std::vector<double> coords3(DIM);
+ std::vector<double> coords4(DIM);
+ std::vector<double> coords5(DIM);
+ std::vector<double> coords6(DIM);
- double *coords1 = new double[dims_n[1]];
- double *coords2 = new double[dims_n[1]];
- double *coords3 = new double[dims_n[1]];
+ std::vector<double> coords_avg(DIM);
- double *coords_avg = new double[dims_n[1]];
+ for (size_t i = 0; i < DIM; i ++) {
+ set_high_coordinate(i, DBL_MIN);
+ set_low_coordinate(i, DBL_MAX);
+ }
- double *coords_max = new double[DIM];
- double *coords_min = new double[DIM];
+ double smallest_distance;
+ double largest_distance;
+
+ for (int i = 0; i < n_cells_subsurface; i ++) {
- for (size_t i = 0; i < DIM; i ++)
- {
- coords_max[i] = DBL_MIN;
- coords_min[i] = DBL_MAX;
- }
+ int n1 = element_node_indices[i][1];
+ int n2 = element_node_indices[i][2];
+ int n3 = element_node_indices[i][3];
+ int n4 = element_node_indices[i][4];
+ int n5 = element_node_indices[i][5];
+ int n6 = element_node_indices[i][6];
- for (size_t i = 0; i < ny_e; i ++)
- {
+ std::vector<double> edge_length(N_VERTICES);
- cells[i].init_neighbors(SURFACE_N_EDGES);
+ for (size_t k = 0; k < N_VERTICES; k ++)
+ edge_length[k] = DBL_MAX;
+
+ for (size_t j = 0; j < DIM; j++) {
- 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;
- }
+ coords1[j] = nodes[n1][j];
+ coords2[j] = nodes[n2][j];
+ coords3[j] = nodes[n3][j];
+ coords4[j] = nodes[n4][j];
+ coords5[j] = nodes[n5][j];
+ coords6[j] = nodes[n6][j];
+
+ coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0;
+
+ // the characteristic length is the radius of the insphere of
+ // the wedge
+ 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]);
+ 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]);
+
+ 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]);
+ }
+
+ }
- 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];
- }
+ subsurface_cells[i].set_characteristic_length(smallest_distance);
+ subsurface_cells[i].set_maximum_length(largest_distance);
+
+ for (size_t k = 0; k < DIM; k ++) {
- cells[i].set_cell_id((int) i);
- cells[i].set_coordinates(coords_avg);
+ // update lowest and highest corner of mesh
+ high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]);
+ low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]);
+ subsurface_cells[i].set_coordinate(k, coords_avg[k]);
+
}
+
+ subsurface_cells[i].set_cell_id((int) i);
+
+ }
+
+}
- 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());
+void Mesh::get_mesh_metadata(H5::H5File fp, int flag)
+{
- #if VERBOSE
+ H5::DataSet dataset_e = fp.openDataSet(MESH_MIXED_ELEMENTS);
+ H5::DataType datatype_e = dataset_e.getDataType();
+ H5::DataSpace dataspace_e = dataset_e.getSpace();
- maxprint = 10;
+ int rank = dataspace_e.getSimpleExtentNdims();
- if (ny_e < maxprint)
- maxprint = ny_e;
+ hsize_t *dims_e = new hsize_t[rank];
+ dataspace_e.getSimpleExtentDims(dims_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());
- }
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
+#endif
- printf("[ok] ...\n");
-
- #endif
+ int *data_e = new int[dims_e[0]*dims_e[1]];
+ dataset_e.read(data_e, H5::PredType::STD_I32LE);
- delete[] coords_avg;
- delete[] coords1;
- delete[] coords2;
- delete[] coords3;
+ size_t nx_e = flag;
+ size_t ny_e = dims_e[0] * dims_e[1] / nx_e;
- printf("[OK] } done\n");
+ element_node_indices.resize(ny_e);
- printf("[OK] initialize edges {\n");
-
- int ns[N_EDGES];
- hsize_t **szudzik_ids = new hsize_t*[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];
+ }
- for (size_t i = 0; i < ny_e; i ++)
- szudzik_ids[i] = new hsize_t[N_EDGES];
+ 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);
+ }
- for (size_t i = 0; i < ny_e; i ++)
- {
- for (size_t j = 0; j < N_EDGES; j ++)
- ns[j] = md_e[i][j+1];
-
- for (size_t j = 0; j < N_EDGES - 1; j ++)
- szudzik_ids[i][j] = szudzik(ns[j], ns[j+1]);
-
- szudzik_ids[i][N_EDGES-1] = szudzik(ns[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 < N_EDGES; j ++)
- {
- printf(" %6lld ", szudzik_ids[i][j]);
- }
- printf("\n");
- }
+void Mesh::read_velocity(H5::H5File fp, std::string groupname, std::string datasetname, int flag, int index)
+{
- #endif
+ 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();
- printf("[OK] } done\n");
-
- printf("[OK] initialize neighbors {\n");
-
- hsize_t e0;
- hsize_t e1;
-
- for (size_t i = 0; i < ny_e; i ++)
- {
- 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) {
- cells[i].set_neighbor(j, (hsize_t) k);
- }
-
- }
- }
- }
- }
+ int rank = dataspace.getSimpleExtentNdims();
- #if VERBOSE
+ hsize_t *dims = new hsize_t[rank];
+ dataspace.getSimpleExtentDims(dims);
- 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
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
- printf("[OK] } done\n");
-
- delete[] szudzik_ids;
- delete[] md_e;
- delete[] data_e;
- delete[] dims_e;
+ double *data = new double[dims[0]*dims[1]];
+ dataset.read(data, H5::PredType::IEEE_F64LE);
- delete[] md_n;
- delete[] data_n;
- delete[] dims_n;
-
- printf("[OK] mesh built\n\n");
+ 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]);
+ }
- return 0;
-
}
-int Mesh::read_surface_data(H5::H5File dp)
+void Mesh::read_water_depth(H5::H5File fp)
{
-
- printf("[OK] build surface environment\n");
-
- // read velocity x
-
- H5::Group group = dp.openGroup("surface-velocity.cell.0");
-
- printf("[OK] reading velocity x { \n");
- printf("[ok] # of objects: %lld\n", group.getNumObjs());
-
- H5::DataSet dataset = group.openDataSet(DATA_GROUP_NAME_H5);
- H5::DataType datatype = dataset.getDataType();
- H5::DataSpace dataspace = dataset.getSpace();
+ 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();
- printf("[ok] rank of data: %d\n", rank);
-
hsize_t *dims = new hsize_t[rank];
dataspace.getSimpleExtentDims(dims);
- for (int i = 0; i < (int) rank; i ++)
- printf("[ok] dims[%d] of data: %lld\n", i, dims[i]);
-
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
+
double *data = new double[dims[0]*dims[1]];
- double **md = new double*[dims[0]];
-
- for (size_t i = 0; i < dims[0]; i ++)
- md[i] = data + i * dims[1];
-
dataset.read(data, H5::PredType::IEEE_F64LE);
-
- #if VERBOSE
- size_t maxprint = (dims[0] > 10) ? 10 : dims[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", data[i]);
- #endif
-
- printf("[OK] } done\n");
- printf("[OK] set velocity x { \n");
-
- for (int i = 0; i < n_cells; i ++)
- cells[i].set_v(0, data[i]);
-
- #if VERBOSE
- maxprint = (dims[0] > 10) ? 10 : dims[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));
+
+ for (int i = 0; i < n_cells_surface; i ++) {
+ surface_cells[i].set_water_depth(data[i]);
+#if VERBOSE
+ if (i%(n_cells_surface/5)==0) {
+ printf("[~~~~>] Surface ponded depth in cell[%d] = %e m\n", i,
+ surface_cells[i].get_water_depth());
}
+#endif
}
- #endif
- printf("[OK] } done\n");
-
- // read velocity y
-
- H5::Group group_vy = dp.openGroup("surface-velocity.cell.1");
-
- printf("[OK] reading velocity y { \n");
- printf("[ok] # of objects: %lld\n", group_vy.getNumObjs());
-
- H5::DataSet dataset_vy = group_vy.openDataSet(DATA_GROUP_NAME_H5);
- H5::DataType datatype_vy = dataset_vy.getDataType();
- H5::DataSpace dataspace_vy = dataset_vy.getSpace();
-
- rank = dataspace_vy.getSimpleExtentNdims();
-
- printf("[ok] rank of data: %d\n", rank);
-
- hsize_t *dims_vy = new hsize_t[rank];
- dataspace_vy.getSimpleExtentDims(dims_vy);
-
- for (int i = 0; i < (int) rank; i ++)
- printf("[ok] dims[%d] of data: %lld\n", i, dims_vy[i]);
-
- double *data_vy = new double[dims_vy[0]*dims_vy[1]];
- double **md_vy = new double*[dims_vy[0]];
+}
- for (size_t i = 0; i < dims_vy[0]; i ++)
- md_vy[i] = data_vy + i * dims_vy[1];
+void Mesh::read_elevation(H5::H5File fp)
+{
- dataset_vy.read(data_vy, H5::PredType::IEEE_F64LE);
+ 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();
- #if VERBOSE
- maxprint = (dims_vy[0] > 10) ? 10 : dims_vy[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", data_vy[i]);
- #endif
-
- printf("[OK] } done\n");
- printf("[OK] set velocity y { \n");
-
- for (int i = 0; i < n_cells; i ++)
- cells[i].set_v(1, data_vy[i]);
-
- #if VERBOSE
- maxprint = (dims_vy[0] > 10) ? 10 : dims_vy[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", cells[i].get_v(1));
- #endif
-
- printf("[OK] } done\n");
-
- // read water depth
- H5::Group group_h = dp.openGroup("surface-ponded_depth.cell.0");
-
- printf("[OK] reading ponded water depth { \n");
- printf("[ok] # of objects: %lld\n", group_vy.getNumObjs());
-
- H5::DataSet dataset_h = group_h.openDataSet(DATA_GROUP_NAME_H5);
- H5::DataType datatype_h = dataset_h.getDataType();
- H5::DataSpace dataspace_h = dataset_h.getSpace();
-
- rank = dataspace_h.getSimpleExtentNdims();
-
- printf("[ok] rank of data: %d\n", rank);
-
- hsize_t *dims_h = new hsize_t[rank];
- dataspace_h.getSimpleExtentDims(dims_h);
-
- for (int i = 0; i < (int) rank; i ++)
- printf("[ok] dims[%d] of data: %lld\n", i, dims_h[i]);
-
- double *data_h = new double[dims_h[0]*dims_h[1]];
- double **md_h = new double*[dims_h[0]];
+ int rank = dataspace.getSimpleExtentNdims();
- for (size_t i = 0; i < dims_h[0]; i ++)
- md_h[i] = data_h + i * dims_h[1];
+ hsize_t *dims = new hsize_t[rank];
+ dataspace.getSimpleExtentDims(dims);
- dataset_h.read(data_h, H5::PredType::IEEE_F64LE);
-
- #if VERBOSE
- maxprint = (dims_h[0] > 10) ? 10 : dims_h[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", data_h[i]);
- #endif
-
- printf("[OK] } done\n");
- printf("[OK] set ponded water depth { \n");
-
- for (int i = 0; i < n_cells; i ++)
- cells[i].set_h(data_h[i]);
-
- #if VERBOSE
- maxprint = (dims_h[0] > 10) ? 10 : dims_h[0];
- for (size_t i = 0; i < maxprint; i ++)
- printf("[ok] %f\n", cells[i].get_h());
- #endif
-
- printf("[OK] } done\n");
-
- delete [] md_h;
- delete [] data_h;
- delete [] md_vy;
- delete [] data_vy;
- delete [] md;
- delete [] data;
- delete [] dims;
+#if VERBOSE
+ for (int i = 0; i < rank; i ++)
+ std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
+#endif
- printf("[OK] data read\n\n");
+ double *data = new double[dims[0]*dims[1]];
+ dataset.read(data, H5::PredType::IEEE_F64LE);
- return 0;
+ for (int i = 0; i < n_cells_surface; i ++) {
+ surface_cells[i].set_coordinate(2, data[i]);
+#if VERBOSE
+ if (i%(n_cells_surface/5)==0) {
+ printf("[~~~~>] Surface elevation in cell[%d] = %e m\n", i,
+ surface_cells[i].get_coordinate(2));
+ }
+#endif
+ }
}
-
diff --git a/mesh.hpp b/mesh.hpp
index 008f302..70f7c34 100644
--- a/mesh.hpp
+++ b/mesh.hpp
@@ -1,115 +1,59 @@
-// -*- mode: c++; c-default-style: "linux" -*-
+// -*- mode: c++ -*-
#ifndef MESH_HPP
#define MESH_HPP
#include "header.hpp"
+#include "cell.hpp"
-// mesh.hpp
-// ========
-// Mesh topology to get environment variables that inform the
-// HydroAgent (agent.hpp).
-// Reads in a HDF5 mesh file and builds the topology. The HydroAgent
-// stores the cell it was last located in. Then, in the next iteration
-// it searches the neighborhood to locate itself. This makes the search
-// more efficient.
-
-class Cell
+class Mesh
{
public:
- Cell();
- Cell(double* coordinates);
- Cell(const Cell &cell);
-
- double get_x();
- double get_y();
- double get_z();
-
- double get_h();
- double get_v(int index);
-
- int get_cell_id();
-
- hsize_t get_neighbor(int index);
+ Mesh();
- int set_coordinates(double *coordinates);
-
- int set_neighbors(hsize_t *value_array);
- int set_neighbor(int index, hsize_t value);
+ void build_surface_mesh(H5::H5File fp, H5::H5File ep);
+ void build_subsurface_mesh(H5::H5File fp);
- int set_h(double value);
- int set_v(int index, double value);
+ void build_environment(H5::H5File ep, int flag);
+ void build_subsurface_environment(H5::H5File ep);
- int set_cell_id(int value);
+ Cell* get_cell(int index, int flag);
- int init_neighbors(size_t n_neighbors);
+ 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);
-protected:
+ int get_n_cells(int flag);
private:
-
- int dim; // n_neighbors = 2 * dim - 1
- int cell_id;
-
- double x;
- double y;
- double z;
-
- double h;
- double v[3];
-
- std::vector<hsize_t> neighbors;
-
-};
-
-class Mesh
-{
-
-public:
-
- Mesh();
- int get_n_cells();
+ int n_cells_subsurface;
+ int n_cells_surface;
- int build_surface(H5::H5File fp);
- int read_surface_data(H5::H5File dp);
-
-#if SUBSURFACE
- int build_subsurface(H5::H5File fp);
- int read_subsurface_data(H5::H5File dp);
-#endif
+ std::vector<Cell> subsurface_cells;
+ std::vector<Cell> surface_cells;
- Cell* get_cells();
+ std::vector<double> high_coordinates;
+ std::vector<double> low_coordinates;
- int set_min_coordinates(double *coordinates_min);
- int set_max_coordinates(double *coordinates_max);
-
- double get_xmax();
- double get_ymax();
- double get_zmax();
-
- double get_xmin();
- double get_ymin();
- double get_zmin();
-
-protected:
+ std::vector<std::vector<double>> nodes;
+ std::vector<std::vector<int>> element_node_indices;
+ //std::vector<std::vector<std::vector<hsize_t>>> edge_ids;
-private:
+ void handle_subsurface_node_data(H5::H5File fp);
+ void build_subsurface_cells();
- int n_cells;
- Cell *cells;
+ void handle_surface_node_data(H5::H5File fp);
+ void build_surface_cells(H5::H5File fp);
- double xmin;
- double xmax;
-
- double ymin;
- double ymax;
-
- double zmin;
- double zmax;
-
+ 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);
+ void read_elevation(H5::H5File fp);
+
};
#endif
diff --git a/new/agent.cpp b/new/agent.cpp
deleted file mode 100644
index 97d2e69..0000000
--- a/new/agent.cpp
+++ /dev/null
@@ -1,245 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "agent.hpp"
-
-// ------------------------------------------------------------------------
-// -- PUBLIC METHODS
-// ------------------------------------------------------------------------
-
-void Agent::advance(Mesh* mesh, size_t step)
-{
-
- switch(domain[step-1]) {
- case SURFACE_FLAG:
- advance_surface(mesh,step);
- break;
- case SUBSURFACE_FLAG:
- advance_subsurface(mesh,step);
- break;
- default:
- set_path(path[step-1][0], path[step-1][1], path[step-1][2], step);
- time[step] = time[step-1];
- domain[step] = OUTSIDE_FLAG;
- break;
- }
-
-}
-
-// ------------------------------------------------------------------------
-// GETTERS AND SETTERS
-// ------------------------------------------------------------------------
-
-Agent::Agent()
-{
- time.resize(N_ITERATION);
- domain.resize(N_ITERATION);
- path.resize(N_ITERATION);
- for (size_t i = 0; i < N_ITERATION; i ++)
- path[i].resize(DIM);
-}
-
-void Agent::set_cell_index(int value)
-{
- cell_index = value;
-}
-
-void Agent::set_domain(int value, size_t step)
-{
- domain[step] = value;
-}
-
-double Agent::get_time(size_t step)
-{
- return time[step];
-}
-
-void Agent::set_path(double x, double y, double z, size_t step)
-{
- path[step][0] = x;
- path[step][1] = y;
- path[step][2] = z;
-}
-
-double Agent::get_timestepsize()
-{
- return time_step_size;
-}
-
-double Agent::get_path(size_t step, size_t index)
-{
- return path[step][index];
-}
-
-int Agent::get_domain(size_t step)
-{
- return domain[step];
-}
-
-// ------------------------------------------------------------------------
-// -- PRIVATE METHODS
-// ------------------------------------------------------------------------
-
-void Agent::advance_surface(Mesh* mesh, size_t step)
-{
- Cell *location = mesh->get_cell(cell_index, SURFACE_FLAG);
-
- std::vector<double> velocity(DIM);
- std::vector<double> coordinates(DIM-1);
- std::vector<double> coordinates_neighbor(DIM-1);
- std::vector<double> distance(DIM-1);
-
- domain[step] = SURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
-
- for (size_t i = 0; i < DIM; i ++)
- velocity[i] = location->get_velocity(i);
-
- // CFL criteria
- double length = location->get_characteristic_length();
- double velocity_norm = sqrt(velocity[0]*velocity[0] + velocity[1]*velocity[1]);
-
- if (velocity_norm > EPS) {
- time_step_size = CFL * length / velocity_norm;
- } else {
- time_step_size = MAXTIMESTEP;
- }
-
- time[step] = time[step-1] + time_step_size;
-
- // check if agent infiltrates
- double h = location->get_water_depth();
- int flag = INFILTRATE_FALSE;
-
- if (h > EPS) { // only wet cell can infiltrate
- double qe = std::min(velocity[2], 0.0);
- double pi = fabs(qe) * time_step_size / h;
- double trial = (double) (rand()%100) / 100.0;
- if (trial < pi)
- flag = INFILTRATE_TRUE;
- }
-
- if (flag == INFILTRATE_FALSE) {
-
- // advect the agent horizontally, which means the last dimension
- // (z-direction) is omitted
- for (size_t i = 0; i < DIM-1; i ++) {
- coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
- distance[i] = location->get_coordinate(i) - coordinates[i];
- }
-
- // horizontal distance to the centroid of current cell
- double distance_norm = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
-
- // check if cell has changed
- size_t n_cells = mesh->get_n_cells(SURFACE_FLAG);
- for (size_t i = 0; i < n_cells; i ++) {
-
- for (size_t j = 0; j < DIM-1; j ++) {
- coordinates_neighbor[j] = mesh->get_cell(i, SURFACE_FLAG)->get_coordinate(j);
- distance[j] = coordinates_neighbor[j] - coordinates[j];
- }
-
- double distance_norm_ = sqrt(distance[0]*distance[0]+distance[1]*distance[1]);
-
- if (distance_norm > distance_norm_) {
- cell_index = i;
- distance_norm = distance_norm_;
- }
-
- }
-
- // check if out of the domain: if the agent travelled thrice the
- // maximum length of the cell and still can't find a neighboring
- // cell, it must be advected out of the domain
- if (distance_norm > 10.0 * location->get_maximum_length())
- domain[step] = OUTSIDE_FLAG;
-
- } else {
- // infiltration detected
- domain[step] = SUBSURFACE_FLAG;
- coordinates[0] = path[step-1][0];
- coordinates[1] = path[step-1][1];
- cell_index = mesh->get_cell(cell_index, SURFACE_FLAG)->get_coupled_neighbor();
- }
-
- set_path(coordinates[0], coordinates[1], location->get_coordinate(2), step);
-
-}
-
-
-void Agent::advance_subsurface(Mesh* mesh, size_t step)
-{
- Cell *location = mesh->get_cell(cell_index, SUBSURFACE_FLAG);
-
- std::vector<double> velocity(DIM);
- std::vector<double> coordinates(DIM);
- std::vector<double> coordinates_neighbor(DIM);
- std::vector<double> distance(DIM);
-
- domain[step] = SUBSURFACE_FLAG; // set domain for this step, may change later if agent infiltrates
-
- for (size_t i = 0; i < DIM; i ++)
- velocity[i] = location->get_velocity(i);
-
- // CFL criteria
- double length = location->get_characteristic_length();
- double velocity_norm = sqrt(velocity[0] * velocity[0] + velocity[1] * velocity[1] + velocity[2] * velocity[2]);
-
- if (velocity_norm > EPS) {
- time_step_size = CFL * length / velocity_norm;
- } else {
- time_step_size = MAXTIMESTEP;
- }
-
- time[step] = time[step-1] + time_step_size;
-
- // advect the agent
- for (size_t i = 0; i < DIM; i ++) {
- coordinates[i] = path[step-1][i] + time_step_size * velocity[i];
- distance[i] = location->get_coordinate(i) - coordinates[i];
- }
-
- set_path(coordinates[0], coordinates[1], coordinates[2], step);
-
- // distance to the centroid of current cell
- 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(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, SUBSURFACE_FLAG)->get_coordinate(j);
- distance[j] = coordinates_neighbor[j] - coordinates[j];
- }
-
- double distance_norm_ = sqrt(distance[0] * distance[0] + distance[1] * distance[1] + distance[2] * distance[2]);
-
- if (distance_norm > distance_norm_) {
- cell_index = i;
- distance_norm = distance_norm_;
- }
-
- }
-
- // check if out of the domain: if the agent travelled thrice the
- // maximum length of the cell and still can't find a neighboring
- // cell, it must be advected out of the domain
- if (distance_norm > 10.0 * location->get_maximum_length()) {
- domain[step] = OUTSIDE_FLAG;
- }
-
- // check if exfiltrates
- int neighbor_index = mesh->get_cell(cell_index, SUBSURFACE_FLAG)->get_coupled_neighbor();
- double z_surface;
- if (neighbor_index > -1) {
- // border to surface detected
- z_surface = mesh->get_cell(neighbor_index, SURFACE_FLAG)->get_coordinate(2);
- if (coordinates[2] >= z_surface) {
- // exfiltration detected
- cell_index = (int) neighbor_index;
- domain[step] = SURFACE_FLAG;
- set_path(coordinates[0], coordinates[1], z_surface, step);
- }
- }
-
-}
diff --git a/new/agent.hpp b/new/agent.hpp
deleted file mode 100644
index f99c2ab..0000000
--- a/new/agent.hpp
+++ /dev/null
@@ -1,41 +0,0 @@
-// -*- mode: c++ -*-
-
-#ifndef AGENT_HPP
-#define AGENT_HPP
-
-#include "header.hpp"
-#include "cell.hpp"
-#include "mesh.hpp"
-
-class Agent
-{
-
-public:
-
- Agent();
-
- void set_cell_index(int value);
- void set_domain(int value, size_t step);
- void set_path(double x, double y, double z, size_t step);
- double get_path(size_t step, size_t index);
- int get_domain(size_t step);
- double get_timestepsize();
- double get_time(size_t step);
-
- void advance(Mesh* mesh, size_t step);
-
-private:
-
- int cell_index;
- double time_step_size;
-
- std::vector<int> domain;
- std::vector<double> time;
- std::vector<std::vector<double>> path;
-
- void advance_surface(Mesh *mesh, size_t step);
- void advance_subsurface(Mesh *mesh, size_t step);
-
-};
-
-#endif
diff --git a/new/header.hpp b/new/header.hpp
deleted file mode 100644
index cff75d2..0000000
--- a/new/header.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-// -*- mode:c++ -*-
-
-#ifndef CONSTANTS_HPP
-#define CONSTANTS_HPP
-
-#include <iostream>
-#include <vector>
-#include <string> // string class
-//#include <omp.h>
-#include <algorithm> // min
-#include <float.h> // DBL_MAX
-#include <math.h> // fabs, sqrt
-#include <mpi.h>
-
-#include <H5Cpp.h>
-#include <H5File.h>
-
-#define MAXPRINT 5
-#define MAXTIMESTEP 1.0e5
-#define EPS 1.0e-10
-
-#define N_AGENT 100
-#define N_ITERATION 2000
-#define SEED 108
-#define CFL 0.5
-
-#define DIM 3
-#define N_NEIGHBORS 5
-#define N_VERTICES 6
-#define N_EDGES 9
-
-#define SURFACE_NCOL 4
-#define SUBSURFACE_NCOL 7
-
-#define INFILTRATE_TRUE 1
-#define INFILTRATE_FALSE 2
-
-#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
-#define SURFACE_FLAG 2
-#define OUTSIDE_FLAG 3
-
-#endif
diff --git a/new/main.cpp b/new/main.cpp
deleted file mode 100644
index 1e45aac..0000000
--- a/new/main.cpp
+++ /dev/null
@@ -1,133 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "mesh.hpp"
-#include "agent.hpp"
-
-int main(int argc, char* argv[])
-{
-
- int rank;
- int size;
-
- MPI_Init(&argc, &argv);
- MPI_Comm_size(MPI_COMM_WORLD, &size);
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-
- Mesh mesh = Mesh();
-
- H5::H5File fp_subsurface = H5::H5File(SUBSURFACE_MESH_FILE, H5F_ACC_RDONLY);
- H5::H5File dp_subsurface = H5::H5File(SUBSURFACE_DATA_FILE, H5F_ACC_RDONLY);
- H5::H5File fp_surface = H5::H5File(SURFACE_MESH_FILE, H5F_ACC_RDONLY);
- H5::H5File dp_surface = H5::H5File(SURFACE_DATA_FILE, H5F_ACC_RDONLY);
-
- if (rank == 0)
- std::cout << "[~>] Building the subsurface mesh" << std::endl;
-
- mesh.build_subsurface_mesh(fp_subsurface);
- fp_subsurface.close();
-
- if (rank == 0) {
- std::cout << "[~~>] Number of cells: " << mesh.get_n_cells(SUBSURFACE_FLAG) << std::endl;
- std::cout << "[~>] Building the surface mesh (includes coupling with subsurface and may take a while)" << std::endl;
- }
-
- mesh.build_surface_mesh(fp_surface, dp_surface);
- fp_surface.close();
-
- if (rank == 0) {
-
- int len_surf = mesh.get_n_cells(SURFACE_FLAG);
-
- std::cout << "[~~>] Number of cells: " << len_surf << std::endl;
-
- 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 < len_surf; n ++) {
- c = mesh.get_cell(n, SURFACE_FLAG);
- fprintf(mesh_fp, "%f, %f, %f, %d\n",
- c->get_coordinate(0),
- c->get_coordinate(1),
- c->get_coordinate(2), n);
- }
-
- fclose(mesh_fp);
-
- std::cout << "[~>] Reading the subsurface environment" << std::endl;
-
- }
-
- mesh.build_environment(dp_subsurface, SUBSURFACE_FLAG);
- dp_subsurface.close();
-
- mesh.build_environment(dp_surface, SURFACE_FLAG);
- dp_subsurface.close();
-
- // Place agents in the domain and advance N_ITERATION steps
-
- // Bootstrap
- srand(SEED + rank);
- int len = mesh.get_n_cells(SUBSURFACE_FLAG);
- int index;
-
- std::cout << "[~>] Releasing " << N_AGENT << " agents" << std::endl;
-
- size_t partition = N_AGENT / size;
- size_t low_bound = rank * partition;
- size_t high_bound = (rank + 1) * partition;
- if (rank == size - 1) // top processor picks up extras
- high_bound = N_AGENT;
-
- for (size_t i = low_bound; i < high_bound; i ++) {
-
- Agent agent;
-
- index = rand()%len;
- agent.set_cell_index(index);
- agent.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);
-
- agent.set_domain(SUBSURFACE_FLAG, 0);
- agent.set_domain(SUBSURFACE_FLAG, 1);
-
- agent.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;
-
- char output_fname[1024];
-
- // move agent
- for (step = 1; step < N_ITERATION; step ++) {
- if (step % (N_ITERATION / 5) == 0 && i % (N_AGENT / 5) == 0) {
- std::cout << "[~~~>] agent[" << i << "] step = " << step << " dt = " << agent.get_timestepsize() << " s" << std::endl;
- }
- agent.advance(&mesh, step);
- }
-
- // write out results
- sprintf(output_fname, "results/agent-%06lu.csv", i);
- FILE *result_fp = fopen(output_fname, "w");
-
- fprintf(result_fp, "time (s), x (m), y (m), z (m), domain\n");
- for (step = 0; step < N_ITERATION; step ++) {
- fprintf(result_fp, "%f, %f, %f, %f, %d\n",
- agent.get_time(step),
- agent.get_path(step, 0),
- agent.get_path(step, 1),
- agent.get_path(step, 2),
- agent.get_domain(step));
- }
- fclose(result_fp);
-
- }
-
- MPI_Finalize();
-
- return MPI_SUCCESS;
-
-}
diff --git a/new/makefile b/new/makefile
deleted file mode 100644
index 4a91e28..0000000
--- a/new/makefile
+++ /dev/null
@@ -1,33 +0,0 @@
-PRG=agent-provocateur
-
-#CC=g++-11
-CC=mpic++
-CXXFLAGS=-Wall -Wextra -pedantic -O0 -fopenmp
-LFLAGS=-I/opt/hdf5/include -L/opt/hdf5/lib
-SRC=cell.cpp mesh.cpp main.cpp agent.cpp
-LIB=-lm -lhdf5 -lhdf5_cpp
-
-OBJ=$(SRC:.cpp=.o)
-
-all:$(PRG)
-
-$(PRG): $(OBJ)
- echo "[OK] linking objects"
- $(CC) $(LFLAGS) $(CXXFLAGS) $(OBJ) -o $@ $(LIB)
-
-%.o: %.cpp
- echo "[OK] building objects"
- $(CC) $(LFLAGS) $(CXXFLAGS) -c $< -o $@
-
-run: $(PRG)
- ./$(PRG)
-
-clean:
- rm -f $(PRG) $(OBJ)
-
-cleanall:
- make clean
- rm -f *.dat; rm -f *.png
-
-.PHONY: all clean cleanall
-
diff --git a/new/mesh.cpp b/new/mesh.cpp
deleted file mode 100644
index 8112ede..0000000
--- a/new/mesh.cpp
+++ /dev/null
@@ -1,546 +0,0 @@
-// -*- mode: c++ -*-
-
-#include "mesh.hpp"
-
-
-// ------------------------------------------------------------------------
-// -- PUBLIC METHODS
-// ------------------------------------------------------------------------
-
-void Mesh::build_surface_mesh(H5::H5File fp, H5::H5File ep)
-{
-
- get_mesh_metadata(fp, SURFACE_NCOL);
-
- handle_surface_node_data(fp);
-
- build_surface_cells(ep);
-
-}
-
-void Mesh::build_subsurface_mesh(H5::H5File fp)
-{
-
- 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) {
- groupname = "darcy_velocity.cell.0";
- datasetname = SUBSURFACE_DATA_GROUP_NAME_H5;
- } else if (flag == SURFACE_FLAG) {
- 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 VERBOSE
- 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));
-#endif
- }
-
- // Read in water depth
- read_water_depth(ep);
-
- }
-#if VERBOSE
- 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));
- }
-
- }
-#endif
-
-}
-
-// ------------------------------------------------------------------------
-// GETTERS AND SETTERS
-// ------------------------------------------------------------------------
-
-Mesh::Mesh()
-{
- n_cells_surface = 0;
- n_cells_subsurface = 0;
- high_coordinates.resize(DIM);
- low_coordinates.resize(DIM);
-}
-
-int Mesh::get_n_cells(int flag)
-{
- int n_cells;
-
- if (flag == SUBSURFACE_FLAG) {
- n_cells = n_cells_subsurface;
- } else {
- n_cells = n_cells_surface;
- }
-
- return n_cells;
-}
-
-Cell* Mesh::get_cell(int index, int flag)
-{
-
- if (flag == SUBSURFACE_FLAG) {
- return &(subsurface_cells[index]);
- } else {
- return &(surface_cells[index]);
- }
-}
-
-void Mesh::set_high_coordinate(int index, double value)
-{
- high_coordinates[index] = value;
-}
-
-double Mesh::get_high_coordinate(int index)
-{
- return high_coordinates[index];
-}
-
-void Mesh::set_low_coordinate(int index, double value)
-{
- low_coordinates[index] = value;
-}
-
-double Mesh::get_low_coordinate(int index)
-{
- return low_coordinates[index];
-}
-
-// ------------------------------------------------------------------------
-// -- PRIVATE METHODS
-// ------------------------------------------------------------------------
-
-void Mesh::handle_surface_node_data(H5::H5File fp)
-{
-
- H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
- H5::DataType datatype_n = dataset_n.getDataType();
- H5::DataSpace dataspace_n = dataset_n.getSpace();
-
- int rank = dataspace_n.getSimpleExtentNdims();
-
- hsize_t *dims_n = new hsize_t[rank];
- dataspace_n.getSimpleExtentDims(dims_n);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
-#endif
-
- double *data_n = new double[dims_n[0]*dims_n[1]];
- dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
-
- nodes.resize(dims_n[0]);
-
- 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];
- }
-
-}
-
-void Mesh::build_surface_cells(H5::H5File ep)
-{
-
- 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 ++) {
-
- int n1 = element_node_indices[i][1];
- int n2 = element_node_indices[i][2];
- int n3 = element_node_indices[i][3];
-
- std::vector<double> edge_length(DIM);
-
- for (size_t k = 0; k < DIM; k ++)
- edge_length[k] = DBL_MAX;
-
-
- 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];
- coords3[j] = nodes[n3][j];
-
- coords_avg[j] = (coords1[j] + coords2[j] + coords3[j]) / 3.0;
-
- 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]);
-
- 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-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);
- 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) - 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 = k;
- if (dz <= 2.0 * subsurface_cells[k].get_characteristic_length())
- break;
- }
-
- }
-
- surface_cells[i].set_coupled_neighbor(coupled_neighbor);
- subsurface_cells[coupled_neighbor].set_coupled_neighbor(i);
-
-#if VERBOSE
- 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;
- }
-#endif
-
- }
-
-}
-
-void Mesh::handle_subsurface_node_data(H5::H5File fp)
-{
-
- H5::DataSet dataset_n = fp.openDataSet(MESH_NODES);
- H5::DataType datatype_n = dataset_n.getDataType();
- H5::DataSpace dataspace_n = dataset_n.getSpace();
-
- int rank = dataspace_n.getSimpleExtentNdims();
-
- hsize_t *dims_n = new hsize_t[rank];
- dataspace_n.getSimpleExtentDims(dims_n);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimension in " << i << " of nodes: " << dims_n[i] << std::endl;
-#endif
-
- double *data_n = new double[dims_n[0]*dims_n[1]];
- dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
-
- nodes.resize(dims_n[0]);
-
- 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];
- }
-
-}
-
-void Mesh::build_subsurface_cells()
-{
-
- std::vector<double> coords1(DIM);
- std::vector<double> coords2(DIM);
- std::vector<double> coords3(DIM);
- std::vector<double> coords4(DIM);
- std::vector<double> coords5(DIM);
- std::vector<double> coords6(DIM);
-
- std::vector<double> coords_avg(DIM);
-
- for (size_t i = 0; i < DIM; i ++) {
- set_high_coordinate(i, DBL_MIN);
- 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];
- int n2 = element_node_indices[i][2];
- int n3 = element_node_indices[i][3];
- int n4 = element_node_indices[i][4];
- int n5 = element_node_indices[i][5];
- int n6 = element_node_indices[i][6];
-
- std::vector<double> edge_length(N_VERTICES);
-
- for (size_t k = 0; k < N_VERTICES; k ++)
- edge_length[k] = DBL_MAX;
-
- for (size_t j = 0; j < DIM; j++) {
-
- coords1[j] = nodes[n1][j];
- coords2[j] = nodes[n2][j];
- coords3[j] = nodes[n3][j];
- coords4[j] = nodes[n4][j];
- coords5[j] = nodes[n5][j];
- coords6[j] = nodes[n6][j];
-
- coords_avg[j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0;
-
- // the characteristic length is the radius of the insphere of
- // the wedge
- 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]);
- 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]);
-
- 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);
- subsurface_cells[i].set_maximum_length(largest_distance);
-
- for (size_t k = 0; k < DIM; k ++) {
-
- // update lowest and highest corner of mesh
- high_coordinates[k] = std::max(coords_avg[k], high_coordinates[k]);
- low_coordinates[k] = std::min(coords_avg[k], low_coordinates[k]);
-
- subsurface_cells[i].set_coordinate(k, coords_avg[k]);
-
- }
-
- 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);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimension in " << i << " of element node map: " << dims_e[i] << std::endl;
-#endif
-
- 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);
- }
-
-}
-
-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);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
-#endif
-
- 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);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
-#endif
-
- 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 VERBOSE
- if (i%(n_cells_surface/5)==0) {
- printf("[~~~~>] Surface ponded depth in cell[%d] = %e m\n", i,
- surface_cells[i].get_water_depth());
- }
-#endif
- }
-
-}
-
-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);
-
-#if VERBOSE
- for (int i = 0; i < rank; i ++)
- std::cout << "[~~~>] Dimensions in " << i << " of velocity in x: " << dims[i] << std::endl;
-#endif
-
- 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 VERBOSE
- if (i%(n_cells_surface/5)==0) {
- printf("[~~~~>] Surface elevation in cell[%d] = %e m\n", i,
- surface_cells[i].get_coordinate(2));
- }
-#endif
- }
-
-}
diff --git a/new/mesh.hpp b/new/mesh.hpp
deleted file mode 100644
index 70f7c34..0000000
--- a/new/mesh.hpp
+++ /dev/null
@@ -1,59 +0,0 @@
-// -*- mode: c++ -*-
-
-#ifndef MESH_HPP
-#define MESH_HPP
-
-#include "header.hpp"
-#include "cell.hpp"
-
-class Mesh
-{
-
-public:
-
- Mesh();
-
- void build_surface_mesh(H5::H5File fp, H5::H5File ep);
- void build_subsurface_mesh(H5::H5File fp);
-
- void build_environment(H5::H5File ep, int flag);
- void build_subsurface_environment(H5::H5File ep);
-
- 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 flag);
-
-private:
-
- 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;
-
- void handle_subsurface_node_data(H5::H5File fp);
- void build_subsurface_cells();
-
- void handle_surface_node_data(H5::H5File fp);
- void build_surface_cells(H5::H5File fp);
-
- 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);
- void read_elevation(H5::H5File fp);
-
-};
-
-#endif
diff --git a/test.cpp b/test.cpp
deleted file mode 100755
index e31010a..0000000
--- a/test.cpp
+++ /dev/null
@@ -1,140 +0,0 @@
-#include <H5Cpp.h>
-#include <H5File.h>
-
-#include <iostream>
-
-#define MESH_FILE_NAME "test1/visdump_mesh.h5"
-#define DATA_FILE_NAME "test1/visdump_data.h5"
-
-int main(void)
-{
- std::cout << "Hello, world!" << std::endl;
-
- H5::H5File fp = H5::H5File(MESH_FILE_NAME, H5F_ACC_RDONLY);
- 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();
-
- std::cout << "Rank of element node map: " << rank << std::endl;
-
- hsize_t *dims_e = new hsize_t[rank];
- dataspace_e.getSimpleExtentDims(dims_e);
-
- for (int i = 0; i < (int) rank; i ++)
- std::cout << "dims[" << i << "] of element node map: " << dims_e[i] << std::endl;
-
- int *data_e = new int[dims_e[0]*dims_e[1]];
-
- size_t nx_e = 7;
- size_t ny_e = dims_e[0] * dims_e[1] / 7;
-
- int **md_e = new int*[ny_e];
-
- for (size_t i = 0; i < ny_e; i ++)
- md_e[i] = data_e + i * nx_e;
-
- dataset_e.read(data_e, H5::PredType::STD_I32LE);
-
- 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");
- }
-
- 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();
-
- hsize_t *dims_n = new hsize_t[rank];
- dataspace_n.getSimpleExtentDims(dims_n);
-
- for (int i = 0; i < (int) rank; i ++)
- std::cout << "dims[" << i << "] of nodes: " << dims_n[i] << std::endl;
-
- 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];
-
- dataset_n.read(data_n, H5::PredType::IEEE_F64LE);
-
- 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");
- }
-
- 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[3 * ny_e];
-
- for (size_t i = 0; i < ny_e; i ++) {
-
- 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];
- coords4[j] = md_n[n4][j];
- coords5[j] = md_n[n5][j];
- coords6[j] = md_n[n6][j];
-
- coords_avg[dims_n[1] * i + j] = (coords1[j] + coords2[j] + coords3[j] + coords4[j] + coords5[j] + coords6[j]) / 6.0;
-
- }
-
- }
-
-#define GNUPLOT "gnuplot -persist"
-
- FILE *gp;
- gp = popen(GNUPLOT, "w");
- fprintf(gp, "splot '-' w p lt -1 pt 6 title 'e and z'\n");
-
- for (size_t i = 0; i < ny_e; i ++)
- fprintf(gp, "%f %f %f\n", coords_avg[dims_n[1] * i], coords_avg[dims_n[1] * i + 1], coords_avg[dims_n[1] * i + 2]);
- fprintf(gp, "\033");
-
- fflush(gp);
- fclose(gp);
-
- fp.close();
-
- int ns[9];
- hsize_t **szudzik_ids = new hsize_t*[ny_e];
-
-
-
- return 0;
-}
diff --git a/tools.cpp b/tools.cpp
deleted file mode 100644
index 4a62950..0000000
--- a/tools.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-// -*- mode: c++; c-default-style: "linux" -*-
-
-#include "header.hpp"
-
-hsize_t szudzik(int n1, int n2)
-{
- hsize_t k1 = (hsize_t) ((n1 < n2) ? n1 : n2);
- hsize_t k2 = (hsize_t) ((n1 < n2) ? n2 : n1);
-
- return k2 * k2 + k1;
-}
diff --git a/tools.hpp b/tools.hpp
deleted file mode 100644
index 9bea0f0..0000000
--- a/tools.hpp
+++ /dev/null
@@ -1,8 +0,0 @@
-// -*- mode: c++; c-default-style: "linux" -*-
-
-#ifndef TOOLS_HPP
-#define TOOLS_HPP
-
-hsize_t szudzik(int n1, int n2);
-
-#endif