GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
N_geom.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: part of the gpde library
9 * allocation, destroing and initializing the geometric struct
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 
19 
20 #include "grass/N_pde.h"
21 
22 /* *************************************************************** *
23  * *********** Konstruktor *************************************** *
24  * *************************************************************** */
31 {
32  N_geom_data *geom = (N_geom_data *) G_calloc(1, sizeof(N_geom_data));
33 
34  geom->area = NULL;
35  geom->planimetric = 1;
36  geom->dim = 0;
37 
38  return geom;
39 }
40 
41 /* *************************************************************** *
42  * *********** Destruktor **************************************** *
43  * *************************************************************** */
51 {
52  if (geom->area != NULL)
53  G_free(geom->area);
54 
55  G_free(geom);
56  return;
57 }
58 
59 /* *************************************************************** *
60  * *************************************************************** *
61  * *************************************************************** */
73 N_geom_data *N_init_geom_data_3d(G3D_Region * region3d, N_geom_data * geodata)
74 {
75  N_geom_data *geom = geodata;
76  struct Cell_head region2d;
77 
78 #pragma omp critical
79  {
80 
81  G_debug(2,
82  "N_init_geom_data_3d: initializing the geometry structure");
83 
84  if (geom == NULL)
85  geom = N_alloc_geom_data();
86 
87  geom->dz = region3d->tb_res * G_database_units_to_meters_factor(); /*this function is not thread safe */
88  geom->depths = region3d->depths;
89  geom->dim = 3;
90 
91  /*convert the 3d into a 2d region and begin the area calculation */
92  G_get_set_window(&region2d); /*this function is not thread safe */
93  G3d_regionToCellHead(region3d, &region2d);
94  }
95 
96  return N_init_geom_data_2d(&region2d, geom);
97 }
98 
99 
100 /* *************************************************************** *
101  * *************************************************************** *
102  * *************************************************************** */
114 N_geom_data *N_init_geom_data_2d(struct Cell_head * region,
115  N_geom_data * geodata)
116 {
117  N_geom_data *geom = geodata;
118  struct Cell_head backup;
119  double meters;
120  short ll = 0;
121  int i;
122 
123 
124  /*create an openmp lock to assure that only one thread at a time will access this function */
125 #pragma omp critical
126  {
127  G_debug(2,
128  "N_init_geom_data_2d: initializing the geometry structure");
129 
130  /*make a backup from this region */
131  G_get_set_window(&backup); /*this function is not thread safe */
132  /*set the current region */
133  G_set_window(region); /*this function is not thread safe */
134 
135  if (geom == NULL)
136  geom = N_alloc_geom_data();
137 
138  meters = G_database_units_to_meters_factor(); /*this function is not thread safe */
139 
140  /*set the dim to 2d if it was not initiated with 3, thats a bit ugly :( */
141  if (geom->dim != 3)
142  geom->dim = 2;
143 
144  geom->planimetric = 1;
145  geom->rows = region->rows;
146  geom->cols = region->cols;
147  geom->dx = region->ew_res * meters;
148  geom->dy = region->ns_res * meters;
149  geom->Az = geom->dy * geom->dx; /*square meters in planimetric proj */
150  /*depths and dz are initialized with a 3d region */
151 
152  /*Begin the area calculation */
153  ll = G_begin_cell_area_calculations(); /*this function is not thread safe */
154 
155  /*if the projection is not planimetric, calc the area for each row */
156  if (ll == 2) {
157  G_debug(2,
158  "N_init_geom_data_2d: calculating the areas for non parametric projection");
159  geom->planimetric = 0;
160 
161  if (geom->area != NULL)
162  G_free(geom->area);
163  else
164  geom->area = G_calloc(geom->rows, sizeof(double));
165 
166  /*fill the area vector */
167  for (i = 0; i < geom->rows; i++) {
168  geom->area[i] = G_area_of_cell_at_row(i); /*square meters */
169  }
170  }
171 
172  /*restore the old region */
173  G_set_window(&backup); /*this function is not thread safe */
174  }
175 
176  return geom;
177 }
178 
179 /* *************************************************************** *
180  * *************************************************************** *
181  * *************************************************************** */
193 {
194  if (geom->planimetric) {
195  G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->Az);
196  return geom->Az;
197  }
198  else {
199  G_debug(6, "N_get_geom_data_area_of_cell: %g", geom->area[row]);
200  return geom->area[row];
201  }
202 
203  return 0.0;
204 }