Loading...
Searching...
No Matches
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
*
* COPYRIGHT: (C) 2000 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
#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;
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
G_area_of_cell_at_row
double G_area_of_cell_at_row(int row)
Cell area in specified row.
Definition area.c:86
G_begin_cell_area_calculations
int G_begin_cell_area_calculations(void)
Begin cell area calculations.
Definition area.c:45
NULL
#define NULL
Definition ccmath.h:32
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
G_get_set_window
void G_get_set_window(struct Cell_head *window)
Get the current working window (region).
Definition gis/set_window.c:35
N_alloc_geom_data
N_geom_data * N_alloc_geom_data(void)
Allocate the pde geometry data structure and return a pointer to the new allocated structure.
Definition n_geom.c:29
N_get_geom_data_area_of_cell
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
Definition n_geom.c:196
N_init_geom_data_3d
N_geom_data * N_init_geom_data_3d(RASTER3D_Region *region3d, N_geom_data *geodata)
Initiate a pde geometry data structure with a 3d region.
Definition n_geom.c:73
N_init_geom_data_2d
N_geom_data * N_init_geom_data_2d(struct Cell_head *region, N_geom_data *geodata)
Initiate a pde geometry data structure with a 2d region.
Definition n_geom.c:117
N_free_geom_data
void N_free_geom_data(N_geom_data *geom)
Release memory of a pde geometry data structure.
Definition n_geom.c:49
G_database_units_to_meters_factor
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition proj3.c:146
N_geom_data
Geometric information about the structured grid.
Definition N_pde.h:101
N_geom_data::dx
double dx
Definition N_pde.h:108
N_geom_data::depths
int depths
Definition N_pde.h:114
N_geom_data::dy
double dy
Definition N_pde.h:109
N_geom_data::area
double * area
Definition N_pde.h:104
N_geom_data::dz
double dz
Definition N_pde.h:110
N_geom_data::dim
int dim
Definition N_pde.h:106
N_geom_data::rows
int rows
Definition N_pde.h:115
N_geom_data::planimetric
int planimetric
Definition N_pde.h:102
N_geom_data::cols
int cols
Definition N_pde.h:116
N_geom_data::Az
double Az
Definition N_pde.h:112