1 #ifndef __CS_CDO_SCHEME_GEOMETRY_H__
2 #define __CS_CDO_SCHEME_GEOMETRY_H__
3
4 /*============================================================================
5 * Geometric computations for building discretization operators which is
6 * shared by several files
7 *============================================================================*/
8
9 /*
10 This file is part of Code_Saturne, a general-purpose CFD tool.
11
12 Copyright (C) 1998-2021 EDF S.A.
13
14 This program is free software; you can redistribute it and/or modify it under
15 the terms of the GNU General Public License as published by the Free Software
16 Foundation; either version 2 of the License, or (at your option) any later
17 version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
22 details.
23
24 You should have received a copy of the GNU General Public License along with
25 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
26 Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 */
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 /*----------------------------------------------------------------------------
34 * Local headers
35 *----------------------------------------------------------------------------*/
36
37 #include "cs_cdo_local.h"
38 #include "cs_math.h"
39
40 /*----------------------------------------------------------------------------*/
41
42 BEGIN_C_DECLS
43
44 /*============================================================================
45 * Macro definitions
46 *============================================================================*/
47
48 /*============================================================================
49 * Type definitions
50 *============================================================================*/
51
52 /*=============================================================================
53 * Inline static function prototypes
54 *============================================================================*/
55
56 /*----------------------------------------------------------------------------*/
57 /*!
58 * \brief Compute the value of the constant gradient of the Lagrange function
59 * attached to xc in p_{f,c} (constant inside this volume)
60 * Cellwise version.
61 *
62 * \param[in] f face number in the cellwise numbering to handle
63 * \param[in] cm pointer to a cell_mesh_t structure
64 * \param[in, out] grd_c gradient of the Lagrange function related to xc
65 */
66 /*----------------------------------------------------------------------------*/
67
68 inline static void
cs_compute_grdfc_cw(short int f,const cs_cell_mesh_t * cm,cs_real_t * grd_c)69 cs_compute_grdfc_cw(short int f,
70 const cs_cell_mesh_t *cm,
71 cs_real_t *grd_c)
72 {
73 /* Sanity checks */
74 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_HFQ));
75
76 const cs_real_t ohf = -cm->f_sgn[f]/cm->hfc[f];
77
78 grd_c[0] = ohf * cm->face[f].unitv[0];
79 grd_c[1] = ohf * cm->face[f].unitv[1];
80 grd_c[2] = ohf * cm->face[f].unitv[2];
81 }
82
83 /*----------------------------------------------------------------------------*/
84 /*!
85 * \brief Compute the value of the constant gradient of the Lagrange function
86 * attached to xc in p_{f,c} (constant inside this volume)
87 * Facewise version.
88 *
89 * \param[in] f face number in the cellwise numbering to handle
90 * \param[in] cm pointer to a cell_mesh_t structure
91 * \param[in, out] grd_c gradient of the Lagrange function related to xc
92 */
93 /*----------------------------------------------------------------------------*/
94
95 inline static void
cs_compute_grdfc_fw(const cs_face_mesh_t * fm,cs_real_t * grd_c)96 cs_compute_grdfc_fw(const cs_face_mesh_t *fm,
97 cs_real_t *grd_c)
98 {
99 const cs_real_t ohf = -fm->f_sgn/fm->hfc;
100
101 grd_c[0] = ohf * fm->face.unitv[0];
102 grd_c[1] = ohf * fm->face.unitv[1];
103 grd_c[2] = ohf * fm->face.unitv[2];
104 }
105
106 /*----------------------------------------------------------------------------*/
107 /*!
108 * \brief Compute the weight wef = |tef|/|f|
109 *
110 * \param[in] f id of the face in the cell-wise numbering
111 * \param[in] cm pointer to a cs_cell_mesh_t structure
112 * \param[in, out] wef pointer to an array storing the weight/vertex
113 */
114 /*----------------------------------------------------------------------------*/
115
116 inline static void
cs_compute_wef(short int f,const cs_cell_mesh_t * cm,cs_real_t * wef)117 cs_compute_wef(short int f,
118 const cs_cell_mesh_t *cm,
119 cs_real_t *wef)
120 {
121 /* Sanity checks */
122 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ));
123
124 const short int *f2e_idx = cm->f2e_idx + f;
125 const double *tef_vals = cm->tef + f2e_idx[0];
126 const double inv_f = 1./cm->face[f].meas;
127
128 /* Compute a weight for each edge of the current face */
129 for (short int e = 0; e < f2e_idx[1] - f2e_idx[0]; e++)
130 wef[e] = tef_vals[e] * inv_f;
131 }
132
133 /*----------------------------------------------------------------------------*/
134 /*!
135 * \brief Compute the volume related to each tetrahedron of base tef and apex
136 * x_c (there are as many tetrahedra as edges in a face)
137 *
138 * \param[in] f id of the face in the cell-wise numbering
139 * \param[in] cm pointer to a cs_cell_mesh_t structure
140 * \param[in, out] pefc pointer to an array storing the volumes
141 */
142 /*----------------------------------------------------------------------------*/
143
144 inline static void
cs_compute_pefc(short int f,const cs_cell_mesh_t * cm,cs_real_t * pefc)145 cs_compute_pefc(short int f,
146 const cs_cell_mesh_t *cm,
147 cs_real_t *pefc)
148 {
149 /* Sanity checks */
150 assert(cs_eflag_test(cm->flag,
151 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC));
152
153 const short int *f2e_idx = cm->f2e_idx + f;
154 const double *tef_vals = cm->tef + f2e_idx[0];
155 const double f_coef = cm->pvol_f[f]/cm->face[f].meas;
156
157 /* Compute a weight for each edge of the current face */
158 for (short int e = 0; e < f2e_idx[1] - f2e_idx[0]; e++)
159 pefc[e] = tef_vals[e] * f_coef;
160 }
161
162 /*----------------------------------------------------------------------------*/
163 /*!
164 * \brief Compute for a face the weight related to each vertex w_{v,f}
165 * This weight is equal to |dc(v) cap f|/|f| so that the sum of the
166 * weights is equal to 1.
167 *
168 * \param[in] f id of the face in the cell-wise numbering
169 * \param[in] cm pointer to a cs_cell_mesh_t structure
170 * \param[in, out] wvf pointer to an array storing the weight/vertex
171 * (allocated to the number of cell vertices)
172 */
173 /*----------------------------------------------------------------------------*/
174
175 inline static void
cs_compute_wvf(short int f,const cs_cell_mesh_t * cm,cs_real_t * wvf)176 cs_compute_wvf(short int f,
177 const cs_cell_mesh_t *cm,
178 cs_real_t *wvf)
179 {
180 /* Sanity checks */
181 assert(cs_eflag_test(cm->flag,
182 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
183
184 /* Reset weights */
185 memset(wvf, 0, cm->n_vc*sizeof(cs_real_t));
186
187 const short int *f2e_idx = cm->f2e_idx + f;
188 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
189 const double *tef_vals = cm->tef + f2e_idx[0];
190 const double inv_f = 1./cm->face[f].meas;
191
192 /* Compute a weight for each vertex of the current face */
193 for (short int e = 0; e < f2e_idx[1] - f2e_idx[0]; e++) {
194
195 const short int *v = cm->e2v_ids + 2*f2e_ids[e];
196 const double w = 0.5*tef_vals[e] * inv_f;
197 wvf[v[0]] += w;
198 wvf[v[1]] += w;
199
200 }
201 }
202
203 /*============================================================================
204 * Public function prototypes
205 *============================================================================*/
206
207 /*----------------------------------------------------------------------------*/
208 /*!
209 * \brief Compute the inertial matrix of a cell with respect to the point
210 * called "center". This computation is performed exactly thanks to
211 * quadrature based on a "tetrahedrization" of the cell.
212 *
213 * \param[in] cm pointer to a cs_cell_mesh_t structure
214 * \param[in] f id of the face in the cell numbering
215 * \param[in] ax main X-axis for the face-related coordinate system
216 * \param[in] ay main Y-axis for the face-related coordinate system
217 * \param[in] center coordinates of the face center
218 * \param[in, out] cov 2x2 symmetric covariance matrix to compute
219 */
220 /*----------------------------------------------------------------------------*/
221
222 void
223 cs_compute_face_covariance_tensor(const cs_cell_mesh_t *cm,
224 short int f,
225 const cs_nvec3_t ax,
226 const cs_nvec3_t ay,
227 const cs_real_t center[3],
228 cs_real_t cov[3]);
229
230 /*----------------------------------------------------------------------------*/
231 /*!
232 * \brief Compute the inertial matrix of a cell with respect to the point
233 * called "center". This computation is performed exactly thanks to
234 * quadrature based on a "tetrahedrization" of the cell.
235 *
236 * \param[in] cm pointer to a cs_cell_mesh_t structure
237 * \param[in] center coordinates of the cell center
238 * \param[in, out] inertia inertia matrix to compute
239 */
240 /*----------------------------------------------------------------------------*/
241
242 void
243 cs_compute_inertia_tensor(const cs_cell_mesh_t *cm,
244 const cs_real_t center[3],
245 cs_real_t inertia[3][3]);
246
247 /*----------------------------------------------------------------------------*/
248 /*!
249 * \brief Compute the gradient of a Lagrange function related to primal
250 * vertices in a p_{ef,c} subvolume of a cell c where e is an edge
251 * belonging to the face f with vertices v1 and v2
252 *
253 * \param[in] v1 number of the first vertex in cell numbering
254 * \param[in] v2 number of the second vertex in cell numbering
255 * \param[in] deq dual edge quantities
256 * \param[in] uvc xc --> xv unit tangent vector
257 * \param[in] lvc xc --> xv vector length
258 * \param[in, out] grd_v1 gradient of Lagrange function related to v1
259 * \param[in, out] grd_v2 gradient of Lagrange function related to v2
260 */
261 /*----------------------------------------------------------------------------*/
262
263 void
264 cs_compute_grd_ve(const short int v1,
265 const short int v2,
266 const cs_nvec3_t deq,
267 const cs_real_3_t uvc[],
268 const cs_real_t lvc[],
269 cs_real_t *grd_v1,
270 cs_real_t *grd_v2);
271
272 /*----------------------------------------------------------------------------*/
273 /*!
274 * \brief Compute for a face the weight related to each vertex w_{v,f} and
275 * the weight related to each edge
276 * w_{v,f} = |dc(v) cap f|/|f|
277 * Sum of w_{v,f} over the face vertices is equal to 1
278 * Sum of w_{e,f} over the face edges is equal to 1
279 *
280 * \param[in] f id of the face in the cell-wise numbering
281 * \param[in] cm pointer to a cs_cell_mesh_t structure
282 * \param[in, out] wvf weights of each face vertex
283 * \param[in, out] wef weights of each face edge
284 */
285 /*----------------------------------------------------------------------------*/
286
287 void
288 cs_compute_wef_wvf(short int f,
289 const cs_cell_mesh_t *cm,
290 cs_real_t *wvf,
291 cs_real_t *wef);
292
293 /*----------------------------------------------------------------------------*/
294
295 END_C_DECLS
296
297 #endif /* __CS_CDO_SCHEME_GEOMETRY_H__ */
298