path: root/mesh.hpp

// -*- mode: c++ -*-

#ifndef MESH_HPP
#define MESH_HPP

#include "header.hpp"
#include "parser.hpp"
#include "subgrid.hpp"
#include "tooling.hpp"
#include "column.hpp"

// Content
// =======
//
// 1. Mesh
// ---
// 1.1 compile
// 1.2 dump
// 1.3 initialise_subgrid
// 1.4 map_surface
// 1.5 map_subsurface
// ---
// 1.6 get_bounding_box
// 1.7 read_node_indices
// 1.8 read_nodes
// 1.9 initialise_surface
// 1.10 initialise_subsurface
//

// ------------------------------------------------------------------------
// Mesh
// ------------------------------------------------------------------------
// The mesh is a collection of columns. The columns needn't to know
// their neighbours.
// ------------------------------------------------------------------------
class Mesh
{

public:

    std::vector<Column> columns;

    Mesh(Parser &value)
    {
        parser = value;
    }

    // ---
    // re-build the mesh and the subgrid from what little
    // information is contained in the ats output files
    // ---
    void compile()
    {

        std::cout << "[msh] compiling the mesh" << std::endl;

        H5::H5File fp_surface    = H5::H5File(parser.surface_mesh_file.c_str(),
                                              H5F_ACC_RDONLY);
        H5::H5File fp_subsurface = H5::H5File(parser.subsurface_mesh_file.c_str(),
                                              H5F_ACC_RDONLY);


        //

        std::cout << "[msh] surface mesh indexing" << std::endl;

        std::vector<std::vector<int>> node_indices;
        read_node_indices(fp_surface, node_indices, SURFACE_MESH_MIXED_ELEMENTS);
        n_vertices = node_indices[0].size();

        std::vector<point> nodes;
        read_nodes(fp_surface, nodes, SURFACE_MESH_NODES);

        initialise_surface(node_indices, nodes);

        //


        //

        std::cout << "[msh] subsurface mesh indexing" << std::endl;

        for (size_t i = 0; i < node_indices.size(); i ++)
            node_indices[i].clear();
        read_node_indices(fp_subsurface, node_indices, SUBSURFACE_MESH_MIXED_ELEMENTS);

        nodes.clear();
        read_nodes(fp_subsurface, nodes, SUBSURFACE_MESH_NODES);

        //


        std::cout << "[msh] bounding box: {";
        std::cout << low_coordinates[0] << ", " << low_coordinates[1];
        std::cout << "} {";
        std::cout << high_coordinates[0] << ", " << high_coordinates[1];
        std::cout << "}" << std::endl;


        //

        int nz = probe_column_depth(node_indices, nodes);
        initialise_subgrid(nz);

        map_surface();

        initialise_subsurface(node_indices, nodes);

        map_subsurface();

        std::cout << std::endl;

    }

