GRASS Programmer's Manual  6.4.3(2013)-r
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Macros Pages
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  N_spvector * spvect, N_array_2d * cell_count,
27  N_array_2d * status, N_array_2d * start_val,
28  double entry, int cell_type);
29 
30 static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
31  int offset_k, int count, int pos, N_les * les,
32  N_spvector * spvect, N_array_3d * cell_count,
33  N_array_3d * status, N_array_3d * start_val,
34  double entry, int cell_type);
35 
36 /* *************************************************************** *
37  * ********************** N_alloc_5star ************************** *
38  * *************************************************************** */
45 {
46  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
47 
48  star->type = N_5_POINT_STAR;
49  star->count = 5;
50  return star;
51 }
52 
53 /* *************************************************************** *
54  * ********************* N_alloc_7star *************************** *
55  * *************************************************************** */
62 {
63  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
64 
65  star->type = N_7_POINT_STAR;
66  star->count = 7;
67  return star;
68 }
69 
70 /* *************************************************************** *
71  * ********************* N_alloc_9star *************************** *
72  * *************************************************************** */
82 {
83  N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
84 
85  star->type = N_9_POINT_STAR;
86  star->count = 9;
87  return star;
88 }
89 
90 /* *************************************************************** *
91  * ********************* N_alloc_27star ************************** *
92  * *************************************************************** */
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  * *************************************************************** */
124 N_data_star *N_create_5star(double C, double W, double E, double N,
125  double S, 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  * *************************************************************** */
159 N_data_star *N_create_7star(double C, double W, double E, double N,
160  double S, 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  * *************************************************************** */
200 N_data_star *N_create_9star(double C, double W, double E, double N,
201  double S, double NW, double SW, double NE,
202  double SE, 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 v %g\n",
221  star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
222  star->SE, star->C, star->V);
223 
224  return star;
225 }
226 
227 /* *************************************************************** *
228  * ************************ N_create_27star *********************** *
229  * *************************************************************** */
263 N_data_star *N_create_27star(double C, double W, double E, double N, double S,
264  double NW, double SW, double NE, double SE,
265  double T, double W_T, double E_T, double N_T,
266  double S_T, double NW_T, double SW_T,
267  double NE_T, double SE_T, double B, double W_B,
268  double E_B, double N_B, double S_B, double NW_B,
269  double SW_B, double NE_B, double SE_B, double V)
270 {
271  N_data_star *star = N_alloc_27star();
272 
273  star->C = C;
274  star->W = W;
275  star->E = E;
276  star->N = N;
277  star->S = S;
278 
279  star->NW = NW;
280  star->SW = SW;
281  star->NE = NE;
282  star->SE = SE;
283 
284  star->T = T;
285  star->W_T = W_T;
286  star->E_T = E_T;
287  star->N_T = N_T;
288  star->S_T = S_T;
289 
290  star->NW_T = NW_T;
291  star->SW_T = SW_T;
292  star->NE_T = NE_T;
293  star->SE_T = SE_T;
294 
295  star->B = B;
296  star->W_B = W_B;
297  star->E_B = E_B;
298  star->N_B = N_B;
299  star->S_B = S_B;
300 
301  star->NW_B = NW_B;
302  star->SW_B = SW_B;
303  star->NE_B = NE_B;
304  star->SE_B = SE_B;
305 
306  star->V = V;
307 
308  G_debug(5,
309  "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
310  star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
311  star->SE, star->C, star->V);
312 
313  G_debug(5,
314  "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",
315  star->W_T, star->E_T, star->N_T, star->S_T, star->NW_T,
316  star->SW_T, star->NE_T, star->SE_T, star->T);
317 
318  G_debug(5,
319  "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",
320  star->W_B, star->E_B, star->N_B, star->S_B, star->NW_B,
321  star->SW_B, star->NE_B, star->SE_B, star->B);
322 
323 
324 
325  return star;
326 }
327 
328 
329 /* *************************************************************** *
330  * ****************** N_set_les_callback_3d_func ***************** *
331  * *************************************************************** */
339 void
341  N_data_star * (*callback_func_3d) ())
342 {
343  data->callback = callback_func_3d;
344 }
345 
346 /* *************************************************************** *
347  * *************** N_set_les_callback_2d_func ******************** *
348  * *************************************************************** */
356 void
358  N_data_star * (*callback_func_2d) ())
359 {
360  data->callback = callback_func_2d;
361 }
362 
363 /* *************************************************************** *
364  * ************** N_alloc_les_callback_3d ************************ *
365  * *************************************************************** */
375 {
377 
378  call = (N_les_callback_3d *) G_calloc(1, sizeof(N_les_callback_3d *));
380 
381  return call;
382 }
383 
384 /* *************************************************************** *
385  * *************** N_alloc_les_callback_2d *********************** *
386  * *************************************************************** */
396 {
398 
399  call = (N_les_callback_2d *) G_calloc(1, sizeof(N_les_callback_2d *));
401 
402  return call;
403 }
404 
405 /* *************************************************************** *
406  * ******************** N_callback_template_3d ******************* *
407  * *************************************************************** */
422 N_data_star *N_callback_template_3d(void *data, N_geom_data * geom, int col,
423  int row, int depth)
424 {
425  N_data_star *star = N_alloc_7star();
426 
427  star->E = 1 / geom->dx;
428  star->W = 1 / geom->dx;
429  star->N = 1 / geom->dy;
430  star->S = 1 / geom->dy;
431  star->T = 1 / geom->dz;
432  star->B = 1 / geom->dz;
433  star->C = -1 * (2 / geom->dx + 2 / geom->dy + 2 / geom->dz);
434  star->V = -1;
435 
436  G_debug(5,
437  "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
438  star->W, star->E, star->N, star->S, star->T, star->B, star->C,
439  star->V);
440 
441 
442  return star;
443 }
444 
445 /* *************************************************************** *
446  * ********************* N_callback_template_2d ****************** *
447  * *************************************************************** */
461 N_data_star *N_callback_template_2d(void *data, N_geom_data * geom, int col,
462  int row)
463 {
464  N_data_star *star = N_alloc_9star();
465 
466  star->E = 1 / geom->dx;
467  star->NE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
468  star->SE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
469  star->W = 1 / geom->dx;
470  star->NW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
471  star->SW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
472  star->N = 1 / geom->dy;
473  star->S = 1 / geom->dy;
474  star->C =
475  -1 * (star->E + star->NE + star->SE + star->W + star->NW + star->SW +
476  star->N + star->S);
477  star->V = 0;
478 
479  return star;
480 }
481 
482 /* *************************************************************** *
483  * ******************** N_assemble_les_2d ************************ *
484  * *************************************************************** */
491 N_les *N_assemble_les_2d(int les_type, N_geom_data * geom,
492  N_array_2d * status, N_array_2d * start_val,
493  void *data, N_les_callback_2d * call)
494 {
495  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
496  call, N_CELL_ACTIVE);
497 }
498 
506  N_array_2d * status, N_array_2d * start_val,
507  void *data, N_les_callback_2d * call)
508 {
509  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
510  call, N_CELL_ACTIVE);
511 }
512 
520  N_array_2d * status,
521  N_array_2d * start_val, void *data,
523 {
524  return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
525  call, N_CELL_DIRICHLET);
526 }
527 
559  N_array_2d * status, N_array_2d * start_val,
560  void *data, N_les_callback_2d * call,
561  int cell_type)
562 {
563  int i, j, count = 0, pos = 0;
564  int cell_type_count = 0;
565  int **index_ij;
566  N_array_2d *cell_count;
567  N_les *les = NULL;
568 
569  G_debug(2,
570  "N_assemble_les_2d: starting to assemble the linear equation system");
571 
572  /* At first count the number of valid cells and save the
573  * each number in a new 2d array. Those numbers are used
574  * to create the linear equation system.
575  * */
576 
577  cell_count = N_alloc_array_2d(geom->cols, geom->rows, 1, CELL_TYPE);
578 
579  /* include dirichlet cells in the les */
580  if (cell_type == N_CELL_DIRICHLET) {
581  for (j = 0; j < geom->rows; j++) {
582  for (i = 0; i < geom->cols; i++) {
583  /*use all non-inactive cells for les creation */
584  if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
586  cell_type_count++;
587  }
588  }
589  }
590  /*use only active cell in the les */
591  if (cell_type == N_CELL_ACTIVE) {
592  for (j = 0; j < geom->rows; j++) {
593  for (i = 0; i < geom->cols; i++) {
594  /*count only active cells */
595  if (N_CELL_ACTIVE == N_get_array_2d_d_value(status, i, j))
596  cell_type_count++;
597  }
598  }
599  }
600 
601  G_debug(2, "N_assemble_les_2d: number of used cells %i\n",
602  cell_type_count);
603 
604  if (cell_type_count == 0)
606  ("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.",
607  cell_type_count);
608 
609  /* Then allocate the memory for the linear equation system (les).
610  * Only valid cells are used to create the les. */
611  index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
612  for (i = 0; i < cell_type_count; i++)
613  index_ij[i] = (int *)G_calloc(2, sizeof(int));
614 
615  les = N_alloc_les_Ax_b(cell_type_count, les_type);
616 
617  count = 0;
618 
619  /*count the number of cells which should be used to create the linear equation system */
620  /*save the i and j indices and create a ordered numbering */
621  for (j = 0; j < geom->rows; j++) {
622  for (i = 0; i < geom->cols; i++) {
623  /*count every non-inactive cell */
624  if (cell_type == N_CELL_DIRICHLET) {
625  if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
626  N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE) {
627  N_put_array_2d_c_value(cell_count, i, j, count);
628  index_ij[count][0] = i;
629  index_ij[count][1] = j;
630  count++;
631  G_debug(5,
632  "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
633  count, i, j);
634  }
635  /*count every active cell */
636  }
637  else if (N_CELL_ACTIVE == N_get_array_2d_c_value(status, i, j)) {
638  N_put_array_2d_c_value(cell_count, i, j, count);
639  index_ij[count][0] = i;
640  index_ij[count][1] = j;
641  count++;
642  G_debug(5,
643  "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
644  count, i, j);
645  }
646  }
647  }
648 
649  G_debug(2, "N_assemble_les_2d: starting the parallel assemble loop");
650 
651  /* Assemble the matrix in parallel */
652 #pragma omp parallel for private(i, j, pos, count) schedule(static)
653  for (count = 0; count < cell_type_count; count++) {
654  i = index_ij[count][0];
655  j = index_ij[count][1];
656 
657  /*create the entries for the */
658  N_data_star *items = call->callback(data, geom, i, j);
659 
660  /* we need a sparse vector pointer anytime */
661  N_spvector *spvect = NULL;
662 
663  /*allocate a sprase vector */
664  if (les_type == N_SPARSE_LES) {
665  spvect = N_alloc_spvector(items->count);
666  }
667  /* initial conditions */
668  les->x[count] = N_get_array_2d_d_value(start_val, i, j);
669 
670  /* the entry in the vector b */
671  les->b[count] = items->V;
672 
673  /* pos describes the position in the sparse vector.
674  * the first entry is always the diagonal entry of the matrix*/
675  pos = 0;
676 
677  if (les_type == N_SPARSE_LES) {
678  spvect->index[pos] = count;
679  spvect->values[pos] = items->C;
680  }
681  else {
682  les->A[count][count] = items->C;
683  }
684  /* western neighbour, entry is col - 1 */
685  if (i > 0) {
686  pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
687  cell_count, status, start_val, items->W,
688  cell_type);
689  }
690  /* eastern neighbour, entry col + 1 */
691  if (i < geom->cols - 1) {
692  pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
693  cell_count, status, start_val, items->E,
694  cell_type);
695  }
696  /* northern neighbour, entry row - 1 */
697  if (j > 0) {
698  pos =
699  make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
700  cell_count, status, start_val, items->N,
701  cell_type);
702  }
703  /* southern neighbour, entry row + 1 */
704  if (j < geom->rows - 1) {
705  pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
706  cell_count, status, start_val, items->S,
707  cell_type);
708  }
709  /*in case of a nine point star, we have additional entries */
710  if (items->type == N_9_POINT_STAR) {
711  /* north-western neighbour, entry is col - 1 row - 1 */
712  if (i > 0 && j > 0) {
713  pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
714  cell_count, status, start_val,
715  items->NW, cell_type);
716  }
717  /* north-eastern neighbour, entry col + 1 row - 1 */
718  if (i < geom->cols - 1 && j > 0) {
719  pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
720  cell_count, status, start_val,
721  items->NE, cell_type);
722  }
723  /* south-western neighbour, entry is col - 1 row + 1 */
724  if (i > 0 && j < geom->rows - 1) {
725  pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect,
726  cell_count, status, start_val,
727  items->SW, cell_type);
728  }
729  /* south-eastern neighbour, entry col + 1 row + 1 */
730  if (i < geom->cols - 1 && j < geom->rows - 1) {
731  pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
732  cell_count, status, start_val,
733  items->SE, cell_type);
734  }
735  }
736 
737  /*How many entries in the les */
738  if (les->type == N_SPARSE_LES) {
739  spvect->cols = pos + 1;
740  N_add_spvector_to_les(les, spvect, count);
741  }
742 
743  if (items)
744  G_free(items);
745  }
746 
747  /*release memory */
748  N_free_array_2d(cell_count);
749 
750  for (i = 0; i < cell_type_count; i++)
751  G_free(index_ij[i]);
752 
753  G_free(index_ij);
754 
755  return les;
756 }
757 
785  N_array_2d * status, N_array_2d * start_val)
786 {
787  int rows, cols;
788  int count = 0;
789  int i, j, x, y, stat;
790  double *dvect1;
791  double *dvect2;
792 
793  G_debug(2,
794  "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition");
795 
796  rows = geom->rows;
797  cols = geom->cols;
798 
799  /*we nned to additional vectors */
800  dvect1 = (double *)G_calloc(les->cols, sizeof(double));
801  dvect2 = (double *)G_calloc(les->cols, sizeof(double));
802 
803  /*fill the first one with the x vector data of Dirichlet cells */
804  count = 0;
805  for (y = 0; y < rows; y++) {
806  for (x = 0; x < cols; x++) {
807  stat = N_get_array_2d_c_value(status, x, y);
808  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
809  dvect1[count] = N_get_array_2d_d_value(start_val, x, y);
810  count++;
811  }
812  else if (stat == N_CELL_ACTIVE) {
813  dvect1[count] = 0.0;
814  count++;
815  }
816  }
817  }
818 
819 #pragma omp parallel default(shared)
820  {
821  /*performe the matrix vector product */
822  if (les->type == N_SPARSE_LES)
823  N_sparse_matrix_vector_product(les, dvect1, dvect2);
824  else
825  N_matrix_vector_product(les, dvect1, dvect2);
826 #pragma omp for schedule (static) private(i)
827  for (i = 0; i < les->cols; i++)
828  les->b[i] = les->b[i] - dvect2[i];
829  }
830 
831  /*now set the Dirichlet cell rows and cols to zero and the
832  * diagonal entry to 1*/
833  count = 0;
834  for (y = 0; y < rows; y++) {
835  for (x = 0; x < cols; x++) {
836  stat = N_get_array_2d_c_value(status, x, y);
837  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
838  if (les->type == N_SPARSE_LES) {
839  /*set the rows to zero */
840  for (i = 0; i < les->Asp[count]->cols; i++)
841  les->Asp[count]->values[i] = 0.0;
842  /*set the cols to zero */
843  for (i = 0; i < les->rows; i++) {
844  for (j = 0; j < les->Asp[i]->cols; j++) {
845  if (les->Asp[i]->index[j] == count)
846  les->Asp[i]->values[j] = 0.0;
847  }
848  }
849 
850  /*entry on the diagonal */
851  les->Asp[count]->values[0] = 1.0;
852 
853  }
854  else {
855  /*set the rows to zero */
856  for (i = 0; i < les->cols; i++)
857  les->A[count][i] = 0.0;
858  /*set the cols to zero */
859  for (i = 0; i < les->rows; i++)
860  les->A[i][count] = 0.0;
861 
862  /*entry on the diagonal */
863  les->A[count][count] = 1.0;
864  }
865  }
866  if (stat >= N_CELL_ACTIVE)
867  count++;
868  }
869  }
870 
871  return 0;
872 
873 }
874 
875 /* **************************************************************** */
876 /* **** make an entry in the les (2d) ***************************** */
877 /* **************************************************************** */
878 int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count,
879  int pos, N_les * les, N_spvector * spvect,
880  N_array_2d * cell_count, N_array_2d * status,
881  N_array_2d * start_val, double entry, int cell_type)
882 {
883  int K;
884  int di = offset_i;
885  int dj = offset_j;
886 
887  K = N_get_array_2d_c_value(cell_count, i + di, j + dj) -
888  N_get_array_2d_c_value(cell_count, i, j);
889 
890  /* active cells build the linear equation system */
891  if (cell_type == N_CELL_ACTIVE) {
892  /* dirichlet or transmission cells must be handled like this */
893  if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_ACTIVE &&
894  N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE)
895  les->b[count] -=
896  N_get_array_2d_d_value(start_val, i + di, j + dj) * entry;
897  else if (N_get_array_2d_c_value(status, i + di, j + dj) ==
898  N_CELL_ACTIVE) {
899  if ((count + K) >= 0 && (count + K) < les->cols) {
900  G_debug(5,
901  " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
902  count, count + K, entry);
903  pos++;
904  if (les->type == N_SPARSE_LES) {
905  spvect->index[pos] = count + K;
906  spvect->values[pos] = entry;
907  }
908  else {
909  les->A[count][count + K] = entry;
910  }
911  }
912  }
913  } /* if dirichlet cells should be used then check for all valid cell neighbours */
914  else if (cell_type == N_CELL_DIRICHLET) {
915  /* all valid cells */
916  if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_INACTIVE
917  && N_get_array_2d_c_value(status, i + di,
918  j + dj) < N_MAX_CELL_STATE) {
919  if ((count + K) >= 0 && (count + K) < les->cols) {
920  G_debug(5,
921  " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
922  count, count + K, entry);
923  pos++;
924  if (les->type == N_SPARSE_LES) {
925  spvect->index[pos] = count + K;
926  spvect->values[pos] = entry;
927  }
928  else {
929  les->A[count][count + K] = entry;
930  }
931  }
932  }
933  }
934 
935  return pos;
936 }
937 
938 
939 /* *************************************************************** *
940  * ******************** N_assemble_les_3d ************************ *
941  * *************************************************************** */
947 N_les *N_assemble_les_3d(int les_type, N_geom_data * geom,
948  N_array_3d * status, N_array_3d * start_val,
949  void *data, N_les_callback_3d * call)
950 {
951  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
952  call, N_CELL_ACTIVE);
953 }
954 
961  N_array_3d * status, N_array_3d * start_val,
962  void *data, N_les_callback_3d * call)
963 {
964  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
965  call, N_CELL_ACTIVE);
966 }
967 
974  N_array_3d * status,
975  N_array_3d * start_val, void *data,
977 {
978  return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
979  call, N_CELL_DIRICHLET);
980 }
981 
1011  N_array_3d * status, N_array_3d * start_val,
1012  void *data, N_les_callback_3d * call,
1013  int cell_type)
1014 {
1015  int i, j, k, count = 0, pos = 0;
1016  int cell_type_count = 0;
1017  N_array_3d *cell_count;
1018  N_les *les = NULL;
1019  int **index_ij;
1020 
1021  G_debug(2,
1022  "N_assemble_les_3d: starting to assemble the linear equation system");
1023 
1024  cell_count =
1025  N_alloc_array_3d(geom->cols, geom->rows, geom->depths, 1, DCELL_TYPE);
1026 
1027  /* First count the number of valid cells and save
1028  * each number in a new 3d array. Those numbers are used
1029  * to create the linear equation system.*/
1030 
1031  if (cell_type == N_CELL_DIRICHLET) {
1032  /* include dirichlet cells in the les */
1033  for (k = 0; k < geom->depths; k++) {
1034  for (j = 0; j < geom->rows; j++) {
1035  for (i = 0; i < geom->cols; i++) {
1036  /*use all non-inactive cells for les creation */
1037  if (N_CELL_INACTIVE <
1038  (int)N_get_array_3d_d_value(status, i, j, k) &&
1039  (int)N_get_array_3d_d_value(status, i, j,
1040  k) < N_MAX_CELL_STATE)
1041  cell_type_count++;
1042  }
1043  }
1044  }
1045  }
1046  else {
1047  /*use only active cell in the les */
1048  for (k = 0; k < geom->depths; k++) {
1049  for (j = 0; j < geom->rows; j++) {
1050  for (i = 0; i < geom->cols; i++) {
1051  /*count only active cells */
1052  if (N_CELL_ACTIVE
1053  == (int)N_get_array_3d_d_value(status, i, j, k))
1054  cell_type_count++;
1055 
1056  }
1057  }
1058  }
1059  }
1060 
1061  G_debug(2,
1062  "N_assemble_les_3d: number of used cells %i\n", cell_type_count);
1063 
1064  if (cell_type_count == 0.0)
1066  ("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.",
1067  cell_type_count);
1068 
1069  /* allocate the memory for the linear equation system (les).
1070  * Only valid cells are used to create the les. */
1071  les = N_alloc_les_Ax_b(cell_type_count, les_type);
1072 
1073  index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
1074  for (i = 0; i < cell_type_count; i++)
1075  index_ij[i] = (int *)G_calloc(3, sizeof(int));
1076 
1077  count = 0;
1078  /*count the number of cells which should be used to create the linear equation system */
1079  /*save the k, i and j indices and create a ordered numbering */
1080  for (k = 0; k < geom->depths; k++) {
1081  for (j = 0; j < geom->rows; j++) {
1082  for (i = 0; i < geom->cols; i++) {
1083  if (cell_type == N_CELL_DIRICHLET) {
1084  if (N_CELL_INACTIVE <
1085  (int)N_get_array_3d_d_value(status, i, j, k) &&
1086  (int)N_get_array_3d_d_value(status, i, j,
1087  k) < N_MAX_CELL_STATE) {
1088  N_put_array_3d_d_value(cell_count, i, j, k, count);
1089  index_ij[count][0] = i;
1090  index_ij[count][1] = j;
1091  index_ij[count][2] = k;
1092  count++;
1093  G_debug(5,
1094  "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
1095  count, i, j, k);
1096  }
1097  }
1098  else if (N_CELL_ACTIVE ==
1099  (int)N_get_array_3d_d_value(status, i, j, k)) {
1100  N_put_array_3d_d_value(cell_count, i, j, k, count);
1101  index_ij[count][0] = i;
1102  index_ij[count][1] = j;
1103  index_ij[count][2] = k;
1104  count++;
1105  G_debug(5,
1106  "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
1107  count, i, j, k);
1108  }
1109  }
1110  }
1111  }
1112 
1113  G_debug(2, "N_assemble_les_3d: starting the parallel assemble loop");
1114 
1115 #pragma omp parallel for private(i, j, k, pos, count) schedule(static)
1116  for (count = 0; count < cell_type_count; count++) {
1117  i = index_ij[count][0];
1118  j = index_ij[count][1];
1119  k = index_ij[count][2];
1120 
1121  /*create the entries for the */
1122  N_data_star *items = call->callback(data, geom, i, j, k);
1123 
1124  N_spvector *spvect = NULL;
1125 
1126  /*allocate a sprase vector */
1127  if (les_type == N_SPARSE_LES)
1128  spvect = N_alloc_spvector(items->count);
1129  /* initial conditions */
1130 
1131  les->x[count] = N_get_array_3d_d_value(start_val, i, j, k);
1132 
1133  /* the entry in the vector b */
1134  les->b[count] = items->V;
1135 
1136  /* pos describes the position in the sparse vector.
1137  * the first entry is always the diagonal entry of the matrix*/
1138  pos = 0;
1139 
1140  if (les_type == N_SPARSE_LES) {
1141  spvect->index[pos] = count;
1142  spvect->values[pos] = items->C;
1143  }
1144  else {
1145  les->A[count][count] = items->C;
1146  }
1147  /* western neighbour, entry is col - 1 */
1148  if (i > 0) {
1149  pos =
1150  make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1151  cell_count, status, start_val, items->W,
1152  cell_type);
1153  }
1154  /* eastern neighbour, entry col + 1 */
1155  if (i < geom->cols - 1) {
1156  pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1157  cell_count, status, start_val, items->E,
1158  cell_type);
1159  }
1160  /* northern neighbour, entry row -1 */
1161  if (j > 0) {
1162  pos =
1163  make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1164  cell_count, status, start_val, items->N,
1165  cell_type);
1166  }
1167  /* southern neighbour, entry row +1 */
1168  if (j < geom->rows - 1) {
1169  pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1170  cell_count, status, start_val, items->S,
1171  cell_type);
1172  }
1173  /*only for a 7 star entry needed */
1174  if (items->type == N_7_POINT_STAR || items->type == N_27_POINT_STAR) {
1175  /* the upper cell (top), entry depth + 1 */
1176  if (k < geom->depths - 1) {
1177  pos =
1178  make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1179  spvect, cell_count, status, start_val,
1180  items->T, cell_type);
1181  }
1182  /* the lower cell (bottom), entry depth - 1 */
1183  if (k > 0) {
1184  pos =
1185  make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1186  spvect, cell_count, status, start_val,
1187  items->B, cell_type);
1188  }
1189  }
1190 
1191  /*How many entries in the les */
1192  if (les->type == N_SPARSE_LES) {
1193  spvect->cols = pos + 1;
1194  N_add_spvector_to_les(les, spvect, count);
1195  }
1196 
1197  if (items)
1198  G_free(items);
1199  }
1200 
1201  N_free_array_3d(cell_count);
1202 
1203  for (i = 0; i < cell_type_count; i++)
1204  G_free(index_ij[i]);
1205 
1206  G_free(index_ij);
1207 
1208  return les;
1209 }
1210 
1238  N_array_3d * status, N_array_3d * start_val)
1239 {
1240  int rows, cols, depths;
1241  int count = 0;
1242  int i, j, x, y, z, stat;
1243  double *dvect1;
1244  double *dvect2;
1245 
1246  G_debug(2,
1247  "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition");
1248 
1249  rows = geom->rows;
1250  cols = geom->cols;
1251  depths = geom->depths;
1252 
1253  /*we nned to additional vectors */
1254  dvect1 = (double *)G_calloc(les->cols, sizeof(double));
1255  dvect2 = (double *)G_calloc(les->cols, sizeof(double));
1256 
1257  /*fill the first one with the x vector data of Dirichlet cells */
1258  count = 0;
1259  for (z = 0; z < depths; z++) {
1260  for (y = 0; y < rows; y++) {
1261  for (x = 0; x < cols; x++) {
1262  stat = (int)N_get_array_3d_d_value(status, x, y, z);
1263  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1264  dvect1[count] =
1265  N_get_array_3d_d_value(start_val, x, y, z);
1266  count++;
1267  }
1268  else if (stat == N_CELL_ACTIVE) {
1269  dvect1[count] = 0.0;
1270  count++;
1271  }
1272  }
1273  }
1274  }
1275 
1276 #pragma omp parallel default(shared)
1277  {
1278  /*performe the matrix vector product and */
1279  if (les->type == N_SPARSE_LES)
1280  N_sparse_matrix_vector_product(les, dvect1, dvect2);
1281  else
1282  N_matrix_vector_product(les, dvect1, dvect2);
1283 #pragma omp for schedule (static) private(i)
1284  for (i = 0; i < les->cols; i++)
1285  les->b[i] = les->b[i] - dvect2[i];
1286  }
1287 
1288  /*now set the Dirichlet cell rows and cols to zero and the
1289  * diagonal entry to 1*/
1290  count = 0;
1291  for (z = 0; z < depths; z++) {
1292  for (y = 0; y < rows; y++) {
1293  for (x = 0; x < cols; x++) {
1294  stat = (int)N_get_array_3d_d_value(status, x, y, z);
1295  if (stat > N_CELL_ACTIVE && stat < N_MAX_CELL_STATE) {
1296  if (les->type == N_SPARSE_LES) {
1297  /*set the rows to zero */
1298  for (i = 0; i < les->Asp[count]->cols; i++)
1299  les->Asp[count]->values[i] = 0.0;
1300  /*set the cols to zero */
1301  for (i = 0; i < les->rows; i++) {
1302  for (j = 0; j < les->Asp[i]->cols; j++) {
1303  if (les->Asp[i]->index[j] == count)
1304  les->Asp[i]->values[j] = 0.0;
1305  }
1306  }
1307 
1308  /*entry on the diagonal */
1309  les->Asp[count]->values[0] = 1.0;
1310 
1311  }
1312  else {
1313  /*set the rows to zero */
1314  for (i = 0; i < les->cols; i++)
1315  les->A[count][i] = 0.0;
1316  /*set the cols to zero */
1317  for (i = 0; i < les->rows; i++)
1318  les->A[i][count] = 0.0;
1319 
1320  /*entry on the diagonal */
1321  les->A[count][count] = 1.0;
1322  }
1323  }
1324  count++;
1325  }
1326  }
1327  }
1328 
1329  return 0;
1330 
1331 }
1332 
1333 /* **************************************************************** */
1334 /* **** make an entry in the les (3d) ***************************** */
1335 /* **************************************************************** */
1336 int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
1337  int offset_k, int count, int pos, N_les * les,
1338  N_spvector * spvect, N_array_3d * cell_count,
1339  N_array_3d * status, N_array_3d * start_val,
1340  double entry, int cell_type)
1341 {
1342  int K;
1343  int di = offset_i;
1344  int dj = offset_j;
1345  int dk = offset_k;
1346 
1347  K = N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) -
1348  N_get_array_3d_d_value(cell_count, i, j, k);
1349 
1350  if (cell_type == N_CELL_ACTIVE) {
1351  if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) >
1352  N_CELL_ACTIVE &&
1353  (int)N_get_array_3d_d_value(status, i + di, j + dj,
1354  k + dk) < N_MAX_CELL_STATE)
1355  les->b[count] -=
1356  N_get_array_3d_d_value(start_val, i + di, j + dj,
1357  k + dk) * entry;
1358  else if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk)
1359  == N_CELL_ACTIVE) {
1360  if ((count + K) >= 0 && (count + K) < les->cols) {
1361  G_debug(5,
1362  " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
1363  count, count + K, entry);
1364  pos++;
1365  if (les->type == N_SPARSE_LES) {
1366  spvect->index[pos] = count + K;
1367  spvect->values[pos] = entry;
1368  }
1369  else {
1370  les->A[count][count + K] = entry;
1371  }
1372  }
1373  }
1374  }
1375  else if (cell_type == N_CELL_DIRICHLET) {
1376  if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk)
1377  != N_CELL_INACTIVE) {
1378  if ((count + K) >= 0 && (count + K) < les->cols) {
1379  G_debug(5,
1380  " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
1381  count, count + K, entry);
1382  pos++;
1383  if (les->type == N_SPARSE_LES) {
1384  spvect->index[pos] = count + K;
1385  spvect->values[pos] = entry;
1386  }
1387  else {
1388  les->A[count][count + K] = entry;
1389  }
1390  }
1391  }
1392  }
1393 
1394  return pos;
1395 }