GRASS GIS 7 Programmer's Manual  7.0.3(2016)-r00000
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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  * *************************************************************** */
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 
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 
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 
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 
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  * *************************************************************** */
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 
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 
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 
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 
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 }
double S_T
Definition: N_pde.h:278
N_data_star *(* callback)()
Definition: N_pde.h:296
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
double NW_B
Definition: N_pde.h:280
void G_math_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:76
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
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
Definition: blas_level_2.c:47
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
int count
#define NE
Definition: dataquad.h:30
#define N_SPARSE_LES
Definition: N_pde.h:27
#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 N_CELL_ACTIVE
Definition: N_pde.h:32
int count
Definition: N_pde.h:275
double S
Definition: N_pde.h:276
double dy
Definition: N_pde.h:110
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
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
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
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
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
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
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
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
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
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
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse matrix at position row.
Definition: sparse_matrix.c:35
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