    // ---
    // dump the subgrid with all topology and geometry information
    // ---
    void dump()
    {

        std::cout << "[msh] dumping the subgrid" << std::endl;

        int nx = sgrid.indices.size();
        int ny = sgrid.indices[0].size();
        int nz = sgrid.indices[0][0].size();

        int *subgrid_indices_array = new int[nx * ny * nz];
        double *subgrid_elevations_array = new double[nx * ny * nz];
        double *depths_array = new double[nz-1];

        for (int i = 0; i < nx; ++i) {
            for (int j = 0; j < ny; ++j) {
                for (int k = 0; k < nz; ++k) {
                    subgrid_indices_array[i + nx * (j + ny * k)] = sgrid.indices[i][j][k];
                    subgrid_elevations_array[i + nx * (j + ny * k)] = sgrid.elevations[i][j][k];
                }
            }
        }

        for (int i = 0; i < nz-1; ++i) {
            depths_array[i] = sgrid.dzs[i];
        }

        // open a hdf5 file and write stuff to file
        H5::H5File file("grid.h5", H5F_ACC_TRUNC);

        hsize_t rank = 1;
        hsize_t dims[1];

        // ---
        // write cellsize
        // ---
        double cellsize[1] = {parser.subgrid_dx};
        dims[0] = 1;

        H5::DataSpace dspace_cellsize(rank, dims);
        H5::FloatType dtype_cellsize(H5::PredType::NATIVE_DOUBLE);
        dtype_cellsize.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_cellsize = file.createDataSet("cellsize", dtype_cellsize, dspace_cellsize);
        dset_cellsize.write(cellsize, H5::PredType::NATIVE_DOUBLE);

        // ---
        // write the lower left corner
        // ---
        double upperleft[2] = {low_coordinates[0], high_coordinates[1]};
        dims[0] = 2;

        H5::DataSpace dspace_upperleft(rank, dims);
        H5::FloatType dtype_upperleft(H5::PredType::NATIVE_DOUBLE);
        dtype_upperleft.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_upperleft = file.createDataSet("upperleft", dtype_upperleft, dspace_upperleft);
        dset_upperleft.write(upperleft, H5::PredType::NATIVE_DOUBLE);

        // ---
        // write depths
        // ---
        dims[0] = nz-1;

        H5::DataSpace dspace_depths(rank, dims);
        H5::FloatType dtype_depths(H5::PredType::NATIVE_DOUBLE);
        dtype_depths.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_depths = file.createDataSet("depths", dtype_depths, dspace_depths);
        dset_depths.write(depths_array, H5::PredType::NATIVE_DOUBLE);

        // ---
        // write dimensions
        // ---
        int dimensions[3] = {nx, ny, nz};
        dims[0] = 3;

        H5::DataSpace dspace_dims(rank, dims);
        H5::IntType dtype_dims(H5::PredType::NATIVE_INT);
        dtype_dims.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_dims = file.createDataSet("dimensions", dtype_dims, dspace_dims);
        dset_dims.write(dimensions, H5::PredType::NATIVE_INT);

        dims[0] = nx * ny * nz;

        // ---
        // write indices
        // ---
        H5::DataSpace dspace_indices(rank, dims);
        H5::IntType dtype_indices(H5::PredType::NATIVE_INT);
        dtype_indices.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_indices = file.createDataSet("indices", dtype_indices, dspace_indices);
        dset_indices.write(subgrid_indices_array, H5::PredType::NATIVE_INT);

        // ---
        // write elevations
        // ---
        H5::DataSpace dspace_elevations(rank, dims);
        H5::FloatType dtype_elevations(H5::PredType::NATIVE_DOUBLE);
        dtype_elevations.setOrder(H5T_ORDER_LE);

        H5::DataSet dset_elevations = file.createDataSet("elevations", dtype_elevations, dspace_elevations);
        dset_elevations.write(subgrid_elevations_array, H5::PredType::NATIVE_DOUBLE);

        file.close();

        // free memory

        delete []subgrid_indices_array;
        delete []subgrid_elevations_array;
        delete []depths_array;

        std::cout << "[msh] subgrid dumped to grid.h5" << std::endl;

    }

    // ---
    // initialise the subgrid, compute its extent and how many grid
    // points are needed
    // ---
    void initialise_subgrid(int nz)
    {

        std::vector<double> depths(nz);
        std::cout << "[msh] initialising the subgrid" << std::endl;

        for (int i = 0; i < nz; i ++)
            depths[i] = columns[0].stack_upper_horizon[i] - columns[0].stack_lower_horizon[i];

        sgrid = Subgrid(parser, low_coordinates, high_coordinates, depths);

    }

