← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to lib/gpde/N_les_assemble.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gpde/N_les_assemble.c
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
 
 * *************************************************************** */
39
 
/*!
40
 
 * \brief allocate a 5 point star data structure
41
 
 *
42
 
 * \return N_data_star *
43
 
 * */
44
 
N_data_star *N_alloc_5star(void)
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
 
 * *************************************************************** */
56
 
/*!
57
 
 * \brief allocate a 7 point star data structure
58
 
 *
59
 
 * \return N_data_star *
60
 
 * */
61
 
N_data_star *N_alloc_7star(void)
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
 
 * *************************************************************** */
73
 
/*!
74
 
 * \brief allocate a 9 point star data structure
75
 
 *
76
 
 * \return N_data_star *
77
 
 *
78
 
 * \attention The 9 point start is not yet implemented in the matrix assembling function
79
 
 *
80
 
 * */
81
 
N_data_star *N_alloc_9star(void)
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
 
 * *************************************************************** */
93
 
/*!
94
 
 * \brief allocate a 27 point star data structure
95
 
 *
96
 
 * \return N_data_star *
97
 
 *
98
 
 * \attention The 27 point start is not yet implemented in the matrix assembling function
99
 
 *
100
 
 * */
101
 
N_data_star *N_alloc_27star(void)
102
 
{
103
 
    N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));
104
 
 
105
 
    star->type = N_27_POINT_STAR;
106
 
    star->count = 27;
107
 
    return star;
108
 
}
109
 
 
110
 
/* *************************************************************** * 
111
 
 * ********************** N_create_5star ************************* * 
112
 
 * *************************************************************** */
113
 
/*!
114
 
 * \brief allocate and initialize a 5 point star data structure
115
 
 *
116
 
 * \param C double
117
 
 * \param W double
118
 
 * \param E double
119
 
 * \param N double
120
 
 * \param S double
121
 
 * \param V double
122
 
 * \return N_data_star *
123
 
 * */
124
 
