n_geom.c

Go to the documentation of this file.

/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      part of the gpde library
 *                 allocation, destroying and initializing the geometric struct
 *
 * SPDX-FileCopyrightText: 2000 GRASS Development Team
 * SPDX-License-Identifier: GPL-2.0-or-later
 *
 *****************************************************************************/
 
#include <grass/N_pde.h>
 
/* *************************************************************** *
 * *********** Konstruktor *************************************** *
 * *************************************************************** */
/*!
 * \brief Allocate the pde geometry data structure and return a pointer to the
 * new allocated structure
 *
 * \return N_geom_data *
 * */
N_geom_data *N_alloc_geom_data(void)
{
    N_geom_data *geom = (N_geom_data *)G_calloc(1, sizeof(N_geom_data));
 
    geom->area = NULL;
    geom->planimetric = 1;
    geom->dim = 0;
 
    return geom;
}
 
/* *************************************************************** *
 * *********** Destructor **************************************** *
 * *************************************************************** */
/*!
 * \brief Release memory of a pde geometry data structure
 *
 * \param  geom N_geom_data *
 * \return void
 * */
void N_free_geom_data(N_geom_data *geom)
{
    if (geom->area != NULL)
        G_free(geom->area);
 
    G_free(geom);
    return;
}
 
/* *************************************************************** *
 * *************************************************************** *
 * *************************************************************** */
/*!
 * \brief Initiate a pde geometry data structure with a 3d region
 *
 * If the projection is not planimetric, a double array will be created based on
 * the number of rows of the provided region
 *
 * \param region3d RASTER3D_Region *
 * \param geodata N_geom_data * - if a NULL pointer is given, a new structure
 * will be allocatet and returned
 *
 * \return N_geom_data *
 * */
N_geom_data *N_init_geom_data_3d(RASTER3D_Region *region3d,
                                 N_geom_data *geodata)
{
    N_geom_data *geom = geodata;
    struct Cell_head region2d;
 
#pragma omp critical
    {
 
        G_debug(2, "N_init_geom_data_3d: initializing the geometry structure");
 
        if (geom == NULL)
            geom = N_alloc_geom_data();
 
        geom->dz = region3d->tb_res *
                   G_database_units_to_meters_factor(); /*this function is not
                                                           thread safe */
        geom->depths = region3d->depths;
        geom->dim = 3;
 
        /*convert the 3d into a 2d region and begin the area calculation */
        G_get_set_window(&region2d); /*this function is not thread safe */
        Rast3d_region_to_cell_head(region3d, &region2d);
    }
 
    return N_init_geom_data_2d(&region2d, geom);
}
 
/* *************************************************************** *
 * *************************************************************** *
 * *************************************************************** */
/*!
 * \brief Initiate a pde geometry data structure with a 2d region
 *
 * If the projection is not planimetric, a double array will be created based on
 * the number of rows of the provided region storing all computed areas for each
 * row
 *
 * \param region struct Cell_head *
 * \param geodata N_geom_data * - if a NULL pointer is given, a new structure
 * will be allocatet and returned
 *
 * \return N_geom_data *
 * */
N_geom_data *N_init_geom_data_2d(struct Cell_head *region, N_geom_data *geodata)
{
    N_geom_data *geom = geodata;
    struct Cell_head backup;
    double meters;
    short ll = 0;
    int i;
 
    /*create an openmp lock to assure that only one thread at a time will access
     * this function */
#pragma omp critical
    {
        G_debug(2, "N_init_geom_data_2d: initializing the geometry structure");
 
        /*make a backup from this region */
        G_get_set_window(&backup); /*this function is not thread safe */
        /*set the current region */
        Rast_set_window(region); /*this function is not thread safe */
 
        if (geom == NULL)
            geom = N_alloc_geom_data();
 
        meters = G_database_units_to_meters_factor(); /*this function is not
                                                         thread safe */
 
        /*set the dim to 2d if it was not initiated with 3, that's a bit ugly :(
         */
        if (geom->dim != 3)
            geom->dim = 2;
 
        geom->planimetric = 1;
        geom->rows = region->rows;
        geom->cols = region->cols;
        geom->dx = region->ew_res * meters;
        geom->dy = region->ns_res * meters;
        geom->Az = geom->dy * geom->dx; /*square meters in planimetric proj */
        /*depths and dz are initialized with a 3d region */
 
        /*Begin the area calculation */
        ll = G_begin_cell_area_calculations(); /*this function is not thread
                                                  safe */
 
        /*if the projection is not planimetric, calc the area for each row */
        if (ll == 2) {
            G_debug(2, "N_init_geom_data_2d: calculating the areas for non "
                       "parametric projection");
            geom->planimetric = 0;
 
            if (geom->area != NULL)
                G_free(geom->area);
            else
                geom->area = G_calloc(geom->rows, sizeof(double));
 
            /*fill the area vector */
            for (i = 0; i < geom->rows; i++) {
                geom->area[i] = G_area_of_cell_at_row(i); /*square meters */
            }
        }
 
        /*restore the old region */
        Rast_set_window(&backup); /*this function is not thread safe */
    }
 
    return geom;
}
 
/* *************************************************************** *
 * *************************************************************** *
 * *************************************************************** */
/*!
 * \brief Get the areay size in square meter of one cell (x*y) at row
 *
 * This function works for two and three dimensions
 *
 * \param geom N_geom_data *
 * \param row int
 * \return area double
 *
 * */
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
{
    if (geom->planimetric) {
        G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->Az);
        return geom->Az;
    }
    else {
        G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->area[row]);
        return geom->area[row];
    }
 
    return 0.0;
}