1
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: functions to assemble a linear equation system
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18
19
20 #include <math.h>
21 #include <grass/N_pde.h>
22
23 /* local protos */
24 static int make_les_entry_2d(int i, int j, int offset_i, int offset_j,
25 int count, int pos, N_les * les,
26 G_math_spvector * spvect,
27 N_array_2d * cell_count, N_array_2d * status,
28 N_array_2d * start_val, double entry,
29 int cell_type);
30
31 static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
32 int offset_k, int count, int pos, N_les * les,
33 G_math_spvector * spvect,
34 N_array_3d * cell_count, N_array_3d * status,
35 N_array_3d * start_val, double entry,
36 int cell_type);
37
38 /* *************************************************************** *
39 * ********************** N_alloc_5star ************************** *
40 * *************************************************************** */
41 /*!
42 * \brief allocate a 5 point star data structure
43 *
44 * \return N_data_star *
45 * */
N_alloc_5star(void)46 N_data_star *N_alloc_5star(void)
47 {
48 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
49
50 star->type = N_5_POINT_STAR;
51 star->count = 5;
52 return star;
53 }
54
55 /* *************************************************************** *
56 * ********************* N_alloc_7star *************************** *
57 * *************************************************************** */
58 /*!
59 * \brief allocate a 7 point star data structure
60 *
61 * \return N_data_star *
62 * */
N_alloc_7star(void)63 N_data_star *N_alloc_7star(void)
64 {
65 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
66
67 star->type = N_7_POINT_STAR;
68 star->count = 7;
69 return star;
70 }
71
72 /* *************************************************************** *
73 * ********************* N_alloc_9star *************************** *
74 * *************************************************************** */
75 /*!
76 * \brief allocate a 9 point star data structure
77 *
78 * \return N_data_star *
79 *
80 * \attention The 9 point start is not yet implemented in the matrix assembling function
81 *
82 * */
N_alloc_9star(void)83 N_data_star *N_alloc_9star(void)
84 {
85 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
86
87 star->type = N_9_POINT_STAR;
88 star->count = 9;
89 return star;
90 }
91
92 /* *************************************************************** *
93 * ********************* N_alloc_27star ************************** *
94 * *************************************************************** */
95 /*!
96 * \brief allocate a 27 point star data structure
97 *
98 * \return N_data_star *
99 *
100 * \attention The 27 point start is not yet implemented in the matrix assembling function
101 *
102 * */
N_alloc_27star(void)103 N_data_star *N_alloc_27star(void)
104 {
105 N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
106
107 star->type = N_27_POINT_STAR;
108 star->count = 27;
109 return star;
110 }
111
112 /* *************************************************************** *
113 * ********************** N_create_5star ************************* *
114 * *************************************************************** */
115 /*!
116 * \brief allocate and initialize a 5 point star data structure
117 *
118 * \param C double
119 * \param W double
120 * \param E double
121 * \param N double
122 * \param S double
123 * \param V double
124 * \return N_data_star *
125 * */
N_create_5star(double C,double W,double E,double N,double S,double V)126 N_data_star *N_create_5star(double C, double W, double E, double N,
127 double S, double V)
128 {
129 N_data_star *star = N_alloc_5star();
130
131 star->C = C;
132 star->W = W;
133 star->E = E;
134 star->N = N;
135 star->S = S;
136
137 star->V = V;
138
139 G_debug(5, "N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->W,
140 star->E, star->N, star->S, star->C, star->V);
141
142 return star;
143 }
144
145 /* *************************************************************** *
146 * ************************* N_create_7star ********************** *
147 * *************************************************************** */
148 /*!
149 * \brief allocate and initialize a 7 point star data structure
150 *
151 * \param C double
152 * \param W double
153 * \param E double
154 * \param N double
155 * \param S double
156 * \param T double
157 * \param B double
158 * \param V double
159 * \return N_data_star *
160 * */
N_create_7star(double C,double W,double E,double N,double S,double T,double B,double V)161 N_data_star *N_create_7star(double C, double W, double E, double N,
162 double S, double T, double B, double V)
163 {
164 N_data_star *star = N_alloc_7star();
165
166 star->C = C;
167 star->W = W;
168 star->E = E;
169 star->N = N;
170 star->S = S;
171
172 star->T = T;
173 star->B = B;
174
175 star->V = V;
176
177 G_debug(5, "N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n",
178 star->W, star->E, star->N, star->S, star->T, star->B, star->C,
179 star->V);
180
181 return star;
182 }
183
184 /* *************************************************************** *
185 * ************************ N_create_9star *********************** *
186 * *************************************************************** */
187 /*!
188 * \brief allocate and initialize a 9 point star data structure
189 *
190 * \param C double
191 * \param W double
192 * \param E double
193 * \param N double
194 * \param S double
195 * \param NW double
196 * \param SW double
197 * \param NE double
198 * \param SE double
199 * \param V double
200 * \return N_data_star *
201 * */
N_create_9star(double C,double W,double E,double N,double S,double NW,double SW,double NE,double SE,double V)202 N_data_star *N_create_9star(double C, double W, double E, double N,
203 double S, double NW, double SW, double NE,
204 double SE, double V)
205 {
206 N_data_star *star = N_alloc_9star();
207
208 star->C = C;
209 star->W = W;
210 star->E = E;
211 star->N = N;
212 star->S = S;
213
214 star->NW = NW;
215 star->SW = SW;
216 star->NE = NE;
217 star->SE = SE;
218
219 star->V = V;
220
221 G_debug(5,
222 "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
223 star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
224 star->SE, star->C, star->V);
225
226 return star;
227 }
228
229 /* *************************************************************** *
230 * ************************ N_create_27star *********************** *
231 * *************************************************************** */
232 /*!
233 * \brief allocate and initialize a 27 point star data structure
234 *
235 * \param C double
236 * \param W double
237 * \param E double
238 * \param N double
239 * \param S double
240 * \param NW double
241 * \param SW double
242 * \param NE double
243 * \param SE double
244 * \param T double
245 * \param W_T double
246 * \param E_T double
247 * \param N_T double
248 * \param S_T double
249 * \param NW_T double
250 * \param SW_T double
251 * \param NE_T double
252 * \param SE_T double
253 * \param B double
254 * \param W_B double
255 * \param E_B double
256 * \param N_B double
257 * \param S_B double
258 * \param NW_B double
259 * \param SW_B double
260 * \param NE_B double
261 * \param SE_B double
262 * \param V double
263 * \return N_data_star *
264 * */
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)265 N_data_star *N_create_27star(double C, double W, double E, double N, double S,
266 double NW, double SW, double NE, double SE,
267 double T, double W_T, double E_T, double N_T,
268 double S_T, double NW_T, double SW_T,
269 double NE_T, double SE_T, double B, double W_B,
270 double E_B, double N_B, double S_B, double NW_B,
271 double SW_B, double NE_B, double SE_B, double V)
272 {
273 N_data_star *star = N_alloc_27star();
274
275 star->C = C;
276 star->W = W;
277 star->E = E;
278 star->N = N;
279 star->S = S;
280
281 star->NW = NW;
282 star->SW = SW;
283 star->NE = NE;
284 star->SE = SE;
285
286 star->T = T;
287 star->W_T = W_T;
288 star->E_T = E_T;
289 star->N_T = N_T;
290 star->S_T = S_T;
291
292 star->NW_T = NW_T;
293 star->SW_T = SW_T;
294 star->NE_T = NE_T;
295 star->SE_T = SE_T;
296
297 star->B = B;
298 star->W_B = W_B;
299 star->E_B = E_B;
300 star->N_B = N_B;
301 star->S_B = S_B;
302
303 star->NW_B = NW_B;
304 star->SW_B = SW_B;
305 star->NE_B = NE_B;
306 star->SE_B = SE_B;
307
308 star->V = V;
309
310 G_debug(5,
311 "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
312 star->W, star->E, star->N, star->S, star->NW, star->SW, star->NE,
313 star->SE, star->C, star->V);
314
315 G_debug(5,
316 "N_create_27star: w_t %g e_t %g n_t %g s_t %g nw_t %g sw_t %g ne_t %g se_t %g t %g \n",
317 star->W_T, star->E_T, star->N_T, star->S_T, star->NW_T,
318 star->SW_T, star->NE_T, star->SE_T, star->T);
319
320 G_debug(5,
321 "N_create_27star: w_b %g e_b %g n_b %g s_b %g nw_b %g sw_b %g ne_b %g se_B %g b %g\n",
322 star->W_B, star->E_B, star->N_B, star->S_B, star->NW_B,
323 star->SW_B, star->NE_B, star->SE_B, star->B);
324
325
326
327 return star;
328 }
329
330
331 /* *************************************************************** *
332 * ****************** N_set_les_callback_3d_func ***************** *
333 * *************************************************************** */
334 /*!
335 * \brief Set the callback function which is called while assembling the les in 3d
336 *
337 * \param data N_les_callback_3d *
338 * \param callback_func_3d N_data_star *
339 * \return void
340 * */
341 void
N_set_les_callback_3d_func(N_les_callback_3d * data,N_data_star * (* callback_func_3d)())342 N_set_les_callback_3d_func(N_les_callback_3d * data,
343 N_data_star * (*callback_func_3d) ())
344 {
345 data->callback = callback_func_3d;
346 }
347
348 /* *************************************************************** *
349 * *************** N_set_les_callback_2d_func ******************** *
350 * *************************************************************** */
351 /*!
352 * \brief Set the callback function which is called while assembling the les in 2d
353 *
354 * \param data N_les_callback_2d *
355 * \param callback_func_2d N_data_star *
356 * \return void
357 * */
358 void
N_set_les_callback_2d_func(N_les_callback_2d * data,N_data_star * (* callback_func_2d)())359 N_set_les_callback_2d_func(N_les_callback_2d * data,
360 N_data_star * (*callback_func_2d) ())
361 {
362 data->callback = callback_func_2d;
363 }
364
365 /* *************************************************************** *
366 * ************** N_alloc_les_callback_3d ************************ *
367 * *************************************************************** */
368 /*!
369 * \brief Allocate the structure holding the callback function
370 *
371 * A template callback is set. Use N_set_les_callback_3d_func
372 * to set up a specific function.
373 *
374 * \return N_les_callback_3d *
375 * */
N_alloc_les_callback_3d(void)376 N_les_callback_3d *N_alloc_les_callback_3d(void)
377 {
378 N_les_callback_3d *call;
379
380 call = (N_les_callback_3d *) G_calloc(1, sizeof(N_les_callback_3d *));
381 call->callback = N_callback_template_3d;
382
383 return call;
384 }
385
386 /* *************************************************************** *
387 * *************** N_alloc_les_callback_2d *********************** *
388 * *************************************************************** */
389 /*!
390 * \brief Allocate the structure holding the callback function
391 *
392 * A template callback is set. Use N_set_les_callback_2d_func
393 * to set up a specific function.
394 *
395 * \return N_les_callback_2d *
396 * */
N_alloc_les_callback_2d(void)397 N_les_callback_2d *N_alloc_les_callback_2d(void)
398 {
399 N_les_callback_2d *call;
400
401 call = (N_les_callback_2d *) G_calloc(1, sizeof(N_les_callback_2d *));
402 call->callback = N_callback_template_2d;
403
404 return call;
405 }
406
407 /* *************************************************************** *
408 * ******************** N_callback_template_3d ******************* *
409 * *************************************************************** */
410 /*!
411 * \brief A callback template creates a 7 point star structure
412 *
413 * This is a template callback for mass balance calculation with 7 point stars
414 * based on 3d data (g3d).
415 *
416 * \param data void *
417 * \param geom N_geom_data *
418 * \param depth int
419 * \param row int
420 * \param col int
421 * \return N_data_star *
422 *
423 * */
N_callback_template_3d(void * data,N_geom_data * geom,int col,int row,int depth)424 N_data_star *N_callback_template_3d(void *data, N_geom_data * geom, int col,
425 int row, int depth)
426 {
427 N_data_star *star = N_alloc_7star();
428
429 star->E = 1 / geom->dx;
430 star->W = 1 / geom->dx;
431 star->N = 1 / geom->dy;
432 star->S = 1 / geom->dy;
433 star->T = 1 / geom->dz;
434 star->B = 1 / geom->dz;
435 star->C = -1 * (2 / geom->dx + 2 / geom->dy + 2 / geom->dz);
436 star->V = -1;
437
438 G_debug(5,
439 "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
440 star->W, star->E, star->N, star->S, star->T, star->B, star->C,
441 star->V);
442
443
444 return star;
445 }
446
447 /* *************************************************************** *
448 * ********************* N_callback_template_2d ****************** *
449 * *************************************************************** */
450 /*!
451 * \brief A callback template creates a 9 point star structure
452 *
453 * This is a template callback for mass balance calculation with 9 point stars
454 * based on 2d data (raster).
455 *
456 * \param data void *
457 * \param geom N_geom_data *
458 * \param row int
459 * \param col int
460 * \return N_data_star *
461 *
462 * */
N_callback_template_2d(void * data,N_geom_data * geom,int col,int row)463 N_data_star *N_callback_template_2d(void *data, N_geom_data * geom, int col,
464 int row)
465 {
466 N_data_star *star = N_alloc_9star();
467
468 star->E = 1 / geom->dx;
469 star->NE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
470 star->SE = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
471 star->W = 1 / geom->dx;
472 star->NW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
473 star->SW = 1 / sqrt(geom->dx * geom->dx + geom->dy * geom->dy);
474 star->N = 1 / geom->dy;
475 star->S = 1 / geom->dy;
476 star->C =
477 -1 * (star->E + star->NE + star->SE + star->W + star->NW + star->SW +
478 star->N + star->S);
479 star->V = 0;
480
481 return star;
482 }
483
484 /* *************************************************************** *
485 * ******************** N_assemble_les_2d ************************ *
486 * *************************************************************** */
487 /*!
488 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
489 *
490 * This function calls #N_assemble_les_2d_param
491 *
492 */
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)493 N_les *N_assemble_les_2d(int les_type, N_geom_data * geom,
494 N_array_2d * status, N_array_2d * start_val,
495 void *data, N_les_callback_2d * call)
496 {
497 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
498 call, N_CELL_ACTIVE);
499 }
500
501 /*!
502 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
503 *
504 * This function calls #N_assemble_les_2d_param
505 *
506 */
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)507 N_les *N_assemble_les_2d_active(int les_type, N_geom_data * geom,
508 N_array_2d * status, N_array_2d * start_val,
509 void *data, N_les_callback_2d * call)
510 {
511 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
512 call, N_CELL_ACTIVE);
513 }
514
515 /*!
516 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet cells
517 *
518 * This function calls #N_assemble_les_2d_param
519 *
520 */
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)521 N_les *N_assemble_les_2d_dirichlet(int les_type, N_geom_data * geom,
522 N_array_2d * status,
523 N_array_2d * start_val, void *data,
524 N_les_callback_2d * call)
525 {
526 return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
527 call, N_CELL_DIRICHLET);
528 }
529
530 /*!
531 * \brief Assemble a linear equation system (les) based on 2d location data (raster)
532 *
533 *
534 * The linear equation system type can be set to N_NORMAL_LES to create a regular
535 * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
536 * a new created linear equation system which can be solved with
537 * linear equation solvers. An 2d array with start values and an 2d status array
538 * must be provided as well as the location geometry and a void pointer to data
539 * passed to the callback which creates the les row entries. This callback
540 * must be defined in the N_les_callback_2d strcuture.
541 *
542 * The creation of the les is parallelized with OpenMP.
543 * If you implement new callbacks, please make sure that the
544 * function calls are thread safe.
545 *
546 *
547 * the les can be created in two ways, with dirichlet and similar cells and without them,
548 * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
549 * must be added.
550 *
551 * \param les_type int
552 * \param geom N_geom_data*
553 * \param status N_array_2d *
554 * \param start_val N_array_2d *
555 * \param data void *
556 * \param cell_type int -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
557 * \param call N_les_callback_2d *
558 * \return N_les *
559 * */
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)560 N_les *N_assemble_les_2d_param(int les_type, N_geom_data * geom,
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) &&
587 N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE)
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)
607 G_fatal_error
608 ("Not enough cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
609 cell_type_count);
610
611 /* Then allocate the memory for the linear equation system (les).
612 * Only valid cells are used to create the les. */
613 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
614 for (i = 0; i < cell_type_count; i++)
615 index_ij[i] = (int *)G_calloc(2, sizeof(int));
616
617 les = N_alloc_les_Ax_b(cell_type_count, les_type);
618
619 count = 0;
620
621 /*count the number of cells which should be used to create the linear equation system */
622 /*save the i and j indices and create a ordered numbering */
623 for (j = 0; j < geom->rows; j++) {
624 for (i = 0; i < geom->cols; i++) {
625 /*count every non-inactive cell */
626 if (cell_type == N_CELL_DIRICHLET) {
627 if (N_CELL_INACTIVE < N_get_array_2d_c_value(status, i, j) &&
628 N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE) {
629 N_put_array_2d_c_value(cell_count, i, j, count);
630 index_ij[count][0] = i;
631 index_ij[count][1] = j;
632 count++;
633 G_debug(5,
634 "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
635 count, i, j);
636 }
637 /*count every active cell */
638 }
639 else if (N_CELL_ACTIVE == N_get_array_2d_c_value(status, i, j)) {
640 N_put_array_2d_c_value(cell_count, i, j, count);
641 index_ij[count][0] = i;
642 index_ij[count][1] = j;
643 count++;
644 G_debug(5,
645 "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
646 count, i, j);
647 }
648 }
649 }
650
651 G_debug(2, "N_assemble_les_2d: starting the parallel assemble loop");
652
653 /* Assemble the matrix in parallel */
654 #pragma omp parallel for private(i, j, pos, count) schedule(static)
655 for (count = 0; count < cell_type_count; count++) {
656 i = index_ij[count][0];
657 j = index_ij[count][1];
658
659 /*create the entries for the */
660 N_data_star *items = call->callback(data, geom, i, j);
661
662 /* we need a sparse vector pointer anytime */
663 G_math_spvector *spvect = NULL;
664
665 /*allocate a sprase vector */
666 if (les_type == N_SPARSE_LES) {
667 spvect = G_math_alloc_spvector(items->count);
668 }
669 /* initial conditions */
670 les->x[count] = N_get_array_2d_d_value(start_val, i, j);
671
672 /* the entry in the vector b */
673 les->b[count] = items->V;
674
675 /* pos describes the position in the sparse vector.
676 * the first entry is always the diagonal entry of the matrix*/
677 pos = 0;
678
679 if (les_type == N_SPARSE_LES) {
680 spvect->index[pos] = count;
681 spvect->values[pos] = items->C;
682 }
683 else {
684 les->A[count][count] = items->C;
685 }
686 /* western neighbour, entry is col - 1 */
687 if (i > 0) {
688 pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
689 cell_count, status, start_val, items->W,
690 cell_type);
691 }
692 /* eastern neighbour, entry col + 1 */
693 if (i < geom->cols - 1) {
694 pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
695 cell_count, status, start_val, items->E,
696 cell_type);
697 }
698 /* northern neighbour, entry row - 1 */
699 if (j > 0) {
700 pos =
701 make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
702 cell_count, status, start_val, items->N,
703 cell_type);
704 }
705 /* southern neighbour, entry row + 1 */
706 if (j < geom->rows - 1) {
707 pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
708 cell_count, status, start_val, items->S,
709 cell_type);
710 }
711 /*in case of a nine point star, we have additional entries */
712 if (items->type == N_9_POINT_STAR) {
713 /* north-western neighbour, entry is col - 1 row - 1 */
714 if (i > 0 && j > 0) {
715 pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
716 cell_count, status, start_val,
717 items->NW, cell_type);
718 }
719 /* north-eastern neighbour, entry col + 1 row - 1 */
720 if (i < geom->cols - 1 && j > 0) {
721 pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
722 cell_count, status, start_val,
723 items->NE, cell_type);
724 }
725 /* south-western neighbour, entry is col - 1 row + 1 */
726 if (i > 0 && j < geom->rows - 1) {
727 pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect,
728 cell_count, status, start_val,
729 items->SW, cell_type);
730 }
731 /* south-eastern neighbour, entry col + 1 row + 1 */
732 if (i < geom->cols - 1 && j < geom->rows - 1) {
733 pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
734 cell_count, status, start_val,
735 items->SE, cell_type);
736 }
737 }
738
739 /*How many entries in the les */
740 if (les->type == N_SPARSE_LES) {
741 spvect->cols = pos + 1;
742 G_math_add_spvector(les->Asp, spvect, count);
743 }
744
745 if (items)
746 G_free(items);
747 }
748
749 /*release memory */
750 N_free_array_2d(cell_count);
751
752 for (i = 0; i < cell_type_count; i++)
753 G_free(index_ij[i]);
754
755 G_free(index_ij);
756
757 return les;
758 }
759
760 /*!
761 * \brief Integrate Dirichlet or Transmission boundary conditions into the les (2s)
762 *
763 * Dirichlet and Transmission boundary conditions will be integrated into
764 * the provided linear equation system. This is meaningful if
765 * the les was created with #N_assemble_les_2d_dirichlet, because in
766 * this case Dirichlet boundary conditions are not automatically included.
767 *
768 * The provided les will be modified:
769 *
770 * Ax = b will be split into Ax_u + Ax_d = b
771 *
772 * x_u - the unknowns
773 * x_d - the Dirichlet cells
774 *
775 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
776 *
777 * | A_u 0 | x_u
778 * | 0 I | x_d
779 *
780 * \param les N_les* -- the linear equation system
781 * \param geom N_geom_data* -- geometrical data information
782 * \param status N_array_2d* -- the status array containing the cell types
783 * \param start_val N_array_2d* -- an array with start values
784 * \return int -- 1 = success, 0 = failure
785 * */
N_les_integrate_dirichlet_2d(N_les * les,N_geom_data * geom,N_array_2d * status,N_array_2d * start_val)786 int N_les_integrate_dirichlet_2d(N_les * les, N_geom_data * geom,
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 /* **************************************************************** */
make_les_entry_2d(int i,int j,int offset_i,int offset_j,int count,int pos,N_les * les,G_math_spvector * spvect,N_array_2d * cell_count,N_array_2d * status,N_array_2d * start_val,double entry,int cell_type)880 int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count,
881 int pos, N_les * les, G_math_spvector * spvect,
882 N_array_2d * cell_count, N_array_2d * status,
883 N_array_2d * start_val, double entry, int cell_type)
884 {
885 int K;
886 int di = offset_i;
887 int dj = offset_j;
888
889 K = N_get_array_2d_c_value(cell_count, i + di, j + dj) -
890 N_get_array_2d_c_value(cell_count, i, j);
891
892 /* active cells build the linear equation system */
893 if (cell_type == N_CELL_ACTIVE) {
894 /* dirichlet or transmission cells must be handled like this */
895 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_ACTIVE &&
896 N_get_array_2d_c_value(status, i + di, j + dj) < N_MAX_CELL_STATE)
897 les->b[count] -=
898 N_get_array_2d_d_value(start_val, i + di, j + dj) * entry;
899 else if (N_get_array_2d_c_value(status, i + di, j + dj) ==
900 N_CELL_ACTIVE) {
901 if ((count + K) >= 0 && (count + K) < les->cols) {
902 G_debug(5,
903 " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
904 count, count + K, entry);
905 pos++;
906 if (les->type == N_SPARSE_LES) {
907 spvect->index[pos] = count + K;
908 spvect->values[pos] = entry;
909 }
910 else {
911 les->A[count][count + K] = entry;
912 }
913 }
914 }
915 } /* if dirichlet cells should be used then check for all valid cell neighbours */
916 else if (cell_type == N_CELL_DIRICHLET) {
917 /* all valid cells */
918 if (N_get_array_2d_c_value(status, i + di, j + dj) > N_CELL_INACTIVE
919 && N_get_array_2d_c_value(status, i + di,
920 j + dj) < N_MAX_CELL_STATE) {
921 if ((count + K) >= 0 && (count + K) < les->cols) {
922 G_debug(5,
923 " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
924 count, count + K, entry);
925 pos++;
926 if (les->type == N_SPARSE_LES) {
927 spvect->index[pos] = count + K;
928 spvect->values[pos] = entry;
929 }
930 else {
931 les->A[count][count + K] = entry;
932 }
933 }
934 }
935 }
936
937 return pos;
938 }
939
940
941 /* *************************************************************** *
942 * ******************** N_assemble_les_3d ************************ *
943 * *************************************************************** */
944 /*!
945 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
946 *
947 * This function calls #N_assemble_les_3d_param
948 * */
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)949 N_les *N_assemble_les_3d(int les_type, N_geom_data * geom,
950 N_array_3d * status, N_array_3d * start_val,
951 void *data, N_les_callback_3d * call)
952 {
953 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
954 call, N_CELL_ACTIVE);
955 }
956
957 /*!
958 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
959 *
960 * This function calls #N_assemble_les_3d_param
961 * */
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)962 N_les *N_assemble_les_3d_active(int les_type, N_geom_data * geom,
963 N_array_3d * status, N_array_3d * start_val,
964 void *data, N_les_callback_3d * call)
965 {
966 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
967 call, N_CELL_ACTIVE);
968 }
969
970 /*!
971 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells
972 *
973 * This function calls #N_assemble_les_3d_param
974 * */
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)975 N_les *N_assemble_les_3d_dirichlet(int les_type, N_geom_data * geom,
976 N_array_3d * status,
977 N_array_3d * start_val, void *data,
978 N_les_callback_3d * call)
979 {
980 return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
981 call, N_CELL_DIRICHLET);
982 }
983
984 /*!
985 * \brief Assemble a linear equation system (les) based on 3d location data (g3d)
986 *
987 * The linear equation system type can be set to N_NORMAL_LES to create a regular
988 * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
989 * a new created linear equation system which can be solved with
990 * linear equation solvers. An 3d array with start values and an 3d status array
991 * must be provided as well as the location geometry and a void pointer to data
992 * passed to the callback which creates the les row entries. This callback
993 * must be defined in the N_les_callback_3d structure.
994 *
995 * The creation of the les is parallelized with OpenMP.
996 * If you implement new callbacks, please make sure that the
997 * function calls are thread safe.
998 *
999 * the les can be created in two ways, with dirichlet and similar cells and without them,
1000 * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
1001 * must be added.
1002 *
1003 * \param les_type int
1004 * \param geom N_geom_data*
1005 * \param status N_array_3d *
1006 * \param start_val N_array_3d *
1007 * \param data void *
1008 * \param call N_les_callback_3d *
1009 * \param cell_type int -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
1010 * \return N_les *
1011 * */
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)1012 N_les *N_assemble_les_3d_param(int les_type, N_geom_data * geom,
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)
1067 G_fatal_error
1068 ("Not enough active cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
1069 cell_type_count);
1070
1071 /* allocate the memory for the linear equation system (les).
1072 * Only valid cells are used to create the les. */
1073 les = N_alloc_les_Ax_b(cell_type_count, les_type);
1074
1075 index_ij = (int **)G_calloc(cell_type_count, sizeof(int *));
1076 for (i = 0; i < cell_type_count; i++)
1077 index_ij[i] = (int *)G_calloc(3, sizeof(int));
1078
1079 count = 0;
1080 /*count the number of cells which should be used to create the linear equation system */
1081 /*save the k, i and j indices and create a ordered numbering */
1082 for (k = 0; k < geom->depths; k++) {
1083 for (j = 0; j < geom->rows; j++) {
1084 for (i = 0; i < geom->cols; i++) {
1085 if (cell_type == N_CELL_DIRICHLET) {
1086 if (N_CELL_INACTIVE <
1087 (int)N_get_array_3d_d_value(status, i, j, k) &&
1088 (int)N_get_array_3d_d_value(status, i, j,
1089 k) < N_MAX_CELL_STATE) {
1090 N_put_array_3d_d_value(cell_count, i, j, k, count);
1091 index_ij[count][0] = i;
1092 index_ij[count][1] = j;
1093 index_ij[count][2] = k;
1094 count++;
1095 G_debug(5,
1096 "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
1097 count, i, j, k);
1098 }
1099 }
1100 else if (N_CELL_ACTIVE ==
1101 (int)N_get_array_3d_d_value(status, i, j, k)) {
1102 N_put_array_3d_d_value(cell_count, i, j, k, count);
1103 index_ij[count][0] = i;
1104 index_ij[count][1] = j;
1105 index_ij[count][2] = k;
1106 count++;
1107 G_debug(5,
1108 "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
1109 count, i, j, k);
1110 }
1111 }
1112 }
1113 }
1114
1115 G_debug(2, "N_assemble_les_3d: starting the parallel assemble loop");
1116
1117 #pragma omp parallel for private(i, j, k, pos, count) schedule(static)
1118 for (count = 0; count < cell_type_count; count++) {
1119 i = index_ij[count][0];
1120 j = index_ij[count][1];
1121 k = index_ij[count][2];
1122
1123 /*create the entries for the */
1124 N_data_star *items = call->callback(data, geom, i, j, k);
1125
1126 G_math_spvector *spvect = NULL;
1127
1128 /*allocate a sprase vector */
1129 if (les_type == N_SPARSE_LES)
1130 spvect = G_math_alloc_spvector(items->count);
1131 /* initial conditions */
1132
1133 les->x[count] = N_get_array_3d_d_value(start_val, i, j, k);
1134
1135 /* the entry in the vector b */
1136 les->b[count] = items->V;
1137
1138 /* pos describes the position in the sparse vector.
1139 * the first entry is always the diagonal entry of the matrix*/
1140 pos = 0;
1141
1142 if (les_type == N_SPARSE_LES) {
1143 spvect->index[pos] = count;
1144 spvect->values[pos] = items->C;
1145 }
1146 else {
1147 les->A[count][count] = items->C;
1148 }
1149 /* western neighbour, entry is col - 1 */
1150 if (i > 0) {
1151 pos =
1152 make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1153 cell_count, status, start_val, items->W,
1154 cell_type);
1155 }
1156 /* eastern neighbour, entry col + 1 */
1157 if (i < geom->cols - 1) {
1158 pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1159 cell_count, status, start_val, items->E,
1160 cell_type);
1161 }
1162 /* northern neighbour, entry row -1 */
1163 if (j > 0) {
1164 pos =
1165 make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1166 cell_count, status, start_val, items->N,
1167 cell_type);
1168 }
1169 /* southern neighbour, entry row +1 */
1170 if (j < geom->rows - 1) {
1171 pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1172 cell_count, status, start_val, items->S,
1173 cell_type);
1174 }
1175 /*only for a 7 star entry needed */
1176 if (items->type == N_7_POINT_STAR || items->type == N_27_POINT_STAR) {
1177 /* the upper cell (top), entry depth + 1 */
1178 if (k < geom->depths - 1) {
1179 pos =
1180 make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1181 spvect, cell_count, status, start_val,
1182 items->T, cell_type);
1183 }
1184 /* the lower cell (bottom), entry depth - 1 */
1185 if (k > 0) {
1186 pos =
1187 make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1188 spvect, cell_count, status, start_val,
1189 items->B, cell_type);
1190 }
1191 }
1192
1193 /*How many entries in the les */
1194 if (les->type == N_SPARSE_LES) {
1195 spvect->cols = pos + 1;
1196 G_math_add_spvector(les->Asp, spvect, count);
1197 }
1198
1199 if (items)
1200 G_free(items);
1201 }
1202
1203 N_free_array_3d(cell_count);
1204
1205 for (i = 0; i < cell_type_count; i++)
1206 G_free(index_ij[i]);
1207
1208 G_free(index_ij);
1209
1210 return les;
1211 }
1212
1213 /*!
1214 * \brief Integrate Dirichlet or Transmission boundary conditions into the les (3d)
1215 *
1216 * Dirichlet and Transmission boundary conditions will be integrated into
1217 * the provided linear equation system. This is meaningful if
1218 * the les was created with #N_assemble_les_2d_dirichlet, because in
1219 * this case Dirichlet boundary conditions are not automatically included.
1220 *
1221 * The provided les will be modified:
1222 *
1223 * Ax = b will be split into Ax_u + Ax_d = b
1224 *
1225 * x_u - the unknowns
1226 * x_d - the Dirichlet cells
1227 *
1228 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
1229 *
1230 * | A_u 0 | x_u
1231 * | 0 I | x_d
1232 *
1233 * \param les N_les* -- the linear equation system
1234 * \param geom N_geom_data* -- geometrical data information
1235 * \param status N_array_2d* -- the status array containing the cell types
1236 * \param start_val N_array_2d* -- an array with start values
1237 * \return int -- 1 = success, 0 = failure
1238 * */
N_les_integrate_dirichlet_3d(N_les * les,N_geom_data * geom,N_array_3d * status,N_array_3d * start_val)1239 int N_les_integrate_dirichlet_3d(N_les * les, N_geom_data * geom,
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 /* **************************************************************** */
make_les_entry_3d(int i,int j,int k,int offset_i,int offset_j,int offset_k,int count,int pos,N_les * les,G_math_spvector * spvect,N_array_3d * cell_count,N_array_3d * status,N_array_3d * start_val,double entry,int cell_type)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 }
1398