1 /*============================================================================
2 * Manage geometrical quantities needed in CDO schemes
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 <limits.h>
34 #include <assert.h>
35 #include <float.h>
36 #include <math.h>
37 #include <string.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_mem.h>
44
45 #include "cs_array_reduce.h"
46 #include "cs_log.h"
47 #include "cs_order.h"
48 #include "cs_parall.h"
49 #include "cs_param_cdo.h"
50 #include "cs_prototypes.h"
51
52 /*----------------------------------------------------------------------------
53 * Header for the current file
54 *----------------------------------------------------------------------------*/
55
56 #include "cs_cdo_quantities.h"
57
58 /*----------------------------------------------------------------------------*/
59
60 BEGIN_C_DECLS
61
62 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
63
64 /*=============================================================================
65 * Local Macro definitions
66 *============================================================================*/
67
68 #define CS_CDO_QUANTITIES_DBG 0 /* Switch off/on debug information */
69
70 /* Redefined names of function from cs_math to get shorter names */
71 #define _n3 cs_math_3_norm
72 #define _dp3 cs_math_3_dot_product
73
74 /*=============================================================================
75 * Local structure definitions
76 *============================================================================*/
77
78 /* Temporary structures to build mesh quantities */
79
80 typedef struct {
81
82 int XYZ[3]; /* Direct permutation of the ref. axis such that nZ
83 is maximal */
84 cs_nvec3_t q; /* face surface and its unit normal */
85 double omega; /* P = Point belonging to the face omega = - < n, P> */
86
87 } _cdo_fspec_t;
88
89 typedef struct { /* Face sub-quantities */
90
91 double F1;
92 double Fa;
93 double Fb;
94 double Fc;
95 double Fa2;
96 double Fb2;
97 double Fc2;
98
99 } _cdo_fsubq_t;
100
101 typedef struct { /* These quantities are the integral of q on the plane
102 (alpha, beta) where the face is projected */
103
104 double p1; /* q = 1 */
105 double pa; /* q = alpha */
106 double pb; /* q = beta */
107 double pc; /* q = gamma */
108 double pab; /* q = alpha * beta */
109 double pa2; /* q = alpha^2 */
110 double pb2; /* q = beta^2 */
111
112 } _cdo_projq_t;
113
114 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
115
116 /*============================================================================
117 * Static global variables
118 *============================================================================*/
119
120 /* Store in a flag which quantities have to be computed */
121 cs_flag_t cs_cdo_quantities_flag = 0;
122
123 /* Algorithm used to compute the cell center */
124 cs_cdo_quantities_cell_center_algo_t
125 cs_cdo_cell_center_algo = CS_CDO_QUANTITIES_BARYC_CENTER;
126
127 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
128
129 /*============================================================================
130 * Private function prototypes
131 *============================================================================*/
132
133 /*----------------------------------------------------------------------------*/
134 /*!
135 * \brief Allocate and initialize a \ref cs_mesh_quantities_t structure
136 *
137 * \return a pointer to the newly allocated structure
138 */
139 /*----------------------------------------------------------------------------*/
140
141 static cs_cdo_quantities_t *
_create_cdo_quantities(void)142 _create_cdo_quantities(void)
143 {
144 cs_cdo_quantities_t *cdoq = NULL;
145
146 /* Build cs_cdo_quantities_t structure */
147 BFT_MALLOC(cdoq, 1, cs_cdo_quantities_t);
148
149 cdoq->remove_boundary_faces = false;
150
151 cdoq->vol_tot = 0.;
152
153 /* Cell-based quantities */
154 cdoq->cell_info.h_min = cdoq->cell_info.meas_min = DBL_MAX;
155 cdoq->cell_info.h_max = cdoq->cell_info.meas_max = -DBL_MAX;
156
157 cdoq->n_cells = 0;
158 cdoq->n_g_cells = 0;
159 cdoq->cell_centers = NULL;
160 cdoq->cell_vol = NULL;
161
162 /* Face-based quantities */
163 cdoq->face_info.h_min = cdoq->face_info.meas_min = DBL_MAX;
164 cdoq->face_info.h_max = cdoq->face_info.meas_max = -DBL_MAX;
165
166 cdoq->n_faces = cdoq->n_i_faces = cdoq->n_b_faces = 0;
167 cdoq->n_g_faces = 0;
168 cdoq->dedge_vector = NULL;
169 cdoq->pvol_fc = NULL;
170
171 /* Edge-based quantities */
172 cdoq->edge_info.h_min = cdoq->edge_info.meas_min = DBL_MAX;
173 cdoq->edge_info.h_max = cdoq->edge_info.meas_max = -DBL_MAX;
174
175 cdoq->n_edges = 0;
176 cdoq->n_g_edges = 0;
177 cdoq->edge_vector = NULL;
178 cdoq->pvol_ec = NULL;
179 cdoq->dface_normal = NULL;
180
181 /* Vertex-based quantities */
182 cdoq->n_vertices = 0;
183 cdoq->n_g_vertices = 0;
184 cdoq->dcell_vol = NULL;
185
186 /* Shared pointers are not initialized at this stage */
187
188 return cdoq;
189 }
190
191 /*----------------------------------------------------------------------------*/
192 /*!
193 * \brief Function related to the Mirtich algorithm
194 * - Define an unitary normal to the current face
195 * - Compute omega = - <n, P> where P belongs to the face
196 * - Choose projection axis in order to maximize the projected area
197 * - Define a direct basis (alpha, beta, gamma) with this choice
198 *
199 * \param[in] f_id id of the face to treat
200 * \param[in] m pointer to a cs_mesh_t structure
201 * \param[in] mq pointer to a cs_mesh_quantities_t structure
202 *
203 * \return a _cdo_fspec_t structure storing the computed quantities
204 */
205 /*----------------------------------------------------------------------------*/
206
207 static _cdo_fspec_t
_get_fspec(cs_lnum_t f_id,const cs_mesh_t * m,const cs_mesh_quantities_t * mq)208 _get_fspec(cs_lnum_t f_id,
209 const cs_mesh_t *m,
210 const cs_mesh_quantities_t *mq)
211 {
212 const int X = 0, Y = 1, Z = 2;
213
214 cs_lnum_t f, j, k, v, e, s;
215 double inv_n, nx, ny, nz;
216 double P[3]; /* Point belonging to the current face */
217
218 _cdo_fspec_t fspec =
219 { .XYZ = {X, Y, Z},
220 .omega = 0.,
221 .q.meas = 0.,
222 .q.unitv[0] = 0, .q.unitv[1] = 0, .q.unitv[2] = 0,
223 };
224
225 /* Treatment according to the kind of face (interior or border) */
226
227 if (f_id < m->n_i_faces) { /* Interior face */
228
229 /* Choose a vertex belonging to this face */
230 f = f_id;
231 s = m->i_face_vtx_idx[f];
232 e = m->i_face_vtx_idx[f+1];
233 inv_n = 1.0 / (e - s);
234
235 for (k = 0; k < 3; k++)
236 P[k] = 0.0;
237
238 for (j = s; j < e; j++) {
239 v = m->i_face_vtx_lst[j];
240 for (k = 0; k < 3; k++)
241 P[k] += m->vtx_coord[3*v+k];
242 }
243
244 for (k = 0; k < 3; k++)
245 P[k] *= inv_n;
246
247 cs_nvec3(&(mq->i_face_normal[3*f]), &(fspec.q));
248
249 }
250 else { /* Border face */
251
252 /* Choose a vertex belonging to this face */
253 f = f_id - m->n_i_faces;
254 s = m->b_face_vtx_idx[f];
255 e = m->b_face_vtx_idx[f+1];
256 inv_n = 1.0 / (e - s);
257
258 for (k = 0; k < 3; k++)
259 P[k] = 0.0;
260
261 for (j = s; j < e; j++) {
262 v = m->b_face_vtx_lst[j];
263 for (k = 0; k < 3; k++)
264 P[k] += m->vtx_coord[3*v+k];
265 }
266
267 for (k = 0; k < 3; k++)
268 P[k] *= inv_n;
269
270 cs_nvec3(&(mq->b_face_normal[3*f]), &(fspec.q));
271
272 }
273
274 /* Define omega = -<n,P>*/
275 fspec.omega = - _dp3(P, fspec.q.unitv);
276
277 /* Define a direct basis such that n[Z] maximal */
278 nx = fabs(fspec.q.unitv[X]);
279 ny = fabs(fspec.q.unitv[Y]);
280 nz = fabs(fspec.q.unitv[Z]);
281
282 if (nx > ny && nx > nz)
283 fspec.XYZ[Z] = X;
284 else
285 fspec.XYZ[Z] = (ny > nz) ? Y : Z;
286
287 fspec.XYZ[X] = (fspec.XYZ[Z] + 1) % 3;
288 fspec.XYZ[Y] = (fspec.XYZ[X] + 1) % 3;
289
290 #if CS_CDO_QUANTITIES_DBG > 1 && defined(DEBUG) && !defined(NDEBUG)
291 printf("\n F: (%d, %d) >> surf: %e; omega: %e; XYZ: %d%d%d; [%e, %e, %e]\n",
292 f_id, f, fspec.q.meas, fspec.omega,
293 fspec.XYZ[0], fspec.XYZ[1], fspec.XYZ[2],
294 fspec.q.unitv[0], fspec.q.unitv[1], fspec.q.unitv[2]);
295 #endif
296
297 return fspec;
298 }
299
300 /*----------------------------------------------------------------------------*/
301 /*!
302 * \brief Function related to the Mirtich algorithm
303 * Compute projected integrals and quantities
304 *
305 * \param[in] f_id id of the face to treat
306 * \param[in] connect pointer to a cs_cdo_connect_t structure
307 * \param[in] coords coordinates of the mesh vertices
308 * \param[in] axis local basis to the given face
309 *
310 * \return a _cdo_projq_t structure storing the computed quantities
311 */
312 /*----------------------------------------------------------------------------*/
313
314 static _cdo_projq_t
_get_proj_quantities(cs_lnum_t f_id,const cs_cdo_connect_t * connect,const cs_real_t coords[],const int axis[])315 _get_proj_quantities(cs_lnum_t f_id,
316 const cs_cdo_connect_t *connect,
317 const cs_real_t coords[],
318 const int axis[])
319 {
320 cs_lnum_t v_id[2];
321 cs_real_t a[2], b[2], a2[2], b2[2];
322
323 const cs_adjacency_t *f2e = connect->f2e;
324 const cs_adjacency_t *e2v = connect->e2v;
325
326 /* Initialize structure */
327 _cdo_projq_t projq;
328
329 /* These quantities are the integral of q on the plane
330 (alpha, beta) where the face is projected */
331
332 projq.p1 = projq.pa = projq.pb = projq.pc = 0.0;
333 projq.pab = projq.pa2 = projq.pb2 = 0.0;
334
335 /* Scan edges which belong to the current face */
336 for (cs_lnum_t i = f2e->idx[f_id]; i < f2e->idx[f_id+1]; i++) {
337
338 short int e_sgn = f2e->sgn[i];
339 cs_lnum_t e_id = f2e->ids[i];
340 cs_lnum_t s = 2*e_id;
341
342 if (e_sgn > 0)
343 v_id[0] = e2v->ids[s], v_id[1] = e2v->ids[s+1];
344 else
345 v_id[0] = e2v->ids[s+1], v_id[1] = e2v->ids[s];
346
347 /* Vertices in the plane (alpha, beta) */
348 for (int k = 0; k < 2; k++) {
349 a[k] = coords[3*v_id[k] + axis[0]];
350 a2[k] = a[k]*a[k];
351 b[k] = coords[3*v_id[k] + axis[1]];
352 b2[k] = b[k]*b[k];
353 }
354
355 /* Related variables */
356 cs_real_t a0_3 = a2[0] * a[0];
357 cs_real_t b0_3 = b2[0] * b[0];
358 cs_real_t da = a[1] - a[0], db = b[1] - b[0];
359 cs_real_t C1 = a[0] + a[1];
360 cs_real_t Ca = C1 * a[1] + a2[0];
361 cs_real_t Cb = b2[1] + b[1]*b[0] + b2[0];
362 cs_real_t Ca2 = a[1] * Ca + a0_3;
363 cs_real_t Cb2 = b[1] * Cb + b0_3;
364 cs_real_t Cab = 3*a2[1] + 2*a[1]*a[0] + a2[0];
365 cs_real_t Kab = a2[1] + 2*a[1]*a[0] + 3*a2[0];
366
367 projq.p1 += db * C1;
368 projq.pa += db * Ca;
369 projq.pb += da * Cb;
370 projq.pa2 += db * Ca2;
371 projq.pb2 += da * Cb2;
372 projq.pab += db * (b[1] * Cab + b[0] * Kab);
373
374 } /* Loop on face edges */
375
376 projq.p1 *= 0.5;
377 projq.pa *= cs_math_1ov6;
378 projq.pb *= -cs_math_1ov6;
379 projq.pab *= cs_math_1ov24;
380 projq.pa2 *= cs_math_1ov12;
381 projq.pb2 *= -cs_math_1ov12;
382
383 return projq;
384 }
385
386 /*----------------------------------------------------------------------------*/
387 /*!
388 * \brief Function related to the Mirtich algorithm
389 * Compute quantities related to the sub-faces of a given face
390 *
391 * \param[in] f_id id of the face to treat
392 * \param[in] connect pointer to a cs_cdo_connect_t structure
393 * \param[in] coord coordinates of the mesh vertices
394 *
395 *
396 * \return a _cdo_fsubq_t structure storing the computed quantities
397 */
398 /*----------------------------------------------------------------------------*/
399
400 static _cdo_fsubq_t
_get_fsub_quantities(cs_lnum_t f_id,const cs_cdo_connect_t * connect,const cs_real_t * coord,const _cdo_fspec_t fspec)401 _get_fsub_quantities(cs_lnum_t f_id,
402 const cs_cdo_connect_t *connect,
403 const cs_real_t *coord,
404 const _cdo_fspec_t fspec)
405 {
406 const double na = fspec.q.unitv[fspec.XYZ[0]];
407 const double nb = fspec.q.unitv[fspec.XYZ[1]];
408 const double nc = fspec.q.unitv[fspec.XYZ[2]];
409 const double k1 = 1./nc, k2 = k1 * k1, k3 = k2 * k1;
410
411 /* Compute projected quantities */
412 const _cdo_projq_t projq = _get_proj_quantities(f_id, connect, coord, fspec.XYZ);
413
414 #if CS_CDO_QUANTITIES_DBG > 1 && defined(DEBUG) && !defined(NDEBUG)
415 printf(" F: %d >> p1: %.4e, pa: %.4e, pb: %.4e, pc: %.4e,"
416 " pab: %.4e, pa2: %.4e, pb2: %.4e\n",
417 f_id, projq.p1, projq.pa, projq.pb, projq.pc,
418 projq.pab, projq.pa2, projq.pb2);
419 #endif
420
421 _cdo_fsubq_t fsubq;
422
423 /* Compute face sub-quantities */
424 fsubq.F1 = k1*projq.p1;
425 fsubq.Fa = k1 * projq.pa;
426 fsubq.Fb = k1 * projq.pb;
427 fsubq.Fc = -k2 * (projq.pa * na + projq.pb * nb + fspec.omega * projq.p1);
428
429 fsubq.Fa2 = k1 * projq.pa2;
430 fsubq.Fb2 = k1 * projq.pb2;
431 fsubq.Fc2 = k3 * (na*na * projq.pa2 + 2*na*nb * projq.pab +
432 nb*nb * projq.pb2 +
433 fspec.omega * (2*na * projq.pa +
434 2*nb * projq.pb +
435 fspec.omega * projq.p1));
436
437 return fsubq;
438 }
439
440 /*----------------------------------------------------------------------------*/
441 /*!
442 * \brief Compute additional quantities related to faces like dual edges
443 * (segment between x_f and x_c and scanned with c2f adjacency)
444 *
445 * \param[in] topo pointer to a cs_cdo_connect_t structure
446 * \param[in, out] cdoq pointer to cs_cdo_quantities_t structure
447 */
448 /*----------------------------------------------------------------------------*/
449
450 static void
_compute_face_based_quantities(const cs_cdo_connect_t * topo,cs_cdo_quantities_t * cdoq)451 _compute_face_based_quantities(const cs_cdo_connect_t *topo,
452 cs_cdo_quantities_t *cdoq)
453 {
454 /* Compute dual edge quantities */
455 const cs_lnum_t n_cells = cdoq->n_cells;
456 const cs_adjacency_t *c2f = topo->c2f;
457
458 BFT_MALLOC(cdoq->dedge_vector, 3*c2f->idx[n_cells], cs_real_t);
459
460 # pragma omp parallel for if (n_cells > CS_THR_MIN)
461 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
462
463 const cs_real_t *xc = cdoq->cell_centers + 3*c_id;
464
465 for (cs_lnum_t i = c2f->idx[c_id]; i < c2f->idx[c_id+1]; i++) {
466
467 const cs_lnum_t f_id = c2f->ids[i];
468 const short int sgn = c2f->sgn[i];
469 const cs_real_t *xf = (f_id < cdoq->n_i_faces) ?
470 cdoq->i_face_center + 3*f_id :
471 cdoq->b_face_center + 3*(f_id - cdoq->n_i_faces);
472
473 for (int k = 0; k < 3; k++)
474 cdoq->dedge_vector[3*i+k] = sgn * (xf[k] - xc[k]);
475
476 } /* Loop on cell faces */
477
478 } /* End of loop on cells */
479
480 /* Compute the volume of the pyramid */
481 cs_flag_t masks[3] = { CS_CDO_QUANTITIES_FB_SCHEME,
482 CS_CDO_QUANTITIES_HHO_SCHEME,
483 CS_CDO_QUANTITIES_VCB_SCHEME };
484 if (cs_flag_at_least(cs_cdo_quantities_flag, 3, masks))
485 cs_cdo_quantities_compute_pvol_fc(cdoq, topo->c2f, &(cdoq->pvol_fc));
486 }
487
488 /*----------------------------------------------------------------------------*/
489 /*!
490 * \brief Compute dual face normals (face crossed by primal edges).
491 * Given a cell c and an edge e, there are two faces attached to the
492 * couple (c, e)
493 * For a face f in c, the triplet (e, f, c) induces a triangle s(e,f,c)
494 * The dual face is the union of these two triangles.
495 * Storage is based on the c2e connectivity
496 *
497 * \param[in] topo pointer to a cs_cdo_connect_t structure
498 * \param[in, out] quant pointer to cs_cdo_quantities_t structure
499 */
500 /*----------------------------------------------------------------------------*/
501
502 static void
_compute_edge_based_quantities(const cs_cdo_connect_t * topo,cs_cdo_quantities_t * quant)503 _compute_edge_based_quantities(const cs_cdo_connect_t *topo,
504 cs_cdo_quantities_t *quant)
505 {
506 /* Sanity check */
507 assert(topo->e2v != NULL);
508 assert(topo->f2e != NULL && topo->f2c != NULL && topo->c2e != NULL);
509
510 const cs_lnum_t n_cells = quant->n_cells;
511 const cs_lnum_t n_edges = quant->n_edges;
512
513 cs_real_t *edge_center = NULL;
514
515 /* Build edge centers and edge vectors */
516 BFT_MALLOC(edge_center, 3*n_edges, cs_real_t);
517 BFT_MALLOC(quant->edge_vector, 3*n_edges, cs_real_t);
518
519 # pragma omp parallel for if (n_edges > CS_THR_MIN)
520 for (cs_lnum_t e_id = 0; e_id < n_edges; e_id++) {
521
522 /* Get the two vertex ids related to the current edge */
523 const cs_lnum_t *v_ids = topo->e2v->ids + 2*e_id;
524 const cs_real_t *xa = quant->vtx_coord + 3*v_ids[0];
525 const cs_real_t *xb = quant->vtx_coord + 3*v_ids[1];
526
527 cs_real_t *xe = edge_center + 3*e_id;
528 cs_real_t *vect = quant->edge_vector + 3*e_id;
529 if (v_ids[1] > v_ids[0]) { /* vb > va */
530
531 for (int k = 0; k < 3; k++) {
532 vect[k] = xb[k] - xa[k];
533 xe[k] = 0.5* (xb[k] + xa[k]) ;
534 }
535
536 }
537 else {
538
539 for (int k = 0; k < 3; k++) {
540 vect[k] = xa[k] - xb[k];
541 xe[k] = 0.5* (xb[k] + xa[k]) ;
542 }
543
544 }
545
546 } /* End of loop on edges */
547
548 cs_flag_t masks[3] = { CS_CDO_QUANTITIES_EB_SCHEME,
549 CS_CDO_QUANTITIES_VB_SCHEME,
550 CS_CDO_QUANTITIES_VCB_SCHEME };
551 if (!cs_flag_at_least(cs_cdo_quantities_flag, 3, masks)) {
552 BFT_FREE(edge_center);
553 return;
554 }
555
556 /* Allocate and initialize array
557 * a) Compute the two vector areas composing each dual face
558 * b) Compute the volume associated to each edge in a cell
559 */
560 BFT_MALLOC(quant->pvol_ec, topo->c2e->idx[n_cells], cs_real_t);
561 BFT_MALLOC(quant->dface_normal, 3*topo->c2e->idx[n_cells], cs_real_t);
562
563 # pragma omp parallel shared(quant, topo, edge_center)
564 { /* OMP Block */
565
566 const cs_adjacency_t *c2f = topo->c2f, *f2e = topo->f2e;
567
568 # pragma omp for CS_CDO_OMP_SCHEDULE
569 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
570
571 const cs_lnum_t *c2e_idx = topo->c2e->idx + c_id;
572 const cs_lnum_t *c2e_ids = topo->c2e->ids + c2e_idx[0];
573 const short int n_ec = c2e_idx[1] - c2e_idx[0];
574
575 /* Initialize cell_dface */
576 cs_real_t *cell_dface = quant->dface_normal + 3*c2e_idx[0];
577 memset(cell_dface, 0, 3*n_ec*sizeof(cs_real_t));
578
579 /* Get the cell center */
580 const cs_real_t *xc = quant->cell_centers + 3*c_id;
581
582 for (cs_lnum_t i = c2f->idx[c_id]; i < c2f->idx[c_id+1]; i++) {
583
584 const cs_lnum_t f_id = c2f->ids[i];
585
586 /* Compute xf -> xc */
587 const cs_real_t *xf = (f_id < quant->n_i_faces) ?
588 quant->i_face_center + 3* f_id :
589 quant->b_face_center + 3*(f_id - quant->n_i_faces);
590 const cs_real_3_t xfxc = { xf[0] - xc[0],
591 xf[1] - xc[1],
592 xf[2] - xc[2] };
593
594 for (cs_lnum_t j = f2e->idx[f_id]; j < f2e->idx[f_id+1]; j++) {
595
596 const cs_lnum_t e_id = topo->f2e->ids[j];
597 const cs_real_t *xe = edge_center + 3*e_id;
598
599 /* Compute the vectorial area for the triangle : xc, xf, xe */
600 cs_real_3_t tria_vect, xexc;
601 for (int k = 0; k < 3; k++)
602 xexc[k] = xc[k] - xe[k];
603 cs_math_3_cross_product(xfxc, xexc, tria_vect);
604
605 /* Find the corresponding local id for this cell edge */
606 short int e = n_ec;
607 for (short int _e = 0; _e < n_ec; _e++) {
608 if (c2e_ids[_e] == e_id) {
609 e = _e;
610 break;
611 }
612 }
613 CS_CDO_OMP_ASSERT(e < n_ec);
614
615 /* One should have (tria.unitv, edge.unitv) > 0 */
616 cs_nvec3_t edge = cs_quant_set_edge_nvec(e_id, quant);
617 cs_nvec3_t tria;
618 cs_nvec3(tria_vect, &tria);
619 const double orient = _dp3(tria.unitv, edge.unitv);
620 CS_CDO_OMP_ASSERT(fabs(orient) > 0);
621
622 /* Store this portion of dual face area at the right place */
623 cs_real_t *_dface = cell_dface + 3*e;
624 if (orient < 0)
625 for (int k = 0; k < 3; k++) _dface[k] -= 0.5 * tria_vect[k];
626 else
627 for (int k = 0; k < 3; k++) _dface[k] += 0.5 * tria_vect[k];
628
629 } /* Loop on face edges */
630
631 } /* Loop on cell faces */
632
633 /* Compute pvol_ec */
634 cs_real_t *_pvol = quant->pvol_ec + c2e_idx[0];
635 for (short int e = 0; e < n_ec; e++) {
636 _pvol[e] = cs_math_1ov3 * _dp3(cell_dface + 3*e,
637 quant->edge_vector + 3*c2e_ids[e]);
638 CS_CDO_OMP_ASSERT(_pvol[e] > 0);
639 }
640
641 } /* End of loop on cells */
642
643 } /* End of OpenMP block */
644
645 BFT_FREE(edge_center);
646 }
647
648 /*----------------------------------------------------------------------------*/
649 /*!
650 * \brief Compute dual cell volumes related to primal vertices
651 * Storage is based on the c2v connectivity
652 *
653 * \param[in] topo pointer to a cs_cdo_connect_t structure
654 * \param[in, out] quant pointer to cs_cdo_quantities_t structure
655 */
656 /*----------------------------------------------------------------------------*/
657
658 static void
_compute_dcell_quantities(const cs_cdo_connect_t * topo,cs_cdo_quantities_t * quant)659 _compute_dcell_quantities(const cs_cdo_connect_t *topo,
660 cs_cdo_quantities_t *quant)
661 {
662 cs_flag_t masks[2] = { CS_CDO_QUANTITIES_VB_SCHEME |
663 CS_CDO_QUANTITIES_VCB_SCHEME };
664 if (!cs_flag_at_least(cs_cdo_quantities_flag, 2, masks))
665 return;
666
667 /* Sanity checks */
668 assert(topo->f2e != NULL && topo->f2c != NULL && topo->c2v != NULL);
669
670 const cs_adjacency_t *c2f = topo->c2f, *f2e = topo->f2e;
671
672 /* Allocate and initialize arrays */
673 BFT_MALLOC(quant->dcell_vol, topo->c2v->idx[quant->n_cells], double);
674
675 # pragma omp parallel for shared(quant, topo, c2f, f2e) \
676 CS_CDO_OMP_SCHEDULE
677 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
678
679 /* Compute part of dual volume related to each primal cell */
680 const cs_lnum_t *c2v_idx = topo->c2v->idx + c_id;
681 const cs_lnum_t *c2v_ids = topo->c2v->ids + c2v_idx[0];
682 const short int n_vc = c2v_idx[1] - c2v_idx[0];
683 const cs_real_t *xc = quant->cell_centers + 3*c_id;
684
685 /* Initialize */
686 double *vol_vc = quant->dcell_vol + c2v_idx[0];
687 for (short int v = 0; v < n_vc; v++) vol_vc[v] = 0.0;
688
689 for (cs_lnum_t jf = c2f->idx[c_id]; jf < c2f->idx[c_id+1]; jf++) {
690
691 const cs_lnum_t f_id = topo->c2f->ids[jf];
692 const cs_lnum_t bf_id = f_id - quant->n_i_faces;
693
694 const cs_real_t *xf;
695 if (bf_id > -1)
696 xf = quant->b_face_center + 3*bf_id;
697 else
698 xf = quant->i_face_center + 3*f_id;
699
700 for (cs_lnum_t je = f2e->idx[f_id]; je < f2e->idx[f_id+1]; je++) {
701
702 const cs_lnum_t e_id = f2e->ids[je];
703 const cs_lnum_t *v_ids = topo->e2v->ids + 2*e_id;
704 const cs_lnum_t v1_id = v_ids[0], v2_id = v_ids[1];
705 const double pvol = 0.5 * cs_math_voltet(quant->vtx_coord + 3*v1_id,
706 quant->vtx_coord + 3*v2_id,
707 xf,
708 xc);
709
710 /* Find the corresponding local cell edge */
711 short int _v1 = n_vc, _v2 = n_vc;
712 for (short int _v = 0; _v < n_vc; _v++) {
713 if (c2v_ids[_v] == v1_id) _v1 = _v;
714 if (c2v_ids[_v] == v2_id) _v2 = _v;
715 }
716 CS_CDO_OMP_ASSERT(_v1 < n_vc && _v2 < n_vc);
717 vol_vc[_v1] += pvol;
718 vol_vc[_v2] += pvol;
719
720 } /* Loop on face edges */
721
722 } /* Loop on cell faces */
723
724 } /* Loop on cells */
725
726 }
727
728
729 /*----------------------------------------------------------------------------*/
730 /*!
731 * \brief Define the cs_quant_info_t structures related to cells, faces and
732 * edges. For edges, it depends on scheme flags
733 *
734 * \param[in, out] quant pointer to a cs_cdo_quantities_t structure
735 */
736 /*----------------------------------------------------------------------------*/
737
738 static void
_compute_quant_info(cs_cdo_quantities_t * quant)739 _compute_quant_info(cs_cdo_quantities_t *quant)
740 {
741 if (quant == NULL)
742 return;
743
744 /* Cell info */
745 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
746
747 const double meas = quant->cell_vol[c_id];
748
749 if (meas > quant->cell_info.meas_max) {
750 quant->cell_info.meas_max = meas;
751 quant->cell_info.h_max = cbrt(meas);
752 }
753 if (meas < quant->cell_info.meas_min) {
754 quant->cell_info.meas_min = meas;
755 quant->cell_info.h_min = cbrt(meas);
756 }
757
758 } /* Loop on cells */
759
760 /* Face info */
761 for (cs_lnum_t f_id = 0; f_id < quant->n_i_faces; f_id++) {
762
763 const cs_real_t meas = quant->i_face_surf[f_id];
764
765 if (meas > quant->face_info.meas_max) {
766 quant->face_info.meas_max = meas;
767 quant->face_info.h_max = sqrt(meas);
768 }
769 if (meas < quant->face_info.meas_min) {
770 quant->face_info.meas_min = meas;
771 quant->face_info.h_min = sqrt(meas);
772 }
773
774 } /* Loop on interior faces */
775
776 for (cs_lnum_t f_id = 0; f_id < quant->n_b_faces; f_id++) {
777
778 const cs_real_t meas = quant->b_face_surf[f_id];
779
780 if (meas > quant->face_info.meas_max) {
781 quant->face_info.meas_max = meas;
782 quant->face_info.h_max = sqrt(meas);
783 }
784 if (meas < quant->face_info.meas_min) {
785 quant->face_info.meas_min = meas;
786 quant->face_info.h_min = sqrt(meas);
787 }
788
789 } /* Loop on border faces */
790
791 /* Edge info */
792 if (quant->edge_vector != NULL) {
793
794 for (cs_lnum_t e_id = 0; e_id < quant->n_edges; e_id++) {
795
796 cs_nvec3_t edge = cs_quant_set_edge_nvec(e_id, quant);
797
798 if (edge.meas > quant->edge_info.meas_max) {
799 quant->edge_info.meas_max = edge.meas;
800 quant->edge_info.h_max = edge.meas;
801 }
802 if (edge.meas < quant->edge_info.meas_min) {
803 quant->edge_info.meas_min = edge.meas;
804 quant->edge_info.h_min = edge.meas;
805 }
806
807 } /* Loop on edges */
808
809 }
810
811 if (cs_glob_n_ranks > 1) { /* Synchronization across ranks */
812
813 double buf[12] =
814 { quant->cell_info.meas_max, quant->cell_info.h_max,
815 -quant->cell_info.meas_min, -quant->cell_info.h_min,
816 quant->face_info.meas_max, quant->face_info.h_max,
817 -quant->face_info.meas_min, -quant->face_info.h_min,
818 quant->edge_info.meas_max, quant->edge_info.h_max,
819 -quant->edge_info.meas_min, -quant->edge_info.h_min };
820
821 cs_parall_max(12, CS_DOUBLE, buf);
822
823 quant->cell_info.meas_max = buf[0];
824 quant->cell_info.h_max = buf[1];
825 quant->cell_info.meas_min = -buf[2];
826 quant->cell_info.h_min = -buf[3];
827 quant->face_info.meas_max = buf[4];
828 quant->face_info.h_max = buf[5];
829 quant->face_info.meas_min = -buf[6];
830 quant->face_info.h_min = -buf[7];
831 quant->edge_info.meas_max = buf[8];
832 quant->edge_info.h_max = buf[9];
833 quant->edge_info.meas_min = -buf[10];
834 quant->edge_info.h_min = -buf[11];
835
836 } /* Synchronization */
837 }
838
839 /*----------------------------------------------------------------------------
840 Algorithm for computing mesh quantities : cell centers are computed as the
841 vertex average over cell vertices. Other quantities are copied from those
842 computed by Code_Saturne (default algorithm)
843 ----------------------------------------------------------------------------*/
844
845 static void
_vtx_algorithm(const cs_cdo_connect_t * connect,cs_cdo_quantities_t * quant)846 _vtx_algorithm(const cs_cdo_connect_t *connect,
847 cs_cdo_quantities_t *quant) /* In/out */
848 {
849 const cs_lnum_t n_cells = quant->n_cells;
850 const cs_adjacency_t *c2v = connect->c2v;
851
852 /* Compute cell centers */
853 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
854
855 const cs_lnum_t vs = c2v->idx[c_id];
856 const cs_lnum_t ve = c2v->idx[c_id+1];
857
858 /* Sanity checks */
859 assert(ve - vs > 0);
860 assert(ve - vs < SHRT_MAX);
861
862 double *xc = quant->cell_centers + 3*c_id;
863
864 xc[0] = xc[1] = xc[2] = 0;
865 for (cs_lnum_t jv = vs; jv < ve; jv++) {
866
867 const cs_real_t *xv = quant->vtx_coord + 3*c2v->ids[jv];
868
869 xc[0] += xv[0];
870 xc[1] += xv[1];
871 xc[2] += xv[2];
872
873 } /* Loop on cell vertices */
874
875 const double coef = 1./(ve-vs);
876 for (int k = 0; k < 3; k++) xc[k] *= coef;
877
878 } /* Loop on cells */
879
880 }
881
882 /*----------------------------------------------------------------------------*/
883 /*!
884 * \brief Algorithm for computing the real cell barycenters inspired from the
885 * article "Fast and accurate computation of polyhedral mass properties"
886 * Journal of Graphics, 1997 by Brian Mirtich
887 *
888 * \param[in] m pointer to a cs_mesh_t structure
889 * \param[in] mq pointer to a cs_mesh_quantities_t structure
890 * \param[in] connect pointer to a cs_cdo_connect_t structure
891 * \param[in,out] quant pointer to a cs_cdo_quantities_t structure
892 */
893 /*----------------------------------------------------------------------------*/
894
895 static void
_mirtich_algorithm(const cs_mesh_t * mesh,const cs_mesh_quantities_t * mq,const cs_cdo_connect_t * connect,cs_cdo_quantities_t * quant)896 _mirtich_algorithm(const cs_mesh_t *mesh,
897 const cs_mesh_quantities_t *mq,
898 const cs_cdo_connect_t *connect,
899 cs_cdo_quantities_t *quant)
900 {
901 const int X = 0, Y = 1, Z = 2;
902 const cs_lnum_t n_cells = quant->n_cells;
903 const cs_lnum_t n_faces = quant->n_faces;
904
905 /* Sanity checks */
906 assert(connect->f2c != NULL);
907 assert(connect->c2f != NULL);
908
909 # pragma omp parallel for if (3*n_cells > CS_THR_MIN)
910 for (cs_lnum_t i = 0; i < 3*n_cells; i++)
911 quant->cell_centers[i] = 0.0;
912
913 for (cs_lnum_t f_id = 0; f_id < n_faces; f_id++) { /* Loop on faces */
914
915 /* Choose gamma to maximize normal according gamma (x, y, or z)
916 Define a direct basis (alpha, beta, gamma) with this choice
917 Compute omega = - <n, P> where P belongs to the face */
918
919 _cdo_fspec_t fspec = _get_fspec(f_id, mesh, mq);
920
921 cs_lnum_t A = fspec.XYZ[X];
922 cs_lnum_t B = fspec.XYZ[Y];
923 cs_lnum_t C = fspec.XYZ[Z];
924
925 _cdo_fsubq_t fsubq = _get_fsub_quantities(f_id,
926 connect, quant->vtx_coord,
927 fspec);
928
929 /* Update cell quantities */
930 const cs_adjacency_t *f2c = connect->f2c;
931 for (cs_lnum_t i = f2c->idx[f_id]; i < f2c->idx[f_id+1]; i++) {
932
933 const cs_lnum_t c_id = f2c->ids[i];
934 const short int sgn = f2c->sgn[i];
935
936 quant->cell_centers[3*c_id + A] += sgn * fspec.q.unitv[A] * fsubq.Fa2;
937 quant->cell_centers[3*c_id + B] += sgn * fspec.q.unitv[B] * fsubq.Fb2;
938 quant->cell_centers[3*c_id + C] += sgn * fspec.q.unitv[C] * fsubq.Fc2;
939
940 } /* End of loop on cell faces */
941
942 } /* End of loop on faces */
943
944 /* Compute cell center of gravity and total volume */
945 # pragma omp parallel for if (n_cells > CS_THR_MIN)
946 for (cs_lnum_t i = 0; i < n_cells; i++) {
947 const double inv_2vol = 0.5 / quant->cell_vol[i];
948 for (int k = 0; k < 3; k++)
949 quant->cell_centers[3*i+k] *= inv_2vol;
950 }
951
952 }
953
954 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
955
956 /*============================================================================
957 * Public function prototypes
958 *============================================================================*/
959
960 /*----------------------------------------------------------------------------*/
961 /*!
962 * \brief Set which quantities have to be computed. Additionnal quantities
963 * are added to cs_cdo_quantities_flag (static variable)
964 *
965 * \param[in] option_flag flag to set geometrical quantities to compute
966 */
967 /*----------------------------------------------------------------------------*/
968
969 void
cs_cdo_quantities_set(cs_flag_t option_flag)970 cs_cdo_quantities_set(cs_flag_t option_flag)
971 {
972 cs_cdo_quantities_flag |= option_flag;
973 }
974
975 /*----------------------------------------------------------------------------*/
976 /*!
977 * \brief Set the type of algorithm to use for computing the cell center
978 *
979 * \param[in] algo type of algorithm
980 */
981 /*----------------------------------------------------------------------------*/
982
983 void
cs_cdo_quantities_set_algo_ccenter(cs_cdo_quantities_cell_center_algo_t algo)984 cs_cdo_quantities_set_algo_ccenter(cs_cdo_quantities_cell_center_algo_t algo)
985 {
986 cs_cdo_cell_center_algo = algo;
987 }
988
989 /*----------------------------------------------------------------------------*/
990 /*!
991 * \brief Build a cs_cdo_quantities_t structure. Some quantities are shared
992 * with the \ref cs_mesh_quantities_t structure and other are not
993 * built according to the given flags in cs_cdo_quantities_flag.
994 *
995 * \param[in] m pointer to a cs_mesh_t structure
996 * \param[in] mq pointer to a cs_mesh_quantities_t structure
997 * \param[in] topo pointer to a cs_cdo_connect_t structure
998 *
999 * \return a new allocated pointer to a cs_cdo_quantities_t structure
1000 */
1001 /*----------------------------------------------------------------------------*/
1002
1003 cs_cdo_quantities_t *
cs_cdo_quantities_build(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,const cs_cdo_connect_t * topo)1004 cs_cdo_quantities_build(const cs_mesh_t *m,
1005 const cs_mesh_quantities_t *mq,
1006 const cs_cdo_connect_t *topo)
1007 {
1008 cs_timer_t t0 = cs_timer_time();
1009
1010 /* Sanity checks */
1011
1012 assert(topo != NULL && topo->c2f != NULL);
1013
1014 cs_cdo_quantities_t *cdoq = _create_cdo_quantities();
1015
1016 /* Compute the volume of the whole domain */
1017
1018 cdoq->vol_tot = mq->tot_vol;
1019
1020 /* Is there a removal of boundary faces to speed-up 2D computations */
1021
1022 if (m->n_g_b_faces_all > m->n_g_b_faces) {
1023
1024 cdoq->remove_boundary_faces = true;
1025
1026 /* The computation of the cell barycenter with the Mirtich algorithm is
1027 false when boundary faces are removed. Switch to a correct one. */
1028
1029 if (cs_cdo_cell_center_algo == CS_CDO_QUANTITIES_BARYC_CENTER)
1030 cs_cdo_cell_center_algo = CS_CDO_QUANTITIES_MEANV_CENTER;
1031
1032 }
1033
1034 /* 1) Initialize shared quantities */
1035 /* ============================ */
1036
1037 /* Face-related quantities */
1038 /* ----------------------- */
1039
1040 /* Shared quantities related to faces (interior and border) */
1041
1042 cdoq->n_i_faces = m->n_i_faces;
1043 cdoq->i_face_normal = mq->i_face_normal;
1044 cdoq->i_face_center = mq->i_face_cog;
1045 cdoq->i_face_surf = mq->i_face_surf;
1046
1047 cdoq->n_b_faces = m->n_b_faces;
1048 cdoq->b_face_normal = mq->b_face_normal;
1049 cdoq->b_face_center = mq->b_face_cog;
1050 cdoq->b_face_surf = mq->b_face_surf;
1051
1052 cdoq->n_faces = m->n_i_faces + m->n_b_faces;
1053 cdoq->n_g_faces = m->n_g_i_faces + m->n_g_b_faces;
1054
1055 /* Vertex-related quantities */
1056 /* ------------------------- */
1057
1058 cdoq->n_vertices = m->n_vertices;
1059 cdoq->n_g_vertices = m->n_g_vertices;
1060 cdoq->vtx_coord = m->vtx_coord;
1061
1062 /* Edge-related quantities */
1063 /* ----------------------- */
1064
1065 cdoq->n_edges = topo->n_edges;
1066 cdoq->n_g_edges = topo->n_g_edges;
1067
1068 /* Cell-related quantities */
1069 /* ----------------------- */
1070
1071 const cs_lnum_t n_cells = m->n_cells;
1072
1073 cdoq->n_cells = n_cells;
1074 cdoq->n_g_cells = m->n_g_cells;
1075 cdoq->cell_vol = mq->cell_vol;
1076
1077 /* Compute the cell centers */
1078
1079 switch (cs_cdo_cell_center_algo) {
1080
1081 case CS_CDO_QUANTITIES_SATURNE_CENTER:
1082 cdoq->cell_centers = mq->cell_cen; /* shared */
1083 break;
1084
1085 case CS_CDO_QUANTITIES_BARYC_CENTER:
1086 BFT_MALLOC(cdoq->cell_centers, 3*n_cells, cs_real_t);
1087 _mirtich_algorithm(m, mq, topo, cdoq);
1088 break;
1089
1090 case CS_CDO_QUANTITIES_MEANV_CENTER:
1091 BFT_MALLOC(cdoq->cell_centers, 3*n_cells, cs_real_t);
1092 _vtx_algorithm(topo, cdoq);
1093 break;
1094
1095 } /* Cell center algorithm */
1096
1097 /* 2) Define quantities available for CDO schemes */
1098 /* =========================================== */
1099
1100 /* Face-related quantities */
1101 /* ----------------------- */
1102
1103 _compute_face_based_quantities(topo, cdoq);
1104
1105 /* Vertex-related quantities */
1106 /* ----------------------- */
1107
1108 /* Compute dual cell volume attached to each vertex in a cell */
1109 _compute_dcell_quantities(topo, cdoq);
1110
1111 /* Edge-related quantities */
1112 /* ----------------------- */
1113
1114 _compute_edge_based_quantities(topo, cdoq);
1115
1116 /* 3) Define metadata */
1117 /* =============== */
1118
1119 /* Define cs_quant_info_t structure */
1120
1121 _compute_quant_info(cdoq);
1122
1123 /* Monitoring */
1124
1125 cs_timer_t t1 = cs_timer_time();
1126 cs_timer_counter_t time_count = cs_timer_diff(&t0, &t1);
1127 cs_log_printf(CS_LOG_PERFORMANCE, " %-35s %9.3f s\n",
1128 "<CDO/Quantities> Runtime", time_count.nsec*1e-9);
1129
1130 return cdoq;
1131 }
1132
1133 /*----------------------------------------------------------------------------*/
1134 /*!
1135 * \brief Destroy a cs_cdo_quantities_t structure
1136 *
1137 * \param[in] cdoq pointer to the cs_cdo_quantities_t struct. to free
1138 *
1139 * \return a NULL pointer
1140 */
1141 /*----------------------------------------------------------------------------*/
1142
1143 cs_cdo_quantities_t *
cs_cdo_quantities_free(cs_cdo_quantities_t * cdoq)1144 cs_cdo_quantities_free(cs_cdo_quantities_t *cdoq)
1145 {
1146 if (cdoq == NULL)
1147 return cdoq;
1148
1149 /* Cell-related quantities */
1150 if (cs_cdo_cell_center_algo != CS_CDO_QUANTITIES_SATURNE_CENTER)
1151 BFT_FREE(cdoq->cell_centers);
1152
1153 /* Face-related quantities */
1154 BFT_FREE(cdoq->dedge_vector);
1155 BFT_FREE(cdoq->pvol_fc);
1156
1157 /* Edge-related quantities */
1158 BFT_FREE(cdoq->edge_vector);
1159 BFT_FREE(cdoq->dface_normal);
1160 BFT_FREE(cdoq->pvol_ec);
1161
1162 /* Vertex-related quantities */
1163 BFT_FREE(cdoq->dcell_vol);
1164
1165 /* vtx_coord is free when the structure cs_mesh_t is destroyed */
1166 BFT_FREE(cdoq);
1167
1168 return NULL;
1169 }
1170
1171 /*----------------------------------------------------------------------------*/
1172 /*!
1173 * \brief Summarize generic information about the cdo mesh quantities
1174 *
1175 * \param[in] quant pointer to cs_cdo_quantities_t structure
1176 */
1177 /*----------------------------------------------------------------------------*/
1178
1179 void
cs_cdo_quantities_summary(const cs_cdo_quantities_t * quant)1180 cs_cdo_quantities_summary(const cs_cdo_quantities_t *quant)
1181 {
1182 cs_log_printf(CS_LOG_SETUP, "\n## CDO quantities settings\n");
1183
1184 switch (cs_cdo_cell_center_algo) {
1185
1186 case CS_CDO_QUANTITIES_SATURNE_CENTER:
1187 cs_log_printf(CS_LOG_SETUP, " * Cell.Center.Algo: Shared with FV\n");
1188 break;
1189
1190 case CS_CDO_QUANTITIES_BARYC_CENTER:
1191 cs_log_printf(CS_LOG_SETUP, " * Cell.Center.Algo: Barycenter (Mirtich)\n");
1192 break;
1193
1194 case CS_CDO_QUANTITIES_MEANV_CENTER:
1195 cs_log_printf(CS_LOG_SETUP, " * Cell.Center.Algo: Mean-value of vertices\n");
1196 break;
1197
1198 } /* Cell center algorithm */
1199
1200 /* Output */
1201 cs_log_printf(CS_LOG_DEFAULT, "\n CDO mesh quantities information:\n");
1202 cs_log_printf(CS_LOG_DEFAULT,
1203 " --cdo-- h_cell %6.4e %6.4e (min/max)\n"
1204 " --cdo-- h_face %6.4e %6.4e (min/max)\n",
1205 quant->cell_info.h_min, quant->cell_info.h_max,
1206 quant->face_info.h_min, quant->face_info.h_max);
1207
1208 if (quant->edge_vector != NULL)
1209 cs_log_printf(CS_LOG_DEFAULT,
1210 " --cdo-- h_edge %6.4e %6.4e (min/max)\n",
1211 quant->edge_info.h_min, quant->edge_info.h_max);
1212 else
1213 cs_log_printf(CS_LOG_DEFAULT, "\n");
1214
1215 #if CS_CDO_QUANTITIES_DBG > 0 && defined(DEBUG) && !defined(NDEBUG)
1216 cs_cdo_quantities_dump(quant);
1217 #endif
1218 }
1219
1220 /*----------------------------------------------------------------------------*/
1221 /*!
1222 * \brief Dump a cs_cdo_quantities_t structure
1223 *
1224 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
1225 *
1226 */
1227 /*----------------------------------------------------------------------------*/
1228
1229 void
cs_cdo_quantities_dump(const cs_cdo_quantities_t * cdoq)1230 cs_cdo_quantities_dump(const cs_cdo_quantities_t *cdoq)
1231 {
1232 int lname = strlen("DumpQuantities.dat") + 1;
1233
1234 /* Define the name of the dump file */
1235 char *fname = NULL;
1236 if (cs_glob_n_ranks > 1) {
1237 lname += 6;
1238 BFT_MALLOC(fname, lname, char);
1239 sprintf(fname, "DumpQuantities.%05d.dat", cs_glob_rank_id);
1240 }
1241 else {
1242 BFT_MALLOC(fname, lname, char);
1243 sprintf(fname, "DumpQuantities.dat");
1244 }
1245 FILE *fdump = fopen(fname, "w");
1246
1247 if (cdoq == NULL) {
1248 fprintf(fdump, "Empty structure.\n");
1249 fclose(fdump);
1250 return;
1251 }
1252
1253 fprintf(fdump, "\n Quantities structure: %p\n\n", (const void *)cdoq);
1254
1255 fprintf(fdump, " -cdoq- n_cells = %ld\n", (long)cdoq->n_cells);
1256 fprintf(fdump, " -cdoq- n_faces = %ld\n", (long)cdoq->n_faces);
1257 fprintf(fdump, " -cdoq- n_edges = %ld\n", (long)cdoq->n_edges);
1258 fprintf(fdump, " -cdoq- n_vertices = %ld\n", (long)cdoq->n_vertices);
1259 fprintf(fdump, " -cdoq- Total volume = %.6e\n\n", cdoq->vol_tot);
1260
1261 fprintf(fdump, "\n *** Cell Quantities ***\n");
1262 fprintf(fdump, "-msg- num.; volume ; center (3)\n");
1263 for (cs_lnum_t i = 0; i < cdoq->n_cells; i++) {
1264 cs_lnum_t p = 3*i;
1265 fprintf(fdump, " [%6ld] | %12.8e | % -12.8e | % -12.8e |% -12.8e\n",
1266 (long)i+1, cdoq->cell_vol[i], cdoq->cell_centers[p],
1267 cdoq->cell_centers[p+1], cdoq->cell_centers[p+2]);
1268 }
1269
1270 fprintf(fdump, "\n\n *** Interior Face Quantities ***\n");
1271 fprintf(fdump, "-msg- num. ; measure ; unitary vector (3) ; center (3)\n");
1272 for (cs_lnum_t f_id = 0; f_id < cdoq->n_i_faces; f_id++) {
1273 cs_quant_t q = cs_quant_set_face(f_id, cdoq);
1274 cs_quant_dump(fdump, f_id+1, q);
1275 }
1276
1277 fprintf(fdump, "\n\n *** Border Face Quantities ***\n");
1278 fprintf(fdump, "-msg- num. ; measure ; unitary vector (3) ; center (3)\n");
1279 for (cs_lnum_t f_id = cdoq->n_i_faces; f_id < cdoq->n_faces; f_id++) {
1280 cs_quant_t q = cs_quant_set_face(f_id, cdoq);
1281 cs_quant_dump(fdump, f_id-cdoq->n_i_faces+1, q);
1282 }
1283
1284 fprintf(fdump, "\n\n *** Edge Quantities ***\n");
1285 fprintf(fdump, "-msg- num. ; measure ; unitary vector (3) ; center (3)\n");
1286 for (cs_lnum_t i = 0; i < cdoq->n_edges; i++) {
1287 const cs_nvec3_t e_vect = cs_quant_set_edge_nvec(i, cdoq);
1288 fprintf(fdump, " -cdoq- [%8ld] | % -10.6e | % -10.6e | % -10.6e |"
1289 " % -10.6e |\n", (long)i+1, e_vect.meas,
1290 e_vect.unitv[0], e_vect.unitv[1], e_vect.unitv[2]);
1291 }
1292
1293 fclose(fdump);
1294 BFT_FREE(fname);
1295 }
1296
1297 /*----------------------------------------------------------------------------*/
1298 /*!
1299 * \brief Compute the portion of volume surrounding each face of a cell.
1300 * This volume corresponds to a pyramid with base f and apex x_f
1301 * The computed quantity is scanned with the c2f adjacency
1302 *
1303 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
1304 * \param[in] c2f pointer to the cell --> edges connectivity
1305 * \param[in, out] p_pvol_fc double pointer to the face volume in each cell
1306 * If not allocated before calling this function,
1307 * one allocates the array storing the volumes
1308 */
1309 /*----------------------------------------------------------------------------*/
1310
1311 void
cs_cdo_quantities_compute_pvol_fc(const cs_cdo_quantities_t * cdoq,const cs_adjacency_t * c2f,cs_real_t ** p_pvol_fc)1312 cs_cdo_quantities_compute_pvol_fc(const cs_cdo_quantities_t *cdoq,
1313 const cs_adjacency_t *c2f,
1314 cs_real_t **p_pvol_fc)
1315 {
1316 if (cdoq == NULL || c2f == NULL)
1317 bft_error(__FILE__, __LINE__, 0,
1318 " %s: A mandatory structure is not allocated.\n", __func__);
1319
1320 const cs_lnum_t n_cells = cdoq->n_cells;
1321
1322 cs_real_t *pvol_fc = *p_pvol_fc;
1323
1324 /* Initialize array */
1325 if (pvol_fc == NULL)
1326 BFT_MALLOC(pvol_fc, c2f->idx[n_cells], cs_real_t);
1327
1328 #if defined(DEBUG) && !defined(NDEBUG)
1329 memset(pvol_fc, 0, c2f->idx[n_cells]*sizeof(cs_real_t));
1330 #endif
1331
1332 # pragma omp parallel for if (n_cells > CS_THR_MIN)
1333 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
1334 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++) {
1335
1336 const cs_lnum_t f_id = c2f->ids[j];
1337 const cs_nvec3_t fp_nvec = cs_quant_set_face_nvec(f_id, cdoq);
1338 const cs_nvec3_t ed_nvec = cs_quant_set_dedge_nvec(j, cdoq);
1339
1340 cs_real_t p_fc = _dp3(fp_nvec.unitv, ed_nvec.unitv);
1341 p_fc *= cs_math_1ov3 * fp_nvec.meas * ed_nvec.meas;
1342
1343 pvol_fc[j] = p_fc;
1344
1345 } /* Loop on cell faces */
1346 } /* Loop on cells */
1347
1348 /* Return pointer */
1349 *p_pvol_fc = pvol_fc;
1350 }
1351
1352 /*----------------------------------------------------------------------------*/
1353 /*!
1354 * \brief Compute the portion of volume surrounding each edge of a cell
1355 * The computed quantity is scanned with the c2e adjacency
1356 *
1357 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
1358 * \param[in] c2e pointer to the cell --> edges connectivity
1359 * \param[in, out] p_pvol_ec double pointer to the edge volume in each cell
1360 * If not allocated before calling this function,
1361 * one allocates the array storing the volumes
1362 */
1363 /*----------------------------------------------------------------------------*/
1364
1365 void
cs_cdo_quantities_compute_pvol_ec(const cs_cdo_quantities_t * cdoq,const cs_adjacency_t * c2e,cs_real_t ** p_pvol_ec)1366 cs_cdo_quantities_compute_pvol_ec(const cs_cdo_quantities_t *cdoq,
1367 const cs_adjacency_t *c2e,
1368 cs_real_t **p_pvol_ec)
1369 {
1370 if (cdoq == NULL || c2e == NULL)
1371 bft_error(__FILE__, __LINE__, 0,
1372 " %s: A mandatory structure is not allocated.\n", __func__);
1373
1374 const cs_lnum_t n_cells = cdoq->n_cells;
1375
1376 cs_real_t *pvol_ec = *p_pvol_ec;
1377
1378 /* Initialize array */
1379 if (pvol_ec == NULL)
1380 BFT_MALLOC(pvol_ec, c2e->idx[n_cells], cs_real_t);
1381
1382 if (cdoq->pvol_ec != NULL)
1383 memcpy(pvol_ec, cdoq->pvol_ec, c2e->idx[n_cells]*sizeof(cs_real_t));
1384
1385 else {
1386
1387 #if defined(DEBUG) && !defined(NDEBUG)
1388 memset(pvol_ec, 0, c2e->idx[n_cells]*sizeof(cs_real_t));
1389 #endif
1390
1391 # pragma omp parallel for if (n_cells > CS_THR_MIN)
1392 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
1393
1394 const cs_lnum_t start = c2e->idx[c_id];
1395 const cs_real_t *sface = cdoq->dface_normal + 3*start;
1396 const cs_lnum_t *c2e_ids = c2e->ids + start;
1397
1398 cs_real_t *_pvol = pvol_ec + start;
1399 for (cs_lnum_t j = 0; j < c2e->idx[c_id+1]-start; j++)
1400 _pvol[j] = cs_math_1ov3 * _dp3(sface + 3*j,
1401 cdoq->edge_vector + 3*c2e_ids[j]);
1402
1403 } /* Loop on cells */
1404
1405 } /* Not already existing */
1406
1407 /* Return pointer */
1408 *p_pvol_ec = pvol_ec;
1409 }
1410
1411 /*----------------------------------------------------------------------------*/
1412 /*!
1413 * \brief Compute the dual volume surrounding each vertex
1414 *
1415 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
1416 * \param[in] c2v pointer to the cell-->vertices connectivity
1417 * \param[in, out] dual_vol dual volumes related to each vertex
1418 */
1419 /*----------------------------------------------------------------------------*/
1420
1421 void
cs_cdo_quantities_compute_dual_volumes(const cs_cdo_quantities_t * cdoq,const cs_adjacency_t * c2v,cs_real_t * dual_vol)1422 cs_cdo_quantities_compute_dual_volumes(const cs_cdo_quantities_t *cdoq,
1423 const cs_adjacency_t *c2v,
1424 cs_real_t *dual_vol)
1425 {
1426 if (dual_vol == NULL)
1427 return;
1428
1429 assert(cdoq != NULL && c2v != NULL);
1430
1431 memset(dual_vol, 0, cdoq->n_vertices*sizeof(cs_real_t));
1432
1433 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++)
1434 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++)
1435 dual_vol[c2v->ids[j]] += cdoq->dcell_vol[j];
1436 }
1437
1438 /*----------------------------------------------------------------------------*/
1439 /*!
1440 * \brief Compute the area of the triangles with basis each edge of the face
1441 * and apex the face center.
1442 * Case of interior faces.
1443 * Storage in agreement with the if2v adjacency structure
1444 *
1445 * \param[in] connect pointer to a cs_cdo_connect_t structure
1446 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1447 * \param[in] f_id interior face id
1448 * \param[in, out] tef quantities to compute (pre-allocated)
1449 */
1450 /*----------------------------------------------------------------------------*/
1451
1452 void
cs_cdo_quantities_compute_i_tef(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_lnum_t f_id,cs_real_t tef[])1453 cs_cdo_quantities_compute_i_tef(const cs_cdo_connect_t *connect,
1454 const cs_cdo_quantities_t *cdoq,
1455 cs_lnum_t f_id,
1456 cs_real_t tef[])
1457 {
1458 if (tef == NULL) return;
1459
1460 const cs_real_t *xf = cdoq->i_face_center + 3*f_id;
1461 const cs_lnum_t *idx = connect->if2v->idx + f_id;
1462 const cs_lnum_t *ids = connect->if2v->ids + idx[0];
1463 const int n_ef = idx[1] - idx[0]; /* n_ef = n_vf */
1464
1465 cs_lnum_t _v0, _v1;
1466 for (int e = 0; e < n_ef; e++) { /* Compute */
1467
1468 if (e < n_ef - 1)
1469 _v0 = ids[e], _v1 = ids[e+1];
1470 else
1471 _v0 = ids[n_ef-1], _v1 = ids[0];
1472
1473 tef[e] = cs_math_surftri(cdoq->vtx_coord + 3*_v0,
1474 cdoq->vtx_coord + 3*_v1,
1475 xf);
1476
1477 } /* Loop on face edges */
1478 }
1479
1480 /*----------------------------------------------------------------------------*/
1481 /*!
1482 * \brief Compute the area of the triangles with basis each edge of the face
1483 * and apex the face center.
1484 * Case of boundary faces.
1485 * Storage in agreement with the bf2v adjacency structure
1486 *
1487 * \param[in] connect pointer to a cs_cdo_connect_t structure
1488 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1489 * \param[in] bf_id border face id
1490 * \param[in, out] tef quantities to compute (pre-allocated)
1491 */
1492 /*----------------------------------------------------------------------------*/
1493
1494 void
cs_cdo_quantities_compute_b_tef(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_lnum_t bf_id,cs_real_t tef[])1495 cs_cdo_quantities_compute_b_tef(const cs_cdo_connect_t *connect,
1496 const cs_cdo_quantities_t *cdoq,
1497 cs_lnum_t bf_id,
1498 cs_real_t tef[])
1499 {
1500 if (tef == NULL) return;
1501
1502 const cs_real_t *xf = cdoq->b_face_center + 3*bf_id;
1503 const cs_lnum_t *idx = connect->bf2v->idx + bf_id;
1504 const cs_lnum_t *ids = connect->bf2v->ids + idx[0];
1505 const int n_ef = idx[1] - idx[0]; /* n_ef = n_vf */
1506
1507 cs_lnum_t _v0, _v1;
1508 for (int e = 0; e < n_ef; e++) { /* Compute */
1509
1510 if (e < n_ef - 1)
1511 _v0 = ids[e], _v1 = ids[e+1];
1512 else
1513 _v0 = ids[n_ef-1], _v1 = ids[0];
1514
1515 tef[e] = cs_math_surftri(cdoq->vtx_coord + 3*_v0,
1516 cdoq->vtx_coord + 3*_v1,
1517 xf);
1518
1519 } /* Loop on face edges */
1520 }
1521
1522 /*----------------------------------------------------------------------------*/
1523 /*!
1524 * \brief Compute the weight related to each vertex of a face. This weight
1525 * ensures a 2nd order approximation if the face center is the face
1526 * barycenter.
1527 * Case of interior faces.
1528 *
1529 * \param[in] connect pointer to a cs_cdo_connect_t structure
1530 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1531 * \param[in] f_id interior face id
1532 * \param[in, out] wvf quantities to compute (pre-allocated)
1533 */
1534 /*----------------------------------------------------------------------------*/
1535
1536 void
cs_cdo_quantities_compute_i_wvf(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_lnum_t f_id,cs_real_t wvf[])1537 cs_cdo_quantities_compute_i_wvf(const cs_cdo_connect_t *connect,
1538 const cs_cdo_quantities_t *cdoq,
1539 cs_lnum_t f_id,
1540 cs_real_t wvf[])
1541 {
1542 if (wvf == NULL) return;
1543
1544 const cs_real_t *xf = cdoq->i_face_center + 3*f_id;
1545 const cs_lnum_t *idx = connect->if2v->idx + f_id;
1546 const cs_lnum_t *ids = connect->if2v->ids + idx[0];
1547 const int n_vf = idx[1] - idx[0];
1548
1549 for (cs_lnum_t v = 0; v < n_vf; v++) wvf[v] = 0.; /* Init */
1550
1551 int _v0, _v1;
1552 for (cs_lnum_t v = 0; v < n_vf; v++) { /* Compute */
1553
1554 if (v < n_vf - 1)
1555 _v0 = v, _v1 = v+1;
1556 else
1557 _v0 = n_vf-1, _v1 = 0;
1558
1559 const double tef = cs_math_surftri(cdoq->vtx_coord + 3*ids[_v0],
1560 cdoq->vtx_coord + 3*ids[_v1],
1561 xf);
1562 wvf[_v0] += tef;
1563 wvf[_v1] += tef;
1564
1565 } /* Loop on face vertices */
1566
1567 const cs_real_t coef = 0.5/cdoq->i_face_surf[f_id];
1568 for (cs_lnum_t v = 0; v < n_vf; v++) wvf[v] *= coef;
1569 }
1570
1571 /*----------------------------------------------------------------------------*/
1572 /*!
1573 * \brief Compute the weight related to each vertex of a face. This weight
1574 * ensures a 2nd order approximation if the face center is the face
1575 * barycenter.
1576 * Case of boundary faces.
1577 *
1578 * \param[in] connect pointer to a cs_cdo_connect_t structure
1579 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1580 * \param[in] bf_id border face id
1581 * \param[in, out] wvf quantities to compute (pre-allocated)
1582 */
1583 /*----------------------------------------------------------------------------*/
1584
1585 void
cs_cdo_quantities_compute_b_wvf(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * cdoq,cs_lnum_t bf_id,cs_real_t wvf[])1586 cs_cdo_quantities_compute_b_wvf(const cs_cdo_connect_t *connect,
1587 const cs_cdo_quantities_t *cdoq,
1588 cs_lnum_t bf_id,
1589 cs_real_t wvf[])
1590 {
1591 if (wvf == NULL) return;
1592
1593 const cs_real_t *xf = cdoq->b_face_center + 3*bf_id;
1594 const cs_lnum_t *idx = connect->bf2v->idx + bf_id;
1595 const cs_lnum_t *ids = connect->bf2v->ids + idx[0];
1596 const int n_vf = idx[1] - idx[0];
1597
1598 for (cs_lnum_t v = 0; v < n_vf; v++) wvf[v] = 0.; /* Init */
1599
1600 int _v0, _v1;
1601 for (cs_lnum_t v = 0; v < n_vf; v++) { /* Compute */
1602
1603 if (v < n_vf - 1)
1604 _v0 = v, _v1 = v+1;
1605 else
1606 _v0 = n_vf-1, _v1 = 0;
1607
1608 const double tef = cs_math_surftri(cdoq->vtx_coord + 3*ids[_v0],
1609 cdoq->vtx_coord + 3*ids[_v1],
1610 xf);
1611 wvf[_v0] += tef;
1612 wvf[_v1] += tef;
1613
1614 } /* Loop on face vertices */
1615
1616 const cs_real_t coef = 0.5/cdoq->b_face_surf[bf_id];
1617 for (cs_lnum_t v = 0; v < n_vf; v++) wvf[v] *= coef;
1618 }
1619
1620 /*----------------------------------------------------------------------------*/
1621 /*!
1622 * \brief Define a cs_quant_t structure for a primal face (interior or border)
1623 *
1624 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
1625 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1626 *
1627 * \return a initialize structure
1628 */
1629 /*----------------------------------------------------------------------------*/
1630
1631 cs_quant_t
cs_quant_set_face(cs_lnum_t f_id,const cs_cdo_quantities_t * cdoq)1632 cs_quant_set_face(cs_lnum_t f_id,
1633 const cs_cdo_quantities_t *cdoq)
1634 {
1635 cs_quant_t q = {.meas = 0.,
1636 .unitv[0] = 0., .unitv[1] = 0., .unitv[2] = 0.,
1637 .center[0] = 0., .center[1] = 0., .center[2] = 0.};
1638 cs_nvec3_t nv;
1639
1640 if (f_id < cdoq->n_i_faces) { /* Interior face */
1641
1642 q.meas = cdoq->i_face_surf[f_id];
1643
1644 cs_nvec3(cdoq->i_face_normal + 3*f_id, &nv);
1645 for (int k = 0; k < 3; k++)
1646 q.unitv[k] = nv.unitv[k];
1647
1648 const cs_real_t *xf = cdoq->i_face_center + 3*f_id;
1649 for (int k = 0; k < 3; k++)
1650 q.center[k] = xf[k];
1651
1652 }
1653 else { /* Border face */
1654
1655 const cs_lnum_t bf_id = f_id - cdoq->n_i_faces;
1656
1657 q.meas = cdoq->b_face_surf[bf_id];
1658 cs_nvec3(cdoq->b_face_normal + 3*bf_id, &nv);
1659 for (int k = 0; k < 3; k++)
1660 q.unitv[k] = nv.unitv[k];
1661
1662 const cs_real_t *xf = cdoq->b_face_center + 3*bf_id;
1663 for (int k = 0; k < 3; k++)
1664 q.center[k] = xf[k];
1665
1666 }
1667
1668 return q;
1669 }
1670
1671 /*----------------------------------------------------------------------------*/
1672 /*!
1673 * \brief Retrieve the face surface and its unit normal vector for a primal
1674 * face (interior or border)
1675 *
1676 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
1677 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1678 *
1679 * \return a pointer to the face normalized vector
1680 */
1681 /*----------------------------------------------------------------------------*/
1682
1683 cs_nvec3_t
cs_quant_set_face_nvec(cs_lnum_t f_id,const cs_cdo_quantities_t * cdoq)1684 cs_quant_set_face_nvec(cs_lnum_t f_id,
1685 const cs_cdo_quantities_t *cdoq)
1686 {
1687 cs_nvec3_t nv;
1688
1689 if (f_id < cdoq->n_i_faces) /* Interior face */
1690 cs_nvec3(cdoq->i_face_normal + 3*f_id, &nv);
1691 else /* Border face */
1692 cs_nvec3(cdoq->b_face_normal + 3*(f_id - cdoq->n_i_faces), &nv);
1693
1694 return nv;
1695 }
1696
1697 /*----------------------------------------------------------------------------*/
1698 /*!
1699 * \brief Get the normalized vector associated to a primal edge
1700 *
1701 * \param[in] e_id id related to an edge
1702 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1703 *
1704 * \return a pointer to the edge normalized vector
1705 */
1706 /*----------------------------------------------------------------------------*/
1707
1708 cs_nvec3_t
cs_quant_set_edge_nvec(cs_lnum_t e_id,const cs_cdo_quantities_t * cdoq)1709 cs_quant_set_edge_nvec(cs_lnum_t e_id,
1710 const cs_cdo_quantities_t *cdoq)
1711 {
1712 cs_nvec3_t nv;
1713 cs_nvec3(cdoq->edge_vector + 3*e_id, &nv);
1714
1715 return nv;
1716 }
1717
1718 /*----------------------------------------------------------------------------*/
1719 /*!
1720 * \brief Get the two normalized vector associated to a dual edge
1721 *
1722 * \param[in] shift position in c2f_idx
1723 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
1724 *
1725 * \return a pointer to the dual edge normalized vector
1726 */
1727 /*----------------------------------------------------------------------------*/
1728
1729 cs_nvec3_t
cs_quant_set_dedge_nvec(cs_lnum_t shift,const cs_cdo_quantities_t * cdoq)1730 cs_quant_set_dedge_nvec(cs_lnum_t shift,
1731 const cs_cdo_quantities_t *cdoq)
1732 {
1733 cs_nvec3_t nv;
1734 cs_nvec3(cdoq->dedge_vector + 3*shift, &nv);
1735
1736 return nv;
1737 }
1738
1739 /*----------------------------------------------------------------------------*/
1740 /*!
1741 * \brief Dump a cs_quant_t structure
1742 *
1743 * \param[in] f FILE struct (stdout if NULL)
1744 * \param[in] num entity number related to this quantity struct.
1745 * \param[in] q cs_quant_t structure to dump
1746 */
1747 /*----------------------------------------------------------------------------*/
1748
1749 void
cs_quant_dump(FILE * f,cs_lnum_t num,const cs_quant_t q)1750 cs_quant_dump(FILE *f,
1751 cs_lnum_t num,
1752 const cs_quant_t q)
1753 {
1754 FILE *_f = f;
1755
1756 if (_f == NULL) _f = stdout;
1757
1758 fprintf(_f, " -cdoq- [%8ld] | % -10.6e | % -10.6e | % -10.6e | % -10.6e |"
1759 " % -10.6e | % -10.6e | % -10.6e\n", (long)num, q.meas,
1760 q.unitv[0], q.unitv[1], q.unitv[2], q.center[0], q.center[1],
1761 q.center[2]);
1762 }
1763
1764 /*----------------------------------------------------------------------------*/
1765
1766 #undef _dp3
1767 #undef _n3
1768
1769 END_C_DECLS
1770