GRASS Programmer's Manual
6.4.2(2012)
|
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(®ion2d); /*this function is not thread safe */ 00093 G3d_regionToCellHead(region3d, ®ion2d); 00094 } 00095 00096 return N_init_geom_data_2d(®ion2d, 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 }