    void get_bounding_box(std::vector<point> &vertices, double &xmin, double &xmax, double &ymin, double &ymax)
    {

        xmin = DBL_MAX;
        ymin = DBL_MAX;
        xmax = DBL_MIN;
        ymax = DBL_MIN;

        for (size_t k = 0; k < vertices.size(); k ++) {

            xmin = std::min(vertices[k][0], xmin);
            xmax = std::max(vertices[k][0], xmax);

            ymin = std::min(vertices[k][1], ymin);
            ymax = std::max(vertices[k][1], ymax);

        }

    }

    // ---
    // assign each subgrid subsurface point the index of the
    // subsurface cell containing it
    // ---
    void map_subsurface()
    {

        size_t i;
        std::cout << "[msh] mapping subsurface cells to subgrid...";

        double offset_x = low_coordinates[0];
        double offset_y = high_coordinates[1];
        double gridsize = parser.subgrid_dx;

        for (i = 0; i < columns.size(); ++ i) {

            // get bounding box
            double xmin = DBL_MAX;
            double ymin = DBL_MAX;
            double xmax = DBL_MIN;
            double ymax = DBL_MIN;

            get_bounding_box(columns[i].vertices, xmin, xmax, ymin, ymax);

            int colmin = (int) floor((xmin - offset_x) / gridsize);
            int rowmin = (int) floor((offset_y - ymax) / gridsize);
            int colmax = (int) floor((xmax - offset_x) / gridsize);
            int rowmax = (int) floor((offset_y - ymin) / gridsize);

            // check limits
            colmax = std::min(colmax, sgrid.nx - 1);
            rowmax = std::min(rowmax, sgrid.ny - 1);

            for (int ix = colmin; ix <= colmax; ++ ix) {
                for (int iy = rowmin; iy <= rowmax; ++ iy) {

                    point p1 { offset_x + 0.5 * gridsize + ix * gridsize,
                               offset_y - 0.5 * gridsize - iy * gridsize, 0.0 };

                    if (is_inside2d(p1, columns[i].centroid, columns[i].vertices)) {

                        for (size_t iz = 1; iz < columns[i].stack_index.size(); ++ iz) {

                            // map subsurface
                            sgrid.indices[ix][iy][iz] = columns[i].stack_index[iz];
                            sgrid.elevations[ix][iy][iz] = 0.5 * (columns[i].stack_lower_horizon[iz] + columns[i].stack_upper_horizon[iz]);

                        }

                    }

                }
            }

        }

        std::cout << " done" << std::endl;

    }

    // ---
    // assign each subgrid surface point the index of the surface cell
    // containing it
    // ---
    void map_surface()
    {

        // the strategy is to assign the top layer of the subgrid to the
        // corresponding cell in a 2d manner.  the underlying cells can be
        // mapped using indices and offset computation.

        size_t i;
        std::cout << "[msh] mapping surface cells to subgrid...";

        // index [0,0] corresponds to upper left corner of the domain.
        // we compute the offset in x and y direction accordingly.
        double offset_x = low_coordinates[0];
        double offset_y = high_coordinates[1];
        double gridsize = parser.subgrid_dx;

        #pragma omp parallel
        {
            #pragma omp for private(i)
            for (i = 0; i < columns.size(); ++ i) {

                // get bounding box
                double xmin = DBL_MAX;
                double ymin = DBL_MAX;
                double xmax = DBL_MIN;
                double ymax = DBL_MIN;

                get_bounding_box(columns[i].vertices, xmin, xmax, ymin, ymax);

                int colmin = (int) floor((xmin - offset_x) / gridsize);
                int rowmin = (int) floor((offset_y - ymax) / gridsize);
                int colmax = (int) floor((xmax - offset_x) / gridsize);
                int rowmax = (int) floor((offset_y - ymin) / gridsize);

                // check limits
                colmax = std::min(colmax, sgrid.nx - 1);
                rowmax = std::min(rowmax, sgrid.ny - 1);

                for (int ix = colmin; ix <= colmax; ++ ix) {
                    for (int iy = rowmin; iy <= rowmax; ++ iy) {

                        point p1 { offset_x + 0.5 * gridsize + ix * gridsize, offset_y - 0.5 * gridsize - iy * gridsize, 0.0 };

                        if (is_inside2d(p1, columns[i].centroid, columns[i].vertices)) {

                            // map surface
                            sgrid.indices[ix][iy][0] = i;
                            sgrid.elevations[ix][iy][0] = columns[i].centroid[2];

                        }

                    }
                }

            }

        } // end of parallel region

        std::cout << " done" << std::endl;

    }

private:

