1 /* ===========================================================================
2 * Functions to handle low-level functions related to CDO local quantities:
3 * - local matrices (stored in dense format),
4 * - local mesh structure related to a cell or to a couple cell/face
5 *============================================================================*/
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
12 This program is free software; you can redistribute it and/or modify it under
13 the terms of the GNU General Public License as published by the Free Software
14 Foundation; either version 2 of the License, or (at your option) any later
15 version.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 details.
21
22 You should have received a copy of the GNU General Public License along with
23 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24 Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29 #include "cs_defs.h"
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <float.h>
37 #include <limits.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_mem.h>
44 #include <bft_printf.h>
45
46 #include "cs_math.h"
47 #include "cs_param_types.h"
48
49 /*----------------------------------------------------------------------------
50 * Header for the current file
51 *----------------------------------------------------------------------------*/
52
53 #include "cs_cdo_local.h"
54
55 /*----------------------------------------------------------------------------*/
56
57 BEGIN_C_DECLS
58
59 /*!
60 \file cs_cdo_local.c
61
62 \brief Functions to handle low-level actions related to CDO local
63 quantities such as cell mesh structures or cellwise systems.
64
65 */
66
67 /*=============================================================================
68 * Local type definitions
69 *============================================================================*/
70
71 #define CS_CDO_LOCAL_DBG 0
72
73 /* Redefined names of function from cs_math to get shorter names */
74 #define _dp3 cs_math_3_dot_product
75
76 /*============================================================================
77 * Global variables
78 *============================================================================*/
79
80 /* Auxiliary buffers for extra-operations related to local problems. These
81 * buffers are also used for computing quantities related to a cs_cell_mesh_t
82 * (there are as many buffers as threads since a call to these buffers can be
83 * inside an OpenMP directive */
84 int cs_cdo_local_d_buffer_size = 0;
85 double **cs_cdo_local_d_buffer = NULL;
86
87 /* Pointer of pointers to global structures */
88 cs_cell_mesh_t **cs_cdo_local_cell_meshes = NULL;
89 cs_face_mesh_t **cs_cdo_local_face_meshes = NULL;
90 cs_face_mesh_light_t **cs_cdo_local_face_meshes_light = NULL;
91
92 /*============================================================================
93 * Local static variables
94 *============================================================================*/
95
96 static const int n_robin_parameters = 3;
97 static int cs_cdo_local_n_structures = 0;
98
99 /* Auxiliary buffers for computing quantities related to a cs_cell_mesh_t
100 (there are as many buffers as threads since a call to these buffers can be
101 inside an OpenMP directive */
102 static short int **cs_cdo_local_kbuf = NULL;
103
104 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
105
106 /*============================================================================
107 * Private function prototypes
108 *============================================================================*/
109
110 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
111
112 /*============================================================================
113 * Public function prototypes
114 *============================================================================*/
115
116 /*----------------------------------------------------------------------------*/
117 /*!
118 * \brief Allocate global structures used for build system with a cellwise or
119 * facewise process
120 *
121 * \param[in] connect pointer to a \ref cs_cdo_connect_t structure
122 */
123 /*----------------------------------------------------------------------------*/
124
125 void
cs_cdo_local_initialize(const cs_cdo_connect_t * connect)126 cs_cdo_local_initialize(const cs_cdo_connect_t *connect)
127 {
128 /* Sanity check */
129 assert(cs_glob_n_threads > 0);
130
131 int nthr = cs_glob_n_threads;
132 int n_vc = connect->n_max_vbyc;
133 int max_ent = 3*CS_MAX(n_vc, CS_MAX(connect->n_max_ebyc,
134 connect->n_max_fbyc));
135
136 cs_cdo_local_d_buffer_size = CS_MAX(n_vc*(n_vc+1)/2, max_ent);
137
138 cs_cdo_local_n_structures = nthr;
139 BFT_MALLOC(cs_cdo_local_cell_meshes, nthr, cs_cell_mesh_t *);
140 BFT_MALLOC(cs_cdo_local_face_meshes, nthr, cs_face_mesh_t *);
141 BFT_MALLOC(cs_cdo_local_face_meshes_light, nthr, cs_face_mesh_light_t *);
142 BFT_MALLOC(cs_cdo_local_d_buffer, nthr, double *);
143 BFT_MALLOC(cs_cdo_local_kbuf, nthr, short int *);
144
145 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
146 #pragma omp parallel
147 {
148 int t_id = omp_get_thread_num();
149 assert(t_id < cs_glob_n_threads);
150
151 cs_cdo_local_cell_meshes[t_id] = cs_cell_mesh_create(connect);
152 cs_cdo_local_face_meshes[t_id] = cs_face_mesh_create(connect->n_max_vbyf);
153 cs_cdo_local_face_meshes_light[t_id] =
154 cs_face_mesh_light_create(connect->n_max_vbyf, connect->n_max_vbyc);
155
156 BFT_MALLOC(cs_cdo_local_d_buffer[t_id], cs_cdo_local_d_buffer_size, double);
157 BFT_MALLOC(cs_cdo_local_kbuf[t_id],
158 CS_MAX(connect->v_max_cell_range, connect->e_max_cell_range)+1,
159 short int);
160 }
161 #else
162
163 assert(cs_glob_n_threads == 1);
164
165 cs_cdo_local_cell_meshes[0] = cs_cell_mesh_create(connect);
166 cs_cdo_local_face_meshes[0] = cs_face_mesh_create(connect->n_max_vbyf);
167 cs_cdo_local_face_meshes_light[0] =
168 cs_face_mesh_light_create(connect->n_max_vbyf, connect->n_max_vbyc);
169
170 BFT_MALLOC(cs_cdo_local_d_buffer[0], cs_cdo_local_d_buffer_size, double);
171 BFT_MALLOC(cs_cdo_local_kbuf[0],
172 CS_MAX(connect->v_max_cell_range, connect->e_max_cell_range)+1,
173 short int);
174
175 #endif /* openMP */
176 }
177
178 /*----------------------------------------------------------------------------*/
179 /*!
180 * \brief Free global structures related to \ref cs_cell_mesh_t and
181 * \ref cs_face_mesh_t structures
182 */
183 /*----------------------------------------------------------------------------*/
184
185 void
cs_cdo_local_finalize(void)186 cs_cdo_local_finalize(void)
187 {
188 if (cs_cdo_local_n_structures < 1)
189 return;
190
191 assert(cs_cdo_local_n_structures == cs_glob_n_threads);
192
193 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
194 #pragma omp parallel
195 {
196 int t_id = omp_get_thread_num();
197 assert(t_id < cs_glob_n_threads);
198
199 cs_cell_mesh_free(&(cs_cdo_local_cell_meshes[t_id]));
200 cs_face_mesh_free(&(cs_cdo_local_face_meshes[t_id]));
201 cs_face_mesh_light_free(&(cs_cdo_local_face_meshes_light[t_id]));
202 BFT_FREE(cs_cdo_local_d_buffer[t_id]);
203 BFT_FREE(cs_cdo_local_kbuf[t_id]);
204
205 }
206 #else
207 assert(cs_glob_n_threads == 1);
208 cs_cell_mesh_free(&(cs_cdo_local_cell_meshes[0]));
209 cs_face_mesh_free(&(cs_cdo_local_face_meshes[0]));
210 cs_face_mesh_light_free(&(cs_cdo_local_face_meshes_light[0]));
211 BFT_FREE(cs_cdo_local_d_buffer[0]);
212 BFT_FREE(cs_cdo_local_kbuf[0]);
213 #endif /* openMP */
214
215 BFT_FREE(cs_cdo_local_cell_meshes);
216 BFT_FREE(cs_cdo_local_face_meshes);
217 BFT_FREE(cs_cdo_local_face_meshes_light);
218 BFT_FREE(cs_cdo_local_d_buffer);
219 BFT_FREE(cs_cdo_local_kbuf);
220 }
221
222 /*----------------------------------------------------------------------------*/
223 /*!
224 * \brief Allocate a \ref cs_cell_sys_t structure
225 *
226 * \param[in] n_max_dofbyc max number of entries
227 * \param[in] n_max_fbyc max number of faces in a cell
228 * \param[in] n_blocks number of blocks in a row/column
229 * \param[in] block_sizes size of each block or NULL.
230 * Specific treatment n_blocks = 1.
231 *
232 * \return a pointer to a new allocated \ref cs_cell_sys_t structure
233 */
234 /*----------------------------------------------------------------------------*/
235
236 cs_cell_sys_t *
cs_cell_sys_create(int n_max_dofbyc,int n_max_fbyc,int n_blocks,int * block_sizes)237 cs_cell_sys_create(int n_max_dofbyc,
238 int n_max_fbyc,
239 int n_blocks,
240 int *block_sizes)
241 {
242 cs_cell_sys_t *csys = NULL;
243
244 BFT_MALLOC(csys, 1, cs_cell_sys_t);
245
246 /* Metadata about DoFs */
247
248 csys->c_id = -1;
249 csys->n_dofs = 0;
250 csys->dof_ids = NULL;
251 csys->dof_flag = NULL;
252
253 /* System and previous values */
254
255 csys->mat = NULL;
256 csys->rhs = NULL;
257 csys->source = NULL;
258 csys->val_n = NULL;
259 csys->val_nm1 = NULL;
260
261 /* Internal enforcement */
262
263 csys->has_internal_enforcement = false;
264 csys->dof_is_forced = NULL;
265
266 if (n_max_dofbyc > 0)
267 BFT_MALLOC(csys->dof_is_forced, n_max_dofbyc, bool);
268
269 /* Boundary conditions */
270
271 csys->n_bc_faces = 0;
272 csys->_f_ids = NULL;
273 csys->bf_ids = NULL;
274 csys->bf_flag = NULL;
275
276 csys->has_dirichlet = false;
277 csys->dir_values = NULL;
278
279 csys->has_nhmg_neumann = false;
280 csys->neu_values = NULL;
281
282 csys->has_robin = false;
283 csys->rob_values = NULL;
284
285 csys->has_sliding = false;
286
287 if (n_max_fbyc > 0) {
288
289 BFT_MALLOC(csys->bf_flag, n_max_fbyc, cs_flag_t);
290 memset(csys->bf_flag, 0, sizeof(cs_flag_t)*n_max_fbyc);
291
292 BFT_MALLOC(csys->_f_ids, n_max_fbyc, short int);
293 memset(csys->_f_ids, 0, sizeof(short int)*n_max_fbyc);
294
295 BFT_MALLOC(csys->bf_ids, n_max_fbyc, cs_lnum_t);
296 memset(csys->bf_ids, 0, sizeof(cs_lnum_t)*n_max_fbyc);
297
298 }
299
300 if (n_max_dofbyc > 0) {
301
302 BFT_MALLOC(csys->dof_flag, n_max_dofbyc, cs_flag_t);
303 memset(csys->dof_flag, 0, sizeof(cs_flag_t)*n_max_dofbyc);
304
305 BFT_MALLOC(csys->dof_ids, n_max_dofbyc, cs_lnum_t);
306 memset(csys->dof_ids, 0, sizeof(cs_lnum_t)*n_max_dofbyc);
307
308 if (block_sizes == NULL)
309 csys->mat = cs_sdm_square_create(n_max_dofbyc);
310
311 else if (n_blocks == 1) {
312
313 assert(block_sizes != NULL);
314 if (block_sizes[0] == 3) {
315 int n_row_blocks = n_max_dofbyc/3;
316 assert(n_max_dofbyc % 3 == 0);
317 csys->mat = cs_sdm_block33_create(n_row_blocks, n_row_blocks);
318 }
319 else
320 bft_error(__FILE__, __LINE__, 0,
321 "%s: Invalid initialization of the cellwise block matrix\n",
322 __func__);
323
324 }
325 else
326 csys->mat = cs_sdm_block_create(n_blocks, n_blocks,
327 block_sizes,
328 block_sizes);
329
330 BFT_MALLOC(csys->rhs , n_max_dofbyc, double);
331 BFT_MALLOC(csys->source , n_max_dofbyc, double);
332 BFT_MALLOC(csys->val_n , n_max_dofbyc, double);
333 BFT_MALLOC(csys->val_nm1 , n_max_dofbyc, double);
334 BFT_MALLOC(csys->dir_values, n_max_dofbyc, double);
335 BFT_MALLOC(csys->neu_values, n_max_dofbyc, double);
336
337 const size_t s = n_max_dofbyc * sizeof(double);
338
339 memset(csys->rhs , 0, s);
340 memset(csys->source , 0, s);
341 memset(csys->val_n , 0, s);
342 memset(csys->val_nm1 , 0, s);
343 memset(csys->dir_values, 0, s);
344 memset(csys->neu_values, 0, s);
345 }
346
347 int n_rob_size = n_robin_parameters*CS_MAX(n_max_dofbyc, n_max_fbyc);
348 BFT_MALLOC(csys->rob_values, n_rob_size, double);
349 memset(csys->rob_values, 0, n_rob_size*sizeof(cs_real_t));
350
351 return csys;
352 }
353
354 /*----------------------------------------------------------------------------*/
355 /*!
356 * \brief Reset all members related to BC and some other ones in a
357 * \ref cs_cell_sys_t structure
358 *
359 * \param[in] n_fbyc number of faces in a cell
360 * \param[in, out] csys pointer to the \ref cs_cell_sys_t struct to reset
361 */
362 /*----------------------------------------------------------------------------*/
363
364 void
cs_cell_sys_reset(int n_fbyc,cs_cell_sys_t * csys)365 cs_cell_sys_reset(int n_fbyc,
366 cs_cell_sys_t *csys)
367 {
368 if (n_fbyc < 1 || csys->n_dofs < 1)
369 return;
370
371 const size_t s = csys->n_dofs * sizeof(double);
372
373 memset(csys->rhs, 0, s);
374 memset(csys->source, 0, s);
375
376 csys->has_internal_enforcement = false;
377 for (int i = 0; i < csys->n_dofs; i++)
378 csys->dof_is_forced[i] = false; /* Not selected */
379
380 memset(csys->dof_flag, 0, sizeof(cs_flag_t)*csys->n_dofs);
381
382 csys->n_bc_faces = 0;
383 csys->has_dirichlet = csys->has_nhmg_neumann = false;
384 csys->has_robin = csys->has_sliding = false;
385
386 memset(csys->bf_flag , 0, sizeof(cs_flag_t)*n_fbyc);
387 memset(csys->_f_ids , 0, sizeof(short int)*n_fbyc);
388 memset(csys->bf_ids , 0, sizeof(cs_lnum_t)*n_fbyc);
389
390 memset(csys->dir_values, 0, s);
391 memset(csys->neu_values, 0, s);
392 memset(csys->rob_values, 0,
393 CS_MAX(n_fbyc, csys->n_dofs)*sizeof(double)*n_robin_parameters);
394 }
395
396 /*----------------------------------------------------------------------------*/
397 /*!
398 * \brief Free a \ref cs_cell_sys_t structure
399 *
400 * \param[in, out] p_csys pointer of pointer to a \ref cs_cell_sys_t struct
401 */
402 /*----------------------------------------------------------------------------*/
403
404 void
cs_cell_sys_free(cs_cell_sys_t ** p_csys)405 cs_cell_sys_free(cs_cell_sys_t **p_csys)
406 {
407 cs_cell_sys_t *csys = *p_csys;
408
409 if (csys == NULL)
410 return;
411
412 BFT_FREE(csys->dof_ids);
413 BFT_FREE(csys->dof_flag);
414
415 csys->mat = cs_sdm_free(csys->mat);
416
417 BFT_FREE(csys->rhs);
418 BFT_FREE(csys->source);
419 BFT_FREE(csys->val_n);
420 BFT_FREE(csys->val_nm1);
421
422 BFT_FREE(csys->_f_ids);
423 BFT_FREE(csys->bf_ids);
424 BFT_FREE(csys->bf_flag);
425 BFT_FREE(csys->dir_values);
426 BFT_FREE(csys->neu_values);
427 BFT_FREE(csys->rob_values);
428
429 BFT_FREE(csys->dof_is_forced);
430
431 BFT_FREE(csys);
432 *p_csys= NULL;
433 }
434
435 /*----------------------------------------------------------------------------*/
436 /*!
437 * \brief Dump a local system for debugging purpose
438 *
439 * \param[in] msg associated message to print
440 * \param[in] csys pointer to a \ref cs_cell_sys_t structure
441 */
442 /*----------------------------------------------------------------------------*/
443
444 void
cs_cell_sys_dump(const char msg[],const cs_cell_sys_t * csys)445 cs_cell_sys_dump(const char msg[],
446 const cs_cell_sys_t *csys)
447 {
448 # pragma omp critical
449 {
450 bft_printf( "[rank:%d] %s\n", cs_glob_rank_id, msg);
451
452 if (csys->has_dirichlet || csys->has_nhmg_neumann || csys->has_robin ||
453 csys->has_sliding) {
454
455 bft_printf(">> dirichlet:%s | nhmg_neumann:%s | robin:%s | sliding:%s\n",
456 cs_base_strtf(csys->has_dirichlet),
457 cs_base_strtf(csys->has_nhmg_neumann),
458 cs_base_strtf(csys->has_robin),
459 cs_base_strtf(csys->has_sliding));
460 if (csys->n_bc_faces > 0) {
461 bft_printf(">> Boundary faces\n"
462 ">> %-8s | %-8s | %-6s\n", "_ID", "ID", "FLAG");
463 for (int i = 0; i < csys->n_bc_faces; i++) {
464 short int f = csys->_f_ids[i];
465 bft_printf(">> %8d | %8ld | %6d\n",
466 f, (long)csys->bf_ids[f], csys->bf_flag[f]);
467 }
468 }
469
470 } /* At least one kind of boundary conditions */
471
472 if (csys->mat->flag & CS_SDM_BY_BLOCK)
473 cs_sdm_block_dump(csys->c_id, csys->mat);
474 else
475 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, csys->mat);
476
477 bft_printf(">> %-8s | %-6s | %-10s | %-10s | %-10s | %-8s |"
478 " %-10s | %-10s\n",
479 "IDS", "FLAG", "RHS", "TS", "DIR_VALS", "ENFORCED", "VAL_N",
480 "VAL_N-1");
481 for (int i = 0; i < csys->n_dofs; i++)
482 bft_printf(">> %8ld | %6d | % -.3e | % -.3e | % -.3e |"
483 " %-8s | % -.3e | % -.3e\n",
484 (long)csys->dof_ids[i], csys->dof_flag[i], csys->rhs[i],
485 csys->source[i], csys->dir_values[i],
486 cs_base_strtf(csys->dof_is_forced[i]),
487 csys->val_n[i], csys->val_nm1[i]);
488 }
489 }
490
491 /*----------------------------------------------------------------------------*/
492 /*!
493 * \brief Allocate \ref cs_cell_builder_t structure
494 *
495 * \return a pointer to the new allocated \ref cs_cell_builder_t structure
496 */
497 /*----------------------------------------------------------------------------*/
498
499 cs_cell_builder_t *
cs_cell_builder_create(void)500 cs_cell_builder_create(void)
501 {
502 cs_cell_builder_t *cb = NULL;
503
504 /* Common part to all discretization */
505
506 BFT_MALLOC(cb, 1, cs_cell_builder_t);
507
508 cb->t_pty_eval = 0.;
509 cb->t_bc_eval = 0.;
510 cb->t_st_eval = 0.;
511
512 cb->cell_flag = 0;
513
514 cb->gpty_val = 1; /* grad-div property */
515 cb->tpty_val = 1; /* time property */
516 cb->rpty_val = 1; /* reaction property */
517
518 for (int r = 0; r < CS_CDO_N_MAX_REACTIONS; r++) cb->rpty_vals[r] = 1;
519
520 cb->adv_fluxes = NULL;
521 cb->ids = NULL;
522 cb->values = NULL;
523 cb->vectors = NULL;
524
525 /* Local matrices */
526
527 cb->loc = cb->aux = NULL;
528
529 return cb;
530 }
531
532 /*----------------------------------------------------------------------------*/
533 /*!
534 * \brief Free a \ref cs_cell_builder_t structure
535 *
536 * \param[in, out] p_cb pointer of pointer to a \ref cs_cell_builder_t struct
537 */
538 /*----------------------------------------------------------------------------*/
539
540 void
cs_cell_builder_free(cs_cell_builder_t ** p_cb)541 cs_cell_builder_free(cs_cell_builder_t **p_cb)
542 {
543 cs_cell_builder_t *cb = *p_cb;
544
545 if (cb == NULL)
546 return;
547
548 BFT_FREE(cb->adv_fluxes);
549 BFT_FREE(cb->ids);
550 BFT_FREE(cb->values);
551 BFT_FREE(cb->vectors);
552
553 cb->loc = cs_sdm_free(cb->loc);
554 cb->aux = cs_sdm_free(cb->aux);
555
556 BFT_FREE(cb);
557 *p_cb = NULL;
558 }
559
560 /*----------------------------------------------------------------------------*/
561 /*!
562 * \brief Allocate and initialize a cs_cell_mesh_t structure
563 *
564 * \param[in] connect pointer to a cs_cdo_connect_t structure
565 *
566 * \return a pointer to a new allocated cs_cell_mesh_t structure
567 */
568 /*----------------------------------------------------------------------------*/
569
570 cs_cell_mesh_t *
cs_cell_mesh_create(const cs_cdo_connect_t * connect)571 cs_cell_mesh_create(const cs_cdo_connect_t *connect)
572 {
573 cs_cell_mesh_t *cm = NULL;
574
575 BFT_MALLOC(cm, 1, cs_cell_mesh_t);
576
577 /* Sizes used to allocate buffers */
578
579 cm->n_max_vbyc = connect->n_max_vbyc;
580 cm->n_max_ebyc = connect->n_max_ebyc;
581 cm->n_max_fbyc = connect->n_max_fbyc;
582
583 cm->flag = 0;
584 cm->n_vc = 0;
585 cm->n_ec = 0;
586 cm->n_fc = 0;
587
588 /* Vertex information */
589
590 BFT_MALLOC(cm->v_ids, cm->n_max_vbyc, cs_lnum_t);
591 BFT_MALLOC(cm->wvc, cm->n_max_vbyc, double);
592 BFT_MALLOC(cm->xv, 3*cm->n_max_vbyc, double);
593
594 /* Edge information */
595
596 BFT_MALLOC(cm->e_ids, cm->n_max_ebyc, cs_lnum_t);
597 BFT_MALLOC(cm->e2v_sgn, cm->n_max_ebyc, short int);
598 BFT_MALLOC(cm->edge, cm->n_max_ebyc, cs_quant_t);
599 BFT_MALLOC(cm->dface, cm->n_max_ebyc, cs_nvec3_t);
600 BFT_MALLOC(cm->pvol_e, cm->n_max_ebyc, double);
601
602 /* Face information */
603
604 BFT_MALLOC(cm->f_ids, cm->n_max_fbyc, cs_lnum_t);
605 BFT_MALLOC(cm->f_sgn, cm->n_max_fbyc, short int);
606 BFT_MALLOC(cm->f_diam, cm->n_max_fbyc, double);
607 BFT_MALLOC(cm->face, cm->n_max_fbyc, cs_quant_t);
608 BFT_MALLOC(cm->dedge, cm->n_max_fbyc, cs_nvec3_t);
609 BFT_MALLOC(cm->hfc, cm->n_max_fbyc, double);
610 BFT_MALLOC(cm->pvol_f, cm->n_max_fbyc, double);
611
612 /* face --> vertices connectivity */
613
614 BFT_MALLOC(cm->f2v_idx, cm->n_max_fbyc + 1, short int);
615 BFT_MALLOC(cm->f2v_ids, 2*cm->n_max_ebyc, short int);
616
617 /* face --> edges connectivity */
618
619 BFT_MALLOC(cm->f2e_idx, cm->n_max_fbyc + 1, short int);
620 BFT_MALLOC(cm->f2e_ids, 2*cm->n_max_ebyc, short int);
621 BFT_MALLOC(cm->f2e_sgn, 2*cm->n_max_ebyc, short int);
622 BFT_MALLOC(cm->tef, 2*cm->n_max_ebyc, double);
623 BFT_MALLOC(cm->sefc, 2*cm->n_max_ebyc, cs_nvec3_t);
624
625 /* edge --> vertices connectivity */
626
627 BFT_MALLOC(cm->e2v_ids, 2*cm->n_max_ebyc, short int);
628
629 /* edge --> face connectivity */
630
631 BFT_MALLOC(cm->e2f_ids, 2*cm->n_max_ebyc, short int);
632
633 cs_cell_mesh_reset(cm);
634
635 return cm;
636 }
637
638 /*----------------------------------------------------------------------------*/
639 /*!
640 * \brief Get a pointer to a cs_cell_mesh_t structure corresponding to mesh id
641 *
642 * \param[in] mesh_id id in the array of pointer to cs_cell_mesh_t struct.
643 *
644 * \return a pointer to a cs_cell_mesh_t structure
645 */
646 /*----------------------------------------------------------------------------*/
647
648 cs_cell_mesh_t *
cs_cdo_local_get_cell_mesh(int mesh_id)649 cs_cdo_local_get_cell_mesh(int mesh_id)
650 {
651 if (mesh_id < 0 || mesh_id >= cs_glob_n_threads)
652 return NULL;
653
654 cs_cell_mesh_t *cm = cs_cdo_local_cell_meshes[mesh_id];
655
656 #if defined(DEBUG) && !defined(NDEBUG)
657 /* This is to check that the mesh flag is correctly set */
658 cs_cell_mesh_reset(cm);
659 #endif
660
661 return cm;
662 }
663
664 /*----------------------------------------------------------------------------*/
665 /*!
666 * \brief Initialize to invalid values a cs_cell_mesh_t structure
667 *
668 * \param[in] cm pointer to a cs_cell_mesh_t structure
669 */
670 /*----------------------------------------------------------------------------*/
671
672 void
cs_cell_mesh_reset(cs_cell_mesh_t * cm)673 cs_cell_mesh_reset(cs_cell_mesh_t *cm)
674 {
675 cm->n_vc = -1;
676 cm->n_ec = -1;
677 cm->n_fc = -1;
678
679 /* Cell information */
680
681 cm->c_id = -1;
682 cm->xc[0] = cm->xc[1] = cm->xc[2] = -DBL_MAX;
683 cm->vol_c = -DBL_MAX;
684 cm->diam_c = -DBL_MAX;
685
686 /* Vertex information */
687
688 for (short int v = 0; v < cm->n_max_vbyc; v++) {
689 cm->v_ids[v] = -1;
690 cm->wvc[v] = -DBL_MAX;
691 cm->xv[3*v] = cm->xv[3*v+1] = cm->xv[3*v+2] = -DBL_MAX;
692 }
693
694 /* Edge information */
695
696 for (short int e = 0; e < cm->n_max_ebyc; e++) {
697 cm->e_ids[e] = -1;
698 cm->e2v_sgn[e] = 0;
699 cm->pvol_e[e] = -DBL_MAX;
700 cm->edge[e].meas = cm->dface[e].meas = -DBL_MAX;
701 cm->edge[e].unitv[0] = cm->dface[e].unitv[0] = -DBL_MAX;
702 cm->edge[e].unitv[1] = cm->dface[e].unitv[1] = -DBL_MAX;
703 cm->edge[e].unitv[2] = cm->dface[e].unitv[2] = -DBL_MAX;
704 cm->edge[e].center[0] = -DBL_MAX;
705 cm->edge[e].center[1] = -DBL_MAX;
706 cm->edge[e].center[2] = -DBL_MAX;
707 }
708
709 /* Face information */
710
711 for (short int f = 0; f < cm->n_max_fbyc; f++) {
712 cm->f_ids[f] = -1;
713 cm->f_sgn[f] = 0;
714 cm->f_diam[f] = -DBL_MAX;
715 cm->hfc[f] = -DBL_MAX;
716 cm->pvol_f[f] = -DBL_MAX;
717 cm->face[f].meas = cm->dedge[f].meas = -DBL_MAX;
718 cm->face[f].unitv[0] = cm->dedge[f].unitv[0] = -DBL_MAX;
719 cm->face[f].unitv[1] = cm->dedge[f].unitv[1] = -DBL_MAX;
720 cm->face[f].unitv[2] = cm->dedge[f].unitv[2] = -DBL_MAX;
721 cm->face[f].center[0] = -DBL_MAX;
722 cm->face[f].center[1] = -DBL_MAX;
723 cm->face[f].center[2] = -DBL_MAX;
724 }
725
726 /* face --> edges connectivity */
727
728 for (short int f = 0; f < cm->n_max_fbyc + 1; f++)
729 cm->f2e_idx[f] = cm->f2v_idx[f] = -1;
730
731 for (int i = 0; i < 2*cm->n_max_ebyc; i++) {
732 cm->e2v_ids[i] = cm->e2f_ids[i] = -1;
733 cm->f2e_ids[i] = cm->f2v_ids[i] = -1;
734 cm->f2e_sgn[i] = 0;
735 cm->tef[i] = cm->sefc[i].meas = -DBL_MAX;
736 cm->sefc[i].unitv[0]=cm->sefc[i].unitv[1]=cm->sefc[i].unitv[2] = -DBL_MAX;
737 }
738 }
739
740 /*----------------------------------------------------------------------------*/
741 /*!
742 * \brief Dump a cs_cell_mesh_t structure
743 *
744 * \param[in] cm pointer to a cs_cell_mesh_t structure
745 */
746 /*----------------------------------------------------------------------------*/
747
748 void
cs_cell_mesh_dump(const cs_cell_mesh_t * cm)749 cs_cell_mesh_dump(const cs_cell_mesh_t *cm)
750 {
751 if (cm == NULL) {
752 bft_printf("\n>> Dump cs_cell_mesh_t %p\n", (const void *)cm);
753 return;
754 }
755
756 bft_printf("\n>> [rank: %d] Dump cs_cell_mesh_t %p; %s; flag: %d\n"
757 " c_id:%ld; vol: %9.6e; xc (% .4e % .4e % .4e); diam: % .4e\n",
758 cs_glob_rank_id, (const void *)cm, fvm_element_type_name[cm->type],
759 cm->flag, (long)cm->c_id, cm->vol_c, cm->xc[0], cm->xc[1],
760 cm->xc[2], cm->diam_c);
761
762 /* Information related to primal vertices */
763
764 if (cm->flag & cs_flag_need_v) {
765
766 bft_printf(" %s | %6s | %35s | %10s\n",
767 "v", "id", "coord", "wvc");
768 for (short int v = 0; v < cm->n_vc; v++)
769 bft_printf("%2d | %6ld | % .4e % .4e % .4e | %.4e\n",
770 v, (long)cm->v_ids[v], cm->xv[3*v], cm->xv[3*v+1],
771 cm->xv[3*v+2], cm->wvc[v]);
772
773 } /* Vertex quantities */
774
775 /* Information related to primal edges */
776
777 if (cm->flag & cs_flag_need_e) {
778
779 bft_printf(" %s | %6s | %3s | %2s | %2s | %9s |"
780 " %35s | %35s | %10s | %35s | %9s\n",
781 "e", "id", "sgn", "v1", "v2", "length", "unit", "coords",
782 "df.meas", "df.unit", "pvol_e");
783 for (short int e = 0; e < cm->n_ec; e++) {
784
785 cs_quant_t peq = cm->edge[e];
786 cs_nvec3_t dfq = cm->dface[e];
787 bft_printf("%2d | %6ld | %3d | %2d | %2d | %.3e |"
788 " % .4e % .4e % .4e | % .4e % .4e % .4e | %.4e |"
789 " % .4e % .4e % .4e | % .4e\n",
790 e, (long)cm->e_ids[e], cm->e2v_sgn[e], cm->e2v_ids[2*e],
791 cm->e2v_ids[2*e+1], peq.meas, peq.unitv[0], peq.unitv[1],
792 peq.unitv[2], peq.center[0], peq.center[1], peq.center[2],
793 dfq.meas, dfq.unitv[0], dfq.unitv[1], dfq.unitv[2],
794 cm->pvol_e[e]);
795
796 } /* Loop on edges */
797
798 } /* Edge quantities */
799
800 /* Information related to primal faces */
801
802 if (cm->flag & cs_flag_need_f) {
803
804 bft_printf(" %s | %6s | %9s | %3s | %35s | %35s |"
805 " %10s | %35s | %11s %11s %11s\n",
806 "f", "id", "surf", "sgn", "unit", "coords", "dlen", "dunitv",
807 "pfc", "hfc", "diam");
808 for (short int f = 0; f < cm->n_fc; f++) {
809 cs_quant_t pfq = cm->face[f];
810 cs_nvec3_t deq = cm->dedge[f];
811 bft_printf("%2d | %6ld | %.3e | %3d | % .4e % .4e % .4e |"
812 " % .4e % .4e % .4e | %.4e | % .4e % .4e % .4e | %.3e |"
813 " %.3e | %.3e\n",
814 f, (long)cm->f_ids[f], pfq.meas, cm->f_sgn[f],
815 pfq.unitv[0], pfq.unitv[1], pfq.unitv[2], pfq.center[0],
816 pfq.center[1], pfq.center[2], deq.meas, deq.unitv[0],
817 deq.unitv[1], deq.unitv[2], cm->pvol_f[f], cm->hfc[f],
818 cm->f_diam[f]);
819 }
820
821 } /* Face quantities */
822
823 if (cm->flag & cs_flag_need_fe) {
824
825 bft_printf(" f | n_ef | e:tef\n");
826 for (short int f = 0; f < cm->n_fc; f++) {
827 bft_printf(" f%2d | %4d |", f, cm->f2e_idx[f+1] - cm->f2e_idx[f]);
828 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++)
829 bft_printf(" e%2d:%.4e (%+1d)",
830 cm->f2e_ids[i], cm->tef[i], cm->f2e_sgn[i]);
831 bft_printf("\n");
832 }
833
834 bft_printf(" e | f0 | sefc ...\n");
835 for (short int e = 0; e < cm->n_ec; e++) {
836 int count = 0;
837 bft_printf(" %2d", e);
838 for (short int f = 0; f < cm->n_fc; f++) {
839 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
840 if (e == cm->f2e_ids[i]) {
841 cs_nvec3_t _s = cm->sefc[i];
842 bft_printf(" | %2d | %.4e (%- .4e %- .4e %- .4e)", f,
843 _s.meas, _s.unitv[0], _s.unitv[1], _s.unitv[2]);
844 count++;
845 }
846 } /* Loop on face edges */
847 if (count == 2)
848 break;
849 } /* Loop on faces */
850
851 bft_printf("\n");
852
853 } /* Loop on edges */
854
855 }
856 }
857
858 /*----------------------------------------------------------------------------*/
859 /*!
860 * \brief Free a cs_cell_mesh_t structure
861 *
862 * \param[in, out] p_cm pointer of pointer to a cs_cell_mesh_t structure
863 */
864 /*----------------------------------------------------------------------------*/
865
866 void
cs_cell_mesh_free(cs_cell_mesh_t ** p_cm)867 cs_cell_mesh_free(cs_cell_mesh_t **p_cm)
868 {
869 cs_cell_mesh_t *cm = *p_cm;
870
871 if (cm == NULL)
872 return;
873
874 BFT_FREE(cm->v_ids);
875 BFT_FREE(cm->wvc);
876 BFT_FREE(cm->xv);
877
878 BFT_FREE(cm->e_ids);
879 BFT_FREE(cm->edge);
880 BFT_FREE(cm->dface);
881 BFT_FREE(cm->pvol_e);
882
883 BFT_FREE(cm->f_ids);
884 BFT_FREE(cm->f_sgn);
885 BFT_FREE(cm->f_diam);
886 BFT_FREE(cm->hfc);
887 BFT_FREE(cm->pvol_f);
888 BFT_FREE(cm->face);
889 BFT_FREE(cm->dedge);
890
891 BFT_FREE(cm->e2v_ids);
892 BFT_FREE(cm->e2v_sgn);
893
894 BFT_FREE(cm->f2v_idx);
895 BFT_FREE(cm->f2v_ids);
896
897 BFT_FREE(cm->f2e_idx);
898 BFT_FREE(cm->f2e_ids);
899 BFT_FREE(cm->f2e_sgn);
900 BFT_FREE(cm->tef);
901
902 BFT_FREE(cm->e2f_ids);
903 BFT_FREE(cm->sefc);
904
905 BFT_FREE(cm);
906 *p_cm = NULL;
907 }
908
909 /*----------------------------------------------------------------------------*/
910 /*!
911 * \brief Define a cs_cell_mesh_t structure for a given cell id. According
912 * to the requested level, some quantities may not be defined;
913 *
914 * \param[in] c_id cell id
915 * \param[in] build_flag indicate which members are really built
916 * \param[in] connect pointer to a cs_cdo_connect_t structure
917 * \param[in] quant pointer to a cs_cdo_quantities_t structure
918 * \param[in, out] cm pointer to a cs_cell_mesh_t structure to set
919 */
920 /*----------------------------------------------------------------------------*/
921
922 void
cs_cell_mesh_build(cs_lnum_t c_id,cs_eflag_t build_flag,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_cell_mesh_t * cm)923 cs_cell_mesh_build(cs_lnum_t c_id,
924 cs_eflag_t build_flag,
925 const cs_cdo_connect_t *connect,
926 const cs_cdo_quantities_t *quant,
927 cs_cell_mesh_t *cm)
928 {
929 if (cm == NULL)
930 return;
931
932 cm->flag = build_flag;
933 cm->type = connect->cell_type[c_id];
934
935 /* Store information related to cell */
936
937 cm->c_id = c_id;
938 cm->vol_c = quant->cell_vol[c_id];
939 for (int k = 0; k < 3; k++)
940 cm->xc[k] = quant->cell_centers[3*c_id+k];
941
942 /* Store the number of cell faces (useful to allocated boundary quantities) */
943
944 const cs_lnum_t *c2f_idx = connect->c2f->idx + c_id;
945
946 cm->n_fc = c2f_idx[1] - c2f_idx[0];
947 cm->bface_shift = quant->n_i_faces; /* Border faces come after interior
948 * faces */
949
950 assert(cm->n_fc > 3); /* A tetrahedron has at least 4 faces */
951
952 if (build_flag == 0)
953 return;
954
955 /* Information related to primal vertices */
956
957 if (build_flag & cs_flag_need_v) {
958
959 const cs_lnum_t *c2v_idx = connect->c2v->idx + c_id;
960 const cs_lnum_t *c2v_ids = connect->c2v->ids + c2v_idx[0];
961
962 cm->n_vc = c2v_idx[1] - c2v_idx[0];
963
964 assert(cm->n_vc > 3); /* A tetrahedron has at least 4 vertices */
965
966 for (short int v = 0; v < cm->n_vc; v++) {
967 const cs_lnum_t v_id = c2v_ids[v];
968 cm->v_ids[v] = v_id;
969 for (int k = 0; k < 3; k++)
970 cm->xv[3*v+k] = quant->vtx_coord[3*v_id+k];
971 }
972
973 /* Primal vertices quantities */
974
975 if (build_flag & CS_FLAG_COMP_PVQ) {
976
977 const double *wvc = quant->dcell_vol + c2v_idx[0];
978 const double invvol = 1/cm->vol_c;
979 for (short int v = 0; v < cm->n_vc; v++)
980 cm->wvc[v] = invvol * wvc[v];
981
982 }
983
984 } /* vertices */
985
986 /* Information related to primal edges */
987
988 if (build_flag & cs_flag_need_e) {
989
990 const cs_lnum_t *c2e_idx = connect->c2e->idx + c_id;
991 const cs_lnum_t *c2e_ids = connect->c2e->ids + c2e_idx[0];
992
993 cm->n_ec = c2e_idx[1] - c2e_idx[0];
994
995 for (short int e = 0; e < cm->n_ec; e++)
996 cm->e_ids[e] = c2e_ids[e];
997
998 if (build_flag & cs_flag_need_peq) {
999
1000 assert(build_flag & CS_FLAG_COMP_PV);
1001
1002 /* Primal edge quantities */
1003
1004 for (short int e = 0; e < cm->n_ec; e++) {
1005
1006 const cs_lnum_t e_id = cm->e_ids[e];
1007 const cs_nvec3_t nv = cs_quant_set_edge_nvec(e_id, quant);
1008 const cs_lnum_t *v_id = connect->e2v->ids + 2*e_id;
1009 const cs_real_t *xv1 = quant->vtx_coord + 3*v_id[0];
1010 const cs_real_t *xv2 = quant->vtx_coord + 3*v_id[1];
1011
1012 cm->edge[e].meas = nv.meas;
1013 for (int k = 0; k < 3; k++) {
1014 cm->edge[e].center[k] = 0.5 * (xv1[k] + xv2[k]);
1015 cm->edge[e].unitv[k] = nv.unitv[k];
1016 }
1017
1018 }
1019
1020 } /* Primal edge quantities */
1021
1022 /* Dual face quantities related to each edge */
1023
1024 if (build_flag & cs_flag_need_dfq) {
1025
1026 const cs_real_t *dface = quant->dface_normal + 3*c2e_idx[0];
1027
1028 for (short int e = 0; e < cm->n_ec; e++) {
1029
1030 cs_nvec3_t df_nvect;
1031 cs_nvec3(dface + 3*e, &df_nvect);
1032 cm->dface[e].meas = df_nvect.meas;
1033 for (int k = 0; k < 3; k++)
1034 cm->dface[e].unitv[k] = df_nvect.unitv[k];
1035
1036 }
1037
1038 } /* Dual face quantities */
1039
1040 if (build_flag & CS_FLAG_COMP_PEC) {
1041
1042 assert(quant->pvol_ec != NULL);
1043 const cs_real_t *_pvol = quant->pvol_ec + c2e_idx[0];
1044 for (short int e = 0; e < cm->n_ec; e++)
1045 cm->pvol_e[e] = _pvol[e];
1046
1047 }
1048
1049 } /* Edge-related quantities */
1050
1051 /* Information related to primal faces */
1052
1053 if (build_flag & cs_flag_need_f) {
1054
1055 const cs_lnum_t *c2f_ids = connect->c2f->ids + c2f_idx[0];
1056 const short int *c2f_sgn = connect->c2f->sgn + c2f_idx[0];
1057
1058 for (short int f = 0; f < cm->n_fc; f++) {
1059 cm->f_ids[f] = c2f_ids[f];
1060 cm->f_sgn[f] = c2f_sgn[f];
1061 }
1062
1063 /* Face related quantities */
1064
1065 if (build_flag & cs_flag_need_pfq) {
1066
1067 for (short int f = 0; f < cm->n_fc; f++) {
1068
1069 const cs_quant_t pfq = cs_quant_set_face(cm->f_ids[f], quant);
1070
1071 cm->face[f].meas = pfq.meas;
1072 for (int k = 0; k < 3; k++) {
1073 cm->face[f].center[k] = pfq.center[k];
1074 cm->face[f].unitv[k] = pfq.unitv[k];
1075 }
1076
1077 }
1078
1079 } /* Primal face quantities */
1080
1081 if (build_flag & cs_flag_need_deq) {
1082
1083 for (short int f = 0; f < cm->n_fc; f++) {
1084
1085 const cs_nvec3_t de_nvect = cs_quant_set_dedge_nvec(c2f_idx[0]+f,
1086 quant);
1087
1088 /* Copy cs_nvec3_t structure */
1089
1090 cm->dedge[f].meas = de_nvect.meas;
1091 for (int k = 0; k < 3; k++)
1092 cm->dedge[f].unitv[k] = de_nvect.unitv[k];
1093
1094 }
1095
1096 } /* Dual edge quantities */
1097
1098 if (build_flag & cs_flag_need_pfc) {
1099
1100 if (quant->pvol_fc != NULL) {
1101
1102 const cs_real_t *_pvol = quant->pvol_fc + c2f_idx[0];
1103 for (short int f = 0; f < cm->n_fc; f++)
1104 cm->pvol_f[f] = _pvol[f];
1105
1106 }
1107 else {
1108
1109 assert(cs_eflag_test(build_flag,
1110 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1111 for (short int f = 0; f < cm->n_fc; f++) {
1112 cm->pvol_f[f] = _dp3(cm->face[f].unitv, cm->dedge[f].unitv);
1113 cm->pvol_f[f] *= cs_math_1ov3 * cm->face[f].meas * cm->dedge[f].meas;
1114 }
1115
1116 }
1117
1118 }
1119
1120 if (build_flag & CS_FLAG_COMP_HFQ) {
1121
1122 /* Compute the height of the pyramid of base f whose apex is
1123 the cell center */
1124
1125 for (short int f = 0; f < cm->n_fc; f++)
1126 cm->hfc[f] = 3 * cm->pvol_f[f] / cm->face[f].meas;
1127
1128 } /* Quantities related to the pyramid of base f */
1129
1130 } /* Face information */
1131
1132 if (build_flag & CS_FLAG_COMP_EV || build_flag & CS_FLAG_COMP_FV) {
1133
1134 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1135 int t_id = omp_get_thread_num();
1136 assert(t_id < cs_glob_n_threads);
1137 #else
1138 int t_id = 0;
1139 #endif /* openMP */
1140
1141 short int *kbuf = cs_cdo_local_kbuf[t_id];
1142
1143 /* Store in a compact way ids for vertices: mesh --> cell mesh */
1144
1145 cs_lnum_t v_shift = cm->v_ids[0];
1146 for (short int v = 1; v < cm->n_vc; v++)
1147 if (cm->v_ids[v] < v_shift)
1148 v_shift = cm->v_ids[v];
1149 for (short int v = 0; v < cm->n_vc; v++)
1150 kbuf[cm->v_ids[v]-v_shift] = v;
1151
1152 if (build_flag & CS_FLAG_COMP_EV) {
1153 for (short int e = 0; e < cm->n_ec; e++) {
1154
1155 const cs_lnum_t e_id = cm->e_ids[e];
1156
1157 /* Store only the sign related to the first vertex since the sign
1158 related to the second one is minus the first one */
1159
1160 cm->e2v_sgn[e] = connect->e2v->sgn[2*e_id];
1161 cm->e2v_ids[2*e] = kbuf[connect->e2v->ids[2*e_id] - v_shift];
1162 cm->e2v_ids[2*e+1] = kbuf[connect->e2v->ids[2*e_id+1] - v_shift];
1163
1164 } /* Loop on cell edges */
1165 } /* edge-vertices information */
1166
1167 if (build_flag & CS_FLAG_COMP_FV) {
1168
1169 /* Build the index and the list of face vertices in the cellwise
1170 * numbering */
1171
1172 cm->f2v_idx[0] = 0;
1173 for (short int f = 0; f < cm->n_fc; f++) {
1174
1175 const cs_lnum_t f_id = cm->f_ids[f];
1176 if (f_id < cm->bface_shift) { /* Interior face */
1177
1178 const cs_lnum_t *idx = connect->if2v->idx + f_id;
1179 const cs_lnum_t n_ent = idx[1] - idx[0];
1180 const cs_lnum_t *v_ids = connect->if2v->ids + idx[0];
1181
1182 cm->f2v_idx[f+1] = cm->f2v_idx[f] + n_ent;
1183
1184 short int *_ids = cm->f2v_ids + cm->f2v_idx[f];
1185 for (cs_lnum_t i = 0; i < n_ent; i++)
1186 _ids[i] = kbuf[v_ids[i] - v_shift];
1187
1188 }
1189 else { /* Border face */
1190
1191 const cs_lnum_t bf_id = f_id - cm->bface_shift;
1192 const cs_lnum_t *idx = connect->bf2v->idx + bf_id;
1193 const cs_lnum_t n_ent = idx[1] - idx[0];
1194 const cs_lnum_t *v_ids = connect->bf2v->ids + idx[0];
1195
1196 cm->f2v_idx[f+1] = cm->f2v_idx[f] + n_ent;
1197
1198 short int *_ids = cm->f2v_ids + cm->f2v_idx[f];
1199 for (cs_lnum_t i = 0; i < n_ent; i++)
1200 _ids[i] = kbuf[v_ids[i] - v_shift];
1201
1202 }
1203
1204 } /* Loop on cell faces */
1205
1206 assert(cm->f2v_idx[cm->n_fc] == 2*cm->n_ec);
1207
1208 } /* face --> vertices connectivity */
1209
1210 } /* EV or FV flag */
1211
1212 if (build_flag & cs_flag_need_fe) {
1213
1214 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1215 int t_id = omp_get_thread_num();
1216 assert(t_id < cs_glob_n_threads);
1217 #else
1218 int t_id = 0;
1219 #endif /* openMP */
1220
1221 short int *kbuf = cs_cdo_local_kbuf[t_id];
1222
1223 /* Store in compact way: mesh --> cell mesh ids for edges */
1224
1225 cs_lnum_t shift = cm->e_ids[0];
1226 for (short int e = 1; e < cm->n_ec; e++)
1227 if (cm->e_ids[e] < shift)
1228 shift = cm->e_ids[e];
1229 for (short int e = 0; e < cm->n_ec; e++)
1230 kbuf[cm->e_ids[e]-shift] = e;
1231
1232 const cs_adjacency_t *f2e = connect->f2e;
1233
1234 cm->f2e_idx[0] = 0;
1235 if (build_flag & CS_FLAG_COMP_FES) {
1236
1237 for (short int f = 0; f < cm->n_fc; f++) {
1238
1239 const cs_lnum_t *idx = f2e->idx + cm->f_ids[f];
1240 const cs_lnum_t *ids = f2e->ids + idx[0];
1241 const short int *sgn = f2e->sgn + idx[0];
1242 const cs_lnum_t n_ef = idx[1]-idx[0];
1243
1244 short int *_ids = cm->f2e_ids + cm->f2e_idx[f];
1245 short int *_sgn = cm->f2e_sgn + cm->f2e_idx[f];
1246
1247 cm->f2e_idx[f+1] = cm->f2e_idx[f] + n_ef;
1248 for (cs_lnum_t i = 0; i < n_ef; i++) {
1249 _ids[i] = kbuf[ids[i] - shift]; /* cellwise numbering */
1250 _sgn[i] = sgn[i];
1251 }
1252
1253 } /* Loop on cell faces */
1254
1255 }
1256 else {
1257
1258 for (short int f = 0; f < cm->n_fc; f++) {
1259
1260 const cs_lnum_t *idx = f2e->idx + cm->f_ids[f];
1261 const cs_lnum_t *ids = f2e->ids + idx[0];
1262 const cs_lnum_t n_ef = idx[1]-idx[0];
1263
1264 short int *_ids = cm->f2e_ids + cm->f2e_idx[f];
1265
1266 cm->f2e_idx[f+1] = cm->f2e_idx[f] + n_ef;
1267 for (cs_lnum_t i = 0; i < n_ef; i++)
1268 _ids[i] = kbuf[ids[i] - shift]; /* cellwise numbering */
1269
1270 } /* Loop on cell faces */
1271
1272 } /* Build f2e_sgn ? */
1273
1274 #if defined(DEBUG) && !defined(NDEBUG)
1275 /* Sanity check */
1276 if (quant->remove_boundary_faces == false)
1277 if (cm->f2e_idx[cm->n_fc] != 2*cm->n_ec)
1278 bft_error(__FILE__, __LINE__, 0,
1279 " %s: Inconsistency detected in f2e_idx for c_id = %d\n"
1280 " cm->f2e_idx[cm->n_fc] = %d and 2*cm->n_ec = %d\n",
1281 __func__, cm->c_id, cm->f2e_idx[cm->n_fc], 2*cm->n_ec);
1282 #endif
1283
1284 if (build_flag & CS_FLAG_COMP_FEQ) {
1285
1286 for (short int f = 0; f < cm->n_fc; f++) {
1287 for (int ie = cm->f2e_idx[f]; ie < cm->f2e_idx[f+1]; ie++)
1288 cm->tef[ie] = cs_compute_area_from_quant(cm->edge[cm->f2e_ids[ie]],
1289 cm->face[f].center);
1290 }
1291
1292 } /* face --> edges quantities */
1293
1294 if (build_flag & CS_FLAG_COMP_SEF) { /* Build cm->sefc */
1295
1296 for (short int f = 0; f < cm->n_fc; f++) {
1297
1298 const cs_quant_t pfq = cm->face[f];
1299 const cs_nvec3_t deq = cm->dedge[f];
1300 const cs_real_t de_coef = 0.5*deq.meas;
1301
1302 /* Loop on face edges */
1303
1304 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1305
1306 const cs_quant_t peq = cm->edge[cm->f2e_ids[i]];
1307
1308 /* Vector area : 0.5 (x_f - x_e) cross.prod (x_f - x_c) */
1309
1310 const cs_real_3_t xexf = { pfq.center[0] - peq.center[0],
1311 pfq.center[1] - peq.center[1],
1312 pfq.center[2] - peq.center[2] };
1313 cs_real_3_t vect_area
1314 = { deq.unitv[1]*xexf[2] - deq.unitv[2]*xexf[1],
1315 deq.unitv[2]*xexf[0] - deq.unitv[0]*xexf[2],
1316 deq.unitv[0]*xexf[1] - deq.unitv[1]*xexf[0] };
1317
1318 if (_dp3(peq.unitv, vect_area) > 0)
1319 for (int k = 0; k < 3; k++) vect_area[k] *= de_coef;
1320 else
1321 for (int k = 0; k < 3; k++) vect_area[k] *= -de_coef;
1322
1323 cs_nvec3(vect_area, cm->sefc + i);
1324
1325 } /* Loop on face edges */
1326
1327 } /* Loop on cell faces */
1328
1329 } /* sefc quantities */
1330
1331 } /* face --> edges connectivity */
1332
1333 if (build_flag & cs_flag_need_ef) {
1334
1335 /* Build the e2f connectivity */
1336
1337 for (short int i = 0; i < 2*cm->n_ec; i++) cm->e2f_ids[i] = -1;
1338
1339 for (short int f = 0; f < cm->n_fc; f++) {
1340 for (short int ie = cm->f2e_idx[f]; ie < cm->f2e_idx[f+1]; ie++) {
1341
1342 short int *ids = cm->e2f_ids + 2*cm->f2e_ids[ie];
1343 if (ids[0] == -1)
1344 ids[0] = f;
1345 else {
1346 assert(ids[1] == -1);
1347 ids[1] = f;
1348 }
1349
1350 } /* Loop on face edges */
1351 } /* Loop on cell faces */
1352
1353 } /* edge-->faces */
1354
1355 if (build_flag & CS_FLAG_COMP_DIAM) {
1356
1357 if (cm->n_fc == 4) { /* Tetrahedron */
1358
1359 assert(cs_flag_test(build_flag, CS_FLAG_COMP_PEQ | CS_FLAG_COMP_FE));
1360
1361 cs_real_t perim_surf = 0.;
1362 for (short int f = 0; f < cm->n_fc; f++) {
1363
1364 cs_real_t perim = 0.;
1365 for (int j = cm->f2e_idx[f]; j < cm->f2e_idx[f+1]; j++) {
1366
1367 short int e = cm->f2e_ids[j];
1368 perim += cm->edge[e].meas;
1369
1370 } /* Loop on edges */
1371
1372 cm->f_diam[f] = 4*cm->face[f].meas/perim;
1373 perim_surf += cm->face[f].meas;
1374
1375 } /* Loop on faces */
1376
1377 cm->diam_c = 6*cm->vol_c/perim_surf;
1378
1379 }
1380 else {
1381
1382 assert(cs_flag_test(build_flag, CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
1383
1384 double cbox[6] = {cm->xv[0], cm->xv[1], cm->xv[2],
1385 cm->xv[0], cm->xv[1], cm->xv[2]};
1386
1387 for (int v = 1; v < cm->n_vc; v++) {
1388 const double *xv = cm->xv + 3*v;
1389 if (xv[0] < cbox[0]) cbox[0] = xv[0];
1390 if (xv[1] < cbox[1]) cbox[1] = xv[1];
1391 if (xv[2] < cbox[2]) cbox[2] = xv[2];
1392 if (xv[0] > cbox[3]) cbox[3] = xv[0];
1393 if (xv[1] > cbox[4]) cbox[4] = xv[1];
1394 if (xv[2] > cbox[5]) cbox[5] = xv[2];
1395 }
1396 cm->diam_c = cs_math_3_distance(cbox, cbox + 3);
1397
1398 /* Now compute an approximation of the diameter for each cell face */
1399
1400 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1401 int t_id = omp_get_thread_num();
1402 assert(t_id < cs_glob_n_threads);
1403 #else
1404 int t_id = 0;
1405 #endif /* openMP */
1406 short int *vtag = cs_cdo_local_kbuf[t_id];
1407
1408 for (short int f = 0; f < cm->n_fc; f++) {
1409
1410 double fbox[6] = { DBL_MAX, DBL_MAX, DBL_MAX,
1411 -DBL_MAX, -DBL_MAX, -DBL_MAX};
1412
1413 /* Reset vtag */
1414
1415 for (short int v = 0; v < cm->n_vc; v++) vtag[v] = -1;
1416
1417 /* Tag face vertices */
1418
1419 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1420 const short int *id = cm->e2v_ids + 2*cm->f2e_ids[i];
1421
1422 if (vtag[id[0]] < 0) {
1423 const double *xv = cm->xv + 3*id[0];
1424 if (xv[0] < fbox[0]) fbox[0] = xv[0];
1425 if (xv[1] < fbox[1]) fbox[1] = xv[1];
1426 if (xv[2] < fbox[2]) fbox[2] = xv[2];
1427 if (xv[0] > fbox[3]) fbox[3] = xv[0];
1428 if (xv[1] > fbox[4]) fbox[4] = xv[1];
1429 if (xv[2] > fbox[5]) fbox[5] = xv[2];
1430
1431 vtag[id[0]] = 1;
1432 }
1433
1434 if (vtag[id[1]] < 0) {
1435 const double *xv = cm->xv + 3*id[1];
1436 if (xv[0] < fbox[0]) fbox[0] = xv[0];
1437 if (xv[1] < fbox[1]) fbox[1] = xv[1];
1438 if (xv[2] < fbox[2]) fbox[2] = xv[2];
1439 if (xv[0] > fbox[3]) fbox[3] = xv[0];
1440 if (xv[1] > fbox[4]) fbox[4] = xv[1];
1441 if (xv[2] > fbox[5]) fbox[5] = xv[2];
1442
1443 vtag[id[1]] = 1;
1444 }
1445
1446 } /* Loop on face edges */
1447
1448 cm->f_diam[f] = cs_math_3_distance(fbox, fbox + 3);
1449
1450 } /* Loop on cell faces */
1451
1452 } /* Not a tetrahedron */
1453
1454 } /* Compute diameters */
1455 }
1456
1457 /*----------------------------------------------------------------------------*/
1458 /*!
1459 * \brief Allocate a cs_face_mesh_t structure
1460 *
1461 * \param[in] n_max_vbyf max. number of vertices fir a face
1462 *
1463 * \return a pointer to a new allocated cs_face_mesh_t structure
1464 */
1465 /*----------------------------------------------------------------------------*/
1466
1467 cs_face_mesh_t *
cs_face_mesh_create(short int n_max_vbyf)1468 cs_face_mesh_create(short int n_max_vbyf)
1469 {
1470 cs_face_mesh_t *fm = NULL;
1471
1472 BFT_MALLOC(fm, 1, cs_face_mesh_t);
1473
1474 fm->n_max_vbyf = n_max_vbyf;
1475
1476 fm->c_id = -1;
1477 fm->xc[0] = fm->xc[1] = fm->xc[2] = 0.;
1478
1479 /* Face-related quantities */
1480
1481 fm->f_id = -1;
1482 fm->f_sgn = 0;
1483 fm->pvol = 0.;
1484 fm->hfc = 0.;
1485
1486 /* Vertex-related quantities */
1487
1488 fm->n_vf = 0;
1489 BFT_MALLOC(fm->v_ids, fm->n_max_vbyf, cs_lnum_t);
1490 BFT_MALLOC(fm->xv, 3*fm->n_max_vbyf, double);
1491 BFT_MALLOC(fm->wvf, fm->n_max_vbyf, double);
1492
1493 /* Edge-related quantities */
1494
1495 fm->n_ef = 0;
1496 BFT_MALLOC(fm->e_ids, fm->n_max_vbyf, cs_lnum_t);
1497 BFT_MALLOC(fm->edge, fm->n_max_vbyf, cs_quant_t);
1498 BFT_MALLOC(fm->e2v_ids, 2*fm->n_max_vbyf, short int);
1499 BFT_MALLOC(fm->tef, fm->n_max_vbyf, double);
1500
1501 return fm;
1502 }
1503
1504 /*----------------------------------------------------------------------------*/
1505 /*!
1506 * \brief Get a pointer to a cs_face_mesh_t structure corresponding to mesh id
1507 *
1508 * \param[in] mesh_id id in the array of pointer to cs_face_mesh_t struct.
1509 *
1510 * \return a pointer to a cs_face_mesh_t structure
1511 */
1512 /*----------------------------------------------------------------------------*/
1513
1514 cs_face_mesh_t *
cs_cdo_local_get_face_mesh(int mesh_id)1515 cs_cdo_local_get_face_mesh(int mesh_id)
1516 {
1517 if (mesh_id < 0 || mesh_id >= cs_glob_n_threads)
1518 return NULL;
1519
1520 return cs_cdo_local_face_meshes[mesh_id];
1521 }
1522
1523 /*----------------------------------------------------------------------------*/
1524 /*!
1525 * \brief Free a cs_face_mesh_t structure
1526 *
1527 * \param[in, out] p_fm pointer of pointer to a cs_face_mesh_t structure
1528 */
1529 /*----------------------------------------------------------------------------*/
1530
1531 void
cs_face_mesh_free(cs_face_mesh_t ** p_fm)1532 cs_face_mesh_free(cs_face_mesh_t **p_fm)
1533 {
1534 cs_face_mesh_t *fm = *p_fm;
1535
1536 if (fm == NULL)
1537 return;
1538
1539 BFT_FREE(fm->v_ids);
1540 BFT_FREE(fm->xv);
1541 BFT_FREE(fm->wvf);
1542
1543 BFT_FREE(fm->e_ids);
1544 BFT_FREE(fm->edge);
1545 BFT_FREE(fm->e2v_ids);
1546 BFT_FREE(fm->tef);
1547
1548 BFT_FREE(fm);
1549 *p_fm = NULL;
1550 }
1551
1552 /*----------------------------------------------------------------------------*/
1553 /*!
1554 * \brief Define a cs_face_mesh_t structure for a given face/cell id.
1555 *
1556 * \param[in] c_id cell id
1557 * \param[in] f_id face id in the mesh structure
1558 * \param[in] connect pointer to a cs_cdo_connect_t structure
1559 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1560 * \param[in, out] fm pointer to a cs_face_mesh_t structure to set
1561 */
1562 /*----------------------------------------------------------------------------*/
1563
1564 void
cs_face_mesh_build(cs_lnum_t c_id,cs_lnum_t f_id,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_face_mesh_t * fm)1565 cs_face_mesh_build(cs_lnum_t c_id,
1566 cs_lnum_t f_id,
1567 const cs_cdo_connect_t *connect,
1568 const cs_cdo_quantities_t *quant,
1569 cs_face_mesh_t *fm)
1570 {
1571 if (fm == NULL)
1572 return;
1573
1574 assert(c_id > -1);
1575 assert(f_id > -1);
1576
1577 fm->c_id = c_id;
1578 const cs_real_t *xc = quant->cell_centers + 3*c_id;
1579 for (int k = 0; k < 3; k++) fm->xc[k] = xc[k];
1580
1581 /* Face-related quantities */
1582
1583 const cs_quant_t pfq = cs_quant_set_face(f_id, quant);
1584
1585 fm->f_id = f_id;
1586 fm->face.meas = pfq.meas;
1587 for (int k = 0; k < 3; k++) {
1588 fm->face.center[k] = pfq.center[k];
1589 fm->face.unitv[k] = pfq.unitv[k];
1590 }
1591
1592 const cs_lnum_t *c2f_idx = connect->c2f->idx + c_id;
1593 const cs_lnum_t *c2f_ids = connect->c2f->ids + c2f_idx[0];
1594 const int n_fc = c2f_idx[1] - c2f_idx[0];
1595
1596 assert(quant->pvol_fc != NULL);
1597
1598 short int _f = n_fc;
1599 for (short int f = 0; f < n_fc; f++) {
1600 if (c2f_ids[f] == f_id) {
1601
1602 const cs_lnum_t f_shift = c2f_idx[0]+f;
1603 const cs_nvec3_t de_nvect = cs_quant_set_dedge_nvec(f_shift, quant);
1604
1605 for (int k = 0; k < 3; k++) fm->dedge.unitv[k] = de_nvect.unitv[k];
1606 fm->dedge.meas = de_nvect.meas;
1607 fm->f_sgn = connect->c2f->sgn[f_shift];
1608 fm->pvol = quant->pvol_fc[f_shift];
1609 fm->hfc = 3*fm->pvol/fm->face.meas;
1610
1611 _f = f;
1612 break;
1613 }
1614 }
1615
1616 if (_f == n_fc) /* Sanity check */
1617 bft_error(__FILE__, __LINE__, 0,
1618 _(" Face %ld not found.\n Stop build a face mesh."), (long)f_id);
1619
1620 const cs_lnum_t *f2e_idx = connect->f2e->idx + f_id;
1621 const cs_lnum_t *f2e_lst = connect->f2e->ids + f2e_idx[0];
1622
1623 fm->n_vf = fm->n_ef = f2e_idx[1] - f2e_idx[0];
1624 short int nv = 0;
1625 for (int i = 0; i < fm->n_vf; i++)
1626 fm->v_ids[i] = -1;
1627
1628 for (short int e = 0; e < fm->n_ef; e++) {
1629
1630 const cs_lnum_t e_id = f2e_lst[e];
1631 const cs_nvec3_t e_nvect = cs_quant_set_edge_nvec(e_id, quant);
1632
1633 fm->e_ids[e] = e_id;
1634 fm->edge[e].meas = e_nvect.meas;
1635 for (int k = 0; k < 3; k++)
1636 fm->edge[e].unitv[k] = e_nvect.unitv[k];
1637 /* Still to handle the edge barycenter */
1638
1639 const cs_lnum_t *e2v_ids = connect->e2v->ids + 2*e_id;
1640 short int v1 = -1, v2 = -1;
1641 for (int v = 0; v < fm->n_vf && fm->v_ids[v] != -1; v++) {
1642 if (fm->v_ids[v] == e2v_ids[0])
1643 v1 = v;
1644 else if (fm->v_ids[v] == e2v_ids[1])
1645 v2 = v;
1646 }
1647
1648 /* Add vertices if not already identified */
1649
1650 if (v1 == -1) /* Not found -> Add v1 */
1651 fm->v_ids[nv] = e2v_ids[0], v1 = nv++;
1652 if (v2 == -1) /* Not found -> Add v2 */
1653 fm->v_ids[nv] = e2v_ids[1], v2 = nv++;
1654
1655 /* Update e2v_ids */
1656
1657 const int _eshft = 2*e;
1658 fm->e2v_ids[_eshft] = v1;
1659 fm->e2v_ids[_eshft+1] = v2;
1660
1661 } /* Loop on face edges */
1662
1663 assert(nv == fm->n_vf); /* Sanity check */
1664
1665 /* Update vertex coordinates */
1666
1667 int shift = 0;
1668 for (short int v = 0; v < fm->n_vf; v++) {
1669 const cs_real_t *xv = quant->vtx_coord + 3*fm->v_ids[v];
1670 for (int k = 0; k < 3; k++)
1671 fm->xv[shift++] = xv[k];
1672 }
1673
1674 /* Update the edge center. Define wvf and tef */
1675
1676 for (int i = 0; i < fm->n_vf; i++)
1677 fm->wvf[i] = 0;
1678
1679 for (short int e = 0; e < fm->n_ef; e++) {
1680
1681 const short int v1 = fm->e2v_ids[2*e];
1682 const short int v2 = fm->e2v_ids[2*e+1];
1683 const cs_real_t *xv1 = fm->xv + 3*v1;
1684 const cs_real_t *xv2 = fm->xv + 3*v2;
1685
1686 /* Update the edge center */
1687
1688 for (int k = 0; k < 3; k++)
1689 fm->edge[e].center[k] = 0.5 * (xv1[k] + xv2[k]);
1690
1691 /* tef = ||(xe -xf) x e||/2 = s(v1,e,f) + s(v2, e, f) */
1692
1693 const double tef = cs_compute_area_from_quant(fm->edge[e], pfq.center);
1694
1695 fm->wvf[v1] += tef;
1696 fm->wvf[v2] += tef;
1697 fm->tef[e] = tef;
1698
1699 } /* Loop on face edges */
1700
1701 const double invf = 0.5/pfq.meas;
1702 for (short int v = 0; v < fm->n_vf; v++) fm->wvf[v] *= invf;
1703 }
1704
1705 /*----------------------------------------------------------------------------*/
1706 /*!
1707 * \brief Define a cs_face_mesh_t structure for a given cell from a
1708 * cs_cell_mesh_t structure
1709 *
1710 * \param[in] cm pointer to the reference cs_cell_mesh_t structure
1711 * \param[in] f face id in the cs_cell_mesh_t structure
1712 * \param[in, out] fm pointer to a cs_face_mesh_t structure to set
1713 */
1714 /*----------------------------------------------------------------------------*/
1715
1716 void
cs_face_mesh_build_from_cell_mesh(const cs_cell_mesh_t * cm,short int f,cs_face_mesh_t * fm)1717 cs_face_mesh_build_from_cell_mesh(const cs_cell_mesh_t *cm,
1718 short int f,
1719 cs_face_mesh_t *fm)
1720 {
1721 if (fm == NULL || cm == NULL)
1722 return;
1723
1724 assert(f > -1 && f < cm->n_fc);
1725 assert(cs_eflag_test(cm->flag,
1726 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
1727 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV |
1728 CS_FLAG_COMP_HFQ));
1729
1730 fm->c_id = cm->c_id;
1731 for (int k = 0; k < 3; k++) fm->xc[k] = cm->xc[k];
1732
1733 /* Face-related quantities */
1734
1735 fm->f_id = f;
1736 fm->f_sgn = cm->f_sgn[f];
1737 fm->pvol = cm->pvol_f[f];
1738 fm->hfc = cm->hfc[f];
1739
1740 const cs_quant_t pfq = cm->face[f];
1741 fm->face.meas = pfq.meas;
1742 for (int k = 0; k < 3; k++) {
1743 fm->face.center[k] = pfq.center[k];
1744 fm->face.unitv[k] = pfq.unitv[k];
1745 }
1746
1747 const cs_nvec3_t deq = cm->dedge[f];
1748 fm->dedge.meas = deq.meas;
1749 for (int k = 0; k < 3; k++)
1750 fm->dedge.unitv[k] = deq.unitv[k];
1751
1752 const short int *f2e_idx = cm->f2e_idx + f;
1753 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
1754 const double *_tef = cm->tef + f2e_idx[0];
1755
1756 fm->n_vf = fm->n_ef = f2e_idx[1] - f2e_idx[0];
1757 short int nv = 0;
1758 for (int i = 0; i < fm->n_vf; i++) fm->v_ids[i] = -1;
1759
1760 for (short int ef = 0; ef < fm->n_ef; ef++) {
1761
1762 const short int ec = f2e_ids[ef];
1763
1764 fm->e_ids[ef] = ec;
1765 fm->tef[ef] = _tef[ef];
1766
1767 const cs_quant_t peq = cm->edge[ec];
1768 fm->edge[ef].meas = peq.meas;
1769 for (int k = 0; k < 3; k++) {
1770 fm->edge[ef].center[k] = peq.center[k];
1771 fm->edge[ef].unitv[k] = peq.unitv[k];
1772 }
1773
1774 const int eshft = 2*ec;
1775 short int v1c_id = cm->e2v_ids[eshft];
1776 short int v2c_id = cm->e2v_ids[eshft+1];
1777
1778 /* Compact vertex numbering to this face */
1779
1780 short int v1 = -1, v2 = -1;
1781 for (int v = 0; v < fm->n_vf && fm->v_ids[v] != -1; v++) {
1782 if (fm->v_ids[v] == v1c_id)
1783 v1 = v;
1784 else if (fm->v_ids[v] == v2c_id)
1785 v2 = v;
1786 }
1787
1788 /* Add vertices if not already identified */
1789
1790 if (v1 == -1) /* Not found -> Add v1 */
1791 fm->v_ids[nv] = v1c_id, v1 = nv++;
1792 if (v2 == -1) /* Not found -> Add v2 */
1793 fm->v_ids[nv] = v2c_id, v2 = nv++;
1794
1795 /* Update e2v_ids */
1796
1797 const int _eshft = 2*ef;
1798 fm->e2v_ids[_eshft] = v1;
1799 fm->e2v_ids[_eshft+1] = v2;
1800
1801 } /* Loop on face edges */
1802
1803 assert(nv == fm->n_vf); /* Sanity check */
1804
1805 /* Update vertex coordinates */
1806
1807 int shift = 0;
1808 for (short int v = 0; v < fm->n_vf; v++) {
1809 const cs_real_t *xv = cm->xv + 3*fm->v_ids[v];
1810 for (int k = 0; k < 3; k++)
1811 fm->xv[shift++] = xv[k];
1812 }
1813
1814 /* Define wvf and tef */
1815
1816 for (int i = 0; i < fm->n_vf; i++)
1817 fm->wvf[i] = 0;
1818
1819 for (short int e = 0; e < fm->n_ef; e++) {
1820
1821 const short int v1 = fm->e2v_ids[2*e];
1822 const short int v2 = fm->e2v_ids[2*e+1];
1823
1824 fm->wvf[v1] += _tef[e];
1825 fm->wvf[v2] += _tef[e];
1826
1827 }
1828
1829 const double invf = 0.5/pfq.meas;
1830 for (short int v = 0; v < fm->n_vf; v++) fm->wvf[v] *= invf;
1831 }
1832
1833 /*----------------------------------------------------------------------------*/
1834 /*!
1835 * \brief Allocate a cs_face_mesh_light_t structure
1836 *
1837 * \param[in] n_max_vbyf max. number of vertices for a face
1838 * \param[in] n_max_vbyc max. number of vertices for a cell
1839 *
1840 * \return a pointer to a new allocated cs_face_mesh_light_t structure
1841 */
1842 /*----------------------------------------------------------------------------*/
1843
1844 cs_face_mesh_light_t *
cs_face_mesh_light_create(short int n_max_vbyf,short int n_max_vbyc)1845 cs_face_mesh_light_create(short int n_max_vbyf,
1846 short int n_max_vbyc)
1847 {
1848 cs_face_mesh_light_t *fm = NULL;
1849
1850 BFT_MALLOC(fm, 1, cs_face_mesh_light_t);
1851
1852 fm->n_max_vbyf = n_max_vbyf;
1853 fm->c_id = -1;
1854 fm->f = -1;
1855
1856 /* Vertex-related quantities */
1857
1858 fm->n_vf = 0;
1859 BFT_MALLOC(fm->v_ids, n_max_vbyc, short int);
1860 BFT_MALLOC(fm->wvf, n_max_vbyc, double);
1861
1862 /* Edge-related quantities */
1863
1864 fm->n_ef = 0;
1865 BFT_MALLOC(fm->e_ids, fm->n_max_vbyf, short int);
1866 BFT_MALLOC(fm->tef, fm->n_max_vbyf, double);
1867
1868 return fm;
1869 }
1870
1871 /*----------------------------------------------------------------------------*/
1872 /*!
1873 * \brief Get a pointer to a cs_face_mesh_light_t structure corresponding to
1874 * mesh id
1875 *
1876 * \param[in] mesh_id id in the cs_face_mesh_light_t array
1877 *
1878 * \return a pointer to a cs_face_mesh_light_t structure
1879 */
1880 /*----------------------------------------------------------------------------*/
1881
1882 cs_face_mesh_light_t *
cs_cdo_local_get_face_mesh_light(int mesh_id)1883 cs_cdo_local_get_face_mesh_light(int mesh_id)
1884 {
1885 if (mesh_id < 0 || mesh_id >= cs_glob_n_threads)
1886 return NULL;
1887
1888 return cs_cdo_local_face_meshes_light[mesh_id];
1889 }
1890
1891 /*----------------------------------------------------------------------------*/
1892 /*!
1893 * \brief Free a cs_face_mesh_light_t structure
1894 *
1895 * \param[in, out] p_fm pointer of pointer to a cs_face_mesh_light_t struct.
1896 */
1897 /*----------------------------------------------------------------------------*/
1898
1899 void
cs_face_mesh_light_free(cs_face_mesh_light_t ** p_fm)1900 cs_face_mesh_light_free(cs_face_mesh_light_t **p_fm)
1901 {
1902 cs_face_mesh_light_t *fm = *p_fm;
1903
1904 if (fm == NULL)
1905 return;
1906
1907 BFT_FREE(fm->v_ids);
1908 BFT_FREE(fm->wvf);
1909 BFT_FREE(fm->e_ids);
1910 BFT_FREE(fm->tef);
1911
1912 BFT_FREE(fm);
1913 *p_fm = NULL;
1914 }
1915
1916 /*----------------------------------------------------------------------------*/
1917 /*!
1918 * \brief Define a cs_face_mesh_light_t structure starting from a
1919 * cs_cell_mesh_t structure.
1920 *
1921 * \param[in] cm pointer to the reference cs_cell_mesh_t structure
1922 * \param[in] f face id in the cs_cell_mesh_t structure
1923 * \param[in, out] fm pointer to a cs_face_mesh_light_t structure to set
1924 */
1925 /*----------------------------------------------------------------------------*/
1926
1927 void
cs_face_mesh_light_build(const cs_cell_mesh_t * cm,short int f,cs_face_mesh_light_t * fm)1928 cs_face_mesh_light_build(const cs_cell_mesh_t *cm,
1929 short int f,
1930 cs_face_mesh_light_t *fm)
1931 {
1932 if (fm == NULL || cm == NULL)
1933 return;
1934
1935 assert(f > -1 && f < cm->n_fc);
1936 assert(cs_eflag_test(cm->flag,
1937 CS_FLAG_COMP_PV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
1938
1939 fm->c_id = cm->c_id;
1940 fm->f = f;
1941
1942 const short int *f2e_idx = cm->f2e_idx + f;
1943 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
1944
1945 /* Initialization */
1946
1947 fm->n_vf = fm->n_ef = f2e_idx[1] - f2e_idx[0];
1948 for (int i = 0; i < cm->n_vc; i++) {
1949 fm->v_ids[i] = -1;
1950 fm->wvf[i] = 0;
1951 }
1952
1953 /* Define wvf from the knowledge of tef */
1954
1955 const double *_tef = cm->tef + f2e_idx[0];
1956
1957 for (short int e = 0; e < fm->n_ef; e++) {
1958
1959 const short int e_cellwise = f2e_ids[e];
1960 const short int *vc = cm->e2v_ids + 2*e_cellwise;
1961
1962 fm->e_ids[e] = e_cellwise;
1963 fm->tef[e] = _tef[e];
1964 fm->v_ids[vc[0]] = 1;
1965 fm->v_ids[vc[1]] = 1;
1966
1967 /* Build wvf */
1968
1969 fm->wvf[vc[0]] += _tef[e];
1970 fm->wvf[vc[1]] += _tef[e];
1971
1972 } /* Loop on face edges */
1973
1974 /* Compact vertex numbering to this face */
1975
1976 short int nv = 0; /* current vertex id in the face numbering */
1977 for (short int v = 0; v < cm->n_vc; v++) {
1978 if (fm->v_ids[v] > 0) {
1979 fm->v_ids[nv] = v;
1980 fm->wvf[nv] = fm->wvf[v];
1981 nv++;
1982 }
1983 }
1984
1985 assert(nv == fm->n_vf); /* Sanity check */
1986 const double invf = 0.5/cm->face[f].meas;
1987 for (short int v = 0; v < fm->n_vf; v++) fm->wvf[v] *= invf;
1988 }
1989
1990 /*----------------------------------------------------------------------------*/
1991
1992 #undef _dp3
1993
1994 END_C_DECLS
1995