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