1 /*============================================================================
2 * Management of mesh quantities
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <math.h>
34 #include <float.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <assert.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45 #include "bft_error.h"
46 #include "bft_printf.h"
47
48 #include "cs_base.h"
49 #include "cs_halo_perio.h"
50 #include "cs_log.h"
51 #include "cs_math.h"
52 #include "cs_mesh.h"
53 #include "cs_mesh_connect.h"
54 #include "cs_parall.h"
55 #include "cs_bad_cells_regularisation.h"
56
57 /*----------------------------------------------------------------------------
58 * Header for the current file
59 *----------------------------------------------------------------------------*/
60
61 #include "cs_mesh_quantities.h"
62
63 /*----------------------------------------------------------------------------*/
64
65 BEGIN_C_DECLS
66
67 /*----------------------------------------------------------------------------*/
68 /*! \file cs_mesh_quantities.c
69 *
70 * \brief Management of mesh quantities.
71 *
72 * Please refer to the
73 * <a href="../../theory.pdf#meshquantities"><b>geometric quantities</b></a>
74 * section of the theory guide for more informations.
75 */
76 /*----------------------------------------------------------------------------*/
77
78 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
79
80 /*=============================================================================
81 * Local Macro definitions
82 *============================================================================*/
83
84 /*============================================================================
85 * Local type definitions
86 *============================================================================*/
87
88 /*============================================================================
89 * Static global variables
90 *============================================================================*/
91
92 /* Pointer to cs_mesh_quantities_t structure for the main mesh */
93
94 cs_mesh_quantities_t *cs_glob_mesh_quantities = NULL;
95
96 /* Choice of the algorithm for computing gravity centers of the cells */
97
98 static int _cell_cen_algorithm = 0;
99 static int _ajust_face_cog_compat_v11_v52 = 0;
100
101 /* Flag (mask) to activate bad cells correction
102 * CS_BAD_CELLS_WARPED_CORRECTION
103 * CS_FACE_DISTANCE_CLIP
104 * are set as default options
105 * */
106 unsigned cs_glob_mesh_quantities_flag =
107 CS_BAD_CELLS_WARPED_CORRECTION + CS_FACE_DISTANCE_CLIP;
108
109 /* Number of computation updates */
110
111 static int _n_computations = 0;
112
113 /*============================================================================
114 * Prototypes for functions intended for use only by Fortran wrappers.
115 * (descriptions follow, with function bodies).
116 *============================================================================*/
117
118 void
119 cs_f_mesh_quantities_fluid_vol_reductions(void);
120
121 /*=============================================================================
122 * Private function definitions
123 *============================================================================*/
124
125 /*----------------------------------------------------------------------------
126 * Build the geometrical matrix linear gradient correction
127 *
128 * parameters:
129 * m <-- mesh
130 * fvq <-> mesh quantities
131 *----------------------------------------------------------------------------*/
132
133 static void
_compute_corr_grad_lin(const cs_mesh_t * m,cs_mesh_quantities_t * fvq)134 _compute_corr_grad_lin(const cs_mesh_t *m,
135 cs_mesh_quantities_t *fvq)
136 {
137 /* Local variables */
138
139 const int n_cells = m->n_cells;
140 const int n_cells_with_ghosts = m->n_cells_with_ghosts;
141 const int n_i_faces = m->n_i_faces;
142 const int n_b_faces = m->n_b_faces;
143
144 const cs_lnum_t *b_face_cells = m->b_face_cells;
145 const cs_lnum_2_t *restrict i_face_cells
146 = (const cs_lnum_2_t *restrict)m->i_face_cells;
147
148 const cs_real_t *restrict cell_vol = fvq->cell_vol;
149 const cs_real_3_t *restrict i_face_normal
150 = (const cs_real_3_t *restrict)fvq->i_face_normal;
151 const cs_real_3_t *restrict b_face_normal
152 = (const cs_real_3_t *restrict)fvq->b_face_normal;
153 const cs_real_3_t *restrict b_face_cog
154 = (const cs_real_3_t *restrict)fvq->b_face_cog;
155 const cs_real_3_t *restrict i_face_cog
156 = (const cs_real_3_t *restrict)fvq->i_face_cog;
157
158 if (fvq->corr_grad_lin_det == NULL)
159 BFT_MALLOC(fvq->corr_grad_lin_det, n_cells_with_ghosts, cs_real_t);
160 if (fvq->corr_grad_lin == NULL)
161 BFT_MALLOC(fvq->corr_grad_lin, n_cells_with_ghosts, cs_real_33_t);
162
163 cs_real_t *restrict corr_grad_lin_det = fvq->corr_grad_lin_det;
164 cs_real_33_t *restrict corr_grad_lin = fvq->corr_grad_lin;
165
166 /* Initialization */
167 for (cs_lnum_t cell_id = 0; cell_id < n_cells_with_ghosts; cell_id++) {
168 for (cs_lnum_t i = 0; i < 3; i++) {
169 for (cs_lnum_t j = 0; j < 3; j++)
170 corr_grad_lin[cell_id][i][j] = 0.;
171 }
172 }
173
174 /* Internal faces contribution */
175 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
176 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
177 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
178
179 for (cs_lnum_t i = 0; i < 3; i++) {
180 for (cs_lnum_t j = 0; j < 3; j++) {
181 cs_real_t flux = i_face_cog[face_id][i] * i_face_normal[face_id][j];
182 corr_grad_lin[cell_id1][i][j] += flux;
183 corr_grad_lin[cell_id2][i][j] -= flux;
184 }
185 }
186 }
187
188 /* Boundary faces contribution */
189 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
190 cs_lnum_t cell_id = b_face_cells[face_id];
191 for (cs_lnum_t i = 0; i < 3; i++) {
192 for (cs_lnum_t j = 0; j < 3; j++) {
193 cs_real_t flux = b_face_cog[face_id][i] * b_face_normal[face_id][j];
194 corr_grad_lin[cell_id][i][j] += flux;
195 }
196 }
197 }
198
199 /* Matrix inversion */
200 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
201 double cocg11 = corr_grad_lin[cell_id][0][0] / cell_vol[cell_id];
202 double cocg12 = corr_grad_lin[cell_id][1][0] / cell_vol[cell_id];
203 double cocg13 = corr_grad_lin[cell_id][2][0] / cell_vol[cell_id];
204 double cocg21 = corr_grad_lin[cell_id][0][1] / cell_vol[cell_id];
205 double cocg22 = corr_grad_lin[cell_id][1][1] / cell_vol[cell_id];
206 double cocg23 = corr_grad_lin[cell_id][2][1] / cell_vol[cell_id];
207 double cocg31 = corr_grad_lin[cell_id][0][2] / cell_vol[cell_id];
208 double cocg32 = corr_grad_lin[cell_id][1][2] / cell_vol[cell_id];
209 double cocg33 = corr_grad_lin[cell_id][2][2] / cell_vol[cell_id];
210
211 double a11 = cocg22 * cocg33 - cocg32 * cocg23;
212 double a12 = cocg32 * cocg13 - cocg12 * cocg33;
213 double a13 = cocg12 * cocg23 - cocg22 * cocg13;
214 double a21 = cocg31 * cocg23 - cocg21 * cocg33;
215 double a22 = cocg11 * cocg33 - cocg31 * cocg13;
216 double a23 = cocg21 * cocg13 - cocg11 * cocg23;
217 double a31 = cocg21 * cocg32 - cocg31 * cocg22;
218 double a32 = cocg31 * cocg12 - cocg11 * cocg32;
219 double a33 = cocg11 * cocg22 - cocg21 * cocg12;
220
221 double det_inv = cocg11 * a11 + cocg21 * a12 + cocg31 * a13;
222
223 if (fabs(det_inv) >= 1.e-15) {
224 det_inv = 1. / det_inv;
225
226 corr_grad_lin[cell_id][0][0] = a11 * det_inv;
227 corr_grad_lin[cell_id][0][1] = a12 * det_inv;
228 corr_grad_lin[cell_id][0][2] = a13 * det_inv;
229 corr_grad_lin[cell_id][1][0] = a21 * det_inv;
230 corr_grad_lin[cell_id][1][1] = a22 * det_inv;
231 corr_grad_lin[cell_id][1][2] = a23 * det_inv;
232 corr_grad_lin[cell_id][2][0] = a31 * det_inv;
233 corr_grad_lin[cell_id][2][1] = a32 * det_inv;
234 corr_grad_lin[cell_id][2][2] = a33 * det_inv;
235
236 double a1 = corr_grad_lin[cell_id][0][0];
237 double a2 = corr_grad_lin[cell_id][0][1];
238 double a3 = corr_grad_lin[cell_id][0][2];
239 double a4 = corr_grad_lin[cell_id][1][0];
240 double a5 = corr_grad_lin[cell_id][1][1];
241 double a6 = corr_grad_lin[cell_id][1][2];
242 double a7 = corr_grad_lin[cell_id][2][0];
243 double a8 = corr_grad_lin[cell_id][2][1];
244 double a9 = corr_grad_lin[cell_id][2][2];
245
246 double determinant = a1 * (a5*a9 - a8*a6)
247 - a2 * (a4*a9 - a7*a6)
248 + a3 * (a4*a8 - a7*a5);
249
250 corr_grad_lin_det[cell_id] = determinant;
251 }
252 else {
253 corr_grad_lin[cell_id][0][0] = 0.;
254 corr_grad_lin[cell_id][0][1] = 0.;
255 corr_grad_lin[cell_id][0][2] = 0.;
256 corr_grad_lin[cell_id][1][0] = 0.;
257 corr_grad_lin[cell_id][1][1] = 0.;
258 corr_grad_lin[cell_id][1][2] = 0.;
259 corr_grad_lin[cell_id][2][0] = 0.;
260 corr_grad_lin[cell_id][2][1] = 0.;
261 corr_grad_lin[cell_id][2][2] = 0.;
262
263 corr_grad_lin_det[cell_id] = 1.;
264 }
265 }
266
267 if (m->halo != NULL) {
268 cs_halo_sync_var(m->halo, CS_HALO_STANDARD, corr_grad_lin_det);
269 cs_halo_sync_var_strided(m->halo, CS_HALO_STANDARD,
270 (cs_real_t *)corr_grad_lin, 9);
271 /* TODO handle rotational periodicity */
272 }
273 }
274
275 /*----------------------------------------------------------------------------
276 * Compute quantities associated to faces (border or internal)
277 *
278 * parameters:
279 * n_faces <-- number of faces
280 * vtx_coord <-- vertex coordinates
281 * face_vtx_idx <-- "face -> vertices" connectivity index
282 * face_vtx_lst <-- "face -> vertices" connectivity list
283 * face_normal --> surface normal of the face
284 *
285 *
286 * Pi+1
287 * *---------* B : barycenter of the polygon
288 * / . . \
289 * / . . \ Pi : vertices of the polygon
290 * / . . \
291 * / . . Ti \ Ti : triangle
292 * *.........B.........* Pi
293 * Pn-1 \ . . /
294 * \ . . /
295 * \ . . /
296 * \ . T0 . /
297 * *---------*
298 * P0
299 *----------------------------------------------------------------------------*/
300
301 static void
_compute_face_normal(cs_lnum_t n_faces,const cs_real_t vtx_coord[][3],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx[],cs_real_t face_normal[][3])302 _compute_face_normal(cs_lnum_t n_faces,
303 const cs_real_t vtx_coord[][3],
304 const cs_lnum_t face_vtx_idx[],
305 const cs_lnum_t face_vtx[],
306 cs_real_t face_normal[][3])
307 {
308 /* Checking */
309
310 assert(face_normal != NULL || n_faces == 0);
311
312 /* Loop on faces */
313
314 # pragma omp parallel for if (n_faces > CS_THR_MIN)
315 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) {
316
317 /* Define the polygon (P) according to the vertices (Pi) of the face */
318
319 cs_lnum_t s_id = face_vtx_idx[f_id];
320 cs_lnum_t e_id = face_vtx_idx[f_id + 1];
321
322 cs_lnum_t n_face_vertices = e_id - s_id;
323
324 if (n_face_vertices == 3) {
325 const cs_lnum_t v0 = face_vtx[s_id];
326 const cs_lnum_t v1 = face_vtx[s_id+1];
327 const cs_lnum_t v2 = face_vtx[s_id+2];
328 cs_real_t v01[3], v02[3], vn[3];
329 for (cs_lnum_t i = 0; i < 3; i++)
330 v01[i] = vtx_coord[v1][i] - vtx_coord[v0][i];
331 for (cs_lnum_t i = 0; i < 3; i++)
332 v02[i] = vtx_coord[v2][i] - vtx_coord[v0][i];
333 cs_math_3_cross_product(v01, v02, vn);
334 for (cs_lnum_t i = 0; i < 3; i++)
335 face_normal[f_id][i] = 0.5*vn[i];
336 }
337
338 else {
339
340 /* Compute approximate face center coordinates for the polygon */
341
342 cs_real_t a_center[3] = {0, 0, 0};
343 cs_real_t f_norm[3] = {0, 0, 0};
344
345 for (cs_lnum_t j = s_id; j < e_id; j++) {
346 const cs_lnum_t v0 = face_vtx[j];
347 for (cs_lnum_t i = 0; i < 3; i++)
348 a_center[i] += vtx_coord[v0][i];
349 }
350
351 for (cs_lnum_t i = 0; i < 3; i++)
352 a_center[i] /= n_face_vertices;
353
354 /* loop on edges, with implied subdivision into triangles
355 defined by edge and cell center */
356
357 cs_real_t vc0[3], vc1[3], vn[3];
358
359 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
360
361 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
362 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
363
364 for (cs_lnum_t i = 0; i < 3; i++) {
365 vc0[i] = vtx_coord[v0][i] - a_center[i];
366 vc1[i] = vtx_coord[v1][i] - a_center[i];
367 }
368
369 cs_math_3_cross_product(vc0, vc1, vn);
370
371 for (cs_lnum_t i = 0; i < 3; i++)
372 f_norm[i] += vn[i];
373
374 }
375
376 for (cs_lnum_t i = 0; i < 3; i++)
377 face_normal[f_id][i] = 0.5*f_norm[i];
378
379 } /* end of test on triangle */
380
381 } /* end of loop on faces */
382 }
383
384 /*----------------------------------------------------------------------------
385 * Compute quantities associated to faces (border or internal)
386 *
387 * parameters:
388 * n_faces <-- number of faces
389 * vtx_coord <-- vertex coordinates
390 * face_vtx_idx <-- "face -> vertices" connectivity index
391 * face_vtx <-- "face -> vertices" connectivity
392 * face_cog --> coordinates of the center of gravity of the faces
393 * face_normal --> face surface normals
394 *
395 * Pi+1
396 * *---------* B : barycenter of the polygon
397 * / . . \
398 * / . . \ Pi : vertices of the polygon
399 * / . . \
400 * / . . Ti \ Ti : triangle
401 * *.........B.........* Pi
402 * Pn-1 \ . . /
403 * \ . . /
404 * \ . . /
405 * \ . T0 . /
406 * *---------*
407 * P0
408 *----------------------------------------------------------------------------*/
409
410 static void
_compute_face_quantities(const cs_lnum_t n_faces,const cs_real_3_t vtx_coord[],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx[],cs_real_3_t face_cog[],cs_real_3_t face_normal[])411 _compute_face_quantities(const cs_lnum_t n_faces,
412 const cs_real_3_t vtx_coord[],
413 const cs_lnum_t face_vtx_idx[],
414 const cs_lnum_t face_vtx[],
415 cs_real_3_t face_cog[],
416 cs_real_3_t face_normal[])
417 {
418 /* Checking */
419
420 assert(face_cog != NULL || n_faces == 0);
421 assert(face_normal != NULL || n_faces == 0);
422
423 const cs_real_t one_third = 1./3.;
424 const cs_real_t s_epsilon = 1.e-32; /* TODO define better "zero" threshold */
425
426 /* Loop on faces */
427
428 # pragma omp parallel for if (n_faces > CS_THR_MIN)
429 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) {
430
431 /* Define the polygon (P) according to the vertices (Pi) of the face */
432
433 cs_lnum_t s_id = face_vtx_idx[f_id];
434 cs_lnum_t e_id = face_vtx_idx[f_id + 1];
435
436 cs_lnum_t n_face_vertices = e_id - s_id;
437
438 if (n_face_vertices == 3) {
439 const cs_lnum_t v0 = face_vtx[s_id];
440 const cs_lnum_t v1 = face_vtx[s_id+1];
441 const cs_lnum_t v2 = face_vtx[s_id+2];
442 cs_real_t v01[3], v02[3], vn[3];
443 for (cs_lnum_t i = 0; i < 3; i++)
444 face_cog[f_id][i] = one_third * ( vtx_coord[v0][i]
445 + vtx_coord[v1][i]
446 + vtx_coord[v2][i]);
447 for (cs_lnum_t i = 0; i < 3; i++)
448 v01[i] = vtx_coord[v1][i] - vtx_coord[v0][i];
449 for (cs_lnum_t i = 0; i < 3; i++)
450 v02[i] = vtx_coord[v2][i] - vtx_coord[v0][i];
451 cs_math_3_cross_product(v01, v02, vn);
452 for (cs_lnum_t i = 0; i < 3; i++)
453 face_normal[f_id][i] = 0.5*vn[i];
454 }
455
456 else { /* For non-triangle faces, assume a division into triangles
457 joining edges and an approximate face center */
458
459 /* Compute approximate face center coordinates for the polygon */
460
461 cs_real_t a_center[3] = {0, 0, 0};
462 cs_real_t f_center[3] = {0, 0, 0}, f_norm[3] = {0, 0, 0};
463
464 for (cs_lnum_t j = s_id; j < e_id; j++) {
465 const cs_lnum_t v0 = face_vtx[j];
466 for (cs_lnum_t i = 0; i < 3; i++)
467 a_center[i] += vtx_coord[v0][i];
468 }
469
470 for (cs_lnum_t i = 0; i < 3; i++)
471 a_center[i] /= n_face_vertices;
472
473 cs_real_t sum_w = 0;
474
475 /* In most cases, the following 2 loops can be merged into a single loop,
476 but for very bad quality faces, some sub-triangles could be oriented
477 differently, so we use 2 passes for safety. */
478
479 if (n_face_vertices < 8) { /* version with local caching for most cases */
480
481 cs_real_t vc0[3], vc1[3], vn[8][3], vtc[8][3];
482
483 /* First pass (face normal) */
484
485 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
486
487 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
488 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
489
490 for (cs_lnum_t i = 0; i < 3; i++) {
491 vc0[i] = vtx_coord[v0][i] - a_center[i];
492 vc1[i] = vtx_coord[v1][i] - a_center[i];
493 vtc[tri_id][i] = vtx_coord[v1][i] + vtx_coord[v0][i] + a_center[i];
494 }
495
496 cs_math_3_cross_product(vc0, vc1, vn[tri_id]);
497
498 for (cs_lnum_t i = 0; i < 3; i++)
499 f_norm[i] += vn[tri_id][i];
500
501 }
502
503 for (cs_lnum_t i = 0; i < 3; i++)
504 f_norm[i] = 0.5*f_norm[i];
505
506 /* Second pass (face center) */
507
508 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
509
510 cs_real_t w = cs_math_3_norm(vn[tri_id]);
511
512 if (cs_math_3_dot_product(vn[tri_id], f_norm) < 0.0)
513 w *= -1.0;
514
515 for (cs_lnum_t i = 0; i < 3; i++)
516 f_center[i] += w*vtc[tri_id][i];
517
518 sum_w += w;
519
520 }
521
522 }
523 else { /* generic version */
524
525 cs_real_t vc0[3], vc1[3], vn[3], vtc[3];
526
527 /* First pass (face normal) */
528
529 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
530
531 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
532 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
533
534 for (cs_lnum_t i = 0; i < 3; i++) {
535 vc0[i] = vtx_coord[v0][i] - a_center[i];
536 vc1[i] = vtx_coord[v1][i] - a_center[i];
537 }
538
539 cs_math_3_cross_product(vc0, vc1, vn);
540
541 for (cs_lnum_t i = 0; i < 3; i++)
542 f_norm[i] += vn[i];
543
544 }
545
546 for (cs_lnum_t i = 0; i < 3; i++)
547 f_norm[i] = 0.5*f_norm[i];
548
549 /* Second pass (face center) */
550
551 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
552
553 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
554 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
555
556 for (cs_lnum_t i = 0; i < 3; i++) {
557 vc0[i] = vtx_coord[v0][i] - a_center[i];
558 vc1[i] = vtx_coord[v1][i] - a_center[i];
559 vtc[i] = vtx_coord[v1][i] + vtx_coord[v0][i] + a_center[i];
560 }
561
562 cs_math_3_cross_product(vc0, vc1, vn);
563
564 cs_real_t w = cs_math_3_norm(vn);
565
566 if (cs_math_3_dot_product(vn, f_norm) < 0.0)
567 w *= -1.0;
568
569 for (cs_lnum_t i = 0; i < 3; i++)
570 f_center[i] += w*vtc[i];
571
572 sum_w += w;
573
574 }
575
576 }
577
578 for (cs_lnum_t i = 0; i < 3; i++)
579 face_normal[f_id][i] = f_norm[i];
580
581 if (sum_w > s_epsilon) {
582 for (cs_lnum_t i = 0; i < 3; i++)
583 face_cog[f_id][i] = one_third * f_center[i]/sum_w;
584 }
585 else {
586 for (cs_lnum_t i = 0; i < 3; i++)
587 face_cog[f_id][i] = a_center[i];
588 }
589
590 } /* end of test on triangle */
591
592 } /* end of loop on faces */
593 }
594
595 /*----------------------------------------------------------------------------
596 * Compute face surfaces based on face norms.
597 *
598 * parameters:
599 * n_faces <-- number of faces
600 * face_norm <-- face surface normals
601 * face_surf --> face surfaces
602 *----------------------------------------------------------------------------*/
603
604 static void
_compute_face_surface(cs_lnum_t n_faces,const cs_real_t face_norm[],cs_real_t face_surf[])605 _compute_face_surface(cs_lnum_t n_faces,
606 const cs_real_t face_norm[],
607 cs_real_t face_surf[])
608 {
609 # pragma omp parallel for if (n_faces > CS_THR_MIN)
610 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++)
611 face_surf[f_id] = cs_math_3_norm(face_norm + f_id*3);
612 }
613
614 /*----------------------------------------------------------------------------
615 * Adjust the position of face centers to obtain a contribution to cell
616 * volumes consistent with that obtained using a splitting of the face
617 * into triangles using face edges and an approxmate center.
618 *
619 * This adjustment was done by default from Code_Saturne 1.1 to 5.2.
620 * To handle warped faces, a correction of the cell center is used so
621 * that the contribution to the cell volume using Green's theorem is
622 * consistent with the volume computed for sub-faces.
623 *
624 * It is based on using the coordinates origin instead of the cell centers
625 * to compute the correction, whereas the volume computation uses the cell
626 * center for better rounding behavior, and uses the cell vertices center
627 * instead of the computed face center, so the consistency gain with
628 * warped faces (where this is supposed to be useful) would need to be
629 * checked.
630 *
631 * parameters:
632 * n_faces <-- number of faces
633 * vtx_coord <-- vertex coordinates
634 * face_vtx_idx <-- "face -> vertices" connectivity index
635 * face_vtx <-- "face -> vertices" connectivity list
636 * face_cog <-> coordinates of the center of gravity of the faces
637 * face_norm <-- face surface normals
638 *
639 * Pi+1
640 * *---------* B : barycenter of the polygon
641 * / . . \
642 * / . . \ Pi : vertices of the polygon
643 * / . . \
644 * / . . Ti \ Ti : triangle
645 * *.........B.........* Pi
646 * Pn-1 \ . . /
647 * \ . . /
648 * \ . . /
649 * \ . T0 . /
650 * *---------*
651 * P0
652 *----------------------------------------------------------------------------*/
653
654 static void
_adjust_face_cog_v11_v52(cs_lnum_t n_faces,const cs_real_t vtx_coord[][3],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx[],cs_real_t face_cog[][3],const cs_real_t face_norm[][3])655 _adjust_face_cog_v11_v52(cs_lnum_t n_faces,
656 const cs_real_t vtx_coord[][3],
657 const cs_lnum_t face_vtx_idx[],
658 const cs_lnum_t face_vtx[],
659 cs_real_t face_cog[][3],
660 const cs_real_t face_norm[][3])
661 {
662 const cs_real_t one_third = 1./3.;
663
664 /* Loop on faces
665 --------------- */
666
667 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) {
668
669 cs_real_t tri_vol_part = 0.;
670
671 /* Define the polygon (P) according to the vertices (Pi) of the face */
672
673 cs_lnum_t s_id = face_vtx_idx[f_id];
674 cs_lnum_t e_id = face_vtx_idx[f_id + 1];
675
676 cs_lnum_t n_face_vertices = e_id - s_id;
677
678 /* No warping - related correction required for triangles*/
679 if (n_face_vertices < 4)
680 continue;
681
682 /* Compute approximate face center coordinates for the polygon */
683
684 cs_real_t a_center[3] = {0, 0, 0}, f_center[3] = {0, 0, 0};
685
686 for (cs_lnum_t j = s_id; j < e_id; j++) {
687 const cs_lnum_t v0 = face_vtx[j];
688 for (cs_lnum_t i = 0; i < 3; i++)
689 a_center[i] += vtx_coord[v0][i];
690 }
691
692 for (cs_lnum_t i = 0; i < 3; i++)
693 a_center[i] /= n_face_vertices;
694
695 /* loop on edges, with implied subdivision into triangles
696 defined by edge and cell center */
697
698 cs_real_t vc0[3], vc1[3], vn[3], vtc[3];
699 cs_real_t sum_w = 0;
700
701 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
702
703 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
704 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
705
706 for (cs_lnum_t i = 0; i < 3; i++) {
707 vc0[i] = vtx_coord[v0][i] - a_center[i];
708 vc1[i] = vtx_coord[v1][i] - a_center[i];
709 vtc[i] = vtx_coord[v1][i] + vtx_coord[v0][i] + a_center[i];
710 }
711
712 cs_math_3_cross_product(vc0, vc1, vn);
713
714 tri_vol_part += one_third * 0.5 * cs_math_3_dot_product(vtc, vn);
715
716 cs_real_t w = 0.5 * cs_math_3_norm(vn);
717
718 if (cs_math_3_dot_product(vn, face_norm[f_id]) < 0.0)
719 w *= -1.0;
720
721 sum_w += w;
722
723 for (cs_lnum_t i = 0; i < 3; i++)
724 f_center[i] += w * vtc[i];
725
726 } /* End of loop on triangles of the face */
727
728 /*
729 * Compute the center of gravity G(P) of the polygon P, and
730 * compute the part of volume of the polygon (before rectification)
731 *
732 * --> ->
733 * OG(P).N(P)
734 *------------ */
735
736 cs_real_t face_vol_part = 0.0;
737
738 for (cs_lnum_t i = 0; i < 3; i++) {
739 f_center[i] *= one_third / sum_w;
740 face_vol_part += (f_center[i] * face_norm[f_id][i]);
741 }
742
743 cs_real_t rectif_cog = (tri_vol_part - face_vol_part)
744 / (sum_w * sum_w);
745
746 for (cs_lnum_t i = 0; i < 3; i++)
747 face_cog[f_id][i] = f_center[i] + rectif_cog * face_norm[f_id][i];
748
749 } /* End of loop on faces */
750 }
751
752 /*----------------------------------------------------------------------------
753 * Refine face center computation based on initial value.
754 *
755 * This may be useful for warped faces; for plance faces, the face center
756 * should remain identical.
757 *
758 * parameters:
759 * n_faces <-- number of faces
760 * vtx_coord <-- vertex coordinates
761 * face_vtx_idx <-- "face -> vertices" connectivity index
762 * face_vtx <-- "face -> vertices" connectivity
763 * face_cog <-> coordinates of the center of gravity of the faces
764 * face_norm <-- face surface normals
765 *
766 * Pi+1
767 * *---------* B : barycenter of the polygon
768 * / . . \
769 * / . . \ Pi : vertices of the polygon
770 * / . . \
771 * / . . Ti \ Ti : triangle
772 * *.........B.........* Pi
773 * Pn-1 \ . . /
774 * \ . . /
775 * \ . . /
776 * \ . T0 . /
777 * *---------*
778 * P0
779 *----------------------------------------------------------------------------*/
780
781 static void
_refine_warped_face_centers(cs_lnum_t n_faces,const cs_real_3_t vtx_coord[],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx[],cs_real_t face_cog[][3],const cs_real_t face_norm[][3])782 _refine_warped_face_centers(cs_lnum_t n_faces,
783 const cs_real_3_t vtx_coord[],
784 const cs_lnum_t face_vtx_idx[],
785 const cs_lnum_t face_vtx[],
786 cs_real_t face_cog[][3],
787 const cs_real_t face_norm[][3])
788 {
789 /* Checking */
790
791 assert(face_cog != NULL || n_faces == 0);
792
793 const cs_real_t one_third = 1./3.;
794 const cs_real_t s_epsilon = 1.e-32; /* TODO define better "zero" threshold */
795
796 /* Loop on faces */
797
798 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) {
799
800 /* Define the polygon (P) according to the vertices (Pi) of the face */
801
802 cs_lnum_t s_id = face_vtx_idx[f_id];
803 cs_lnum_t e_id = face_vtx_idx[f_id + 1];
804
805 cs_lnum_t n_face_vertices = e_id - s_id;
806
807 if (n_face_vertices > 3) {
808
809 cs_real_t ref_size = sqrt(cs_math_3_norm(face_norm[f_id]));
810
811 for (int ite = 0; ite < 5; ite++) {
812
813 /* Use previous face center coordinates for the polygon */
814
815 cs_real_t a_center[3];
816 cs_real_t f_center[3] = {0, 0, 0};
817
818 for (cs_lnum_t i = 0; i < 3; i++)
819 a_center[i] = face_cog[f_id][i];
820
821 /* loop on edges, with implied subdivision into triangles
822 defined by edge and cell center */
823
824 cs_real_t vc0[3], vc1[3], vn[3], vtc[3];
825 cs_real_t sum_w = 0;
826
827 for (cs_lnum_t tri_id = 0; tri_id < n_face_vertices; tri_id++) {
828
829 const cs_lnum_t v0 = face_vtx[s_id + tri_id];
830 const cs_lnum_t v1 = face_vtx[s_id + (tri_id+1)%n_face_vertices];
831
832 for (cs_lnum_t i = 0; i < 3; i++) {
833 vc0[i] = vtx_coord[v0][i] - a_center[i];
834 vc1[i] = vtx_coord[v1][i] - a_center[i];
835 vtc[i] = vtx_coord[v1][i] + vtx_coord[v0][i] + a_center[i];
836 }
837
838 cs_math_3_cross_product(vc0, vc1, vn);
839
840 cs_real_t w = cs_math_3_norm(vn);
841
842 if (cs_math_3_dot_product(vn, face_norm[f_id]) < 0.0)
843 w *= -1.0;
844
845 sum_w += w;
846
847 for (cs_lnum_t i = 0; i < 3; i++)
848 f_center[i] += w*vtc[i];
849
850 }
851
852 if (sum_w > s_epsilon) {
853 for (cs_lnum_t i = 0; i < 3; i++)
854 face_cog[f_id][i] = one_third * f_center[i]/sum_w;
855 if ( cs_math_3_distance(face_cog[f_id], a_center)/ref_size < 1e-8)
856 break;
857 }
858 else
859 break;
860
861 } /* end of face iterations */
862
863 } /* end of test on triangle */
864
865 } /* end of loop on faces */
866 }
867
868 /*----------------------------------------------------------------------------
869 * Recompute quantities associated to faces (border or internal) when
870 * the quality of the mesh is not good enough
871 *----------------------------------------------------------------------------*/
872
873 static void
_correct_cell_face_center(const cs_mesh_t * mesh,const cs_lnum_t n_cells_with_ghosts,const cs_lnum_t n_i_faces,const cs_lnum_t n_b_faces,const cs_lnum_2_t i_face_cells[],const cs_lnum_t b_face_cells[],cs_real_3_t cell_cen[],cs_real_3_t i_face_cog[],cs_real_3_t b_face_cog[],cs_real_3_t i_face_normal[],cs_real_3_t b_face_normal[])874 _correct_cell_face_center(const cs_mesh_t *mesh,
875 const cs_lnum_t n_cells_with_ghosts,
876 const cs_lnum_t n_i_faces,
877 const cs_lnum_t n_b_faces,
878 const cs_lnum_2_t i_face_cells[],
879 const cs_lnum_t b_face_cells[],
880 cs_real_3_t cell_cen[],
881 cs_real_3_t i_face_cog[],
882 cs_real_3_t b_face_cog[],
883 cs_real_3_t i_face_normal[],
884 cs_real_3_t b_face_normal[])
885 {
886 int nitmax = 500;
887 cs_real_3_t *i_face_cog0, *b_face_cog0;
888 cs_real_3_t *i_face_cen, *b_face_cen;
889
890 cs_real_t *relaxf;
891 cs_real_t *relaxb;
892
893 cs_real_33_t *dxidxj;
894 cs_real_t *determinant;
895
896 BFT_MALLOC(i_face_cog0, n_i_faces, cs_real_3_t);
897 BFT_MALLOC(b_face_cog0, n_b_faces, cs_real_3_t);
898 BFT_MALLOC(i_face_cen, n_i_faces, cs_real_3_t);
899 BFT_MALLOC(b_face_cen, n_b_faces, cs_real_3_t);
900 BFT_MALLOC(relaxf, n_i_faces, cs_real_t);
901 BFT_MALLOC(relaxb, n_b_faces, cs_real_t);
902 BFT_MALLOC(dxidxj, n_cells_with_ghosts, cs_real_33_t);
903 BFT_MALLOC(determinant, n_cells_with_ghosts, cs_real_t);
904
905 /* Iterative process */
906 for (int sweep = 0; sweep < 1; sweep++) {
907
908 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++)
909 for (int i = 0; i < 3; i++)
910 i_face_cog0[face_id][i] = i_face_cog[face_id][i];
911
912 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++)
913 for (int i = 0; i < 3; i++)
914 b_face_cog0[face_id][i] = b_face_cog[face_id][i];
915
916 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
917 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
918 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
919
920 /* IF0 . S */
921 double ps1 = cs_math_3_distance_dot_product(cell_cen[cell_id1],
922 i_face_cog0[face_id],
923 i_face_normal[face_id]);
924
925 /* IJ . S */
926 double ps2 = cs_math_3_distance_dot_product(cell_cen[cell_id1],
927 cell_cen[cell_id2],
928 i_face_normal[face_id]);
929
930 double lambda = 0.5;
931 if (CS_ABS(ps2) > 1.e-20)
932 lambda = ps1 / ps2;
933
934 lambda = CS_MAX(lambda, 1./3.);
935 lambda = CS_MIN(lambda, 2./3.);
936
937 /* F = I + lambda * IJ */
938 for (int i = 0; i < 3; i++)
939 i_face_cen[face_id][i] = cell_cen[cell_id1][i]
940 + lambda * ( cell_cen[cell_id2][i]
941 - cell_cen[cell_id1][i]);
942 }
943
944 /* Compute the projection of I on the boundary face: b_face_cen */
945 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
946 cs_lnum_t cell_id = b_face_cells[face_id];
947
948 cs_real_3_t normal;
949 /* Normal is vector 0 if the b_face_normal norm is too small */
950 cs_math_3_normalise(b_face_normal[face_id], normal);
951
952 if (cs_math_3_norm(normal) > 0.) {
953 double lambda = cs_math_3_distance_dot_product(cell_cen[cell_id],
954 b_face_cog[face_id],
955 normal);
956
957 for (int i = 0; i < 3; i++)
958 b_face_cen[face_id][i] = cell_cen[cell_id][i] + lambda * normal[i];
959 } else {
960 for (int i = 0; i < 3; i++)
961 b_face_cen[face_id][i] = b_face_cog[face_id][i];
962 }
963 }
964
965 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++)
966 relaxf[face_id] = 1.;
967 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++)
968 relaxb[face_id] = 1.;
969
970 int iiter = 0;
971 cs_gnum_t irelax = 0;
972
973 do
974 {
975 iiter +=1;
976
977 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
978 for (int i = 0; i < 3; i++)
979 i_face_cog[face_id][i] = (1. - relaxf[face_id]) * i_face_cog0[face_id][i]
980 + relaxf[face_id] * i_face_cen[face_id][i];
981 }
982
983 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
984 for (int i = 0; i < 3; i++)
985 b_face_cog[face_id][i] = (1. - relaxb[face_id]) * b_face_cog0[face_id][i]
986 + relaxb[face_id] * b_face_cen[face_id][i];
987 }
988
989 for (cs_lnum_t cell_id = 0; cell_id < n_cells_with_ghosts; cell_id++)
990 for (int i = 0; i < 3; i++)
991 for (int j = 0; j < 3; j++)
992 dxidxj[cell_id][i][j] = 0.;
993
994 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
995 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
996 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
997
998 for (int i = 0; i < 3; i++)
999 for (int j = 0; j < 3; j++) {
1000 double fluxij = i_face_cog[face_id][i]//TODO minus celli
1001 * i_face_normal[face_id][j];
1002 dxidxj[cell_id1][i][j] += fluxij;
1003 dxidxj[cell_id2][i][j] -= fluxij;
1004 }
1005 }
1006
1007 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1008 cs_lnum_t cell_id = b_face_cells[face_id];
1009
1010 for (int i = 0; i < 3; i++)
1011 for (int j = 0; j < 3; j++) {
1012 double fluxij = b_face_cog[face_id][i]
1013 * b_face_normal[face_id][j];
1014 dxidxj[cell_id][i][j] += fluxij;
1015 }
1016 }
1017
1018 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
1019 double vol = ( dxidxj[cell_id][0][0]
1020 + dxidxj[cell_id][1][1]
1021 + dxidxj[cell_id][2][2] ) / 3.;
1022
1023 //FIXME
1024 if (vol >= 0)
1025 vol = CS_MAX(vol,1.e-20);
1026
1027 double a1 = dxidxj[cell_id][0][0] / vol;
1028 double a2 = dxidxj[cell_id][0][1] / vol;
1029 double a3 = dxidxj[cell_id][0][2] / vol;
1030 double a4 = dxidxj[cell_id][1][0] / vol;
1031 double a5 = dxidxj[cell_id][1][1] / vol;
1032 double a6 = dxidxj[cell_id][1][2] / vol;
1033 double a7 = dxidxj[cell_id][2][0] / vol;
1034 double a8 = dxidxj[cell_id][2][1] / vol;
1035 double a9 = dxidxj[cell_id][2][2] / vol;
1036
1037 determinant[cell_id] = fabs( a1 * (a5*a9 - a8*a6)
1038 - a2 * (a4*a9 - a7*a6)
1039 + a3 * (a4*a8 - a7*a5) );
1040
1041 //FIXME
1042 determinant[cell_id] = CS_MAX(determinant[cell_id],1.e-20);
1043 determinant[cell_id] = CS_MIN(determinant[cell_id],1./determinant[cell_id]);
1044
1045 /* M-matrix structure control */
1046 double mmatrice1a = a1 - CS_ABS(a2) - CS_ABS(a3);
1047 double mmatrice1b = a1 - CS_ABS(a4) - CS_ABS(a7);
1048 double mmatrice2a = a5 - CS_ABS(a4) - CS_ABS(a6);
1049 double mmatrice2b = a5 - CS_ABS(a2) - CS_ABS(a8);
1050 double mmatrice3a = a9 - CS_ABS(a7) - CS_ABS(a8);
1051 double mmatrice3b = a9 - CS_ABS(a6) - CS_ABS(a3);
1052
1053 if ( mmatrice1a <= 0. || mmatrice1b <= 0.
1054 || mmatrice2a <= 0. || mmatrice2b <= 0.
1055 || mmatrice3a <= 0. || mmatrice3b <= 0.)
1056 determinant[cell_id] = 0.;
1057 }
1058
1059 if (mesh->halo != NULL)
1060 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, determinant);
1061
1062 irelax = 0;
1063
1064 //FIXME test was 0.001
1065 cs_real_t threshold = 0.1;
1066 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1067 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1068 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1069
1070 if (determinant[cell_id1] < threshold || determinant[cell_id2] < threshold) {
1071 irelax +=1;
1072 relaxf[face_id] *= 0.95;
1073 }
1074 }
1075
1076 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1077 cs_lnum_t cell_id = b_face_cells[face_id];
1078
1079 if (determinant[cell_id] < threshold) {
1080 irelax += 1;
1081 relaxb[face_id] *= 0.95;
1082 }
1083 }
1084
1085 cs_parall_counter(&irelax, 1);
1086
1087 } while (iiter < nitmax && irelax > 0);
1088
1089 }
1090 BFT_FREE(i_face_cog0);
1091 BFT_FREE(b_face_cog0);
1092 BFT_FREE(i_face_cen);
1093 BFT_FREE(b_face_cen);
1094 BFT_FREE(relaxf);
1095 BFT_FREE(relaxb);
1096 BFT_FREE(dxidxj);
1097 BFT_FREE(determinant);
1098 }
1099
1100 /*----------------------------------------------------------------------------*/
1101 /*!
1102 * \brief Compute cell centers and volumes.
1103 *
1104 * \param[in] mesh pointer to mesh structure
1105 * \param[in] i_face_norm surface normal of internal faces
1106 * \param[in] i_face_cog center of gravity of internal faces
1107 * \param[in] b_face_norm surface normal of border faces
1108 * \param[in] b_face_cog center of gravity of border faces
1109 * \param[out] cell_cen cell centers
1110 * \param[out] cell_vol cell volumes
1111 */
1112 /*----------------------------------------------------------------------------*/
1113
1114 static void
_compute_cell_quantities(const cs_mesh_t * mesh,const cs_real_3_t i_face_norm[],const cs_real_3_t i_face_cog[],const cs_real_3_t b_face_norm[],const cs_real_3_t b_face_cog[],cs_real_3_t cell_cen[restrict],cs_real_t cell_vol[restrict])1115 _compute_cell_quantities(const cs_mesh_t *mesh,
1116 const cs_real_3_t i_face_norm[],
1117 const cs_real_3_t i_face_cog[],
1118 const cs_real_3_t b_face_norm[],
1119 const cs_real_3_t b_face_cog[],
1120 cs_real_3_t cell_cen[restrict],
1121 cs_real_t cell_vol[restrict])
1122 {
1123 /* Mesh connectivity */
1124
1125 const cs_lnum_t n_i_faces = mesh->n_i_faces;
1126 const cs_lnum_t n_b_faces = CS_MAX(mesh->n_b_faces, mesh->n_b_faces_all);
1127 const cs_lnum_t n_cells = mesh->n_cells;
1128 const cs_lnum_t n_cells_ext = mesh->n_cells_with_ghosts;
1129 const cs_lnum_2_t *i_face_cells
1130 = (const cs_lnum_2_t *)(mesh->i_face_cells);
1131 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
1132
1133 /* Checking */
1134
1135 assert(cell_cen != NULL);
1136 assert(cell_vol != NULL);
1137
1138 /* Compute approximate cell center using face centers */
1139
1140 cs_real_3_t *a_cell_cen;
1141 BFT_MALLOC(a_cell_cen, n_cells_ext, cs_real_3_t);
1142
1143 cs_mesh_quantities_cell_faces_cog(mesh,
1144 (const cs_real_t *)i_face_norm,
1145 (const cs_real_t *)i_face_cog,
1146 (const cs_real_t *)b_face_norm,
1147 (const cs_real_t *)b_face_cog,
1148 (cs_real_t *)a_cell_cen);
1149
1150 /* Initialization */
1151
1152 for (cs_lnum_t j = 0; j < n_cells_ext; j++)
1153 cell_vol[j] = 0.;
1154
1155 for (cs_lnum_t j = 0; j < n_cells_ext; j++) {
1156 for (cs_lnum_t i = 0; i < 3; i++)
1157 cell_cen[j][i] = 0.;
1158 }
1159
1160 /* Loop on interior faces
1161 ---------------------- */
1162
1163 for (cs_lnum_t f_id = 0; f_id < n_i_faces; f_id++) {
1164
1165 /* For each cell sharing the internal face, we update
1166 * cell_cen and cell_area */
1167
1168 cs_lnum_t c_id1 = i_face_cells[f_id][0];
1169 cs_lnum_t c_id2 = i_face_cells[f_id][1];
1170
1171 /* Implicit subdivision of cell into face vertices-cell-center pyramids */
1172
1173 if (c_id1 > -1) {
1174 cs_real_t pyra_vol_3 = cs_math_3_distance_dot_product(a_cell_cen[c_id1],
1175 i_face_cog[f_id],
1176 i_face_norm[f_id]);
1177 for (cs_lnum_t i = 0; i < 3; i++)
1178 cell_cen[c_id1][i] += pyra_vol_3 *( 0.75*i_face_cog[f_id][i]
1179 + 0.25*a_cell_cen[c_id1][i]);
1180 cell_vol[c_id1] += pyra_vol_3;
1181 }
1182 if (c_id2 > -1) {
1183 cs_real_t pyra_vol_3 = cs_math_3_distance_dot_product(i_face_cog[f_id],
1184 a_cell_cen[c_id2],
1185 i_face_norm[f_id]);
1186 for (cs_lnum_t i = 0; i < 3; i++)
1187 cell_cen[c_id2][i] += pyra_vol_3 *( 0.75*i_face_cog[f_id][i]
1188 + 0.25*a_cell_cen[c_id2][i]);
1189 cell_vol[c_id2] += pyra_vol_3;
1190 }
1191
1192 } /* End of loop on interior faces */
1193
1194 /* Loop on boundary faces
1195 --------------------- */
1196
1197 for (cs_lnum_t f_id = 0; f_id < n_b_faces; f_id++) {
1198
1199 /* For each cell sharing a border face, we update the numerator
1200 * of cell_cen and cell_area */
1201
1202 cs_lnum_t c_id1 = b_face_cells[f_id];
1203
1204 /* Computation of the area of the face
1205 (note that c_id1 == -1 may happen for isolated faces,
1206 which are cleaned afterwards) */
1207
1208 if (c_id1 > -1) {
1209 cs_real_t pyra_vol_3 = cs_math_3_distance_dot_product(a_cell_cen[c_id1],
1210 b_face_cog[f_id],
1211 b_face_norm[f_id]);
1212 for (cs_lnum_t i = 0; i < 3; i++)
1213 cell_cen[c_id1][i] += pyra_vol_3 *( 0.75*b_face_cog[f_id][i]
1214 + 0.25*a_cell_cen[c_id1][i]);
1215 cell_vol[c_id1] += pyra_vol_3;
1216 }
1217
1218 } /* End of loop on boundary faces */
1219
1220 BFT_FREE(a_cell_cen);
1221
1222 /* Loop on cells to finalize the computation
1223 ----------------------------------------- */
1224
1225 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
1226
1227 for (cs_lnum_t i = 0; i < 3; i++)
1228 cell_cen[c_id][i] /= cell_vol[c_id];
1229
1230 cell_vol[c_id] /= 3.0;
1231
1232 }
1233 }
1234
1235 /*----------------------------------------------------------------------------*
1236 * Compute new cell centers by minimizing the distance to faces
1237 *
1238 * parameters:
1239 * mesh <-- pointer to mesh structure
1240 * i_face_normal <-- surface normal of internal faces
1241 * i_face_cog <-- center of gravity of internal faces
1242 * b_face_normal <-- surface normal of border faces
1243 * b_face_cog <-- center of gravity of border faces
1244 * cell_cen --> center of gravity of cells
1245 *----------------------------------------------------------------------------*/
1246
1247 static void
_recompute_cell_cen_face(const cs_mesh_t * mesh,const cs_real_3_t i_face_normal[],const cs_real_3_t i_face_cog[],const cs_real_3_t b_face_normal[],const cs_real_3_t b_face_cog[],cs_real_3_t cell_cen[])1248 _recompute_cell_cen_face(const cs_mesh_t *mesh,
1249 const cs_real_3_t i_face_normal[],
1250 const cs_real_3_t i_face_cog[],
1251 const cs_real_3_t b_face_normal[],
1252 const cs_real_3_t b_face_cog[],
1253 cs_real_3_t cell_cen[])
1254 {
1255 const cs_lnum_t n_i_faces = mesh->n_i_faces;
1256 const cs_lnum_t n_b_faces = CS_MAX(mesh->n_b_faces, mesh->n_b_faces_all);
1257
1258 const cs_lnum_t n_cells_with_ghosts = mesh->n_cells_with_ghosts;
1259
1260 const cs_lnum_2_t *i_face_cells
1261 = (const cs_lnum_2_t *)(mesh->i_face_cells);
1262 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
1263
1264 /* First pass of verification */
1265 int *pb1;
1266 BFT_MALLOC(pb1, n_cells_with_ghosts, int);
1267
1268 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++)
1269 pb1[cell_id] = 0;
1270
1271 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1272
1273 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1274 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1275
1276 /* IF . S */
1277 double psi1 = cs_math_3_distance_dot_product(cell_cen[cell_id1],
1278 i_face_cog[face_id],
1279 i_face_normal[face_id]);
1280 /* JF . S */
1281 double psj1 = cs_math_3_distance_dot_product(cell_cen[cell_id2],
1282 i_face_cog[face_id],
1283 i_face_normal[face_id]);
1284 if (psi1 < 0.)
1285 pb1[cell_id1]++;
1286 if (psj1 > 0.)
1287 pb1[cell_id2]++;
1288 }
1289
1290 cs_gnum_t cpt1 = 0;
1291 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++)
1292 if (pb1[cell_id] > 0) cpt1++;
1293 cs_parall_counter(&cpt1, 1);
1294
1295 if (cpt1 > 0) {
1296 bft_printf("Total number of cell centers on the other side of a face "
1297 "(before correction) = %llu / %ld\n",
1298 (unsigned long long)cpt1,
1299 (long)mesh->n_cells);
1300
1301 /* Second pass */
1302 cs_real_33_t *a;
1303 cs_real_3_t *b;
1304 cs_real_3_t *cdgbis;
1305
1306 BFT_MALLOC(a, n_cells_with_ghosts, cs_real_33_t);
1307 BFT_MALLOC(b, n_cells_with_ghosts, cs_real_3_t);
1308 BFT_MALLOC(cdgbis, n_cells_with_ghosts, cs_real_3_t);
1309
1310 /* init matrice et second membre */
1311 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++) {
1312 for (int i = 0; i < 3; i++) {
1313 b[cell_id][i] = 0.;
1314 for (int j = 0; j < 3; j++)
1315 a[cell_id][i][j] = 0.;
1316 }
1317 }
1318
1319 /* Contribution from interior faces */
1320 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1321
1322 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1323 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1324 double surfn = cs_math_3_norm(i_face_normal[face_id]);
1325
1326 for (int i = 0; i < 3; i++)
1327 for (int j = 0; j < 3; j++) {
1328 a[cell_id1][i][j] += i_face_normal[face_id][i]
1329 * i_face_normal[face_id][j] / surfn;
1330 a[cell_id2][i][j] += i_face_normal[face_id][i]
1331 * i_face_normal[face_id][j] / surfn;
1332 }
1333
1334 double ps = cs_math_3_dot_product(i_face_normal[face_id],
1335 i_face_cog[face_id]);
1336
1337 for (int i = 0; i < 3; i++) {
1338 b[cell_id1][i] += ps * i_face_normal[face_id][i] / surfn;
1339 b[cell_id2][i] += ps * i_face_normal[face_id][i] / surfn;
1340 }
1341
1342 }
1343
1344 /* Contribution from boundary faces */
1345 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1346
1347 cs_lnum_t cell_id = b_face_cells[face_id];
1348 double surfn = cs_math_3_norm(b_face_normal[face_id]);
1349
1350 for (int i = 0; i < 3; i++)
1351 for (int j = 0; j < 3; j++) {
1352 a[cell_id][i][j] += b_face_normal[face_id][i]
1353 * b_face_normal[face_id][j] / surfn;
1354 }
1355
1356 double ps = cs_math_3_dot_product(b_face_normal[face_id],
1357 b_face_cog[face_id]);
1358
1359 for (int i = 0; i < 3; i++) {
1360 b[cell_id][i] += ps * b_face_normal[face_id][i] / surfn;
1361 }
1362
1363 }
1364
1365 /* invert system */
1366 double aainv[3][3];
1367 double bb[3];
1368 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
1369
1370 cdgbis[cell_id][0] = cell_cen[cell_id][0];
1371 cdgbis[cell_id][1] = cell_cen[cell_id][1];
1372 cdgbis[cell_id][2] = cell_cen[cell_id][2];
1373
1374 if (pb1[cell_id] > 0) {
1375
1376 double adim = a[cell_id][0][0] + a[cell_id][1][1] + a[cell_id][2][2];
1377
1378 if (adim > 0.) {
1379 bb[0] = b[cell_id][0] / adim;
1380 bb[1] = b[cell_id][1] / adim;
1381 bb[2] = b[cell_id][2] / adim;
1382
1383 /* Matrix inversion */
1384 double cocg11 = a[cell_id][0][0] / adim;
1385 double cocg12 = a[cell_id][0][1] / adim;
1386 double cocg13 = a[cell_id][0][2] / adim;
1387 double cocg21 = a[cell_id][1][0] / adim;
1388 double cocg22 = a[cell_id][1][1] / adim;
1389 double cocg23 = a[cell_id][1][2] / adim;
1390 double cocg31 = a[cell_id][2][0] / adim;
1391 double cocg32 = a[cell_id][2][1] / adim;
1392 double cocg33 = a[cell_id][2][2] / adim;
1393
1394 double a11 = cocg22 * cocg33 - cocg32 * cocg23;
1395 double a12 = cocg32 * cocg13 - cocg12 * cocg33;
1396 double a13 = cocg12 * cocg23 - cocg22 * cocg13;
1397 double a21 = cocg31 * cocg23 - cocg21 * cocg33;
1398 double a22 = cocg11 * cocg33 - cocg31 * cocg13;
1399 double a23 = cocg21 * cocg13 - cocg11 * cocg23;
1400 double a31 = cocg21 * cocg32 - cocg31 * cocg22;
1401 double a32 = cocg31 * cocg12 - cocg11 * cocg32;
1402 double a33 = cocg11 * cocg22 - cocg21 * cocg12;
1403
1404 double det_inv = cocg11 * a11 + cocg21 * a12 + cocg31 * a13;
1405
1406 if (CS_ABS(det_inv) >= 1.e-15) {
1407 det_inv = 1. / det_inv;
1408
1409 aainv[0][0] = a11 * det_inv;
1410 aainv[0][1] = a12 * det_inv;
1411 aainv[0][2] = a13 * det_inv;
1412 aainv[1][0] = a21 * det_inv;
1413 aainv[1][1] = a22 * det_inv;
1414 aainv[1][2] = a23 * det_inv;
1415 aainv[2][0] = a31 * det_inv;
1416 aainv[2][1] = a32 * det_inv;
1417 aainv[2][2] = a33 * det_inv;
1418
1419 for (int i = 0; i < 3; i++)
1420 cdgbis[cell_id][i] = aainv[i][0] * bb[0]
1421 + aainv[i][1] * bb[1]
1422 + aainv[i][2] * bb[2];
1423 }
1424 }
1425 }
1426 }
1427
1428 if (mesh->halo != NULL) {
1429 cs_halo_sync_var_strided(mesh->halo, CS_HALO_EXTENDED,
1430 (cs_real_t *)cdgbis, 3);
1431 if (mesh->n_init_perio > 0)
1432 cs_halo_perio_sync_coords(mesh->halo, CS_HALO_EXTENDED,
1433 (cs_real_t *)cdgbis);
1434 }
1435
1436 /* Second verification */
1437
1438 int *pb2;
1439 BFT_MALLOC(pb2, n_cells_with_ghosts, int);
1440
1441 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++)
1442 pb2[cell_id] = 0;
1443
1444 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1445
1446 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1447 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1448
1449 /* IF . S */
1450 double psi1 = cs_math_3_distance_dot_product(cdgbis[cell_id1],
1451 i_face_cog[face_id],
1452 i_face_normal[face_id]);
1453 /* JF . S */
1454 double psj1 = cs_math_3_distance_dot_product(cdgbis[cell_id2],
1455 i_face_cog[face_id],
1456 i_face_normal[face_id]);
1457 if (psi1 < 0.)
1458 pb2[cell_id1]++;
1459 if (psj1 > 0.)
1460 pb2[cell_id2]++;
1461 }
1462
1463 cs_gnum_t cpt2 = 0;
1464 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++)
1465 if (pb2[cell_id] > 0)
1466 cpt2++;
1467 cs_parall_counter(&cpt2, 1);
1468
1469 bft_printf("Total number of cell centers on the other side of a face "
1470 "(after correction) = %llu / %ld\n",
1471 (unsigned long long)cpt2,
1472 (long)mesh->n_cells);
1473
1474 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
1475 if (pb1[cell_id] > 0 && pb2[cell_id] == 0) {
1476
1477 cell_cen[cell_id][0] = cdgbis[cell_id][0];
1478 cell_cen[cell_id][1] = cdgbis[cell_id][1];
1479 cell_cen[cell_id][2] = cdgbis[cell_id][2];
1480 }
1481 }
1482
1483 if (mesh->halo != NULL) {
1484 cs_halo_sync_var_strided(mesh->halo, CS_HALO_EXTENDED,
1485 (cs_real_t *)cell_cen, 3);
1486 if (mesh->n_init_perio > 0)
1487 cs_halo_perio_sync_coords(mesh->halo, CS_HALO_EXTENDED,
1488 (cs_real_t *)cell_cen);
1489 }
1490
1491 BFT_FREE(a);
1492 BFT_FREE(b);
1493 BFT_FREE(cdgbis);
1494 BFT_FREE(pb2);
1495 }
1496
1497 /* Free memory */
1498 BFT_FREE(pb1);
1499 }
1500
1501 /*----------------------------------------------------------------------------
1502 * Compute the volume of cells C from their n faces F(i) and their center of
1503 * gravity G(Fi) where i=0, n-1
1504 *
1505 * 1 n-1
1506 * G(C) = - . Sum Surf(Fi) G(Fi)
1507 * 3 i=0
1508 *
1509 * parameters:
1510 * mesh <-- pointer to mesh structure
1511 * i_face_norm <-- surface normal of internal faces
1512 * i_face_cog <-- center of gravity of internal faces
1513 * b_face_norm <-- surface normal of border faces
1514 * b_face_cog <-- center of gravity of border faces
1515 * cell_cen <-- center of gravity of cells
1516 * cell_vol --> cells volume
1517 *----------------------------------------------------------------------------*/
1518
1519 static void
_compute_cell_volume(const cs_mesh_t * mesh,const cs_real_3_t i_face_norm[],const cs_real_3_t i_face_cog[],const cs_real_3_t b_face_norm[],const cs_real_3_t b_face_cog[],const cs_real_3_t cell_cen[],cs_real_t cell_vol[])1520 _compute_cell_volume(const cs_mesh_t *mesh,
1521 const cs_real_3_t i_face_norm[],
1522 const cs_real_3_t i_face_cog[],
1523 const cs_real_3_t b_face_norm[],
1524 const cs_real_3_t b_face_cog[],
1525 const cs_real_3_t cell_cen[],
1526 cs_real_t cell_vol[])
1527 {
1528 const cs_lnum_t n_b_faces = CS_MAX(mesh->n_b_faces, mesh->n_b_faces_all);
1529 const cs_real_t a_third = 1.0/3.0;
1530
1531 /* Initialization */
1532
1533 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++)
1534 cell_vol[cell_id] = 0;
1535
1536 /* Loop on internal faces */
1537
1538 for (cs_lnum_t fac_id = 0; fac_id < mesh->n_i_faces; fac_id++) {
1539
1540 cs_lnum_t cell_id1 = mesh->i_face_cells[fac_id][0];
1541 cs_lnum_t cell_id2 = mesh->i_face_cells[fac_id][1];
1542
1543 cell_vol[cell_id1] += cs_math_3_distance_dot_product(cell_cen[cell_id1],
1544 i_face_cog[fac_id],
1545 i_face_norm[fac_id]);
1546 cell_vol[cell_id2] -= cs_math_3_distance_dot_product(cell_cen[cell_id2],
1547 i_face_cog[fac_id],
1548 i_face_norm[fac_id]);
1549
1550 }
1551
1552 /* Loop on border faces */
1553
1554 for (cs_lnum_t fac_id = 0; fac_id < n_b_faces; fac_id++) {
1555
1556 cs_lnum_t cell_id1 = mesh->b_face_cells[fac_id];
1557
1558 cell_vol[cell_id1] += cs_math_3_distance_dot_product(cell_cen[cell_id1],
1559 b_face_cog[fac_id],
1560 b_face_norm[fac_id]);
1561 }
1562
1563 /* First computation of the volume */
1564
1565 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++)
1566 cell_vol[cell_id] *= a_third;
1567 }
1568
1569 /*----------------------------------------------------------------------------
1570 * Correction of small or negative volumes.
1571 *
1572 * Based on smoothing with neighbors; does not conserve the total volume.
1573 *----------------------------------------------------------------------------*/
1574
1575 static void
_cell_bad_volume_correction(const cs_mesh_t * mesh,cs_real_t cell_vol[])1576 _cell_bad_volume_correction(const cs_mesh_t *mesh,
1577 cs_real_t cell_vol[])
1578 {
1579 if (mesh->halo != NULL)
1580 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, cell_vol);
1581
1582 /* Iterations in order to get vol_I / max(vol_J) > critmin */
1583
1584 double *vol_neib_max;
1585 BFT_MALLOC(vol_neib_max, mesh->n_cells_with_ghosts, double);
1586
1587 for (int iter = 0; iter < 10; iter++) {
1588
1589 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++)
1590 vol_neib_max[cell_id] = 0.;
1591
1592 for (cs_lnum_t fac_id = 0; fac_id < mesh->n_i_faces; fac_id++) {
1593
1594 cs_lnum_t cell_id1 = mesh->i_face_cells[fac_id][0];
1595 cs_lnum_t cell_id2 = mesh->i_face_cells[fac_id][1];
1596 double vol1 = cell_vol[cell_id1];
1597 double vol2 = cell_vol[cell_id2];
1598
1599 if (vol2 > 0.)
1600 vol_neib_max[cell_id1] = CS_MAX(vol_neib_max[cell_id1], vol2);
1601
1602 if (vol1 > 0.)
1603 vol_neib_max[cell_id2] = CS_MAX(vol_neib_max[cell_id2], vol1);
1604 }
1605
1606 /* Previous value of 0.2 sometimes leads to computation divergence */
1607 /* 0.01 seems better and safer for the moment */
1608 double critmin = 0.01;
1609
1610 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++)
1611 cell_vol[cell_id] = CS_MAX(cell_vol[cell_id],
1612 critmin * vol_neib_max[cell_id]);
1613
1614 if (mesh->halo != NULL)
1615 cs_halo_sync_var(mesh->halo, CS_HALO_STANDARD, cell_vol);
1616 }
1617
1618 BFT_FREE(vol_neib_max);
1619 }
1620
1621 /*----------------------------------------------------------------------------
1622 * Compute the total, min, and max volumes of cells.
1623 *
1624 * parameters:
1625 * mesh <-- pointer to mesh structure
1626 * cell_vol <-- cells volume
1627 * min_vol --> minimum control volume
1628 * max_vol --> maximum control volume
1629 * tot_vol --> total control volume
1630 *----------------------------------------------------------------------------*/
1631
1632 static void
_cell_volume_reductions(const cs_mesh_t * mesh,cs_real_t cell_vol[],cs_real_t * min_vol,cs_real_t * max_vol,cs_real_t * tot_vol)1633 _cell_volume_reductions(const cs_mesh_t *mesh,
1634 cs_real_t cell_vol[],
1635 cs_real_t *min_vol,
1636 cs_real_t *max_vol,
1637 cs_real_t *tot_vol)
1638 {
1639 *min_vol = cs_math_infinite_r;
1640 *max_vol = -cs_math_infinite_r;
1641 *tot_vol = 0.;
1642
1643 for (cs_lnum_t cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
1644 *min_vol = CS_MIN(*min_vol, cell_vol[cell_id]);
1645 *max_vol = CS_MAX(*max_vol, cell_vol[cell_id]);
1646 *tot_vol = *tot_vol + cell_vol[cell_id];
1647 }
1648 }
1649
1650 /*----------------------------------------------------------------------------
1651 * Compute some distances relative to faces and associated weighting.
1652 *
1653 * parameters:
1654 * n_i_faces <-- number of interior faces
1655 * n_b_faces <-- number of border faces
1656 * i_face_cells <-- interior "faces -> cells" connectivity
1657 * b_face_cells <-- border "faces -> cells" connectivity
1658 * i_face_norm <-- surface normal of interior faces
1659 * b_face_norm <-- surface normal of border faces
1660 * i_face_cog <-- center of gravity of interior faces
1661 * b_face_cog <-- center of gravity of border faces
1662 * cell_cen <-- cell center
1663 * cell_vol <-- cell volume
1664 * i_dist --> distance IJ.Nij for interior faces
1665 * b_dist --> likewise for border faces
1666 * weight --> weighting factor (Aij=pond Ai+(1-pond)Aj)
1667 *----------------------------------------------------------------------------*/
1668
1669 static void
_compute_face_distances(cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_2_t i_face_cells[],const cs_lnum_t b_face_cells[],const cs_real_t i_face_normal[][3],const cs_real_t b_face_normal[][3],const cs_real_t i_face_cog[][3],const cs_real_t b_face_cog[][3],const cs_real_t cell_cen[][3],const cs_real_t cell_vol[],cs_real_t i_dist[],cs_real_t b_dist[],cs_real_t weight[])1670 _compute_face_distances(cs_lnum_t n_i_faces,
1671 cs_lnum_t n_b_faces,
1672 const cs_lnum_2_t i_face_cells[],
1673 const cs_lnum_t b_face_cells[],
1674 const cs_real_t i_face_normal[][3],
1675 const cs_real_t b_face_normal[][3],
1676 const cs_real_t i_face_cog[][3],
1677 const cs_real_t b_face_cog[][3],
1678 const cs_real_t cell_cen[][3],
1679 const cs_real_t cell_vol[],
1680 cs_real_t i_dist[],
1681 cs_real_t b_dist[],
1682 cs_real_t weight[])
1683 {
1684 cs_gnum_t w_count = 0;
1685
1686 /* Interior faces */
1687
1688 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1689
1690 const cs_real_t *face_nomal = i_face_normal[face_id];
1691 cs_real_t normal[3];
1692 cs_math_3_normalise(face_nomal, normal);
1693
1694 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1695 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1696
1697 /* Distance between the neighbor cell centers
1698 * and dot-product with the normal */
1699 i_dist[face_id] = cs_math_3_distance_dot_product(cell_cen[cell_id1],
1700 cell_cen[cell_id2],
1701 normal);
1702
1703 if (CS_ABS(i_dist[face_id]) > 1e-20) {
1704 /* Distance between the face center of gravity
1705 and the neighbor cell center
1706 and dot-product with the normal */
1707 cs_real_t dist2f = cs_math_3_distance_dot_product(i_face_cog[face_id],
1708 cell_cen[cell_id2],
1709 normal);
1710 weight[face_id] = dist2f / i_dist[face_id];
1711 }
1712 else {
1713 weight[face_id] = 0.5;
1714 }
1715
1716 /* Clipping of cell cell distances */
1717 if (cs_glob_mesh_quantities_flag & CS_FACE_DISTANCE_CLIP) {
1718
1719 /* Min value between IJ and
1720 * (Omega_i+Omega_j)/S_ij which is exactly the distance for tetras */
1721 cs_real_t distmax
1722 = cs_math_fmin(cs_math_3_distance(cell_cen[cell_id1],
1723 cell_cen[cell_id2]),
1724 ( (cell_vol[cell_id1] + cell_vol[cell_id2])
1725 / cs_math_3_norm(face_nomal)));
1726
1727 /* Previous value of 0.2 sometimes leads to computation divergence */
1728 /* 0.01 seems better and safer for the moment */
1729 double critmin = 0.01;
1730 if (i_dist[face_id] < critmin * distmax) {
1731 w_count++;
1732 i_dist[face_id] = cs_math_fmax(i_dist[face_id], critmin * distmax);
1733 }
1734
1735 /* Clipping of weighting */
1736 weight[face_id] = cs_math_fmax(weight[face_id], 0.001);
1737 weight[face_id] = cs_math_fmin(weight[face_id], 0.999);
1738 }
1739 }
1740
1741 cs_parall_counter(&w_count, 1);
1742
1743 if (w_count > 0)
1744 bft_printf(_("\n"
1745 "%llu faces have a too small distance between centers.\n"
1746 "For these faces, the weight may be clipped.\n"),
1747 (unsigned long long)w_count);
1748
1749 /* Boundary faces */
1750
1751 w_count = 0;
1752
1753 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1754
1755 const cs_real_t *face_nomal = b_face_normal[face_id];
1756 cs_real_t normal[3];
1757 cs_math_3_normalise(face_nomal, normal);
1758
1759 cs_lnum_t cell_id = b_face_cells[face_id];
1760
1761 /* Distance between the face center of gravity
1762 and the neighbor cell center */
1763 b_dist[face_id] = cs_math_3_distance_dot_product(cell_cen[cell_id],
1764 b_face_cog[face_id],
1765 normal);
1766 /* Clipping of cell boundary distances */
1767 if (cs_glob_mesh_quantities_flag & CS_FACE_DISTANCE_CLIP) {
1768
1769 /* Min value between IF and
1770 * (Omega_i)/S which is exactly the distance for tetras */
1771 double distmax = CS_MIN(
1772 cs_math_3_distance(cell_cen[cell_id], b_face_cog[face_id]),
1773 cell_vol[cell_id]/cs_math_3_norm(face_nomal));
1774
1775 double critmin = 0.01;
1776 if (b_dist[face_id] < critmin * distmax) {
1777 w_count++;
1778 b_dist[face_id] = CS_MAX(b_dist[face_id], critmin * distmax);
1779 }
1780
1781 }
1782 }
1783
1784 cs_parall_counter(&w_count, 1);
1785
1786 if (w_count > 0)
1787 bft_printf(_("\n"
1788 "%llu boundary faces have a too small distance between\n"
1789 "cell centre and face centre.\n"),
1790 (unsigned long long)w_count);
1791 }
1792
1793 /*----------------------------------------------------------------------------
1794 * Compute some vectors to handle non-orthogonalities.
1795 *
1796 * Let a face and I, J the centers of neighboring cells
1797 * (only I is defined for a border face)
1798 *
1799 * The face is oriented from I to J, with Nij its normal.
1800 * (border faces are oriented towards the exterior)
1801 * The norm of Nij is 1.
1802 * The face surface is Sij.
1803 *
1804 * I' and J' are defined as the orthogonal projection of I and J on the line
1805 * orthogonal to the face passing through the center of gravity F of the face.
1806 * (only I' is defined for a border face)
1807 *
1808 * We compute here the vector I'J' for interior faces (dijpf)
1809 * the vector II' for border faces (diipb)
1810 * the vector OF for interior faces (dofij)
1811 *
1812 * We also have the following formulae
1813 * II' = IG - (IG.Nij)Nij
1814 * JJ' = JG - (JG.Nij)Nij
1815 *
1816 * parameters:
1817 * dim <-- dimension
1818 * n_i_faces <-- number of interior faces
1819 * n_b_faces <-- number of border faces
1820 * i_face_cells <-- interior "faces -> cells" connectivity
1821 * b_face_cells <-- border "faces -> cells" connectivity
1822 * i_face_norm <-- surface normal of interior faces
1823 * b_face_norm <-- surface normal of border faces
1824 * i_face_cog <-- center of gravity of interior faces
1825 * b_face_cog <-- center of gravity of border faces
1826 * i_face_surf <-- interior faces surface
1827 * cell_cen <-- cell center
1828 * weight <-- weighting factor (Aij=pond Ai+(1-pond)Aj)
1829 * dijpf --> vector i'j' for interior faces
1830 * diipb --> vector ii' for border faces
1831 * dofij --> vector OF for interior faces
1832 *----------------------------------------------------------------------------*/
1833
1834 static void
_compute_face_vectors(int dim,const cs_lnum_t n_i_faces,const cs_lnum_t n_b_faces,const cs_lnum_2_t i_face_cells[],const cs_lnum_t b_face_cells[],const cs_real_t i_face_normal[],const cs_real_t b_face_normal[],const cs_real_t i_face_cog[],const cs_real_t b_face_cog[],const cs_real_t i_face_surf[],const cs_real_t cell_cen[],const cs_real_t weight[],const cs_real_t b_dist[],cs_real_t dijpf[],cs_real_t diipb[],cs_real_t dofij[])1835 _compute_face_vectors(int dim,
1836 const cs_lnum_t n_i_faces,
1837 const cs_lnum_t n_b_faces,
1838 const cs_lnum_2_t i_face_cells[],
1839 const cs_lnum_t b_face_cells[],
1840 const cs_real_t i_face_normal[],
1841 const cs_real_t b_face_normal[],
1842 const cs_real_t i_face_cog[],
1843 const cs_real_t b_face_cog[],
1844 const cs_real_t i_face_surf[],
1845 const cs_real_t cell_cen[],
1846 const cs_real_t weight[],
1847 const cs_real_t b_dist[],
1848 cs_real_t dijpf[],
1849 cs_real_t diipb[],
1850 cs_real_t dofij[])
1851 {
1852 cs_lnum_t face_id;
1853 cs_lnum_t cell_id;
1854
1855 cs_real_t dipjp, pond;
1856 cs_real_t surfnx, surfny, surfnz;
1857 cs_real_t vecijx, vecijy, vecijz;
1858
1859 /* Interior faces */
1860
1861 for (face_id = 0; face_id < n_i_faces; face_id++) {
1862
1863 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
1864 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
1865
1866 /* Normalized normal */
1867 surfnx = i_face_normal[face_id*dim] / i_face_surf[face_id];
1868 surfny = i_face_normal[face_id*dim + 1] / i_face_surf[face_id];
1869 surfnz = i_face_normal[face_id*dim + 2] / i_face_surf[face_id];
1870
1871 /* ---> IJ */
1872 vecijx = cell_cen[cell_id2*dim] - cell_cen[cell_id1*dim];
1873 vecijy = cell_cen[cell_id2*dim + 1] - cell_cen[cell_id1*dim + 1];
1874 vecijz = cell_cen[cell_id2*dim + 2] - cell_cen[cell_id1*dim + 2];
1875
1876 /* ---> DIJPP = IJ.NIJ */
1877 dipjp = vecijx*surfnx + vecijy*surfny + vecijz*surfnz;
1878
1879 /* ---> DIJPF = (IJ.NIJ).NIJ */
1880 dijpf[face_id*dim] = dipjp*surfnx;
1881 dijpf[face_id*dim + 1] = dipjp*surfny;
1882 dijpf[face_id*dim + 2] = dipjp*surfnz;
1883
1884 pond = weight[face_id];
1885
1886 /* ---> DOFIJ = OF */
1887 dofij[face_id*dim] = i_face_cog[face_id*dim]
1888 - ( pond *cell_cen[cell_id1*dim]
1889 + (1. - pond)*cell_cen[cell_id2*dim]);
1890
1891 dofij[face_id*dim + 1] = i_face_cog[face_id*dim + 1]
1892 - ( pond *cell_cen[cell_id1*dim + 1]
1893 + (1. - pond)*cell_cen[cell_id2*dim + 1]);
1894
1895 dofij[face_id*dim + 2] = i_face_cog[face_id*dim + 2]
1896 - ( pond *cell_cen[cell_id1*dim + 2]
1897 + (1. - pond)*cell_cen[cell_id2*dim + 2]);
1898 }
1899
1900 /* Boundary faces */
1901 cs_gnum_t w_count = 0;
1902
1903 for (face_id = 0; face_id < n_b_faces; face_id++) {
1904
1905 cell_id = b_face_cells[face_id];
1906
1907 cs_real_3_t normal;
1908 /* Normal is vector 0 if the b_face_normal norm is too small */
1909 cs_math_3_normalise(&b_face_normal[face_id*dim], normal);
1910
1911 /* ---> IF */
1912 cs_real_t vec_if[3] = {
1913 b_face_cog[face_id*dim] - cell_cen[cell_id*dim],
1914 b_face_cog[face_id*dim + 1] - cell_cen[cell_id*dim + 1],
1915 b_face_cog[face_id*dim + 2] - cell_cen[cell_id*dim + 2]};
1916
1917 /* ---> diipb = IF - (IF.NIJ)NIJ */
1918 cs_math_3_orthogonal_projection(normal, vec_if, &diipb[face_id*dim]);
1919
1920 /* Limiter on boundary face reconstruction */
1921 if (cs_glob_mesh_quantities_flag & CS_FACE_RECONSTRUCTION_CLIP) {
1922 cs_real_t iip = cs_math_3_norm(&diipb[face_id*dim]);
1923
1924 bool is_clipped = false;
1925 cs_real_t corri = 1.;
1926
1927 if (iip > 0.5 * b_dist[face_id]) {
1928 is_clipped = true;
1929 corri = 0.5 * b_dist[face_id] / iip;
1930 }
1931
1932 diipb[face_id*dim] *= corri;
1933 diipb[face_id*dim +1] *= corri;
1934 diipb[face_id*dim +2] *= corri;
1935
1936 if (is_clipped)
1937 w_count++;
1938 }
1939 }
1940
1941 cs_parall_counter(&w_count, 1);
1942
1943 if (w_count > 0)
1944 bft_printf(_("\n"
1945 "%llu boundary faces have a too large reconstruction distance.\n"
1946 "For these faces, reconstruction are limited.\n"),
1947 (unsigned long long)w_count);
1948
1949 }
1950
1951 /*----------------------------------------------------------------------------
1952 * Compute some vectors to handle non-orthogonalities.
1953 *
1954 * Let a face and I, J the centers of neighboring cells
1955 * (only I is defined for a border face)
1956 *
1957 * The face is oriented from I to J, with Nij its normal.
1958 * (border faces are oriented towards the exterior)
1959 * The norm of Nij is 1.
1960 * The face surface is Sij.
1961 *
1962 * I' and J' are defined as the orthogonal projection of I and J on the line
1963 * orthogonal to the face passing through the center of gravity F of the face.
1964 * (only I' is defined for a border face)
1965 *
1966 * We compute here the vector II' for interior faces (diipf)
1967 * the vector JJ' for interior faces (djjpf)
1968 *
1969 * We also have the following formulae
1970 * II' = IF - (IF.Nij)Nij
1971 * JJ' = JF - (JF.Nij)Nij
1972 *
1973 * parameters:
1974 * n_cells <-- number of cells
1975 * n_i_faces <-- number of interior faces
1976 * i_face_cells <-- interior "faces -> cells" connectivity
1977 * i_face_norm <-- surface normal of interior faces
1978 * i_face_cog <-- center of gravity of interior faces
1979 * i_face_surf <-- interior faces surface
1980 * cell_cen <-- cell center
1981 * cell_vol <-- cell volume
1982 * dist <-- interior distance
1983 * diipf --> vector ii' for interior faces
1984 * djjpf --> vector jj' for interior faces
1985 *----------------------------------------------------------------------------*/
1986
1987 static void
_compute_face_sup_vectors(const cs_lnum_t n_cells,const cs_lnum_t n_i_faces,const cs_lnum_2_t i_face_cells[],const cs_real_t i_face_normal[][3],const cs_real_t i_face_cog[][3],const cs_real_t cell_cen[][3],const cs_real_t cell_vol[],const cs_real_t dist[],cs_real_t diipf[][3],cs_real_t djjpf[][3])1988 _compute_face_sup_vectors(const cs_lnum_t n_cells,
1989 const cs_lnum_t n_i_faces,
1990 const cs_lnum_2_t i_face_cells[],
1991 const cs_real_t i_face_normal[][3],
1992 const cs_real_t i_face_cog[][3],
1993 const cs_real_t cell_cen[][3],
1994 const cs_real_t cell_vol[],
1995 const cs_real_t dist[],
1996 cs_real_t diipf[][3],
1997 cs_real_t djjpf[][3])
1998 {
1999 cs_gnum_t w_count = 0;
2000
2001 /* Interior faces */
2002
2003 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
2004
2005 cs_lnum_t cell_id1 = i_face_cells[face_id][0];
2006 cs_lnum_t cell_id2 = i_face_cells[face_id][1];
2007
2008 /* Normalized normal */
2009 cs_real_3_t normal;
2010 cs_math_3_normalise(i_face_normal[face_id], normal);
2011
2012 /* ---> IF and JF */
2013 cs_real_t vec_if[3] = {
2014 i_face_cog[face_id][0] - cell_cen[cell_id1][0],
2015 i_face_cog[face_id][1] - cell_cen[cell_id1][1],
2016 i_face_cog[face_id][2] - cell_cen[cell_id1][2]};
2017
2018 cs_real_t vec_jf[3] = {
2019 i_face_cog[face_id][0] - cell_cen[cell_id2][0],
2020 i_face_cog[face_id][1] - cell_cen[cell_id2][1],
2021 i_face_cog[face_id][2] - cell_cen[cell_id2][2]};
2022
2023 /* ---> diipf = IF - (IF.Nij)Nij */
2024 cs_math_3_orthogonal_projection(normal, vec_if, diipf[face_id]);
2025
2026 /* ---> djjpf = JF - (JF.Nij)Nij */
2027 cs_math_3_orthogonal_projection(normal, vec_jf, djjpf[face_id]);
2028
2029 /* Limiter on interior face reconstruction */
2030 if (cs_glob_mesh_quantities_flag & CS_FACE_RECONSTRUCTION_CLIP) {
2031
2032 cs_real_t surfn = cs_math_3_norm(i_face_normal[face_id]);
2033 cs_real_t iip = cs_math_3_norm(diipf[face_id]);
2034
2035 cs_real_t corri = 1.;
2036 bool is_clipped = false;
2037
2038 if (0.5 * dist[face_id] < iip) {
2039 is_clipped = true;
2040 corri = 0.5 * dist[face_id] / iip;
2041 }
2042
2043 diipf[face_id][0] *= corri;
2044 diipf[face_id][1] *= corri;
2045 diipf[face_id][2] *= corri;
2046
2047 iip = cs_math_3_norm(diipf[face_id]);
2048
2049 corri = 1.;
2050
2051 if (0.9 * cell_vol[cell_id1] < surfn * iip) {
2052 is_clipped = true;
2053 corri = 0.9 * cell_vol[cell_id1] / (surfn * iip);
2054 }
2055
2056 diipf[face_id][0] *= corri;
2057 diipf[face_id][1] *= corri;
2058 diipf[face_id][2] *= corri;
2059
2060 cs_real_t jjp = cs_math_3_norm(djjpf[face_id]);
2061
2062 cs_real_t corrj = 1.;
2063
2064 if (0.5 * dist[face_id] < jjp) {
2065 is_clipped = true;
2066 corrj = 0.5 * dist[face_id] / jjp;
2067 }
2068
2069 djjpf[face_id][0] *= corrj;
2070 djjpf[face_id][1] *= corrj;
2071 djjpf[face_id][2] *= corrj;
2072
2073 jjp = cs_math_3_norm(djjpf[face_id]);
2074
2075 corrj = 1.;
2076 if (0.9 * cell_vol[cell_id2] < surfn * jjp) {
2077 is_clipped = true;
2078 corrj = 0.9 * cell_vol[cell_id2] / (surfn * jjp);
2079 }
2080
2081 if (is_clipped && cell_id1 < n_cells)
2082 w_count++;
2083
2084 djjpf[face_id][0] *= corrj;
2085 djjpf[face_id][1] *= corrj;
2086 djjpf[face_id][2] *= corrj;
2087
2088 }
2089 }
2090
2091 cs_parall_counter(&w_count, 1);
2092
2093 if (w_count > 0)
2094 bft_printf(_("\n"
2095 "%llu internal faces have a too large reconstruction distance.\n"
2096 "For these faces, reconstruction are limited.\n"
2097 "\n"),
2098 (unsigned long long)w_count);
2099
2100 }
2101
2102 /*----------------------------------------------------------------------------
2103 * Evaluate boundary thickness.
2104 *
2105 * parameters:
2106 * m <-- pointer to mesh structure
2107 * m_quantities <-- pointer to mesh quantities structures.
2108 * b_thickness --> boundary thickness
2109 *----------------------------------------------------------------------------*/
2110
2111 static void
_b_thickness(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,cs_real_t b_thickness[])2112 _b_thickness(const cs_mesh_t *m,
2113 const cs_mesh_quantities_t *mq,
2114 cs_real_t b_thickness[])
2115 {
2116 const cs_real_3_t *cell_cen
2117 = (const cs_real_3_t *)(mq->cell_cen);
2118 const cs_real_3_t *b_face_cog
2119 = (const cs_real_3_t *)(mq->b_face_cog);
2120 const cs_real_3_t *b_face_normal
2121 = (const cs_real_3_t *)(mq->b_face_normal);
2122 const cs_real_t *b_face_surf
2123 = (const cs_real_t *)(mq->b_face_surf);
2124
2125 for (cs_lnum_t f_id = 0; f_id < m->n_b_faces; f_id++) {
2126 cs_lnum_t c_id = m->b_face_cells[f_id];
2127 b_thickness[f_id]
2128 = ( (b_face_cog[f_id][0] - cell_cen[c_id][0])*b_face_normal[f_id][0]
2129 + (b_face_cog[f_id][1] - cell_cen[c_id][1])*b_face_normal[f_id][1]
2130 + (b_face_cog[f_id][2] - cell_cen[c_id][2])*b_face_normal[f_id][2])
2131 * 2.0 / b_face_surf[f_id];
2132 }
2133 }
2134
2135 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2136
2137 /*============================================================================
2138 * Fortran wrapper function definitions
2139 *============================================================================*/
2140
2141 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
2142
2143 /*----------------------------------------------------------------------------
2144 * Compute the total, min, and max fluid volumes of cells
2145 *----------------------------------------------------------------------------*/
2146
2147 void
cs_f_mesh_quantities_fluid_vol_reductions(void)2148 cs_f_mesh_quantities_fluid_vol_reductions(void)
2149 {
2150 const cs_mesh_t *m = cs_glob_mesh;
2151 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
2152
2153 cs_mesh_quantities_fluid_vol_reductions(m,
2154 mq);
2155 }
2156
2157 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2158
2159 /*=============================================================================
2160 * Public function definitions
2161 *============================================================================*/
2162
2163 /*----------------------------------------------------------------------------*/
2164 /*!
2165 * \brief Query or modification of the option for computing cell centers.
2166 *
2167 * \param[in] algo_choice < 0 : query
2168 * 0 : computation based on face centers (default)
2169 * 1 : computation by cell sub-volumes
2170 *
2171 * \return 0 or 1 according to the selected algorithm
2172 */
2173 /*----------------------------------------------------------------------------*/
2174
2175 int
cs_mesh_quantities_cell_cen_choice(int algo_choice)2176 cs_mesh_quantities_cell_cen_choice(int algo_choice)
2177 {
2178 if (algo_choice > -1 && algo_choice < 2)
2179 _cell_cen_algorithm = algo_choice;
2180
2181 return _cell_cen_algorithm;
2182 }
2183
2184 /*----------------------------------------------------------------------------*/
2185 /*!
2186 * \brief Query or modification of the option for computing face centers.
2187 *
2188 * \param[in] algo_choice < 0 : query
2189 * 0 : standard computation
2190 * 1 : use adjustment for volume
2191 * from versions 1.1 to 5.3
2192 *
2193 * \return 0 or 1 according to the selected algorithm
2194 */
2195 /*----------------------------------------------------------------------------*/
2196
2197 int
cs_mesh_quantities_face_cog_choice(int algo_choice)2198 cs_mesh_quantities_face_cog_choice(int algo_choice)
2199 {
2200 if (algo_choice > -1 && algo_choice < 2)
2201 _ajust_face_cog_compat_v11_v52 = algo_choice;
2202
2203 return _ajust_face_cog_compat_v11_v52;
2204 }
2205
2206 /*----------------------------------------------------------------------------*/
2207 /*!
2208 * \brief Create a mesh quantities structure.
2209 *
2210 * \return pointer to created cs_mesh_quantities_t structure
2211 */
2212 /*----------------------------------------------------------------------------*/
2213
2214 cs_mesh_quantities_t *
cs_mesh_quantities_create(void)2215 cs_mesh_quantities_create(void)
2216 {
2217 cs_mesh_quantities_t *mesh_quantities = NULL;
2218
2219 BFT_MALLOC(mesh_quantities, 1, cs_mesh_quantities_t);
2220
2221 mesh_quantities->cell_cen = NULL;
2222 mesh_quantities->cell_vol = NULL;
2223 mesh_quantities->cell_f_vol = NULL;
2224 mesh_quantities->i_face_normal = NULL;
2225 mesh_quantities->b_face_normal = NULL;
2226 mesh_quantities->i_f_face_normal = NULL;
2227 mesh_quantities->b_f_face_normal = NULL;
2228 mesh_quantities->i_face_cog = NULL;
2229 mesh_quantities->b_face_cog = NULL;
2230 mesh_quantities->i_face_surf = NULL;
2231 mesh_quantities->b_face_surf = NULL;
2232 mesh_quantities->i_f_face_surf = NULL;
2233 mesh_quantities->b_f_face_surf = NULL;
2234 mesh_quantities->i_f_face_factor = NULL;
2235 mesh_quantities->b_f_face_factor = NULL;
2236 mesh_quantities->i_dist = NULL;
2237 mesh_quantities->b_dist = NULL;
2238 mesh_quantities->weight = NULL;
2239 mesh_quantities->dijpf = NULL;
2240 mesh_quantities->diipb = NULL;
2241 mesh_quantities->dofij = NULL;
2242 mesh_quantities->diipf = NULL;
2243 mesh_quantities->djjpf = NULL;
2244 mesh_quantities->corr_grad_lin_det = NULL;
2245 mesh_quantities->corr_grad_lin = NULL;
2246 mesh_quantities->b_sym_flag = NULL;
2247 mesh_quantities->has_disable_flag = 0;
2248 mesh_quantities->c_disable_flag = NULL;
2249 mesh_quantities->bad_cell_flag = NULL;
2250
2251 return (mesh_quantities);
2252 }
2253
2254 /*----------------------------------------------------------------------------*/
2255 /*!
2256 * \brief Destroy a mesh quantities structure.
2257 *
2258 * \param[in] mq pointer to mesh quantities structures
2259 *
2260 * \return NULL
2261 */
2262 /*----------------------------------------------------------------------------*/
2263
2264 cs_mesh_quantities_t *
cs_mesh_quantities_destroy(cs_mesh_quantities_t * mq)2265 cs_mesh_quantities_destroy(cs_mesh_quantities_t *mq)
2266 {
2267 cs_mesh_quantities_free_all(mq);
2268
2269 BFT_FREE(mq);
2270
2271 return (mq);
2272 }
2273
2274 /*----------------------------------------------------------------------------*/
2275 /*!
2276 * \brief Reset a mesh quantities structure to its empty initial state.
2277 *
2278 * \param[in] mq pointer to mesh quantities structures
2279 */
2280 /*----------------------------------------------------------------------------*/
2281
2282 void
cs_mesh_quantities_free_all(cs_mesh_quantities_t * mq)2283 cs_mesh_quantities_free_all(cs_mesh_quantities_t *mq)
2284 {
2285 CS_FREE_HD(mq->cell_cen);
2286 BFT_FREE(mq->cell_vol);
2287 BFT_FREE(mq->i_face_normal);
2288 CS_FREE_HD(mq->b_face_normal);
2289 BFT_FREE(mq->i_face_cog);
2290 CS_FREE_HD(mq->b_face_cog);
2291 BFT_FREE(mq->i_face_surf);
2292 CS_FREE_HD(mq->b_face_surf);
2293 BFT_FREE(mq->i_dist);
2294 CS_FREE_HD(mq->b_dist);
2295 CS_FREE_HD(mq->weight);
2296 BFT_FREE(mq->dijpf);
2297 CS_FREE_HD(mq->diipb);
2298 BFT_FREE(mq->dofij);
2299 BFT_FREE(mq->diipf);
2300 BFT_FREE(mq->djjpf);
2301 BFT_FREE(mq->corr_grad_lin_det);
2302 BFT_FREE(mq->corr_grad_lin);
2303 BFT_FREE(mq->b_sym_flag);
2304 BFT_FREE(mq->c_disable_flag);
2305 BFT_FREE(mq->bad_cell_flag);
2306 }
2307
2308 /*----------------------------------------------------------------------------*/
2309 /*!
2310 * \brief Compute mesh quantities needed for preprocessing.
2311 *
2312 * \param[in] m pointer to mesh structure
2313 * \param[in, out] mq pointer to mesh quantities structures.
2314 */
2315 /*----------------------------------------------------------------------------*/
2316
2317 void
cs_mesh_quantities_compute_preprocess(const cs_mesh_t * m,cs_mesh_quantities_t * mq)2318 cs_mesh_quantities_compute_preprocess(const cs_mesh_t *m,
2319 cs_mesh_quantities_t *mq)
2320 {
2321 cs_lnum_t n_i_faces = m->n_i_faces;
2322 cs_lnum_t n_b_faces = CS_MAX(m->n_b_faces, m->n_b_faces_all);
2323 cs_lnum_t n_cells_with_ghosts = m->n_cells_with_ghosts;
2324
2325 /* If this is not an update, allocate members of the structure */
2326
2327 if (mq->cell_cen == NULL)
2328 CS_MALLOC_HD(mq->cell_cen, n_cells_with_ghosts*3, cs_real_t, cs_alloc_mode);
2329
2330 if (mq->cell_vol == NULL)
2331 BFT_MALLOC(mq->cell_vol, n_cells_with_ghosts, cs_real_t);
2332
2333 if (mq->i_face_normal == NULL)
2334 BFT_MALLOC(mq->i_face_normal, n_i_faces*3, cs_real_t);
2335
2336 if (mq->b_face_normal == NULL)
2337 CS_MALLOC_HD(mq->b_face_normal, n_b_faces*3, cs_real_t, cs_alloc_mode);
2338
2339 if (mq->i_face_cog == NULL)
2340 BFT_MALLOC(mq->i_face_cog, n_i_faces*3, cs_real_t);
2341
2342 if (mq->b_face_cog == NULL)
2343 CS_MALLOC_HD(mq->b_face_cog, n_b_faces*3, cs_real_t, cs_alloc_mode);
2344
2345 if (mq->i_face_surf == NULL)
2346 BFT_MALLOC(mq->i_face_surf, n_i_faces, cs_real_t);
2347
2348 if (mq->b_face_surf == NULL)
2349 CS_MALLOC_HD(mq->b_face_surf, n_b_faces, cs_real_t, cs_alloc_mode);
2350
2351 /* Compute face centers of gravity, normals, and surfaces */
2352
2353 _compute_face_quantities(n_i_faces,
2354 (const cs_real_3_t *)m->vtx_coord,
2355 m->i_face_vtx_idx,
2356 m->i_face_vtx_lst,
2357 (cs_real_3_t *)mq->i_face_cog,
2358 (cs_real_3_t *)mq->i_face_normal);
2359
2360 _compute_face_surface(n_i_faces,
2361 mq->i_face_normal,
2362 mq->i_face_surf);
2363
2364 _compute_face_quantities(n_b_faces,
2365 (const cs_real_3_t *)m->vtx_coord,
2366 m->b_face_vtx_idx,
2367 m->b_face_vtx_lst,
2368 (cs_real_3_t *)mq->b_face_cog,
2369 (cs_real_3_t *)mq->b_face_normal);
2370
2371 _compute_face_surface(n_b_faces,
2372 mq->b_face_normal,
2373 mq->b_face_surf);
2374
2375 if (cs_glob_mesh_quantities_flag & CS_FACE_CENTER_REFINE) {
2376 _refine_warped_face_centers
2377 (n_i_faces,
2378 (const cs_real_3_t *)m->vtx_coord,
2379 m->i_face_vtx_idx,
2380 m->i_face_vtx_lst,
2381 (cs_real_3_t *)mq->i_face_cog,
2382 (const cs_real_3_t *)mq->i_face_normal);
2383
2384 _refine_warped_face_centers
2385 (n_b_faces,
2386 (const cs_real_3_t *)m->vtx_coord,
2387 m->b_face_vtx_idx,
2388 m->b_face_vtx_lst,
2389 (cs_real_3_t *)mq->b_face_cog,
2390 (const cs_real_3_t *)mq->b_face_normal);
2391 }
2392
2393 if (_ajust_face_cog_compat_v11_v52) {
2394 _adjust_face_cog_v11_v52
2395 (n_i_faces,
2396 (const cs_real_3_t *)m->vtx_coord,
2397 m->i_face_vtx_idx,
2398 m->i_face_vtx_lst,
2399 (cs_real_3_t *)mq->i_face_cog,
2400 (const cs_real_3_t *)mq->i_face_normal);
2401
2402 _adjust_face_cog_v11_v52
2403 (n_b_faces,
2404 (const cs_real_3_t *)m->vtx_coord,
2405 m->b_face_vtx_idx,
2406 m->b_face_vtx_lst,
2407 (cs_real_3_t *)mq->b_face_cog,
2408 (const cs_real_3_t *)mq->b_face_normal);
2409 }
2410
2411 /* Compute cell centers from face barycenters or vertices */
2412
2413 bool volume_computed = false;
2414
2415 switch (_cell_cen_algorithm) {
2416
2417 case 0:
2418 cs_mesh_quantities_cell_faces_cog(m,
2419 mq->i_face_normal,
2420 mq->i_face_cog,
2421 mq->b_face_normal,
2422 mq->b_face_cog,
2423 mq->cell_cen);
2424
2425 break;
2426 case 1:
2427 _compute_cell_quantities(m,
2428 (const cs_real_3_t *)mq->i_face_normal,
2429 (const cs_real_3_t *)mq->i_face_cog,
2430 (const cs_real_3_t *)mq->b_face_normal,
2431 (const cs_real_3_t *)mq->b_face_cog,
2432 (cs_real_3_t *)mq->cell_cen,
2433 mq->cell_vol);
2434 volume_computed = true;
2435 break;
2436
2437 default:
2438 assert(0);
2439
2440 }
2441
2442 if (cs_glob_mesh_quantities_flag & CS_CELL_CENTER_CORRECTION) {
2443 _recompute_cell_cen_face(m,
2444 (const cs_real_3_t *)(mq->i_face_normal),
2445 (const cs_real_3_t *)(mq->i_face_cog),
2446 (const cs_real_3_t *)(mq->b_face_normal),
2447 (const cs_real_3_t *)(mq->b_face_cog),
2448 (cs_real_3_t *)(mq->cell_cen));
2449 volume_computed = false; /* should not be different with plane faces,
2450 not sure with warped faces */
2451 }
2452
2453 /* Recompute face centers as the middle of two cell centers if possible */
2454 if (cs_glob_mesh_quantities_flag & CS_CELL_FACE_CENTER_CORRECTION) {
2455
2456 if (m->halo != NULL) {
2457 cs_halo_sync_var_strided(m->halo, CS_HALO_EXTENDED,
2458 mq->cell_cen, 3);
2459 if (m->n_init_perio > 0)
2460 cs_halo_perio_sync_coords(m->halo, CS_HALO_EXTENDED,
2461 mq->cell_cen);
2462 }
2463
2464 _correct_cell_face_center(m,
2465 n_cells_with_ghosts,
2466 n_i_faces,
2467 n_b_faces,
2468 (const cs_lnum_2_t *)(m->i_face_cells),
2469 m->b_face_cells,
2470 (cs_real_3_t *)(mq->cell_cen),
2471 (cs_real_3_t *)(mq->i_face_cog),
2472 (cs_real_3_t *)(mq->b_face_cog),
2473 (cs_real_3_t *)(mq->i_face_normal),
2474 (cs_real_3_t *)(mq->b_face_normal));
2475
2476 volume_computed = false;
2477
2478 }
2479
2480 /* Compute the volume of cells */
2481
2482 if (volume_computed == false)
2483 _compute_cell_volume(m,
2484 (const cs_real_3_t *)(mq->i_face_normal),
2485 (const cs_real_3_t *)(mq->i_face_cog),
2486 (const cs_real_3_t *)(mq->b_face_normal),
2487 (const cs_real_3_t *)(mq->b_face_cog),
2488 (const cs_real_3_t *)(mq->cell_cen),
2489 mq->cell_vol);
2490
2491 /* Correction of small or negative volumes
2492 (doesn't conserve the total volume) */
2493
2494 if (cs_glob_mesh_quantities_flag & CS_CELL_VOLUME_RATIO_CORRECTION)
2495 _cell_bad_volume_correction(m,
2496 mq->cell_vol);
2497
2498 /* Synchronize geometric quantities */
2499
2500 if (m->halo != NULL) {
2501
2502 cs_halo_sync_var_strided(m->halo, CS_HALO_EXTENDED,
2503 mq->cell_cen, 3);
2504 if (m->n_init_perio > 0)
2505 cs_halo_perio_sync_coords(m->halo, CS_HALO_EXTENDED,
2506 mq->cell_cen);
2507
2508 cs_halo_sync_var(m->halo, CS_HALO_EXTENDED, mq->cell_vol);
2509
2510 }
2511
2512 _cell_volume_reductions(m,
2513 mq->cell_vol,
2514 &(mq->min_vol),
2515 &(mq->max_vol),
2516 &(mq->tot_vol));
2517
2518 #if defined(HAVE_MPI)
2519 if (cs_glob_n_ranks > 1) {
2520
2521 cs_real_t _min_vol, _max_vol, _tot_vol;
2522
2523 MPI_Allreduce(&(mq->min_vol), &_min_vol, 1, CS_MPI_REAL,
2524 MPI_MIN, cs_glob_mpi_comm);
2525
2526 MPI_Allreduce(&(mq->max_vol), &_max_vol, 1, CS_MPI_REAL,
2527 MPI_MAX, cs_glob_mpi_comm);
2528
2529 MPI_Allreduce(&(mq->tot_vol), &_tot_vol, 1, CS_MPI_REAL,
2530 MPI_SUM, cs_glob_mpi_comm);
2531
2532 mq->min_vol = _min_vol;
2533 mq->max_vol = _max_vol;
2534 mq->tot_vol = _tot_vol;
2535
2536 }
2537 #endif
2538 }
2539
2540 /*----------------------------------------------------------------------------*/
2541 /*!
2542 * \brief Compute mesh quantities.
2543 *
2544 * \param[in] m pointer to mesh structure
2545 * \param[in, out] mq pointer to mesh quantities structures.
2546 */
2547 /*----------------------------------------------------------------------------*/
2548
2549 void
cs_mesh_quantities_compute(const cs_mesh_t * m,cs_mesh_quantities_t * mq)2550 cs_mesh_quantities_compute(const cs_mesh_t *m,
2551 cs_mesh_quantities_t *mq)
2552 {
2553 cs_lnum_t dim = m->dim;
2554 cs_lnum_t n_i_faces = m->n_i_faces;
2555 cs_lnum_t n_b_faces = m->n_b_faces;
2556 cs_lnum_t n_cells_with_ghosts = m->n_cells_with_ghosts;
2557
2558 /* Update the number of passes */
2559
2560 _n_computations++;
2561
2562 cs_mesh_quantities_compute_preprocess(m, mq);
2563
2564 /* Fluid surfaces and volumes: point to standard quantities and
2565 * may be modified afterwards */
2566 mq->i_f_face_normal = mq->i_face_normal;
2567 mq->b_f_face_normal = mq->b_face_normal;
2568 mq->i_f_face_surf = mq->i_face_surf;
2569 mq->b_f_face_surf = mq->b_face_surf;
2570
2571 mq->cell_f_vol = mq->cell_vol;
2572
2573 /* Porous models */
2574 if (mq->c_disable_flag == NULL) {
2575 if (mq->has_disable_flag == 1) {
2576 cs_lnum_t n_cells_ext = n_cells_with_ghosts;
2577 BFT_MALLOC(mq->c_disable_flag, n_cells_ext, int);
2578 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
2579 mq->c_disable_flag[cell_id] = 0;
2580 }
2581 else {
2582 BFT_MALLOC(mq->c_disable_flag, 1, int);
2583 mq->c_disable_flag[0] = 0;
2584 }
2585 }
2586
2587 mq->min_f_vol = mq->min_vol;
2588 mq->max_f_vol = mq->max_vol;
2589 mq->tot_f_vol = mq->tot_vol;
2590
2591 if (mq->i_dist == NULL)
2592 BFT_MALLOC(mq->i_dist, n_i_faces, cs_real_t);
2593
2594 if (mq->b_dist == NULL)
2595 CS_MALLOC_HD(mq->b_dist, n_b_faces, cs_real_t, cs_alloc_mode);
2596
2597 if (mq->weight == NULL)
2598 CS_MALLOC_HD(mq->weight, n_i_faces, cs_real_t, cs_alloc_mode);
2599
2600 if (mq->dijpf == NULL)
2601 BFT_MALLOC(mq->dijpf, n_i_faces*dim, cs_real_t);
2602
2603 if (mq->diipb == NULL)
2604 CS_MALLOC_HD(mq->diipb, n_b_faces*dim, cs_real_t, cs_alloc_mode);
2605
2606 if (mq->dofij == NULL)
2607 BFT_MALLOC(mq->dofij, n_i_faces*dim, cs_real_t);
2608
2609 if (mq->diipf == NULL)
2610 BFT_MALLOC(mq->diipf, n_i_faces*dim, cs_real_t);
2611
2612 if (mq->djjpf == NULL)
2613 BFT_MALLOC(mq->djjpf, n_i_faces*dim, cs_real_t);
2614
2615 if (mq->b_sym_flag == NULL) {
2616 BFT_MALLOC(mq->b_sym_flag, n_b_faces, int);
2617 for (cs_lnum_t i = 0; i < n_b_faces; i++)
2618 mq->b_sym_flag[i] = 1;
2619 }
2620
2621 /* Compute some distances relative to faces and associated weighting */
2622
2623 _compute_face_distances(m->n_i_faces,
2624 m->n_b_faces,
2625 (const cs_lnum_2_t *)(m->i_face_cells),
2626 m->b_face_cells,
2627 (const cs_real_3_t *)(mq->i_face_normal),
2628 (const cs_real_3_t *)(mq->b_face_normal),
2629 (const cs_real_3_t *)(mq->i_face_cog),
2630 (const cs_real_3_t *)(mq->b_face_cog),
2631 (const cs_real_3_t *)(mq->cell_cen),
2632 (const cs_real_t *)(mq->cell_vol),
2633 mq->i_dist,
2634 mq->b_dist,
2635 mq->weight);
2636
2637 /* Compute some vectors relative to faces to handle non-orthogonalities */
2638
2639 _compute_face_vectors(dim,
2640 m->n_i_faces,
2641 m->n_b_faces,
2642 (const cs_lnum_2_t *)(m->i_face_cells),
2643 m->b_face_cells,
2644 mq->i_face_normal,
2645 mq->b_face_normal,
2646 mq->i_face_cog,
2647 mq->b_face_cog,
2648 mq->i_face_surf,
2649 mq->cell_cen,
2650 mq->weight,
2651 mq->b_dist,
2652 mq->dijpf,
2653 mq->diipb,
2654 mq->dofij);
2655
2656 /* Compute additional vectors relative to faces to handle non-orthogonalities */
2657
2658 _compute_face_sup_vectors
2659 (m->n_cells,
2660 m->n_i_faces,
2661 (const cs_lnum_2_t *)(m->i_face_cells),
2662 (const cs_real_3_t *)(mq->i_face_normal),
2663 (const cs_real_3_t *)(mq->i_face_cog),
2664 (const cs_real_3_t *)(mq->cell_cen),
2665 mq->cell_vol,
2666 mq->i_dist,
2667 (cs_real_3_t *)(mq->diipf),
2668 (cs_real_3_t *)(mq->djjpf));
2669
2670 /* Build the geometrical matrix linear gradient correction */
2671 if (cs_glob_mesh_quantities_flag & CS_BAD_CELLS_WARPED_CORRECTION)
2672 _compute_corr_grad_lin(m, mq);
2673
2674 /* Print some information on the control volumes, and check min volume */
2675
2676 if (_n_computations == 1)
2677 bft_printf(_(" --- Information on the volumes\n"
2678 " Minimum control volume = %14.7e\n"
2679 " Maximum control volume = %14.7e\n"
2680 " Total volume for the domain = %14.7e\n"),
2681 mq->min_vol, mq->max_vol,
2682 mq->tot_vol);
2683 else {
2684 if (mq->min_vol <= 0.) {
2685 bft_printf(_(" --- Information on the volumes\n"
2686 " Minimum control volume = %14.7e\n"
2687 " Maximum control volume = %14.7e\n"
2688 " Total volume for the domain = %14.7e\n"),
2689 mq->min_vol, mq->max_vol,
2690 mq->tot_vol);
2691 bft_printf(_("\nAbort due to the detection of a negative control "
2692 "volume.\n"));
2693 }
2694 }
2695 }
2696
2697 /*----------------------------------------------------------------------------
2698 * Compute min, max, and total
2699 *
2700 * parameters:
2701 * mesh <-- pointer to a cs_mesh_t structure
2702 * mesh_quantities <-> pointer to a cs_mesh_quantities_t structure
2703 *----------------------------------------------------------------------------*/
2704
2705 void
cs_mesh_quantities_fluid_compute(const cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities)2706 cs_mesh_quantities_fluid_compute(const cs_mesh_t *mesh,
2707 cs_mesh_quantities_t *mesh_quantities)
2708 {
2709 CS_UNUSED(mesh);
2710 CS_UNUSED(mesh_quantities);
2711 }
2712
2713 /*----------------------------------------------------------------------------
2714 * Compute the total, min, and max fluid volumes of cells
2715 *
2716 * parameters:
2717 * mesh <-- pointer to mesh structure
2718 * mesh_quantities <-> pointer to a mesh quantities structure
2719 *----------------------------------------------------------------------------*/
2720
2721 void
cs_mesh_quantities_fluid_vol_reductions(const cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities)2722 cs_mesh_quantities_fluid_vol_reductions(const cs_mesh_t *mesh,
2723 cs_mesh_quantities_t *mesh_quantities)
2724 {
2725 _cell_volume_reductions(mesh,
2726 mesh_quantities->cell_f_vol,
2727 &(mesh_quantities->min_f_vol),
2728 &(mesh_quantities->max_f_vol),
2729 &(mesh_quantities->tot_f_vol));
2730
2731 #if defined(HAVE_MPI)
2732 if (cs_glob_n_ranks > 1) {
2733
2734 cs_real_t _min_f_vol, _max_f_vol, _tot_f_vol;
2735
2736 MPI_Allreduce(&(mesh_quantities->min_f_vol), &_min_f_vol, 1, CS_MPI_REAL,
2737 MPI_MIN, cs_glob_mpi_comm);
2738
2739 MPI_Allreduce(&(mesh_quantities->max_f_vol), &_max_f_vol, 1, CS_MPI_REAL,
2740 MPI_MAX, cs_glob_mpi_comm);
2741
2742 MPI_Allreduce(&(mesh_quantities->tot_f_vol), &_tot_f_vol, 1, CS_MPI_REAL,
2743 MPI_SUM, cs_glob_mpi_comm);
2744
2745 mesh_quantities->min_f_vol = _min_f_vol;
2746 mesh_quantities->max_f_vol = _max_f_vol;
2747 mesh_quantities->tot_f_vol = _tot_f_vol;
2748
2749 }
2750 #endif
2751 }
2752
2753 /*----------------------------------------------------------------------------
2754 * Compute fluid section mesh quantities at the initial step
2755 *
2756 * parameters:
2757 * mesh <-- pointer to a cs_mesh_t structure
2758 * mesh_quantities <-> pointer to a cs_mesh_quantities_t structure
2759 *----------------------------------------------------------------------------*/
2760
2761 void
cs_mesh_init_fluid_sections(const cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities)2762 cs_mesh_init_fluid_sections(const cs_mesh_t *mesh,
2763 cs_mesh_quantities_t *mesh_quantities)
2764 {
2765 cs_lnum_t n_i_faces = mesh->n_i_faces;
2766 cs_lnum_t n_b_faces = mesh->n_b_faces;
2767
2768 cs_real_3_t *restrict i_face_normal
2769 = (cs_real_3_t *restrict)mesh_quantities->i_face_normal;
2770 cs_real_3_t *restrict b_face_normal
2771 = (cs_real_3_t *restrict)mesh_quantities->b_face_normal;
2772 cs_real_3_t *restrict i_f_face_normal
2773 = (cs_real_3_t *restrict)mesh_quantities->i_f_face_normal;
2774 cs_real_3_t *restrict b_f_face_normal
2775 = (cs_real_3_t *restrict)mesh_quantities->b_f_face_normal;
2776
2777 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
2778 mesh_quantities->i_f_face_surf[face_id] =
2779 mesh_quantities->i_face_surf[face_id];
2780
2781 for (int i = 0; i < 3; i++)
2782 i_f_face_normal[face_id][i] = i_face_normal[face_id][i];
2783
2784 mesh_quantities->i_f_face_factor[face_id][0] = 1.;
2785 mesh_quantities->i_f_face_factor[face_id][1] = 1.;
2786 }
2787
2788 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
2789 mesh_quantities->b_f_face_surf[face_id]
2790 = mesh_quantities->b_face_surf[face_id];
2791
2792 for (int i = 0; i < 3; i++)
2793 b_f_face_normal[face_id][i] = b_face_normal[face_id][i];
2794
2795 mesh_quantities->b_f_face_factor[face_id] = 1.;
2796 }
2797 }
2798
2799 /*----------------------------------------------------------------------------
2800 * Compute mesh quantities -> vectors II' and JJ'
2801 *
2802 * parameters:
2803 * mesh <-- pointer to a cs_mesh_t structure
2804 * mesh_quantities <-> pointer to a cs_mesh_quantities_t structure
2805 *----------------------------------------------------------------------------*/
2806
2807 void
cs_mesh_quantities_sup_vectors(const cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities)2808 cs_mesh_quantities_sup_vectors(const cs_mesh_t *mesh,
2809 cs_mesh_quantities_t *mesh_quantities)
2810 {
2811 cs_lnum_t dim = mesh->dim;
2812 cs_lnum_t n_i_faces = mesh->n_i_faces;
2813
2814 if (mesh_quantities->diipf == NULL)
2815 BFT_MALLOC(mesh_quantities->diipf, n_i_faces*dim, cs_real_t);
2816
2817 if (mesh_quantities->djjpf == NULL)
2818 BFT_MALLOC(mesh_quantities->djjpf, n_i_faces*dim, cs_real_t);
2819
2820 _compute_face_sup_vectors
2821 (mesh->n_cells,
2822 mesh->n_i_faces,
2823 (const cs_lnum_2_t *)(mesh->i_face_cells),
2824 (const cs_real_3_t *)(mesh_quantities->i_face_normal),
2825 (const cs_real_3_t *)(mesh_quantities->i_face_cog),
2826 (const cs_real_3_t *)(mesh_quantities->cell_cen),
2827 mesh_quantities->cell_vol,
2828 mesh_quantities->i_dist,
2829 (cs_real_3_t *)(mesh_quantities->diipf),
2830 (cs_real_3_t *)(mesh_quantities->djjpf));
2831 }
2832
2833 /*----------------------------------------------------------------------------
2834 * Compute internal and border face normal.
2835 *
2836 * parameters:
2837 * mesh <-- pointer to a cs_mesh_t structure
2838 * p_i_face_normal <-> pointer to the internal face normal array
2839 * p_b_face_normal <-> pointer to the border face normal array
2840 *----------------------------------------------------------------------------*/
2841
2842 void
cs_mesh_quantities_face_normal(const cs_mesh_t * mesh,cs_real_t * p_i_face_normal[],cs_real_t * p_b_face_normal[])2843 cs_mesh_quantities_face_normal(const cs_mesh_t *mesh,
2844 cs_real_t *p_i_face_normal[],
2845 cs_real_t *p_b_face_normal[])
2846 {
2847 cs_real_t *i_face_normal = NULL, *b_face_normal = NULL;
2848
2849 const cs_lnum_t n_b_faces = mesh->n_b_faces;
2850 const cs_lnum_t n_i_faces = mesh->n_i_faces;
2851
2852 /* Internal face treatment */
2853
2854 BFT_MALLOC(i_face_normal, n_i_faces * 3, cs_real_t);
2855
2856 _compute_face_normal(mesh->n_i_faces,
2857 (const cs_real_3_t *)mesh->vtx_coord,
2858 mesh->i_face_vtx_idx,
2859 mesh->i_face_vtx_lst,
2860 (cs_real_3_t *)i_face_normal);
2861
2862 *p_i_face_normal = i_face_normal;
2863
2864 /* Boundary face treatment */
2865
2866 BFT_MALLOC(b_face_normal, n_b_faces * 3, cs_real_t);
2867
2868 _compute_face_normal(mesh->n_b_faces,
2869 (const cs_real_3_t *)mesh->vtx_coord,
2870 mesh->b_face_vtx_idx,
2871 mesh->b_face_vtx_lst,
2872 (cs_real_3_t *)b_face_normal);
2873
2874 *p_b_face_normal = b_face_normal;
2875 }
2876
2877 /*----------------------------------------------------------------------------
2878 * Compute interior face centers and normals.
2879 *
2880 * The corresponding arrays are allocated by this function, and it is the
2881 * caller's responsibility to free them when they are no longer needed.
2882 *
2883 * parameters:
2884 * mesh <-- pointer to a cs_mesh_t structure
2885 * p_i_face_cog <-> pointer to the interior face center array
2886 * p_i_face_normal <-> pointer to the interior face normal array
2887 *----------------------------------------------------------------------------*/
2888
2889 void
cs_mesh_quantities_i_faces(const cs_mesh_t * mesh,cs_real_t * p_i_face_cog[],cs_real_t * p_i_face_normal[])2890 cs_mesh_quantities_i_faces(const cs_mesh_t *mesh,
2891 cs_real_t *p_i_face_cog[],
2892 cs_real_t *p_i_face_normal[])
2893 {
2894 cs_real_t *i_face_cog = NULL, *i_face_normal = NULL;
2895
2896 BFT_MALLOC(i_face_cog, mesh->n_i_faces * mesh->dim, cs_real_t);
2897 BFT_MALLOC(i_face_normal, mesh->n_i_faces * mesh->dim, cs_real_t);
2898
2899 _compute_face_quantities(mesh->n_i_faces,
2900 (const cs_real_3_t *)mesh->vtx_coord,
2901 mesh->i_face_vtx_idx,
2902 mesh->i_face_vtx_lst,
2903 (cs_real_3_t *)i_face_cog,
2904 (cs_real_3_t *)i_face_normal);
2905
2906 *p_i_face_cog = i_face_cog;
2907 *p_i_face_normal = i_face_normal;
2908 }
2909
2910 /*----------------------------------------------------------------------------
2911 * Compute border face centers and normals.
2912 *
2913 * The corresponding arrays are allocated by this function, and it is the
2914 * caller's responsibility to free them when they are no longer needed.
2915 *
2916 * parameters:
2917 * mesh <-- pointer to a cs_mesh_t structure
2918 * p_b_face_cog <-> pointer to the border face center array
2919 * p_b_face_normal <-> pointer to the border face normal array
2920 *----------------------------------------------------------------------------*/
2921
2922 void
cs_mesh_quantities_b_faces(const cs_mesh_t * mesh,cs_real_t * p_b_face_cog[],cs_real_t * p_b_face_normal[])2923 cs_mesh_quantities_b_faces(const cs_mesh_t *mesh,
2924 cs_real_t *p_b_face_cog[],
2925 cs_real_t *p_b_face_normal[])
2926 {
2927 cs_real_t *b_face_cog = NULL, *b_face_normal = NULL;
2928
2929 BFT_MALLOC(b_face_cog, mesh->n_b_faces * mesh->dim, cs_real_t);
2930 BFT_MALLOC(b_face_normal, mesh->n_b_faces * mesh->dim, cs_real_t);
2931
2932 _compute_face_quantities(mesh->n_b_faces,
2933 (const cs_real_3_t *)mesh->vtx_coord,
2934 mesh->b_face_vtx_idx,
2935 mesh->b_face_vtx_lst,
2936 (cs_real_3_t *)b_face_cog,
2937 (cs_real_3_t *)b_face_normal);
2938
2939 *p_b_face_cog = b_face_cog;
2940 *p_b_face_normal = b_face_normal;
2941 }
2942
2943 /*----------------------------------------------------------------------------*/
2944 /*!
2945 * \brief Compute approximate cells centers as the mean of the given face
2946 * centers weighted by the associated surfaces.
2947 *
2948 * n-1
2949 * Sum Surf(Fi) G(Fi)
2950 * i=0
2951 * G(C) = -----------------------
2952 * n-1
2953 * Sum Surf(Fi)
2954 * i=0
2955 *
2956 * \param[in] mesh pointer to mesh structure
2957 * \param[in] i_face_norm surface normal of internal faces
2958 * \param[in] i_face_cog center of gravity of internal faces
2959 * \param[in] b_face_norm surface normal of border faces
2960 * \param[in] b_face_cog center of gravity of border faces
2961 * \param[out] cell_cen cell centers
2962 */
2963 /*----------------------------------------------------------------------------*/
2964
2965 void
cs_mesh_quantities_cell_faces_cog(const cs_mesh_t * mesh,const cs_real_t i_face_norm[],const cs_real_t i_face_cog[],const cs_real_t b_face_norm[],const cs_real_t b_face_cog[],cs_real_t cell_cen[])2966 cs_mesh_quantities_cell_faces_cog(const cs_mesh_t *mesh,
2967 const cs_real_t i_face_norm[],
2968 const cs_real_t i_face_cog[],
2969 const cs_real_t b_face_norm[],
2970 const cs_real_t b_face_cog[],
2971 cs_real_t cell_cen[])
2972 {
2973 cs_real_t *cell_area = NULL;
2974
2975 /* Mesh connectivity */
2976
2977 const cs_lnum_t n_i_faces = mesh->n_i_faces;
2978 const cs_lnum_t n_b_faces = mesh->n_b_faces;
2979 const cs_lnum_t n_cells = mesh->n_cells;
2980 const cs_lnum_t n_cells_with_ghosts = mesh->n_cells_with_ghosts;
2981 const cs_lnum_2_t *i_face_cells
2982 = (const cs_lnum_2_t *)(mesh->i_face_cells);
2983 const cs_lnum_t *b_face_cells = mesh->b_face_cells;
2984
2985 /* Return if ther is not enough data (Solcom case except rediative module
2986 or Pre-processor 1.2.d without option "-n") */
2987
2988 if (mesh->i_face_vtx_lst == NULL && mesh->b_face_vtx_lst == NULL)
2989 return;
2990
2991 /* Checking */
2992
2993 assert(cell_cen != NULL);
2994
2995 /* Initialization */
2996
2997 BFT_MALLOC(cell_area, n_cells_with_ghosts, cs_real_t);
2998
2999 for (cs_lnum_t j = 0; j < n_cells_with_ghosts; j++) {
3000
3001 cell_area[j] = 0.;
3002
3003 for (cs_lnum_t i = 0; i < 3; i++)
3004 cell_cen[3*j + i] = 0.;
3005
3006 }
3007
3008 /* Loop on interior faces
3009 ---------------------- */
3010
3011 for (cs_lnum_t f_id = 0; f_id < n_i_faces; f_id++) {
3012
3013 /* For each cell sharing the internal face, we update
3014 * cell_cen and cell_area */
3015
3016 cs_lnum_t c_id1 = i_face_cells[f_id][0];
3017 cs_lnum_t c_id2 = i_face_cells[f_id][1];
3018
3019 /* Computation of the area of the face */
3020
3021 cs_real_t area = cs_math_3_norm(i_face_norm + 3*f_id);
3022
3023 if (c_id1 > -1) {
3024 cell_area[c_id1] += area;
3025 for (cs_lnum_t i = 0; i < 3; i++)
3026 cell_cen[3*c_id1 + i] += i_face_cog[3*f_id + i]*area;
3027 }
3028 if (c_id2 > -1) {
3029 cell_area[c_id2] += area;
3030 for (cs_lnum_t i = 0; i < 3; i++)
3031 cell_cen[3*c_id2 + i] += i_face_cog[3*f_id + i]*area;
3032 }
3033
3034 } /* End of loop on interior faces */
3035
3036 /* Loop on boundary faces
3037 --------------------- */
3038
3039 for (cs_lnum_t f_id = 0; f_id < n_b_faces; f_id++) {
3040
3041 /* For each cell sharing a border face, we update the numerator
3042 * of cell_cen and cell_area */
3043
3044 cs_lnum_t c_id1 = b_face_cells[f_id];
3045
3046 /* Computation of the area of the face
3047 (note that c_id1 == -1 may happen for isolated faces,
3048 which are cleaned afterwards) */
3049
3050 if (c_id1 > -1) {
3051
3052 cs_real_t area = cs_math_3_norm(b_face_norm + 3*f_id);
3053
3054 cell_area[c_id1] += area;
3055
3056 /* Computation of the numerator */
3057
3058 for (cs_lnum_t i = 0; i < 3; i++)
3059 cell_cen[3*c_id1 + i] += b_face_cog[3*f_id + i]*area;
3060
3061 }
3062
3063 } /* End of loop on boundary faces */
3064
3065 /* Loop on cells to finalize the computation of center of gravity
3066 -------------------------------------------------------------- */
3067
3068 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
3069
3070 for (cs_lnum_t i = 0; i < 3; i++)
3071 cell_cen[c_id*3 + i] /= cell_area[c_id];
3072
3073 }
3074
3075 /* Free memory */
3076
3077 BFT_FREE(cell_area);
3078 }
3079
3080 /*----------------------------------------------------------------------------
3081 * Compute cell volumes.
3082 *
3083 * The corresponding array is allocated by this function, and it is the
3084 * caller's responsability to free it when they are no longer needed.
3085 *
3086 * parameters:
3087 * mesh <-- pointer to a cs_mesh_t structure
3088 *
3089 * return:
3090 * pointer to newly allocated cell volumes array
3091 *----------------------------------------------------------------------------*/
3092
3093 cs_real_t *
cs_mesh_quantities_cell_volume(const cs_mesh_t * mesh)3094 cs_mesh_quantities_cell_volume(const cs_mesh_t *mesh)
3095 {
3096 cs_real_t *cell_vol;
3097 BFT_MALLOC(cell_vol, mesh->n_cells_with_ghosts, cs_real_t);
3098
3099 cs_real_3_t *cell_cen;
3100 BFT_MALLOC(cell_cen, mesh->n_cells_with_ghosts, cs_real_3_t);
3101
3102 cs_real_t *i_face_cog = NULL, *i_face_normal = NULL;
3103 cs_real_t *b_face_cog = NULL, *b_face_normal = NULL;
3104
3105 cs_mesh_quantities_i_faces(mesh, &i_face_cog, &i_face_normal);
3106 cs_mesh_quantities_b_faces(mesh, &b_face_cog, &b_face_normal);
3107
3108 _compute_cell_quantities(mesh,
3109 (const cs_real_3_t *)i_face_normal,
3110 (const cs_real_3_t *)i_face_cog,
3111 (const cs_real_3_t *)b_face_normal,
3112 (const cs_real_3_t *)b_face_cog,
3113 (cs_real_3_t *)cell_cen,
3114 cell_vol);
3115
3116 BFT_FREE(cell_cen);
3117 BFT_FREE(b_face_normal);
3118 BFT_FREE(b_face_cog);
3119 BFT_FREE(i_face_normal);
3120 BFT_FREE(i_face_cog);
3121
3122 return cell_vol;
3123 }
3124
3125 /*----------------------------------------------------------------------------
3126 * Check that no negative volumes are present, and exit on error otherwise.
3127 *
3128 * parameters:
3129 * mesh <-- pointer to mesh structure
3130 * mesh_quantities <-- pointer to mesh quantities structure
3131 * allow_error <-- 1 if errors are allowed, 0 otherwise
3132 *----------------------------------------------------------------------------*/
3133
3134 void
cs_mesh_quantities_check_vol(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities,int allow_error)3135 cs_mesh_quantities_check_vol(const cs_mesh_t *mesh,
3136 const cs_mesh_quantities_t *mesh_quantities,
3137 int allow_error)
3138 {
3139 cs_lnum_t cell_id;
3140
3141 cs_gnum_t error_count = 0;
3142
3143 for (cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
3144 if (mesh_quantities->cell_vol[cell_id] < 0.0)
3145 error_count += 1;
3146 }
3147
3148 #if defined(HAVE_MPI)
3149 if (cs_glob_n_ranks > 1) {
3150 cs_gnum_t tot_error_count = 0;
3151 MPI_Allreduce(&error_count, &tot_error_count, 1, CS_MPI_GNUM, MPI_SUM,
3152 cs_glob_mpi_comm);
3153 error_count = tot_error_count;
3154 }
3155 #endif
3156
3157 /* Exit with error */
3158
3159 if (error_count > 0) {
3160 const char fmt[]
3161 = N_(" %llu cells have a Negative volume.\n"
3162 " Run mesh quality check for post-processing output.\n"
3163 " In case of mesh joining, this may be due to overly "
3164 " agressive joining parameters.");
3165
3166 if (allow_error) {
3167 cs_base_warn(__FILE__, __LINE__);
3168 bft_printf(_(fmt), (unsigned long long)error_count);
3169 bft_printf("\n\n");
3170 }
3171 else
3172 bft_error(__FILE__, __LINE__, 0,
3173 _(fmt), (unsigned long long)error_count);
3174 }
3175 }
3176
3177 /*----------------------------------------------------------------------------*/
3178 /*!
3179 * \brief Compute the bounding box for cells.
3180 *
3181 * The corresponding array is allocated by this function, and it is the
3182 * caller's responsability to free it when they are no longer needed.
3183 *
3184 * \param[in] m pointer to mesh structure
3185 * \param[in] tolerance addition to local extents of each element:
3186 * extent = base_extent * (1 + tolerance)
3187 *
3188 * \return pointer to newly allocated cell volumes array
3189 */
3190 /*----------------------------------------------------------------------------*/
3191
3192 cs_real_6_t *
cs_mesh_quantities_cell_extents(const cs_mesh_t * m,cs_real_t tolerance)3193 cs_mesh_quantities_cell_extents(const cs_mesh_t *m,
3194 cs_real_t tolerance)
3195 {
3196 cs_real_6_t *bbox;
3197
3198 BFT_MALLOC(bbox, m->n_cells_with_ghosts, cs_real_6_t);
3199
3200 for (cs_lnum_t i = 0; i < m->n_cells_with_ghosts; i++) {
3201 bbox[i][0] = HUGE_VAL;
3202 bbox[i][1] = HUGE_VAL;
3203 bbox[i][2] = HUGE_VAL;
3204 bbox[i][3] = -HUGE_VAL;
3205 bbox[i][4] = -HUGE_VAL;
3206 bbox[i][5] = -HUGE_VAL;
3207 }
3208
3209 const cs_lnum_t n_i_faces = m->n_i_faces;
3210
3211 for (cs_lnum_t i = 0; i < n_i_faces; i++) {
3212 cs_lnum_t c_id_0 = m->i_face_cells[i][0];
3213 cs_lnum_t c_id_1 = m->i_face_cells[i][1];
3214 cs_lnum_t s_id = m->i_face_vtx_idx[i];
3215 cs_lnum_t e_id = m->i_face_vtx_idx[i+1];
3216 for (cs_lnum_t j = s_id; j < e_id; j++) {
3217 cs_lnum_t vtx_id = m->i_face_vtx_lst[j];
3218 const cs_real_t *coo = m->vtx_coord + vtx_id*3;
3219 if (c_id_0 > -1) {
3220 for (cs_lnum_t k = 0; k < 3; k++) {
3221 bbox[c_id_0][k] = CS_MIN(bbox[c_id_0][k], coo[k]);
3222 bbox[c_id_0][k+3] = CS_MAX(bbox[c_id_0][k+3], coo[k]);
3223 }
3224 }
3225 if (c_id_1 > -1) {
3226 for (cs_lnum_t k = 0; k < 3; k++) {
3227 bbox[c_id_1][k] = CS_MIN(bbox[c_id_1][k], coo[k]);
3228 bbox[c_id_1][k+3] = CS_MAX(bbox[c_id_1][k+3], coo[k]);
3229 }
3230 }
3231 }
3232 }
3233
3234 const cs_lnum_t n_b_faces = m->n_b_faces;
3235
3236 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
3237 cs_lnum_t c_id = m->b_face_cells[i];
3238 cs_lnum_t s_id = m->b_face_vtx_idx[i];
3239 cs_lnum_t e_id = m->b_face_vtx_idx[i+1];
3240 for (cs_lnum_t j = s_id; j < e_id; j++) {
3241 cs_lnum_t vtx_id = m->b_face_vtx_lst[j];
3242 const cs_real_t *coo = m->vtx_coord + vtx_id*3;
3243 if (c_id > -1) {
3244 for (cs_lnum_t k = 0; k < 3; k++) {
3245 bbox[c_id][k] = CS_MIN(bbox[c_id][k], coo[k]);
3246 bbox[c_id][k+3] = CS_MAX(bbox[c_id][k+3], coo[k]);
3247 }
3248 }
3249 }
3250 }
3251
3252 {
3253 const cs_lnum_t n_cells = m->n_cells;
3254
3255 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
3256
3257 cs_real_t delta[3];
3258
3259 for (cs_lnum_t i = 0; i < 3; i++)
3260 delta[i] = (bbox[c_id][3+i] - bbox[c_id][i]) * tolerance;
3261
3262 for (cs_lnum_t i = 0; i < 3; i++) {
3263 bbox[c_id][i] = bbox[c_id][i] - delta[i];
3264 bbox[c_id][3+i] = bbox[c_id][3+i] + delta[i];
3265 }
3266
3267 }
3268 }
3269
3270 return bbox;
3271 }
3272
3273 /*----------------------------------------------------------------------------
3274 * Return the number of times mesh quantities have been computed.
3275 *
3276 * returns:
3277 * number of times mesh quantities have been computed
3278 *----------------------------------------------------------------------------*/
3279
3280 int
cs_mesh_quantities_compute_count(void)3281 cs_mesh_quantities_compute_count(void)
3282 {
3283 return _n_computations;
3284 }
3285
3286 /*----------------------------------------------------------------------------*/
3287 /*!
3288 * \brief Determine local boundary thickness around each vertex.
3289 *
3290 * \param[in] m pointer to mesh structure
3291 * \param[in] mq pointer to mesh quantities structures.
3292 * \param[in] n_passes number of smoothing passes
3293 * \param[out] b_thickness thickness for each mesh vertex
3294 * (0 at non-boundary vertices)
3295 */
3296 /*----------------------------------------------------------------------------*/
3297
3298 void
cs_mesh_quantities_b_thickness_v(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,int n_passes,cs_real_t b_thickness[])3299 cs_mesh_quantities_b_thickness_v(const cs_mesh_t *m,
3300 const cs_mesh_quantities_t *mq,
3301 int n_passes,
3302 cs_real_t b_thickness[])
3303 {
3304 cs_real_t *v_sum = NULL;
3305 cs_real_t *f_b_thickness = NULL;
3306
3307 BFT_MALLOC(v_sum, m->n_vertices*2, cs_real_t);
3308
3309 BFT_MALLOC(f_b_thickness, m->n_b_faces*2, cs_real_t);
3310 _b_thickness(m, mq, f_b_thickness);
3311
3312 if (n_passes < 1)
3313 n_passes = 1;
3314
3315 for (int i = 0; i < n_passes; i++) {
3316
3317 for (cs_lnum_t j = 0; j < m->n_vertices*2; j++)
3318 v_sum[j] = 0.;
3319
3320 for (cs_lnum_t f_id = 0; f_id < m->n_b_faces; f_id++) {
3321 cs_lnum_t s_id = m->b_face_vtx_idx[f_id];
3322 cs_lnum_t e_id = m->b_face_vtx_idx[f_id+1];
3323 const cs_real_t f_s = mq->b_face_surf[f_id];
3324 for (cs_lnum_t k = s_id; k < e_id; k++) {
3325 cs_lnum_t v_id = m->b_face_vtx_lst[k];
3326 v_sum[v_id*2] += f_s * f_b_thickness[f_id];
3327 v_sum[v_id*2+1] += f_s;
3328 }
3329 }
3330
3331 if (m->vtx_interfaces != NULL)
3332 cs_interface_set_sum(m->vtx_interfaces,
3333 m->n_vertices,
3334 2,
3335 true,
3336 CS_REAL_TYPE,
3337 v_sum);
3338
3339 /* Prepare for next smoothing */
3340
3341 if (i < n_passes-1) {
3342
3343 for (cs_lnum_t j = 0; j < m->n_b_faces*2; j++)
3344 f_b_thickness[j] = 0.;
3345
3346 for (cs_lnum_t f_id = 0; f_id < m->n_b_faces; f_id++) {
3347 cs_lnum_t s_id = m->b_face_vtx_idx[f_id];
3348 cs_lnum_t e_id = m->b_face_vtx_idx[f_id+1];
3349 for (cs_lnum_t k = s_id; k < e_id; k++) {
3350 cs_lnum_t v_id = m->b_face_vtx_lst[k];
3351 f_b_thickness[f_id] += v_sum[v_id*2];
3352 f_b_thickness[f_id + m->n_b_faces] += v_sum[v_id*2 + 1];
3353 }
3354 }
3355
3356 for (cs_lnum_t f_id = 0; f_id < m->n_b_faces; f_id++) {
3357 if (f_b_thickness[f_id + m->n_b_faces] > 0)
3358 f_b_thickness[f_id] /= f_b_thickness[f_id + m->n_b_faces];
3359 }
3360
3361 }
3362
3363 }
3364
3365 BFT_FREE(f_b_thickness);
3366
3367 for (cs_lnum_t j = 0; j < m->n_vertices; j++) {
3368 if (v_sum[j*2+1] > 0)
3369 b_thickness[j] = v_sum[j*2] / v_sum[j*2+1];
3370 else
3371 b_thickness[j] = 0;
3372 }
3373
3374 BFT_FREE(v_sum);
3375 }
3376
3377 /*----------------------------------------------------------------------------*/
3378 /*!
3379 * \brief Determine local boundary thickness around each boundary face.
3380 *
3381 * \param[in] m pointer to mesh structure
3382 * \param[in] mq pointer to mesh quantities structures.
3383 * \param[in] n_passes number of optional smoothing passes
3384 * \param[out] b_thickness thickness for each mesh boundary face
3385 */
3386 /*----------------------------------------------------------------------------*/
3387
3388 void
cs_mesh_quantities_b_thickness_f(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,int n_passes,cs_real_t b_thickness[])3389 cs_mesh_quantities_b_thickness_f(const cs_mesh_t *m,
3390 const cs_mesh_quantities_t *mq,
3391 int n_passes,
3392 cs_real_t b_thickness[])
3393 {
3394 if (n_passes < 1)
3395 _b_thickness(m, mq, b_thickness);
3396
3397 else {
3398
3399 cs_real_t *v_b_thickness = NULL;
3400
3401 BFT_MALLOC(v_b_thickness, m->n_vertices, cs_real_t);
3402
3403 cs_mesh_quantities_b_thickness_v(m,
3404 mq,
3405 n_passes,
3406 v_b_thickness);
3407
3408 for (cs_lnum_t f_id = 0; f_id < m->n_b_faces; f_id++) {
3409 b_thickness[f_id] = 0;
3410 cs_lnum_t s_id = m->b_face_vtx_idx[f_id];
3411 cs_lnum_t e_id = m->b_face_vtx_idx[f_id+1];
3412 for (cs_lnum_t k = s_id; k < e_id; k++) {
3413 cs_lnum_t v_id = m->b_face_vtx_lst[k];
3414 b_thickness[f_id] += v_b_thickness[v_id];
3415 }
3416 b_thickness[f_id] /= (e_id - s_id);
3417 }
3418
3419 BFT_FREE(v_b_thickness);
3420
3421 }
3422 }
3423
3424 /*----------------------------------------------------------------------------*/
3425 /*!
3426 * \brief Log mesh quantities options to setup file.
3427 */
3428 /*----------------------------------------------------------------------------*/
3429
3430 void
cs_mesh_quantities_log_setup(void)3431 cs_mesh_quantities_log_setup(void)
3432 {
3433 if (cs_glob_mesh_quantities_flag != 0 || _cell_cen_algorithm != 1)
3434 cs_log_printf(CS_LOG_SETUP, _("\n"
3435 "Mesh quantity computation options\n"
3436 "---------------------------------\n\n"));
3437
3438 const char *cen_type_name[] = {N_("weighted center of face centers"),
3439 N_("center of mass")};
3440 cs_log_printf(CS_LOG_SETUP,
3441 _(" Cell centers: %s\n"),
3442 _(cen_type_name[_cell_cen_algorithm]));
3443
3444 if (cs_glob_mesh_quantities_flag != 0) {
3445
3446 const char *correction_name[] = {"CS_BAD_CELLS_WARPED_CORRECTION",
3447 "CS_BAD_CELLS_REGULARISATION",
3448 "CS_CELL_FACE_CENTER_CORRECTION",
3449 "CS_CELL_CENTER_CORRECTION",
3450 "CS_FACE_DISTANCE_CLIP",
3451 "CS_FACE_RECONSTRUCTION_CLIP",
3452 "CS_CELL_VOLUME_RATIO_CORRECTION",
3453 "CS_FACE_CENTER_REFINE"};
3454
3455 cs_log_printf(CS_LOG_SETUP,
3456 ("\n"
3457 " Mesh quantity corrections:\n"));
3458
3459 for (int i = 0; i < 8; i++) {
3460 if (cs_glob_mesh_quantities_flag & (1 << i))
3461 cs_log_printf(CS_LOG_SETUP, " %s\n", correction_name[i]);
3462 }
3463
3464 }
3465 }
3466
3467 /*----------------------------------------------------------------------------
3468 * Dump a cs_mesh_quantities_t structure
3469 *
3470 * parameters:
3471 * mesh <-- pointer to a cs_mesh_t structure
3472 * mesh_quantities <-- pointer to a cs_mesh_quantities_t structure
3473 *----------------------------------------------------------------------------*/
3474
3475 void
cs_mesh_quantities_dump(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mesh_quantities)3476 cs_mesh_quantities_dump(const cs_mesh_t *mesh,
3477 const cs_mesh_quantities_t *mesh_quantities)
3478 {
3479 cs_lnum_t i;
3480
3481 const cs_lnum_t n_cells = mesh->n_cells_with_ghosts;
3482 const cs_lnum_t n_i_faces = mesh->n_i_faces;
3483 const cs_lnum_t n_b_faces = mesh->n_b_faces;
3484
3485 const cs_real_t *cell_cen = mesh_quantities->cell_cen;
3486 const cs_real_t *cell_vol = mesh_quantities->cell_vol;
3487 const cs_real_t *i_fac_norm = mesh_quantities->i_face_normal;
3488 const cs_real_t *b_fac_norm = mesh_quantities->b_face_normal;
3489 const cs_real_t *i_fac_cog = mesh_quantities->i_face_cog;
3490 const cs_real_t *b_fac_cog = mesh_quantities->b_face_cog;
3491 const cs_real_t *i_fac_surf = mesh_quantities->i_face_surf;
3492 const cs_real_t *b_fac_surf = mesh_quantities->b_face_surf;
3493
3494 bft_printf("\n\nDUMP OF A MESH QUANTITIES STRUCTURE: %p\n\n",
3495 (const void *)mesh_quantities);
3496
3497 if (mesh_quantities == NULL)
3498 return;
3499
3500 /* Cell data */
3501
3502 bft_printf("\n\n"
3503 " ---------------"
3504 " Cell quantities"
3505 " ---------------\n\n");
3506
3507 bft_printf("Cell center coordinates:\n");
3508 for (i = 0; i < n_cells; i++)
3509 bft_printf(" < %ld > %.3f %.3f %.3f\n", (long)i+1,
3510 cell_cen[3*i], cell_cen[3*i+1], cell_cen[3*i+2]);
3511
3512 bft_printf("\nCell volume:\n");
3513 for (i = 0; i < n_cells; i++)
3514 bft_printf(" < %ld > %.3f\n", (long)i+1, cell_vol[i]);
3515
3516 /* Internal faces data */
3517
3518 bft_printf("\n\n"
3519 " ------------------------"
3520 " Interior face quantities"
3521 " ------------------------\n\n");
3522
3523 bft_printf("\nInterior face normals\n");
3524 for (i = 0; i < n_i_faces; i++)
3525 bft_printf(" < %ld > %.3f %.3f %.3f\n", (long)i+1,
3526 i_fac_norm[3*i], i_fac_norm[3*i+1], i_fac_norm[3*i+2]);
3527
3528 bft_printf("\nInterior face centers\n");
3529 for (i = 0; i < n_i_faces; i++)
3530 bft_printf(" < %ld > %.3f %.3f %.3f\n", (long)i+1,
3531 i_fac_cog[3*i], i_fac_cog[3*i+1], i_fac_cog[3*i+2]);
3532
3533 bft_printf("\nInterior face surfaces\n");
3534 for (i = 0; i < n_i_faces; i++)
3535 bft_printf(" < %ld > %.3f\n", (long)i+1, i_fac_surf[i]);
3536
3537 /* Border faces data */
3538
3539 bft_printf("\n\n"
3540 " ------------------------"
3541 " Boundary face quantities"
3542 " ------------------------\n\n");
3543
3544 bft_printf("\nBoundary face normals\n");
3545 for (i = 0; i < n_b_faces; i++)
3546 bft_printf(" < %ld > %.3f %.3f %.3f\n", (long)i+1,
3547 b_fac_norm[3*i], b_fac_norm[3*i+1], b_fac_norm[3*i+2]);
3548
3549 bft_printf("\nBoundary faces centers\n");
3550 for (i = 0; i < n_b_faces; i++)
3551 bft_printf(" < %ld > %.3f %.3f %.3f\n", (long)i+1,
3552 b_fac_cog[3*i], b_fac_cog[3*i+1], b_fac_cog[3*i+2]);
3553
3554 bft_printf("\nBoundary face surfaces\n");
3555 for (i = 0; i < n_b_faces; i++)
3556 bft_printf(" < %ld > %.3f\n", (long)i+1, b_fac_surf[i]);
3557
3558 bft_printf("\n\nEND OF DUMP OF MESH QUANTITIES STRUCTURE\n\n");
3559 bft_printf_flush();
3560 }
3561
3562 /*----------------------------------------------------------------------------*/
3563
3564 #if 0 /* Test if face orientation is OK */
3565
3566 cs_lnum_t i, fac_id, cell_id;
3567 cs_real_t cogfac[3];
3568 cs_real_t cogcel[3];
3569 cs_real_t normal[3];
3570 cs_real_t pscal;
3571
3572 for (fac_id = 0 ; fac_id < mesh->n_b_faces ; fac_id++) {
3573
3574 cell_id = mesh->b_face_cells[fac_id];
3575 pscal = 0;
3576
3577 for (i = 0 ; i < 3 ; i++) {
3578 cogcel[i] = cs_glob_mesh_quantities->cell_cen[cell_id*3 + i];
3579 cogfac[i] = cs_glob_mesh_quantities->b_face_cog[fac_id*3 + i];
3580 normal[i] = cs_glob_mesh_quantities->b_face_normal[fac_id*3 + i];
3581 pscal += normal[i] * (cogfac[i] - cogcel[i]);
3582 }
3583
3584 if (pscal < 0.0)
3585 printf("num_fac_brd = %d, num_cel = %d, pscal = %f\n",
3586 fac_id + 1, cell_id + 1, pscal);
3587 }
3588
3589 #endif
3590
3591 /*----------------------------------------------------------------------------*/
3592
3593 END_C_DECLS
3594