N_data_star *N_create_5star(double C, double W, double E, double N,
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
 
 * *************************************************************** */
146
 
/*!
147
 
 * \brief allocate and initialize a 7 point star data structure
148
 
 *
149
 
 * \param C double
150
 
 * \param W double
151
 
 * \param E double
152
 
 * \param N double
153
 
 * \param S double
154
 
 * \param T double
155
 
 * \param B double
156
 
 * \param V double
157
 
 * \return N_data_star *
158
 
 * */
159
 
N_data_star *N_create_7star(double C, double W, double E, double N,
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
 
 * *************************************************************** */
185
 
/*!
186
 
 * \brief allocate and initialize a 9 point star data structure
187
 
 *
188
 
 * \param C  double
189
 
 * \param W  double
190
 
 * \param E  double
191
 
 * \param N  double
192
 
 * \param S  double
193
 
 * \param NW double
194
 
 * \param SW double
195
 
 * \param NE double
196
 
 * \param SE double
197
 
 * \param V  double
198
 
 * \return N_data_star *
199
 
 * */
200
 
N_data_star *N_create_9star(double C, double W, double E, double N,
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
 
 * *************************************************************** */
230
 
/*!
231
 
 * \brief allocate and initialize a 27 point star data structure
232
 
 *
233
 
 * \param C  double
234
 
 * \param W  double
235
 
 * \param E  double
236
 
 * \param N  double
237
 
 * \param S  double
238
 
 * \param NW double
239
 
 * \param SW double
240
 
 * \param NE double
241
 
 * \param SE double
242
 
 * \param T  double
243
 
 * \param W_T  double
244
 
 * \param E_T  double
245
 
 * \param N_T  double
246
 
 * \param S_T  double
247
 
 * \param NW_T double
248
 
 * \param SW_T double
249
 
 * \param NE_T double
250
 
 * \param SE_T double
251
 
 * \param B  double
252
 
 * \param W_B  double
253
 
 * \param E_B  double
254
 
 * \param N_B  double
255
 
 * \param S_B  double
256
 
 * \param NW_B double
257
 
 * \param SW_B double
258
 
 * \param NE_B double
259
 
 * \param SE_B double
260
 
 * \param V  double
261
 
 * \return N_data_star *
262
 
 * */
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
 
 * *************************************************************** */
332
 
/*!
333
 
 * \brief Set the callback function which is called while assembling the les in 3d
334
 
 *
335
 
 * \param data N_les_callback_3d *
336
 
 * \param callback_func_3d N_data_star *
337
 
 * \return void
338
 
 * */
339
 
void
340
 
N_set_les_callback_3d_func(N_les_callback_3d * data,
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
 
 * *************************************************************** */
349
 
/*!
350
 
 * \brief Set the callback function which is called while assembling the les in 2d
351
 
 *
352
 
 * \param data N_les_callback_2d *
353
 
 * \param callback_func_2d N_data_star * 
354
 
 * \return void
355
 
 * */
356
 
void
357
 
N_set_les_callback_2d_func(N_les_callback_2d * data,
358
 
                           N_data_star * (*callback_func_2d) ())
359
 
{
360
 
    data->callback = callback_func_2d;
361
 
}
362
 
 
363
 
/* *************************************************************** * 
364
 
 * ************** N_alloc_les_callback_3d ************************ * 
365
 
 * *************************************************************** */
366
 
/*!
367
 
 * \brief Allocate the structure holding the callback function
368
 
 *
369
 
 * A template callback is set. Use N_set_les_callback_3d_func
370
 
 * to set up a specific function.
371
 
 *
372
 
 * \return N_les_callback_3d *
373
 
 * */
374
 
N_les_callback_3d *N_alloc_les_callback_3d(void)
375
 
{
376
 
    N_les_callback_3d *call;
377
 
 
378
 
    call = (N_les_callback_3d *) G_calloc(1, sizeof(N_les_callback_3d *));
379
 
    call->callback = N_callback_template_3d;
380
 
 
381
 
    return call;
382
 
}
383
 
 
384
 
/* *************************************************************** * 
385
 
 * *************** N_alloc_les_callback_2d *********************** * 
386
 
 * *************************************************************** */
387
 
/*!
388
 
 * \brief Allocate the structure holding the callback function
389
 
 *
390
 
 * A template callback is set. Use N_set_les_callback_2d_func
391
 
 * to set up a specific function.
392
 
 *
393
 
 * \return N_les_callback_2d *
394
 
 * */
395
 
N_les_callback_2d *N_alloc_les_callback_2d(void)
396
 
{
397
 
    N_les_callback_2d *call;
398
 
 
399
 
    call = (N_les_callback_2d *) G_calloc(1, sizeof(N_les_callback_2d *));
400
 
    call->callback = N_callback_template_2d;
401
 
 
402
 
    return call;
403
 
}
404
 
 
405
 
/* *************************************************************** * 
406
 
 * ******************** N_callback_template_3d ******************* * 
407
 
 * *************************************************************** */
408
 
/*!
409
 
 * \brief A callback template creates a 7 point star structure
410
 
 *
411
 
 * This is a template callback for mass balance calculation with 7 point stars
412
 
 * based on 3d data (g3d).
413
 
 *
414
 
 * \param data void *
415
 
 * \param geom N_geom_data *
416
 
 * \param depth int
417
 
 * \param row   int
418
 
 * \param col   int
419
 
 * \return N_data_star *
420
 
 *
421
 
 * */
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
 
 * *************************************************************** */
448
 
/*!
449
 
 * \brief A callback template creates a 9 point star structure
450
 
 *
451
 
 * This is a template callback for mass balance calculation with 9 point stars
452
 
 * based on 2d data (raster).
453
 
 *
454
 
 * \param data void *
455
 
 * \param geom N_geom_data *
456
 
 * \param row int
457
 
 * \param col int
458
 
 * \return N_data_star *
459
 
 *
460
 
 * */
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
 
 * *************************************************************** */
485
 
/*!
486
 
 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
487
 
 *
488
 
 * This function calls #N_assemble_les_2d_param
489
 
 *
490
 
 */
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
 
 
499
 
/*!
500
 
 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active cells
501
 
 *
502
 
 * This function calls #N_assemble_les_2d_param
503
 
 *
504
 
 */
505
 
N_les *N_assemble_les_2d_active(int les_type, N_geom_data * geom,
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
 
 
513
 
/*!
514
 
 * \brief Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet cells
515
 
 *
516
 
 * This function calls #N_assemble_les_2d_param
517
 
 *
518
 
 */
519
 
N_les *N_assemble_les_2d_dirichlet(int les_type, N_geom_data * geom,
520
 
                                   N_array_2d * status,
521
 
                                   N_array_2d * start_val, void *data,
522
 
                                   N_les_callback_2d * call)
523
 
{
524
 
    return N_assemble_les_2d_param(les_type, geom, status, start_val, data,
525
 
                                   call, N_CELL_DIRICHLET);
526
 
}
527
 
 
528
 
/*!
529
 
 * \brief Assemble a linear equation system (les) based on 2d location data  (raster)
530
 
 *
531
 
 * 
532
 
 * The linear equation system type can be set to N_NORMAL_LES to create a regular
533
 
 * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
534
 
 * a new created linear equation system which can be solved with 
535
 
 * linear equation solvers. An 2d array with start values and an 2d status array
536
 
 * must be provided as well as the location geometry and a void pointer to data 
537
 
 * passed to the callback which creates the les row entries. This callback
538
 
 * must be defined in the N_les_callback_2d strcuture.
539
 
 *
540
 
 * The creation of the les is parallelized with OpenMP. 
541
 
 * If you implement new callbacks, please make sure that the 
542
 
 * function calls are thread safe.
543
 
 *
544
 
 *
545
 
 * the les can be created in two ways, with dirichlet and similar cells and without them,
546
 
 * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
547
 
 * must be added.
548
 
 *
549
 
 * \param les_type int
550
 
 * \param geom      N_geom_data*
551
 
 * \param status    N_array_2d *
552
 
 * \param start_val N_array_2d *
553
 
 * \param data void *
554
 
 * \param cell_type int  -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
555
 
 * \param call N_les_callback_2d *
556
 
 * \return N_les *
557
 
 * */
558
 
N_les *N_assemble_les_2d_param(int les_type, N_geom_data * geom,
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) &&
585
 
                    N_get_array_2d_c_value(status, i, j) < N_MAX_CELL_STATE)
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)
605
 
        G_fatal_error
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
 
 
758
 
/*!
759
 
 * \brief Integrate Dirichlet or Transmission boundary conditions into the les (2s)
760
 
 *
761
 
 * Dirichlet and Transmission boundary conditions will be integrated into
762
 
 * the provided linear equation system. This is meaningfull if
763
 
 * the les was created with #N_assemble_les_2d_dirichlet, because in
764
 
 * this case Dirichlet boundary conditions are not automatically included.
765
 
 *
766
 
 * The provided les will be modified:
767
 
 *
768
 
 * Ax = b will be splitted into Ax_u + Ax_d = b
769
 
 *
770
 
 * x_u - the unknowns
771
 
 * x_d - the Dirichlet cells
772
 
 *
773
 
 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
774
 
 *
775
 
 * | A_u  0 | x_u
776
 
 * |  0   I | x_d
777
 
 *
778
 
 * \param les N_les* -- the linear equation system
779
 
 * \param geom N_geom_data* -- geometrical data information
780
 
 * \param status N_array_2d* -- the status array containing the cell types
781
 
 * \param start_val N_array_2d* -- an array with start values
782
 
 * \return int -- 1 = success, 0 = failure
783
 
 * */
784
 
int N_les_integrate_dirichlet_2d(N_les * les, N_geom_data * geom,
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
 
 * *************************************************************** */
942
 
/*!
943
 
 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
944
 
 *
945
 
 * This function calls #N_assemble_les_3d_param
946
 
 * */
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
 
 
955
 
/*!
956
 
 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active cells
957
 
 *
958
 
 * This function calls #N_assemble_les_3d_param
959
 
 * */
960
 
N_les *N_assemble_les_3d_active(int les_type, N_geom_data * geom,
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
 
 
968
 
/*!
969
 
 * \brief Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells
970
 
 *
971
 
 * This function calls #N_assemble_les_3d_param
972
 
 * */
973
 
N_les *N_assemble_les_3d_dirichlet(int les_type, N_geom_data * geom,
974
 
                                   N_array_3d * status,
975
 
                                   N_array_3d * start_val, void *data,
976
 
                                   N_les_callback_3d * call)
977
 
{
978
 
    return N_assemble_les_3d_param(les_type, geom, status, start_val, data,
979
 
                                   call, N_CELL_DIRICHLET);
980
 
}
981
 
 
982
 
/*!
983
 
 * \brief Assemble a linear equation system (les) based on 3d location data (g3d)
984
 
 *
985
 
 * The linear equation system type can be set to N_NORMAL_LES to create a regular
986
 
 * matrix, or to N_SPARSE_LES to create a sparse matrix. This function returns
987
 
 * a new created linear equation system which can be solved with 
988
 
 * linear equation solvers. An 3d array with start values and an 3d status array
989
 
 * must be provided as well as the location geometry and a void pointer to data 
990
 
 * passed to the callback which creates the les row entries. This callback
991
 
 * must be defined in the N_les_callback_3d structure.
992
 
 * 
993
 
 * The creation of the les is parallelized with OpenMP. 
994
 
 * If you implement new callbacks, please make sure that the 
995
 
 * function calls are thread safe.
996
 
 *
997
 
 * the les can be created in two ways, with dirichlet and similar cells and without them,
998
 
 * to spare some memory. If the les is created with dirichlet cell, the dirichlet boundary condition
999
 
 * must be added.
1000
 
 *
1001
 
 * \param les_type int
1002
 
 * \param geom      N_geom_data*
1003
 
 * \param status    N_array_3d *
1004
 
 * \param start_val N_array_3d *
1005
 
 * \param data void *
1006
 
 * \param call N_les_callback_3d *
1007
 
 * \param cell_type int  -- les assemble based on N_CELL_ACTIVE or N_CELL_DIRICHLET
1008
 
 * \return N_les *
1009
 
 * */
1010
 
N_les *N_assemble_les_3d_param(int les_type, N_geom_data * geom,
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)
1065
 
        G_fatal_error
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
 
 
1211
 
/*!
1212
 
 * \brief Integrate Dirichlet or Transmission boundary conditions into the les (3d)
1213
 
 *
1214
 
 * Dirichlet and Transmission boundary conditions will be integrated into
1215
 
 * the provided linear equation system. This is meaningfull if
1216
 
 * the les was created with #N_assemble_les_2d_dirichlet, because in
1217
 
 * this case Dirichlet boundary conditions are not automatically included.
1218
 
 *
1219
 
 * The provided les will be modified:
1220
 
 *
1221
 
 * Ax = b will be splitted into Ax_u + Ax_d = b
1222
 
 *
1223
 
 * x_u - the unknowns
1224
 
 * x_d - the Dirichlet cells
1225
 
 *
1226
 
 * Ax_u = b -Ax_d will be computed. Then the matrix A will be modified to
1227
 
 *
1228
 
 * | A_u  0 | x_u
1229
 
 * |  0   I | x_d
1230
 
 *
1231
 
 * \param les N_les* -- the linear equation system
1232
 
 * \param geom N_geom_data* -- geometrical data information
1233
 
 * \param status N_array_2d* -- the status array containing the cell types
1234
 
 * \param start_val N_array_2d* -- an array with start values
1235
 
 * \return int -- 1 = success, 0 = failure
1236
 
 * */
1237
 
int N_les_integrate_dirichlet_3d(N_les * les, N_geom_data * geom,
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
 
}