1 /*============================================================================
2 * Radiation solver operations.
3 *============================================================================*/
4
5 /* This file is part of Code_Saturne, a general-purpose CFD tool.
6
7 Copyright (C) 1998-2021 EDF S.A.
8
9 This program is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the Free Software
11 Foundation; either version 2 of the License, or (at your option) any later
12 version.
13
14 This program is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
17 details.
18
19 You should have received a copy of the GNU General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
21 Street, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 /*----------------------------------------------------------------------------*/
24
25 #include "cs_defs.h"
26 #include "cs_math.h"
27
28 /*----------------------------------------------------------------------------
29 * Standard C library headers
30 *----------------------------------------------------------------------------*/
31
32 #include <assert.h>
33 #include <errno.h>
34 #include <stdio.h>
35 #include <stdarg.h>
36 #include <string.h>
37 #include <math.h>
38 #include <float.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include "bft_error.h"
49 #include "bft_mem.h"
50 #include "bft_printf.h"
51
52 #include "cs_log.h"
53 #include "cs_mesh.h"
54 #include "cs_mesh_quantities.h"
55 #include "cs_parall.h"
56 #include "cs_parameters.h"
57 #include "cs_sles.h"
58 #include "cs_sles_it.h"
59 #include "cs_timer.h"
60
61 #include "cs_rad_transfer.h"
62
63 /*----------------------------------------------------------------------------
64 * Header for the current file
65 *----------------------------------------------------------------------------*/
66
67 #include "cs_rad_transfer_dir.h"
68
69 /*----------------------------------------------------------------------------*/
70
71 BEGIN_C_DECLS
72
73 /*=============================================================================
74 * Additional Doxygen documentation
75 *============================================================================*/
76
77 /*! \file cs_rad_transfer_dir.c */
78
79 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
80
81 /*=============================================================================
82 * Private variable
83 *============================================================================*/
84
85 const int x = 0;
86 const int y = 1;
87 const int z = 2;
88
89 /*============================================================================
90 * Private function definitions
91 *============================================================================*/
92
93 /*----------------------------------------------------------------------------
94 * Initialize quadratures.
95 *----------------------------------------------------------------------------*/
96
97 static void
_init_quadrature(void)98 _init_quadrature(void)
99 {
100 if (cs_glob_rad_transfer_params->vect_s == NULL)
101 BFT_MALLOC(cs_glob_rad_transfer_params->vect_s,
102 cs_glob_rad_transfer_params->ndirs,
103 cs_real_3_t);
104
105 if (cs_glob_rad_transfer_params->angsol == NULL)
106 BFT_MALLOC(cs_glob_rad_transfer_params->angsol,
107 cs_glob_rad_transfer_params->ndirs,
108 cs_real_t);
109 }
110
111 /*----------------------------------------------------------------------------
112 * Vector normalization
113 *
114 * parameters:
115 * vec <-> Vector to be normalized
116 *----------------------------------------------------------------------------*/
117
118 inline static void
_normve(cs_real_t vec[3])119 _normve(cs_real_t vec[3])
120 {
121 cs_real_t norm;
122
123 norm = cs_math_3_norm(vec);
124
125 for (int i = 0; i < 3; i++)
126 vec[i] = vec[i] / norm;
127 }
128
129 /*----------------------------------------------------------------------------
130 * Compute the solid angle associated to a given direction of the Tn quadrature
131 *
132 * parameters:
133 * posnod <-- Node coordinate
134 *----------------------------------------------------------------------------*/
135
136 inline static cs_real_t
_lhuilier(cs_real_t posnod[3][3])137 _lhuilier(cs_real_t posnod[3][3])
138 {
139 cs_real_t a, b, c, p;
140 cs_real_t sol_angle;
141
142 /* Renormalisation for dot products */
143
144 _normve(posnod[0]);
145 _normve(posnod[1]);
146 _normve(posnod[2]);
147
148 /* segment lengths of the curvilinear triangle (R=1) */
149
150 a = acos(cs_math_3_dot_product(posnod[0], posnod[1]));
151 b = acos(cs_math_3_dot_product(posnod[1], posnod[2]));
152 c = acos(cs_math_3_dot_product(posnod[2], posnod[0]));
153
154 /* perimeter */
155
156 p = 0.5 * (a + b + c);
157
158 /* Solid angle */
159
160 sol_angle = 4.0 * atan(sqrt( tan(p / 2.0)
161 * tan((p - a) / 2.0)
162 * tan((p - b) / 2.0)
163 * tan((p - c) / 2.0)
164 ));
165
166 return sol_angle;
167 }
168
169 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
170
171 /*=============================================================================
172 * Public function definitions
173 *============================================================================*/
174
175 /*----------------------------------------------------------------------------*/
176 /*!
177 * \brief Compute a quadrature Sn or Tn
178 */
179 /*----------------------------------------------------------------------------*/
180
181 void
cs_rad_transfer_dir(void)182 cs_rad_transfer_dir(void)
183 {
184 cs_rad_transfer_params_t *cs_rp = cs_glob_rad_transfer_params;
185
186 /* Initializations */
187
188 switch(cs_rp->i_quadrature){
189
190 case CS_RAD_QUADRATURE_S4: /* quadrature S4 : 24 directions */
191 cs_rp->ndirs = 3;
192 break;
193
194 case CS_RAD_QUADRATURE_S6: /* quadrature S6 : 48 directions */
195 cs_rp->ndirs = 6;
196 break;
197
198 case CS_RAD_QUADRATURE_S8: /* quadrature S8 : 80 directions */
199 cs_rp->ndirs = 10;
200 break;
201
202 case 4: /* Quadrature T2 : 32 directions */
203 cs_rp->ndirs = 4;
204 break;
205
206 case 5: /* Quadrature T4 : 128 directions */
207 cs_rp->ndirs = 16;
208 break;
209
210 case 6: /* Quadrature Tn : 8 n^2 directions */
211 cs_rp->ndirs = cs_rp->ndirec * cs_rp->ndirec;
212 break;
213
214 case 7: /* Quadrature 120 directions (LC11) */
215 cs_rp->ndirs = 15;
216 break;
217
218 case 8: /* Quadrature 48 directions (DCT020-2468) */
219 cs_rp->ndirs = 6;
220 break;
221
222 default:
223 assert(0);
224 break;
225
226 }
227
228 _init_quadrature();
229
230 cs_real_t vec[10] = {0.0};
231 cs_real_t weight[5] = {0.0};
232
233 /* =================================================================
234 * Quadrature Sn : n(n+2) directions
235 * ================================================================= */
236
237 /* quadrature S4 : 24 directions */
238
239 if (cs_rp->i_quadrature == 1) {
240 vec[0] = 0.2958759;
241 vec[1] = 0.9082483;
242 weight[0] = 0.5235987;
243 cs_rp->vect_s[0][0] = vec[0];
244 cs_rp->vect_s[1][0] = vec[0];
245 cs_rp->vect_s[2][0] = vec[1];
246 cs_rp->vect_s[0][1] = vec[0];
247 cs_rp->vect_s[1][1] = vec[1];
248 cs_rp->vect_s[2][1] = vec[0];
249 cs_rp->vect_s[0][2] = vec[1];
250 cs_rp->vect_s[1][2] = vec[0];
251 cs_rp->vect_s[2][2] = vec[0];
252 cs_rp->angsol[0] = weight[0];
253 cs_rp->angsol[1] = weight[0];
254 cs_rp->angsol[2] = weight[0];
255 }
256
257 /* quadrature S6 : 48 directions */
258
259 else if (cs_rp->i_quadrature == 2) {
260 vec[0] = 0.183867;
261 vec[1] = 0.6950514;
262 vec[2] = 0.9656013;
263 weight[0] = 0.1609517;
264 weight[1] = 0.3626469;
265
266 cs_rp->vect_s[0][0] = vec[0];
267 cs_rp->vect_s[1][0] = vec[0];
268 cs_rp->vect_s[2][0] = vec[0];
269 cs_rp->vect_s[3][0] = vec[1];
270 cs_rp->vect_s[4][0] = vec[1];
271 cs_rp->vect_s[5][0] = vec[2];
272 cs_rp->vect_s[0][1] = vec[0];
273 cs_rp->vect_s[1][1] = vec[1];
274 cs_rp->vect_s[2][1] = vec[2];
275 cs_rp->vect_s[3][1] = vec[0];
276 cs_rp->vect_s[4][1] = vec[1];
277 cs_rp->vect_s[5][1] = vec[0];
278 cs_rp->vect_s[0][2] = vec[2];
279 cs_rp->vect_s[1][2] = vec[1];
280 cs_rp->vect_s[2][2] = vec[0];
281 cs_rp->vect_s[3][2] = vec[1];
282 cs_rp->vect_s[4][2] = vec[0];
283 cs_rp->vect_s[5][2] = vec[0];
284 cs_rp->angsol[0] = weight[0];
285 cs_rp->angsol[1] = weight[1];
286 cs_rp->angsol[2] = weight[0];
287 cs_rp->angsol[3] = weight[1];
288 cs_rp->angsol[4] = weight[1];
289 cs_rp->angsol[5] = weight[0];
290 }
291
292 /* quadrature S8 : 80 directions */
293
294 else if (cs_rp->i_quadrature == 3) {
295 vec[0] = 0.1422555;
296 vec[1] = 0.5773503;
297 vec[2] = 0.8040087;
298 vec[3] = 0.9795543;
299 weight[0] = 0.0992284;
300 weight[1] = 0.1712359;
301 weight[2] = 0.4617179;
302
303 cs_rp->vect_s[0][0] = vec[0];
304 cs_rp->vect_s[1][0] = vec[0];
305 cs_rp->vect_s[2][0] = vec[0];
306 cs_rp->vect_s[3][0] = vec[0];
307 cs_rp->vect_s[4][0] = vec[1];
308 cs_rp->vect_s[5][0] = vec[1];
309 cs_rp->vect_s[6][0] = vec[1];
310 cs_rp->vect_s[7][0] = vec[2];
311 cs_rp->vect_s[8][0] = vec[2];
312 cs_rp->vect_s[9][0] = vec[3];
313 cs_rp->vect_s[0][1] = vec[0];
314 cs_rp->vect_s[1][1] = vec[1];
315 cs_rp->vect_s[2][1] = vec[2];
316 cs_rp->vect_s[3][1] = vec[3];
317 cs_rp->vect_s[4][1] = vec[0];
318 cs_rp->vect_s[5][1] = vec[1];
319 cs_rp->vect_s[6][1] = vec[2];
320 cs_rp->vect_s[7][1] = vec[0];
321 cs_rp->vect_s[8][1] = vec[1];
322 cs_rp->vect_s[9][1] = vec[0];
323 cs_rp->vect_s[0][2] = vec[3];
324 cs_rp->vect_s[1][2] = vec[2];
325 cs_rp->vect_s[2][2] = vec[1];
326 cs_rp->vect_s[3][2] = vec[0];
327 cs_rp->vect_s[4][2] = vec[2];
328 cs_rp->vect_s[5][2] = vec[1];
329 cs_rp->vect_s[6][2] = vec[0];
330 cs_rp->vect_s[7][2] = vec[1];
331 cs_rp->vect_s[8][2] = vec[0];
332 cs_rp->vect_s[9][2] = vec[0];
333 cs_rp->angsol[0] = weight[1];
334 cs_rp->angsol[1] = weight[0];
335 cs_rp->angsol[2] = weight[0];
336 cs_rp->angsol[3] = weight[1];
337 cs_rp->angsol[4] = weight[0];
338 cs_rp->angsol[5] = weight[2];
339 cs_rp->angsol[6] = weight[0];
340 cs_rp->angsol[7] = weight[0];
341 cs_rp->angsol[8] = weight[0];
342 cs_rp->angsol[9] = weight[1];
343 }
344
345 /* =================================================================
346 * Quadrature Tn : 8n^2 directions
347 * ================================================================= */
348
349 /* Quadrature T2 : 32 directions */
350
351 else if (cs_rp->i_quadrature == 4) {
352 vec[0] = 0.2357022604;
353 vec[1] = 0.9428090416;
354 vec[2] = 0.5773502692;
355 weight[0] = 0.5512855984;
356 weight[1] = 0.3398369095;
357
358 cs_rp->vect_s[0][0] = vec[0];
359 cs_rp->vect_s[1][0] = vec[1];
360 cs_rp->vect_s[2][0] = vec[2];
361 cs_rp->vect_s[3][0] = vec[0];
362 cs_rp->vect_s[0][1] = vec[0];
363 cs_rp->vect_s[1][1] = vec[0];
364 cs_rp->vect_s[2][1] = vec[2];
365 cs_rp->vect_s[3][1] = vec[1];
366 cs_rp->vect_s[0][2] = vec[1];
367 cs_rp->vect_s[1][2] = vec[0];
368 cs_rp->vect_s[2][2] = vec[2];
369 cs_rp->vect_s[3][2] = vec[0];
370 cs_rp->angsol[0] = weight[1];
371 cs_rp->angsol[1] = weight[1];
372 cs_rp->angsol[2] = weight[0];
373 cs_rp->angsol[3] = weight[1];
374 }
375
376 /* Quadrature T4 : 128 directions */
377
378 else if (cs_rp->i_quadrature == 5) {
379 vec[0] = 0.0990147543;
380 vec[1] = 0.4923659639;
381 vec[2] = 0.2357022604;
382 vec[3] = 0.123091491;
383 vec[4] = 0.8616404369;
384 vec[5] = 0.6804138174;
385 vec[6] = 0.5773502692;
386 vec[7] = 0.272165527;
387 vec[8] = 0.990147543;
388 vec[9] = 0.9428090416;
389 weight[0] = 0.0526559083;
390 weight[1] = 0.0995720042;
391 weight[2] = 0.0880369928;
392 weight[3] = 0.1320249278;
393 weight[4] = 0.155210815;
394
395 cs_rp->vect_s[0][0] = vec[0];
396 cs_rp->vect_s[1][0] = vec[1];
397 cs_rp->vect_s[2][0] = vec[2];
398 cs_rp->vect_s[3][0] = vec[3];
399 cs_rp->vect_s[4][0] = vec[4];
400 cs_rp->vect_s[5][0] = vec[5];
401 cs_rp->vect_s[6][0] = vec[6];
402 cs_rp->vect_s[7][0] = vec[7];
403 cs_rp->vect_s[8][0] = vec[3];
404 cs_rp->vect_s[9][0] = vec[8];
405 cs_rp->vect_s[10][0] = vec[9];
406 cs_rp->vect_s[11][0] = vec[4];
407 cs_rp->vect_s[12][0] = vec[5];
408 cs_rp->vect_s[13][0] = vec[1];
409 cs_rp->vect_s[14][0] = vec[2];
410 cs_rp->vect_s[15][0] = vec[0];
411 cs_rp->vect_s[0][1] = vec[0];
412 cs_rp->vect_s[1][1] = vec[3];
413 cs_rp->vect_s[2][1] = vec[2];
414 cs_rp->vect_s[3][1] = vec[1];
415 cs_rp->vect_s[4][1] = vec[3];
416 cs_rp->vect_s[5][1] = vec[7];
417 cs_rp->vect_s[6][1] = vec[6];
418 cs_rp->vect_s[7][1] = vec[5];
419 cs_rp->vect_s[8][1] = vec[4];
420 cs_rp->vect_s[9][1] = vec[0];
421 cs_rp->vect_s[10][1] = vec[2];
422 cs_rp->vect_s[11][1] = vec[1];
423 cs_rp->vect_s[12][1] = vec[5];
424 cs_rp->vect_s[13][1] = vec[4];
425 cs_rp->vect_s[14][1] = vec[9];
426 cs_rp->vect_s[15][1] = vec[8];
427 cs_rp->vect_s[0][2] = vec[8];
428 cs_rp->vect_s[1][2] = vec[4];
429 cs_rp->vect_s[2][2] = vec[9];
430 cs_rp->vect_s[3][2] = vec[4];
431 cs_rp->vect_s[4][2] = vec[1];
432 cs_rp->vect_s[5][2] = vec[5];
433 cs_rp->vect_s[6][2] = vec[6];
434 cs_rp->vect_s[7][2] = vec[5];
435 cs_rp->vect_s[8][2] = vec[1];
436 cs_rp->vect_s[9][2] = vec[0];
437 cs_rp->vect_s[10][2] = vec[2];
438 cs_rp->vect_s[11][2] = vec[3];
439 cs_rp->vect_s[12][2] = vec[7];
440 cs_rp->vect_s[13][2] = vec[3];
441 cs_rp->vect_s[14][2] = vec[2];
442 cs_rp->vect_s[15][2] = vec[0];
443 cs_rp->angsol[0] = weight[0];
444 cs_rp->angsol[1] = weight[1];
445 cs_rp->angsol[2] = weight[2];
446 cs_rp->angsol[3] = weight[1];
447 cs_rp->angsol[4] = weight[1];
448 cs_rp->angsol[5] = weight[3];
449 cs_rp->angsol[6] = weight[4];
450 cs_rp->angsol[7] = weight[3];
451 cs_rp->angsol[8] = weight[1];
452 cs_rp->angsol[9] = weight[0];
453 cs_rp->angsol[10] = weight[2];
454 cs_rp->angsol[11] = weight[1];
455 cs_rp->angsol[12] = weight[3];
456 cs_rp->angsol[13] = weight[1];
457 cs_rp->angsol[14] = weight[2];
458 cs_rp->angsol[15] = weight[0];
459 }
460
461 /* Quadrature Tn: 8*n^2 directions */
462
463 else if (cs_rp->i_quadrature == 6) {
464
465 int nquad = cs_rp->ndirec;
466
467 /* Compute X position and Z position of the center of all the sub
468 * triangles of the main 2D triangle
469
470 * Here for T2: z
471 * /\
472 * / \
473 * / 1 \ ! level 1 (1 triangle)
474 * /______\
475 * /\ /\
476 * / \ 2 / \ ! level 2 (3 triangles)
477 * / 1 \ / 3 \
478 * /______\/______\
479 * x y
480 */
481
482 /* Max number of points of a level */
483
484 int max_npt = 2 * nquad - 1;
485
486 cs_real_3_t *xyz3d;
487 BFT_MALLOC(xyz3d, nquad * max_npt, cs_real_3_t);
488
489 /* z position */
490
491 for (int ii = 0; ii < nquad; ii++) {
492
493 int lev = nquad - ii;
494 int lev_npt = 2 * lev - 1;
495
496 for (int jj = 0; jj <= lev_npt; jj += 2)
497 xyz3d[max_npt * (lev - 1) + jj][2] = (1.0 + 3.0 * ii) / (3.0 * nquad);
498
499 for (int jj = 1; jj <= lev_npt; jj += 2)
500 xyz3d[max_npt * (lev - 1) + jj][2] = (2.0 + 3.0 * ii) / (3.0 * nquad);
501
502 }
503
504 /* y position
505 * y position y for each point of a given level increases of 1/(3*(nbpoint+1)/2)
506 * in the same column of a triangle. So we fill y3d by column (diagonal going from
507 * the top to the bottom left) */
508
509 int jj = 0;
510 int int1 = 0;
511
512 for (int tri = 1; tri <= max_npt; tri++) {
513
514 for (int lev = jj; lev < nquad; lev++)
515 xyz3d[lev * max_npt + tri - 1][1] = (tri + jj) / (3.0 * (max_npt + 1.0) / 2.0);
516
517 int1++;
518
519 if (int1 == 2) {
520 int1 = 0;
521 jj++;
522
523 }
524
525 }
526
527 /* x position */
528 for (int ii = 0; ii < nquad * max_npt; ii++)
529 xyz3d[ii][0] = 1.0 - xyz3d[ii][2] - xyz3d[ii][1];
530
531 /* write xyz3d in cs_rp->vect_s */
532
533 jj = 0;
534
535 for (int lev = 1; lev <= nquad; lev++) {
536
537 int lev_npt = 2 * lev - 1;
538
539 for (int ii = 0; ii < lev_npt; ii++) {
540 cs_rp->vect_s[jj][0] = xyz3d[(lev - 1) * max_npt + ii][0];
541 cs_rp->vect_s[jj][1] = xyz3d[(lev - 1) * max_npt + ii][1];
542 cs_rp->vect_s[jj][2] = xyz3d[(lev - 1) * max_npt + ii][2];
543 jj++;
544 }
545
546 }
547
548 /* Vectors normalisation */
549
550 for (int ii = 0; ii < cs_rp->ndirs; ii++)
551 _normve(cs_rp->vect_s[ii]);
552
553 BFT_FREE(xyz3d);
554 xyz3d = NULL;
555
556 /* All the directions are now known, weights will be computed */
557
558 /* Compute nodes positions (the number at node are the number of the node of a
559 * given level)
560 * z
561 * 1 ! level 1 (1 point)
562 * /\
563 * / \
564 * 1/____\2 ! level 2 (2 points)
565 * /\ /\
566 * / \ / \
567 * 1/____\/____\3 ! level 3 (3 points)
568 * x 2 y */
569
570 /* Max number of points for a level */
571
572 max_npt = nquad + 1;
573
574 BFT_MALLOC(xyz3d, max_npt * (nquad + 1), cs_real_3_t);
575
576 /* z position */
577
578 int lev = 0;
579 int ipt = 0;
580 xyz3d[max_npt * lev + ipt][2] = 1.0;
581
582 for (int iquad = 0; iquad < nquad; iquad++) {
583 lev++;
584
585 /* There are "lev + 1" points at level "lev". Beware, C language is 0-based */
586 int lev_npt = lev + 1;
587 for (ipt = 0; ipt < lev_npt; ipt++)
588 xyz3d[max_npt * lev + ipt][2] = xyz3d[max_npt * (lev - 1)][2] - 1.0 / nquad;
589 }
590
591 /* y position */
592 /* We fill points in y by column (diagonal going from the top to the bottom left) */
593 lev = 0;
594 ipt = 0;
595 xyz3d[max_npt * lev + ipt][1] = 0.0;
596
597 for (int iquad = 0; iquad < nquad; iquad++) {
598 lev++;
599
600 /* There are "lev + 1" points at level "lev". Beware, C language is 0-based */
601 int lev_npt = lev + 1;
602 for (ipt = 0; ipt < lev_npt; ipt++)
603 xyz3d[max_npt * lev + ipt][1] = ipt * (1.0 - xyz3d[max_npt * lev][2]) / lev;
604
605 }
606
607 /* x position (every points exist on the plane x+y+z=1) */
608
609 for (int ii = 0; ii < nquad * max_npt; ii++)
610 xyz3d[ii][0] = 1.0 - xyz3d[ii][2] - xyz3d[ii][1];
611
612 /* Compute the surface and the weight of each triangle */
613
614 int ntri = 0;
615
616 cs_real_33_t posnod;
617
618 /* Number of triangle levels */
619
620 for (lev = 0; lev < nquad; lev++) {
621
622 /* First triangle of the level */
623 posnod[0][0] = xyz3d[max_npt * lev][0];
624 posnod[0][1] = xyz3d[max_npt * lev][1];
625 posnod[0][2] = xyz3d[max_npt * lev][2];
626 posnod[1][0] = xyz3d[max_npt * (lev + 1)][0];
627 posnod[1][1] = xyz3d[max_npt * (lev + 1)][1];
628 posnod[1][2] = xyz3d[max_npt * (lev + 1)][2];
629 posnod[2][0] = xyz3d[max_npt * (lev + 1) + 1][0];
630 posnod[2][1] = xyz3d[max_npt * (lev + 1) + 1][1];
631 posnod[2][2] = xyz3d[max_npt * (lev + 1) + 1][2];
632 cs_rp->angsol[ntri] = _lhuilier(posnod);
633 ntri++;
634
635 /* Number of double triangle for this level */
636
637 for (int ii = 0; ii < lev; ii++) {
638 posnod[0][0] = xyz3d[max_npt * lev + ii][0];
639 posnod[0][1] = xyz3d[max_npt * lev + ii][1];
640 posnod[0][2] = xyz3d[max_npt * lev + ii][2];
641 posnod[1][0] = xyz3d[max_npt * lev + ii + 1][0];
642 posnod[1][1] = xyz3d[max_npt * lev + ii + 1][1];
643 posnod[1][2] = xyz3d[max_npt * lev + ii + 1][2];
644 posnod[2][0] = xyz3d[max_npt * (lev + 1) + ii + 1][0];
645 posnod[2][1] = xyz3d[max_npt * (lev + 1) + ii + 1][1];
646 posnod[2][2] = xyz3d[max_npt * (lev + 1) + ii + 1][2];
647 cs_rp->angsol[ntri] = _lhuilier(posnod);
648 ntri++;
649
650 posnod[0][0] = xyz3d[max_npt * lev + ii + 1][0];
651 posnod[0][1] = xyz3d[max_npt * lev + ii + 1][1];
652 posnod[0][2] = xyz3d[max_npt * lev + ii + 1][2];
653 posnod[1][0] = xyz3d[max_npt * (lev + 1) + ii + 1][0];
654 posnod[1][1] = xyz3d[max_npt * (lev + 1) + ii + 1][1];
655 posnod[1][2] = xyz3d[max_npt * (lev + 1) + ii + 1][2];
656 posnod[2][0] = xyz3d[max_npt * (lev + 1) + ii + 2][0];
657 posnod[2][1] = xyz3d[max_npt * (lev + 1) + ii + 2][1];
658 posnod[2][2] = xyz3d[max_npt * (lev + 1) + ii + 2][2];
659 cs_rp->angsol[ntri] = _lhuilier(posnod);
660 ntri++;
661 }
662
663 }
664
665 BFT_FREE(xyz3d);
666
667 }
668
669 else if (cs_rp->i_quadrature == 7) {
670 vec[0] = 0.963560905;
671 vec[1] = 0.189143308;
672 vec[2] = 0.772965714;
673 vec[3] = 0.448622338;
674 vec[4] = 0.686596043;
675 vec[5] = 0.239106143;
676 vec[6] = 0.879538138;
677 vec[7] = 0.475828397;
678 vec[8] = 0.0;
679 weight[0] = 16 * atan(1.0) / 96.0;
680 weight[1] = 8 * atan(1.0) / 96.0;
681 cs_rp->vect_s[0][0] = vec[0];
682 cs_rp->vect_s[1][0] = vec[1];
683 cs_rp->vect_s[2][0] = vec[1];
684 cs_rp->vect_s[3][0] = vec[2];
685 cs_rp->vect_s[4][0] = vec[3];
686 cs_rp->vect_s[5][0] = vec[3];
687 cs_rp->vect_s[6][0] = vec[4];
688 cs_rp->vect_s[7][0] = vec[4];
689 cs_rp->vect_s[8][0] = vec[5];
690 cs_rp->vect_s[9][0] = vec[6];
691 cs_rp->vect_s[10][0] = vec[7];
692 cs_rp->vect_s[11][0] = vec[8];
693 cs_rp->vect_s[12][0] = vec[8];
694 cs_rp->vect_s[13][0] = vec[7];
695 cs_rp->vect_s[14][0] = vec[6];
696 cs_rp->vect_s[0][1] = vec[1];
697 cs_rp->vect_s[1][1] = vec[0];
698 cs_rp->vect_s[2][1] = vec[1];
699 cs_rp->vect_s[3][1] = vec[3];
700 cs_rp->vect_s[4][1] = vec[2];
701 cs_rp->vect_s[5][1] = vec[2];
702 cs_rp->vect_s[6][1] = vec[5];
703 cs_rp->vect_s[7][1] = vec[4];
704 cs_rp->vect_s[8][1] = vec[4];
705 cs_rp->vect_s[9][1] = vec[7];
706 cs_rp->vect_s[10][1] = vec[6];
707 cs_rp->vect_s[11][1] = vec[6];
708 cs_rp->vect_s[12][1] = vec[7];
709 cs_rp->vect_s[13][1] = vec[8];
710 cs_rp->vect_s[14][1] = vec[8];
711 cs_rp->vect_s[0][2] = vec[1];
712 cs_rp->vect_s[1][2] = vec[1];
713 cs_rp->vect_s[2][2] = vec[0];
714 cs_rp->vect_s[3][2] = vec[3];
715 cs_rp->vect_s[4][2] = vec[3];
716 cs_rp->vect_s[5][2] = vec[2];
717 cs_rp->vect_s[6][2] = vec[4];
718 cs_rp->vect_s[7][2] = vec[5];
719 cs_rp->vect_s[8][2] = vec[4];
720 cs_rp->vect_s[9][2] = vec[8];
721 cs_rp->vect_s[10][2] = vec[8];
722 cs_rp->vect_s[11][2] = vec[7];
723 cs_rp->vect_s[12][2] = vec[6];
724 cs_rp->vect_s[13][2] = vec[6];
725 cs_rp->vect_s[14][2] = vec[7];
726 for (int ii = 0; ii < 9; ii++)
727 cs_rp->angsol[ii] = weight[0];
728 for (int ii = 9; ii < 15; ii++)
729 cs_rp->angsol[ii] = weight[1];
730 }
731
732 else if (cs_rp->i_quadrature == 8) {
733 vec[0] = 0.939848342;
734 vec[1] = 0.241542019;
735 vec[2] = 0.6817799;
736 vec[3] = 0.265240148;
737 weight[0] = 0.2437531319;
738 weight[1] = 0.2798456437;
739
740 cs_rp->vect_s[0][0] = vec[0];
741 cs_rp->vect_s[1][0] = vec[2];
742 cs_rp->vect_s[2][0] = vec[2];
743 cs_rp->vect_s[3][0] = vec[1];
744 cs_rp->vect_s[4][0] = vec[3];
745 cs_rp->vect_s[5][0] = vec[1];
746 cs_rp->vect_s[0][1] = vec[1];
747 cs_rp->vect_s[1][1] = vec[3];
748 cs_rp->vect_s[2][1] = vec[2];
749 cs_rp->vect_s[3][1] = vec[1];
750 cs_rp->vect_s[4][1] = vec[2];
751 cs_rp->vect_s[5][1] = vec[0];
752 cs_rp->vect_s[0][2] = vec[1];
753 cs_rp->vect_s[1][2] = vec[2];
754 cs_rp->vect_s[2][2] = vec[3];
755 cs_rp->vect_s[3][2] = vec[0];
756 cs_rp->vect_s[4][2] = vec[2];
757 cs_rp->vect_s[5][2] = vec[1];
758 cs_rp->angsol[0] = weight[0];
759 cs_rp->angsol[1] = weight[1];
760 cs_rp->angsol[2] = weight[1];
761 cs_rp->angsol[3] = weight[0];
762 cs_rp->angsol[4] = weight[1];
763 cs_rp->angsol[5] = weight[0];
764 }
765 }
766
767 /*----------------------------------------------------------------------------*/
768
769 END_C_DECLS
770