    Parser parser;
    Subgrid sgrid;
    point high_coordinates;
    point low_coordinates;

    int n_columns;
    int n_vertices;

    // ---
    // read the node indices from hdf5 file
    // ---
    void read_node_indices(H5::H5File fp, std::vector<std::vector<int>> &node_indices, std::string node_index_str)
    {

        H5::DataSet dset = fp.openDataSet(node_index_str.c_str());
        H5::DataType dtype = dset.getDataType();
        H5::DataSpace dspace = dset.getSpace();

        int rank = dspace.getSimpleExtentNdims();
        hsize_t* dims = new hsize_t[rank];
        dspace.getSimpleExtentDims(dims);

        int n_dims = 1;
        for (int i = 0; i < rank; i++)
            n_dims *= dims[i];

        int* data = new int[n_dims];
        dset.read(data, H5::PredType::STD_I32LE);

        // columns and vertices change from surface/subsurface
        size_t ncols = data[0];
        size_t nvers = 0;
        if (ncols == 4) {
            nvers = 3;
        } else if (ncols == 8) {
            ncols = 7;
            nvers = 6;
        }
        size_t nrows = n_dims / ncols;
        //

        std::cout << "[msh] number of vertices per cell: " << nvers << std::endl;

        node_indices.resize(nrows);

        for (size_t i = 0; i < nrows; i ++) {
            node_indices[i].resize(nvers);
            for (size_t j = 0; j < nvers; j ++)
                node_indices[i][j] = data[i * ncols + j + 1];
        }

        n_columns = nrows;

    }

    // ---
    // read in the coordinates from the hdf5 file
    // ---
    void read_nodes(H5::H5File fp, std::vector<point> &nodes, std::string mesh_nodes_str)
    {

        H5::DataSet dset = fp.openDataSet(mesh_nodes_str.c_str());
        H5::DataType dtype = dset.getDataType();
        H5::DataSpace dspace = dset.getSpace();

        int rank = dspace.getSimpleExtentNdims();

        hsize_t *dims = new hsize_t[rank];
        dspace.getSimpleExtentDims(dims);

        int n_dims = 1;
        for (int i = 0; i < rank; i++)
            n_dims *= dims[i];

        double *data = new double[n_dims];
        dset.read(data, H5::PredType::IEEE_F64LE);

        nodes.resize(dims[0]);

        for (size_t i = 0; i < dims[0]; i ++) {
            for (size_t j = 0; j < dims[1]; j ++)
                nodes[i][j] = data[i * dims[1] + j];
        }

    }

