GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
n_les_assemble.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: functions to assemble a linear equation system
9 * part of the gpde library
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 <math.h>
21 #include <grass/N_pde.h>
22 
23 /* local protos */
24 static int make_les_entry_2d(int i, int j, int offset_i, int offset_j,
25  int count, int pos, N_les * les,
26  G_math_spvector * spvect,
27  N_array_2d * cell_count, N_array_2d * status,
28  N_array_2d * start_val, double entry,
29  int cell_type);
30 
31 static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
32  int offset_k, int count, int pos, N_les * les,
33  G_math_spvector * spvect,
34  N_array_3d * cell_count, N_array_3d * status,
35  N_array_3d * start_val, double entry,
36  int cell_type);
37 
38 /* *************************************************************** *
39  * ********************** N_alloc_5star ************************** *
40  * *************************************************************** */
41 /*!
42  * \brief allocate a 5 point star data structure
43  *
44  * \return N_data_star *
45  * */
47 {
48  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
49 
50  star->type = N_5_POINT_STAR;
51  star->count = 5;
52  return star;
53 }
54 
55 /* *************************************************************** *
56  * ********************* N_alloc_7star *************************** *
57  * *************************************************************** */
58 /*!
59  * \brief allocate a 7 point star data structure
60  *
61  * \return N_data_star *
62  * */
64 {
65  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
66 
67  star->type = N_7_POINT_STAR;
68  star->count = 7;
69  return star;
70 }
71 
72 /* *************************************************************** *
73  * ********************* N_alloc_9star *************************** *
74  * *************************************************************** */
75 /*!
76  * \brief allocate a 9 point star data structure
77  *
78  * \return N_data_star *
79  *
80  * \attention The 9 point start is not yet implemented in the matrix assembling function
81  *
82  * */
84 {
85  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
86 
87  star->type = N_9_POINT_STAR;
88  star->count = 9;
89  return star;
90 }
91 
92 /* *************************************************************** *
93  * ********************* N_alloc_27star ************************** *
94  * *************************************************************** */
95 /*!
96  * \brief allocate a 27 point star data structure
97  *
98  * \return N_data_star *
99  *
100  * \attention The 27 point start is not yet implemented in the matrix assembling function
101  *
102  * */
104 {
105  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
106 
107  star->type = N_27_POINT_STAR;
108  star->count = 27;
109  return star;
110 }
111 
112 /* *************************************************************** *
113  * ********************** N_create_5star ************************* *
114  * *************************************************************** */
115 /*!
116  * \brief allocate and initialize a 5 point star data structure
117  *
118  * \param C double
119  * \param W double
120  * \param E double
121  * \param N double
122  * \param S double
123  * \param V double
124  * \return N_data_star *
125  * */
126 N_data_star *N_create_5star(double C, double W, double E, double N,
127  double S, double V)
128 {
129  N_data_star *star = N_alloc_5star();
130 
131  star->C = C;
132  star->W = W;
133  star->E = E;
134  star->N = N;
135  star->S = S;
136 
137  star->V = V;
138 
139  G_debug(5, "N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->W,
140  star->E, star->N, star->S, star->C, star->V);
141 
142  return star;
143 }
144 
145 /* *************************************************************** *
146  * ************************* N_create_7star ********************** *
147  * *************************************************************** */
148 /*!
149  * \brief allocate and initialize a 7 point star data structure
150  *
151  * \param C double
152  * \param W double
153  * \param E double
154  * \param N double
155  * \param S double
156  * \param T double
157  * \param B double
158  * \param V double
159  * \return N_data_star *
160  * */
161 N_data_star *N_create_7star(double C, double W, double E, double N,
162  double S, double T, double B, double V)
163 {
164  N_data_star *star = N_alloc_7star();
165 
166  star->C = C;
167  star->W = W;
168  star->E = E;
169  star->N = N;
170  star->S = S;
171 
172  star->T = T;
173  star->B = B;
174 
175  star->V = V;
176 
177  G_debug(5, "N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n",
178  star->W, star->E, star->N, star->S, star->T, star->B, star->C,
179  star->V);
180 
181  return star;
182 }
183 
184 /* *************************************************************** *
185  * ************************ N_create_9star *********************** *
186  * *************************************************************** */
187 /*!
188  * \brief allocate and initialize a 9 point star data structure
189  *
190  * \param C double
191  * \param W double
192  * \param E double
193  * \param N double
194  * \param S double
195  * \param NW double
196  * \param SW double
197  * \param NE double
198  * \param SE double
199  * \param V double
200  * \return N_data_star *
201  * */
202 N_data_star *N_create_9star(double C, double W, double E, double N,
203  double S, double NW, double SW, double NE,
204  double SE, double V)
205 {
206  N_data_star *star = N_alloc_9star();
207 
208  star->C = C;
209  star->W = W;
210  star->E = E;
211  star->N = N;
212  star->S = S;
213 
214  star->NW = NW;
215  star->SW = SW;
216  star->NE = NE;
217  star->SE = SE;
218 
219  star->V = V;
220 
221  G_debug(5,
222  "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
223  star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
224  star->SE, star->C, star->V);
225 
226  return star;
227 }
228 
229 /* *************************************************************** *
230  * ************************ N_create_27star *********************** *
231  * *************************************************************** */
232 /*!
233  * \brief allocate and initialize a 27 point star data structure
234  *
235  * \param C double
236  * \param W double
237  * \param E double
238  * \param N double
239  * \param S double
240  * \param NW double
241  * \param SW double
242  * \param NE double
243  * \param SE double
244  * \param T double
245  * \param W_T double
246  * \param E_T double
247  * \param N_T double
248  * \param S_T double
249  * \param NW_T double
250  * \param SW_T double
251  * \param NE_T double
252  * \param SE_T double
253  * \param B double
254  * \param W_B double
255  * \param E_B double
256  * \param N_B double
257  * \param S_B double
258  * \param NW_B double
259  * \param SW_B double
260  * \param NE_B double
261  * \param SE_B double
262  * \param V double
263  * \return N_data_star *
264  * */
265 N_data_star *N_create_27star(double C, double W, double E, double N, double S,
266  double NW, double SW, double NE, double SE,
267  double T, double W_T, double E_T, double N_T,
268  double S_T, double NW_T, double SW_T,
269  double NE_T, double SE_T, double B, double W_B,
270  double E_B, double N_B, double S_B, double NW_B,
271  double SW_B, double NE_B, double SE_B, double V)
272 {
273  N_data_star *star = N_alloc_27star();
274 
275  star->C = C;
276  star->W = W;
277  star->E = E;
278  star->N = N;
279  star->S = S;
280 
281  star->NW = NW;
282  star->SW = SW;
283  star->NE = NE;
284  star->SE = SE;
285 
286  star->T = T;
287  star->W_T = W_T;
288  star->E_T = E_T;
289  star->N_T = N_T;
290  star->S_T = S_T;
291 
292  star->NW_T = NW_T;
293  star->SW_T = SW_T;
294  star->NE_T = NE_T;
295  star->SE_T = SE_T;
296 
297  star->B = B;
298  star->W_B = W_B;
299  star->E_B = E_B;
300  star->N_B = N_B;
301  star->S_B = S_B;
302 
303  star->NW_B = NW_B;
304  star->SW_B = SW_B;
305  star->NE_B = NE_B;
306  star->SE_B = SE_B;
307 
308  star->V = V;
309 
310  G_debug(5,
311  "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
312  star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
313  star->SE, star->C, star->V);
314 
315  G_debug(5,
316  "N_create_27star: w_t %g e_t %g n_t %g s_t %g nw_t %g sw_t %g ne_t %g se_t %g t %g \n",
317  star->W_T, star->E_T, star->N_T, star->S_T, star->NW_T,
318  star->SW_T, star->NE_T, star->SE_T, star->T);
319 
320  G_debug(5,
321  "N_create_27star: w_b %g e_b %g n_b %g s_b %g nw_b %g sw_b %g ne_b %g se_B %g b %g\n",
322  star->W_B, star->E_B, star->N_B, star->S_B, star->NW_B,
323  star->SW_B, star->NE_B, star->SE_B, star->B);
324 
325 
326 
327  return star;
328 }
329 
330 
331 /* *************************************************************** *
332  * ****************** N_set_les_callback_3d_func ***************** *
333  * *************************************************************** */
334 /*!
335  * \brief Set the callback function which is called while assembling the les in 3d
336  *
337  * \param data N_les_callback_3d *
338  * \param callback_func_3d N_data_star *
339  * \return void
340  * */
341 void
343  N_data_star * (*callback_func_3d) ())
344 {
345  data->callback = callback_func_3d;
346 }
347 
348 /* *************************************************************** *
349  * *************** N_set_les_callback_2d_func ******************** *
350  * *************************************************************** */
351 /*!
352  * \brief Set the callback function which is called while assembling the les in 2d
353  *
354  * \param data N_les_callback_2d *
355  * \param callback_func_2d N_data_star *
356  * \return void
357  * */
358 void
360  N_data_star * (*callback_func_2d) ())
361 {
362  data->callback = callback_func_2d;
363 }
364 
365 /* *************************************************************** *
366  * ************** N_alloc_les_callback_3d ************************ *
367  * *************************************************************** */
368 /*!
369  * \brief Allocate the structure holding the callback function
370  *
371  * A template callback is set. Use N_set_les_callback_3d_func
372  * to set up a specific function.
373  *
374  * \return N_les_callback_3d *
375  * */
377 {
378  N_les_callback_3d *call;
379 
380  call = (N_les_callback_3d *) G_calloc(1, sizeof(N_les_callback_3d *));
382 
383  return call;
384 }
385 
386 /* *************************************************************** *
387  * *************** N_alloc_les_callback_2d *********************** *
388  * *************************************************************** */
389 /*!
390  * \brief Allocate the structure holding the callback function
391  *
392  * A template callback is set. Use N_set_les_callback_2d_func
393  * to set up a specific function.
394  *
395  * \return N_les_callback_2d *
396  * */
398 {
399  N_les_callback_2d *call;
400 
401  call = (N_les_callback_2d *) G_calloc(1, sizeof(N_les_callback_2d *));
403 
404  return call;
405 }
406 
407 /* *************************************************************** *
408  * ******************** N_callback_template_3d ******************* *
409  * *************************************************************** */
410 /*!
411  * \brief A callback template creates a 7 point star structure
412  *
413  * This is a template callback for mass balance calculation with 7 point stars
414  * based on 3d data (g3d).
415  *
416  * \param data void *
417  * \param geom N_geom_data *
418  * \param depth int
419  * \param row int
420  * \param col int
421  * \return N_data_star *
422  *
423  * */
424 N_data_star *N_callback_template_3d(void *data, N_geom_data * geom, int col,
425  int row, int depth)
426 {
427  N_data_star *star = N_alloc_7star();
428 
429  star->E = 1 / geom->dx;
430  star->W = 1 / geom->dx;
431  star->N = 1 / geom->dy;
432  star->S = 1 / geom->dy;
433  star->T = 1 / geom->dz;
434  star->B = 1 / geom->dz;
435  star->C = -1 * (2 / geom->dx + 2 / geom->dy + 2 / geom->dz);
436  star->V = -1;
437 
438  G_debug(5,
439  "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
440  star->W, star->E, star->N, star->S, star->T, star->B, star->C,
441  star->V);
442 
443 
444  return star;
445 }
446 
447 /* *************************************************************** *
448  * ********************* N_callback_template_2d ****************** *
449  * *************************************************************** */
450 /*!
451  * \brief A callback template creates a 9 point star structure
452  *
453  * This is a template callback for mass balance calculation with 9 point stars
454  * based on 2d data (raster).
455  *
456  * \param data void *
457  * \param geom N_geom_data *
458  * \param row int
459  * \param col int
460  * \return N_data_star *
461  *
462  * */
463 N_data_star *N_callback_template_2d(void *data, N_geom_data * geom, int col,
464  int row)
465 {
466  N_data_star *star = N_alloc_9star();
467 
468  star->E = 1 / geom->dx;
469  star->NE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
470  star->SE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
471  star->W = 1 / geom->dx;
472  star->NW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
473  star->SW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
474  star->N = 1 / geom->dy;
475  star->S = 1 / geom->dy;
476  star->C =
477  -1 * (star->E + star->NE + star->SE + star->W + star->NW + star->SW +
478  star->N + star->S);
479  star->V = 0;
480 
481  return star;
482 }
483 
484 /* *************************************************************** *
485  * ******************** N_assemble_les_2d ************************ *
486  * *************************************************************** */
487 /*!
488  * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
489  *
490  * This function calls #N_assemble_les_2d_param
491  *
492  */
493 N_les *N_assemble_les_2d(int les_type, N_geom_data * geom,
494  N_array_2d * status, N_array_2d * start_val,
495  void *data, N_les_callback_2d * call)
496 {
497  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
498  call, N_CELL_ACTIVE);
499 }
500 
501 /*!
502  * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
503  *
504  * This function calls #N_assemble_les_2d_param
505  *
506  */
508  N_array_2d * status, N_array_2d * start_val,
509  void *data, N_les_callback_2d * call)
510 {
511  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
512  call, N_CELL_ACTIVE);
513 }
514 
515 /*!
516  * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet cells
517  *
518  * This function calls #N_assemble_les_2d_param
519  *
520  */
522  N_array_2d * status,
523  N_array_2d * start_val, void *data,
524  N_les_callback_2d * call)
525 {
526  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
527  call, N_CELL_DIRICHLET);
528 }
529 
530 /*!
531  * \brief Assemble a linear equation system (les) based on 2d location data (raster)
532  *
533  *
534  * The linear equation system type can be set to N_NORMAL_LES to create a regular
535  * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
536  * a new created linear equation system which can be solved with
537  * linear equation solvers. An 2d array with start values and an 2d status array
538  * must be provided as well as the location geometry and a void pointer to data
539  * passed to the callback which creates the les row entries. This callback
540  * must be defined in the N_les_callback_2d strcuture.
541  *
542  * The creation of the les is parallelized with OpenMP.
543  * If you implement new callbacks, please make sure that the
544  * function calls are thread safe.
545  *
546  *
547  * the les can be created in two ways, with dirichlet and similar cells and without them,
548  * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
549  * must be added.
550  *
551  * \param les_type int
552  * \param geom N_geom_data*
553  * \param status N_array_2d *
554  * \param start_val N_array_2d *
555  * \param data void *
556  * \param cell_type int -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
557  * \param call N_les_callback_2d *
558  * \return N_les *
559  * */
561  N_array_2d * status, N_array_2d * start_val,
562  void *data, N_les_callback_2d * call,
563  int cell_type)
564 {
565  int i, j, count = 0, pos = 0;
566  int cell_type_count = 0;
567  int **index_ij;
568  N_array_2d *cell_count;
569  N_les *les = NULL;
570 
571  G_debug(2,
572  "N_assemble_les_2d: starting to assemble the linear equation system");
573 
574  /* At first count the number of valid cells and save
575  * each number in a new 2d array. Those numbers are used
576  * to create the linear equation system.
577  * */
578 
579  cell_count = N_alloc_array_2d(geom->cols, geom->rows, 1, CELL_TYPE);
580 
581  /* include dirichlet cells in the les */
582  if (cell_type == N_CELL_DIRICHLET) {
583  for (j = 0; j < geom->rows; j++) {
584  for (i = 0; i < geom->cols; i++) {
585  /*use all non-inactive cells for les creation */
586  if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
588  cell_type_count++;
589  }
590  }
591  }
592  /*use only active cell in the les */
593  if (cell_type == N_CELL_ACTIVE) {
594  for (j = 0; j < geom->rows; j++) {
595  for (i = 0; i < geom->cols; i++) {
596  /*count only active cells */
597  if (N_CELL_ACTIVE == N_get_array_2d_d_value(status, i, j))
598  cell_type_count++;
599  }
600  }
601  }
602 
603  G_debug(2, "N_assemble_les_2d: number of used cells %i\n",
604  cell_type_count);
605 
606  if (cell_type_count == 0)
608  ("Not enough cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
609  cell_type_count);
610 
611  /* Then allocate the memory for the linear equation system (les).
612  * Only valid cells are used to create the les. */
613  index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
614  for (i = 0; i < cell_type_count; i++)
615  index_ij[i] = (int *)G_calloc(2, sizeof(int));
616 
617  les = N_alloc_les_Ax_b(cell_type_count, les_type);
618 
619  count = 0;
620 
621  /*count the number of cells which should be used to create the linear equation system */
622  /*save the i and j indices and create a ordered numbering */
623  for (j = 0; j < geom->rows; j++) {
624  for (i = 0; i < geom->cols; i++) {
625  /*count every non-inactive cell */
626  if (cell_type == N_CELL_DIRICHLET) {
627  if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
628  N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE) {
629  N_put_array_2d_c_value(cell_count, i, j, count);
630  index_ij[count][0] = i;
631  index_ij[count][1] = j;
632  count++;
633  G_debug(5,
634  "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
635  count, i, j);
636  }
637  /*count every active cell */
638  }
639  else if (N_CELL_ACTIVE == N_get_array_2d_c_value(status, i, j)) {
640  N_put_array_2d_c_value(cell_count, i, j, count);
641  index_ij[count][0] = i;
642  index_ij[count][1] = j;
643  count++;
644  G_debug(5,
645  "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
646  count, i, j);
647  }
648  }
649  }
650 
651  G_debug(2, "N_assemble_les_2d: starting the parallel assemble loop");
652 
653  /* Assemble the matrix in parallel */
654 #pragma omp parallel for private(i, j, pos, count) schedule(static)
655  for (count = 0; count < cell_type_count; count++) {
656  i = index_ij[count][0];
657  j = index_ij[count][1];
658 
659  /*create the entries for the */
660  N_data_star *items = call->callback(data, geom, i, j);
661 
662  /* we need a sparse vector pointer anytime */
663  G_math_spvector *spvect = NULL;
664 
665  /*allocate a sprase vector */
666  if (les_type == N_SPARSE_LES) {
667  spvect = G_math_alloc_spvector(items->count);
668  }
669  /* initial conditions */
670  les->x[count] = N_get_array_2d_d_value(start_val, i, j);
671 
672  /* the entry in the vector b */
673  les->b[count] = items->V;
674 
675  /* pos describes the position in the sparse vector.
676  * the first entry is always the diagonal entry of the matrix*/
677  pos = 0;
678 
679  if (les_type == N_SPARSE_LES) {
680  spvect->index[pos] = count;
681  spvect->values[pos] = items->C;
682  }
683  else {
684  les->A[count][count] = items->C;
685  }
686  /* western neighbour, entry is col - 1 */
687  if (i > 0) {
688  pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
689  cell_count, status, start_val, items->W,
690  cell_type);
691  }
692  /* eastern neighbour, entry col + 1 */
693  if (i < geom->cols - 1) {
694  pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
695  cell_count, status, start_val, items->E,
696  cell_type);
697  }
698  /* northern neighbour, entry row - 1 */
699  if (j > 0) {
700  pos =
701  make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
702  cell_count, status, start_val, items->N,
703  cell_type);
704  }
705  /* southern neighbour, entry row + 1 */
706  if (j < geom->rows - 1) {
707  pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
708  cell_count, status, start_val, items->S,
709  cell_type);
710  }
711  /*in case of a nine point star, we have additional entries */
712  if (items->type == N_9_POINT_STAR) {
713  /* north-western neighbour, entry is col - 1 row - 1 */
714  if (i > 0 && j > 0) {
715  pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
716  cell_count, status, start_val,
717  items->NW, cell_type);
718  }
719  /* north-eastern neighbour, entry col + 1 row - 1 */
720  if (i < geom->cols - 1 && j > 0) {
721  pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
722  cell_count, status, start_val,
723  items->NE, cell_type);
724  }
725  /* south-western neighbour, entry is col - 1 row + 1 */
726  if (i > 0 && j < geom->rows - 1) {
727  pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect,
728  cell_count, status, start_val,
729  items->SW, cell_type);
730  }
731  /* south-eastern neighbour, entry col + 1 row + 1 */
732  if (i < geom->cols - 1 && j < geom->rows - 1) {
733  pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
734  cell_count, status, start_val,
735  items->SE, cell_type);
736  }
737  }
738 
739  /*How many entries in the les */
740  if (les->type == N_SPARSE_LES) {
741  spvect->cols = pos + 1;
742  G_math_add_spvector(les->Asp, spvect, count);
743  }
744 
745  if (items)
746  G_free(items);
747  }
748 
749  /*release memory */
750  N_free_array_2d(cell_count);
751 
752  for (i = 0; i < cell_type_count; i++)
753  G_free(index_ij[i]);
754 
755  G_free(index_ij);
756 
757  return les;
758 }
759 
760 /*!
761  * \brief Integrate Dirichlet or Transmission boundary conditions into the les (2s)
762  *
763  * Dirichlet and Transmission boundary conditions will be integrated into
764  * the provided linear equation system. This is meaningful if
765  * the les was created with #N_assemble_les_2d_dirichlet, because in
766  * this case Dirichlet boundary conditions are not automatically included.
767  *
768  * The provided les will be modified:
769  *
770  * Ax = b will be split into Ax_u + Ax_d = b
771  *
772  * x_u - the unknowns
773  * x_d - the Dirichlet cells
774  *
775  * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
776  *
777  * | A_u 0 | x_u
778  * | 0 I | x_d
779  *
780  * \param les N_les* -- the linear equation system
781  * \param geom N_geom_data* -- geometrical data information
782  * \param status N_array_2d* -- the status array containing the cell types
783  * \param start_val N_array_2d* -- an array with start values
784  * \return int -- 1 = success, 0 = failure
785  * */
787  N_array_2d * status, N_array_2d * start_val)
788 {
789  int rows, cols;
790  int count = 0;
791  int i, j, x, y, stat;
792  double *dvect1;
793  double *dvect2;
794 
795  G_debug(2,
796  "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition");
797 
798  rows = geom->rows;
799  cols = geom->cols;
800 
801  /*we nned to additional vectors */
802  dvect1 = (double *)G_calloc(les->cols, sizeof(double));
803  dvect2 = (double *)G_calloc(les->cols, sizeof(double));
804 
805  /*fill the first one with the x vector data of Dirichlet cells */
806  count = 0;
807  for (y = 0; y < rows; y++) {
808  for (x = 0; x < cols; x++) {
809  stat = N_get_array_2d_c_value(status, x, y);
810  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
811  dvect1[count] = N_get_array_2d_d_value(start_val, x, y);
812  count++;
813  }
814  else if (stat == N_CELL_ACTIVE) {
815  dvect1[count] = 0.0;
816  count++;
817  }
818  }
819  }
820 
821 #pragma omp parallel default(shared)
822  {
823  /*perform the matrix vector product and */
824  if (les->type == N_SPARSE_LES)
825  G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
826  else
827  G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
828 #pragma omp for schedule (static) private(i)
829  for (i = 0; i < les->cols; i++)
830  les->b[i] = les->b[i] - dvect2[i];
831  }
832 
833  /*now set the Dirichlet cell rows and cols to zero and the
834  * diagonal entry to 1*/
835  count = 0;
836  for (y = 0; y < rows; y++) {
837  for (x = 0; x < cols; x++) {
838  stat = N_get_array_2d_c_value(status, x, y);
839  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
840  if (les->type == N_SPARSE_LES) {
841  /*set the rows to zero */
842  for (i = 0; i < les->Asp[count]->cols; i++)
843  les->Asp[count]->values[i] = 0.0;
844  /*set the cols to zero */
845  for (i = 0; i < les->rows; i++) {
846  for (j = 0; j < les->Asp[i]->cols; j++) {
847  if (les->Asp[i]->index[j] == count)
848  les->Asp[i]->values[j] = 0.0;
849  }
850  }
851 
852  /*entry on the diagonal */
853  les->Asp[count]->values[0] = 1.0;
854 
855  }
856  else {
857  /*set the rows to zero */
858  for (i = 0; i < les->cols; i++)
859  les->A[count][i] = 0.0;
860  /*set the cols to zero */
861  for (i = 0; i < les->rows; i++)
862  les->A[i][count] = 0.0;
863 
864  /*entry on the diagonal */
865  les->A[count][count] = 1.0;
866  }
867  }
868  if (stat >= N_CELL_ACTIVE)
869  count++;
870  }
871  }
872 
873  return 0;
874 
875 }
876 
877 /* **************************************************************** */
878 /* **** make an entry in the les (2d) ***************************** */
879 /* **************************************************************** */
880 int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count,
881  int pos, N_les * les, G_math_spvector * spvect,
882  N_array_2d * cell_count, N_array_2d * status,
883  N_array_2d * start_val, double entry, int cell_type)
884 {
885  int K;
886  int di = offset_i;
887  int dj = offset_j;
888 
889  K = N_get_array_2d_c_value(cell_count, i + di, j + dj) -
890  N_get_array_2d_c_value(cell_count, i, j);
891 
892  /* active cells build the linear equation system */
893  if (cell_type == N_CELL_ACTIVE) {
894  /* dirichlet or transmission cells must be handled like this */
895  if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_ACTIVE &&
896  N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE)
897  les->b[count] -=
898  N_get_array_2d_d_value(start_val, i + di, j + dj) * entry;
899  else if (N_get_array_2d_c_value(status, i + di, j + dj) ==
900  N_CELL_ACTIVE) {
901  if ((count + K) >= 0 && (count + K) < les->cols) {
902  G_debug(5,
903  " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
904  count, count + K, entry);
905  pos++;
906  if (les->type == N_SPARSE_LES) {
907  spvect->index[pos] = count + K;
908  spvect->values[pos] = entry;
909  }
910  else {
911  les->A[count][count + K] = entry;
912  }
913  }
914  }
915  } /* if dirichlet cells should be used then check for all valid cell neighbours */
916  else if (cell_type == N_CELL_DIRICHLET) {
917  /* all valid cells */
918  if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_INACTIVE
919  && N_get_array_2d_c_value(status, i + di,
920  j + dj) < N_MAX_CELL_STATE) {
921  if ((count + K) >= 0 && (count + K) < les->cols) {
922  G_debug(5,
923  " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
924  count, count + K, entry);
925  pos++;
926  if (les->type == N_SPARSE_LES) {
927  spvect->index[pos] = count + K;
928  spvect->values[pos] = entry;
929  }
930  else {
931  les->A[count][count + K] = entry;
932  }
933  }
934  }
935  }
936 
937  return pos;
938 }
939 
940 
941 /* *************************************************************** *
942  * ******************** N_assemble_les_3d ************************ *
943  * *************************************************************** */
944 /*!
945  * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
946  *
947  * This function calls #N_assemble_les_3d_param
948  * */
949 N_les *N_assemble_les_3d(int les_type, N_geom_data * geom,
950  N_array_3d * status, N_array_3d * start_val,
951  void *data, N_les_callback_3d * call)
952 {
953  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
954  call, N_CELL_ACTIVE);
955 }
956 
957 /*!
958  * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
959  *
960  * This function calls #N_assemble_les_3d_param
961  * */
963  N_array_3d * status, N_array_3d * start_val,
964  void *data, N_les_callback_3d * call)
965 {
966  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
967  call, N_CELL_ACTIVE);
968 }
969 
970 /*!
971  * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells
972  *
973  * This function calls #N_assemble_les_3d_param
974  * */
976  N_array_3d * status,
977  N_array_3d * start_val, void *data,
978  N_les_callback_3d * call)
979 {
980  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
981  call, N_CELL_DIRICHLET);
982 }
983 
984 /*!
985  * \brief Assemble a linear equation system (les) based on 3d location data (g3d)
986  *
987  * The linear equation system type can be set to N_NORMAL_LES to create a regular
988  * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
989  * a new created linear equation system which can be solved with
990  * linear equation solvers. An 3d array with start values and an 3d status array
991  * must be provided as well as the location geometry and a void pointer to data
992  * passed to the callback which creates the les row entries. This callback
993  * must be defined in the N_les_callback_3d structure.
994  *
995  * The creation of the les is parallelized with OpenMP.
996  * If you implement new callbacks, please make sure that the
997  * function calls are thread safe.
998  *
999  * the les can be created in two ways, with dirichlet and similar cells and without them,
1000  * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
1001  * must be added.
1002  *
1003  * \param les_type int
1004  * \param geom N_geom_data*
1005  * \param status N_array_3d *
1006  * \param start_val N_array_3d *
1007  * \param data void *
1008  * \param call N_les_callback_3d *
1009  * \param cell_type int -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
1010  * \return N_les *
1011  * */
1013  N_array_3d * status, N_array_3d * start_val,
1014  void *data, N_les_callback_3d * call,
1015  int cell_type)
1016 {
1017  int i, j, k, count = 0, pos = 0;
1018  int cell_type_count = 0;
1019  N_array_3d *cell_count;
1020  N_les *les = NULL;
1021  int **index_ij;
1022 
1023  G_debug(2,
1024  "N_assemble_les_3d: starting to assemble the linear equation system");
1025 
1026  cell_count =
1027  N_alloc_array_3d(geom->cols, geom->rows, geom->depths, 1, DCELL_TYPE);
1028 
1029  /* First count the number of valid cells and save
1030  * each number in a new 3d array. Those numbers are used
1031  * to create the linear equation system.*/
1032 
1033  if (cell_type == N_CELL_DIRICHLET) {
1034  /* include dirichlet cells in the les */
1035  for (k = 0; k < geom->depths; k++) {
1036  for (j = 0; j < geom->rows; j++) {
1037  for (i = 0; i < geom->cols; i++) {
1038  /*use all non-inactive cells for les creation */
1039  if (N_CELL_INACTIVE <
1040  (int)N_get_array_3d_d_value(status, i, j, k) &&
1041  (int)N_get_array_3d_d_value(status, i, j,
1042  k) < N_MAX_CELL_STATE)
1043  cell_type_count++;
1044  }
1045  }
1046  }
1047  }
1048  else {
1049  /*use only active cell in the les */
1050  for (k = 0; k < geom->depths; k++) {
1051  for (j = 0; j < geom->rows; j++) {
1052  for (i = 0; i < geom->cols; i++) {
1053  /*count only active cells */
1054  if (N_CELL_ACTIVE
1055  == (int)N_get_array_3d_d_value(status, i, j, k))
1056  cell_type_count++;
1057 
1058  }
1059  }
1060  }
1061  }
1062 
1063  G_debug(2,
1064  "N_assemble_les_3d: number of used cells %i\n", cell_type_count);
1065 
1066  if (cell_type_count == 0.0)
1068  ("Not enough active cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
1069  cell_type_count);
1070 
1071  /* allocate the memory for the linear equation system (les).
1072  * Only valid cells are used to create the les. */
1073  les = N_alloc_les_Ax_b(cell_type_count, les_type);
1074 
1075  index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
1076  for (i = 0; i < cell_type_count; i++)
1077  index_ij[i] = (int *)G_calloc(3, sizeof(int));
1078 
1079  count = 0;
1080  /*count the number of cells which should be used to create the linear equation system */
1081  /*save the k, i and j indices and create a ordered numbering */
1082  for (k = 0; k < geom->depths; k++) {
1083  for (j = 0; j < geom->rows; j++) {
1084  for (i = 0; i < geom->cols; i++) {
1085  if (cell_type == N_CELL_DIRICHLET) {
1086  if (N_CELL_INACTIVE <
1087  (int)N_get_array_3d_d_value(status, i, j, k) &&
1088  (int)N_get_array_3d_d_value(status, i, j,
1089  k) < N_MAX_CELL_STATE) {
1090  N_put_array_3d_d_value(cell_count, i, j, k, count);
1091  index_ij[count][0] = i;
1092  index_ij[count][1] = j;
1093  index_ij[count][2] = k;
1094  count++;
1095  G_debug(5,
1096  "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
1097  count, i, j, k);
1098  }
1099  }
1100  else if (N_CELL_ACTIVE ==
1101  (int)N_get_array_3d_d_value(status, i, j, k)) {
1102  N_put_array_3d_d_value(cell_count, i, j, k, count);
1103  index_ij[count][0] = i;
1104  index_ij[count][1] = j;
1105  index_ij[count][2] = k;
1106  count++;
1107  G_debug(5,
1108  "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
1109  count, i, j, k);
1110  }
1111  }
1112  }
1113  }
1114 
1115  G_debug(2, "N_assemble_les_3d: starting the parallel assemble loop");
1116 
1117 #pragma omp parallel for private(i, j, k, pos, count) schedule(static)
1118  for (count = 0; count < cell_type_count; count++) {
1119  i = index_ij[count][0];
1120  j = index_ij[count][1];
1121  k = index_ij[count][2];
1122 
1123  /*create the entries for the */
1124  N_data_star *items = call->callback(data, geom, i, j, k);
1125 
1126  G_math_spvector *spvect = NULL;
1127 
1128  /*allocate a sprase vector */
1129  if (les_type == N_SPARSE_LES)
1130  spvect = G_math_alloc_spvector(items->count);
1131  /* initial conditions */
1132 
1133  les->x[count] = N_get_array_3d_d_value(start_val, i, j, k);
1134 
1135  /* the entry in the vector b */
1136  les->b[count] = items->V;
1137 
1138  /* pos describes the position in the sparse vector.
1139  * the first entry is always the diagonal entry of the matrix*/
1140  pos = 0;
1141 
1142  if (les_type == N_SPARSE_LES) {
1143  spvect->index[pos] = count;
1144  spvect->values[pos] = items->C;
1145  }
1146  else {
1147  les->A[count][count] = items->C;
1148  }
1149  /* western neighbour, entry is col - 1 */
1150  if (i > 0) {
1151  pos =
1152  make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1153  cell_count, status, start_val, items->W,
1154  cell_type);
1155  }
1156  /* eastern neighbour, entry col + 1 */
1157  if (i < geom->cols - 1) {
1158  pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1159  cell_count, status, start_val, items->E,
1160  cell_type);
1161  }
1162  /* northern neighbour, entry row -1 */
1163  if (j > 0) {
1164  pos =
1165  make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1166  cell_count, status, start_val, items->N,
1167  cell_type);
1168  }
1169  /* southern neighbour, entry row +1 */
1170  if (j < geom->rows - 1) {
1171  pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1172  cell_count, status, start_val, items->S,
1173  cell_type);
1174  }
1175  /*only for a 7 star entry needed */
1176  if (items->type == N_7_POINT_STAR || items->type == N_27_POINT_STAR) {
1177  /* the upper cell (top), entry depth + 1 */
1178  if (k < geom->depths - 1) {
1179  pos =
1180  make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1181  spvect, cell_count, status, start_val,
1182  items->T, cell_type);
1183  }
1184  /* the lower cell (bottom), entry depth - 1 */
1185  if (k > 0) {
1186  pos =
1187  make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1188  spvect, cell_count, status, start_val,
1189  items->B, cell_type);
1190  }
1191  }
1192 
1193  /*How many entries in the les */
1194  if (les->type == N_SPARSE_LES) {
1195  spvect->cols = pos + 1;
1196  G_math_add_spvector(les->Asp, spvect, count);
1197  }
1198 
1199  if (items)
1200  G_free(items);
1201  }
1202 
1203  N_free_array_3d(cell_count);
1204 
1205  for (i = 0; i < cell_type_count; i++)
1206  G_free(index_ij[i]);
1207 
1208  G_free(index_ij);
1209 
1210  return les;
1211 }
1212 
1213 /*!
1214  * \brief Integrate Dirichlet or Transmission boundary conditions into the les (3d)
1215  *
1216  * Dirichlet and Transmission boundary conditions will be integrated into
1217  * the provided linear equation system. This is meaningful if
1218  * the les was created with #N_assemble_les_2d_dirichlet, because in
1219  * this case Dirichlet boundary conditions are not automatically included.
1220  *
1221  * The provided les will be modified:
1222  *
1223  * Ax = b will be split into Ax_u + Ax_d = b
1224  *
1225  * x_u - the unknowns
1226  * x_d - the Dirichlet cells
1227  *
1228  * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
1229  *
1230  * | A_u 0 | x_u
1231  * | 0 I | x_d
1232  *
1233  * \param les N_les* -- the linear equation system
1234  * \param geom N_geom_data* -- geometrical data information
1235  * \param status N_array_2d* -- the status array containing the cell types
1236  * \param start_val N_array_2d* -- an array with start values
1237  * \return int -- 1 = success, 0 = failure
1238  * */
1240  N_array_3d * status, N_array_3d * start_val)
1241 {
1242  int rows, cols, depths;
1243  int count = 0;
1244  int i, j, x, y, z, stat;
1245  double *dvect1;
1246  double *dvect2;
1247 
1248  G_debug(2,
1249  "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition");
1250 
1251  rows = geom->rows;
1252  cols = geom->cols;
1253  depths = geom->depths;
1254 
1255  /*we nned to additional vectors */
1256  dvect1 = (double *)G_calloc(les->cols, sizeof(double));
1257  dvect2 = (double *)G_calloc(les->cols, sizeof(double));
1258 
1259  /*fill the first one with the x vector data of Dirichlet cells */
1260  count = 0;
1261  for (z = 0; z < depths; z++) {
1262  for (y = 0; y < rows; y++) {
1263  for (x = 0; x < cols; x++) {
1264  stat = (int)N_get_array_3d_d_value(status, x, y, z);
1265  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1266  dvect1[count] =
1267  N_get_array_3d_d_value(start_val, x, y, z);
1268  count++;
1269  }
1270  else if (stat == N_CELL_ACTIVE) {
1271  dvect1[count] = 0.0;
1272  count++;
1273  }
1274  }
1275  }
1276  }
1277 
1278 #pragma omp parallel default(shared)
1279  {
1280  /*perform the matrix vector product and */
1281  if (les->type == N_SPARSE_LES)
1282  G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
1283  else
1284  G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
1285 #pragma omp for schedule (static) private(i)
1286  for (i = 0; i < les->cols; i++)
1287  les->b[i] = les->b[i] - dvect2[i];
1288  }
1289 
1290  /*now set the Dirichlet cell rows and cols to zero and the
1291  * diagonal entry to 1*/
1292  count = 0;
1293  for (z = 0; z < depths; z++) {
1294  for (y = 0; y < rows; y++) {
1295  for (x = 0; x < cols; x++) {
1296  stat = (int)N_get_array_3d_d_value(status, x, y, z);
1297  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1298  if (les->type == N_SPARSE_LES) {
1299  /*set the rows to zero */
1300  for (i = 0; i < les->Asp[count]->cols; i++)
1301  les->Asp[count]->values[i] = 0.0;
1302  /*set the cols to zero */
1303  for (i = 0; i < les->rows; i++) {
1304  for (j = 0; j < les->Asp[i]->cols; j++) {
1305  if (les->Asp[i]->index[j] == count)
1306  les->Asp[i]->values[j] = 0.0;
1307  }
1308  }
1309 
1310  /*entry on the diagonal */
1311  les->Asp[count]->values[0] = 1.0;
1312 
1313  }
1314  else {
1315  /*set the rows to zero */
1316  for (i = 0; i < les->cols; i++)
1317  les->A[count][i] = 0.0;
1318  /*set the cols to zero */
1319  for (i = 0; i < les->rows; i++)
1320  les->A[i][count] = 0.0;
1321 
1322  /*entry on the diagonal */
1323  les->A[count][count] = 1.0;
1324  }
1325  }
1326  count++;
1327  }
1328  }
1329  }
1330 
1331  return 0;
1332 
1333 }
1334 
1335 /* **************************************************************** */
1336 /* **** make an entry in the les (3d) ***************************** */
1337 /* **************************************************************** */
1338 int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
1339  int offset_k, int count, int pos, N_les * les,
1340  G_math_spvector * spvect, N_array_3d * cell_count,
1341  N_array_3d * status, N_array_3d * start_val,
1342  double entry, int cell_type)
1343 {
1344  int K;
1345  int di = offset_i;
1346  int dj = offset_j;
1347  int dk = offset_k;
1348 
1349  K = (int)N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) -
1350  (int)N_get_array_3d_d_value(cell_count, i, j, k);
1351 
1352  if (cell_type == N_CELL_ACTIVE) {
1353  if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) >
1354  N_CELL_ACTIVE &&
1355  (int)N_get_array_3d_d_value(status, i + di, j + dj,
1356  k + dk) < N_MAX_CELL_STATE)
1357  les->b[count] -=
1358  N_get_array_3d_d_value(start_val, i + di, j + dj,
1359  k + dk) * entry;
1360  else if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk)
1361  == N_CELL_ACTIVE) {
1362  if ((count + K) >= 0 && (count + K) < les->cols) {
1363  G_debug(5,
1364  " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
1365  count, count + K, entry);
1366  pos++;
1367  if (les->type == N_SPARSE_LES) {
1368  spvect->index[pos] = count + K;
1369  spvect->values[pos] = entry;
1370  }
1371  else {
1372  les->A[count][count + K] = entry;
1373  }
1374  }
1375  }
1376  }
1377  else if (cell_type == N_CELL_DIRICHLET) {
1378  if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk)
1379  != N_CELL_INACTIVE) {
1380  if ((count + K) >= 0 && (count + K) < les->cols) {
1381  G_debug(5,
1382  " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
1383  count, count + K, entry);
1384  pos++;
1385  if (les->type == N_SPARSE_LES) {
1386  spvect->index[pos] = count + K;
1387  spvect->values[pos] = entry;
1388  }
1389  else {
1390  les->A[count][count + K] = entry;
1391  }
1392  }
1393  }
1394  }
1395 
1396  return pos;
1397 }
#define CELL_TYPE
Definition: raster.h:11
double S_T
Definition: N_pde.h:278
N_data_star *(* callback)()
Definition: N_pde.h:296
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
double SW_T
Definition: N_pde.h:278
double N
Definition: N_pde.h:276
double S_B
Definition: N_pde.h:280
N_data_star * N_create_9star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double V)
allocate and initialize a 9 point star data structure
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:272
N_les * N_alloc_les_Ax_b(int rows, int type)
Allocate memory for a quadratic linear equation system which includes the Matrix A, vector x and vector b.
Definition: n_les.c:145
#define N_27_POINT_STAR
Definition: N_pde.h:43
N_les * N_assemble_les_2d_dirichlet(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet c...
double E
Definition: N_pde.h:276
G_math_spvector * G_math_alloc_spvector(int)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:76
The row vector of the sparse matrix.
Definition: gmath.h:54
double NW_B
Definition: N_pde.h:280
int N_les_integrate_dirichlet_2d(N_les *les, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (2s)
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: n_arrays.c:990
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: n_arrays.c:726
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: n_arrays.c:72
N_les_callback_3d * N_alloc_les_callback_3d(void)
Allocate the structure holding the callback function.
double C
Definition: N_pde.h:276
N_data_star * N_alloc_5star(void)
allocate a 5 point star data structure
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int count
#define NE
Definition: dataquad.h:30
#define N_SPARSE_LES
Definition: N_pde.h:27
#define N
Definition: e_intersect.c:917
#define N_MAX_CELL_STATE
the maximum number of available cell states (eg: boundary condition, inactiven active) ...
Definition: N_pde.h:38
int N_les_integrate_dirichlet_3d(N_les *les, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (3d)
N_les * N_assemble_les_2d(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
double T
Definition: N_pde.h:278
#define NW
Definition: dataquad.h:29
#define NULL
Definition: ccmath.h:32
double * x
Definition: N_pde.h:74
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row.
Definition: n_arrays.c:524
int depths
Definition: N_pde.h:115
callback structure for 2d matrix assembling
Definition: N_pde.h:294
N_data_star * N_callback_template_2d(void *data, N_geom_data *geom, int col, int row)
A callback template creates a 9 point star structure.
#define x
#define N_CELL_ACTIVE
Definition: N_pde.h:32
int count
Definition: N_pde.h:275
void G_math_Ax_sparse(G_math_spvector **, double *, double *, int)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
#define G_calloc(m, n)
Definition: defs/gis.h:113
double S
Definition: N_pde.h:276
double dy
Definition: N_pde.h:110
#define W
Definition: ogsf.h:140
double W_B
Definition: N_pde.h:280
N_les * N_assemble_les_2d_param(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call, int cell_type)
Assemble a linear equation system (les) based on 2d location data (raster)
#define SW
Definition: dataquad.h:31
unsigned int * index
Definition: gmath.h:58
N_les * N_assemble_les_3d_dirichlet(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells...
Geometric information about the structured grid.
Definition: N_pde.h:103
double dz
Definition: N_pde.h:111
#define DCELL_TYPE
Definition: raster.h:13
N_les * N_assemble_les_3d(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells.
int rows
Definition: N_pde.h:116
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition: n_arrays.c:130
double ** A
Definition: N_pde.h:76
double E_T
Definition: N_pde.h:278
N_data_star * N_alloc_9star(void)
allocate a 9 point star data structure
N_les_callback_2d * N_alloc_les_callback_2d(void)
Allocate the structure holding the callback function.
double V
Definition: N_pde.h:276
N_data_star * N_create_27star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double T, double W_T, double E_T, double N_T, double S_T, double NW_T, double SW_T, double NE_T, double SE_T, double B, double W_B, double E_B, double N_B, double S_B, double NW_B, double SW_B, double NE_B, double SE_B, double V)
allocate and initialize a 27 point star data structure
double W_T
Definition: N_pde.h:278
double N_B
Definition: N_pde.h:280
double SE_B
Definition: N_pde.h:280
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
Definition: n_arrays.c:1175
double * values
Definition: gmath.h:56
N_data_star * N_alloc_7star(void)
allocate a 7 point star data structure
G_math_spvector ** Asp
Definition: N_pde.h:77
void N_set_les_callback_2d_func(N_les_callback_2d *data, N_data_star *(*callback_func_2d)())
Set the callback function which is called while assembling the les in 2d.
#define N_CELL_DIRICHLET
Definition: N_pde.h:33
#define N_CELL_INACTIVE
Definition: N_pde.h:31
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure
N_data_star * N_alloc_27star(void)
allocate a 27 point star data structure
unsigned int cols
Definition: gmath.h:57
int depths
number of depths for 3D data
Definition: gis.h:446
N_les * N_assemble_les_3d_active(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells.
int cols
Definition: N_pde.h:117
#define N_5_POINT_STAR
Definition: N_pde.h:40
int type
Definition: N_pde.h:274
#define SE
Definition: dataquad.h:32
double NE_B
Definition: N_pde.h:280
double SE
Definition: N_pde.h:276
double SW
Definition: N_pde.h:276
double NW_T
Definition: N_pde.h:278
int cols
Number of columns for 2D data.
Definition: gis.h:442
double NW
Definition: N_pde.h:276
int cols
Definition: N_pde.h:79
double N_T
Definition: N_pde.h:278
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row.
Definition: n_arrays.c:314
double SW_B
Definition: N_pde.h:280
void N_set_les_callback_3d_func(N_les_callback_3d *data, N_data_star *(*callback_func_3d)())
Set the callback function which is called while assembling the les in 3d.
double NE
Definition: N_pde.h:276
double E_B
Definition: N_pde.h:280
N_data_star *(* callback)()
Definition: N_pde.h:288
double NE_T
Definition: N_pde.h:278
double * b
Definition: N_pde.h:75
N_les * N_assemble_les_2d_active(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
double B
Definition: N_pde.h:280
N_les * N_assemble_les_3d_param(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call, int cell_type)
Assemble a linear equation system (les) based on 3d location data (g3d)
#define N_7_POINT_STAR
Definition: N_pde.h:41
void G_math_d_Ax(double **, double *, double *, int, int)
Compute the matrix - vector product of matrix A and vector x.
Definition: blas_level_2.c:47
callback structure for 3d matrix assembling
Definition: N_pde.h:286
double SE_T
Definition: N_pde.h:278
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: n_arrays.c:378
double W
Definition: N_pde.h:276
double dx
Definition: N_pde.h:109
int rows
Definition: N_pde.h:78
#define N_9_POINT_STAR
Definition: N_pde.h:42
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: n_arrays.c:780
int rows
Number of rows for 2D data.
Definition: gis.h:438
int G_math_add_spvector(G_math_spvector **, G_math_spvector *, int)
Adds a sparse vector to a sparse matrix at position row.
Definition: sparse_matrix.c:35
int G_debug(int, const char *,...) __attribute__((format(printf
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure
N_data_star * N_callback_template_3d(void *data, N_geom_data *geom, int col, int row, int depth)
A callback template creates a 7 point star structure.
The linear equation system (les) structure.
Definition: N_pde.h:72
int type
Definition: N_pde.h:81