1 #ifndef __CS_CDO_QUANTITIES_H__
2 #define __CS_CDO_QUANTITIES_H__
3
4 /*============================================================================
5 * Manage geometrical quantities needed in CDO schemes
6 *============================================================================*/
7
8 /*
9 This file is part of Code_Saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2021 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27
28 /*----------------------------------------------------------------------------*/
29
30 /*----------------------------------------------------------------------------
31 * Local headers
32 *----------------------------------------------------------------------------*/
33
34 #include "cs_base.h"
35 #include "cs_cdo_connect.h"
36 #include "cs_flag.h"
37 #include "cs_math.h"
38 #include "cs_mesh.h"
39 #include "cs_mesh_quantities.h"
40
41 /*----------------------------------------------------------------------------*/
42
43 BEGIN_C_DECLS
44
45 /*============================================================================
46 * Macro definitions
47 *============================================================================*/
48
49 /*============================================================================
50 * Type definitions
51 *============================================================================*/
52
53 /*! \enum cs_cdo_quantities_cell_center_algo_t
54 * \brief Type of algorithm used to compute the cell centers
55 *
56 * \var CS_CDO_QUANTITIES_MEANV_CENTER
57 * Center is computed as the mean of cell vertices
58 *
59 * \var CS_CDO_QUANTITIES_BARYC_CENTER
60 * Center is computed as the real cell barycenter (default choice)
61 *
62 * \var CS_CDO_QUANTITIES_SATURNE_CENTER
63 * Center is the one defined in cs_mesh_quantities_t (i.e. the one
64 * used in the legacy Finite Volume scheme).
65 */
66
67 typedef enum {
68
69 CS_CDO_QUANTITIES_MEANV_CENTER,
70 CS_CDO_QUANTITIES_BARYC_CENTER,
71 CS_CDO_QUANTITIES_SATURNE_CENTER
72
73 } cs_cdo_quantities_cell_center_algo_t;
74
75 /*! \enum cs_cdo_quantities_bit_t
76 * \brief Bit values for setting which quantities to compute
77 *
78 * \var CS_CDO_QUANTITIES_EB_SCHEME
79 * Geometrical quantities related to edge-based schemes
80 *
81 * \var CS_CDO_QUANTITIES_FB_SCHEME
82 * Geometrical quantities related to face-based schemes
83 *
84 * \var CS_CDO_QUANTITIES_HHO_SCHEME
85 * Geometrical quantities related to HHO schemes
86 *
87 * \var CS_CDO_QUANTITIES_VB_SCHEME
88 * Geometrical quantities related to vertex-based schemes
89 *
90 * \var CS_CDO_QUANTITIES_VCB_SCHEME
91 * Geometrical quantities related to vertex+cell-based schemes
92 */
93
94 typedef enum {
95
96 /* Set of geometrical quantities related to CDO schemes */
97
98 CS_CDO_QUANTITIES_EB_SCHEME = 1<<0, /* = 1 */
99 CS_CDO_QUANTITIES_FB_SCHEME = 1<<1, /* = 2 */
100 CS_CDO_QUANTITIES_HHO_SCHEME = 1<<2, /* = 4 */
101 CS_CDO_QUANTITIES_VB_SCHEME = 1<<3, /* = 8 */
102 CS_CDO_QUANTITIES_VCB_SCHEME = 1<<4, /* = 16 */
103
104 } cs_cdo_quantities_bit_t;
105
106
107 /* Structure storing information about variation of entities across the
108 mesh for a given type of entity (cell, face and edge) */
109 typedef struct {
110
111 /* Measure is either a volume for cells, a surface for faces or a length
112 for edges */
113 double meas_min; /* Min. value of the entity measure */
114 double meas_max; /* Max. value of the entity measure */
115 double h_min; /* Estimation of the min. value of the diameter */
116 double h_max; /* Estimation of the max. value of the diameter */
117
118 } cs_quant_info_t;
119
120 /* For primal vector quantities (edge or face) */
121 typedef struct {
122
123 double meas; /* length or area */
124 double unitv[3]; /* unitary vector: tangent or normal to the element */
125 double center[3];
126
127 } cs_quant_t;
128
129 typedef struct { /* Specific mesh quantities */
130
131 /* Keep the information about the removal of boundary faces in case of 2D
132 computations */
133 bool remove_boundary_faces;
134
135 /* Global mesh quantities */
136 double vol_tot;
137
138 /* Cell-based quantities */
139 /* ===================== */
140
141 cs_lnum_t n_cells; /* Local number of cells */
142 cs_gnum_t n_g_cells; /* Global number of cells */
143 cs_real_t *cell_centers; /* May be shared according to options */
144 const cs_real_t *cell_vol; /* Shared with cs_mesh_quantities_t */
145
146 cs_quant_info_t cell_info;
147
148 /* Face-based quantities */
149 /* ===================== */
150
151 cs_lnum_t n_faces; /* n_i_faces + n_b_faces */
152 cs_lnum_t n_i_faces; /* Local number of interior faces */
153 cs_lnum_t n_b_faces; /* Local number of border faces */
154 cs_gnum_t n_g_faces; /* Global number of faces */
155
156 /* Remark: cs_quant_t structure attached to a face (interior or border) can
157 be built on-the-fly calling the function cs_quant_set_face(f_id, cdoq).
158 See \ref cs_quant_set_face for more details.
159
160 In order to reduce the memory consumption one shares face quantities with
161 the ones defined in the legacy part and stored in the cs_mesh_quantities_t
162 structure that's why a distinction is made between interior and border
163 faces.
164
165 cs_nvec3_t structure associated to a face can also be built on-the-fly
166 using cs_quant_set_face_nvec(f_id, cdoq).
167 See \ref cs_quant_set_face_nvec for more details.
168 */
169
170 const cs_real_t *i_face_normal; /* Shared with cs_mesh_quantities_t */
171 const cs_real_t *i_face_center; /* Shared with cs_mesh_quantities_t */
172 const cs_real_t *i_face_surf; /* Shared with cs_mesh_quantities_t */
173
174 const cs_real_t *b_face_normal; /* Shared with cs_mesh_quantities_t */
175 const cs_real_t *b_face_center; /* Shared with cs_mesh_quantities_t */
176 const cs_real_t *b_face_surf; /* Shared with cs_mesh_quantities_t */
177
178 /* Remark: cs_nvec3_t structure attached to a dual edge can be built
179 on-the-fly to access to its length and its unit tangential vector using
180 the function cs_quant_set_dedge_nvec(shift, cdoq)
181
182 One recalls that a dual edge is associated to a primal face and is shared
183 with two cells for an interior face and shared with one cell for a
184 boundary face. Scan this quantity with the c2f connectivity.
185 */
186
187 cs_real_t *dedge_vector; /* Allocation to 3*c2f->idx[n_faces] */
188
189 cs_real_t *pvol_fc; /* Portion of volume surrounding a face
190 * in each cell. This is a pyramid of
191 * base the face and apex the cell center
192 * Scanned with the c2f adjacency.
193 * Not always allocated.
194 */
195 cs_quant_info_t face_info;
196
197 /* Edge-based quantities */
198 /* ===================== */
199
200 cs_lnum_t n_edges; /* Local number of edges */
201 cs_gnum_t n_g_edges; /* Global number of edges */
202
203 cs_real_t *edge_vector; /* Allocation to 3*n_edges
204 Norm of the vector is equal to the
205 distance between two vertices.
206 Unit vector is the tangential direction
207 attached to the edge */
208
209 /* For each edge e belonging to a cell c, the dual face is built from the
210 contributions of two triangles s(x_c, x_f, x_e) and s(x_c, x_f', x_e) with
211 the faces f and f' belonging to F_e \cap F_c
212 Scan this quantity with the c2e connectivity */
213
214 cs_real_t *dface_normal; /* Vector-valued normal for each dual face
215 * inside a cell associated to an edge */
216 cs_real_t *pvol_ec; /* Portion of volume surrounding an edge
217 * in each cell. Scanned with the c2e
218 * adjacency.
219 * Not always allocated. */
220
221 cs_quant_info_t edge_info;
222
223 /* Vertex-based quantities */
224 /* ======================= */
225
226 cs_lnum_t n_vertices; /* Local number of vertices */
227 cs_gnum_t n_g_vertices; /* Global number of vertices */
228
229 cs_real_t *dcell_vol; /* Dual volume related to each vertex.
230 * Scanned with the c2v adjacency.
231 * Not always allocated.
232 */
233 const cs_real_t *vtx_coord; /* Shared with the cs_mesh_t structure */
234
235 } cs_cdo_quantities_t;
236
237 /*============================================================================
238 * Global variables
239 *============================================================================*/
240
241 /*============================================================================
242 * Public function prototypes
243 *============================================================================*/
244
245 /*----------------------------------------------------------------------------*/
246 /*!
247 * \brief Compute the area of the triangle of base given by q (related to a
248 * segment) with apex located at xa
249 *
250 * \param[in] qa pointer to a cs_quant_t structure related to a segment
251 * \param[in] xb coordinates of the apex to consider
252 *
253 * \return the value the area of the triangle
254 */
255 /*----------------------------------------------------------------------------*/
256
257 static inline double
cs_compute_area_from_quant(const cs_quant_t qa,const cs_real_t * xb)258 cs_compute_area_from_quant(const cs_quant_t qa,
259 const cs_real_t *xb)
260 {
261 const double xab[3] = {xb[0] - qa.center[0],
262 xb[1] - qa.center[1],
263 xb[2] - qa.center[2]};
264 const double cp[3] = {qa.unitv[1]*xab[2] - qa.unitv[2]*xab[1],
265 qa.unitv[2]*xab[0] - qa.unitv[0]*xab[2],
266 qa.unitv[0]*xab[1] - qa.unitv[1]*xab[0]};
267
268 return 0.5 * qa.meas * cs_math_3_norm(cp);
269 }
270
271 /*----------------------------------------------------------------------------*/
272 /*!
273 * \brief Set which quantities have to be computed. Additionnal quantities
274 * are added to cs_cdo_quantities_flag (static variable)
275 *
276 * \param[in] option_flag flag to set geometrical quantities to compute
277 */
278 /*----------------------------------------------------------------------------*/
279
280 void
281 cs_cdo_quantities_set(cs_flag_t option_flag);
282
283 /*----------------------------------------------------------------------------*/
284 /*!
285 * \brief Set the type of algorithm to use for computing the cell center
286 *
287 * \param[in] algo type of algorithm
288 */
289 /*----------------------------------------------------------------------------*/
290
291 void
292 cs_cdo_quantities_set_algo_ccenter(cs_cdo_quantities_cell_center_algo_t algo);
293
294 /*----------------------------------------------------------------------------*/
295 /*!
296 * \brief Build a cs_cdo_quantities_t structure. Some quantities are shared
297 * with the \ref cs_mesh_quantities_t structure and other are not
298 * built according to the given flags in cs_cdo_quantities_flag.
299 *
300 * \param[in] m pointer to a cs_mesh_t structure
301 * \param[in] mq pointer to a cs_mesh_quantities_t structure
302 * \param[in] topo pointer to a cs_cdo_connect_t structure
303 *
304 * \return a new allocated pointer to a cs_cdo_quantities_t structure
305 */
306 /*----------------------------------------------------------------------------*/
307
308 cs_cdo_quantities_t *
309 cs_cdo_quantities_build(const cs_mesh_t *m,
310 const cs_mesh_quantities_t *mq,
311 const cs_cdo_connect_t *topo);
312
313 /*----------------------------------------------------------------------------*/
314 /*!
315 * \brief Destroy a cs_cdo_quantities_t structure
316 *
317 * \param[in] q pointer to the cs_cdo_quantities_t struct. to free
318 *
319 * \return a NULL pointer
320 */
321 /*----------------------------------------------------------------------------*/
322
323 cs_cdo_quantities_t *
324 cs_cdo_quantities_free(cs_cdo_quantities_t *q);
325
326 /*----------------------------------------------------------------------------*/
327 /*!
328 * \brief Summarize generic information about the cdo mesh quantities
329 *
330 * \param[in] quant pointer to cs_cdo_quantities_t structure
331 *
332 */
333 /*----------------------------------------------------------------------------*/
334
335 void
336 cs_cdo_quantities_summary(const cs_cdo_quantities_t *quant);
337
338 /*----------------------------------------------------------------------------*/
339 /*!
340 * \brief Dump a cs_cdo_quantities_t structure
341 *
342 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
343 *
344 */
345 /*----------------------------------------------------------------------------*/
346
347 void
348 cs_cdo_quantities_dump(const cs_cdo_quantities_t *cdoq);
349
350 /*----------------------------------------------------------------------------*/
351 /*!
352 * \brief Compute the portion of volume surrounding each face of a cell.
353 * This volume corresponds to a pyramid with base f and apex x_f
354 * The computed quantity is scanned with the c2f adjacency
355 *
356 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
357 * \param[in] c2f pointer to the cell --> edges connectivity
358 * \param[in, out] p_pvol_fc double pointer to the face volume in each cell
359 * If not allocated before calling this function,
360 * one allocates the array storing the volumes
361 */
362 /*----------------------------------------------------------------------------*/
363
364 void
365 cs_cdo_quantities_compute_pvol_fc(const cs_cdo_quantities_t *cdoq,
366 const cs_adjacency_t *c2f,
367 cs_real_t **p_pvol_fc);
368
369 /*----------------------------------------------------------------------------*/
370 /*!
371 * \brief Compute the portion of volume surrounding each edge of a cell
372 * The computed quantity is scanned with the c2e adjacency
373 *
374 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
375 * \param[in] c2e pointer to the cell --> edges connectivity
376 * \param[in, out] p_pvol_ec double pointer to the edge volume in each cell
377 * If not allocated before calling this function,
378 * one allocates the array storing the volumes
379 */
380 /*----------------------------------------------------------------------------*/
381
382 void
383 cs_cdo_quantities_compute_pvol_ec(const cs_cdo_quantities_t *cdoq,
384 const cs_adjacency_t *c2e,
385 cs_real_t **p_pvol_ec);
386
387 /*----------------------------------------------------------------------------*/
388 /*!
389 * \brief Compute the dual volume surrounding each vertex
390 *
391 * \param[in] cdoq pointer to cs_cdo_quantities_t structure
392 * \param[in] c2v pointer to the cell-->vertices connectivity
393 * \param[in, out] dual_vol dual volumes related to each vertex
394 */
395 /*----------------------------------------------------------------------------*/
396
397 void
398 cs_cdo_quantities_compute_dual_volumes(const cs_cdo_quantities_t *cdoq,
399 const cs_adjacency_t *c2v,
400 cs_real_t *dual_vol);
401
402 /*----------------------------------------------------------------------------*/
403 /*!
404 * \brief Compute the area of the triangles with basis each edge of the face
405 * and apex the face center.
406 * Case of interior faces.
407 * Storage in agreement with the bf2v adjacency structure
408 *
409 * \param[in] connect pointer to a cs_cdo_connect_t structure
410 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
411 * \param[in] f_id interior face id
412 * \param[in, out] tef quantities to compute (pre-allocated)
413 */
414 /*----------------------------------------------------------------------------*/
415
416 void
417 cs_cdo_quantities_compute_i_tef(const cs_cdo_connect_t *connect,
418 const cs_cdo_quantities_t *cdoq,
419 cs_lnum_t f_id,
420 cs_real_t tef[]);
421
422 /*----------------------------------------------------------------------------*/
423 /*!
424 * \brief Compute the area of the triangles with basis each edge of the face
425 * and apex the face center.
426 * Case of boundary faces.
427 * Storage in agreement with the bf2v adjacency structure
428 *
429 * \param[in] connect pointer to a cs_cdo_connect_t structure
430 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
431 * \param[in] bf_id border face id
432 * \param[in, out] tef quantities to compute (pre-allocated)
433 */
434 /*----------------------------------------------------------------------------*/
435
436 void
437 cs_cdo_quantities_compute_b_tef(const cs_cdo_connect_t *connect,
438 const cs_cdo_quantities_t *cdoq,
439 cs_lnum_t bf_id,
440 cs_real_t tef[]);
441
442 /*----------------------------------------------------------------------------*/
443 /*!
444 * \brief Compute the weight related to each vertex of a face. This weight
445 * ensures a 2nd order approximation if the face center is the face
446 * barycenter.
447 * Case of interior faces.
448 *
449 * \param[in] connect pointer to a cs_cdo_connect_t structure
450 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
451 * \param[in] f_id interior face id
452 * \param[in, out] wvf quantities to compute (pre-allocated)
453 */
454 /*----------------------------------------------------------------------------*/
455
456 void
457 cs_cdo_quantities_compute_i_wvf(const cs_cdo_connect_t *connect,
458 const cs_cdo_quantities_t *cdoq,
459 cs_lnum_t f_id,
460 cs_real_t wvf[]);
461
462 /*----------------------------------------------------------------------------*/
463 /*!
464 * \brief Compute the weight related to each vertex of a face. This weight
465 * ensures a 2nd order approximation if the face center is the face
466 * barycenter.
467 * Case of boundary faces.
468 *
469 * \param[in] connect pointer to a cs_cdo_connect_t structure
470 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
471 * \param[in] bf_id border face id
472 * \param[in, out] wvf quantities to compute (pre-allocated)
473 */
474 /*----------------------------------------------------------------------------*/
475
476 void
477 cs_cdo_quantities_compute_b_wvf(const cs_cdo_connect_t *connect,
478 const cs_cdo_quantities_t *cdoq,
479 cs_lnum_t bf_id,
480 cs_real_t wvf[]);
481
482 /*----------------------------------------------------------------------------*/
483 /*!
484 * \brief Retrieve the face vector which the face_area * face_normal for a
485 * primal face (interior or border)
486 *
487 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
488 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
489 *
490 * \return a pointer to the face vector
491 */
492 /*----------------------------------------------------------------------------*/
493
494 inline static const cs_real_t *
cs_quant_get_face_vector_area(cs_lnum_t f_id,const cs_cdo_quantities_t * cdoq)495 cs_quant_get_face_vector_area(cs_lnum_t f_id,
496 const cs_cdo_quantities_t *cdoq)
497 {
498 if (f_id < cdoq->n_i_faces) /* Interior face */
499 return cdoq->i_face_normal + 3*f_id;
500 else /* Border face */
501 return cdoq->b_face_normal + 3*(f_id - cdoq->n_i_faces);
502 }
503
504 /*----------------------------------------------------------------------------*/
505 /*!
506 * \brief Retrieve the face center for a primal face (interior or border)
507 *
508 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
509 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
510 *
511 * \return a pointer to the face center coordinates
512 */
513 /*----------------------------------------------------------------------------*/
514
515 inline static const cs_real_t *
cs_quant_get_face_center(cs_lnum_t f_id,const cs_cdo_quantities_t * cdoq)516 cs_quant_get_face_center(cs_lnum_t f_id,
517 const cs_cdo_quantities_t *cdoq)
518 {
519 if (f_id < cdoq->n_i_faces) /* Interior face */
520 return cdoq->i_face_center + 3*f_id;
521 else /* Border face */
522 return cdoq->b_face_center + 3*(f_id - cdoq->n_i_faces);
523 }
524
525 /*----------------------------------------------------------------------------*/
526 /*!
527 * \brief Define a cs_quant_t structure for a primal face (interior or border)
528 *
529 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
530 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
531 *
532 * \return a initialize structure
533 */
534 /*----------------------------------------------------------------------------*/
535
536 cs_quant_t
537 cs_quant_set_face(cs_lnum_t f_id,
538 const cs_cdo_quantities_t *cdoq);
539
540 /*----------------------------------------------------------------------------*/
541 /*!
542 * \brief Retrieve the face surface and its unit normal vector for a primal
543 * face (interior or border)
544 *
545 * \param[in] f_id id related to the face (f_id > n_i_face -> border face)
546 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
547 *
548 * \return a pointer to the face normalized vector
549 */
550 /*----------------------------------------------------------------------------*/
551
552 cs_nvec3_t
553 cs_quant_set_face_nvec(cs_lnum_t f_id,
554 const cs_cdo_quantities_t *cdoq);
555
556 /*----------------------------------------------------------------------------*/
557 /*!
558 * \brief Get the normalized vector associated to a primal edge
559 *
560 * \param[in] e_id id related to an edge
561 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
562 *
563 * \return a pointer to the edge normalized vector
564 */
565 /*----------------------------------------------------------------------------*/
566
567 cs_nvec3_t
568 cs_quant_set_edge_nvec(cs_lnum_t e_id,
569 const cs_cdo_quantities_t *cdoq);
570
571 /*----------------------------------------------------------------------------*/
572 /*!
573 * \brief Get the two normalized vector associated to a dual edge
574 *
575 * \param[in] shift position in c2f_idx
576 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
577 *
578 * \return a pointer to the dual edge normalized vector
579 */
580 /*----------------------------------------------------------------------------*/
581
582 cs_nvec3_t
583 cs_quant_set_dedge_nvec(cs_lnum_t shift,
584 const cs_cdo_quantities_t *cdoq);
585
586 /*----------------------------------------------------------------------------*/
587 /*!
588 * \brief Dump a cs_quant_t structure
589 *
590 * \param[in] f FILE struct (stdout if NULL)
591 * \param[in] num entity number related to this quantity struct.
592 * \param[in] q cs_quant_t structure to dump
593 */
594 /*----------------------------------------------------------------------------*/
595
596 void
597 cs_quant_dump(FILE *f,
598 cs_lnum_t num,
599 const cs_quant_t q);
600
601 /*----------------------------------------------------------------------------*/
602
603 END_C_DECLS
604
605 #endif /* __CS_CDO_QUANTITIES_H__ */
606