    // ---
    // initialise the top cell (vertices and centroids)
    // ---
    void initialise_surface(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
    {

        columns.resize(node_indices.size());

        for (int i = 0; i < DIM; i ++) {
            high_coordinates[i] = DBL_MIN;
            low_coordinates[i] = DBL_MAX;
        }

        for (size_t i = 0; i < node_indices.size(); i ++) {

            columns[i].vertices.resize(n_vertices);

            for (int j = 0; j < n_vertices; j ++) {

                int index = node_indices[i][j];

                for (size_t k = 0; k < DIM; k ++) {

                    double val = nodes[index][k];

                    // bounding box for the entire domain
                    high_coordinates[k] = std::max(val, high_coordinates[k]);
                    low_coordinates[k]  = std::min(val, low_coordinates[k]);
                    //

                    columns[i].vertices[j][k] = val;
                    columns[i].centroid[k] += val / (double) n_vertices;

                }
            }
        }

    }

    int probe_column_depth(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
    {

        int maxstack = 0;

        for (size_t i = 0; i < node_indices.size(); ++ i) {

            point centroid;
            for (size_t j = 0; j < DIM; ++ j) {
                centroid[j] = 0.0;
            }

            size_t index_count = node_indices[i].size();

            for (size_t j = 0; j < index_count; ++ j) {
                int index = node_indices[i][j];
                for (size_t k = 0; k < DIM; ++ k) {
                    double val = nodes[index][k];
                    centroid[k] += val / (double) index_count;
                }
            }

            double distance = get_distance2d(columns[0].centroid, centroid);
            if (distance <= parser.fuzz) {

                columns[0].stack_index.push_back(i);

                // determine horizons
                point p = nodes[node_indices[i][0]];
                double zmax = p[2];
                double zmin;

                for (size_t k = 1; k < node_indices[i].size(); ++k) {
                    point q = nodes[node_indices[i][k]];
                    if (get_distance2d(p, q) < parser.fuzz) {
                        if (zmax < q[2]) {
                            zmin = zmax;
                            zmax = q[2];
                        } else {
                            zmin = q[2];
                        }
                        break;
                    }
                }

                columns[0].stack_upper_horizon.push_back(zmax);
                columns[0].stack_lower_horizon.push_back(zmin);

                maxstack += 1;
            }

        }

        columns[0].sort();

        return maxstack;

    }

    // ---
    // initialise the stack, calculating the elevations
    // and node indices in the stack
    // ---
    void initialise_subsurface(std::vector<std::vector<int>> &node_indices, std::vector<point> &nodes)
    {

        size_t i;

        double offset_x = low_coordinates[0];
        double offset_y = high_coordinates[1];
        double gridsize = parser.subgrid_dx;

        #pragma omp parallel
        {
            #pragma omp for private(i)
            for (i = 0; i < node_indices.size(); ++ i) {

                // compute the centroid
                point centroid = {0.0, 0.0, 0.0};
                size_t index_count = node_indices[i].size();

                // vertical bounds
                double zmax = DBL_MIN;
                double zmin = DBL_MAX;

                for (size_t j = 0; j < index_count; ++ j) {

                    int index = node_indices[i][j];

                    double ztmp = nodes[index][2];
                    zmax = std::max(zmax, ztmp);
                    zmin = std::min(zmin, ztmp);

                    for (size_t k = 0; k < DIM; ++ k) {
                        double val = nodes[index][k];
                        centroid[k] += val / (double) index_count;
                    }

                }

                //

                // convert into indices and add to stack of column `index`
                int ix = (int) floor((centroid[0] - offset_x) / gridsize);
                int iy = (int) floor((offset_y - centroid[1]) / gridsize);
                int index = sgrid.indices[ix][iy][0];

                if (index < 0) {

                    printf("[msh] process %d: cannot find a column for cell %ld (%d, %d), centroid {%.12f %.12f}\n",
                           omp_get_thread_num(),
                           i, ix, iy, centroid[0], centroid[1]);
                    exit(EXIT_FAILURE);

                }

                if (index > 0) {
                    columns[index].stack_index.push_back(i);
                    columns[index].stack_upper_horizon.push_back(zmax);
                    columns[index].stack_lower_horizon.push_back(zmin);
                }
                //

            }

            // sanity check and sorting
            #pragma omp for private(i)
            for (i = 1; i < columns.size(); ++ i) {

                if ((int) columns[i].stack_index.size() != sgrid.nz) {
                    printf("[msh] process %d: incomplete column %ld with %ld/%d subsurface cells\n",
                           omp_get_thread_num(),
                           i,
                           columns[i].stack_index.size(), sgrid.nz);
                    exit(EXIT_FAILURE);
                }

                columns[i].sort();
            }


        }

    }

};

#endif