GRASS Programmer's Manual  6.4.2(2012)
N_geom.c
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * MODULE:       Grass PDE Numerical Library
00005 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
00006 *               soerengebbert <at> gmx <dot> de
00007 *               
00008 * PURPOSE:      part of the gpde library
00009 *               allocation, destroing and initializing the geometric struct
00010 *
00011 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *
00017 *****************************************************************************/
00018 
00019 
00020 #include "grass/N_pde.h"
00021 
00022 /* *************************************************************** * 
00023  * *********** Konstruktor *************************************** * 
00024  * *************************************************************** */
00030 N_geom_data *N_alloc_geom_data(void)
00031 {
00032     N_geom_data *geom = (N_geom_data *) G_calloc(1, sizeof(N_geom_data));
00033 
00034     geom->area = NULL;
00035     geom->planimetric = 1;
00036     geom->dim = 0;
00037 
00038     return geom;
00039 }
00040 
00041 /* *************************************************************** * 
00042  * *********** Destruktor **************************************** * 
00043  * *************************************************************** */
00050 void N_free_geom_data(N_geom_data * geom)
00051 {
00052     if (geom->area != NULL)
00053         G_free(geom->area);
00054 
00055     G_free(geom);
00056     return;
00057 }
00058 
00059 /* *************************************************************** * 
00060  * *************************************************************** * 
00061  * *************************************************************** */
00073 N_geom_data *N_init_geom_data_3d(G3D_Region * region3d, N_geom_data * geodata)
00074 {
00075     N_geom_data *geom = geodata;
00076     struct Cell_head region2d;
00077 
00078 #pragma omp critical
00079     {
00080 
00081         G_debug(2,
00082                 "N_init_geom_data_3d: initializing the geometry structure");
00083 
00084         if (geom == NULL)
00085             geom = N_alloc_geom_data();
00086 
00087         geom->dz = region3d->tb_res * G_database_units_to_meters_factor();      /*this function is not thread safe */
00088         geom->depths = region3d->depths;
00089         geom->dim = 3;
00090 
00091         /*convert the 3d into a 2d region and begin the area calculation */
00092         G_get_set_window(&region2d);    /*this function is not thread safe */
00093         G3d_regionToCellHead(region3d, &region2d);
00094     }
00095 
00096     return N_init_geom_data_2d(&region2d, geom);
00097 }
00098 
00099 
00100 /* *************************************************************** * 
00101  * *************************************************************** * 
00102  * *************************************************************** */
00114 N_geom_data *N_init_geom_data_2d(struct Cell_head * region,
00115                                  N_geom_data * geodata)
00116 {
00117     N_geom_data *geom = geodata;
00118     struct Cell_head backup;
00119     double meters;
00120     short ll = 0;
00121     int i;
00122 
00123 
00124     /*create an openmp lock to assure that only one thread at a time will access this function */
00125 #pragma omp critical
00126     {
00127         G_debug(2,
00128                 "N_init_geom_data_2d: initializing the geometry structure");
00129 
00130         /*make a backup from this region */
00131         G_get_set_window(&backup);      /*this function is not thread safe */
00132         /*set the current region */
00133         G_set_window(region);   /*this function is not thread safe */
00134 
00135         if (geom == NULL)
00136             geom = N_alloc_geom_data();
00137 
00138         meters = G_database_units_to_meters_factor();   /*this function is not thread safe */
00139 
00140         /*set the dim to 2d if it was not initiated with 3, thats a bit ugly :( */
00141         if (geom->dim != 3)
00142             geom->dim = 2;
00143 
00144         geom->planimetric = 1;
00145         geom->rows = region->rows;
00146         geom->cols = region->cols;
00147         geom->dx = region->ew_res * meters;
00148         geom->dy = region->ns_res * meters;
00149         geom->Az = geom->dy * geom->dx; /*square meters in planimetric proj */
00150         /*depths and dz are initialized with a 3d region */
00151 
00152         /*Begin the area calculation */
00153         ll = G_begin_cell_area_calculations();  /*this function is not thread safe */
00154 
00155         /*if the projection is not planimetric, calc the area for each row */
00156         if (ll == 2) {
00157             G_debug(2,
00158                     "N_init_geom_data_2d: calculating the areas for non parametric projection");
00159             geom->planimetric = 0;
00160 
00161             if (geom->area != NULL)
00162                 G_free(geom->area);
00163             else
00164                 geom->area = G_calloc(geom->rows, sizeof(double));
00165 
00166             /*fill the area vector */
00167             for (i = 0; i < geom->rows; i++) {
00168                 geom->area[i] = G_area_of_cell_at_row(i);       /*square meters */
00169             }
00170         }
00171 
00172         /*restore the old region */
00173         G_set_window(&backup);  /*this function is not thread safe */
00174     }
00175 
00176     return geom;
00177 }
00178 
00179 /* *************************************************************** * 
00180  * *************************************************************** * 
00181  * *************************************************************** */
00192 double N_get_geom_data_area_of_cell(N_geom_data * geom, int row)
00193 {
00194     if (geom->planimetric) {
00195         G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->Az);
00196         return geom->Az;
00197     }
00198     else {
00199         G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->area[row]);
00200         return geom->area[row];
00201     }
00202 
00203     return 0.0;
00204 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines