1 /*============================================================================
2 * Build discrete stiffness matrices and handled boundary conditions for the
3 * diffusion term in CDO schemes
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <assert.h>
37 #include <math.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include <bft_mem.h>
44
45 #include "cs_hodge.h"
46 #include "cs_log.h"
47 #include "cs_math.h"
48 #include "cs_property.h"
49 #include "cs_reco.h"
50 #include "cs_scheme_geometry.h"
51
52 #if defined(DEBUG) && !defined(NDEBUG)
53 #include "cs_dbg.h"
54 #endif
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_cdo_diffusion.h"
61
62 /*----------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*=============================================================================
67 * Additional doxygen documentation
68 *============================================================================*/
69
70 /*!
71 \file cs_cdo_diffusion.c
72
73 \brief Build discrete stiffness matrices and handled boundary conditions for
74 diffusion term in CDO vertex-based and vertex+cell schemes.
75
76 */
77
78 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
79
80 /*=============================================================================
81 * Local Macro definitions
82 *============================================================================*/
83
84 #define CS_CDO_DIFFUSION_DBG 0
85
86 /* Redefined the name of functions from cs_math to get shorter names */
87 #define _dp3 cs_math_3_dot_product
88
89 /*=============================================================================
90 * Local structure definitions
91 *============================================================================*/
92
93 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
94
95 /*============================================================================
96 * Local variables
97 *============================================================================*/
98
99 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
100
101 /*============================================================================
102 * Private function prototypes
103 *============================================================================*/
104
105 /*----------------------------------------------------------------------------*/
106 /*!
107 * \brief Pre-compute the product between a property and the face vector areas
108 *
109 * \param[in] pty pointer to a \ref cs_property_data_t structure
110 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
111 * \param[in, out] kappa_f property against face vector for all faces
112 */
113 /*----------------------------------------------------------------------------*/
114
115 static inline void
_compute_kappa_f(const cs_property_data_t * pty,const cs_cell_mesh_t * cm,cs_real_3_t * kappa_f)116 _compute_kappa_f(const cs_property_data_t *pty,
117 const cs_cell_mesh_t *cm,
118 cs_real_3_t *kappa_f)
119 {
120 if (pty->is_unity) {
121 for (short int f = 0; f < cm->n_fc; f++) {
122 for (short int k = 0; k < 3; k++)
123 kappa_f[f][k] = cm->face[f].meas*cm->face[f].unitv[k];
124 }
125 }
126 else if (pty->is_iso) {
127 for (short int f = 0; f < cm->n_fc; f++) {
128 const cs_real_t coef = cm->face[f].meas*pty->value;
129 for (short int k = 0; k < 3; k++)
130 kappa_f[f][k] = coef * cm->face[f].unitv[k];
131 }
132 }
133 else {
134 for (short int f = 0; f < cm->n_fc; f++) {
135 cs_math_33_3_product(pty->tensor, cm->face[f].unitv, kappa_f[f]);
136 for (short int k = 0; k < 3; k++) kappa_f[f][k] *= cm->face[f].meas;
137 }
138 }
139 }
140
141 /*----------------------------------------------------------------------------*/
142 /*!
143 * \brief Compute \f$ \int_{fb} \nabla (u) \cdot \nu_{fb} v \f$ where \p fb
144 * is a boundary face (Co+St algorithm)
145 *
146 * \param[in] fb index of the boundary face on the local numbering
147 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
148 * \param[in] hodgep pointer to a \ref cs_hodge_param_t structure
149 * \param[in] kappa_f property against face vector for all faces
150 * \param[in, out] ntrgrd pointer to a local matrix structure
151 */
152 /*----------------------------------------------------------------------------*/
153
154 static void
_cdofb_normal_flux_reco(short int fb,const cs_cell_mesh_t * cm,const cs_hodge_param_t * hodgep,const cs_real_3_t * kappa_f,cs_sdm_t * ntrgrd)155 _cdofb_normal_flux_reco(short int fb,
156 const cs_cell_mesh_t *cm,
157 const cs_hodge_param_t *hodgep,
158 const cs_real_3_t *kappa_f,
159 cs_sdm_t *ntrgrd)
160 {
161 assert(hodgep->type == CS_HODGE_TYPE_EDFP);
162 assert(hodgep->algo == CS_HODGE_ALGO_COST);
163 assert(cs_eflag_test(cm->flag,
164 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFC));
165 assert(cm->f_sgn[fb] == 1); /* +1 because it's a boundary face */
166
167 const short int n_fc = cm->n_fc;
168 const cs_quant_t pfbq = cm->face[fb];
169 const cs_nvec3_t debq = cm->dedge[fb];
170 const double inv_volc = 1./cm->vol_c;
171
172 /* beta * |fb|^2 * nu_{fb}^T.kappa.nu_{fb} / |pvol_fb| */
173
174 const double stab_scaling =
175 hodgep->coef * pfbq.meas * _dp3(kappa_f[fb], pfbq.unitv) / cm->pvol_f[fb];
176
177 /* Get the 'fb' row */
178
179 cs_real_t *ntrgrd_fb = ntrgrd->val + fb * (n_fc + 1);
180 double row_sum = 0.0;
181
182 for (short int f = 0; f < n_fc; f++) {
183
184 const cs_quant_t pfq = cm->face[f];
185
186 /* consistent part */
187
188 const double consist_scaling = cm->f_sgn[f] * pfq.meas * inv_volc;
189 const double consist_part = consist_scaling * _dp3(kappa_f[fb], pfq.unitv);
190
191 /* stabilization part */
192
193 double stab_part = -consist_scaling*debq.meas*_dp3(debq.unitv, pfq.unitv);
194 if (f == fb) stab_part += 1;
195 stab_part *= stab_scaling;
196
197 const double fb_f_part = consist_part + stab_part;
198
199 ntrgrd_fb[f] -= fb_f_part; /* Minus because -du/dn */
200 row_sum += fb_f_part;
201
202 } /* Loop on f */
203
204 /* Cell column */
205
206 ntrgrd_fb[n_fc] += row_sum;
207 }
208
209 /*----------------------------------------------------------------------------*/
210 /*!
211 * \brief Compute the cellwise consistent gradient for Vb schemes
212 *
213 * \param[in] cm pointer to a cs_cell_mesh_t structure
214 * \param[in, out] grd_cell cellwise gradient vector attached to each vertex
215 * degree of freedom
216 */
217 /*----------------------------------------------------------------------------*/
218
219 static void
_svb_cellwise_grd(const cs_cell_mesh_t * cm,cs_real_3_t grd_cell[])220 _svb_cellwise_grd(const cs_cell_mesh_t *cm,
221 cs_real_3_t grd_cell[])
222 {
223 /* Reset array */
224
225 for (short int v = 0; v < cm->n_vc; v++)
226 grd_cell[v][0] = grd_cell[v][1] = grd_cell[v][2] = 0;
227
228 for (short int e = 0; e < cm->n_ec; e++) {
229
230 const short int *_v = cm->e2v_ids + 2*e; /* v1 = _v[0], v2 = _v[1] */
231 const short int sgn_v1 = cm->e2v_sgn[e]; /* sgn_v2 = - sgn_v1 */
232 const cs_nvec3_t dfq = cm->dface[e];
233
234 grd_cell[_v[0]][0] += sgn_v1 * dfq.meas * dfq.unitv[0];
235 grd_cell[_v[0]][1] += sgn_v1 * dfq.meas * dfq.unitv[1];
236 grd_cell[_v[0]][2] += sgn_v1 * dfq.meas * dfq.unitv[2];
237 grd_cell[_v[1]][0] -= sgn_v1 * dfq.meas * dfq.unitv[0];
238 grd_cell[_v[1]][1] -= sgn_v1 * dfq.meas * dfq.unitv[1];
239 grd_cell[_v[1]][2] -= sgn_v1 * dfq.meas * dfq.unitv[2];
240
241 } /* Loop on cell edges */
242
243 const cs_real_t over_vol_c = 1/cm->vol_c;
244 for (short int v = 0; v < cm->n_vc; v++) {
245 grd_cell[v][0] *= over_vol_c;
246 grd_cell[v][1] *= over_vol_c;
247 grd_cell[v][2] *= over_vol_c;
248 }
249 }
250
251 /*----------------------------------------------------------------------------*/
252 /*!
253 * \brief Compute the normal diffusive flux reconstruction for Vb schemes
254 * using the COST/VORONOI/BUBBLE algorithm. This normal flux is
255 * computed in each t_{ek,f}: the triangle with base e_k and apex xf
256 *
257 * \param[in] cm pointer to a cs_cell_mesh_t structure
258 * \param[in] dbeta d times the value of the stabilization coef.
259 * \param[in] ek edge id to handle (related to p_{ek,c} and t_{ek,f}
260 * \param[in] pekq edge quantities for ek
261 * \param[in] dfkq dual face quantities for ek
262 * \param[in] grd_cell cellwise gradient vector attached to each vertex
263 * degree of freedom
264 * \param[in] mnu diff. property times the face normal of t_{ek,f}
265 * \param[in, out] nflux reconstructed normal flux related to each vertex
266 */
267 /*----------------------------------------------------------------------------*/
268
269 static void
_nflux_reco_svb_ocs_in_pec(const cs_cell_mesh_t * cm,const double dbeta,const short int ek,const cs_quant_t pekq,const cs_nvec3_t dfkq,const cs_real_3_t grd_cell[],const cs_real_3_t mnu,cs_real_t nflux[])270 _nflux_reco_svb_ocs_in_pec(const cs_cell_mesh_t *cm,
271 const double dbeta,
272 const short int ek,
273 const cs_quant_t pekq,
274 const cs_nvec3_t dfkq,
275 const cs_real_3_t grd_cell[],
276 const cs_real_3_t mnu,
277 cs_real_t nflux[])
278 {
279 /* Reset the normal flux for each vertex of the cell */
280
281 for (short int v = 0; v < cm->n_vc; v++) nflux[v] = 0;
282
283 /* Compute the gradient reconstruction for the potential arrays of the
284 * canonical basis (size = number of cell vertices)
285 * L_{E_c}(GRAD(p)) on t_{e_k,f} for each vertex of the cell
286 * The restriction to t_{e_k,f} is identical to the restriction to p_{e_k,f}
287 *
288 * L_{E_c}(GRAD(p))|t_{e_k,f}
289 * = Gc(p) ──> Consist. part
290 * + beta*fd_ek(c)/(p_{ek,c} * (GRAD(p)|ek - ek.Gc(p) ) ──> Stabili. part
291 */
292
293 const cs_real_t stab_ek = dbeta/(_dp3(pekq.unitv, dfkq.unitv));
294 const short int *_vk = cm->e2v_ids + 2*ek; /* vk0 = _v[0], vk1 = _v[1] */
295 const short int sgn_vk0 = cm->e2v_sgn[ek]; /* sgn_vk1 = - sgn_vk0 */
296
297 /* Consistent + Stabilization part */
298
299 for (short int v = 0; v < cm->n_vc; v++) {
300
301 const cs_real_t ekgc = _dp3(pekq.unitv, grd_cell[v]);
302
303 /* Initialize the reconstruction of the gradient with the consistent part */
304
305 cs_real_3_t le_grd = {grd_cell[v][0], grd_cell[v][1], grd_cell[v][2]};
306
307 if (v == _vk[0]) { /* v = vk0 belonging to edge ek */
308
309 const cs_real_t stab_coef = stab_ek * (sgn_vk0/pekq.meas - ekgc);
310 for (int p = 0; p < 3; p++)
311 le_grd[p] += stab_coef * dfkq.unitv[p]; /* Stabilization part */
312
313 }
314 else if (v == _vk[1]) { /* v = vk1 belonging to edge ek */
315
316 const cs_real_t stab_coef = -stab_ek * (sgn_vk0/pekq.meas + ekgc);
317 for (int p = 0; p < 3; p++)
318 le_grd[p] += stab_coef * dfkq.unitv[p]; /* Stabilization part */
319
320 }
321 else { /* Other cell vertices */
322
323 const cs_real_t stab_coef = -stab_ek * ekgc;
324 for (int p = 0; p < 3; p++)
325 le_grd[p] += stab_coef * dfkq.unitv[p]; /* Stabilization part */
326
327 }
328
329 nflux[v] = -_dp3(mnu, le_grd); /* Minus because -du/dn */
330
331 } /* Loop on cell vertices */
332 }
333
334 /*----------------------------------------------------------------------------*/
335 /*!
336 * \brief Compute the normal trace operator for a given border face when a
337 * COST/BUBLLE or Voronoï algo. is used for reconstructing a flux
338 * from the degrees of freedom
339 *
340 * \param[in] fm pointer to a cs_face_mesh_t structure
341 * \param[in] cm pointer to a cs_cell_mesh_t structure
342 * \param[in] mnu property tensor times the face normal
343 * \param[in] dbeta d times the value of the stabilization coef.
344 * \param[in, out] cb pointer to a cell builder structure
345 * \param[in, out] ntrgrd local matrix related to the normal trace op. i.e.
346 * the flux operator
347 */
348 /*----------------------------------------------------------------------------*/
349
350 static void
_vb_ocs_normal_flux_op(const short int f,const cs_cell_mesh_t * cm,const cs_real_3_t mnu,double dbeta,cs_cell_builder_t * cb,cs_sdm_t * ntrgrd)351 _vb_ocs_normal_flux_op(const short int f,
352 const cs_cell_mesh_t *cm,
353 const cs_real_3_t mnu,
354 double dbeta,
355 cs_cell_builder_t *cb,
356 cs_sdm_t *ntrgrd)
357 {
358 cs_real_t *nflux = cb->values; /* size = cm->n_vc */
359 cs_real_3_t *grd_cell = cb->vectors; /* size = cm->n_vc */
360
361 /* Compute the cellwise-constant gradient for each array of the canonical
362 * basis of the potential DoFs (size = cm->n_vc)
363 */
364 _svb_cellwise_grd(cm, grd_cell);
365
366 /* Loop on border face edges */
367 for (int fe_idx = cm->f2e_idx[f]; fe_idx < cm->f2e_idx[f+1]; fe_idx++) {
368
369 const short int ek = cm->f2e_ids[fe_idx];
370 const cs_real_t tekf = cm->tef[fe_idx];
371 const cs_quant_t pekq = cm->edge[ek];
372 const cs_nvec3_t dfkq = cm->dface[ek];
373 const short int *_vk = cm->e2v_ids + 2*ek; /* v1 = _v[0], v2 = _v[1] */
374
375 /* Compute the reconstructed normal flux */
376 _nflux_reco_svb_ocs_in_pec(cm, dbeta, ek, pekq, dfkq,
377 (const cs_real_3_t *)grd_cell,
378 mnu,
379 nflux);
380
381 for (short int v = 0; v < cm->n_vc; v++) {
382
383 const double contrib_v = 0.5 * nflux[v] * tekf;
384
385 ntrgrd->val[_vk[0]*cm->n_vc + v] += contrib_v;
386 ntrgrd->val[_vk[1]*cm->n_vc + v] += contrib_v;
387
388 } /* Loop on cell vertices */
389
390 } /* Loop on face edges */
391
392 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 2
393 cs_log_printf(CS_LOG_DEFAULT,
394 ">> Flux.Op (NTRGRD) matrix (c_id: %d,f_id: %d)",
395 cm->c_id, cm->f_ids[f]);
396 cs_sdm_dump(cm->c_id, NULL, NULL, ntrgrd);
397 #endif
398 }
399
400 /*----------------------------------------------------------------------------*/
401 /*!
402 * \brief Compute the transposed of the normal trace operator and the full
403 * boundary operator to enforce weakly a Robin BC for a given border
404 * face when a COST algo. is used for reconstructing the degrees of
405 * freedom
406 *
407 * \param[in] fm pointer to a cs_face_mesh_t structure
408 * \param[in] cm pointer to a cs_cell_mesh_t structure
409 * \param[in] mnu property tensor times the face normal
410 * \param[in] beta value of the stabilization coef. related to reco.
411 * \param[in, out] cb pointer to a cell builder structure
412 * \param[in, out] bc_op local matrix storing du/dn.dv/dn
413 * \param[in, out] ntrgrd_tr local matrix related to the transposed version of
414 * the normal trace op. (i.e. the flux operator)
415 */
416 /*----------------------------------------------------------------------------*/
417
418 static void
_vb_cost_full_flux_op(const short int f,const cs_cell_mesh_t * cm,const cs_real_3_t mnu,double beta,cs_cell_builder_t * cb,cs_sdm_t * bc_op,cs_sdm_t * ntrgrd_tr)419 _vb_cost_full_flux_op(const short int f,
420 const cs_cell_mesh_t *cm,
421 const cs_real_3_t mnu,
422 double beta,
423 cs_cell_builder_t *cb,
424 cs_sdm_t *bc_op,
425 cs_sdm_t *ntrgrd_tr)
426 {
427 cs_real_t *nflux = cb->values; /* size = cm->n_vc */
428 cs_real_3_t *grd_cell = cb->vectors; /* size = cm->n_vc */
429
430 /* Initialize the local operators */
431 cs_sdm_square_init(cm->n_vc, bc_op);
432 cs_sdm_square_init(cm->n_vc, ntrgrd_tr);
433
434 /* Compute the cellwise-constant gradient for each array of the canonical
435 * basis of the potential DoFs (size = cm->n_vc)
436 */
437 _svb_cellwise_grd(cm, grd_cell);
438
439 /* Loop on border face edges */
440 for (int fe_idx = cm->f2e_idx[f]; fe_idx < cm->f2e_idx[f+1]; fe_idx++) {
441
442 const short int ek = cm->f2e_ids[fe_idx];
443 const cs_real_t tekf = cm->tef[fe_idx];
444 const cs_quant_t pekq = cm->edge[ek];
445 const cs_nvec3_t dfkq = cm->dface[ek];
446 const short int *_vk = cm->e2v_ids + 2*ek; /* v1 = _v[0], v2 = _v[1] */
447
448 /* Compute the reconstructed normal flux */
449 _nflux_reco_svb_ocs_in_pec(cm, beta, ek, pekq, dfkq,
450 (const cs_real_3_t *)grd_cell,
451 mnu,
452 nflux);
453
454 for (short int vi = 0; vi < cm->n_vc; vi++) {
455
456 const double contrib_vi = nflux[vi] * tekf;
457
458 ntrgrd_tr->val[vi*cm->n_vc + _vk[0]] += contrib_vi;
459 ntrgrd_tr->val[vi*cm->n_vc + _vk[1]] += contrib_vi;
460
461 for (short int vj = 0; vj < cm->n_vc; vj++) {
462 bc_op->val[vi*cm->n_vc + vj] += contrib_vi * nflux[vj];
463 } /* Loop on cells vertices (vj) */
464
465 } /* Loop on cell vertices (vi) */
466
467 } /* Loop on face edges */
468
469 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 2
470 cs_log_printf(CS_LOG_DEFAULT,
471 ">> Full flux.Op (NGRD.NGRD) matrix (c_id: %d,f_id: %d)",
472 cm->c_id, cm->f_ids[f]);
473 cs_log_printf(CS_LOG_DEFAULT, " grd_cell:\n>> ");
474 for (int k = 0; k < 3; k++) {
475 for (int i = 0; i < cm->n_vc; i++)
476 cs_log_printf(CS_LOG_DEFAULT, "% .4e ", grd_cell[i][k]);
477 if (k < 2)
478 cs_log_printf(CS_LOG_DEFAULT, "\n>> ");
479 else
480 cs_log_printf(CS_LOG_DEFAULT, "\n");
481 }
482 cs_sdm_dump(cm->c_id, NULL, NULL, bc_op);
483 cs_sdm_dump(cm->c_id, NULL, NULL, ntrgrd_tr);
484 #endif
485 }
486
487 /*----------------------------------------------------------------------------*/
488 /*!
489 * \brief Compute the normal trace operator for a given border face when a
490 * WBS algo. is used for reconstructing the degrees of freedom
491 *
492 * \param[in] fm pointer to a cs_face_mesh_t structure
493 * \param[in] cm pointer to a cs_cell_mesh_t structure
494 * \param[in] pty_nuf property tensor times the face normal
495 * \param[in, out] cb pointer to a cell builder structure
496 * \param[in, out] ntrgrd local matrix related to the normal trace op. i.e.
497 * the flux operator
498 */
499 /*----------------------------------------------------------------------------*/
500
501 static void
_vb_wbs_normal_flux_op(const cs_face_mesh_t * fm,const cs_cell_mesh_t * cm,const cs_real_3_t pty_nuf,cs_cell_builder_t * cb,cs_sdm_t * ntrgrd)502 _vb_wbs_normal_flux_op(const cs_face_mesh_t *fm,
503 const cs_cell_mesh_t *cm,
504 const cs_real_3_t pty_nuf,
505 cs_cell_builder_t *cb,
506 cs_sdm_t *ntrgrd)
507 {
508 cs_real_3_t grd_f, grd_v1, grd_v2, grd_c;
509
510 /* Useful quantities are stored in cb->values and cb->vectors */
511 cs_real_t *l_vc = cb->values;
512 cs_real_3_t *mng_ef = cb->vectors;
513 cs_real_3_t *u_vc = cb->vectors + fm->n_vf;
514
515 /* Initialize the local operator */
516 cs_sdm_square_init(cm->n_vc, ntrgrd);
517
518 /* Compute the gradient of the Lagrange function related to xc which is
519 constant inside p_{f,c} */
520 cs_compute_grdfc_fw(fm, grd_c);
521
522 /* f_coef = (pty_tensor*nu_f).grd_c */
523 const cs_real_t f_coef = fm->face.meas * _dp3(pty_nuf, grd_c);
524
525 /* Compute xc --> xv length and unit vector for all face vertices */
526 for (short int v = 0; v < fm->n_vf; v++)
527 cs_math_3_length_unitv(fm->xc, fm->xv + 3*v, l_vc + v, u_vc[v]);
528
529 /* Compute a weight for each vertex of the current face */
530 double sum_ef = 0.;
531 for (short int e = 0; e < fm->n_ef; e++) {
532
533 /* Gradient of the Lagrange function related to v1 and v2 */
534 cs_compute_grd_ve(fm->e2v_ids[2*e],
535 fm->e2v_ids[2*e+1],
536 fm->dedge,
537 (const cs_real_t (*)[3])u_vc, l_vc,
538 grd_v1, grd_v2);
539
540 /* Gradient of the Lagrange function related to a face.
541 This formula is a consequence of the Partition of the Unity */
542 for (int k = 0; k < 3; k++)
543 grd_f[k] = -(grd_c[k] + grd_v1[k] + grd_v2[k]);
544
545 const double tef_coef = fm->tef[e] * cs_math_1ov3;
546 mng_ef[e][0] = _dp3(pty_nuf, grd_v1) * tef_coef;
547 mng_ef[e][1] = _dp3(pty_nuf, grd_v2) * tef_coef;
548 mng_ef[e][2] = _dp3(pty_nuf, grd_f) * tef_coef;
549
550 sum_ef += mng_ef[e][2];
551
552 } /* End of loop on face edges */
553
554 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
555
556 const short int vi = fm->v_ids[vfi];
557 double *ntrgrd_vi = ntrgrd->val + vi*cm->n_vc;
558
559 /* Default contribution for this line */
560 const double default_coef = f_coef * fm->wvf[vfi];
561 for (short int vj = 0; vj < cm->n_vc; vj++)
562 ntrgrd_vi[vj] = default_coef * cm->wvc[vj]; /* two contributions */
563
564 /* Block Vf x Vf */
565 for (short int vfj = 0; vfj < fm->n_vf; vfj++) {
566
567 short int vj = fm->v_ids[vfj];
568 ntrgrd_vi[vj] += sum_ef * fm->wvf[vfi] * fm->wvf[vfj];
569
570 double entry_ij = 0.;
571 for (short int e = 0; e < fm->n_ef; e++) {
572
573 const short int v1 = fm->e2v_ids[2*e];
574 const short int v2 = fm->e2v_ids[2*e+1];
575
576 if (vfj == v1)
577 entry_ij += mng_ef[e][0] * fm->wvf[vfi];
578 else if (vfj == v2)
579 entry_ij += mng_ef[e][1] * fm->wvf[vfi];
580
581 if (vfi == v1 || vfi == v2) { /* i in e */
582 entry_ij += fm->wvf[vfj] * mng_ef[e][2];
583 if (vfj == v1) /* j is also in e */
584 entry_ij += mng_ef[e][0];
585 if (vfj == v2)
586 entry_ij += mng_ef[e][1];
587 }
588
589 } /* Loop on face edges */
590
591 ntrgrd_vi[vj] += entry_ij;
592
593 } /* Loop on face vertices (vj) */
594
595 } /* Loop on face vertices (vi) */
596
597 /* The flux is -1 times the reconstruction of the gradient */
598 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
599
600 const short int vi = fm->v_ids[vfi];
601 double *ntrgrd_vi = ntrgrd->val + vi*cm->n_vc;
602 for (short int vj = 0; vj < cm->n_vc; vj++)
603 ntrgrd_vi[vj] *= -1;
604
605 }
606
607 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 2
608 cs_log_printf(CS_LOG_DEFAULT,
609 ">> Flux.Op (NTRGRD) matrix (c_id: %d,f_id: %d)",
610 cm->c_id,cm->f_ids[fm->f_id]);
611 cs_sdm_dump(cm->c_id, NULL, NULL, ntrgrd);
612 #endif
613 }
614
615 /*----------------------------------------------------------------------------*/
616 /*!
617 * \brief Compute the normal trace operator for a given border face when a
618 * WBS algo. is used for reconstructing the degrees of freedom
619 * Specific to CDO-V+C schemes
620 *
621 * \param[in] fm pointer to a cs_face_mesh_t structure
622 * \param[in] cm pointer to a cs_cell_mesh_t structure
623 * \param[in] pty_nuf property tensor times the face normal
624 * \param[in, out] cb pointer to a cell builder structure
625 * \param[in, out] ntrgrd local matrix related to the normal trace op. i.e.
626 * the flux operator
627 */
628 /*----------------------------------------------------------------------------*/
629
630 static void
_vcb_wbs_normal_flux_op(const cs_face_mesh_t * fm,const cs_cell_mesh_t * cm,const cs_real_3_t pty_nuf,cs_cell_builder_t * cb,cs_sdm_t * ntrgrd)631 _vcb_wbs_normal_flux_op(const cs_face_mesh_t *fm,
632 const cs_cell_mesh_t *cm,
633 const cs_real_3_t pty_nuf,
634 cs_cell_builder_t *cb,
635 cs_sdm_t *ntrgrd)
636 {
637 cs_real_3_t grd_f, grd_v1, grd_v2, grd_c;
638
639 /* Useful quantities are stored in cb->values and cb->tmp-vect */
640 cs_real_t *l_vc = cb->values;
641 cs_real_3_t *mng_ef = cb->vectors;
642 cs_real_3_t *u_vc = cb->vectors + fm->n_vf;
643
644 /* Initialize the local operator */
645 cs_sdm_square_init(cm->n_vc + 1, ntrgrd);
646
647 /* Compute the gradient of the Lagrange function related to xc which is
648 constant inside p_{f,c} */
649 cs_compute_grdfc_fw(fm, grd_c);
650
651 const cs_real_t mng_cf = _dp3(pty_nuf, grd_c); /* (pty_tensor*nu_f).grd_c */
652
653 /* Compute xc --> xv length and unit vector for all face vertices */
654 for (short int v = 0; v < fm->n_vf; v++)
655 cs_math_3_length_unitv(fm->xc, fm->xv + 3*v, l_vc + v, u_vc[v]);
656
657 /* Compute a weight for each vertex of the current face */
658 for (short int e = 0; e < fm->n_ef; e++) {
659
660 const short int v1 = fm->e2v_ids[2*e];
661 const short int v2 = fm->e2v_ids[2*e+1];
662
663 /* Gradient of the Lagrange function related to v1 and v2 */
664 cs_compute_grd_ve(v1, v2, fm->dedge, (const cs_real_t (*)[3])u_vc, l_vc,
665 grd_v1, grd_v2);
666
667 /* Gradient of the Lagrange function related to a face.
668 This formula is a consequence of the Partition of the Unity */
669 for (int k = 0; k < 3; k++)
670 grd_f[k] = -(grd_c[k] + grd_v1[k] + grd_v2[k]);
671
672 const double tef_coef = fm->tef[e] * cs_math_1ov3;
673 mng_ef[e][0] = _dp3(pty_nuf, grd_v1) * tef_coef;
674 mng_ef[e][1] = _dp3(pty_nuf, grd_v2) * tef_coef;
675 mng_ef[e][2] = _dp3(pty_nuf, grd_f) * tef_coef;
676
677 } /* End of loop on face edges */
678
679 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
680
681 short int vi = fm->v_ids[vfi];
682 double *ntrgrd_vi = ntrgrd->val + vi*ntrgrd->n_rows;
683
684 /* Contribution to the cell column */
685 ntrgrd_vi[cm->n_vc] = fm->wvf[vfi] * fm->face.meas * mng_cf;
686
687 /* Block Vf x Vf */
688 for (short int vfj = 0; vfj < fm->n_vf; vfj++) {
689
690 double entry_ij = 0.;
691 for (short int e = 0; e < fm->n_ef; e++) {
692
693 const short int v1 = fm->e2v_ids[2*e];
694 const short int v2 = fm->e2v_ids[2*e+1];
695
696 double coef_i = fm->wvf[vfi];
697 if (vfi == v1 || vfi == v2)
698 coef_i += 1;
699
700 double coef_j = fm->wvf[vfj] * mng_ef[e][2];
701 if (vfj == v1)
702 coef_j += mng_ef[e][0];
703 else if (vfj == v2)
704 coef_j += mng_ef[e][1];
705
706 entry_ij += coef_i * coef_j;
707
708 } /* Loop on face edges */
709
710 ntrgrd_vi[fm->v_ids[vfj]] += entry_ij;
711
712 } /* Loop on face vertices (vj) */
713
714 } /* Loop on face vertices (vi) */
715
716 /* The flux is -1 times the reconstruction of the gradient */
717 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
718
719 const short int vi = fm->v_ids[vfi];
720 double *ntrgrd_vi = ntrgrd->val + vi*ntrgrd->n_rows;
721 for (short int vj = 0; vj < ntrgrd->n_rows; vj++)
722 ntrgrd_vi[vj] *= -1;
723
724 }
725
726 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 2
727 cs_log_printf(CS_LOG_DEFAULT,
728 ">> Flux.Op (NTRGRD) matrix (c_id: %d,f_id: %d)",
729 cm->c_id,cm->f_ids[fm->f_id]);
730 cs_sdm_dump(cm->c_id, NULL, NULL, ntrgrd);
731 #endif
732 }
733
734 /*----------------------------------------------------------------------------*/
735 /*!
736 * \brief Weak enforcement of Dirichlet BCs using a Nitsche technique.
737 * Scalar-valued case for CDO-VB or CDO-VCb schemes
738 *
739 * \param[in] pcoef value of the penalization coefficient
740 * \param[in] fm pointer to a cs_face_mesh_t structure
741 * \param[in, out] ntrgrd pointer to a local matrix structure
742 * \param[in, out] csys structure storing the cell-wise system
743 */
744 /*----------------------------------------------------------------------------*/
745
746 static inline void
_svb_nitsche(const double pcoef,const cs_face_mesh_t * fm,cs_sdm_t * ntrgrd,cs_cell_sys_t * csys)747 _svb_nitsche(const double pcoef,
748 const cs_face_mesh_t *fm,
749 cs_sdm_t *ntrgrd,
750 cs_cell_sys_t *csys)
751 {
752 /* Sanity checks */
753 assert(pcoef > 0);
754 assert(csys->mat->n_rows == ntrgrd->n_rows);
755
756 /* Update local RHS and system */
757 for (short int v = 0; v < fm->n_vf; v++) {
758
759 /* Set the penalty diagonal coefficient */
760 const double pcoef_v = pcoef * fm->wvf[v];
761 const short int vi = fm->v_ids[v];
762
763 ntrgrd->val[vi*(1 + ntrgrd->n_rows)] += pcoef_v;
764 csys->rhs[vi] += pcoef_v * csys->dir_values[vi];
765
766 } /* Dirichlet or homogeneous Dirichlet */
767
768 }
769
770 #if 0 /* Lack of robustness w.r.t the linear algebra */
771 /*----------------------------------------------------------------------------*/
772 /*!
773 * \brief Using the local (cellwise) "normal trace gradient" matrix takes
774 * into account Dirichlet BCs by a weak enforcement using Nitsche
775 * technique.
776 * Case of WBS algorithm either for CDO-VB or CDO-VCb schemes
777 *
778 * \param[in] pcoef value of the penalization coefficient
779 * \param[in] fm pointer to a cs_face_mesh_t structure
780 * \param[in, out] ntrgrd pointer to a local matrix structure
781 * \param[in, out] cb pointer to a cs_cell_builder_t structure
782 * \param[in, out] csys structure storing the cell-wise system
783 */
784 /*----------------------------------------------------------------------------*/
785
786 static void
787 _wbs_nitsche(const double pcoef,
788 const cs_face_mesh_t *fm,
789 cs_sdm_t *ntrgrd,
790 cs_cell_builder_t *cb,
791 cs_cell_sys_t *csys)
792 {
793 /* Sanity checks */
794 assert(pcoef > 0);
795 assert(csys->mat->n_rows == ntrgrd->n_rows);
796
797 const short int n_csys = csys->mat->n_rows;
798
799 /* Build the related border Hodge operator */
800 cs_sdm_t *hloc = cb->aux;
801
802 cs_hodge_compute_wbs_surfacic(fm, hloc); /* hloc is of size n_vf */
803
804 /* Add the border Hodge op. to the normal trace op.
805 Update RHS whith H*p^{dir} */
806 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
807
808 const double *hi = hloc->val + vfi*fm->n_vf;
809 const short int vi = fm->v_ids[vfi];
810
811 double *ntrg_vi = ntrgrd->val + vi*n_csys;
812
813 for (short int vfj = 0; vfj < fm->n_vf; vfj++) {
814
815 const double pcoef_ij = pcoef * hi[vfj];
816 const short int vj = fm->v_ids[vfj];
817
818 ntrg_vi[vj] += pcoef_ij;
819 csys->rhs[vi] += pcoef_ij * csys->dir_values[vj];
820
821 } /* Loop on face vertices vj */
822 } /* Loop on face vertices vi */
823
824 }
825 #endif
826
827 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
828
829 /*============================================================================
830 * Public function prototypes
831 *============================================================================*/
832
833 /*----------------------------------------------------------------------------*/
834 /*!
835 * \brief Take into account Dirichlet BCs by a weak enforcement by a
836 * penalization technique with a huge value
837 * Predefined prototype to match the function pointer
838 * cs_cdo_enforce_bc_t
839 *
840 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
841 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
842 * \param[in, out] fm pointer to a cs_face_mesh_t structure
843 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
844 * \param[in, out] cb pointer to a cs_cell_builder_t structure
845 * \param[in, out] csys structure storing the cell-wise system
846 */
847 /*----------------------------------------------------------------------------*/
848
849 void
cs_cdo_diffusion_pena_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)850 cs_cdo_diffusion_pena_dirichlet(const cs_equation_param_t *eqp,
851 const cs_cell_mesh_t *cm,
852 cs_face_mesh_t *fm,
853 cs_hodge_t *hodge,
854 cs_cell_builder_t *cb,
855 cs_cell_sys_t *csys)
856 {
857 /* Prototype common to cs_cdo_enforce_bc_t. Hence the unused parameters */
858 CS_UNUSED(hodge);
859 CS_UNUSED(fm);
860 CS_UNUSED(cm);
861 CS_UNUSED(cb);
862 assert(csys != NULL); /* Sanity checks */
863
864 /* Enforcement of the Dirichlet BCs */
865 if (csys->has_dirichlet == false)
866 return; /* Nothing to do */
867
868 /* Penalize diagonal entry (and its rhs if needed) */
869 for (short int i = 0; i < csys->n_dofs; i++) {
870
871 if (csys->dof_flag[i] & CS_CDO_BC_HMG_DIRICHLET) {
872 csys->mat->val[i + csys->n_dofs*i] += eqp->strong_pena_bc_coeff;
873 }
874 else if (csys->dof_flag[i] & CS_CDO_BC_DIRICHLET) {
875 csys->mat->val[i + csys->n_dofs*i] += eqp->strong_pena_bc_coeff;
876 csys->rhs[i] += csys->dir_values[i] * eqp->strong_pena_bc_coeff;
877 }
878
879 } /* Loop on degrees of freedom */
880
881 }
882
883 /*----------------------------------------------------------------------------*/
884 /*!
885 * \brief Take into account Dirichlet BCs by a weak enforcement by a
886 * penalization technique with a huge value.
887 * Case of a cellwise system defined by block.
888 * Predefined prototype to match the function pointer
889 * cs_cdo_enforce_bc_t
890 *
891 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
892 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
893 * \param[in, out] fm pointer to a cs_face_mesh_t structure
894 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
895 * \param[in, out] cb pointer to a cs_cell_builder_t structure
896 * \param[in, out] csys structure storing the cell-wise system
897 */
898 /*----------------------------------------------------------------------------*/
899
900 void
cs_cdo_diffusion_pena_block_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)901 cs_cdo_diffusion_pena_block_dirichlet(const cs_equation_param_t *eqp,
902 const cs_cell_mesh_t *cm,
903 cs_face_mesh_t *fm,
904 cs_hodge_t *hodge,
905 cs_cell_builder_t *cb,
906 cs_cell_sys_t *csys)
907 {
908 /* Prototype common to cs_cdo_enforce_bc_t. Hence the unused parameters */
909 CS_UNUSED(hodge);
910 CS_UNUSED(fm);
911 CS_UNUSED(cm);
912 CS_UNUSED(cb);
913 assert(csys != NULL); /* Sanity checks */
914
915 /* Enforcement of the Dirichlet BCs */
916 if (csys->has_dirichlet == false)
917 return; /* Nothing to do */
918
919 cs_sdm_t *m = csys->mat;
920 cs_sdm_block_t *bd = m->block_desc;
921 assert(bd != NULL);
922
923 /* Penalize diagonal entry (and its rhs if needed) */
924 int shift = 0;
925 for (short int bi = 0; bi < bd->n_row_blocks; bi++) {
926
927 cs_sdm_t *mII = cs_sdm_get_block(m, bi, bi);
928 assert(mII->n_rows == mII->n_cols);
929 cs_real_t *_rhs = csys->rhs + shift;
930 const cs_flag_t *_flag = csys->dof_flag + shift;
931 const cs_real_t *_dir_val = csys->dir_values + shift;
932
933 for (int i = 0; i < mII->n_rows; i++) {
934
935 if (_flag[i] & CS_CDO_BC_HMG_DIRICHLET) {
936 mII->val[i + mII->n_rows*i] += eqp->strong_pena_bc_coeff;
937 }
938 else if (_flag[i] & CS_CDO_BC_DIRICHLET) {
939 mII->val[i + mII->n_rows*i] += eqp->strong_pena_bc_coeff;
940 _rhs[i] += _dir_val[i] * eqp->strong_pena_bc_coeff;
941 }
942
943 }
944
945 shift += mII->n_rows;
946
947 } /* Block bi */
948
949 }
950
951 /*----------------------------------------------------------------------------*/
952 /*!
953 * \brief Take into account Dirichlet BCs by keeping the DoFs related to
954 * Dirichlet BCs in the algebraic system (i.e. a weak enforcement)
955 * The corresponding DoFs are algebraically "removed" of the system
956 *
957 * | | | | | | | | | |
958 * | Aii | Aid | | Aii | 0 | |bi| |bi-Aid.bd |
959 * |------------| --> |------------| and |--| --> |----------|
960 * | | | | | | | | | |
961 * | Adi | Add | | 0 | Id | |bd| | xd |
962 *
963 * where xd is the value of the Dirichlet BC
964 * Predefined prototype to match the function pointer
965 * cs_cdo_enforce_bc_t
966 *
967 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
968 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
969 * \param[in, out] fm pointer to a cs_face_mesh_t structure
970 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
971 * \param[in, out] cb pointer to a cs_cell_builder_t structure
972 * \param[in, out] csys structure storing the cell-wise system
973 */
974 /*----------------------------------------------------------------------------*/
975
976 void
cs_cdo_diffusion_alge_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)977 cs_cdo_diffusion_alge_dirichlet(const cs_equation_param_t *eqp,
978 const cs_cell_mesh_t *cm,
979 cs_face_mesh_t *fm,
980 cs_hodge_t *hodge,
981 cs_cell_builder_t *cb,
982 cs_cell_sys_t *csys)
983 {
984 /* Prototype common to cs_cdo_enforce_bc_t Hence the unused parameters */
985 CS_UNUSED(eqp);
986 CS_UNUSED(fm);
987 CS_UNUSED(cm);
988 CS_UNUSED(hodge);
989 assert(csys != NULL); /* Sanity checks */
990
991 /* Enforcement of the Dirichlet BCs */
992 if (csys->has_dirichlet == false)
993 return; /* Nothing to do */
994
995 double *x_dir = cb->values;
996 double *ax_dir = cb->values + csys->n_dofs;
997
998 memset(cb->values, 0, 2*csys->n_dofs*sizeof(double));
999
1000 /* Build x_dir */
1001 for (short int i = 0; i < csys->n_dofs; i++)
1002 if (csys->dof_flag[i] & CS_CDO_BC_DIRICHLET) /* Only non-homogeneous */
1003 x_dir[i] = csys->dir_values[i];
1004
1005 /* Contribution of the Dirichlet conditions */
1006 cs_sdm_matvec(csys->mat, x_dir, ax_dir);
1007
1008 /* Second pass: Replace the Dirichlet block by a diagonal block */
1009 for (short int i = 0; i < csys->n_dofs; i++) {
1010
1011 if (cs_cdo_bc_is_dirichlet(csys->dof_flag[i])) { /* All Dirichlet:
1012 homogeneous or not */
1013 /* Reset row i */
1014 memset(csys->mat->val + csys->n_dofs*i, 0, csys->n_dofs*sizeof(double));
1015 /* Reset column i */
1016 for (short int j = 0; j < csys->n_dofs; j++)
1017 csys->mat->val[i + csys->n_dofs*j] = 0;
1018 csys->mat->val[i*(1 + csys->n_dofs)] = 1;
1019
1020 /* Set the RHS */
1021 csys->rhs[i] = csys->dir_values[i];
1022
1023 } /* DoF associated to a Dirichlet BC */
1024 else
1025 csys->rhs[i] -= ax_dir[i]; /* Update RHS */
1026
1027 } /* Loop on degrees of freedom */
1028
1029 }
1030
1031 /*----------------------------------------------------------------------------*/
1032 /*!
1033 * \brief Take into account Dirichlet BCs by keeping the DoFs related to
1034 * Dirichlet BCs in the algebraic system (i.e. a weak enforcement)
1035 * The corresponding DoFs are algebraically "removed" of the system
1036 * Block version.
1037 *
1038 * | | | | | | | | | |
1039 * | Aii | Aid | | Aii | 0 | |bi| |bi-Aid.xd |
1040 * |------------| --> |------------| and |--| --> |----------|
1041 * | | | | | | | | | |
1042 * | Adi | Add | | 0 | Id | |bd| | xd |
1043 *
1044 * where xd is the value of the Dirichlet BC
1045 * Predefined prototype to match the function pointer
1046 * cs_cdo_enforce_bc_t
1047 *
1048 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1049 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1050 * \param[in, out] fm pointer to a cs_face_mesh_t structure
1051 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1052 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1053 * \param[in, out] csys structure storing the cell-wise system
1054 */
1055 /*----------------------------------------------------------------------------*/
1056
1057 void
cs_cdo_diffusion_alge_block_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1058 cs_cdo_diffusion_alge_block_dirichlet(const cs_equation_param_t *eqp,
1059 const cs_cell_mesh_t *cm,
1060 cs_face_mesh_t *fm,
1061 cs_hodge_t *hodge,
1062 cs_cell_builder_t *cb,
1063 cs_cell_sys_t *csys)
1064 {
1065 /* Prototype common to cs_cdo_enforce_bc_t. Hence the unused parameters */
1066 CS_UNUSED(eqp);
1067 CS_UNUSED(fm);
1068 CS_UNUSED(cm);
1069 CS_UNUSED(hodge);
1070 assert(csys != NULL); /* Sanity checks */
1071
1072 /* Enforcement of the Dirichlet BCs */
1073 if (csys->has_dirichlet == false)
1074 return; /* Nothing to do */
1075
1076 double *x_dir = cb->values;
1077 double *ax_dir = cb->values + csys->n_dofs;
1078 cs_sdm_t *m = csys->mat;
1079 cs_sdm_block_t *bd = m->block_desc;
1080 assert(bd != NULL);
1081
1082 memset(cb->values, 0, 2*csys->n_dofs*sizeof(double));
1083
1084 /* Build x_dir */
1085 for (short int i = 0; i < csys->n_dofs; i++)
1086 if (csys->dof_flag[i] & CS_CDO_BC_DIRICHLET) /* Only non-homogeneous */
1087 x_dir[i] = csys->dir_values[i];
1088
1089 /* Contribution of the Dirichlet conditions */
1090 cs_sdm_block_matvec(csys->mat, x_dir, ax_dir);
1091
1092 /* Second pass: Replace the Dirichlet block by a diagonal block */
1093 int s = 0;
1094 for (int bi = 0; bi < bd->n_row_blocks; bi++) {
1095
1096 cs_real_t *_rhs = csys->rhs + s;
1097 cs_real_t *_dir = csys->dir_values + s;
1098 cs_flag_t *_flg = csys->dof_flag + s;
1099 cs_sdm_t *mII = cs_sdm_get_block(m, bi, bi);
1100 assert(mII->n_rows == mII->n_cols);
1101
1102 /* Is the current block associated to a Dirichlet BC ? */
1103 int n_dir = 0;
1104 for (int i = 0; i < mII->n_rows; i++)
1105 if (cs_cdo_bc_is_dirichlet(_flg[i]))
1106 n_dir += 1;
1107
1108 if (n_dir > 0) {
1109
1110 if (n_dir != mII->n_rows)
1111 bft_error(__FILE__, __LINE__, 0,
1112 "%s: Partial Dirichlet block not yet implemented", __func__);
1113
1114 /* Reset block row bi and block colum bi */
1115 for (int bj = 0; bj < bd->n_col_blocks; bj++) {
1116
1117 if (bj != bi) {
1118
1119 cs_sdm_t *mIJ = cs_sdm_get_block(m, bi, bj);
1120 cs_sdm_t *mJI = cs_sdm_get_block(m, bj, bi);
1121
1122 memset(mIJ->val, 0, mIJ->n_rows*mIJ->n_cols*sizeof(double));
1123 memset(mJI->val, 0, mJI->n_rows*mJI->n_cols*sizeof(double));
1124
1125 }
1126 else { /* mII block --> diagonal */
1127
1128 memset(mII->val, 0, mII->n_rows*mII->n_rows*sizeof(double));
1129
1130 for (int i = 0; i < mII->n_rows; i++) {
1131 mII->val[i + mII->n_rows*i] = 1;
1132 _rhs[i] = _dir[i]; /* Set the RHS */
1133 }
1134
1135 }
1136
1137 } /* Loop on row/columnn blocks */
1138
1139 } /* DoF associated to a Dirichlet BC */
1140 else {
1141
1142 for (int i = 0; i < mII->n_rows; i++)
1143 _rhs[i] -= ax_dir[s+i]; /* Update RHS */
1144
1145 } /* Not a Dirichlet block */
1146
1147 s += mII->n_rows;
1148
1149 } /* Block bi */
1150
1151 }
1152
1153 /*----------------------------------------------------------------------------*/
1154 /*!
1155 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
1156 * technique.
1157 * Case of scalar-valued CDO Face-based schemes
1158 * Predefined prototype to match the function pointer
1159 * cs_cdo_enforce_bc_t
1160 *
1161 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1162 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1163 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1164 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1165 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1166 * \param[in, out] csys structure storing the cellwise system
1167 */
1168 /*----------------------------------------------------------------------------*/
1169
1170 void
cs_cdo_diffusion_sfb_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1171 cs_cdo_diffusion_sfb_weak_dirichlet(const cs_equation_param_t *eqp,
1172 const cs_cell_mesh_t *cm,
1173 cs_face_mesh_t *fm,
1174 cs_hodge_t *hodge,
1175 cs_cell_builder_t *cb,
1176 cs_cell_sys_t *csys)
1177 {
1178 CS_UNUSED(fm);
1179 assert(csys != NULL); /* Sanity checks */
1180
1181 /* Enforcement of the Dirichlet BCs */
1182
1183 if (csys->has_dirichlet == false)
1184 return; /* Nothing to do */
1185
1186 assert(cm != NULL && cb != NULL);
1187
1188 const cs_property_data_t *pdata = hodge->pty_data;
1189 const double chi =
1190 eqp->weak_pena_bc_coeff * fabs(pdata->eigen_ratio)*pdata->eigen_max;
1191
1192 /* First step: pre-compute the product between diffusion property and the
1193 face vector areas */
1194
1195 cs_real_3_t *kappa_f = cb->vectors;
1196 _compute_kappa_f(pdata, cm, kappa_f);
1197
1198 /* Initialize the matrix related to this flux reconstruction operator */
1199
1200 const short int n_dofs = cm->n_fc + 1;
1201 cs_sdm_t *bc_op = cb->loc;
1202 cs_sdm_square_init(n_dofs, bc_op);
1203
1204 /* First pass: build the bc_op matrix */
1205
1206 for (short int i = 0; i < csys->n_bc_faces; i++) {
1207
1208 /* Get the boundary face in the cell numbering */
1209
1210 const short int f = csys->_f_ids[i];
1211
1212 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1213
1214 /* Compute \int_f du/dn v and update the matrix */
1215
1216 _cdofb_normal_flux_reco(f, cm, hodge->param,
1217 (const cs_real_t (*)[3])kappa_f,
1218 bc_op);
1219
1220 } /* If Dirichlet */
1221
1222 } /* Loop boundary faces */
1223
1224 /* Second pass: Update the cell system with the bc_op matrix and the Dirichlet
1225 values. Avoid a truncation error if the arbitrary coefficient of the
1226 Nitsche algorithm is large
1227 */
1228
1229 for (short int i = 0; i < csys->n_bc_faces; i++) {
1230
1231 /* Get the boundary face in the cell numbering */
1232 const short int f = csys->_f_ids[i];
1233
1234 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1235
1236 /* chi * \meas{f} / h_f */
1237
1238 const cs_real_t pcoef = chi * sqrt(cm->face[f].meas);
1239
1240 /* Diagonal term */
1241
1242 bc_op->val[f*(n_dofs + 1)] += pcoef;
1243
1244 /* rhs */
1245
1246 csys->rhs[f] += pcoef * csys->dir_values[f];
1247
1248 } /* If Dirichlet */
1249
1250 } /* Loop on boundary faces */
1251
1252 /* Update the local system matrix */
1253
1254 cs_sdm_add(csys->mat, bc_op);
1255 }
1256
1257 /*----------------------------------------------------------------------------*/
1258 /*!
1259 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
1260 * technique.
1261 * Case of vector-valued CDO Face-based schemes
1262 * The idea is to compute the scalar version and dispatch it three
1263 * times, one for each Cartesian components
1264 * Predefined prototype to match the function pointer
1265 * cs_cdo_enforce_bc_t
1266 *
1267 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1268 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1269 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1270 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1271 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1272 * \param[in, out] csys structure storing the cellwise system
1273 */
1274 /*----------------------------------------------------------------------------*/
1275
1276 void
cs_cdo_diffusion_vfb_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1277 cs_cdo_diffusion_vfb_weak_dirichlet(const cs_equation_param_t *eqp,
1278 const cs_cell_mesh_t *cm,
1279 cs_face_mesh_t *fm,
1280 cs_hodge_t *hodge,
1281 cs_cell_builder_t *cb,
1282 cs_cell_sys_t *csys)
1283 {
1284 CS_UNUSED(fm);
1285 assert(csys != NULL); /* Sanity checks */
1286
1287 /* Enforcement of the Dirichlet BCs */
1288 if (csys->has_dirichlet == false)
1289 return; /* Nothing to do */
1290
1291 assert(cm != NULL && cb != NULL && csys != NULL && hodge);
1292
1293 const cs_property_data_t *pty = hodge->pty_data;
1294 const double chi =
1295 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
1296
1297 /* First step: pre-compute the product between diffusion property and the
1298 face vector areas */
1299 cs_real_3_t *kappa_f = cb->vectors;
1300 _compute_kappa_f(pty, cm, kappa_f);
1301
1302 /* Initialize the matrix related this flux reconstruction operator */
1303 const short int n_dofs = cm->n_fc + 1; /* n_blocks or n_scalar_dofs */
1304 cs_sdm_t *bc_op = cb->loc;
1305 cs_sdm_square_init(n_dofs, bc_op);
1306
1307 /* First pass: build the bc_op matrix */
1308 for (short int i = 0; i < csys->n_bc_faces; i++) {
1309
1310 /* Get the boundary face in the cell numbering */
1311 const short int f = csys->_f_ids[i];
1312
1313 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1314
1315 /* Compute \int_f du/dn v and update the matrix */
1316 _cdofb_normal_flux_reco(f, cm, hodge->param,
1317 (const cs_real_t (*)[3])kappa_f,
1318 bc_op);
1319
1320 } /* If Dirichlet */
1321
1322 } /* Loop boundary faces */
1323
1324 /* Second pass: Update the cell system with the bc_op matrix and the Dirichlet
1325 * values. Avoid a truncation error if the arbitrary coefficient of the
1326 * Nitsche algorithm is large
1327 */
1328 for (short int i = 0; i < csys->n_bc_faces; i++) {
1329
1330 /* Get the boundary face in the cell numbering */
1331 const short int f = csys->_f_ids[i];
1332
1333 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1334
1335 /* chi * \meas{f} / h_f */
1336 const cs_real_t pcoef = chi * sqrt(cm->face[f].meas);
1337
1338 /* Diagonal term */
1339 bc_op->val[f*(n_dofs + 1)] += pcoef;
1340
1341 /* rhs */
1342 for (short int k = 0; k < 3; k++)
1343 csys->rhs[3*f + k] += pcoef * csys->dir_values[3*f + k];
1344
1345 } /* If Dirichlet */
1346
1347 } /* Loop on boundary faces */
1348
1349 assert(pty->is_iso == true); /* if not the case something else TODO ? */
1350
1351 /* Update the local system matrix */
1352 for (int bi = 0; bi < n_dofs; bi++) { /* n_(scalar)_dofs == n_blocks */
1353 for (int bj = 0; bj < n_dofs; bj++) {
1354
1355 /* Retrieve the 3x3 matrix */
1356 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
1357 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1358
1359 const cs_real_t _val = bc_op->val[n_dofs*bi + bj];
1360 /* Update diagonal terms only */
1361 bij->val[0] += _val;
1362 bij->val[4] += _val;
1363 bij->val[8] += _val;
1364
1365 }
1366 }
1367
1368 }
1369
1370 /*----------------------------------------------------------------------------*/
1371 /*!
1372 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
1373 * technique plus a symmetric treatment
1374 * Case of scalar-valued CDO Face-based schemes
1375 * Predefined prototype to match the function pointer
1376 * cs_cdo_enforce_bc_t
1377 *
1378 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1379 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1380 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1381 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1382 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1383 * \param[in, out] csys structure storing the cellwise system
1384 */
1385 /*----------------------------------------------------------------------------*/
1386
1387 void
cs_cdo_diffusion_sfb_wsym_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1388 cs_cdo_diffusion_sfb_wsym_dirichlet(const cs_equation_param_t *eqp,
1389 const cs_cell_mesh_t *cm,
1390 cs_face_mesh_t *fm,
1391 cs_hodge_t *hodge,
1392 cs_cell_builder_t *cb,
1393 cs_cell_sys_t *csys)
1394 {
1395 CS_UNUSED(fm);
1396 assert(csys != NULL); /* Sanity checks */
1397
1398 /* Enforcement of the Dirichlet BCs */
1399 if (csys->has_dirichlet == false)
1400 return; /* Nothing to do */
1401
1402 assert(cm != NULL && cb != NULL && csys != NULL);
1403
1404 const cs_property_data_t *pty = hodge->pty_data;
1405 const double chi =
1406 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
1407
1408 /* First step: pre-compute the product between diffusion property and the
1409 face vector areas */
1410 cs_real_3_t *kappa_f = cb->vectors;
1411 _compute_kappa_f(pty, cm, kappa_f);
1412
1413 const short int n_dofs = cm->n_fc + 1, n_f = cm->n_fc;
1414 cs_sdm_t *bc_op = cb->loc, *bc_op_t = cb->aux;
1415 cs_sdm_square_init(n_dofs, bc_op);
1416
1417 /* First pass: build the bc_op matrix */
1418 for (short int i = 0; i < csys->n_bc_faces; i++) {
1419
1420 /* Get the boundary face in the cell numbering */
1421 const short int f = csys->_f_ids[i];
1422
1423 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1424
1425 /* Compute \int_f du/dn v and update the matrix */
1426 _cdofb_normal_flux_reco(f, cm, hodge->param,
1427 (const cs_real_t (*)[3])kappa_f,
1428 bc_op);
1429
1430 } /* If Dirichlet */
1431
1432 } /* Loop boundary faces */
1433
1434 /* Symmetrize the Nitsche contribution */
1435 cs_real_t *dir_val = cb->values, *u0_trgradv = cb->values + n_dofs;
1436
1437 /* Putting the face DoFs of the BC, into a face- and cell-DoFs array */
1438 memcpy(dir_val, csys->dir_values, n_f*sizeof(cs_real_t));
1439 dir_val[n_f] = 0.;
1440
1441 /* Update bc_op = bc_op + transp and transp = transpose(bc_op) cb->loc
1442 plays the role of the flux operator */
1443 cs_sdm_square_add_transpose(bc_op, bc_op_t);
1444 cs_sdm_square_matvec(bc_op_t, dir_val, u0_trgradv);
1445
1446 /* Second pass: Update the cell system with the bc_op matrix and the Dirichlet
1447 values. Avoid a truncation error if the arbitrary coefficient of the
1448 Nitsche algorithm is large
1449 */
1450 for (short int i = 0; i < n_dofs; i++) /* Cell too! */
1451 csys->rhs[i] += u0_trgradv[i];
1452
1453 for (short int i = 0; i < csys->n_bc_faces; i++) {
1454
1455 /* Get the boundary face in the cell numbering */
1456 const short int f = csys->_f_ids[i];
1457
1458 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1459
1460 /* chi * \meas{f} / h_f */
1461 const cs_real_t pcoef = chi * sqrt(cm->face[f].meas);
1462
1463 /* Diagonal term */
1464 bc_op->val[f*(n_dofs + 1)] += pcoef;
1465
1466 /* rhs */
1467 csys->rhs[f] += pcoef * csys->dir_values[f];
1468
1469 } /* If Dirichlet */
1470
1471 } /* Loop on boundary faces */
1472
1473 /* Update the local system matrix */
1474 cs_sdm_add(csys->mat, bc_op);
1475
1476 }
1477
1478 /*----------------------------------------------------------------------------*/
1479 /*!
1480 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
1481 * technique plus a symmetric treatment.
1482 * Case of vector-valued CDO Face-based schemes
1483 * The idea is to compute the scalar version and dispatch it three
1484 * times, one for each Cartesian components
1485 * Predefined prototype to match the function pointer
1486 * cs_cdo_enforce_bc_t
1487 *
1488 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1489 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1490 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1491 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1492 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1493 * \param[in, out] csys structure storing the cellwise system
1494 */
1495 /*----------------------------------------------------------------------------*/
1496
1497 void
cs_cdo_diffusion_vfb_wsym_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1498 cs_cdo_diffusion_vfb_wsym_dirichlet(const cs_equation_param_t *eqp,
1499 const cs_cell_mesh_t *cm,
1500 cs_face_mesh_t *fm,
1501 cs_hodge_t *hodge,
1502 cs_cell_builder_t *cb,
1503 cs_cell_sys_t *csys)
1504 {
1505 CS_UNUSED(fm);
1506 assert(csys != NULL); /* Sanity checks */
1507
1508 /* Enforcement of the Dirichlet BCs */
1509 if (csys->has_dirichlet == false)
1510 return; /* Nothing to do */
1511
1512 assert(cm != NULL && cb != NULL && csys != NULL && hodge != NULL);
1513
1514 const cs_property_data_t *pty = hodge->pty_data;
1515 const double chi =
1516 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
1517
1518 /* First step: pre-compute the product between diffusion property and the
1519 face vector areas */
1520 cs_real_3_t *kappa_f = cb->vectors;
1521 _compute_kappa_f(pty, cm, kappa_f);
1522
1523 /* Initialize the matrix related this flux reconstruction operator */
1524 const short int n_dofs = cm->n_fc + 1; /* n_blocks or n_scalar_dofs */
1525 cs_sdm_t *bc_op = cb->loc, *bc_op_t = cb->aux;
1526 cs_sdm_square_init(n_dofs, bc_op);
1527
1528 /* First pass: build the bc_op matrix */
1529 for (short int i = 0; i < csys->n_bc_faces; i++) {
1530
1531 /* Get the boundary face in the cell numbering */
1532 const short int f = csys->_f_ids[i];
1533
1534 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1535
1536 /* Compute \int_f du/dn v and update the matrix */
1537 _cdofb_normal_flux_reco(f, cm, hodge->param,
1538 (const cs_real_t (*)[3])kappa_f,
1539 bc_op);
1540
1541 } /* If Dirichlet */
1542
1543 } /* Loop boundary faces */
1544
1545 /* Second pass: add the bc_op matrix, add the BC */
1546
1547 /* Update bc_op = bc_op + transp and bc_op_t = transpose(bc_op) */
1548 cs_sdm_square_add_transpose(bc_op, bc_op_t);
1549
1550 /* Fill a face- and cell-DoFs array with the face DoFs values associated to
1551 * Dirichlet BC. These Dirichlet BC values are interlaced. One first
1552 * de-interlaces this array to perform the local matrix-vector product:
1553 * bc_op_t dir_val = u0_trgradv */
1554 cs_real_t *dir_val = cb->values, *u0_trgradv = cb->values + n_dofs;
1555
1556 dir_val[cm->n_fc] = 0.; /* cell DoF is not attached to a Dirichlet */
1557
1558 for (int k = 0; k < 3; k++) {
1559
1560 /* Build dir_val */
1561 for (short int f = 0; f < cm->n_fc; f++)
1562 dir_val[f] = csys->dir_values[3*f+k];
1563
1564 cs_sdm_square_matvec(bc_op_t, dir_val, u0_trgradv);
1565
1566 for (short int i = 0; i < n_dofs; i++) /* Cell too! */
1567 csys->rhs[3*i+k] += u0_trgradv[i];
1568
1569 } /* Loop on components */
1570
1571 /* Third pass: Update the cell system with the bc_op matrix and the Dirichlet
1572 values. Avoid a truncation error if the arbitrary coefficient of the
1573 Nitsche algorithm is large
1574 */
1575 for (short int i = 0; i < csys->n_bc_faces; i++) {
1576
1577 /* Get the boundary face in the cell numbering */
1578 const short int f = csys->_f_ids[i];
1579
1580 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
1581
1582 /* chi * \meas{f} / h_f */
1583 const cs_real_t pcoef = chi * sqrt(cm->face[f].meas);
1584
1585 /* Diagonal term */
1586 bc_op->val[f*(n_dofs + 1)] += pcoef;
1587
1588 /* rhs */
1589 for (short int k = 0; k < 3; k++)
1590 csys->rhs[3*f + k] += pcoef * csys->dir_values[3*f + k];
1591
1592 } /* If Dirichlet */
1593
1594 } /* Loop on boundary faces */
1595
1596 assert(pty->is_iso == true); /* if not the case something else TODO ? */
1597
1598 /* Update the local system matrix */
1599 for (int bi = 0; bi < n_dofs; bi++) { /* n_(scalar)_dofs == n_blocks */
1600 for (int bj = 0; bj < n_dofs; bj++) {
1601
1602 /* Retrieve the 3x3 matrix */
1603 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
1604 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1605
1606 const cs_real_t _val = bc_op->val[n_dofs*bi + bj];
1607 /* Update diagonal terms only */
1608 bij->val[0] += _val;
1609 bij->val[4] += _val;
1610 bij->val[8] += _val;
1611
1612 }
1613 }
1614
1615 }
1616
1617 /*----------------------------------------------------------------------------*/
1618 /*!
1619 * \brief Take into account sliding BCs by a weak enforcement using Nitsche
1620 * technique plus a symmetric treatment.
1621 * Case of vector-valued CDO Face-based schemes
1622 * Predefined prototype to match the function pointer
1623 * cs_cdo_enforce_bc_t
1624 *
1625 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1626 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1627 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1628 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1629 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1630 * \param[in, out] csys structure storing the cellwise system
1631 */
1632 /*----------------------------------------------------------------------------*/
1633
1634 void
cs_cdo_diffusion_vfb_wsym_sliding(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1635 cs_cdo_diffusion_vfb_wsym_sliding(const cs_equation_param_t *eqp,
1636 const cs_cell_mesh_t *cm,
1637 cs_face_mesh_t *fm,
1638 cs_hodge_t *hodge,
1639 cs_cell_builder_t *cb,
1640 cs_cell_sys_t *csys)
1641 {
1642 CS_UNUSED(fm);
1643 assert(csys != NULL);
1644
1645 /* Enforcement of the Dirichlet BCs */
1646 if (csys->has_sliding == false)
1647 return; /* Nothing to do */
1648
1649 const cs_property_data_t *pty = hodge->pty_data;
1650
1651 /* Sanity checks */
1652 assert(cm != NULL && cb != NULL);
1653 assert(pty->is_iso == true); /* if not the case something else TODO ? */
1654
1655 const double chi =
1656 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
1657 const short int n_f = cm->n_fc;
1658 const short int n_dofs = n_f + 1; /* n_blocks or n_scalar_dofs */
1659
1660 /* First step: pre-compute the product between diffusion property and the
1661 face vector areas */
1662 cs_real_3_t *kappa_f = cb->vectors;
1663 _compute_kappa_f(pty, cm, kappa_f);
1664
1665 /* Initialize the matrix related this flux reconstruction operator */
1666 cs_sdm_t *bc_op = cb->loc;
1667 cs_sdm_square_init(n_dofs, bc_op);
1668
1669 /* First pass: build the bc_op matrix */
1670 for (short int i = 0; i < csys->n_bc_faces; i++) {
1671
1672 /* Get the boundary face in the cell numbering */
1673 const short int f = csys->_f_ids[i];
1674
1675 if (csys->bf_flag[f] & CS_CDO_BC_SLIDING) {
1676
1677 /* Compute \int_f du/dn v and update the matrix */
1678 _cdofb_normal_flux_reco(f, cm, hodge->param,
1679 (const cs_real_t (*)[3])kappa_f,
1680 bc_op);
1681
1682 } /* If sliding */
1683
1684 } /* Loop boundary faces */
1685
1686 /* Second pass: add the bc_op matrix, add the BC */
1687 for (short int i = 0; i < csys->n_bc_faces; i++) {
1688
1689 /* Get the boundary face in the cell numbering */
1690 const short int fi = csys->_f_ids[i];
1691
1692 if (csys->bf_flag[fi] & CS_CDO_BC_SLIDING) {
1693
1694 const cs_quant_t pfq = cm->face[fi];
1695 const cs_real_t *ni = pfq.unitv;
1696 const cs_real_t ni_ni[9] =
1697 { ni[0]*ni[0], ni[0]*ni[1], ni[0]*ni[2],
1698 ni[1]*ni[0], ni[1]*ni[1], ni[1]*ni[2],
1699 ni[2]*ni[0], ni[2]*ni[1], ni[2]*ni[2]};
1700
1701 /* chi * \meas{f} / h_f */
1702 const cs_real_t pcoef = chi * sqrt(pfq.meas);
1703
1704 for (short int xj = 0; xj < n_dofs; xj++) {
1705 /* It should be done both for face- and cell-defined DoFs */
1706
1707 if ( xj == fi ) {
1708 /* Retrieve the 3x3 matrix */
1709 cs_sdm_t *bii = cs_sdm_get_block(csys->mat, fi, fi);
1710 assert(bii->n_rows == bii->n_cols && bii->n_rows == 3);
1711
1712 const cs_real_t _val = pcoef + 2 * bc_op->val[n_dofs*fi + fi];
1713
1714 for (short int k = 0; k < 9; k++)
1715 bii->val[k] += ni_ni[k] * _val;
1716
1717 }
1718 else { /* xj != fi */
1719
1720 /* Retrieve the 3x3 matrix */
1721 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, fi, xj);
1722 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1723 cs_sdm_t *bji = cs_sdm_get_block(csys->mat, xj, fi);
1724 assert(bji->n_rows == bji->n_cols && bji->n_rows == 3);
1725
1726 const cs_real_t _val = bc_op->val[n_dofs*fi + xj];
1727
1728 for (short int k = 0; k < 9; k++) {
1729 bij->val[k] += ni_ni[k] * _val;
1730 /* ni_ni is symmetric */
1731 bji->val[k] += ni_ni[k] * _val;
1732 }
1733
1734 } /* End if */
1735
1736 } /* Loop on xj */
1737
1738 } /* If sliding */
1739
1740 } /* Loop boundary faces */
1741
1742 }
1743
1744 /*----------------------------------------------------------------------------*/
1745 /*!
1746 * \brief Take into account Robin BCs.
1747 * Case of scalar-valued CDO-Fb schemes with a CO+ST algorithm.
1748 * Predefined prototype to match the function pointer
1749 * cs_cdo_enforce_bc_t
1750 *
1751 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1752 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1753 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1754 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1755 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1756 * \param[in, out] csys structure storing the cellwise system
1757 */
1758 /*----------------------------------------------------------------------------*/
1759
1760 void
cs_cdo_diffusion_sfb_cost_robin(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1761 cs_cdo_diffusion_sfb_cost_robin(const cs_equation_param_t *eqp,
1762 const cs_cell_mesh_t *cm,
1763 cs_face_mesh_t *fm,
1764 cs_hodge_t *hodge,
1765 cs_cell_builder_t *cb,
1766 cs_cell_sys_t *csys)
1767 {
1768 CS_UNUSED(eqp);
1769 CS_UNUSED(fm);
1770 CS_UNUSED(hodge);
1771 CS_UNUSED(cb);
1772
1773 /* Sanity checks */
1774 assert(cm != NULL && csys != NULL);
1775
1776 /* Enforcement of the Robin BCs */
1777 if (csys->has_robin == false)
1778 return; /* Nothing to do */
1779
1780 /* Robin BC expression: K du/dn + alpha*(u - u0) = g */
1781 for (short int i = 0; i < csys->n_bc_faces; i++) {
1782
1783 /* Get the boundary face in the cell numbering */
1784 const short int f = csys->_f_ids[i];
1785
1786 if (csys->bf_flag[f] & CS_CDO_BC_ROBIN) {
1787
1788 /* Robin BC expression: K du/dn + alpha*(u - u0) = g */
1789 /* ------------------------------------------------- */
1790
1791 const double alpha = csys->rob_values[3*f];
1792 const double u0 = csys->rob_values[3*f+1];
1793 const double g = csys->rob_values[3*f+2];
1794 const cs_real_t f_meas = cm->face[f].meas;
1795
1796 /* Update the RHS and the local system */
1797 csys->rhs[f] += (alpha*u0 + g)*f_meas;
1798
1799 cs_real_t *row_f_val = csys->mat->val + f*csys->n_dofs;
1800 row_f_val[f] += alpha*f_meas;
1801
1802 } /* Robin face */
1803
1804 } /* Loop on boundary faces */
1805
1806 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
1807 if (cs_dbg_cw_test(eqp, cm, csys))
1808 cs_cell_sys_dump(">> Cell %d, scalar Fb: After Robin", csys);
1809 #endif
1810 }
1811
1812 /*----------------------------------------------------------------------------*/
1813 /*!
1814 * \brief Take into account Robin BCs.
1815 * Case of scalar-valued CDO-Vb schemes with a CO+ST algorithm.
1816 * Predefined prototype to match the function pointer
1817 * cs_cdo_enforce_bc_t
1818 *
1819 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1820 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1821 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1822 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1823 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1824 * \param[in, out] csys structure storing the cellwise system
1825 */
1826 /*----------------------------------------------------------------------------*/
1827
1828 void
cs_cdo_diffusion_svb_cost_robin(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1829 cs_cdo_diffusion_svb_cost_robin(const cs_equation_param_t *eqp,
1830 const cs_cell_mesh_t *cm,
1831 cs_face_mesh_t *fm,
1832 cs_hodge_t *hodge,
1833 cs_cell_builder_t *cb,
1834 cs_cell_sys_t *csys)
1835 {
1836 CS_UNUSED(eqp);
1837 CS_UNUSED(hodge);
1838
1839 /* Sanity checks */
1840 assert(cm != NULL && cb != NULL && csys != NULL);
1841
1842 /* Enforcement of the Robin BCs */
1843 if (csys->has_robin == false)
1844 return; /* Nothing to do */
1845
1846 /* Robin BC expression: K du/dn + alpha*(u - u0) = beta */
1847 cs_sdm_t *bc_op = cb->loc;
1848
1849 /* Reset local operator */
1850 cs_sdm_square_init(cm->n_vc, bc_op);
1851
1852 for (short int i = 0; i < csys->n_bc_faces; i++) {
1853
1854 /* Get the boundary face in the cell numbering */
1855 const short int f = csys->_f_ids[i];
1856
1857 if (csys->bf_flag[f] & CS_CDO_BC_ROBIN) {
1858
1859 /* Compute the face-view of the mesh */
1860 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
1861
1862 /* Robin BC expression: K du/dn + alpha*(u - u0) = beta */
1863 /* ---------------------------------------------------- */
1864
1865 const double alpha = csys->rob_values[3*f];
1866 const double u0 = csys->rob_values[3*f+1];
1867 const double beta = csys->rob_values[3*f+2];
1868 const double g = beta + alpha*u0;
1869
1870 /* Update the RHS and the local system */
1871 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
1872
1873 const double f_cap_dualv = fm->face.meas * fm->wvf[vfi];
1874 const short int vi = fm->v_ids[vfi];
1875
1876 csys->rhs[vi] += g*f_cap_dualv;
1877 bc_op->val[vi*(1+cm->n_vc)] += alpha*f_cap_dualv;
1878
1879 }
1880
1881 } /* Robin face */
1882 } /* Loop on boundary faces */
1883
1884 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
1885 if (cs_dbg_cw_test(eqp, cm, csys)) {
1886 cs_log_printf(CS_LOG_DEFAULT, ">> Cell %d Vb COST Robin bc matrix",
1887 cm->c_id);
1888 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, bc_op);
1889 }
1890 #endif
1891
1892 /* Add contribution to the linear system */
1893 cs_sdm_add(csys->mat, bc_op);
1894 }
1895
1896 /*----------------------------------------------------------------------------*/
1897 /*!
1898 * \brief Take into account generic BCs by a weak enforcement using Nitsche
1899 * technique. According to the settings one can apply Neumann BCs if
1900 * alpha = 0, Dirichlet BCs if alpha >> 1 or Robin BCs
1901 * Case of scalar-valued CDO-Vb schemes with a CO+ST algorithm.
1902 * Predefined prototype to match the function pointer
1903 * cs_cdo_enforce_bc_t
1904 *
1905 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
1906 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1907 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1908 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
1909 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1910 * \param[in, out] csys structure storing the cellwise system
1911 */
1912 /*----------------------------------------------------------------------------*/
1913
1914 void
cs_cdo_diffusion_svb_cost_generic(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)1915 cs_cdo_diffusion_svb_cost_generic(const cs_equation_param_t *eqp,
1916 const cs_cell_mesh_t *cm,
1917 cs_face_mesh_t *fm,
1918 cs_hodge_t *hodge,
1919 cs_cell_builder_t *cb,
1920 cs_cell_sys_t *csys)
1921 {
1922 CS_UNUSED(eqp);
1923
1924 /* Sanity checks */
1925 assert(csys != NULL);
1926
1927 /* Enforcement of the Robin BCs */
1928 if (csys->has_robin == false)
1929 return; /* Nothing to do */
1930
1931 assert(cm != NULL && cb != NULL && hodge != NULL);
1932
1933 const double cs_gamma_generic = 1e-2;
1934 const cs_hodge_param_t *hodgep = hodge->param;
1935 const cs_property_data_t *pty = hodge->pty_data;
1936
1937 assert(pty->need_tensor);
1938
1939 /* Robin BC expression: K du/dn + alpha*(u - u0) = g
1940 * Store x = u0 + g/alpha
1941 */
1942 cs_real_t *x = cb->values;
1943 cs_real_t *ax = cb->values + cm->n_vc;
1944 cs_sdm_t *bc_op = cb->loc;
1945 cs_sdm_t *ntrgrd_tr = cb->aux;
1946
1947 for (short int i = 0; i < csys->n_bc_faces; i++) {
1948
1949 /* Get the boundary face in the cell numbering */
1950 const short int f = csys->_f_ids[i];
1951
1952 if (csys->bf_flag[f] & CS_CDO_BC_ROBIN) {
1953
1954 /* Compute the face-view of the mesh */
1955 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
1956
1957 /* Compute the product: matpty*face unit normal */
1958 cs_real_3_t pty_nuf;
1959 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
1960
1961 /* Compute the flux operator related to the trace on the current face
1962 of the normal gradient */
1963 _vb_cost_full_flux_op(f, cm, pty_nuf, hodgep->coef, cb, bc_op, ntrgrd_tr);
1964
1965 /* Robin BC expression: K du/dn + alpha*(u - u0) = g */
1966 /* ------------------------------------------------- */
1967 const double alpha = csys->rob_values[3*f];
1968 const double u0 = csys->rob_values[3*f+1];
1969 const double g = csys->rob_values[3*f+2];
1970 assert(fabs(alpha) > 0);
1971 const double epsilon = 1./alpha;
1972
1973 const double hf = sqrt(cm->face[f].meas);
1974 const double pena_coef = 1./(epsilon + cs_gamma_generic*hf);
1975 const double flux_coef = cs_gamma_generic*hf*pena_coef;
1976
1977 memset(x, 0, sizeof(cs_real_t)*cm->n_vc);
1978 for (short int v = 0; v < fm->n_vf; v++) {
1979 const short int vi = fm->v_ids[v];
1980 x[vi] = u0 + epsilon*g;
1981 }
1982
1983 /* Update the RHS */
1984 cs_sdm_square_matvec(ntrgrd_tr, x, ax);
1985 for (short int v = 0; v < fm->n_vf; v++) {
1986 const double pcoef_v = pena_coef * fm->face.meas * fm->wvf[v];
1987 const short int vi = fm->v_ids[v];
1988 csys->rhs[vi] += pcoef_v*x[vi] - flux_coef*ax[vi];
1989 }
1990
1991 /* Update the local system matrix:
1992 *
1993 * Initially bc_op = du/dn.dv/dn
1994 * 1) bc_op *= -epsilon*gamma*hf/(epsilon + gamma*hf)
1995 * 2) bc_op += (u,v) * 1/(epsilon + gamma*hf)
1996 * 3) bc_op += (-gamma*hf)/(epsilon + gamma*hf) [ntrgrd + ntrgrd_tr]
1997 */
1998
1999 /* Step 1 */
2000 for (int ij = 0; ij < cm->n_vc*cm->n_vc; ij++)
2001 bc_op->val[ij] *= -epsilon*flux_coef;
2002
2003 /* Step 2 */
2004 for (short int v = 0; v < fm->n_vf; v++) {
2005 const double pcoef_v = pena_coef * fm->face.meas * fm->wvf[v];
2006 const short int vi = fm->v_ids[v];
2007 bc_op->val[vi*(1 + bc_op->n_rows)] += pcoef_v;
2008 }
2009
2010 /* Step 3 */
2011 cs_sdm_square_2symm(ntrgrd_tr); /* ntrgrd_tr is now a symmetric matrix */
2012 cs_sdm_add_mult(bc_op, -flux_coef, ntrgrd_tr);
2013
2014 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2015 if (cs_dbg_cw_test(eqp, cm, csys)) {
2016 cs_log_printf(CS_LOG_DEFAULT,
2017 ">> Cell Vb COST generic bc matrix (f_id: %d)", fm->f_id);
2018 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, bc_op);
2019 }
2020 #endif
2021
2022 /* Add contribution to the linear system */
2023 cs_sdm_add(csys->mat, bc_op);
2024
2025 } /* Robin face */
2026 } /* Loop on boundary faces */
2027
2028 }
2029
2030 /*----------------------------------------------------------------------------*/
2031 /*!
2032 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2033 * technique. Case of scalar-valued CDO-Vb schemes with an orthogonal
2034 * splitting between the consistency/stabilization parts (OCS)
2035 * Predefined prototype to match the function pointer
2036 * cs_cdo_enforce_bc_t
2037 *
2038 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2039 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2040 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2041 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2042 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2043 * \param[in, out] csys structure storing the cellwise system
2044 */
2045 /*----------------------------------------------------------------------------*/
2046
2047 void
cs_cdo_diffusion_svb_ocs_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2048 cs_cdo_diffusion_svb_ocs_weak_dirichlet(const cs_equation_param_t *eqp,
2049 const cs_cell_mesh_t *cm,
2050 cs_face_mesh_t *fm,
2051 cs_hodge_t *hodge,
2052 cs_cell_builder_t *cb,
2053 cs_cell_sys_t *csys)
2054 {
2055 assert(csys != NULL); /* Sanity checks */
2056
2057 /* Enforcement of the Dirichlet BCs */
2058 if (csys->has_dirichlet == false)
2059 return; /* Nothing to do */
2060
2061 assert(cm != NULL && cb != NULL && hodge != NULL);
2062
2063 const cs_property_data_t *pty = hodge->pty_data;
2064 const cs_hodge_param_t *hodgep = hodge->param;
2065 const double chi =
2066 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2067 const cs_real_t dbeta =
2068 (hodgep->algo == CS_HODGE_ALGO_BUBBLE) ? hodgep->coef : 3*hodgep->coef;
2069
2070 assert(pty->need_tensor);
2071
2072 cs_sdm_t *ntrgrd = cb->loc;
2073
2074 /* Initialize the local operator */
2075 cs_sdm_square_init(cm->n_vc, ntrgrd);
2076
2077 for (short int i = 0; i < csys->n_bc_faces; i++) {
2078
2079 /* Get the boundary face in the cell numbering */
2080 const short int f = csys->_f_ids[i];
2081
2082 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2083
2084 /* Compute the face-view of the mesh */
2085 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2086
2087 /* Compute the product: matpty*face unit normal */
2088 cs_real_3_t pty_nuf;
2089 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2090
2091 /* Compute the flux operator related to the trace on the current face
2092 of the normal gradient */
2093 _vb_ocs_normal_flux_op(f, cm, pty_nuf, dbeta, cb, ntrgrd);
2094
2095 /* Update the RHS and the local system matrix */
2096 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2097
2098 } /* Dirichlet face */
2099 } /* Loop on boundary faces */
2100
2101 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2102 if (cs_dbg_cw_test(eqp, cm, csys)) {
2103 cs_log_printf(CS_LOG_DEFAULT, ">> Cell Vb.COST Weak bc matrix");
2104 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2105 }
2106 #endif
2107
2108 /* Add contribution to the linear system */
2109 cs_sdm_add(csys->mat, ntrgrd);
2110 }
2111
2112 /*----------------------------------------------------------------------------*/
2113 /*!
2114 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2115 * technique. A Dirichlet is set for the three components of the
2116 * vector. Case of vector-valued CDO-Vb schemes with an orthogonal
2117 * splitting between the consistency/stabilization parts (OCS)
2118 * Predefined prototype to match the function pointer
2119 * cs_cdo_enforce_bc_t
2120 *
2121 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2122 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2123 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2124 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2125 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2126 * \param[in, out] csys structure storing the cellwise system
2127 */
2128 /*----------------------------------------------------------------------------*/
2129
2130 void
cs_cdo_diffusion_vvb_ocs_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2131 cs_cdo_diffusion_vvb_ocs_weak_dirichlet(const cs_equation_param_t *eqp,
2132 const cs_cell_mesh_t *cm,
2133 cs_face_mesh_t *fm,
2134 cs_hodge_t *hodge,
2135 cs_cell_builder_t *cb,
2136 cs_cell_sys_t *csys)
2137 {
2138 assert(csys != NULL);
2139
2140 /* Enforcement of the Dirichlet BCs */
2141 if (csys->has_dirichlet == false)
2142 return; /* Nothing to do */
2143
2144 assert(cm != NULL && cb != NULL && hodge != NULL);
2145
2146 const cs_property_data_t *pty = hodge->pty_data;
2147 const double chi =
2148 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2149 const cs_hodge_param_t *hodgep = hodge->param;
2150 const cs_real_t dbeta =
2151 (hodgep->algo == CS_HODGE_ALGO_BUBBLE) ? hodgep->coef : 3*hodgep->coef;
2152
2153 assert(pty->need_tensor);
2154 cs_sdm_t *ntrgrd = cb->loc;
2155
2156 /* Initialize the local operator */
2157 cs_sdm_square_init(cm->n_vc, ntrgrd);
2158
2159 for (short int i = 0; i < csys->n_bc_faces; i++) {
2160
2161 /* Get the boundary face in the cell numbering */
2162 const short int f = csys->_f_ids[i];
2163
2164 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2165
2166 /* Compute the face-view of the mesh */
2167 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2168
2169 /* Compute the product: matpty*face unit normal */
2170 cs_real_3_t pty_nuf;
2171 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2172
2173 /* Compute the flux operator related to the trace on the current face
2174 of the normal gradient */
2175 _vb_ocs_normal_flux_op(f, cm, pty_nuf, dbeta, cb, ntrgrd);
2176
2177 /* Update the RHS and the local system matrix */
2178 const double pcoef = chi/sqrt(cm->face[f].meas);
2179
2180 for (short int v = 0; v < fm->n_vf; v++) {
2181
2182 /* Set the penalty diagonal coefficient */
2183 const short int vi = fm->v_ids[v];
2184 const double pcoef_v = pcoef * fm->wvf[v];
2185
2186 ntrgrd->val[vi*(1 + ntrgrd->n_rows)] += pcoef_v;
2187
2188 /* Update the RHS */
2189 for (int k = 0; k < 3; k++)
2190 csys->rhs[3*vi+k] += pcoef_v * csys->dir_values[3*vi+k];
2191
2192 } /* Dirichlet or homogeneous Dirichlet */
2193
2194 } /* Dirichlet face */
2195 } /* Loop on boundary faces */
2196
2197 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2198 if (cs_dbg_cw_test(eqp, cm, csys)) {
2199 cs_log_printf(CS_LOG_DEFAULT, ">> Cell Vb.COST Weak bc matrix");
2200 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2201 }
2202 #endif
2203
2204 /* Add contribution to the linear system */
2205 assert(pty->is_iso == true); /* if not the case something else TODO ? */
2206
2207 for (int bvi = 0; bvi < cm->n_vc; bvi++) {
2208 for (int bvj = 0; bvj < cm->n_vc; bvj++) {
2209
2210 /* Retrieve the 3x3 matrix */
2211 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bvi, bvj);
2212 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
2213
2214 const cs_real_t _val = ntrgrd->val[cm->n_vc*bvi + bvj];
2215 /* Update diagonal terms only */
2216 bij->val[0] += _val;
2217 bij->val[4] += _val;
2218 bij->val[8] += _val;
2219
2220 }
2221 }
2222
2223 }
2224
2225 /*----------------------------------------------------------------------------*/
2226 /*!
2227 * \brief Take into account a sliding BCs.
2228 * Case of vector-valued CDO-Vb schemes with a OCS algorithm.
2229 * Orthogonal splitting between Consistency/Stabilization parts.
2230 * Predefined prototype to match the function pointer
2231 * cs_cdo_enforce_bc_t
2232 *
2233 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2234 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2235 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2236 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2237 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2238 * \param[in, out] csys structure storing the cellwise system
2239 */
2240 /*----------------------------------------------------------------------------*/
2241
2242 void
cs_cdo_diffusion_vvb_ocs_sliding(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2243 cs_cdo_diffusion_vvb_ocs_sliding(const cs_equation_param_t *eqp,
2244 const cs_cell_mesh_t *cm,
2245 cs_face_mesh_t *fm,
2246 cs_hodge_t *hodge,
2247 cs_cell_builder_t *cb,
2248 cs_cell_sys_t *csys)
2249 {
2250 /* Enforcement of the sliding BCs */
2251 if (csys->has_sliding == false)
2252 return; /* Nothing to do */
2253
2254 const cs_hodge_param_t *hodgep = hodge->param;
2255 const cs_real_t dbeta =
2256 (hodgep->algo == CS_HODGE_ALGO_BUBBLE) ? hodgep->coef : 3*hodgep->coef;
2257 const cs_property_data_t *pty = hodge->pty_data;
2258
2259 /* Sanity checks */
2260 assert(pty->is_iso == true); /* if not the case something else TODO ? */
2261
2262 cs_sdm_t *ntrgrd = cb->loc;
2263
2264 /* Initialize the local operator (same as the scalar-valued one) */
2265 cs_sdm_square_init(cm->n_vc, ntrgrd);
2266
2267 for (short int i = 0; i < csys->n_bc_faces; i++) {
2268
2269 /* Get the boundary face in the cell numbering */
2270 const short int f = csys->_f_ids[i];
2271
2272 if (csys->bf_flag[f] & CS_CDO_BC_SLIDING) {
2273
2274 /* Compute the face-view of the mesh */
2275 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2276
2277 const cs_real_t nf[3] = {fm->face.unitv[0],
2278 fm->face.unitv[1],
2279 fm->face.unitv[2]};
2280
2281 /* Compute the product: matpty*face unit normal */
2282 cs_real_t pty_nuf[3] = {pty->value * nf[0],
2283 pty->value * nf[1],
2284 pty->value * nf[2]};
2285
2286 /* Compute the flux operator related to the trace on the current face
2287 * of the normal gradient
2288 * cb->values (size = cm->n_vc) is used inside
2289 */
2290 _vb_ocs_normal_flux_op(f, cm, pty_nuf, dbeta, cb, ntrgrd);
2291
2292 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2293 if (cs_dbg_cw_test(eqp, cm, csys)) {
2294 cs_log_printf(CS_LOG_DEFAULT, ">> Cell Vb.COST Sliding bc matrix");
2295 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2296 }
2297 #endif
2298
2299 /* Add contribution to the linear system of the penalized part */
2300 const double pcoef = eqp->weak_pena_bc_coeff/sqrt(cm->face[f].meas);
2301
2302 for (int vfi = 0; vfi < fm->n_vf; vfi++) {
2303
2304 const short int bvi = fm->v_ids[vfi];
2305
2306 /* Retrieve the diagonal 3x3 block matrix */
2307 cs_sdm_t *bii = cs_sdm_get_block(csys->mat, bvi, bvi);
2308
2309 /* Penalized part (u.n)(v.n) and tangential flux part */
2310 const cs_real_t vcoef = pcoef * fm->wvf[vfi];
2311 const cs_real_t nii = ntrgrd->val[cm->n_vc*bvi + bvi];
2312 const cs_real_t bval = vcoef + 2*nii;
2313
2314 for (int k = 0; k < 3; k++) {
2315 bii->val[3*k ] += bval * nf[0]*nf[k];
2316 bii->val[3*k+1] += bval * nf[1]*nf[k];
2317 bii->val[3*k+2] += bval * nf[2]*nf[k];
2318 }
2319
2320 } /* Loop on face vertices */
2321
2322 /* Update the system matrix on the extra-diagonal blocks */
2323 for (int bvi = 0; bvi < cm->n_vc; bvi++) {
2324 for (int bvj = 0; bvj < cm->n_vc; bvj++) {
2325
2326 if (bvi == bvj)
2327 continue;
2328
2329 /* Retrieve the 3x3 matrix */
2330 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bvi, bvj);
2331 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
2332
2333 const cs_real_t nij = ntrgrd->val[cm->n_vc*bvi + bvj];
2334 const cs_real_t nji = ntrgrd->val[cm->n_vc*bvj + bvi];
2335
2336 for (int k = 0; k < 3; k++) {
2337 bij->val[3*k ] += nf[0]*nf[k] * (nij + nji);
2338 bij->val[3*k+1] += nf[1]*nf[k] * (nij + nji);
2339 bij->val[3*k+2] += nf[2]*nf[k] * (nij + nji);
2340 }
2341
2342 } /* vfj */
2343 } /* vfi */
2344
2345 } /* Sliding face */
2346 } /* Loop on boundary faces */
2347
2348 }
2349
2350 /*----------------------------------------------------------------------------*/
2351 /*!
2352 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2353 * technique plus a symmetric treatment. Case of CDO-Vb schemes with a
2354 * COST/Bubble or Voronoi algorithm. One assumes an Orthogonal
2355 * splitting between Consistency/Stabilization parts (OCS).
2356 * Predefined prototype to match the function pointer
2357 * cs_cdo_enforce_bc_t
2358 *
2359 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2360 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2361 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2362 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2363 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2364 * \param[in, out] csys structure storing the cellwise system
2365 */
2366 /*----------------------------------------------------------------------------*/
2367
2368 void
cs_cdo_diffusion_svb_ocs_wsym_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2369 cs_cdo_diffusion_svb_ocs_wsym_dirichlet(const cs_equation_param_t *eqp,
2370 const cs_cell_mesh_t *cm,
2371 cs_face_mesh_t *fm,
2372 cs_hodge_t *hodge,
2373 cs_cell_builder_t *cb,
2374 cs_cell_sys_t *csys)
2375 {
2376 assert(csys != NULL); /* Sanity checks */
2377
2378 /* Enforcement of the Dirichlet BCs */
2379 if (csys->has_dirichlet == false)
2380 return; /* Nothing to do */
2381
2382 assert(cm != NULL && cb != NULL && hodge != NULL);
2383
2384 const cs_hodge_param_t *hodgep = hodge->param;
2385 const cs_property_data_t *pty = hodge->pty_data;
2386 const double chi =
2387 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2388
2389 assert(pty->need_tensor);
2390
2391 cs_sdm_t *ntrgrd = cb->loc;
2392
2393 for (short int i = 0; i < csys->n_bc_faces; i++) {
2394
2395 /* Get the boundary face in the cell numbering */
2396 const short int f = csys->_f_ids[i];
2397
2398 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2399
2400 /* Compute the face-view of the mesh */
2401 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2402
2403 /* Compute the product: matpty*face unit normal */
2404 cs_real_3_t pty_nuf;
2405 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2406
2407 /* Initialize the local operator */
2408 cs_sdm_square_init(cm->n_vc, ntrgrd);
2409
2410 /* Compute the flux operator related to the trace on the current face
2411 of the normal gradient */
2412 _vb_ocs_normal_flux_op(f, cm, pty_nuf, hodgep->coef, cb, ntrgrd);
2413
2414 /* Update ntrgrd = ntrgrd + transp and transp = transpose(ntrgrd) */
2415 cs_sdm_t *ntrgrd_tr = cb->aux;
2416 cs_sdm_square_add_transpose(ntrgrd, ntrgrd_tr);
2417
2418 /* Update RHS according to the add of ntrgrd_tr (cb->aux) */
2419 cs_sdm_square_matvec(ntrgrd_tr, csys->dir_values, cb->values);
2420 for (short int v = 0; v < csys->n_dofs; v++)
2421 csys->rhs[v] += cb->values[v];
2422
2423 /* Update the RHS and the local system matrix */
2424 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2425
2426 /* Add contribution to the linear system */
2427 cs_sdm_add(csys->mat, ntrgrd);
2428
2429 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2430 if (cs_dbg_cw_test(eqp, cm, csys)) {
2431 cs_log_printf(CS_LOG_DEFAULT,
2432 ">> Cell Vb COST WeakSym bc matrix (f_id: %d)", fm->f_id);
2433 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2434 }
2435 #endif
2436 } /* Dirichlet face */
2437 } /* Loop on boundary faces */
2438
2439 }
2440
2441 /*----------------------------------------------------------------------------*/
2442 /*!
2443 * \brief Take into account Robin BCs.
2444 * Case of scalar-valued CDO-Vb schemes with a WBS algorithm.
2445 * Predefined prototype to match the function pointer
2446 * cs_cdo_enforce_bc_t
2447 *
2448 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2449 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2450 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2451 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2452 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2453 * \param[in, out] csys structure storing the cellwise system
2454 */
2455 /*----------------------------------------------------------------------------*/
2456
2457 void
cs_cdo_diffusion_svb_wbs_robin(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2458 cs_cdo_diffusion_svb_wbs_robin(const cs_equation_param_t *eqp,
2459 const cs_cell_mesh_t *cm,
2460 cs_face_mesh_t *fm,
2461 cs_hodge_t *hodge,
2462 cs_cell_builder_t *cb,
2463 cs_cell_sys_t *csys)
2464 {
2465 CS_UNUSED(eqp);
2466 CS_UNUSED(hodge);
2467
2468 assert(csys != NULL); /* Sanity checks */
2469
2470 /* Enforcement of the Robin BCs */
2471 if (csys->has_robin == false)
2472 return; /* Nothing to do */
2473
2474 assert(cm != NULL && cb != NULL);
2475
2476 /* Robin BC expression: K du/dn + alpha*(u - u0) = g
2477 * Store x = alpha*u0 + g
2478 */
2479 cs_real_t *x = cb->values;
2480 cs_sdm_t *bc_op = cb->loc;
2481 cs_sdm_t *hloc = cb->aux; /* 2D Hodge operator on a boundary face */
2482
2483 for (short int i = 0; i < csys->n_bc_faces; i++) {
2484
2485 /* Get the boundary face in the cell numbering */
2486 const short int f = csys->_f_ids[i];
2487
2488 if (csys->bf_flag[f] & CS_CDO_BC_ROBIN) {
2489
2490 /* Reset local operator */
2491 cs_sdm_square_init(csys->n_dofs, bc_op);
2492
2493 /* Compute the face-view of the mesh */
2494 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2495
2496 cs_hodge_compute_wbs_surfacic(fm, hloc); /* hloc is of size n_vf */
2497
2498 /* Robin BC expression: K du/dn + alpha*(u - u0) = g */
2499 /* ------------------------------------------------- */
2500
2501 const double alpha = csys->rob_values[3*f];
2502 const double u0 = csys->rob_values[3*f+1];
2503 const double g = csys->rob_values[3*f+2];
2504
2505 memset(x, 0, sizeof(cs_real_t)*cm->n_vc);
2506 for (short int v = 0; v < fm->n_vf; v++) {
2507 const short int vi = fm->v_ids[v];
2508 x[vi] = alpha*u0 + g;
2509 }
2510
2511 /* Update the RHS and the local system */
2512 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
2513
2514 const short int vi = fm->v_ids[vfi];
2515 const cs_real_t *hfi = hloc->val + vfi*fm->n_vf;
2516 cs_real_t *opi = bc_op->val + vi*bc_op->n_rows;
2517
2518 for (short int vfj = 0; vfj < fm->n_vf; vfj++) {
2519
2520 const short int vj = fm->v_ids[vfj];
2521 csys->rhs[vi] += hfi[vfj]*x[vj];
2522 opi[vj] += alpha * hfi[vfj];
2523
2524 }
2525
2526 }
2527
2528 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2529 if (cs_dbg_cw_test(eqp, cm, csys)) {
2530 cs_log_printf(CS_LOG_DEFAULT,
2531 ">> Cell Vb WBS Robin bc matrix (f_id: %d)", fm->f_id);
2532 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, bc_op);
2533 }
2534 #endif
2535
2536 /* Add contribution to the linear system */
2537 cs_sdm_add(csys->mat, bc_op);
2538
2539 } /* Robin face */
2540 } /* Loop on boundary faces */
2541
2542 }
2543
2544 /*----------------------------------------------------------------------------*/
2545 /*!
2546 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2547 * technique. Case of CDO-Vb schemes with a WBS algorithm.
2548 * Predefined prototype to match the function pointer
2549 * cs_cdo_enforce_bc_t
2550 *
2551 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2552 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2553 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2554 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2555 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2556 * \param[in, out] csys structure storing the cellwise system
2557 */
2558 /*----------------------------------------------------------------------------*/
2559
2560 void
cs_cdo_diffusion_svb_wbs_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2561 cs_cdo_diffusion_svb_wbs_weak_dirichlet(const cs_equation_param_t *eqp,
2562 const cs_cell_mesh_t *cm,
2563 cs_face_mesh_t *fm,
2564 cs_hodge_t *hodge,
2565 cs_cell_builder_t *cb,
2566 cs_cell_sys_t *csys)
2567 {
2568 assert(csys != NULL); /* Sanity checks */
2569
2570 /* Enforcement of the Dirichlet BCs */
2571 if (csys->has_dirichlet == false)
2572 return; /* Nothing to do */
2573
2574 assert(cm != NULL && cb != NULL && hodge != NULL);
2575
2576 const cs_property_data_t *pty = hodge->pty_data;
2577 const double chi =
2578 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2579
2580 assert(pty->need_tensor);
2581
2582 cs_sdm_t *ntrgrd = cb->loc;
2583
2584 for (short int i = 0; i < csys->n_bc_faces; i++) {
2585
2586 /* Get the boundary face in the cell numbering */
2587 const short int f = csys->_f_ids[i];
2588
2589 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2590
2591 /* Compute the face-view of the mesh */
2592 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2593
2594 /* Compute the product: matpty*face unit normal */
2595 cs_real_3_t pty_nuf;
2596 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2597
2598 /* Compute the flux operator related to the trace on the current face
2599 of the normal gradient */
2600 _vb_wbs_normal_flux_op(fm, cm, pty_nuf, cb, ntrgrd);
2601
2602 /* Update the RHS and the local system matrix */
2603 #if 1 /* Default choice */
2604 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2605 #else /* This option seems less robust w.r.t the linear algebra */
2606 _wbs_nitsche(chi/sqrt(cm->face[f].meas), fm, ntrgrd, cb, csys);
2607 #endif
2608
2609 /* Add contribution to the linear system */
2610 cs_sdm_add(csys->mat, ntrgrd);
2611
2612 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2613 if (cs_dbg_cw_test(eqp, cm, csys)) {
2614 cs_log_printf(CS_LOG_DEFAULT,
2615 ">> Cell Vb WBS Weak bc matrix (f_id: %d)", fm->f_id);
2616 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2617 }
2618 #endif
2619 } /* Dirichlet face */
2620 } /* Loop on boundary faces */
2621
2622 }
2623
2624 /*----------------------------------------------------------------------------*/
2625 /*!
2626 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2627 * technique plus a symmetric treatment. Case of CDO-Vb schemes with a
2628 * WBS algorithm
2629 * Predefined prototype to match the function pointer
2630 * cs_cdo_enforce_bc_t
2631 *
2632 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2633 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2634 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2635 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2636 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2637 * \param[in, out] csys structure storing the cellwise system
2638 */
2639 /*----------------------------------------------------------------------------*/
2640
2641 void
cs_cdo_diffusion_svb_wbs_wsym_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2642 cs_cdo_diffusion_svb_wbs_wsym_dirichlet(const cs_equation_param_t *eqp,
2643 const cs_cell_mesh_t *cm,
2644 cs_face_mesh_t *fm,
2645 cs_hodge_t *hodge,
2646 cs_cell_builder_t *cb,
2647 cs_cell_sys_t *csys)
2648 {
2649 assert(csys != NULL); /* Sanity checks */
2650
2651 /* Enforcement of the Dirichlet BCs */
2652 if (csys->has_dirichlet == false)
2653 return; /* Nothing to do */
2654
2655 assert(cm != NULL && cb != NULL && hodge != NULL);
2656
2657 const cs_property_data_t *pty = hodge->pty_data;
2658 const double chi =
2659 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2660
2661 assert(pty->need_tensor);
2662
2663 cs_sdm_t *ntrgrd = cb->loc;
2664
2665 for (short int i = 0; i < csys->n_bc_faces; i++) {
2666
2667 /* Get the boundary face in the cell numbering */
2668 const short int f = csys->_f_ids[i];
2669
2670 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2671
2672 /* Compute the face-view of the mesh */
2673 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2674
2675 /* Compute the product: matpty*face unit normal */
2676 cs_real_3_t pty_nuf;
2677 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2678
2679 /* Initialize the local operator */
2680 cs_sdm_square_init(cm->n_vc, ntrgrd);
2681
2682 /* Compute the flux operator related to the trace on the current face
2683 of the normal gradient */
2684 _vb_wbs_normal_flux_op(fm, cm, pty_nuf, cb, ntrgrd);
2685
2686 /* Update ntrgrd = ntrgrd + transp and transp = transpose(ntrgrd) cb->loc
2687 plays the role of the flux operator */
2688 cs_sdm_t *ntrgrd_tr = cb->aux;
2689 cs_sdm_square_add_transpose(ntrgrd, ntrgrd_tr);
2690
2691 /* Update RHS according to the add of ntrgrd_tr (cb->aux) */
2692 cs_sdm_square_matvec(ntrgrd_tr, csys->dir_values, cb->values);
2693 for (short int v = 0; v < csys->n_dofs; v++)
2694 csys->rhs[v] += cb->values[v];
2695
2696 /* Update the RHS and the local system matrix */
2697 #if 1 /* Default choice */
2698 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2699 #else /* This option seems less robust w.r.t the linear algebra */
2700 _wbs_nitsche(chi/sqrt(cm->face[f].meas), fm, ntrgrd, cb, csys);
2701 #endif
2702
2703 /* Add contribution to the linear system */
2704 cs_sdm_add(csys->mat, ntrgrd);
2705
2706 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2707 if (cs_dbg_cw_test(eqp, cm, csys)) {
2708 cs_log_printf(CS_LOG_DEFAULT,
2709 ">> Cell Vb WBS WeakSym bc matrix (f_id: %d)", fm->f_id);
2710 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2711 }
2712 #endif
2713 } /* Dirichlet face */
2714 } /* Loop on boundary faces */
2715
2716 }
2717
2718 /*----------------------------------------------------------------------------*/
2719 /*!
2720 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2721 * technique. Case of CDO-VCb schemes with a WBS algorithm.
2722 * Predefined prototype to match the function pointer
2723 * cs_cdo_enforce_bc_t
2724 *
2725 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2726 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2727 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2728 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2729 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2730 * \param[in, out] csys structure storing the cellwise system
2731 */
2732 /*----------------------------------------------------------------------------*/
2733
2734 void
cs_cdo_diffusion_vcb_weak_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2735 cs_cdo_diffusion_vcb_weak_dirichlet(const cs_equation_param_t *eqp,
2736 const cs_cell_mesh_t *cm,
2737 cs_face_mesh_t *fm,
2738 cs_hodge_t *hodge,
2739 cs_cell_builder_t *cb,
2740 cs_cell_sys_t *csys)
2741 {
2742 assert(csys != NULL); /* Sanity checks */
2743
2744 /* Enforcement of the Dirichlet BCs */
2745 if (csys->has_dirichlet == false)
2746 return; /* Nothing to do */
2747
2748 assert(cm != NULL && cb != NULL && hodge != NULL);
2749
2750 const cs_property_data_t *pty = hodge->pty_data;
2751 const double chi =
2752 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2753
2754 assert(pty->need_tensor);
2755
2756 cs_sdm_t *ntrgrd = cb->loc;
2757
2758 for (short int i = 0; i < csys->n_bc_faces; i++) {
2759
2760 /* Get the boundary face in the cell numbering */
2761 const short int f = csys->_f_ids[i];
2762
2763 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2764
2765 /* Compute the face-view of the mesh */
2766 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2767
2768 /* Compute the product: matpty*face unit normal */
2769 cs_real_3_t pty_nuf;
2770 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2771
2772 /* Compute the flux operator related to the trace on the current face
2773 of the normal gradient */
2774 _vcb_wbs_normal_flux_op(fm, cm, pty_nuf, cb, ntrgrd);
2775
2776 /* Update the RHS and the local system matrix */
2777 #if 1 /* Default choice */
2778 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2779 #else /* This option seems less robust w.r.t the linear algebra */
2780 _wbs_nitsche(chi/sqrt(cm->face[f].meas), fm, ntrgrd, cb, csys);
2781 #endif
2782
2783 /* Add contribution to the linear system */
2784 cs_sdm_add(csys->mat, ntrgrd);
2785
2786 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2787 if (cs_dbg_cw_test(eqp, cm, csys)) {
2788 cs_log_printf(CS_LOG_DEFAULT,
2789 ">> Cell VCb Weak bc matrix (f_id: %d)", fm->f_id);
2790 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2791 }
2792 #endif
2793 } /* Dirichlet face */
2794 } /* Loop on boundary faces */
2795
2796 }
2797
2798 /*----------------------------------------------------------------------------*/
2799 /*!
2800 * \brief Take into account Dirichlet BCs by a weak enforcement using Nitsche
2801 * technique plus a symmetric treatment. Case of CDO-VCb schemes with
2802 * a WBS algorithm
2803 * Predefined prototype to match the function pointer
2804 * cs_cdo_enforce_bc_t
2805 *
2806 * \param[in] eqp pointer to a \ref cs_equation_param_t struct.
2807 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2808 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2809 * \param[in, out] hodge pointer to a \ref cs_hodge_t structure
2810 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2811 * \param[in, out] csys structure storing the cellwise system
2812 */
2813 /*----------------------------------------------------------------------------*/
2814
2815 void
cs_cdo_diffusion_vcb_wsym_dirichlet(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_face_mesh_t * fm,cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2816 cs_cdo_diffusion_vcb_wsym_dirichlet(const cs_equation_param_t *eqp,
2817 const cs_cell_mesh_t *cm,
2818 cs_face_mesh_t *fm,
2819 cs_hodge_t *hodge,
2820 cs_cell_builder_t *cb,
2821 cs_cell_sys_t *csys)
2822 {
2823 assert(csys != NULL); /* Sanity checks */
2824
2825 /* Enforcement of the Dirichlet BCs */
2826 if (csys->has_dirichlet == false)
2827 return; /* Nothing to do */
2828
2829 assert(cm != NULL && cb != NULL && hodge != NULL);
2830
2831 const cs_property_data_t *pty = hodge->pty_data;
2832 const double chi =
2833 eqp->weak_pena_bc_coeff * fabs(pty->eigen_ratio)*pty->eigen_max;
2834
2835 assert(pty->need_tensor);
2836
2837 cs_sdm_t *ntrgrd = cb->loc;
2838
2839 for (short int i = 0; i < csys->n_bc_faces; i++) {
2840
2841 /* Get the boundary face in the cell numbering */
2842 const short int f = csys->_f_ids[i];
2843
2844 if (cs_cdo_bc_is_dirichlet(csys->bf_flag[f])) {
2845
2846 /* Compute the face-view of the mesh */
2847 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2848
2849 /* Compute the product: matpty*face unit normal */
2850 cs_real_3_t pty_nuf;
2851 cs_math_33_3_product(pty->tensor, fm->face.unitv, pty_nuf);
2852
2853 /* Compute the flux operator related to the trace on the current face
2854 of the normal gradient */
2855 _vcb_wbs_normal_flux_op(fm, cm, pty_nuf, cb, ntrgrd);
2856
2857 /* Update ntrgrd = ntrgrd + transp and transp = transpose(ntrgrd) cb->loc
2858 plays the role of the flux operator */
2859 cs_sdm_t *ntrgrd_tr = cb->aux;
2860 cs_sdm_square_add_transpose(ntrgrd, ntrgrd_tr);
2861
2862 /* Update RHS according to the add of ntrgrd_tr (cb->aux) */
2863 cs_sdm_square_matvec(ntrgrd_tr, csys->dir_values, cb->values);
2864 for (short int v = 0; v < csys->n_dofs; v++)
2865 csys->rhs[v] += cb->values[v];
2866
2867 /* Update the RHS and the local system matrix */
2868 #if 1 /* Default choice */
2869 _svb_nitsche(chi/sqrt(fm->face.meas), fm, ntrgrd, csys);
2870 #else /* This option seems less robust w.r.t the linear algebra */
2871 _wbs_nitsche(chi/sqrt(cm->face[f].meas), fm, ntrgrd, cb, csys);
2872 #endif
2873
2874 /* Add contribution to the linear system */
2875 cs_sdm_add(csys->mat, ntrgrd);
2876
2877 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_DIFFUSION_DBG > 0
2878 if (cs_dbg_cw_test(eqp, cm, csys)) {
2879 cs_log_printf(CS_LOG_DEFAULT,
2880 ">> Cell VCb WeakSym bc matrix (f_id: %d)", fm->f_id);
2881 cs_sdm_dump(csys->c_id, csys->dof_ids, csys->dof_ids, ntrgrd);
2882 }
2883 #endif
2884 } /* Dirichlet face */
2885 } /* Loop on boundary faces */
2886
2887 }
2888
2889 /*----------------------------------------------------------------------------*/
2890 /*!
2891 * \brief Compute the diffusive flux across dual faces for a given cell.
2892 * Use the same consistent approximation as in the discrete Hodge op.
2893 * for this computation.
2894 * This function is dedicated to vertex-based schemes.
2895 * Flux = -Consistent(Hdg) * GRAD(pot)
2896 * Predefined prototype to match the function pointer
2897 * cs_cdo_diffusion_cw_flux_t
2898 *
2899 * \param[in] cm pointer to a cs_cell_mesh_t structure
2900 * \param[in] pot values of the potential fields at specific locations
2901 * \param[in] hodge pointer to a \ref cs_hodge_t structure
2902 * \param[in, out] cb auxiliary structure for computing the flux
2903 * \param[in, out] flx values of the flux across specific entities
2904 */
2905 /*----------------------------------------------------------------------------*/
2906
2907 void
cs_cdo_diffusion_svb_get_dfbyc_flux(const cs_cell_mesh_t * cm,const double * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,double * flx)2908 cs_cdo_diffusion_svb_get_dfbyc_flux(const cs_cell_mesh_t *cm,
2909 const double *pot,
2910 const cs_hodge_t *hodge,
2911 cs_cell_builder_t *cb,
2912 double *flx)
2913 {
2914 if (flx == NULL)
2915 return;
2916
2917 /* Sanity checks */
2918 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_EV));
2919
2920 /* Cellwise DoFs related to the discrete gradient (size: n_ec) */
2921 double *gec = cb->values;
2922 for (short int e = 0; e < cm->n_ec; e++) {
2923
2924 const short int *v = cm->e2v_ids + 2*e;
2925 /* sgn_v2 = -sgn_v1; flux = - HDG * GRAD(P) */
2926 gec[e] = cm->e2v_sgn[e]*(pot[v[1]] - pot[v[0]]);
2927
2928 } /* Loop on cell edges */
2929
2930 /* Store the local fluxes. flux = -Hdg * grd_c(pdi_c)
2931 * hodge->matrix has been computed just before the call to this function */
2932 cs_sdm_square_matvec(hodge->matrix, gec, flx);
2933 }
2934
2935 /*----------------------------------------------------------------------------*/
2936 /*!
2937 * \brief Compute the constant approximation of the diffusive flux inside a
2938 * (primal) cell. Use the same consistent approximation as in the
2939 * discrete Hodge op. for this computation. This function is dedicated
2940 * to vertex-based schemes.
2941 * Flux = -Hdg * GRAD(pot)
2942 * Predefined prototype to match the function pointer
2943 * cs_cdo_diffusion_cw_flux_t
2944 *
2945 * \param[in] cm pointer to a cs_cell_mesh_t structure
2946 * \param[in] pot values of the potential fields at specific locations
2947 * \param[in] hodge pointer to a \ref cs_hodge_t structure
2948 * \param[in, out] cb auxiliary structure for computing the flux
2949 * \param[in, out] flx values of the flux inside the cell
2950 */
2951 /*----------------------------------------------------------------------------*/
2952
2953 void
cs_cdo_diffusion_svb_get_cell_flux(const cs_cell_mesh_t * cm,const double * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,double * flx)2954 cs_cdo_diffusion_svb_get_cell_flux(const cs_cell_mesh_t *cm,
2955 const double *pot,
2956 const cs_hodge_t *hodge,
2957 cs_cell_builder_t *cb,
2958 double *flx)
2959 {
2960 CS_UNUSED(cb);
2961
2962 if (flx == NULL)
2963 return;
2964
2965 const cs_property_data_t *pty = hodge->pty_data;
2966
2967 /* Sanity checks */
2968 assert(pty->need_tensor);
2969 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_EV | CS_FLAG_COMP_DFQ));
2970
2971 cs_real_t grd[3] = {0., 0., 0.};
2972
2973 /* Cellwise DoFs related to the discrete gradient (size: n_ec) */
2974 for (short int e = 0; e < cm->n_ec; e++) {
2975
2976 const short int *v = cm->e2v_ids + 2*e;
2977
2978 /* sgn_v1 = -sgn_v0; flux = - Kc * GRAD(P) */
2979 const double ge = cm->e2v_sgn[e]*(pot[v[1]] - pot[v[0]]);
2980 const double contrib = ge * cm->dface[e].meas;
2981 for (int k = 0; k < 3; k++)
2982 grd[k] += contrib * cm->dface[e].unitv[k];
2983
2984 } /* Loop on cell edges */
2985
2986 cs_math_33_3_product(pty->tensor, grd, flx);
2987
2988 const double invvol = 1/cm->vol_c;
2989 for (int k = 0; k < 3; k++) flx[k] *= invvol;
2990 }
2991
2992 /*----------------------------------------------------------------------------*/
2993 /*!
2994 * \brief Compute the normal flux for a face assuming only the knowledge of
2995 * the potential at cell vertices. Valid for algorithm relying on a
2996 * spliting of the consistency/stabilization part as in OCS (also
2997 * called CO+ST) or Bubble algorithm. This is used for reconstructing
2998 * the normal flux from the degrees of freedom. The contribution for
2999 * each vertex of the face is then computed.
3000 *
3001 * \param[in] f face id in the cell mesh
3002 * \param[in] cm pointer to a cs_cell_mesh_t structure
3003 * \param[in] pot array of values of the potential (all the mesh)
3004 * \param[in] hodge pointer to a \ref cs_hodge_t structure
3005 * \param[in, out] cb auxiliary structure for building the flux
3006 * \param[in, out] flux array of values to set (size: n_vc)
3007 */
3008 /*----------------------------------------------------------------------------*/
3009
3010 void
cs_cdo_diffusion_svb_vbyf_flux(short int f,const cs_cell_mesh_t * cm,const cs_real_t * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_real_t * flux)3011 cs_cdo_diffusion_svb_vbyf_flux(short int f,
3012 const cs_cell_mesh_t *cm,
3013 const cs_real_t *pot,
3014 const cs_hodge_t *hodge,
3015 cs_cell_builder_t *cb,
3016 cs_real_t *flux)
3017 {
3018 if (flux == NULL)
3019 return;
3020
3021 const cs_property_data_t *pty = hodge->pty_data;
3022 const cs_hodge_param_t *hodgep = hodge->param;
3023
3024 assert(pty->need_tensor);
3025 assert(cs_eflag_test(cm->flag,
3026 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_EV |
3027 CS_FLAG_COMP_DFQ | CS_FLAG_COMP_FE));
3028
3029 const cs_real_t beta =
3030 (hodgep->algo == CS_HODGE_ALGO_BUBBLE) ? hodgep->coef : 3*hodgep->coef;
3031 const cs_quant_t pfq = cm->face[f];
3032
3033 /* Reset the fluxes */
3034 memset(flux, 0, cm->n_vc*sizeof(cs_real_t));
3035
3036 /* Compute the product: matpty*face unit normal */
3037 cs_real_t pty_nuf[3] = {0, 0, 0};;
3038 cs_math_33_3_product(pty->tensor, pfq.unitv, pty_nuf);
3039
3040 /* Cellwise constant and consistent gradient */
3041 cs_real_t grd_cc[3] = {0, 0, 0};
3042 cs_real_t *g = cb->values;
3043
3044 /* Cellwise DoFs related to the discrete gradient (size: n_ec) */
3045 for (short int e = 0; e < cm->n_ec; e++) {
3046
3047 const short int *v = cm->e2v_ids + 2*e;
3048
3049 /* sgn_v1 = -sgn_v0; GRAD(P) */
3050 g[e] = cm->e2v_sgn[e]*(pot[v[0]] - pot[v[1]]);
3051
3052 const double ge_coef = g[e] * cm->dface[e].meas;
3053 for (int k = 0; k < 3; k++)
3054 grd_cc[k] += ge_coef * cm->dface[e].unitv[k];
3055
3056 } /* Loop on cell edges */
3057
3058 const double invvol = 1/cm->vol_c;
3059 for (int k = 0; k < 3; k++) grd_cc[k] *= invvol;
3060
3061 /* Add the stabilisation part which is constant on p_{e,c} --> t_ef if one
3062 restricts to the face f */
3063 for (int ie = cm->f2e_idx[f]; ie < cm->f2e_idx[f+1]; ie++) {
3064
3065 cs_real_t grd_tef[3] = {0, 0, 0};
3066
3067 const short int e = cm->f2e_ids[ie];
3068 const short int *v = cm->e2v_ids + 2*e;
3069 const cs_quant_t peq = cm->edge[e];
3070 const cs_nvec3_t dfq = cm->dface[e];
3071 const cs_real_t pec_coef = beta/(peq.meas*_dp3(peq.unitv, dfq.unitv));
3072 const cs_real_t delta = g[e] - peq.meas*_dp3(peq.unitv, grd_cc);
3073 const cs_real_t stab_coef = pec_coef * delta;
3074
3075 for (int k = 0; k < 3; k++)
3076 grd_tef[k] = grd_cc[k] + stab_coef * dfq.unitv[k];
3077
3078 const cs_real_t tef = (cs_eflag_test(cm->flag, CS_FLAG_COMP_FEQ)) ?
3079 cm->tef[ie] : cs_compute_area_from_quant(peq, pfq.center);
3080
3081 const double _flx = -0.5 * tef * _dp3(grd_tef, pty_nuf);
3082
3083 flux[v[0]] += _flx;
3084 flux[v[1]] += _flx;
3085
3086 } /* Loop on face edges */
3087
3088 }
3089
3090 /*----------------------------------------------------------------------------*/
3091 /*!
3092 * \brief Compute the diffusive flux across dual faces for a given cell
3093 * Use the WBS algo. for approximating the gradient
3094 * The computation takes into account a subdivision into tetrahedra of
3095 * the current cell based on p_{ef,c}
3096 * Predefined prototype to match the function pointer
3097 * cs_cdo_diffusion_cw_flux_t
3098 *
3099 * \param[in] cm pointer to a cs_cell_mesh_t structure
3100 * \param[in] pot values of the potential fields at vertices
3101 * \param[in] hodge pointer to a \ref cs_hodge_t structure
3102 * \param[in, out] cb auxiliary structure for computing the flux
3103 * \param[in, out] flx flux across dual faces inside this cell
3104 */
3105 /*----------------------------------------------------------------------------*/
3106
3107 void
cs_cdo_diffusion_wbs_get_dfbyc_flux(const cs_cell_mesh_t * cm,const cs_real_t * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_real_t * flx)3108 cs_cdo_diffusion_wbs_get_dfbyc_flux(const cs_cell_mesh_t *cm,
3109 const cs_real_t *pot,
3110 const cs_hodge_t *hodge,
3111 cs_cell_builder_t *cb,
3112 cs_real_t *flx)
3113 {
3114 /* Temporary buffers */
3115 cs_real_3_t *u_vc = cb->vectors;
3116 double *l_vc = cb->values;
3117
3118 const double *p_v = pot;
3119 const double p_c = pot[cm->n_vc];
3120 const cs_property_data_t *pty = hodge->pty_data;
3121
3122 /* Sanity checks */
3123 assert(cs_eflag_test(cm->flag,
3124 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
3125 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_SEF));
3126 assert(pty->need_tensor);
3127
3128 /* Reset local fluxes */
3129 for (short int e = 0; e < cm->n_ec; e++) flx[e] = 0.;
3130
3131 /* Store segments xv --> xc for this cell */
3132 for (short int v = 0; v < cm->n_vc; v++)
3133 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, l_vc + v, u_vc[v]);
3134
3135 /* Loop on cell faces */
3136 for (short int f = 0; f < cm->n_fc; f++) {
3137
3138 cs_real_3_t grd_c, grd_v1, grd_v2, grd_pef, mgrd;
3139
3140 /* Compute for the current face:
3141 - the area of each triangle defined by a base e and an apex f
3142 - the gradient of the Lagrange function related xc in p_{f,c} */
3143 cs_compute_grdfc_cw(f, cm, grd_c);
3144
3145 /* Compute the reconstructed value of the potential at p_f */
3146 double p_f = 0.;
3147 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
3148
3149 const short int e = cm->f2e_ids[i];
3150
3151 p_f += cm->tef[i]*( p_v[cm->e2v_ids[2*e]] /* p_v1 */
3152 + p_v[cm->e2v_ids[2*e+1]] ); /* p_v2 */
3153 }
3154 p_f *= 0.5/cm->face[f].meas;
3155
3156 const double dp_cf = p_c - p_f;
3157
3158 /* Loop on face edges to scan p_{ef,c} subvolumes */
3159 const cs_nvec3_t deq = cm->dedge[f];
3160 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
3161
3162 const short int e = cm->f2e_ids[i];
3163 const short int ee = 2*e;
3164 const short int v1 = cm->e2v_ids[ee];
3165 const short int v2 = cm->e2v_ids[ee+1];
3166
3167 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])u_vc, l_vc,
3168 grd_v1, grd_v2);
3169
3170 /* Gradient of the lagrange function related to a face.
3171 grd_f = -(grd_c + grd_v1 + grd_v2)
3172 This formula is a consequence of the Partition of the Unity.
3173 This yields the following formula for grd(Lv^conf)|_p_{ef,c} */
3174 const double dp1f = p_v[v1] - p_f, dp2f = p_v[v2] - p_f;
3175 for (int k = 0; k < 3; k++)
3176 grd_pef[k] = dp_cf*grd_c[k] + dp1f*grd_v1[k] + dp2f*grd_v2[k];
3177
3178 cs_math_33_3_product(pty->tensor, grd_pef, mgrd);
3179
3180 flx[e] -= cm->sefc[i].meas * _dp3(cm->sefc[i].unitv, mgrd);
3181
3182 } /* Loop on face edges */
3183
3184 } /* Loop on cell faces */
3185
3186 }
3187
3188 /*----------------------------------------------------------------------------*/
3189 /*!
3190 * \brief Compute the diffusive flux inside a given primal cell
3191 * Use the WBS algo. for approximating the gradient
3192 * The computation takes into account a subdivision into tetrahedra of
3193 * the current cell based on p_{ef,c}
3194 * Predefined prototype to match the function pointer
3195 * cs_cdo_diffusion_cw_flux_t
3196 *
3197 * \param[in] cm pointer to a cs_cell_mesh_t structure
3198 * \param[in] pot values of the potential fields at vertices
3199 * \param[in] hodge pointer to a \ref cs_hodge_t structure
3200 * \param[in, out] cb auxiliary structure for computing the flux
3201 * \param[in, out] flx flux vector inside this cell
3202 */
3203 /*----------------------------------------------------------------------------*/
3204
3205 void
cs_cdo_diffusion_wbs_get_cell_flux(const cs_cell_mesh_t * cm,const cs_real_t * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_real_t * flx)3206 cs_cdo_diffusion_wbs_get_cell_flux(const cs_cell_mesh_t *cm,
3207 const cs_real_t *pot,
3208 const cs_hodge_t *hodge,
3209 cs_cell_builder_t *cb,
3210 cs_real_t *flx)
3211 {
3212 const cs_property_data_t *pty = hodge->pty_data;
3213
3214 /* Sanity checks */
3215 assert(pty->need_tensor);
3216 assert(cs_eflag_test(cm->flag,
3217 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
3218 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ));
3219
3220 cs_real_t cgrd[3] = {0, 0, 0};
3221
3222 /* Compute the mean-value of the cell gradient */
3223 cs_reco_cw_cgrd_wbs_from_pvc(cm, pot, cb, cgrd);
3224
3225 cs_math_33_3_product(pty->tensor, cgrd, flx);
3226 for (int k = 0; k < 3; k++) flx[k] *= -1; /* Flux = - tensor * grd */
3227 }
3228
3229 /*----------------------------------------------------------------------------*/
3230 /*!
3231 * \brief Compute the normal flux for a face assuming only the knowledge
3232 * of the potential at cell vertices (and at cell center).
3233 * WBS algorithm is used for reconstructing the normal flux from the
3234 * degrees of freedom.
3235 *
3236 * \param[in] f face id in the cell mesh
3237 * \param[in] cm pointer to a cs_cell_mesh_t structure
3238 * \param[in] pot array of values of the potential (all the mesh)
3239 * \param[in] hodge pointer to a \ref cs_hodge_t structure
3240 * \param[in, out] cb auxiliary structure dedicated to diffusion
3241 * \param[in, out] vf_flux array of values to set (size: n_vc)
3242 */
3243 /*----------------------------------------------------------------------------*/
3244
3245 void
cs_cdo_diffusion_wbs_vbyf_flux(short int f,const cs_cell_mesh_t * cm,const cs_real_t * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_real_t * flux)3246 cs_cdo_diffusion_wbs_vbyf_flux(short int f,
3247 const cs_cell_mesh_t *cm,
3248 const cs_real_t *pot,
3249 const cs_hodge_t *hodge,
3250 cs_cell_builder_t *cb,
3251 cs_real_t *flux)
3252 {
3253 if (flux == NULL)
3254 return;
3255
3256 const cs_property_data_t *pty = hodge->pty_data;
3257
3258 /* Sanity checks */
3259 assert(pty->need_tensor);
3260 assert(hodge->param->algo == CS_HODGE_ALGO_WBS);
3261 assert(cs_eflag_test(cm->flag,
3262 CS_FLAG_COMP_FV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ |
3263 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_FE));
3264
3265 /* Reset the fluxes */
3266 memset(flux, 0, cm->n_vc*sizeof(cs_real_t));
3267
3268 /* Compute the product: matpty*face unit normal */
3269 cs_real_t mnuf[3] = {0, 0, 0};
3270 cs_math_33_3_product(pty->tensor, cm->face[f].unitv, mnuf);
3271
3272 /* Compute xc --> xv length and unit vector for all face vertices */
3273 double *l_vc = cb->values;
3274 cs_real_3_t *u_vc = cb->vectors;
3275 for (int i = cm->f2v_idx[f]; i < cm->f2v_idx[f+1]; i++) {
3276 short int v = cm->f2v_ids[i];
3277 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, l_vc + v, u_vc[v]);
3278 }
3279
3280 /* Compute for the current face, the gradient of the Lagrange function
3281 related to xc in p_{f,c} */
3282 cs_real_t grd_c[3], grd_pef[3], grd_v0[3], grd_v1[3];
3283 cs_compute_grdfc_cw(f, cm, grd_c);
3284
3285 const cs_real_t *p_v = pot;
3286 const cs_real_t p_f = cs_reco_cw_scalar_pv_at_face_center(f, cm, p_v);
3287
3288 /* Compute p_c - p_f (where p_c is the reconstructed values at the
3289 cell center */
3290 const double dp_cf = pot[cm->n_vc] - p_f;
3291
3292 /* Loop on face edges to scan p_{ef,c} subvolumes */
3293 for (int ie = cm->f2e_idx[f]; ie < cm->f2e_idx[f+1]; ie++) {
3294
3295 const short int *v = cm->e2v_ids + 2*cm->f2e_ids[ie];
3296
3297 /* Compute the gradient of the Lagrange function related xv0, xv1
3298 in each p_{ef,c} for e in E_f */
3299 cs_compute_grd_ve(v[0], v[1], cm->dedge[f],
3300 (const cs_real_t (*)[3])u_vc, l_vc,
3301 grd_v0, grd_v1);
3302
3303 /* Gradient of the lagrange function related to a face.
3304 grd_f = -(grd_c + grd_v1 + grd_v2)
3305 This formula is a consequence of the Partition of the Unity.
3306 This yields the following formula for grd(Lv^conf)|_p_{ef,c} */
3307 const double dp_v0f = p_v[v[0]] - p_f;
3308 const double dp_v1f = p_v[v[1]] - p_f;
3309 for (int k = 0; k < 3; k++)
3310 grd_pef[k] = dp_cf*grd_c[k] + dp_v0f*grd_v0[k] + dp_v1f*grd_v1[k];
3311
3312 /* tef: area of the triangle defined by the base e and the apex f */
3313 const double _flx = -0.5*cm->tef[ie] * _dp3(mnuf, grd_pef);
3314
3315 flux[v[0]] += _flx;
3316 flux[v[1]] += _flx;
3317
3318 } /* Loop on face edges */
3319
3320 }
3321
3322 /*----------------------------------------------------------------------------*/
3323 /*!
3324 * \brief Compute the diffusive flux across a face (based on a subdivision
3325 * into tetrahedra of the volume p_{f,c})
3326 *
3327 * \param[in] fm pointer to a cs_face_mesh_t structure
3328 * \param[in] pty_tens 3x3 matrix related to the diffusion property
3329 * \param[in] p_v array of values attached to face vertices
3330 * \param[in] p_f value attached to the face
3331 * \param[in] p_c value attached to the cell
3332 * \param[in, out] cb auxiliary structure dedicated to diffusion
3333 *
3334 * \return the value of the diffusive flux across the current face
3335 */
3336 /*----------------------------------------------------------------------------*/
3337
3338 double
cs_cdo_diffusion_wbs_face_flux(const cs_face_mesh_t * fm,const cs_real_t pty_tens[3][3],const double * p_v,const double p_f,const double p_c,cs_cell_builder_t * cb)3339 cs_cdo_diffusion_wbs_face_flux(const cs_face_mesh_t *fm,
3340 const cs_real_t pty_tens[3][3],
3341 const double *p_v,
3342 const double p_f,
3343 const double p_c,
3344 cs_cell_builder_t *cb)
3345 {
3346 cs_real_t grd_c[3], grd_pef[3], grd_v1[3], grd_v2[3], mnuf[3] = {0, 0, 0};
3347 double f_flux = 0.;
3348
3349 /* Retrieve temporary buffers */
3350 double *l_vc = cb->values;
3351 cs_real_3_t *u_vc = cb->vectors;
3352
3353 cs_math_33_3_product((const cs_real_t (*)[3])pty_tens, fm->face.unitv,
3354 mnuf);
3355
3356 /* Compute xc --> xv length and unit vector for all face vertices */
3357 for (short int v = 0; v < fm->n_vf; v++)
3358 cs_math_3_length_unitv(fm->xc, fm->xv + 3*v, l_vc + v, u_vc[v]);
3359
3360 /* Compute for the current face, the gradient of the Lagrange function
3361 related to xc in p_{f,c} */
3362 cs_compute_grdfc_fw(fm, grd_c);
3363
3364 /* Compute p_c - p_f (where p_c is the reconstructed values at the
3365 cell center */
3366 const double dp_cf = p_c - p_f;
3367
3368 /* Loop on face edges to scan p_{ef,c} subvolumes */
3369 for (int e = 0; e < fm->n_ef; e++) {
3370
3371 const short int v1 = fm->e2v_ids[2*e];
3372 const short int v2 = fm->e2v_ids[2*e+1];
3373
3374 /* Compute the gradient of the Lagrange function related xv1, xv2
3375 in each p_{ef,c} for e in E_f */
3376 cs_compute_grd_ve(v1, v2, fm->dedge, (const cs_real_t (*)[3])u_vc, l_vc,
3377 grd_v1, grd_v2);
3378
3379 /* Gradient of the lagrange function related to a face.
3380 grd_f = -(grd_c + grd_v1 + grd_v2)
3381 This formula is a consequence of the Partition of the Unity.
3382 This yields the following formula for grd(Lv^conf)|_p_{ef,c} */
3383 const double dp1f = p_v[v1] - p_f, dp2f = p_v[v2] - p_f;
3384 for (int k = 0; k < 3; k++)
3385 grd_pef[k] = dp_cf*grd_c[k] + dp1f*grd_v1[k] + dp2f*grd_v2[k];
3386
3387 /* Area of the triangle defined by the base e and the apex f */
3388 f_flux -= fm->tef[e] * _dp3(mnuf, grd_pef);
3389
3390 } /* Loop on face edges */
3391
3392 return f_flux;
3393 }
3394
3395 /*----------------------------------------------------------------------------*/
3396 /*!
3397 * \brief Compute the normal diffusive flux for a face assuming only the
3398 * knowledge of the potential at faces and cell.
3399 * CO+ST algorithm is used for reconstructing the normal flux from
3400 * the degrees of freedom.
3401 *
3402 * \param[in] f face id in the cell mesh
3403 * \param[in] cm pointer to a cs_cell_mesh_t structure
3404 * \param[in] pot array of values of the potential (all the mesh)
3405 * \param[in] hodge pointer to a \ref cs_hodge_t structure
3406 * \param[in, out] cb auxiliary structure dedicated to diffusion
3407 * \param[out] flux pointer to the value to set
3408 */
3409 /*----------------------------------------------------------------------------*/
3410
3411 void
cs_cdo_diffusion_sfb_cost_flux(short int f,const cs_cell_mesh_t * cm,const cs_real_t * pot,const cs_hodge_t * hodge,cs_cell_builder_t * cb,cs_real_t * flux)3412 cs_cdo_diffusion_sfb_cost_flux(short int f,
3413 const cs_cell_mesh_t *cm,
3414 const cs_real_t *pot,
3415 const cs_hodge_t *hodge,
3416 cs_cell_builder_t *cb,
3417 cs_real_t *flux)
3418 {
3419 if (flux == NULL)
3420 return;
3421
3422 const cs_property_data_t *pty = hodge->pty_data;
3423
3424 /* Sanity checks */
3425 assert(pty->need_tensor);
3426 assert(hodge->param->algo == CS_HODGE_ALGO_COST);
3427 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3428
3429 const cs_real_t beta = hodge->param->coef;
3430 const cs_quant_t pfq = cm->face[f];
3431
3432 /* Compute the product: matpty*face unit normal */
3433 cs_real_t pty_nuf[3] = {0, 0, 0};
3434 cs_math_33_3_product(pty->tensor, pfq.unitv, pty_nuf);
3435
3436 /* Define a cellwise constant and consistent gradient */
3437 /* -------------------------------------------------- */
3438
3439 cs_real_t grd_cc[3] = {0, 0, 0};
3440 cs_real_t *g = cb->values;
3441
3442 /* Cellwise DoFs related to the discrete gradient (size: n_fc) */
3443 for (short int ff = 0; ff < cm->n_fc; ff++) {
3444
3445 /* Delta of the potential along the dual edge (x_f - x_c) */
3446 g[ff] = cm->f_sgn[ff]*(pot[ff] - pot[cm->n_fc]);
3447
3448 const double gcoef = g[ff] * cm->face[ff].meas;
3449 for (int k = 0; k < 3; k++)
3450 grd_cc[k] += gcoef * cm->face[ff].unitv[k];
3451
3452 } /* Loop on cell faces */
3453
3454 const double invvol = 1/cm->vol_c;
3455 for (int k = 0; k < 3; k++) grd_cc[k] *= invvol;
3456
3457 /* Add the stabilisation part which is constant on p_{f,c} */
3458 const cs_nvec3_t *deq = cm->dedge + f;
3459 const cs_real_t pfc_coef = 3*beta/(_dp3(pfq.unitv, deq->unitv));
3460 const cs_real_t delta = g[f] - deq->meas*_dp3(deq->unitv, grd_cc);
3461 const cs_real_t stab_coef = pfc_coef * delta;
3462
3463 cs_real_t grd_reco[3] = {0., 0., 0.};
3464 for (int k = 0; k < 3; k++)
3465 grd_reco[k] = grd_cc[k] + stab_coef * pfq.unitv[k];
3466
3467 *flux = -pfq.meas * _dp3(grd_reco, pty_nuf);
3468 }
3469
3470 /*----------------------------------------------------------------------------*/
3471 /*!
3472 * \brief Compute the normal flux for a face assuming only the knowledge
3473 * of the potential at cell vertices. COST algorithm is used for
3474 * reconstructing a piecewise constant gradient from the degrees of
3475 * freedom.
3476 *
3477 * \param[in] f face id in the cell mesh
3478 * \param[in] cm pointer to a cs_cell_mesh_t structure
3479 * \param[in] diff_tensor property tensor times the face normal
3480 * \param[in] pot_values array of values of the potential (all the mesh)
3481 * \param[in, out] fluxes values of the fluxes related to each vertex
3482 *
3483 */
3484 /*----------------------------------------------------------------------------*/
3485
3486 void
cs_cdovb_diffusion_p0_face_flux(const short int f,const cs_cell_mesh_t * cm,const cs_real_t diff_tensor[3][3],const cs_real_t * pot_values,cs_real_t * fluxes)3487 cs_cdovb_diffusion_p0_face_flux(const short int f,
3488 const cs_cell_mesh_t *cm,
3489 const cs_real_t diff_tensor[3][3],
3490 const cs_real_t *pot_values,
3491 cs_real_t *fluxes)
3492 {
3493 assert(cs_eflag_test(cm->flag,
3494 CS_FLAG_COMP_PV | CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ));
3495
3496 cs_real_3_t mnuf;
3497 cs_real_3_t gc = {0, 0, 0};
3498
3499 cs_math_33_3_product(diff_tensor, cm->face[f].unitv, mnuf);
3500
3501 cs_reco_dfbyc_in_cell(cm, pot_values, gc);
3502
3503 for (short int v = 0; v < cm->n_vc; v++) fluxes[v] = 0;
3504
3505 const cs_real_t flux_coef = 0.5 * _dp3(gc, mnuf);
3506 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
3507 const cs_real_t _flx = flux_coef * cm->tef[i];
3508 const int eshft = 2*cm->f2e_ids[i];
3509
3510 fluxes[cm->e2v_ids[eshft]] -= _flx;
3511 fluxes[cm->e2v_ids[eshft+1]] -= _flx;
3512 }
3513 }
3514
3515 /*----------------------------------------------------------------------------*/
3516
3517 #undef _dp3
3518
3519 END_C_DECLS
3520