1 /*============================================================================
2 * Manage discrete Hodge operators and closely related operators
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <math.h>
36 #include <float.h>
37 #include <limits.h>
38 #include <assert.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include "bft_mem.h"
45
46 #include "cs_log.h"
47 #include "cs_math.h"
48 #include "cs_scheme_geometry.h"
49
50 /*----------------------------------------------------------------------------
51 * Header for the current file
52 *----------------------------------------------------------------------------*/
53
54 #include "cs_hodge.h"
55
56 /*----------------------------------------------------------------------------*/
57
58 BEGIN_C_DECLS
59
60 /*=============================================================================
61 * Additional doxygen documentation
62 *============================================================================*/
63
64 /*!
65 \file cs_hodge.c
66
67 \brief Build discrete Hodge operators
68
69 */
70
71 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
72
73 /*=============================================================================
74 * Local Macro definitions
75 *============================================================================*/
76
77 #define CS_HODGE_DBG 0
78 #define CS_HODGE_MODULO 1
79
80 /* Redefined the name of functions from cs_math to get shorter names */
81 #define _dp3 cs_math_3_dot_product
82
83 /*=============================================================================
84 * Local structure definitions
85 *============================================================================*/
86
87 /*============================================================================
88 * Local variables
89 *============================================================================*/
90
91 /*============================================================================
92 * Private constant variables
93 *============================================================================*/
94
95 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
96
97 static const double cs_hodge_vc_coef = 3./20;
98
99 static const char
100 cs_hodge_type_desc[CS_HODGE_N_TYPES][CS_BASE_STRING_LEN] =
101 { N_("VpCd"),
102 N_("EpFd"),
103 N_("FpEd"),
104 N_("EdFp"),
105 N_("CpVd") };
106
107 static const char
108 cs_hodge_algo_desc[CS_HODGE_N_ALGOS][CS_BASE_STRING_LEN] =
109 { N_("Voronoi"),
110 N_("Whitney on the Barycentric Subdivision (WBS)"),
111 N_("Orthogonal Consistency/Stabilization (OCS)"),
112 N_("Orthogonal Consistency/Sub-Stabilization (OCS2)"),
113 N_("Orthogonal Consistency/Bubble-Stabilization (BUBBLE)"),
114 N_("Automatic switch") };
115
116 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
117
118 /*============================================================================
119 * Private function prototypes for debugging purpose
120 *============================================================================*/
121
122 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
123 /*----------------------------------------------------------------------------*/
124 /*!
125 * \brief Check the coherency of the values of a stiffness matrix
126 *
127 * \param[in] c_id current cell id
128 * \param[in] sloc pointer to a cs_sdm_t struct.
129 */
130 /*----------------------------------------------------------------------------*/
131
132 static void
_check_stiffness(cs_lnum_t c_id,const cs_sdm_t * sloc)133 _check_stiffness(cs_lnum_t c_id,
134 const cs_sdm_t *sloc)
135 {
136 if (c_id % CS_HODGE_MODULO != 0 || sloc == NULL)
137 return;
138
139 cs_log_printf(CS_LOG_DEFAULT, ">> Local stiffness matrix");
140 cs_sdm_dump(c_id, NULL, NULL, sloc);
141
142 double print_val = 0.;
143
144 for (int i = 0 ; i < sloc->n_rows; i++) {
145 double rsum = 0.;
146 const cs_real_t *rval = sloc->val + i*sloc->n_rows;
147 for (cs_lnum_t j = 0; j < sloc->n_rows; j++)
148 rsum += rval[j];
149
150 print_val += fabs(rsum);
151
152 if (rsum > 100*cs_math_get_machine_epsilon()) {
153 cs_base_warn(__FILE__, __LINE__);
154 cs_log_printf(CS_LOG_DEFAULT, " %s: row %d Row sum = %5.3e > 0 !\n",
155 __func__, i, rsum);
156 }
157
158 }
159 cs_log_printf(CS_LOG_DEFAULT, " %s: err = % -5.3e\n", __func__, print_val);
160 }
161
162 /*----------------------------------------------------------------------------*/
163 /*!
164 * \brief Check the coherency of the values of a discrete Hodge operator
165 *
166 * \param[in] c_id current cell id
167 * \param[in] vec vectors of quantities to test against a hodge
168 * \param[in] res vectors of quantities to compare with
169 * \param[in] hodge pointer to a discrete Hodge operator
170 * \param[in, out] cb pointer to a cell builder structure
171 * buffers to store temporary values
172 */
173 /*----------------------------------------------------------------------------*/
174
175 static void
_check_vector_hodge(cs_lnum_t c_id,const cs_real_3_t * vec,const cs_real_3_t * res,cs_hodge_t * hodge,cs_cell_builder_t * cb)176 _check_vector_hodge(cs_lnum_t c_id,
177 const cs_real_3_t *vec,
178 const cs_real_3_t *res,
179 cs_hodge_t *hodge,
180 cs_cell_builder_t *cb)
181
182 {
183 if (c_id % CS_HODGE_MODULO != 0 || hmat == NULL)
184 return;
185
186 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
187 cs_sdm_dump(c_id, NULL, NULL, hodge->matrix);
188
189 const cs_hodge_param_t *hodgep = hodge->param;
190 const cs_property_data_t *ptyd = hodge->pty_data;
191 const cs_sdm_t *hmat = hodge->matrix;
192
193 cs_real_t *in = cb->values;
194 cs_real_t *h_in = cb->values + hmat->n_rows;
195 cs_real_t *ref = cb->values + 2*hmat->n_rows;
196 double print_val = 0.;
197
198 if (ptyd->is_unity) {
199 ptyd->tensor[0][0] = ptyd->tensor[1][1] = ptyd->tensor[2][2] = 1;
200 ptyd->tensor[0][1] = ptyd->tensor[1][0] = ptyd->tensor[2][0] = 0;
201 ptyd->tensor[0][2] = ptyd->tensor[1][2] = ptyd->tensor[2][1] = 0;
202
203 }
204 else if (ptyd->is_iso) {
205 ptyd->tensor[0][0] = ptyd->tensor[1][1] = ptyd->tensor[2][2] = ptyd->value;
206 ptyd->tensor[0][1] = ptyd->tensor[1][0] = ptyd->tensor[2][0] = 0;
207 ptyd->tensor[0][2] = ptyd->tensor[1][2] = ptyd->tensor[2][1] = 0;
208 }
209
210 const cs_real_3_t a[3] = { {1, 0, 0}, {0, 1, 0}, {0, 0, 1} };
211
212 for (int dim = 0; dim < 3; dim++) {
213
214 cs_real_3_t pty_a;
215 cs_math_33_3_product((const cs_real_t (*)[3])ptyd->tensor, a[dim], pty_a);
216
217 for (int i = 0; i < hmat->n_rows; i++) {
218 in[i] = vec[i][dim];
219 ref[i] = _dp3(pty_a, res[i]);
220 }
221
222 cs_sdm_square_matvec(hmat, in, h_in);
223
224 double err = 0.;
225 for (int i = 0; i < hmat->n_rows; i++)
226 err += fabs(ref[i] - h_in[i]);
227 print_val += err;
228 if (err > 100 * cs_math_get_machine_epsilon()) {
229 cs_base_warn(__FILE__, __LINE__);
230 cs_log_printf(CS_LOG_DEFAULT,
231 " %s: err = %5.3e, dim = %d\n", __func__, err, dim);
232 }
233
234 }
235
236 cs_log_printf(CS_LOG_DEFAULT, "%s: err = % -5.3e\n", __func__, print_val);
237
238 }
239 #endif /* DEBUG */
240
241 /*============================================================================
242 * Private function prototypes
243 *============================================================================*/
244
245 /*----------------------------------------------------------------------------*/
246 /*!
247 * \brief Perform a part of a matrix-vector product
248 *
249 * \param[in] n_ent size of the matrix
250 * \param[in] i starting index
251 * \param[in] dq_pq matrix values
252 * \param[in] vec vector
253 * \param[in, out] mvec resulting vector
254 */
255 /*----------------------------------------------------------------------------*/
256
257 inline static void
_partial_matvec(const int i,const cs_sdm_t * dq_pq,const double * restrict vec,double * mvec)258 _partial_matvec(const int i,
259 const cs_sdm_t *dq_pq,
260 const double *restrict vec,
261 double *mvec)
262
263 {
264 const int n_ent = dq_pq->n_rows;
265
266 for (int irow = i; irow < n_ent; irow++) {
267 const double *restrict m_i = dq_pq->val + irow*n_ent;
268 double s = 0;
269 for (int j = 0; j < n_ent; j++) s += m_i[j]*vec[j];
270 mvec[irow] = s;
271 }
272 }
273
274 /*----------------------------------------------------------------------------*/
275 /*!
276 * \brief Initialize the geometrical quantities for EpFd Hodge operators
277 * (Cost and bubble algorithm)
278 *
279 * \param[in] cm pointer to a cs_cell_mesh_t structure
280 * \param[in, out] pq primal vector-valued quantities
281 * \param[in, out] dq dual vector-valued quantities
282 */
283 /*----------------------------------------------------------------------------*/
284
285 static inline void
_init_vb_geom_quant(const cs_cell_mesh_t * cm,cs_real_3_t * pq,cs_real_3_t * dq)286 _init_vb_geom_quant(const cs_cell_mesh_t *cm,
287 cs_real_3_t *pq,
288 cs_real_3_t *dq)
289 {
290 for (int ii = 0; ii < cm->n_ec; ii++) {
291
292 cs_nvec3_t dfq = cm->dface[ii];
293 cs_quant_t peq = cm->edge[ii];
294
295 for (int k = 0; k < 3; k++) {
296 dq[ii][k] = dfq.meas * dfq.unitv[k];
297 pq[ii][k] = peq.meas * peq.unitv[k];
298 }
299
300 } /* Loop on cell edges */
301 }
302
303 /*----------------------------------------------------------------------------*/
304 /*!
305 * \brief Initialize the geometrical quantities for FpEd Hodge operators
306 * (Cost and bubble algorithm)
307 *
308 * \param[in] cm pointer to a cs_cell_mesh_t structure
309 * \param[in, out] pq primal vector-valued quantities
310 * \param[in, out] dq dual vector-valued quantities
311 */
312 /*----------------------------------------------------------------------------*/
313
314 static inline void
_init_fb_geom_quant(const cs_cell_mesh_t * cm,cs_real_3_t * pq,cs_real_3_t * dq)315 _init_fb_geom_quant(const cs_cell_mesh_t *cm,
316 cs_real_3_t *pq,
317 cs_real_3_t *dq)
318 {
319 for (short int f = 0; f < cm->n_fc; f++) {
320
321 const cs_nvec3_t deq = cm->dedge[f];
322 const cs_quant_t pfq = cm->face[f];
323
324 for (int k = 0; k < 3; k++) {
325 dq[f][k] = deq.meas * deq.unitv[k];
326 pq[f][k] = pfq.meas * pfq.unitv[k];
327 }
328
329 } /* Loop on cell faces */
330 }
331
332 /*----------------------------------------------------------------------------*/
333 /*!
334 * \brief Initialize the local builder structure used for building Hodge op.
335 * cellwise
336 *
337 * \param[in] space_scheme discretization scheme
338 * \param[in] connect pointer to a cs_cdo_connect_t structure
339 *
340 * \return a pointer to a new allocated cs_cell_builder_t structure
341 */
342 /*----------------------------------------------------------------------------*/
343
344 static cs_cell_builder_t *
_cell_builder_create(cs_param_space_scheme_t space_scheme,const cs_cdo_connect_t * connect)345 _cell_builder_create(cs_param_space_scheme_t space_scheme,
346 const cs_cdo_connect_t *connect)
347 {
348 int size;
349
350 const int n_vc = connect->n_max_vbyc;
351 const int n_ec = connect->n_max_ebyc;
352 const int n_fc = connect->n_max_fbyc;
353
354 cs_cell_builder_t *cb = cs_cell_builder_create();
355
356 switch (space_scheme) {
357
358 case CS_SPACE_SCHEME_CDOVB:
359 size = CS_MAX(4*n_ec + 3*n_vc, n_ec*(n_ec+1));
360 BFT_MALLOC(cb->values, size, double);
361 memset(cb->values, 0, size*sizeof(cs_real_t));
362
363 size = 2*n_ec;
364 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
365 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
366 break;
367
368 case CS_SPACE_SCHEME_CDOVCB:
369 size = 2*n_vc + 3*n_ec + n_fc;
370 BFT_MALLOC(cb->values, size, double);
371 memset(cb->values, 0, size*sizeof(cs_real_t));
372
373 size = 2*n_ec + n_vc;
374 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
375 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
376 break;
377
378 case CS_SPACE_SCHEME_CDOFB:
379 size = n_fc*(n_fc+1);
380 BFT_MALLOC(cb->values, size, double);
381 memset(cb->values, 0, size*sizeof(cs_real_t));
382
383 size = 2*n_fc;
384 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
385 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
386 break;
387
388 case CS_SPACE_SCHEME_CDOEB:
389 {
390 int n_ent = CS_MAX(n_fc, n_ec);
391
392 size = n_ent*(n_ent+1);
393 BFT_MALLOC(cb->values, size, double);
394 memset(cb->values, 0, size*sizeof(cs_real_t));
395
396 size = 2*n_ent;
397 BFT_MALLOC(cb->vectors, size, cs_real_3_t);
398 memset(cb->vectors, 0, size*sizeof(cs_real_3_t));
399 }
400 break;
401
402 default:
403 bft_error(__FILE__, __LINE__, 0, _("Invalid space scheme."));
404
405 } /* End of switch on space scheme */
406
407 return cb;
408 }
409
410 /*----------------------------------------------------------------------------*/
411 /*!
412 * \brief Compute quantities used for defining the entries of the discrete
413 * Hodge for COST algo. when the property is isotropic
414 * Initialize the local discrete Hodge op. with the consistency part
415 *
416 * \param[in] cm pointer to a cs_cell_mesh_t structure
417 * \param[in] hmat pointer to a local Hodge matrix
418 * \param[in, out] sloc pointer to the local stiffness matrix
419 */
420 /*----------------------------------------------------------------------------*/
421
422 static void
_define_vb_stiffness(const cs_cell_mesh_t * cm,const cs_sdm_t * hmat,cs_sdm_t * sloc)423 _define_vb_stiffness(const cs_cell_mesh_t *cm,
424 const cs_sdm_t *hmat,
425 cs_sdm_t *sloc)
426 {
427 /* Initialize the local stiffness matrix */
428 cs_sdm_square_init(cm->n_vc, sloc);
429
430 for (int ei = 0; ei < cm->n_ec; ei++) { /* Loop on cell edges I */
431
432 const short int *v_sgn = cm->e2v_sgn + ei;
433 const short int *v_ids = cm->e2v_ids + 2*ei;
434 const short int i1 = v_ids[0], i2 = v_ids[1];
435 assert(i1 < i2);
436
437 double *si1 = sloc->val + i1*sloc->n_rows;
438 double *si2 = sloc->val + i2*sloc->n_rows;
439
440 /* Diagonal value: consistency part has already been computed */
441 const double *hii = hmat->val + ei*(1 + cm->n_ec);
442 const double dval = hii[0];
443
444 si1[i1] += dval;
445 si1[i2] -= dval;
446 si2[i2] += dval;
447
448 /* Compute extra-diag entries */
449 for (int _j = 1; _j < cm->n_ec-ei; _j++) { /* Loop on cell entities J */
450
451 const short int j1 = v_ids[2*_j], j2 = v_ids[2*_j+1];
452 assert(j1 < j2);
453
454 double *sj1 = sloc->val + j1*sloc->n_rows;
455 double *sj2 = sloc->val + j2*sloc->n_rows;
456
457 /* Add contribution from the stabilization part for each sub-volume
458 related to a primal entity */
459 const double xval = hii[_j] * v_sgn[0]*v_sgn[_j];
460
461 if (i2 < j1) { /* i1 < i2 < j1 < j2 */
462 si1[j1] += xval;
463 si1[j2] -= xval;
464 si2[j1] -= xval;
465 si2[j2] += xval;
466 }
467 else if (i2 == j1) { /* i1 < i2 = j1 < j2 */
468 si1[j1] += xval;
469 si1[j2] -= xval;
470 si2[j1] -= 2*xval;
471 si2[j2] += xval;
472 }
473 else if (i2 < j2) {
474
475 assert(i2 > j1);
476 if (i1 < j1) /* i1 < j1 < i2 < j2 */
477 si1[j1] += xval;
478 else if (i1 == j1) /* i1 = j1 < i2 < j2 */
479 si1[j1] += 2*xval;
480 else /* j1 < i1 < i2 < j2 */
481 sj1[i1] += xval;
482
483 si1[j2] -= xval;
484 sj1[i2] -= xval;
485 si2[j2] += xval;
486
487 }
488 else if (i2 == j2) {
489
490 if (i1 < j1) /* i1 < j1 < i2 = j2 */
491 si1[j1] += xval;
492 else /* j1 < i1 < i2 = j2 */
493 sj1[i1] += xval;
494
495 /* Remark: the case i1 == j1 is not possible since ei != ej */
496 si1[j2] -= xval;
497 sj1[i2] -= xval;
498 si2[j2] += 2*xval;
499
500 }
501 else { /* i2 > j2 */
502
503 if (i1 < j1) { /* i1 < j1 < j2 < i2 */
504 si1[j1] += xval;
505 si1[j2] -= xval;
506 }
507 else if (i1 == j1) { /* i1 = j1 < j2 < i2 */
508 si1[j1] += 2*xval;
509 si1[j2] -= xval;
510 }
511 else if (i1 < j2) { /* j1 < i1 < j2 < i2 */
512 sj1[i1] += xval;
513 si1[j2] -= xval;
514 }
515 else if (i1 == j2) { /* j1 < i1 = j2 < i2 */
516 sj1[i1] += xval;
517 si1[j2] -= 2*xval;
518 }
519 else { /* j1 < j2 < i1 < i2 */
520 sj1[i1] += xval;
521 sj2[i1] -= xval;
522 }
523
524 assert(i2 > j2);
525 sj1[i2] -= xval;
526 sj2[i2] += xval;
527
528 } /* End of tests */
529
530 } /* End of loop on J entities */
531
532 } /* End of loop on I entities */
533
534 /* Stiffness matrix is symmetric by construction */
535 cs_sdm_symm_ur(sloc);
536
537 }
538
539 /*----------------------------------------------------------------------------*/
540 /*!
541 * \brief Compute quantities used for defining the entries of the discrete
542 * Hodge for COST algo. when the property is isotropic
543 * Initialize the local discrete Hodge op. with the consistency part
544 *
545 * \param[in] n_ent number of local entities
546 * \param[in] invcvol 1/|c|
547 * \param[in] ptyval values of property inside this cell
548 * \param[in] pq pointer to the first set of quantities
549 * \param[in] dq pointer to the second set of quantities
550 * \param[in, out] alpha geometrical quantity
551 * \param[in, out] kappa geometrical quantity
552 * \param[in, out] hmat pointer to a cs_sdm_t struct.
553 */
554 /*----------------------------------------------------------------------------*/
555
556 static void
_compute_cost_quant_iso(const int n_ent,const double invcvol,const double ptyval,const cs_real_3_t * pq,const cs_real_3_t * dq,double * alpha,double * kappa,cs_sdm_t * hmat)557 _compute_cost_quant_iso(const int n_ent,
558 const double invcvol,
559 const double ptyval,
560 const cs_real_3_t *pq,
561 const cs_real_3_t *dq,
562 double *alpha,
563 double *kappa,
564 cs_sdm_t *hmat)
565 {
566 /* Compute several useful quantities
567 alpha_ij = delta_ij - pq_j.Consist_i where Consist_i = 1/|c| dq_i
568 qmq_ii = dq_i.ptyval.dq_i
569 kappa_i = qmq_ii / |subvol_i|
570 */
571
572 for (int i = 0; i < n_ent; i++) {
573
574 const double dsvol_i = _dp3(dq[i], pq[i]);
575
576 double *alpha_i = alpha + i*n_ent;
577 double *mi = hmat->val + i*n_ent;
578
579 alpha_i[i] = 1 - invcvol * dsvol_i;
580
581 const double qmq_ii = ptyval * _dp3(dq[i], dq[i]);
582
583 mi[i] = invcvol * qmq_ii;
584 kappa[i] = 3. * qmq_ii / dsvol_i;
585
586 for (int j = i+1; j < n_ent; j++) {
587
588 /* Initialize the upper right part of hmat with the consistency part */
589 mi[j] = invcvol * ptyval * _dp3(dq[j], dq[i]);
590
591 /* Compute the alpha matrix (not symmetric) */
592 alpha_i[j] = -invcvol * _dp3(pq[j], dq[i]);
593 alpha[j*n_ent + i] = -invcvol * _dp3(pq[i], dq[j]);
594
595 } /* Loop on entities (J) */
596
597 } /* Loop on entities (I) */
598
599 }
600
601 /*----------------------------------------------------------------------------*/
602 /*!
603 * \brief Compute quantities used for defining the entries of the discrete
604 * Hodge for COST algo.
605 * Initialize the local discrete Hodge op. with the consistency part
606 *
607 * \param[in] n_ent number of local entities
608 * \param[in] ovc 1/|c| where |c| is the cell volume
609 * \param[in] pty values of the tensor related to the material pty
610 * \param[in] pq pointer to the first set of quantities
611 * \param[in] dq pointer to the second set of quantities
612 * \param[in, out] alpha geometrical quantity
613 * \param[in, out] kappa geometrical quantity
614 * \param[in, out] hmat pointer to a cs_sdm_t structure
615 */
616 /*----------------------------------------------------------------------------*/
617
618 static void
_compute_cost_quant(const int n_ent,const double ovc,const cs_real_t pty[3][3],const cs_real_3_t * pq,const cs_real_3_t * dq,double * alpha,double * kappa,cs_sdm_t * hmat)619 _compute_cost_quant(const int n_ent,
620 const double ovc,
621 const cs_real_t pty[3][3],
622 const cs_real_3_t *pq,
623 const cs_real_3_t *dq,
624 double *alpha,
625 double *kappa,
626 cs_sdm_t *hmat)
627 {
628 /* Compute several useful quantities:
629 alpha_ij = delta_ij - 1/|c|*pq_j.dq_i
630 qmq_ii = dq_i.mat.dq_i
631 kappa_i = qmq_ii / |subvol_i| */
632
633 for (int i = 0; i < n_ent; i++) {
634
635 const double mdq_i[3]
636 = { pty[0][0] * dq[i][0] + pty[0][1] * dq[i][1] + pty[0][2] * dq[i][2],
637 pty[1][0] * dq[i][0] + pty[1][1] * dq[i][1] + pty[1][2] * dq[i][2],
638 pty[2][0] * dq[i][0] + pty[2][1] * dq[i][1] + pty[2][2] * dq[i][2] };
639
640 double *h_i = hmat->val + i*n_ent;
641
642 h_i[i] = _dp3(dq[i], mdq_i); /* Add the consistent part */
643
644 double *alpha_i = alpha + i*n_ent;
645
646 alpha_i[i] = _dp3(dq[i], pq[i]);
647 kappa[i] = 3. * h_i[i] / alpha_i[i];
648 alpha_i[i] = 1 - ovc*alpha_i[i];
649 h_i[i] *= ovc;
650
651 for (int j = i+1; j < n_ent; j++) {
652
653 /* Initialize the upper right part of hmat with the consistency part */
654 h_i[j] = ovc * _dp3(dq[j], mdq_i);
655
656 /* Compute the alpha matrix (not symmetric) */
657 alpha_i[j] = -ovc * _dp3(pq[j], dq[i]);
658
659 }
660
661 const double opq_i[3] = { -ovc*pq[i][0], -ovc*pq[i][1], -ovc*pq[i][2] };
662 for (int j = i+1; j < n_ent; j++)
663 alpha[j*n_ent + i] = _dp3(opq_i, dq[j]);
664
665 } /* Loop on entities (I) */
666
667 }
668
669 /*----------------------------------------------------------------------------*/
670 /*!
671 * \brief Compute the discrete EpFd Hodge operator (the upper right part).
672 * Co+St algo. in case of isotropic material property.
673 *
674 * \param[in] n_ent number of local entities
675 * \param[in] dbeta2 space dim * squared value of the stabilization coef.
676 * \param[in] ovc reciprocal of the cell volume
677 * \param[in] pty values of the the material pty in this cell
678 * \param[in] pq pointer to the first set of quantities
679 * \param[in] dq pointer to the second set of quantities
680 * \param[in, out] cb temporary buffers
681 * \param[in, out] hmat pointer to a cs_sdm_t structure
682 */
683 /*----------------------------------------------------------------------------*/
684
685 static void
_compute_iso_hodge_ur(const int n_ent,const double dbeta2,const double ovc,const cs_real_t pty,const cs_real_3_t * pq,const cs_real_3_t * dq,cs_cell_builder_t * cb,cs_sdm_t * hmat)686 _compute_iso_hodge_ur(const int n_ent,
687 const double dbeta2,
688 const double ovc,
689 const cs_real_t pty,
690 const cs_real_3_t *pq,
691 const cs_real_3_t *dq,
692 cs_cell_builder_t *cb,
693 cs_sdm_t *hmat)
694 {
695 const double ptyc = pty*ovc;
696
697 double *kappa = cb->values; /* size = n_ent */
698 double *kappa_pq_dqi = cb->values + n_ent; /* size = n_ent */
699 double *stab = cb->values + 2*n_ent; /* size = n_ent */
700 double *dq_pq = cb->aux->val; /* size = n_ent*n_ent */
701
702 cs_sdm_square_init(n_ent, cb->aux);
703
704 const double dbetac = dbeta2*ovc;
705 const double dbetac2 = dbetac*ovc;
706 const double beta_coef = ovc * (1 - 2*dbeta2);
707
708 /* Initialize the upper right part of the discrete Hodge op and store useful
709 quantities */
710 for (int i = 0; i < n_ent; i++) {
711
712 const double dqi[3] = {dq[i][0], dq[i][1], dq[i][2]};
713
714 double *dqi_pq = dq_pq + i*n_ent;
715 for (int j = 0; j < n_ent; j++)
716 dqi_pq[j] = _dp3(dqi, pq[j]);
717
718 const double dqi_m_dqi = pty * _dp3(dqi, dqi);
719 kappa[i] = dqi_m_dqi/dqi_pq[i];
720
721 double *hi = hmat->val + i*n_ent;
722
723 hi[i] = dqi_m_dqi*beta_coef + dbeta2*kappa[i];
724 for (int j = i+1; j < n_ent; j++)
725 hi[j] = ptyc * _dp3(dq[j], dqi);
726
727 }
728
729 /* Build the upper right part of the discrete Hodge operator. */
730 for (int i = 0; i < n_ent; i++) {
731
732 const double *dqi_pq = dq_pq + n_ent*i;
733
734 for (int k = 0; k < n_ent;k++)
735 kappa_pq_dqi[k] = kappa[k] * dqi_pq[k];
736
737 _partial_matvec(i, cb->aux, kappa_pq_dqi, stab);
738
739 double *hi = hmat->val + i*n_ent;
740
741 hi[i] += dbetac2 * stab[i];
742
743 /* Compute the extra-diagonal terms */
744 const double kai = kappa[i];
745 for (int j = i+1; j < n_ent; j++) {
746 double contrib = ovc * stab[j];
747 contrib -= kai*dq_pq[j*n_ent+i] + kappa[j]*dqi_pq[j];
748 contrib *= dbetac;
749 hi[j] += contrib;
750 }
751
752 } /* Loop on rows (entities i) */
753 }
754
755 /*----------------------------------------------------------------------------*/
756 /*!
757 * \brief Compute the discrete EpFd Hodge operator (the upper right part).
758 * Co+St algo. in case of anisotropic material property.
759 *
760 * \param[in] n_ent number of local entities
761 * \param[in] dbeta2 space dim * squared value of the stabilization coef.
762 * \param[in] vol_c cell volume
763 * \param[in] pty values of the tensor related to the material pty
764 * \param[in] pq pointer to the first set of quantities
765 * \param[in] dq pointer to the second set of quantities
766 * \param[in, out] cb temporary buffers
767 * \param[in, out] hmat pointer to a cs_sdm_t structure
768 */
769 /*----------------------------------------------------------------------------*/
770
771 static void
_compute_aniso_hodge_ur(const int n_ent,const double dbeta2,const double ovc,const cs_real_t pty[3][3],const cs_real_3_t * pq,const cs_real_3_t * dq,cs_cell_builder_t * cb,cs_sdm_t * hmat)772 _compute_aniso_hodge_ur(const int n_ent,
773 const double dbeta2,
774 const double ovc,
775 const cs_real_t pty[3][3],
776 const cs_real_3_t *pq,
777 const cs_real_3_t *dq,
778 cs_cell_builder_t *cb,
779 cs_sdm_t *hmat)
780 {
781 double *kappa = cb->values; /* size = n_ent */
782 double *kappa_pq_dqi = cb->values + n_ent; /* size = n_ent */
783 double *stab = cb->values + 2*n_ent; /* size = n_ent */
784 double *dq_pq = cb->aux->val; /* size = n_ent*n_ent */
785
786 cs_sdm_square_init(n_ent, cb->aux);
787
788 const double dbetac = dbeta2*ovc;
789 const double dbetac2 = dbetac*ovc;
790 const double beta_coef = ovc * (1 - 2*dbeta2);
791
792 /* Initialize the upper right part of the discrete Hodge op and store useful
793 quantities */
794 for (int i = 0; i < n_ent; i++) {
795
796 const double dqi[3] = {dq[i][0], dq[i][1], dq[i][2]};
797
798 double *dqi_pq = dq_pq + i*n_ent;
799 for (int j = 0; j < n_ent; j++)
800 dqi_pq[j] = _dp3(dqi, pq[j]);
801
802 const double mdqi[3]
803 = { pty[0][0] * dqi[0] + pty[0][1] * dqi[1] + pty[0][2] * dqi[2],
804 pty[1][0] * dqi[0] + pty[1][1] * dqi[1] + pty[1][2] * dqi[2],
805 pty[2][0] * dqi[0] + pty[2][1] * dqi[1] + pty[2][2] * dqi[2] };
806
807 const double dqi_m_dqi = _dp3(dqi, mdqi);
808
809 kappa[i] = dqi_m_dqi/dqi_pq[i];
810
811 double *hi = hmat->val + i*n_ent;
812
813 hi[i] = dqi_m_dqi*beta_coef + dbeta2*kappa[i];
814 for (int j = i+1; j < n_ent; j++)
815 hi[j] = ovc * _dp3(dq[j], mdqi);
816
817 }
818
819 /* Build the upper right part of the discrete Hodge operator. */
820 for (int i = 0; i < n_ent; i++) {
821
822 const double *dqi_pq = dq_pq + n_ent*i;
823
824 for (int k = 0; k < n_ent; k++)
825 kappa_pq_dqi[k] = kappa[k] * dqi_pq[k];
826
827 _partial_matvec(i, cb->aux, kappa_pq_dqi, stab);
828
829 double *hi = hmat->val + i*n_ent;
830
831 hi[i] += dbetac2 * stab[i];
832
833 /* Compute the extra-diagonal terms */
834 const double kai = kappa[i];
835 for (int j = i+1; j < n_ent; j++) {
836 double contrib = ovc*stab[j];
837 contrib -= kai*dq_pq[j*n_ent+i] + kappa[j]*dqi_pq[j];
838 contrib *= dbetac;
839 hi[j] += contrib;
840 }
841
842 } /* Loop on rows (entities i) */
843 }
844
845 /*----------------------------------------------------------------------------*/
846 /*!
847 * \brief Compute the discrete Hodge operator (the upper right part).
848 * Co+St algo. with bubble stabilization in case of isotropic
849 * material property.
850 *
851 * \param[in] n_ent number of local entities
852 * \param[in] beta the stabilization coef.
853 * \param[in] ovc reciprocal of the cell volume
854 * \param[in] pty_val value of the property inside this cell
855 * \param[in] pq pointer to the first set of quantities
856 * \param[in] dq pointer to the second set of quantities
857 * \param[in, out] cb temporary buffers
858 * \param[in, out] hmat pointer to a cs_sdm_t structure
859 */
860 /*----------------------------------------------------------------------------*/
861
862 static void
_compute_iso_bubble_hodge_ur(const int n_ent,const double beta,const double ovc,const cs_real_t pty_val,const cs_real_3_t * pq,const cs_real_3_t * dq,cs_cell_builder_t * cb,cs_sdm_t * hmat)863 _compute_iso_bubble_hodge_ur(const int n_ent,
864 const double beta,
865 const double ovc,
866 const cs_real_t pty_val,
867 const cs_real_3_t *pq,
868 const cs_real_3_t *dq,
869 cs_cell_builder_t *cb,
870 cs_sdm_t *hmat)
871 {
872 double *kappa = cb->values; /* size = n_ent */
873 double *dq_pq = cb->aux->val; /* size = n_ent*n_ent */
874
875 /* Initialize the upper right part of the discrete Hodge op and store useful
876 quantities */
877 for (int i = 0; i < n_ent; i++) {
878
879 const double dqi[3] = {dq[i][0], dq[i][1], dq[i][2]};
880 const double dqi_m_dqi = pty_val * _dp3(dqi, dqi);
881 const double dqi_pqi = _dp3(dqi, pq[i]);
882
883 kappa[i] = dqi_m_dqi/dqi_pqi;
884
885 double *dqi_pq = dq_pq + i*n_ent;
886 for (int j = 0; j < i; j++)
887 dqi_pq[j] = _dp3(dqi, pq[j])*ovc;
888 dqi_pq[i] = dqi_pqi*ovc;
889 for (int j = i+1; j < n_ent; j++)
890 dqi_pq[j] = _dp3(dqi, pq[j])*ovc;
891
892 /* Consistent part */
893 double *hi = hmat->val + i*n_ent;
894 hi[i] = dqi_m_dqi*ovc;
895 const double coef = ovc * pty_val;
896 for (int j = i+1; j < n_ent; j++)
897 hi[j] = coef * _dp3(dq[j], dqi);
898
899 }
900
901 /* Add the stabilization part */
902
903 /* \int_{p_{ec}} \theta_e*\theta_e = 0.1*|p_{ec}| and d=3 */
904 const double stab_coef = 0.3*beta*beta;
905
906 /* Build the upper right part of the discrete Hodge operator. */
907 for (int i = 0; i < n_ent; i++) {
908
909 const double *dqi_pq = dq_pq + n_ent*i;
910 double *hi = hmat->val + i*n_ent;
911
912 /* Diagonal term */
913 double stab = 0;
914 for (int k = 0; k < n_ent; k++) {
915 const double a_ik = (i == k) ? 1 - dqi_pq[k] : -dqi_pq[k];
916 stab += kappa[k] * a_ik * a_ik;
917 }
918 hi[i] += stab_coef * stab;
919
920 /* Extra-diag term */
921 for (int j = i+1; j < n_ent; j++) {
922
923 const double *dqj_pq = dq_pq + n_ent*j;
924
925 stab = 0;
926 for (int k = 0; k < n_ent; k++) {
927 const double a_ik = (i == k) ? 1 - dqi_pq[k] : -dqi_pq[k];
928 const double a_jk = (j == k) ? 1 - dqj_pq[k] : -dqj_pq[k];
929 stab += kappa[k] * a_ik * a_jk;
930 }
931 hi[j] += stab_coef * stab;
932
933 } /* Loop on row elements */
934
935 } /* Loop on partition elements */
936
937 }
938
939 /*----------------------------------------------------------------------------*/
940 /*!
941 * \brief Compute the discrete Hodge operator (the upper right part).
942 * Co+St algo. with bubble stabilization in case of anisotropic
943 * material property.
944 *
945 * \param[in] n_ent number of local entities
946 * \param[in] beta the stabilization coef.
947 * \param[in] ovc reciprocal of the cell volume
948 * \param[in] pty value of the property inside this cell
949 * \param[in] pq pointer to the first set of quantities
950 * \param[in] dq pointer to the second set of quantities
951 * \param[in, out] cb temporary buffers
952 * \param[in, out] hmat pointer to a cs_sdm_t structure
953 */
954 /*----------------------------------------------------------------------------*/
955
956 static void
_compute_aniso_bubble_hodge_ur(const int n_ent,const double beta,const double ovc,const cs_real_t pty[3][3],const cs_real_3_t * pq,const cs_real_3_t * dq,cs_cell_builder_t * cb,cs_sdm_t * hmat)957 _compute_aniso_bubble_hodge_ur(const int n_ent,
958 const double beta,
959 const double ovc,
960 const cs_real_t pty[3][3],
961 const cs_real_3_t *pq,
962 const cs_real_3_t *dq,
963 cs_cell_builder_t *cb,
964 cs_sdm_t *hmat)
965 {
966 double *kappa = cb->values; /* size = n_ent */
967 double *dq_pq = cb->aux->val; /* size = n_ent*n_ent */
968
969 /* Initialize the upper right part of the discrete Hodge op and store useful
970 quantities */
971 for (int i = 0; i < n_ent; i++) {
972
973 const double dqi[3] = {dq[i][0], dq[i][1], dq[i][2]};
974 const double mdqi[3]
975 = { pty[0][0] * dqi[0] + pty[0][1] * dqi[1] + pty[0][2] * dqi[2],
976 pty[1][0] * dqi[0] + pty[1][1] * dqi[1] + pty[1][2] * dqi[2],
977 pty[2][0] * dqi[0] + pty[2][1] * dqi[1] + pty[2][2] * dqi[2] };
978 const double dqi_m_dqi = _dp3(dqi, mdqi);
979 const double dqi_pqi = _dp3(dqi, pq[i]);
980
981 kappa[i] = dqi_m_dqi/dqi_pqi;
982
983 double *dqi_pq = dq_pq + i*n_ent;
984 for (int j = 0; j < i; j++)
985 dqi_pq[j] = _dp3(dqi, pq[j])*ovc;
986 dqi_pq[i] = dqi_pqi*ovc;
987 for (int j = i+1; j < n_ent; j++)
988 dqi_pq[j] = _dp3(dqi, pq[j])*ovc;
989
990 /* Consistent part */
991 double *hi = hmat->val + i*n_ent;
992 hi[i] = dqi_m_dqi*ovc;
993 for (int j = i+1; j < n_ent; j++)
994 hi[j] = ovc * _dp3(dq[j], mdqi);
995
996 }
997
998 /* Add the stabilization part */
999
1000 /* \int_{p_{ec}} \theta_e*\theta_e = 0.1*|p_{ec}| and d=3 */
1001 const double stab_coef = 0.3*beta*beta;
1002
1003 /* Build the upper right part of the discrete Hodge operator. */
1004 for (int i = 0; i < n_ent; i++) {
1005
1006 const double *dqi_pq = dq_pq + n_ent*i;
1007 double *hi = hmat->val + i*n_ent;
1008
1009 /* Diagonal term */
1010 double stab = 0;
1011 for (int k = 0; k < n_ent; k++) {
1012 const double a_ik = (i == k) ? 1 - dqi_pq[k] : -dqi_pq[k];
1013 stab += kappa[k] * a_ik * a_ik;
1014 }
1015 hi[i] += stab_coef * stab;
1016
1017 /* Extra-diag term */
1018 for (int j = i+1; j < n_ent; j++) {
1019
1020 const double *dqj_pq = dq_pq + n_ent*j;
1021
1022 stab = 0;
1023 for (int k = 0; k < n_ent; k++) {
1024 const double a_ik = (i == k) ? 1 - dqi_pq[k] : -dqi_pq[k];
1025 const double a_jk = (j == k) ? 1 - dqj_pq[k] : -dqj_pq[k];
1026 stab += kappa[k] * a_ik * a_jk;
1027 }
1028 hi[j] += stab_coef * stab;
1029
1030 } /* Loop on row elements */
1031
1032 } /* Loop on partition elements */
1033
1034 }
1035
1036 /*----------------------------------------------------------------------------*/
1037 /*!
1038 * \brief Compute the discrete EpFd Hodge operator (the upper right part)
1039 * from primal edges to dual faces with the algorithm called:
1040 * Orthogonal Consistent/Sub-Stabilization decomposition (OCS2) with a
1041 * subdivision of pvol_{e,c}
1042 * Case of anisotropic material property.
1043 *
1044 * \param[in] dbeta2 space dim * squared value of the stabilization coef.
1045 * \param[in] pty values of the tensor related to the material pty
1046 * \param[in] cm pointer to a cs_cell_mesh_t structure
1047 * \param[in, out] cb temporary buffers
1048 * \param[in, out] hmat pointer to a cs_sdm_t structure
1049 */
1050 /*----------------------------------------------------------------------------*/
1051
1052 static void
_compute_aniso_hepfd_ocs2_ur(const double dbeta2,const cs_real_t pty[3][3],const cs_cell_mesh_t * cm,cs_cell_builder_t * cb,cs_sdm_t * hmat)1053 _compute_aniso_hepfd_ocs2_ur(const double dbeta2,
1054 const cs_real_t pty[3][3],
1055 const cs_cell_mesh_t *cm,
1056 cs_cell_builder_t *cb,
1057 cs_sdm_t *hmat)
1058 {
1059 const double ovc = 1./cm->vol_c;
1060 const int n_ent = cm->n_ec;
1061
1062 /* Store the consistent part of the reconstruction of the basis element */
1063 cs_real_3_t *consist = cb->vectors;
1064 for (int i = 0; i < n_ent; i++) {
1065 const double fd_coef = ovc * cm->dface[i].meas;
1066 for (int k = 0; k < 3; k++)
1067 consist[i][k] = fd_coef * cm->dface[i].unitv[k];
1068 }
1069
1070 /* Initialize the upper right part of the discrete Hodge op and store useful
1071 quantities */
1072
1073 /* Consistency part */
1074 for (int i = 0; i < n_ent; i++) {
1075
1076 double pty_fd_i[3] = {0, 0, 0};
1077 for (int k = 0; k < 3; k++) {
1078 pty_fd_i[0] += pty[0][k] * cm->dface[i].unitv[k];
1079 pty_fd_i[1] += pty[1][k] * cm->dface[i].unitv[k];
1080 pty_fd_i[2] += pty[2][k] * cm->dface[i].unitv[k];
1081 }
1082 for (int k = 0; k < 3; k++)
1083 pty_fd_i[k] *= cm->dface[i].meas;
1084
1085 double *h_i = hmat->val + i*n_ent;
1086 for (int j = i; j < n_ent; j++)
1087 h_i[j] = _dp3(consist[j], pty_fd_i);
1088
1089 }
1090
1091 /* Compute the contribution part of each edge for the stabilization term */
1092 cs_real_t *contrib_pe = cb->values;
1093 memset(contrib_pe, 0, n_ent*sizeof(cs_real_t));
1094
1095 for (int f = 0; f < cm->n_fc; f++) {
1096
1097 for (int i = cm->f2e_idx[f]; f < cm->f2e_idx[f+1]; f++) {
1098
1099 const short int k = cm->f2e_ids[i]; /* edge k */
1100 const cs_real_t ep_k[3] = { cm->edge[k].meas * cm->edge[k].unitv[0],
1101 cm->edge[k].meas * cm->edge[k].unitv[1],
1102 cm->edge[k].meas * cm->edge[k].unitv[2] };
1103
1104 cs_real_3_t pty_fdk = {0, 0, 0};
1105 for (int kk = 0; kk < 3; kk++) {
1106 pty_fdk[0] += pty[0][kk] * cm->sefc[i].unitv[kk];
1107 pty_fdk[1] += pty[1][kk] * cm->sefc[i].unitv[kk];
1108 pty_fdk[2] += pty[2][kk] * cm->sefc[i].unitv[kk];
1109 }
1110
1111 const cs_real_t coef = cm->sefc[i].meas * dbeta2;
1112 contrib_pe[k] +=
1113 coef * _dp3(cm->sefc[i].unitv, pty_fdk)/_dp3(cm->sefc[i].unitv, ep_k) ;
1114
1115 } /* Loop on face edges */
1116
1117 } /* Loop on cell faces */
1118
1119 /* Stabilization part */
1120 for (int k = 0; k < n_ent; k++) {
1121
1122 const cs_real_t contrib_pek = contrib_pe[k];
1123 const cs_real_t ep_k[3] = { cm->edge[k].meas * cm->edge[k].unitv[0],
1124 cm->edge[k].meas * cm->edge[k].unitv[1],
1125 cm->edge[k].meas * cm->edge[k].unitv[2] };
1126
1127
1128 for (int i = 0; i < n_ent; i++) {
1129
1130 double contrib_i = -_dp3(consist[i], ep_k);
1131 if (i == k) contrib_i += 1;
1132 const double contrib_ik = contrib_pek * contrib_i;
1133
1134 double *h_i = hmat->val + i*n_ent;
1135 h_i[i] += contrib_ik * contrib_i;
1136
1137 for (int j = i+1; j < n_ent; j++) {
1138 if (j != k)
1139 h_i[j] += - _dp3(consist[j], ep_k) * contrib_ik;
1140 else
1141 h_i[j] += (1 -_dp3(consist[j], ep_k)) * contrib_ik;
1142 }
1143
1144 } /* i */
1145
1146 } /* Loop on sub-volume k */
1147
1148 }
1149
1150 /*----------------------------------------------------------------------------*/
1151 /*!
1152 * \brief Build a local discrete Hodge operator using the generic COST algo.
1153 * and a cellwise view of the mesh. Specific for EpFd Hodge operator.
1154 * COST means COnsistency + STabilization
1155 *
1156 * \param[in] n_ent number of local entities
1157 * \param[in] beta2 squared value of the stabilization coefficient
1158 * \param[in] alpha precomputed quantities
1159 * \param[in] kappa precomputed quantities
1160 * \param[in, out] hval values of the coefficient of the Hodge matrix
1161 */
1162 /*----------------------------------------------------------------------------*/
1163
1164 static void
_compute_hodge_cost(const int n_ent,const double beta2,const double alpha[],const double kappa[],double hval[])1165 _compute_hodge_cost(const int n_ent,
1166 const double beta2,
1167 const double alpha[],
1168 const double kappa[],
1169 double hval[])
1170 {
1171 for (int i = 0; i < n_ent; i++) { /* Loop on cell entities I */
1172
1173 const int shift_i = i*n_ent;
1174 const double *alpha_i = alpha + shift_i;
1175
1176 double *mi = hval + shift_i;
1177
1178 /* Add contribution from the stabilization part for
1179 each sub-volume related to a primal entity */
1180 double stab_part = 0;
1181 for (int k = 0; k < n_ent; k++) /* Loop over sub-volumes */
1182 stab_part += kappa[k] * alpha_i[k] * alpha_i[k];
1183
1184 mi[i] += beta2 * stab_part; /* Consistency part has already been computed */
1185
1186 /* Compute extra-diag entries */
1187 for (int j = i + 1; j < n_ent; j++) { /* Loop on cell entities J */
1188
1189 const int shift_j = j*n_ent;
1190 const double *alpha_j = alpha + shift_j;
1191 double *mj = hval + shift_j;
1192
1193 /* Add contribution from the stabilization part for
1194 each sub-volume related to a primal entity */
1195 stab_part = 0;
1196 for (int k = 0; k < n_ent; k++) /* Loop over sub-volumes */
1197 stab_part += kappa[k] * alpha_i[k] * alpha_j[k];
1198
1199 mi[j] += beta2 * stab_part;
1200 mj[i] = mi[j]; /* Symmetric by construction */
1201
1202 } /* End of loop on J entities */
1203
1204 } /* End of loop on I entities */
1205
1206 }
1207
1208 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1209
1210 /*============================================================================
1211 * Public function prototypes
1212 *============================================================================*/
1213
1214 /*----------------------------------------------------------------------------*/
1215 /*!
1216 * \brief Create and initialize a pointer to a cs_hodge_t structure
1217 *
1218 * \param[in] connect pointer to cs_cdo_connect_t structure
1219 * \param[in] property pointer to a property structure
1220 * \param[in] hp pointer to a cs_hodge_param_t structure
1221 * \param[in] need_tensor true if one needs a tensor otherwise false
1222 * \param[in] need_eigen true if one needs to compute eigen valuese
1223 *
1224 * \return a pointer to the new allocated cs_hodge_t structure
1225 */
1226 /*----------------------------------------------------------------------------*/
1227
1228 cs_hodge_t *
cs_hodge_create(const cs_cdo_connect_t * connect,const cs_property_t * property,const cs_hodge_param_t * hp,bool need_tensor,bool need_eigen)1229 cs_hodge_create(const cs_cdo_connect_t *connect,
1230 const cs_property_t *property,
1231 const cs_hodge_param_t *hp,
1232 bool need_tensor,
1233 bool need_eigen)
1234 {
1235 cs_hodge_t *hdg = NULL;
1236
1237 BFT_MALLOC(hdg, 1, cs_hodge_t);
1238
1239 hdg->param = hp;
1240
1241 switch (hp->type) {
1242
1243 case CS_HODGE_TYPE_VPCD:
1244 hdg->matrix = cs_sdm_square_create(connect->n_max_vbyc);
1245 break;
1246 case CS_HODGE_TYPE_EPFD:
1247 hdg->matrix = cs_sdm_square_create(connect->n_max_ebyc);
1248 break;
1249 case CS_HODGE_TYPE_FPED:
1250 case CS_HODGE_TYPE_EDFP:
1251 hdg->matrix = cs_sdm_square_create(connect->n_max_fbyc);
1252 break;
1253 case CS_HODGE_TYPE_CPVD:
1254 hdg->matrix = cs_sdm_square_create(1);
1255 break;
1256 case CS_HODGE_TYPE_FB:
1257 hdg->matrix = cs_sdm_square_create(connect->n_max_fbyc + 1);
1258 break;
1259 case CS_HODGE_TYPE_VC:
1260 hdg->matrix = cs_sdm_square_create(connect->n_max_vbyc + 1);
1261 break;
1262
1263 default:
1264 hdg->matrix = NULL;
1265 break;
1266 }
1267
1268 BFT_MALLOC(hdg->pty_data, 1, cs_property_data_t);
1269 cs_property_data_init(need_tensor, need_eigen, property, hdg->pty_data);
1270 if (hdg->pty_data->is_unity == false && connect->n_cells > 0)
1271 cs_hodge_set_property_value(0, 0, 0, hdg);
1272
1273 return hdg;
1274 }
1275
1276 /*----------------------------------------------------------------------------*/
1277 /*!
1278 * \brief Create and initialize an array of pointers to a cs_hodge_t
1279 * structures. This array is of size the number of OpenMP threads.
1280 * Only the one associated to the current thread is set.
1281 *
1282 * \param[in] connect pointer to cs_cdo_connect_t structure
1283 * \param[in] property pointer to a property structure
1284 * \param[in] hp pointer to a cs_hodge_param_t structure
1285 * \param[in] need_tensor true if one needs a tensor otherwise false
1286 * \param[in] need_eigen true if one needs to compute eigen valuese
1287 *
1288 * \return an array of pointers of cs_hodge_t structures
1289 */
1290 /*----------------------------------------------------------------------------*/
1291
1292 cs_hodge_t **
cs_hodge_init_context(const cs_cdo_connect_t * connect,const cs_property_t * property,const cs_hodge_param_t * hp,bool need_tensor,bool need_eigen)1293 cs_hodge_init_context(const cs_cdo_connect_t *connect,
1294 const cs_property_t *property,
1295 const cs_hodge_param_t *hp,
1296 bool need_tensor,
1297 bool need_eigen)
1298 {
1299 cs_hodge_t **hodge_array = NULL;
1300
1301 BFT_MALLOC(hodge_array, cs_glob_n_threads, cs_hodge_t *);
1302 for (int i = 0; i < cs_glob_n_threads; i++)
1303 hodge_array[i] = NULL;
1304
1305 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
1306 #pragma omp parallel
1307 {
1308 int t_id = omp_get_thread_num();
1309 assert(t_id < cs_glob_n_threads);
1310
1311 hodge_array[t_id] = cs_hodge_create(connect, property, hp,
1312 need_tensor, need_eigen);
1313 }
1314 #else
1315 assert(cs_glob_n_threads == 1);
1316 hodge_array[0] = cs_hodge_create(connect, property, hp,
1317 need_tensor, need_eigen);
1318 #endif /* openMP */
1319
1320 return hodge_array;
1321 }
1322
1323 /*----------------------------------------------------------------------------*/
1324 /*!
1325 * \brief Free a cs_hodge_t structure
1326 *
1327 * \param[in, out] p_hodge double pointer to a cs_hodge_t structure
1328 */
1329 /*----------------------------------------------------------------------------*/
1330
1331 void
cs_hodge_free(cs_hodge_t ** p_hodge)1332 cs_hodge_free(cs_hodge_t **p_hodge)
1333 {
1334 cs_hodge_t *hdg = *p_hodge;
1335
1336 if (hdg == NULL)
1337 return;
1338
1339 /* Other pointers are shared */
1340 hdg->matrix = cs_sdm_free(hdg->matrix);
1341
1342 BFT_FREE(hdg->pty_data);
1343
1344 BFT_FREE(hdg);
1345 *p_hodge = NULL;
1346 }
1347
1348 /*----------------------------------------------------------------------------*/
1349 /*!
1350 * \brief Free a set of cs_hodge_t structures
1351 *
1352 * \param[in, out] p_hodges triple pointer to cs_hodge_t structures
1353 */
1354 /*----------------------------------------------------------------------------*/
1355
1356 void
cs_hodge_free_context(cs_hodge_t *** p_hodges)1357 cs_hodge_free_context(cs_hodge_t ***p_hodges)
1358 {
1359 cs_hodge_t **hodge_array = *p_hodges;
1360
1361 if (hodge_array == NULL)
1362 return;
1363
1364 #if defined(HAVE_OPENMP)
1365 # pragma omp parallel
1366 {
1367 int t_id = omp_get_thread_num();
1368 cs_hodge_free(&(hodge_array[t_id]));
1369 }
1370 #else /* No OpenMP */
1371 cs_hodge_free(&(hodge_array[0]));
1372 #endif /* openMP */
1373
1374 BFT_FREE(hodge_array);
1375 *p_hodges = NULL;
1376 }
1377
1378 /*----------------------------------------------------------------------------*/
1379 /*!
1380 * \brief Retrieve a function pointer to compute a discrete Hodge operator
1381 *
1382 * \param[in] calling_func name of the calling function
1383 * \param[in] hp a cs_hodge_param_t structure
1384 *
1385 * \return a pointer to the corresponding function
1386 */
1387 /*----------------------------------------------------------------------------*/
1388
1389 cs_hodge_compute_t *
cs_hodge_get_func(const char * calling_func,const cs_hodge_param_t hp)1390 cs_hodge_get_func(const char *calling_func,
1391 const cs_hodge_param_t hp)
1392 {
1393 cs_hodge_compute_t *hdg_func = NULL;
1394
1395 switch (hp.type) {
1396
1397 case CS_HODGE_TYPE_VPCD:
1398 {
1399 switch (hp.algo) {
1400
1401 case CS_HODGE_ALGO_COST:
1402 case CS_HODGE_ALGO_OCS2:
1403 case CS_HODGE_ALGO_BUBBLE:
1404 case CS_HODGE_ALGO_VORONOI:
1405 return cs_hodge_vpcd_voro_get;
1406 case CS_HODGE_ALGO_WBS:
1407 return cs_hodge_vpcd_wbs_get;
1408
1409 default:
1410 bft_error(__FILE__, __LINE__, 0,
1411 "%s: Invalid algorithm to compute a Vp-Cd Hodge operator.\n"
1412 " The calling function is %s\n", __func__, calling_func);
1413 break;
1414 }
1415 }
1416 break;
1417
1418 case CS_HODGE_TYPE_EPFD:
1419 {
1420 switch (hp.algo) {
1421
1422 case CS_HODGE_ALGO_COST:
1423 return cs_hodge_epfd_cost_get;
1424 case CS_HODGE_ALGO_OCS2:
1425 return cs_hodge_epfd_ocs2_get;
1426 case CS_HODGE_ALGO_BUBBLE:
1427 case CS_HODGE_ALGO_WBS: /* By default, one should define a specific
1428 algorithm for WBS */
1429 return cs_hodge_epfd_bubble_get;
1430 case CS_HODGE_ALGO_VORONOI:
1431 return cs_hodge_epfd_voro_get;
1432
1433 default:
1434 bft_error(__FILE__, __LINE__, 0,
1435 "%s: Invalid algorithm to compute a Ep-Fd Hodge operator.\n"
1436 " The calling function is %s\n", __func__, calling_func);
1437 break;
1438 }
1439 }
1440 break;
1441
1442 case CS_HODGE_TYPE_EDFP:
1443 {
1444 switch (hp.algo) {
1445
1446 case CS_HODGE_ALGO_COST:
1447 return cs_hodge_edfp_cost_get_opt;
1448 case CS_HODGE_ALGO_BUBBLE:
1449 case CS_HODGE_ALGO_WBS: /* By default, one should define a specific
1450 algorithm for WBS */
1451 return cs_hodge_edfp_bubble_get;
1452 case CS_HODGE_ALGO_VORONOI:
1453 return cs_hodge_edfp_voro_get;
1454
1455 default:
1456 bft_error(__FILE__, __LINE__, 0,
1457 "%s: Invalid algorithm to compute a Ed-Fp Hodge operator.\n"
1458 " The calling function is %s\n", __func__, calling_func);
1459 break;
1460 }
1461 }
1462 break;
1463
1464 case CS_HODGE_TYPE_FPED:
1465 {
1466 switch (hp.algo) {
1467
1468 case CS_HODGE_ALGO_COST:
1469 return cs_hodge_fped_cost_get;
1470 case CS_HODGE_ALGO_BUBBLE:
1471 case CS_HODGE_ALGO_WBS: /* By default, one should define a specific
1472 algorithm for WBS */
1473 return cs_hodge_fped_bubble_get;
1474 case CS_HODGE_ALGO_VORONOI:
1475 return cs_hodge_fped_voro_get;
1476
1477 default:
1478 bft_error(__FILE__, __LINE__, 0,
1479 "%s: Invalid algorithm to compute a Fp-Ed Hodge operator.\n"
1480 " The calling function is %s\n", __func__, calling_func);
1481 break;
1482 }
1483 }
1484 break;
1485
1486 case CS_HODGE_TYPE_FB:
1487 return cs_hodge_fb_get;
1488
1489 case CS_HODGE_TYPE_VC:
1490 switch (hp.algo) {
1491
1492 case CS_HODGE_ALGO_VORONOI:
1493 return cs_hodge_vcb_voro_get;
1494 case CS_HODGE_ALGO_WBS:
1495 return cs_hodge_vcb_wbs_get;
1496 default:
1497 bft_error(__FILE__, __LINE__, 0,
1498 "%s: Invalid algorithm to compute a Fp-Ed Hodge operator.\n"
1499 " The calling function is %s\n", __func__, calling_func);
1500 break;
1501 }
1502 break;
1503
1504 default:
1505 bft_error(__FILE__, __LINE__, 0,
1506 " %s: Invalid type of Hodge operator called by %s\n",
1507 __func__, calling_func);
1508 break;
1509
1510 } /* Switch on the type of Hodge operator */
1511
1512 return hdg_func;
1513 }
1514
1515 /*----------------------------------------------------------------------------*/
1516 /*!
1517 * \brief Check the consistency of the settings between terms related to a
1518 * mass matrix and define the common algorithm to use.
1519 * If a term should not be considered, set the algorithm to
1520 * CS_HODGE_N_ALGOS
1521 *
1522 * \param[in] eqname name of the equation to check
1523 * \param[in] reac_algo optional algo. used for the reaction term
1524 * \param[in] time_algo optional algo. used for the unsteady term
1525 * \param[in] srct_algo optional algo. used for the source term
1526 *
1527 * \return the common algorithm to use
1528 */
1529 /*----------------------------------------------------------------------------*/
1530
1531 cs_hodge_algo_t
cs_hodge_set_mass_algo(const char * eqname,cs_hodge_algo_t reac_algo,cs_hodge_algo_t time_algo,cs_hodge_algo_t srct_algo)1532 cs_hodge_set_mass_algo(const char *eqname,
1533 cs_hodge_algo_t reac_algo,
1534 cs_hodge_algo_t time_algo,
1535 cs_hodge_algo_t srct_algo)
1536 {
1537 cs_hodge_algo_t return_algo = CS_HODGE_ALGO_VORONOI;
1538
1539 if (reac_algo != CS_HODGE_N_ALGOS) { /* Hodge algo. is set for reaction */
1540
1541 return_algo = reac_algo;
1542
1543 if (time_algo != CS_HODGE_N_ALGOS) {
1544
1545 if (reac_algo != time_algo)
1546 bft_error(__FILE__, __LINE__, 0,
1547 " %s: The configuration of the Hodge algorithm between the"
1548 " reaction and unsteady term is not consistent.\n"
1549 " Please check your settings for equation \"%s\"\n",
1550 __func__, eqname);
1551
1552 if (srct_algo != CS_HODGE_N_ALGOS)
1553 if (time_algo != srct_algo)
1554 bft_error(__FILE__, __LINE__, 0,
1555 " %s: The configuration of the Hodge algorithm between the"
1556 " source term and unsteady term is not consistent.\n"
1557 " Please check your settings for equation \"%s\"\n",
1558 __func__, eqname);
1559
1560 }
1561 else { /* Hodge algo not set for the unsteady term */
1562
1563 if (srct_algo != CS_HODGE_N_ALGOS)
1564 if (reac_algo != srct_algo)
1565 bft_error(__FILE__, __LINE__, 0,
1566 " %s: The configuration of the Hodge algorithm between the"
1567 " reaction and source term is not consistent.\n"
1568 " Please check your settings for equation \"%s\"\n",
1569 __func__, eqname);
1570
1571 }
1572
1573 }
1574 else { /* Hodge algo not set for the reaction term */
1575
1576 if (time_algo != CS_HODGE_N_ALGOS) {
1577
1578 return_algo = time_algo;
1579
1580 if (srct_algo != CS_HODGE_N_ALGOS)
1581 if (time_algo != srct_algo)
1582 bft_error(__FILE__, __LINE__, 0,
1583 " %s: The configuration of the Hodge algorithm between the"
1584 " source term and unsteady term is not consistent.\n"
1585 " Please check your settings for equation \"%s\"\n",
1586 __func__, eqname);
1587
1588 }
1589 else { /* Neither time_algo nor reac_algo is set */
1590
1591 if (srct_algo != CS_HODGE_N_ALGOS)
1592 return_algo = srct_algo;
1593
1594 }
1595
1596 }
1597
1598 return return_algo;
1599 }
1600
1601 /*----------------------------------------------------------------------------*/
1602 /*!
1603 * \brief Output the settings related to a cs_hodge_param_t structure
1604 *
1605 * \param[in] prefix optional string
1606 * \param[in] property optional pointer to a property structure
1607 * \param[in] hp a cs_hodge_param_t structure
1608 */
1609 /*----------------------------------------------------------------------------*/
1610
1611 void
cs_hodge_param_log(const char * prefix,const cs_property_t * property,const cs_hodge_param_t hp)1612 cs_hodge_param_log(const char *prefix,
1613 const cs_property_t *property,
1614 const cs_hodge_param_t hp)
1615 {
1616 const char *_p;
1617 const char _empty_prefix[2] = "";
1618 if (prefix == NULL)
1619 _p = _empty_prefix;
1620 else
1621 _p = prefix;
1622
1623 cs_log_printf(CS_LOG_SETUP, "%s | Type: %s\n",
1624 _p, cs_hodge_type_desc[hp.type]);
1625 cs_log_printf(CS_LOG_SETUP, "%s | Algo: %s\n",
1626 _p, cs_hodge_algo_desc[hp.algo]);
1627 if (hp.algo == CS_HODGE_ALGO_COST ||
1628 hp.algo == CS_HODGE_ALGO_OCS2 ||
1629 hp.algo == CS_HODGE_ALGO_BUBBLE)
1630 cs_log_printf(CS_LOG_SETUP, "%s | Algo.Coef: %.3e\n",
1631 _p, hp.coef);
1632
1633 if (property != NULL)
1634 cs_log_printf(CS_LOG_SETUP, "%s | Associated property: %s\n",
1635 _p, cs_property_get_name(property));
1636 cs_log_printf(CS_LOG_SETUP, "%s | Property inversion: %s\n",
1637 _p, cs_base_strtf(hp.inv_pty));
1638 }
1639
1640 /*----------------------------------------------------------------------------*/
1641 /*!
1642 * \brief Copy the set of parameters associated to a discrete Hodge operator
1643 * to another one
1644 *
1645 * \param[in] h_ref reference set of parameters
1646 * \param[in, out] h_cpy set of parameters to update
1647 */
1648 /*----------------------------------------------------------------------------*/
1649
1650 void
cs_hodge_copy_parameters(const cs_hodge_param_t * h_ref,cs_hodge_param_t * h_cpy)1651 cs_hodge_copy_parameters(const cs_hodge_param_t *h_ref,
1652 cs_hodge_param_t *h_cpy)
1653 {
1654 if (h_ref == NULL || h_cpy == NULL)
1655 return;
1656
1657 h_cpy->inv_pty = h_ref->inv_pty;
1658 h_cpy->type = h_ref->type;
1659 h_cpy->algo = h_ref->algo;
1660 h_cpy->coef = h_ref->coef;
1661 }
1662
1663 /*----------------------------------------------------------------------------*/
1664 /*!
1665 * \brief Set the property value (scalar- or tensor-valued) related to a
1666 * discrete Hodge operator inside a cell and if needed other related
1667 * quantities
1668 *
1669 * \param[in] c_id id of the cell to deal with
1670 * \param[in] t_eval time at which one performs the evaluation
1671 * \param[in] c_flag flag related to this cell
1672 * \param[in, out] hodge pointer to a cs_hodge_t structure
1673 */
1674 /*----------------------------------------------------------------------------*/
1675
1676 void
cs_hodge_set_property_value(const cs_lnum_t c_id,const cs_real_t t_eval,const cs_flag_t c_flag,cs_hodge_t * hodge)1677 cs_hodge_set_property_value(const cs_lnum_t c_id,
1678 const cs_real_t t_eval,
1679 const cs_flag_t c_flag,
1680 cs_hodge_t *hodge)
1681 {
1682 assert(hodge != NULL);
1683
1684 cs_property_data_t *ptyd = hodge->pty_data;
1685
1686 if (ptyd->property == NULL)
1687 return; /* The default initialization corresponds to what is needed */
1688
1689 if (ptyd->need_tensor) {
1690
1691 cs_property_get_cell_tensor(c_id,
1692 t_eval,
1693 ptyd->property,
1694 hodge->param->inv_pty,
1695 ptyd->tensor);
1696
1697 if (ptyd->is_iso)
1698 ptyd->value = ptyd->tensor[0][0];
1699
1700 }
1701 else {
1702
1703 if (ptyd->is_iso) {
1704
1705 ptyd->value = cs_property_get_cell_value(c_id,
1706 t_eval,
1707 ptyd->property);
1708
1709 if (hodge->param->inv_pty) {
1710 assert(fabs(ptyd->value) > FLT_MIN);
1711 ptyd->value = 1/ptyd->value;
1712 }
1713
1714 }
1715 else { /* Anisotropic so a tensor is needed */
1716
1717 ptyd->need_tensor = true;
1718 cs_property_get_cell_tensor(c_id,
1719 t_eval,
1720 ptyd->property,
1721 hodge->param->inv_pty,
1722 ptyd->tensor);
1723
1724 }
1725
1726 } /* Tensor-valued evaluation of the property is not needed */
1727
1728 /* Set additional quantities in case of more advanced way of enforcing the
1729 essential BCs */
1730
1731 if (c_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
1732 if (ptyd->need_eigen) {
1733
1734 if (ptyd->need_tensor)
1735 cs_math_33_eigen(ptyd->tensor,
1736 &(ptyd->eigen_ratio), &(ptyd->eigen_max));
1737 else
1738 ptyd->eigen_ratio = 1.0, ptyd->eigen_max = ptyd->value;
1739
1740 }
1741 }
1742
1743 }
1744
1745 /*----------------------------------------------------------------------------*/
1746 /*!
1747 * \brief Set the property value (scalar- or tensor-valued) related to a
1748 * discrete Hodge operator inside a cell and if needed ohter related
1749 * quantities.
1750 * Cell-wise variant (usage of cs_cell_mesh_t structure)
1751 *
1752 * \param[in] cm pointer to a cs_cell_mesh_t structure
1753 * \param[in] t_eval time at which one performs the evaluation
1754 * \param[in] c_flag flag related to this cell
1755 * \param[in, out] hodge pointer to a cs_hodge_t structure
1756 */
1757 /*----------------------------------------------------------------------------*/
1758
1759 void
cs_hodge_set_property_value_cw(const cs_cell_mesh_t * cm,const cs_real_t t_eval,const cs_flag_t c_flag,cs_hodge_t * hodge)1760 cs_hodge_set_property_value_cw(const cs_cell_mesh_t *cm,
1761 const cs_real_t t_eval,
1762 const cs_flag_t c_flag,
1763 cs_hodge_t *hodge)
1764 {
1765 assert(hodge != NULL);
1766
1767 cs_property_data_t *ptyd = hodge->pty_data;
1768
1769 if (ptyd->property == NULL)
1770 return; /* The default initialization corresponds to what is needed */
1771
1772 if (ptyd->need_tensor) {
1773
1774 cs_property_tensor_in_cell(cm,
1775 ptyd->property,
1776 t_eval,
1777 hodge->param->inv_pty,
1778 ptyd->tensor);
1779
1780 if (ptyd->is_iso)
1781 ptyd->value = ptyd->tensor[0][0];
1782
1783 }
1784 else {
1785
1786 if (ptyd->is_iso) {
1787
1788 ptyd->value = cs_property_value_in_cell(cm, ptyd->property, t_eval);
1789
1790 if (hodge->param->inv_pty) {
1791 assert(fabs(ptyd->value) > FLT_MIN);
1792 ptyd->value = 1/ptyd->value;
1793 }
1794
1795 }
1796 else { /* Anisotropic so a tensor is needed */
1797
1798 ptyd->need_tensor = true;
1799 cs_property_tensor_in_cell(cm,
1800 ptyd->property,
1801 t_eval,
1802 hodge->param->inv_pty,
1803 ptyd->tensor);
1804
1805 }
1806
1807 } /* Tensor-valued evaluation of the property is not needed */
1808
1809 /* Set additional quantities in case of more advanced way of enforcing the
1810 essential BCs */
1811
1812 if (c_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) {
1813 if (ptyd->need_eigen) {
1814
1815 if (ptyd->need_tensor)
1816 cs_math_33_eigen(ptyd->tensor,
1817 &(ptyd->eigen_ratio), &(ptyd->eigen_max));
1818 else
1819 ptyd->eigen_ratio = 1.0, ptyd->eigen_max = ptyd->value;
1820
1821 }
1822 }
1823
1824 }
1825
1826 /*----------------------------------------------------------------------------*/
1827 /*!
1828 * \brief Build a local stiffness matrix using the generic COST algo.
1829 * The computed matrix is stored in cb->loc and the related discrete
1830 * hodge operator in hodge->matrix
1831 * Case of CDO face-based schemes
1832 *
1833 * \param[in] cm pointer to a cs_cell_mesh_t structure
1834 * \param[in, out] hodge pointer to a cs_hodge_t structure
1835 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1836 */
1837 /*----------------------------------------------------------------------------*/
1838
1839 void
cs_hodge_fb_cost_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)1840 cs_hodge_fb_cost_get_stiffness(const cs_cell_mesh_t *cm,
1841 cs_hodge_t *hodge,
1842 cs_cell_builder_t *cb)
1843 {
1844 /* Sanity checks */
1845 assert(cs_eflag_test(cm->flag,
1846 CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1847
1848 /* Compute the local discrete Hodge operator */
1849 cs_hodge_edfp_cost_get_opt(cm, hodge, cb);
1850
1851 const cs_sdm_t *hmat = hodge->matrix;
1852
1853 /* Initialize the local stiffness matrix */
1854 cs_sdm_t *sloc = cb->loc;
1855 cs_sdm_square_init(cm->n_fc + 1, sloc);
1856
1857 double *sval_crow = sloc->val + cm->n_fc*sloc->n_rows;
1858 double full_sum = 0.;
1859
1860 for (int i = 0; i < hmat->n_rows; i++) {
1861
1862 const short int fi_sgn = cm->f_sgn[i];
1863 const double *hval_i = hmat->val + i*hmat->n_rows;
1864
1865 double *sval_i = sloc->val + i*sloc->n_rows;
1866 double row_sum = 0.;
1867 for (int j = 0; j < hmat->n_rows; j++) {
1868 const cs_real_t hsgn_coef = fi_sgn * cm->f_sgn[j] * hval_i[j];
1869 row_sum += hsgn_coef;
1870 sval_i[j] = hsgn_coef;
1871 }
1872
1873 sval_i[cm->n_fc] = -row_sum;
1874 sval_crow[i] = -row_sum;
1875 full_sum += row_sum;
1876
1877 }
1878
1879 /* (c, c) diagonal entry */
1880 sval_crow[cm->n_fc] = full_sum;
1881
1882 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
1883 _check_stiffness(cm->c_id, sloc);
1884 #endif
1885 }
1886
1887 /*----------------------------------------------------------------------------*/
1888 /*!
1889 * \brief Build a local stiffness matrix using the generic COST algo.
1890 * with the usage of bubble stabilization.
1891 * The computed matrix is stored in cb->loc and the related discrete
1892 * hodge operator in hodge->matrix
1893 * Case of CDO face-based schemes
1894 *
1895 * \param[in] cm pointer to a cs_cell_mesh_t structure
1896 * \param[in, out] hodge pointer to a cs_hodge_t structure
1897 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1898 */
1899 /*----------------------------------------------------------------------------*/
1900
1901 void
cs_hodge_fb_bubble_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)1902 cs_hodge_fb_bubble_get_stiffness(const cs_cell_mesh_t *cm,
1903 cs_hodge_t *hodge,
1904 cs_cell_builder_t *cb)
1905 {
1906 /* Sanity checks */
1907 assert(cs_eflag_test(cm->flag,
1908 CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1909
1910 /* Compute the local discrete Hodge operator */
1911 cs_hodge_edfp_bubble_get(cm, hodge, cb);
1912
1913 const cs_sdm_t *hmat = hodge->matrix;
1914
1915 /* Initialize the local stiffness matrix */
1916 cs_sdm_t *sloc = cb->loc;
1917 cs_sdm_square_init(cm->n_fc + 1, sloc);
1918
1919 double *sval_crow = sloc->val + cm->n_fc*sloc->n_rows;
1920 double full_sum = 0.;
1921
1922 for (int i = 0; i < hmat->n_rows; i++) {
1923
1924 const short int fi_sgn = cm->f_sgn[i];
1925 const double *hval_i = hmat->val + i*hmat->n_rows;
1926
1927 double *sval_i = sloc->val + i*sloc->n_rows;
1928 double row_sum = 0.;
1929 for (int j = 0; j < hmat->n_rows; j++) {
1930 const cs_real_t hsgn_coef = fi_sgn * cm->f_sgn[j] * hval_i[j];
1931 row_sum += hsgn_coef;
1932 sval_i[j] = hsgn_coef;
1933 }
1934
1935 sval_i[cm->n_fc] = -row_sum;
1936 sval_crow[i] = -row_sum;
1937 full_sum += row_sum;
1938
1939 }
1940
1941 /* (c, c) diagonal entry */
1942 sval_crow[cm->n_fc] = full_sum;
1943
1944 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
1945 _check_stiffness(cm->c_id, sloc);
1946 #endif
1947 }
1948
1949 /*----------------------------------------------------------------------------*/
1950 /*!
1951 * \brief Build a local stiffness matrix using the Voronoi algorithm
1952 * The computed matrix is stored in cb->loc and the related discrete
1953 * hodge operator in hodge->matrix
1954 * Case of CDO face-based schemes
1955 *
1956 * \param[in] cm pointer to a cs_cell_mesh_t structure
1957 * \param[in, out] hodge pointer to a cs_hodge_t structure
1958 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1959 */
1960 /*----------------------------------------------------------------------------*/
1961
1962 void
cs_hodge_fb_voro_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)1963 cs_hodge_fb_voro_get_stiffness(const cs_cell_mesh_t *cm,
1964 cs_hodge_t *hodge,
1965 cs_cell_builder_t *cb)
1966 {
1967 /* Sanity checks */
1968 assert(hodge->param->type == CS_HODGE_TYPE_EDFP);
1969 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
1970 assert(cs_eflag_test(cm->flag,
1971 CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1972
1973 /* Compute the local discrete Hodge operator */
1974 cs_hodge_edfp_voro_get(cm, hodge, cb);
1975
1976 const cs_sdm_t *hmat = hodge->matrix;
1977
1978 /* Initialize the local stiffness matrix */
1979 cs_sdm_t *sloc = cb->loc;
1980 cs_sdm_square_init(cm->n_fc + 1, sloc);
1981
1982 double full_sum = 0.;
1983 double *sval_crow = sloc->val + cm->n_fc*sloc->n_rows;
1984
1985 for (int i = 0; i < hmat->n_rows; i++) {
1986
1987 /* Hodge operator is diagonal */
1988 const double *hval_i = hmat->val + i*hmat->n_rows;
1989 const double row_sum = hval_i[i];
1990
1991 double *sval_i = sloc->val + i*sloc->n_rows;
1992
1993 sval_i[i] = hval_i[i];
1994 sval_i[cm->n_fc] = -row_sum;
1995 sval_crow[i] = -row_sum;
1996 full_sum += row_sum;
1997
1998 }
1999
2000 /* (c, c) diagonal entry */
2001 sval_crow[cm->n_fc] = full_sum;
2002
2003 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2004 _check_stiffness(cm->c_id, sloc);
2005 #endif
2006 }
2007
2008 /*----------------------------------------------------------------------------*/
2009 /*!
2010 * \brief Build a local stiffness matrix using the generic COST algo.
2011 * The computed matrix is stored in cb->loc and the related discrete
2012 * hodge operator in hodge->matrix
2013 * Case of CDO vertex-based schemes and an isotropic property
2014 *
2015 * \param[in] cm pointer to a cs_cell_mesh_t structure
2016 * \param[in, out] hodge pointer to a cs_hodge_t structure
2017 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2018 */
2019 /*----------------------------------------------------------------------------*/
2020
2021 void
cs_hodge_vb_cost_get_iso_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2022 cs_hodge_vb_cost_get_iso_stiffness(const cs_cell_mesh_t *cm,
2023 cs_hodge_t *hodge,
2024 cs_cell_builder_t *cb)
2025 {
2026 const cs_hodge_param_t *hodgep = hodge->param;
2027 const cs_property_data_t *ptyd = hodge->pty_data;
2028
2029 /* Sanity checks */
2030 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2031 assert(hodgep->algo == CS_HODGE_ALGO_COST);
2032 assert(cs_eflag_test(cm->flag,
2033 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2034 CS_FLAG_COMP_EV));
2035
2036 /* Set numbering and geometrical quantities Hodge builder */
2037 cs_real_3_t *pq = cb->vectors;
2038 cs_real_3_t *dq = cb->vectors + cm->n_ec;
2039
2040 _init_vb_geom_quant(cm, pq, dq);
2041
2042 /* Compute the upper right part of the local Hodge matrix
2043 * Rk: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
2044 * or DUAL->PRIMAL space */
2045 cs_sdm_square_init(cm->n_ec, hodge->matrix);
2046
2047 _compute_iso_hodge_ur(cm->n_ec,
2048 3*hodgep->coef*hodgep->coef,
2049 1./cm->vol_c,
2050 ptyd->value,
2051 (const cs_real_t (*)[3])pq,
2052 (const cs_real_t (*)[3])dq,
2053 cb, hodge->matrix);
2054
2055 _define_vb_stiffness(cm, hodge->matrix, cb->loc);
2056
2057 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2058 _check_stiffness(cm->c_id, cb->loc);
2059 #endif
2060 }
2061
2062 /*----------------------------------------------------------------------------*/
2063 /*!
2064 * \brief Build a local stiffness matrix using the generic COST algo.
2065 * The computed matrix is stored in cb->loc and the related discrete
2066 * hodge operator in hodge->matrix
2067 * Case of CDO vertex-based schemes and an anistropic property
2068 *
2069 * \param[in] cm pointer to a cs_cell_mesh_t structure
2070 * \param[in, out] hodge pointer to a cs_hodge_t structure
2071 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2072 */
2073 /*----------------------------------------------------------------------------*/
2074
2075 void
cs_hodge_vb_cost_get_aniso_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2076 cs_hodge_vb_cost_get_aniso_stiffness(const cs_cell_mesh_t *cm,
2077 cs_hodge_t *hodge,
2078 cs_cell_builder_t *cb)
2079 {
2080 const cs_hodge_param_t *hodgep = hodge->param;
2081 const cs_property_data_t *ptyd = hodge->pty_data;
2082
2083 /* Sanity checks */
2084 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2085 assert(hodgep->algo == CS_HODGE_ALGO_COST);
2086 assert(cs_eflag_test(cm->flag,
2087 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2088 CS_FLAG_COMP_EV));
2089
2090 /* Set numbering and geometrical quantities Hodge builder */
2091 cs_real_3_t *pq = cb->vectors;
2092 cs_real_3_t *dq = cb->vectors + cm->n_ec;
2093
2094 _init_vb_geom_quant(cm, pq, dq);
2095
2096 /* Compute the upper right part of the local Hodge matrix
2097 * Remark: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
2098 * or DUAL->PRIMAL space */
2099 cs_sdm_square_init(cm->n_ec, hodge->matrix);
2100
2101 _compute_aniso_hodge_ur(cm->n_ec,
2102 3*hodgep->coef*hodgep->coef,
2103 1./cm->vol_c,
2104 ptyd->tensor,
2105 (const cs_real_t (*)[3])pq,
2106 (const cs_real_t (*)[3])dq,
2107 cb, hodge->matrix);
2108
2109 _define_vb_stiffness(cm, hodge->matrix, cb->loc);
2110
2111 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2112 _check_stiffness(cm->c_id, cb->loc);
2113 #endif
2114 }
2115
2116 /*----------------------------------------------------------------------------*/
2117 /*!
2118 * \brief Build a local stiffness matrix using the generic Bubble algo.
2119 * The computed matrix is stored in cb->loc and the related discrete
2120 * hodge operator in hodge->matrix
2121 * Case of CDO vertex-based schemes and isotropic material property
2122 *
2123 * \param[in] cm pointer to a cs_cell_mesh_t structure
2124 * \param[in, out] hodge pointer to a cs_hodge_t structure
2125 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2126 */
2127 /*----------------------------------------------------------------------------*/
2128
2129 void
cs_hodge_vb_bubble_get_iso_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2130 cs_hodge_vb_bubble_get_iso_stiffness(const cs_cell_mesh_t *cm,
2131 cs_hodge_t *hodge,
2132 cs_cell_builder_t *cb)
2133 {
2134 const cs_hodge_param_t *hodgep = hodge->param;
2135 const cs_property_data_t *ptyd = hodge->pty_data;
2136
2137 /* Sanity checks */
2138 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2139 assert(hodgep->algo == CS_HODGE_ALGO_BUBBLE);
2140 assert(cs_eflag_test(cm->flag,
2141 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2142 CS_FLAG_COMP_EV));
2143
2144 /* Set numbering and geometrical quantities Hodge builder */
2145 cs_real_3_t *pq = cb->vectors;
2146 cs_real_3_t *dq = cb->vectors + cm->n_ec;
2147
2148 _init_vb_geom_quant(cm, pq, dq);
2149
2150 /* Compute the upper right part of the local Hodge matrix
2151 * Rk: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
2152 * or DUAL->PRIMAL space */
2153 cs_sdm_square_init(cm->n_ec, hodge->matrix);
2154
2155 _compute_iso_bubble_hodge_ur(cm->n_ec,
2156 hodgep->coef,
2157 1./cm->vol_c,
2158 ptyd->value,
2159 (const cs_real_t (*)[3])pq,
2160 (const cs_real_t (*)[3])dq,
2161 cb, hodge->matrix);
2162
2163 _define_vb_stiffness(cm, hodge->matrix, cb->loc);
2164
2165 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2166 _check_stiffness(cm->c_id, cb->loc);
2167 #endif
2168 }
2169
2170 /*----------------------------------------------------------------------------*/
2171 /*!
2172 * \brief Build a local stiffness matrix using the generic Bubble algo.
2173 * The computed matrix is stored in cb->loc and the related discrete
2174 * hodge operator in hodge->matrix
2175 * Case of CDO vertex-based schemes and anisotropic material property
2176 *
2177 * \param[in] cm pointer to a cs_cell_mesh_t structure
2178 * \param[in, out] hodge pointer to a cs_hodge_t structure
2179 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2180 */
2181 /*----------------------------------------------------------------------------*/
2182
2183 void
cs_hodge_vb_bubble_get_aniso_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2184 cs_hodge_vb_bubble_get_aniso_stiffness(const cs_cell_mesh_t *cm,
2185 cs_hodge_t *hodge,
2186 cs_cell_builder_t *cb)
2187 {
2188 const cs_hodge_param_t *hodgep = hodge->param;
2189 const cs_property_data_t *ptyd = hodge->pty_data;
2190
2191 /* Sanity checks */
2192 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2193 assert(hodgep->algo == CS_HODGE_ALGO_BUBBLE);
2194 assert(cs_eflag_test(cm->flag,
2195 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2196 CS_FLAG_COMP_EV));
2197
2198 /* Set numbering and geometrical quantities Hodge builder */
2199 cs_real_3_t *pq = cb->vectors;
2200 cs_real_3_t *dq = cb->vectors + cm->n_ec;
2201
2202 _init_vb_geom_quant(cm, pq, dq);
2203
2204 /* Compute the upper right part of the local Hodge matrix
2205 * Rk: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
2206 * or DUAL->PRIMAL space */
2207 cs_sdm_square_init(cm->n_ec, hodge->matrix);
2208
2209 _compute_aniso_bubble_hodge_ur(cm->n_ec,
2210 hodgep->coef,
2211 1./cm->vol_c,
2212 ptyd->tensor,
2213 (const cs_real_t (*)[3])pq,
2214 (const cs_real_t (*)[3])dq,
2215 cb, hodge->matrix);
2216
2217 _define_vb_stiffness(cm, hodge->matrix, cb->loc);
2218
2219 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2220 _check_stiffness(cm->c_id, cb->loc);
2221 #endif
2222 }
2223
2224 /*----------------------------------------------------------------------------*/
2225 /*!
2226 * \brief Build a local stiffness matrix using the Orthogonal
2227 * Consistent/Sub-Stabilization decomposition (OCS2) with a
2228 * subdivision of pvol_{e,c}.
2229 * The computed matrix is stored in cb->loc and the related discrete
2230 * hodge operator in hodge->matrix
2231 * Case Vb schemes and an anisotropic material property
2232 *
2233 * \param[in] cm pointer to a cs_cell_mesh_t structure
2234 * \param[in, out] hodge pointer to a cs_hodge_t structure
2235 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2236 */
2237 /*----------------------------------------------------------------------------*/
2238
2239 void
cs_hodge_vb_ocs2_get_aniso_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2240 cs_hodge_vb_ocs2_get_aniso_stiffness(const cs_cell_mesh_t *cm,
2241 cs_hodge_t *hodge,
2242 cs_cell_builder_t *cb)
2243 {
2244 const cs_hodge_param_t *hodgep = hodge->param;
2245 const cs_property_data_t *ptyd = hodge->pty_data;
2246
2247 /* Sanity checks */
2248 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2249 assert(hodgep->algo == CS_HODGE_ALGO_OCS2);
2250 assert(cs_eflag_test(cm->flag,
2251 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2252 CS_FLAG_COMP_EV | CS_FLAG_COMP_SEF));
2253
2254 /* Initialize the hodge matrix */
2255 cs_sdm_square_init(cm->n_ec, hodge->matrix);
2256
2257 /* Compute the upper right part of the local Hodge matrix
2258 * Rk: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
2259 * or DUAL->PRIMAL space */
2260 _compute_aniso_hepfd_ocs2_ur(3*hodgep->coef*hodgep->coef, ptyd->tensor, cm,
2261 cb, hodge->matrix);
2262
2263 _define_vb_stiffness(cm, hodge->matrix, cb->loc);
2264
2265 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2266 _check_stiffness(cm->c_id, cb->loc);
2267 #endif
2268 }
2269
2270 /*----------------------------------------------------------------------------*/
2271 /*!
2272 * \brief Build a local stiffness matrix using the generic COST algo.
2273 * The computed matrix is stored in cb->loc and the related discrete
2274 * hodge operator in hodge->matrix
2275 * Case of CDO vertex-based schemes
2276 *
2277 * \param[in] cm pointer to a cs_cell_mesh_t structure
2278 * \param[in, out] hodge pointer to a cs_hodge_t structure
2279 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2280 */
2281 /*----------------------------------------------------------------------------*/
2282
2283 void
cs_hodge_vb_cost_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2284 cs_hodge_vb_cost_get_stiffness(const cs_cell_mesh_t *cm,
2285 cs_hodge_t *hodge,
2286 cs_cell_builder_t *cb)
2287 {
2288 const cs_hodge_param_t *hodgep = hodge->param;
2289 const cs_property_data_t *ptyd = hodge->pty_data;
2290
2291 /* Sanity checks */
2292 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
2293 assert(hodgep->algo == CS_HODGE_ALGO_COST);
2294 assert(cs_eflag_test(cm->flag,
2295 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2296 CS_FLAG_COMP_EV));
2297
2298 /* Set numbering and geometrical quantities Hodge builder */
2299 cs_real_3_t *pq = cb->vectors;
2300 cs_real_3_t *dq = cb->vectors + cm->n_ec;
2301
2302 _init_vb_geom_quant(cm, pq, dq);
2303
2304 /* Compute additional geometrical quantities.
2305 Initial the local Hodge matrix with the consistency part which is
2306 constant over a cell.
2307 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
2308 and discrete Hodge operator from DUAL->PRIMAL space */
2309
2310 double *kappa = cb->values;
2311 double *alpha = cb->values + cm->n_ec;
2312
2313 /* Initialize the hodge matrix */
2314 cs_sdm_t *hmat = hodge->matrix;
2315 cs_sdm_square_init(cm->n_ec, hmat);
2316
2317 const double invcvol = 1/cm->vol_c;
2318 const double beta2 = hodgep->coef*hodgep->coef;
2319
2320 if (ptyd->is_iso || ptyd->is_unity) {
2321 _compute_cost_quant_iso(cm->n_ec, invcvol, ptyd->value,
2322 (const cs_real_t (*)[3])pq,
2323 (const cs_real_t (*)[3])dq,
2324 alpha, kappa, hmat);
2325
2326 }
2327 else
2328 _compute_cost_quant(cm->n_ec, invcvol, ptyd->tensor,
2329 (const cs_real_t (*)[3])pq,
2330 (const cs_real_t (*)[3])dq,
2331 alpha, kappa, hmat);
2332
2333 /* Initialize the local stiffness matrix */
2334 cs_sdm_t *sloc = cb->loc;
2335 cs_sdm_square_init(cm->n_vc, sloc);
2336
2337 for (int ei = 0; ei < cm->n_ec; ei++) { /* Loop on cell edges I */
2338
2339 const int shift_i = ei*cm->n_ec;
2340 const double *alpha_i = alpha + shift_i;
2341 const short int i1ei = cm->e2v_sgn[ei];
2342 const short int i1 = cm->e2v_ids[2*ei];
2343 const short int i2 = cm->e2v_ids[2*ei+1];
2344 const double *hi = hmat->val + shift_i;
2345
2346 assert(i1 < i2);
2347
2348 double *si1 = sloc->val + i1*sloc->n_rows;
2349 double *si2 = sloc->val + i2*sloc->n_rows;
2350
2351 /* Add contribution from the stabilization part for each sub-volume
2352 related to a primal entity */
2353 double stab_part = 0;
2354 for (int ek = 0; ek < cm->n_ec; ek++) /* Loop over sub-volumes */
2355 stab_part += kappa[ek] * alpha_i[ek] * alpha_i[ek];
2356
2357 /* Diagonal value: consistency part has already been computed */
2358 const double dval = hi[ei] + beta2 * stab_part;
2359
2360 si1[i1] += dval;
2361 si1[i2] -= dval;
2362 si2[i2] += dval;
2363
2364 /* Compute extra-diag entries */
2365 for (int ej = ei + 1; ej < cm->n_ec; ej++) { /* Loop on cell entities J */
2366
2367 const int shift_j = ej*cm->n_ec;
2368 const double *alpha_j = alpha + shift_j;
2369 const short int j1ej = cm->e2v_sgn[ej];
2370 const short int j1 = cm->e2v_ids[2*ej];
2371 const short int j2 = cm->e2v_ids[2*ej+1];
2372
2373 assert(j1 < j2);
2374
2375 double *sj1 = sloc->val + j1*sloc->n_rows;
2376 double *sj2 = sloc->val + j2*sloc->n_rows;
2377
2378 /* Add contribution from the stabilization part for each sub-volume
2379 related to a primal entity */
2380 stab_part = 0;
2381 for (int ek = 0; ek < cm->n_ec; ek++) /* Loop over sub-volumes */
2382 stab_part += kappa[ek] * alpha_i[ek] * alpha_j[ek];
2383
2384 /* Extra-diagonal value */
2385 const double xval = (hi[ej] + beta2 * stab_part) * i1ei * j1ej;
2386
2387 if (i2 < j1) { /* i1 < i2 < j1 < j2 */
2388 si1[j1] += xval;
2389 si1[j2] -= xval;
2390 si2[j1] -= xval;
2391 si2[j2] += xval;
2392 }
2393 else if (i2 == j1) { /* i1 < i2 = j1 < j2 */
2394 si1[j1] += xval;
2395 si1[j2] -= xval;
2396 si2[j1] -= 2*xval;
2397 si2[j2] += xval;
2398 }
2399 else if (i2 < j2) {
2400
2401 assert(i2 > j1);
2402 if (i1 < j1) /* i1 < j1 < i2 < j2 */
2403 si1[j1] += xval;
2404 else if (i1 == j1) /* i1 = j1 < i2 < j2 */
2405 si1[j1] += 2*xval;
2406 else /* j1 < i1 < i2 < j2 */
2407 sj1[i1] += xval;
2408
2409 si1[j2] -= xval;
2410 sj1[i2] -= xval;
2411 si2[j2] += xval;
2412
2413 }
2414 else if (i2 == j2) {
2415
2416 if (i1 < j1) /* i1 < j1 < i2 = j2 */
2417 si1[j1] += xval;
2418 else if (i1 == j1) /* i1 = j1 < i2 = j2 */
2419 si1[j1] += 2*xval;
2420 else /* j1 < i1 < i2 = j2 */
2421 sj1[i1] += xval;
2422
2423 si1[j2] -= xval;
2424 sj1[i2] -= xval;
2425 si2[j2] += 2*xval;
2426
2427 }
2428 else { /* i2 > j2 */
2429
2430 if (i1 < j1) { /* i1 < j1 < j2 < i2 */
2431 si1[j1] += xval;
2432 si1[j2] -= xval;
2433 }
2434 else if (i1 == j1) { /* i1 = j1 < j2 < i2 */
2435 si1[j1] += 2*xval;
2436 si1[j2] -= xval;
2437 }
2438 else if (i1 < j2) { /* j1 < i1 < j2 < i2 */
2439 sj1[i1] += xval;
2440 si1[j2] -= xval;
2441 }
2442 else if (i1 == j2) { /* j1 < i1 = j2 < i2 */
2443 sj1[i1] += xval;
2444 si1[j2] -= 2*xval;
2445 }
2446 else { /* j1 < j2 < i1 < i2 */
2447 sj1[i1] += xval;
2448 sj2[i1] -= xval;
2449 }
2450
2451 assert(i2 > j2);
2452 sj1[i2] -= xval;
2453 sj2[i2] += xval;
2454
2455 } /* End of tests */
2456
2457 } /* End of loop on J entities */
2458
2459 } /* End of loop on I entities */
2460
2461 /* Stiffness matrix is symmetric by construction */
2462 for (int ei = 0; ei < sloc->n_rows; ei++) {
2463 double *si = sloc->val + ei*sloc->n_rows;
2464 for (int ej = 0; ej < ei; ej++)
2465 si[ej] = sloc->val[ej*sloc->n_rows + ei];
2466 }
2467
2468 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2469 _check_stiffness(cm->c_id, sloc);
2470 #endif
2471 }
2472
2473 /*----------------------------------------------------------------------------*/
2474 /*!
2475 * \brief Build a local stiffness matrix using the Voronoi algorithm
2476 * The computed matrix is stored in cb->loc and the related discrete
2477 * hodge operator in hodge->matrix
2478 * Case of CDO vertex-based schemes
2479 *
2480 * \param[in] cm pointer to a cs_cell_mesh_t structure
2481 * \param[in, out] hodge pointer to a cs_hodge_t structure
2482 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2483 */
2484 /*----------------------------------------------------------------------------*/
2485
2486 void
cs_hodge_vb_voro_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2487 cs_hodge_vb_voro_get_stiffness(const cs_cell_mesh_t *cm,
2488 cs_hodge_t *hodge,
2489 cs_cell_builder_t *cb)
2490 {
2491 const cs_property_data_t *ptyd = hodge->pty_data;
2492
2493 /* Sanity checks */
2494 assert(hodge->param->type == CS_HODGE_TYPE_EPFD);
2495 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
2496 assert(cs_eflag_test(cm->flag,
2497 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
2498 CS_FLAG_COMP_EV));
2499
2500 /* Initialize the local stiffness matrix */
2501 cs_sdm_t *sloc = cb->loc;
2502 cs_sdm_square_init(cm->n_vc, sloc);
2503
2504 if (ptyd->is_iso || ptyd->is_unity) {
2505
2506 double dpty_val = 1.0; /* is_unity */
2507 if (ptyd->is_iso)
2508 dpty_val = ptyd->value;
2509
2510 /* Loop on cell edges */
2511 for (int ii = 0; ii < cm->n_ec; ii++) {
2512
2513 cs_nvec3_t dfq = cm->dface[ii];
2514 cs_quant_t peq = cm->edge[ii];
2515
2516 /* Only a diagonal term */
2517 const double dval = dpty_val * dfq.meas/peq.meas;
2518 const short int vi = cm->e2v_ids[2*ii];
2519 const short int vj = cm->e2v_ids[2*ii+1];
2520
2521 double *si = sloc->val + vi*sloc->n_rows;
2522 double *sj = sloc->val + vj*sloc->n_rows;
2523
2524 si[vi] += dval;
2525 sj[vj] += dval;
2526 si[vj] = sj[vi] = -dval; /* sgn_i * sgn_j = -1 */
2527
2528 } /* End of loop on cell edges */
2529
2530 }
2531 else { /* Diffusion property is anisotropic */
2532
2533 cs_real_3_t mv;
2534
2535 /* Loop on cell edges */
2536 for (int ii = 0; ii < cm->n_ec; ii++) {
2537
2538 cs_nvec3_t dfq = cm->dface[ii];
2539 cs_quant_t peq = cm->edge[ii];
2540
2541 cs_math_33_3_product((const cs_real_3_t *)ptyd->tensor, dfq.unitv, mv);
2542
2543 /* Only a diagonal term */
2544 const double dval = _dp3(mv, dfq.unitv) * dfq.meas/peq.meas;
2545 const short int vi = cm->e2v_ids[2*ii];
2546 const short int vj = cm->e2v_ids[2*ii+1];
2547
2548 double *si = sloc->val + vi*sloc->n_rows;
2549 double *sj = sloc->val + vj*sloc->n_rows;
2550
2551 si[vi] += dval;
2552 sj[vj] += dval;
2553 si[vj] = sj[vi] = -dval; /* sgn_j * sgn_i = -1 */
2554
2555 } /* End of loop on cell edges */
2556
2557 } /* Tensor-valued property */
2558
2559 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2560 _check_stiffness(cm->c_id, sloc);
2561 #endif
2562 }
2563
2564 /*----------------------------------------------------------------------------*/
2565 /*!
2566 * \brief Build a local stiffness matrix using the generic WBS algo.
2567 * WBS standing for Whitney Barycentric Subdivision (WBS)
2568 * algo.
2569 * The computed matrix is stored in cb->loc
2570 *
2571 * \param[in] cm pointer to a cs_cell_mesh_t structure
2572 * \param[in, out] hodge pointer to a cs_hodge_t structure
2573 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2574 */
2575 /*----------------------------------------------------------------------------*/
2576
2577 void
cs_hodge_vb_wbs_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2578 cs_hodge_vb_wbs_get_stiffness(const cs_cell_mesh_t *cm,
2579 cs_hodge_t *hodge,
2580 cs_cell_builder_t *cb)
2581 {
2582 const cs_property_data_t *ptyd = hodge->pty_data;
2583
2584 /* Sanity checks */
2585 assert(ptyd->need_tensor);
2586 assert(hodge->param->type == CS_HODGE_TYPE_EPFD);
2587 assert(hodge->param->algo == CS_HODGE_ALGO_WBS);
2588 assert(cs_eflag_test(cm->flag,
2589 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
2590 CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ | CS_FLAG_COMP_FEQ));
2591
2592 cs_real_3_t grd_c, grd_f, grd_v1, grd_v2, matg;
2593
2594 cs_real_3_t *uvc = cb->vectors;
2595 cs_real_3_t *glv = cb->vectors + cm->n_vc;
2596 cs_real_t *lvc = cb->values;
2597 cs_real_t *wvf = cb->values + cm->n_vc;
2598 cs_real_t *wef = cb->values + 2*cm->n_vc;
2599
2600 /* Initialize the local stiffness matrix */
2601 cs_sdm_t *sloc = cb->loc;
2602 cs_sdm_square_init(cm->n_vc, sloc);
2603
2604 /* Define the length and unit vector of the segment x_c --> x_v */
2605 for (short int v = 0; v < cm->n_vc; v++)
2606 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, lvc + v, uvc[v]);
2607
2608 /* Loop on cell faces */
2609 for (short int f = 0; f < cm->n_fc; f++) {
2610
2611 const cs_nvec3_t deq = cm->dedge[f];
2612
2613 /* Compute the gradient of the lagrange function related to a cell
2614 in each p_{f,c} and the weights for each vertex related to this face */
2615 cs_compute_grdfc_cw(f, cm, grd_c);
2616 cs_compute_wef_wvf(f, cm, wvf, wef);
2617
2618 /* Loop on face edges to scan p_{e,f,c} subvolumes */
2619 const short int *f2e_idx = cm->f2e_idx + f;
2620 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
2621 for (int i = 0; i < f2e_idx[1] - f2e_idx[0]; i++) {
2622
2623 const short int ee = 2*f2e_ids[i];
2624 const double subvol = wef[i]*cm->pvol_f[f];
2625 const short int v1 = cm->e2v_ids[ee];
2626 const short int v2 = cm->e2v_ids[ee+1];
2627
2628 /* Gradient of the lagrange function related to v1 and v2 */
2629 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])uvc, lvc,
2630 grd_v1, grd_v2);
2631
2632 /* Gradient of the lagrange function related to a face.
2633 This formula is a consequence of the Partition of the Unity */
2634 for (int k = 0; k < 3; k++)
2635 grd_f[k] = -(grd_c[k] + grd_v1[k] + grd_v2[k]);
2636
2637 /* Compute the gradient of the conforming reconstruction functions for
2638 each vertex of the cell in this subvol (pefc) */
2639 for (int si = 0; si < sloc->n_rows; si++) {
2640
2641 for (int k = 0; k < 3; k++)
2642 glv[si][k] = cm->wvc[si]*grd_c[k];
2643
2644 if (wvf[si] > 0) /* Face contrib. */
2645 for (int k = 0; k < 3; k++)
2646 glv[si][k] += wvf[si]*grd_f[k];
2647
2648 if (si == v1) /* Vertex 1 contrib */
2649 for (int k = 0; k < 3; k++)
2650 glv[si][k] += grd_v1[k];
2651
2652 if (si == v2) /* Vertex 2 contrib */
2653 for (int k = 0; k < 3; k++)
2654 glv[si][k] += grd_v2[k];
2655
2656 } /* Loop on cell vertices */
2657
2658 /* Build the upper right part */
2659 for (int si = 0; si < sloc->n_rows; si++) {
2660
2661 cs_math_33_3_product(ptyd->tensor, glv[si], matg);
2662
2663 /* Diagonal contribution */
2664 double *mi = sloc->val + si*sloc->n_rows;
2665 mi[si] += subvol * _dp3(matg, glv[si]);
2666
2667 /* Loop on vertices v_j (j > i) */
2668 for (int sj = si+1; sj < sloc->n_rows; sj++)
2669 mi[sj] += subvol * _dp3(matg, glv[sj]);
2670
2671 } /* Loop on vertices v_i */
2672
2673 }
2674
2675 } /* Loop on cell faces */
2676
2677 /* Matrix is symmetric by construction */
2678 cs_sdm_symm_ur(sloc);
2679
2680 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2681 _check_stiffness(cm->c_id, sloc);
2682 #endif
2683 }
2684
2685 /*----------------------------------------------------------------------------*/
2686 /*!
2687 * \brief Build a local stiffness matrix using the generic WBS algo.
2688 * WBS standing for Whitney Barycentric Subdivision (WBS) algo.
2689 * The computed matrix is stored in cb->loc
2690 *
2691 * \param[in] cm pointer to a cs_cell_mesh_t structure
2692 * \param[in, out] hodge pointer to a cs_hodge_t structure
2693 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2694 */
2695 /*----------------------------------------------------------------------------*/
2696
2697 void
cs_hodge_vcb_get_stiffness(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2698 cs_hodge_vcb_get_stiffness(const cs_cell_mesh_t *cm,
2699 cs_hodge_t *hodge,
2700 cs_cell_builder_t *cb)
2701 {
2702 const cs_property_data_t *ptyd = hodge->pty_data;
2703
2704 /* Sanity checks */
2705 assert(ptyd->need_tensor);
2706 assert(cs_eflag_test(cm->flag,
2707 CS_FLAG_COMP_PV | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
2708 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC));
2709
2710 cs_real_3_t grd_c, grd_f, grd_v1, grd_v2, matg, matg_c;
2711
2712 cs_real_3_t *uvc = cb->vectors;
2713 cs_real_3_t *glv = cb->vectors + cm->n_vc;
2714 cs_real_t *lvc = cb->values;
2715 cs_real_t *wvf = cb->values + cm->n_vc;
2716 cs_real_t *wef = cb->values + 2*cm->n_vc;
2717
2718
2719 /* Initialize the local stiffness matrix */
2720 const int nc_dofs = cm->n_vc + 1;
2721 /* index to the (cell,cell) entry */
2722 const int cc = nc_dofs*cm->n_vc + cm->n_vc;
2723
2724 cs_sdm_t *sloc = cb->loc;
2725 cs_sdm_square_init(nc_dofs, sloc);
2726
2727 /* Define the length and unit vector of the segment x_c --> x_v */
2728 for (short int v = 0; v < cm->n_vc; v++)
2729 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, lvc + v, uvc[v]);
2730
2731 /* Loop on cell faces */
2732 for (short int f = 0; f < cm->n_fc; f++) {
2733
2734 const cs_nvec3_t deq = cm->dedge[f];
2735
2736 /* Compute for the current face:
2737 - the gradient of the Lagrange function related xc in p_{f,c}
2738 - weights related to vertices
2739 - subvolume p_{ef,c} related to edges
2740 */
2741 cs_compute_grdfc_cw(f, cm, grd_c);
2742 cs_compute_wef_wvf(f, cm, wvf, wef);
2743
2744 /* Compute the contribution to the entry A(c,c) */
2745 cs_math_33_3_product(ptyd->tensor, grd_c, matg_c);
2746 sloc->val[cc] += cm->pvol_f[f] * _dp3(grd_c, matg_c);
2747
2748 /* Loop on face edges to scan p_{e,f,c} subvolumes */
2749 const short int *f2e_idx = cm->f2e_idx + f;
2750 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
2751 for (short int i = 0; i < f2e_idx[1] - f2e_idx[0]; i++) {
2752
2753 const short int ee = 2*f2e_ids[i];
2754 const double subvol = wef[i]*cm->pvol_f[f];
2755 const short int v1 = cm->e2v_ids[ee];
2756 const short int v2 = cm->e2v_ids[ee+1];
2757
2758 /* Gradient of the lagrange function related to v1 and v2 */
2759 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])uvc, lvc,
2760 grd_v1, grd_v2);
2761
2762 /* Gradient of the Lagrange function related to a face.
2763 This formula is a consequence of the Partition of the Unity */
2764 for (int k = 0; k < 3; k++)
2765 grd_f[k] = -(grd_c[k] + grd_v1[k] + grd_v2[k]);
2766
2767 /* Compute the gradient of the conforming reconstruction functions for
2768 each vertex of the cell in this subvol (pefc) */
2769 for (int si = 0; si < cm->n_vc; si++) {
2770
2771 for (int k = 0; k < 3; k++)
2772 glv[si][k] = 0;
2773
2774 if (wvf[si] > 0) /* Face contrib. */
2775 for (int k = 0; k < 3; k++)
2776 glv[si][k] += wvf[si]*grd_f[k];
2777
2778 if (si == v1) /* Vertex 1 contrib */
2779 for (int k = 0; k < 3; k++)
2780 glv[si][k] += grd_v1[k];
2781
2782 if (si == v2) /* Vertex 2 contrib */
2783 for (int k = 0; k < 3; k++)
2784 glv[si][k] += grd_v2[k];
2785
2786 } /* Loop on cell vertices */
2787
2788 /* Build the upper right part (v-v and v-c)
2789 Be careful: sloc->n_rows = cm->n_vc + 1 */
2790 for (int si = 0; si < cm->n_vc; si++) {
2791
2792 double *mi = sloc->val + si*sloc->n_rows;
2793
2794 /* Add v-c contribution */
2795 mi[cm->n_vc] += subvol * _dp3(matg_c, glv[si]);
2796
2797 /* Add v-v contribution */
2798 cs_math_33_3_product(ptyd->tensor, glv[si], matg);
2799
2800 /* Loop on vertices v_j (j >= i) */
2801 for (int sj = si; sj < cm->n_vc; sj++)
2802 mi[sj] += subvol * _dp3(matg, glv[sj]);
2803
2804 } /* Loop on vertices v_i */
2805
2806 }
2807
2808 } /* Loop on cell faces */
2809
2810 /* Matrix is symmetric by construction */
2811 cs_sdm_symm_ur(sloc);
2812
2813 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
2814 _check_stiffness(cm->c_id, sloc);
2815 #endif
2816 }
2817
2818 /*----------------------------------------------------------------------------*/
2819 /*!
2820 * \brief Build a local Hodge operator on a given cell which is equivalent of
2821 * a mass matrix. It relies on a CO+ST algo. and is specific to CDO-Fb
2822 * schemes.
2823 * The discrete Hodge operator is stored in hodge->matrix
2824 *
2825 * \param[in] cm pointer to a cs_cell_mesh_t structure
2826 * \param[in, out] hodge pointer to a cs_hodge_t structure
2827 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2828 */
2829 /*----------------------------------------------------------------------------*/
2830
2831 void
cs_hodge_fb_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2832 cs_hodge_fb_get(const cs_cell_mesh_t *cm,
2833 cs_hodge_t *hodge,
2834 cs_cell_builder_t *cb)
2835 {
2836 CS_UNUSED(cb);
2837
2838 /* Sanity check */
2839 assert(hodge != NULL);
2840 assert(hodge->param->type == CS_HODGE_TYPE_FB);
2841 assert(hodge->param->algo == CS_HODGE_ALGO_COST);
2842 assert(hodge->pty_data->is_unity);
2843 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ));
2844
2845 const int n_cols = cm->n_fc + 1;
2846 const cs_real_t over_cell = 1./(cm->vol_c*cm->vol_c);
2847
2848 /* Initialize the local matrix related to this discrete Hodge operator */
2849 cs_sdm_t *hmat = hodge->matrix;
2850 cs_sdm_square_init(n_cols, hmat);
2851
2852 /* cell-cell entry (cell-face and face-cell block are null) */
2853 hmat->val[cm->n_fc*(n_cols + 1)] = cm->vol_c;
2854
2855 /* Compute the inertia (useful for the face-face block) */
2856 cs_real_t inertia_tensor[3][3];
2857 cs_compute_inertia_tensor(cm, cm->xc, inertia_tensor);
2858
2859 for (short int fi = 0; fi < cm->n_fc; fi++) {
2860
2861 const cs_quant_t pfq_i = cm->face[fi];
2862 const short int sgn_i = cm->f_sgn[fi];
2863
2864 /* Diagonal entry (a bit more optimized) */
2865 cs_real_t dval = 0;
2866 for (int k = 0; k < 3; k++) {
2867 dval += pfq_i.unitv[k]*pfq_i.unitv[k]*inertia_tensor[k][k];
2868 for (int l = k+1; l < 3; l++)
2869 dval += 2*pfq_i.unitv[k]*pfq_i.unitv[l]*inertia_tensor[k][l];
2870 }
2871
2872 hmat->val[fi*(n_cols+1)] = over_cell * pfq_i.meas*pfq_i.meas * dval;
2873
2874 /* Extra-diag entry */
2875 for (short int fj = fi + 1; fj < cm->n_fc; fj++) {
2876
2877 const cs_quant_t pfq_j = cm->face[fj];
2878
2879 cs_real_t eval = 0;
2880 for (int k = 0; k < 3; k++) {
2881 for (int l = 0; l < 3; l++)
2882 eval += pfq_i.unitv[k]*pfq_j.unitv[l]*inertia_tensor[k][l];
2883 }
2884 eval *= over_cell * sgn_i*pfq_i.meas * cm->f_sgn[fj]*pfq_j.meas;
2885
2886 /* Symmetric by construction */
2887 hmat->val[fi*n_cols+fj] = eval;
2888 hmat->val[fj*n_cols+fi] = eval;
2889
2890 } /* Loop on cell faces (fj) */
2891
2892 } /* Loop on cell faces (fi) */
2893 }
2894
2895 /*----------------------------------------------------------------------------*/
2896 /*!
2897 * \brief Build a local Hodge operator for a given cell using the Voronoi
2898 * algo. This leads to a diagonal operator.
2899 * This function is specific for vertex+cell-based schemes
2900 * The discrete Hodge operator is stored in hodge->matrix
2901 *
2902 * \param[in] cm pointer to a cs_cell_mesh_t structure
2903 * \param[in, out] hodge pointer to a cs_hodge_t structure
2904 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2905 */
2906 /*----------------------------------------------------------------------------*/
2907
2908 void
cs_hodge_vcb_voro_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2909 cs_hodge_vcb_voro_get(const cs_cell_mesh_t *cm,
2910 cs_hodge_t *hodge,
2911 cs_cell_builder_t *cb)
2912 {
2913 CS_UNUSED(cb);
2914
2915 const cs_property_data_t *ptyd = hodge->pty_data;
2916
2917 /* Sanity check */
2918 assert(cm != NULL && hodge != NULL);
2919 assert(hodge->param->type == CS_HODGE_TYPE_VC);
2920 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
2921 assert(ptyd->is_iso);
2922 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
2923
2924 cs_sdm_t *hmat = hodge->matrix;
2925
2926 /* Initialize the local matrix related to this discrete Hodge operator */
2927 cs_sdm_square_init(cm->n_vc + 1, hmat);
2928
2929 const int msize = cm->n_vc + 1;
2930
2931 if (ptyd->is_unity) {
2932
2933 /* H(c,c) = 0.25*|c| */
2934 hmat->val[msize*cm->n_vc] = 0.25*cm->vol_c;
2935
2936 /* H(c,c) = 0.75*|dcell(v) \cap c| */
2937 const double vol_coef = 0.75*cm->vol_c;
2938 for (short int vi = 0; vi < cm->n_vc; vi++)
2939 hmat->val[msize*vi] = vol_coef*cm->wvc[vi];
2940
2941 }
2942 else {
2943
2944 /* H(c,c) = 0.25*|c| */
2945 hmat->val[msize*cm->n_vc] = ptyd->value*0.25*cm->vol_c;
2946
2947 /* H(c,c) = 0.75*|dcell(v) \cap c| */
2948 const double vol_coef = 0.75*cm->vol_c*ptyd->value;
2949 for (short int vi = 0; vi < cm->n_vc; vi++)
2950 hmat->val[msize*vi] = vol_coef*cm->wvc[vi];
2951
2952 }
2953
2954 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
2955 if (cm->c_id % CS_HODGE_MODULO == 0) {
2956 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
2957 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
2958 }
2959 #endif
2960 }
2961
2962 /*----------------------------------------------------------------------------*/
2963 /*!
2964 * \brief Build a local Hodge operator for a given cell using the WBS algo.
2965 * This function is specific for vertex+cell-based schemes
2966 * The discrete Hodge operator is stored in hodge->matrix
2967 *
2968 * \param[in] cm pointer to a cs_cell_mesh_t structure
2969 * \param[in, out] hodge pointer to a cs_hodge_t structure
2970 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2971 */
2972 /*----------------------------------------------------------------------------*/
2973
2974 void
cs_hodge_vcb_wbs_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)2975 cs_hodge_vcb_wbs_get(const cs_cell_mesh_t *cm,
2976 cs_hodge_t *hodge,
2977 cs_cell_builder_t *cb)
2978 {
2979 const cs_property_data_t *ptyd = hodge->pty_data;
2980
2981 /* Sanity check */
2982 assert(cb != NULL && hodge != NULL);
2983 assert(hodge->param->type == CS_HODGE_TYPE_VC);
2984 assert(hodge->param->algo == CS_HODGE_ALGO_WBS);
2985 assert(ptyd->is_iso);
2986 assert(cs_eflag_test(cm->flag,
2987 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
2988 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_HFQ));
2989
2990 cs_sdm_t *hmat = hodge->matrix;
2991
2992 /* Initialize the local matrix related to this discrete Hodge operator */
2993 cs_sdm_square_init(cm->n_vc + 1, hmat);
2994
2995 double *wvf = cb->values;
2996 double *wef = cb->values + cm->n_vc;
2997
2998 const int msize = cm->n_vc + 1;
2999 const double c_coef1 = 0.2*cm->vol_c;
3000 const double c_coef2 = cs_hodge_vc_coef * cm->vol_c;
3001
3002 /* H(c,c) = 0.1*|c| */
3003 hmat->val[msize*cm->n_vc + cm->n_vc] = 0.1*cm->vol_c;
3004
3005 /* Initialize the upper part of the local Hodge matrix
3006 diagonal and cell column entries */
3007 for (short int vi = 0; vi < cm->n_vc; vi++) {
3008
3009 double *mi = hmat->val + vi*msize;
3010
3011 mi[vi] = c_coef1 * cm->wvc[vi]; /* Diagonal entry */
3012 for (short int vj = vi+1; vj < cm->n_vc; vj++)
3013 mi[vj] = 0.;
3014 mi[cm->n_vc] = c_coef2 * cm->wvc[vi]; /* Cell column */
3015
3016 } /* Loop on cell vertices */
3017
3018 /* Loop on each pef and add the contribution */
3019 for (short int f = 0; f < cm->n_fc; f++) {
3020
3021 /* Define useful quantities for WBS algo. */
3022 cs_compute_wef_wvf(f, cm, wvf, wef);
3023
3024 const double f_coef = 0.3 * cm->pvol_f[f];
3025 const double ef_coef = 0.05 * cm->pvol_f[f];
3026
3027 /* Add face contribution:
3028 Diagonal entry H(i,i) += 0.3*wif*wif*pfc_vol
3029 Extra-diag. entry H(i,j) += 0.3*wjf*wif*pfc_vol */
3030 for (short int vi = 0; vi < cm->n_vc; vi++) {
3031
3032 const double coef_if = f_coef * wvf[vi];
3033 double *mi = hmat->val + vi*msize;
3034
3035 /* Diagonal and Extra-diagonal entries: Add face contribution */
3036 for (short int vj = vi; vj < cm->n_vc; vj++)
3037 mi[vj] += coef_if * wvf[vj];
3038
3039 } /* Extra-diag entries */
3040
3041 /* Add edge-face contribution (only extra-diag) = 0.05 * |p_{ef,c}| */
3042 const short int *f2e_idx = cm->f2e_idx + f;
3043 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
3044 for (short int i = 0; i < f2e_idx[1] - f2e_idx[0]; i++) {
3045
3046 const short int ee = 2*f2e_ids[i];
3047 const short int v1 = cm->e2v_ids[ee];
3048 const short int v2 = cm->e2v_ids[ee+1];
3049
3050 /* Sanity check */
3051 assert(v1 > -1 && v2 > -1);
3052 if (v1 < v2)
3053 hmat->val[v1*msize+v2] += ef_coef * wef[i];
3054 else
3055 hmat->val[v2*msize+v1] += ef_coef * wef[i];
3056
3057 } /* Loop on face edges */
3058
3059 } /* Loop on cell faces */
3060
3061 /* Take into account the value of the associated property */
3062 if (!ptyd->is_unity) {
3063 for (short int vi = 0; vi < msize; vi++) {
3064 double *mi = hmat->val + vi*msize;
3065 for (short int vj = vi; vj < msize; vj++)
3066 mi[vj] *= ptyd->value;
3067 }
3068 }
3069
3070 /* Local matrix is symmetric by construction. Set the lower part. */
3071 for (short int vj = 0; vj < msize; vj++) {
3072 double *mj = hmat->val + vj*msize;
3073 for (short int vi = vj+1; vi < msize; vi++)
3074 hmat->val[vi*msize + vj] = mj[vi];
3075 }
3076
3077 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3078 if (cm->c_id % CS_HODGE_MODULO == 0) {
3079 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3080 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
3081 }
3082 #endif
3083 }
3084
3085 /*----------------------------------------------------------------------------*/
3086 /*!
3087 * \brief Build a local Hodge operator for a given cell using WBS algo.
3088 * Hodge op. from primal vertices to dual cells.
3089 * This function is specific for vertex-based schemes
3090 * The discrete Hodge operator is stored in hodge->matrix
3091 *
3092 * \param[in] cm pointer to a cs_cell_mesh_t structure
3093 * \param[in, out] hodge pointer to a cs_hodge_t structure
3094 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3095 */
3096 /*----------------------------------------------------------------------------*/
3097
3098 void
cs_hodge_vpcd_wbs_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3099 cs_hodge_vpcd_wbs_get(const cs_cell_mesh_t *cm,
3100 cs_hodge_t *hodge,
3101 cs_cell_builder_t *cb)
3102 {
3103 const cs_property_data_t *ptyd = hodge->pty_data;
3104
3105 /* Sanity checks */
3106 assert(cb != NULL);
3107 assert(hodge->param->type == CS_HODGE_TYPE_VPCD);
3108 assert(hodge->param->algo == CS_HODGE_ALGO_WBS);
3109 assert(ptyd->is_iso);
3110 assert(cs_eflag_test(cm->flag,
3111 CS_FLAG_COMP_PVQ |CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
3112 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFC));
3113
3114 double *wvf = cb->values;
3115 double *wef = cb->values + cm->n_vc;
3116
3117 cs_sdm_t *hmat = hodge->matrix;
3118 cs_sdm_square_init(cm->n_vc, hmat);
3119
3120 const double c_coef = 0.1*cm->vol_c;
3121
3122 /* Initialize the upper part of the local Hodge matrix */
3123 for (short int vi = 0; vi < cm->n_vc; vi++) {
3124
3125 double *mi = hmat->val + vi*cm->n_vc;
3126 const double vi_coef = 4 * c_coef * cm->wvc[vi];
3127
3128 /* Diag. entry has an additional contrib */
3129 mi[vi] = vi_coef * (0.5 + cm->wvc[vi]);
3130 for (short int vj = vi + 1; vj < cm->n_vc; vj++)
3131 mi[vj] = vi_coef * cm->wvc[vj]; /* Extra-diagonal entries */
3132
3133 } /* Loop on cell vertices */
3134
3135 /* Loop on each pef and add the contribution */
3136 for (short int f = 0; f < cm->n_fc; f++) {
3137
3138 /* Define useful quantities for WBS algo. */
3139 cs_compute_wef_wvf(f, cm, wvf, wef);
3140
3141 const double f_coef = 0.3 * cm->pvol_f[f];
3142 const double ef_coef = 0.05 * cm->pvol_f[f];
3143
3144 /* Add face contribution:
3145 Diagonal entry H(i,i) += 0.3*wif*wif*pfc_vol
3146 Extra-diag. entry H(i,j) += 0.3*wjf*wif*pfc_vol */
3147 for (short int vi = 0; vi < cm->n_vc; vi++) {
3148
3149 double *mi = hmat->val + vi*cm->n_vc;
3150
3151 /* Diagonal and Extra-diagonal entries: Add face contribution */
3152 const double coef_if = f_coef * wvf[vi];
3153 for (short int vj = vi; vj < cm->n_vc; vj++)
3154 mi[vj] += coef_if * wvf[vj];
3155
3156 } /* Face contribution */
3157
3158 /* Add edge-face contribution (only extra-diag) = 0.05 * |p_{ef,c}| */
3159 const short int *f2e_idx = cm->f2e_idx + f;
3160 const short int *f2e_ids = cm->f2e_ids + f2e_idx[0];
3161 for (short int i = 0; i < f2e_idx[1] - f2e_idx[0]; i++) {
3162
3163 const short int ee = 2*f2e_ids[i];
3164 const short int v1 = cm->e2v_ids[ee];
3165 const short int v2 = cm->e2v_ids[ee+1];
3166
3167 /* Sanity check */
3168 assert(v1 > -1 && v2 > -1);
3169 if (v1 < v2)
3170 hmat->val[v1*cm->n_vc+v2] += ef_coef * wef[i];
3171 else
3172 hmat->val[v2*cm->n_vc+v1] += ef_coef * wef[i];
3173
3174 } /* Loop on face edges */
3175
3176 } /* Loop on cell faces */
3177
3178 /* Take into account the value of the associated property */
3179 if (!ptyd->is_unity) {
3180 for (short int vi = 0; vi < cm->n_vc; vi++) {
3181 double *mi = hmat->val + vi*cm->n_vc;
3182 for (short int vj = vi; vj < cm->n_vc; vj++)
3183 mi[vj] *= ptyd->value;
3184 }
3185 }
3186
3187 /* Local matrix is symmetric by construction. Set the lower part. */
3188 for (short int vj = 0; vj < cm->n_vc; vj++) {
3189 double *mj = hmat->val + vj*cm->n_vc;
3190 for (short int vi = vj+1; vi < cm->n_vc; vi++)
3191 hmat->val[vi*cm->n_vc + vj] = mj[vi];
3192 }
3193
3194 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3195 if (cm->c_id % CS_HODGE_MODULO == 0) {
3196 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3197 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
3198 }
3199 #endif
3200 }
3201
3202 /*----------------------------------------------------------------------------*/
3203 /*!
3204 * \brief Build a local Hodge operator for a given cell using VORONOI algo.
3205 * Hodge op. from primal vertices to dual cells.
3206 * This function is specific for vertex-based schemes
3207 * The discrete Hodge operator is stored in hodge->matrix
3208 *
3209 * \param[in] cm pointer to a cs_cell_mesh_t structure
3210 * \param[in, out] hodge pointer to a cs_hodge_t structure
3211 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3212 */
3213 /*----------------------------------------------------------------------------*/
3214
3215 void
cs_hodge_vpcd_voro_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3216 cs_hodge_vpcd_voro_get(const cs_cell_mesh_t *cm,
3217 cs_hodge_t *hodge,
3218 cs_cell_builder_t *cb)
3219 {
3220 CS_UNUSED(cb);
3221
3222 const cs_property_data_t *ptyd = hodge->pty_data;
3223
3224 /* Sanity check */
3225 assert(hodge->param->type == CS_HODGE_TYPE_VPCD);
3226 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
3227 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
3228
3229 cs_sdm_t *hmat = hodge->matrix;
3230 cs_sdm_square_init(cm->n_vc, hmat);
3231
3232 const int stride = cm->n_vc + 1;
3233 if (ptyd->is_unity) {
3234
3235 for (int v = 0; v < cm->n_vc; v++)
3236 hmat->val[v*stride] = cm->wvc[v] * cm->vol_c;
3237
3238 }
3239 else {
3240
3241 const cs_real_t coef = ptyd->value * cm->vol_c;
3242 for (int v = 0; v < cm->n_vc; v++)
3243 hmat->val[v*stride] = coef * cm->wvc[v];
3244
3245 }
3246
3247 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3248 if (cm->c_id % CS_HODGE_MODULO == 0) {
3249 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3250 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
3251 }
3252 #endif
3253 }
3254
3255 /*----------------------------------------------------------------------------*/
3256 /*!
3257 * \brief Build a local Hodge operator for a given cell using VORONOI algo.
3258 * Hodge op. from primal edges to dual faces.
3259 * The discrete Hodge operator is stored in hodge->matrix
3260 * This function is specific for vertex-based schemes
3261 *
3262 * \param[in] cm pointer to a cs_cell_mesh_t structure
3263 * \param[in, out] hodge pointer to a cs_hodge_t structure
3264 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3265 */
3266 /*----------------------------------------------------------------------------*/
3267
3268 void
cs_hodge_epfd_voro_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3269 cs_hodge_epfd_voro_get(const cs_cell_mesh_t *cm,
3270 cs_hodge_t *hodge,
3271 cs_cell_builder_t *cb)
3272 {
3273 CS_UNUSED(cb);
3274
3275 const cs_property_data_t *ptyd = hodge->pty_data;
3276
3277 /* Sanity check */
3278 assert(hodge->param->type == CS_HODGE_TYPE_EPFD);
3279 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
3280 assert(cs_eflag_test(cm->flag,
3281 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ | CS_FLAG_COMP_SEF));
3282
3283 /* Initialize the local matrix related to this discrete Hodge operator */
3284 cs_sdm_t *hmat = hodge->matrix;
3285 cs_sdm_square_init(cm->n_ec, hmat);
3286
3287 /* Set the diagonal entries */
3288 if (ptyd->is_iso) {
3289
3290 for (short int e = 0; e < cm->n_ec; e++)
3291 hmat->val[e*cm->n_ec+e] =
3292 ptyd->value * cm->dface[e].meas/cm->edge[e].meas;
3293
3294 }
3295 else {
3296
3297 cs_real_3_t mv;
3298 for (short int f = 0; f < cm->n_fc; f++) {
3299 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
3300 const short int e = cm->f2e_ids[i];
3301 const cs_nvec3_t *sefc = cm->sefc + i;
3302 cs_math_33_3_product(ptyd->tensor, sefc->unitv, mv);
3303 hmat->val[e*cm->n_ec+e] += sefc->meas * _dp3(mv, sefc->unitv);
3304 }
3305 }
3306
3307 for (short int e = 0; e < cm->n_ec; e++)
3308 hmat->val[e*cm->n_ec+e] /= cm->edge[e].meas;
3309
3310 } /* anisotropic */
3311
3312 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3313 if (cm->c_id % CS_HODGE_MODULO == 0) {
3314 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3315 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
3316 }
3317 #endif
3318 }
3319
3320 /*----------------------------------------------------------------------------*/
3321 /*!
3322 * \brief Build a local Hodge operator for a given cell using the COST algo.
3323 * Hodge op. from primal edges to dual faces.
3324 * The discrete Hodge operator is stored in hodge->matrix
3325 * This function is specific for vertex-based schemes
3326 *
3327 * \param[in] cm pointer to a cs_cell_mesh_t structure
3328 * \param[in, out] hodge pointer to a cs_hodge_t structure
3329 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3330 */
3331 /*----------------------------------------------------------------------------*/
3332
3333 void
cs_hodge_epfd_cost_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3334 cs_hodge_epfd_cost_get(const cs_cell_mesh_t *cm,
3335 cs_hodge_t *hodge,
3336 cs_cell_builder_t *cb)
3337 {
3338 const cs_hodge_param_t *hodgep = hodge->param;
3339 const cs_property_data_t *ptyd = hodge->pty_data;
3340
3341 /* Sanity check */
3342 assert(cb != NULL && hodgep != NULL);
3343 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
3344 assert(hodgep->algo == CS_HODGE_ALGO_COST);
3345 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ));
3346
3347 /* Initialize the local matrix related to this discrete Hodge operator */
3348 cs_sdm_t *hmat = hodge->matrix;
3349 cs_sdm_square_init(cm->n_ec, hmat);
3350
3351 /* Set numbering and geometrical quantities Hodge builder */
3352 cs_real_3_t *pq = cb->vectors;
3353 cs_real_3_t *dq = cb->vectors + cm->n_ec;
3354
3355 _init_vb_geom_quant(cm, pq, dq);
3356
3357 /* Compute additional geometrical quantities: qmq and T
3358 Initial the local Hodge matrix with the consistency part which is
3359 constant over a cell.
3360 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3361 and discrete Hodge operator from DUAL->PRIMAL space */
3362
3363 double *kappa = cb->values;
3364 double *alpha = cb->values + cm->n_ec;
3365
3366 if (ptyd->is_unity)
3367 _compute_cost_quant_iso(cm->n_ec, 1/cm->vol_c, 1.0,
3368 (const cs_real_t (*)[3])pq,
3369 (const cs_real_t (*)[3])dq,
3370 alpha, kappa, hmat);
3371 else if (ptyd->is_iso)
3372 _compute_cost_quant_iso(cm->n_ec, 1/cm->vol_c, ptyd->value,
3373 (const cs_real_t (*)[3])pq,
3374 (const cs_real_t (*)[3])dq,
3375 alpha, kappa, hmat);
3376 else
3377 _compute_cost_quant(cm->n_ec, 1/cm->vol_c,
3378 (const cs_real_3_t *)ptyd->tensor,
3379 (const cs_real_t (*)[3])pq,
3380 (const cs_real_t (*)[3])dq,
3381 alpha, kappa, hmat);
3382
3383 double beta2 = hodgep->coef*hodgep->coef;
3384 _compute_hodge_cost(cm->n_ec, beta2, alpha, kappa, hmat->val);
3385
3386 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3387 _check_vector_hodge(cm->c_id,
3388 (const cs_real_t (*)[3])pq, (const cs_real_t (*)[3])dq,
3389 hodge, cb);
3390 #endif
3391 }
3392
3393 /*----------------------------------------------------------------------------*/
3394 /*!
3395 * \brief Build a local Hodge operator for a given cell using the COST algo.
3396 * with a bubble stabilization.
3397 * The discrete Hodge operator is stored in hodge->matrix
3398 * Hodge op. from primal edges to dual faces. This function is
3399 * specific for vertex-based schemes
3400 *
3401 * \param[in] cm pointer to a cs_cell_mesh_t structure
3402 * \param[in, out] hodge pointer to a cs_hodge_t structure
3403 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3404 */
3405 /*----------------------------------------------------------------------------*/
3406
3407 void
cs_hodge_epfd_bubble_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3408 cs_hodge_epfd_bubble_get(const cs_cell_mesh_t *cm,
3409 cs_hodge_t *hodge,
3410 cs_cell_builder_t *cb)
3411 {
3412 const cs_hodge_param_t *hodgep = hodge->param;
3413 const cs_property_data_t *ptyd = hodge->pty_data;
3414
3415 /* Sanity check */
3416 assert(cb != NULL && hodgep != NULL && ptyd != NULL);
3417 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
3418 assert(hodgep->algo == CS_HODGE_ALGO_BUBBLE);
3419 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ));
3420
3421 /* Set numbering and geometrical quantities Hodge builder */
3422 cs_real_3_t *pq = cb->vectors;
3423 cs_real_3_t *dq = cb->vectors + cm->n_ec;
3424
3425 _init_vb_geom_quant(cm, pq, dq);
3426
3427 /* Compute the upper right part of the local Hodge matrix
3428 * Rk: Switch arguments between discrete Hodge operator from PRIMAL->DUAL
3429 * or DUAL->PRIMAL space */
3430 cs_sdm_t *hmat = hodge->matrix;
3431 cs_sdm_square_init(cm->n_ec, hmat);
3432
3433 if (ptyd->is_iso || ptyd->is_unity)
3434 _compute_iso_bubble_hodge_ur(cm->n_ec,
3435 hodgep->coef,
3436 1./cm->vol_c,
3437 ptyd->value,
3438 (const cs_real_t (*)[3])pq,
3439 (const cs_real_t (*)[3])dq,
3440 cb, hmat);
3441 else
3442 _compute_aniso_bubble_hodge_ur(cm->n_ec,
3443 hodgep->coef,
3444 1./cm->vol_c,
3445 ptyd->tensor,
3446 (const cs_real_t (*)[3])pq,
3447 (const cs_real_t (*)[3])dq,
3448 cb, hmat);
3449
3450 /* Hodge matrix is symmetric by construction */
3451 cs_sdm_symm_ur(hmat);
3452
3453 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3454 _check_vector_hodge(cm->c_id,
3455 (const cs_real_t (*)[3])pq, (const cs_real_t (*)[3])dq,
3456 hodge, cb);
3457 #endif
3458 }
3459
3460 /*----------------------------------------------------------------------------*/
3461 /*!
3462 * \brief Build a local Hodge operator for a given cell using the Orthogonal
3463 * Consistent/Sub-Stabilization decomposition (OCS2) with a
3464 * subdivision of pvol_{e,c}.
3465 * The discrete Hodge operator is stored in hodge->matrix
3466 * Hodge op. from primal edges to dual faces.
3467 * This function is specific for vertex-based schemes
3468 *
3469 * \param[in] cm pointer to a cs_cell_mesh_t structure
3470 * \param[in, out] hodge pointer to a cs_hodge_t structure
3471 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3472 */
3473 /*----------------------------------------------------------------------------*/
3474
3475 void
cs_hodge_epfd_ocs2_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3476 cs_hodge_epfd_ocs2_get(const cs_cell_mesh_t *cm,
3477 cs_hodge_t *hodge,
3478 cs_cell_builder_t *cb)
3479 {
3480 const cs_hodge_param_t *hodgep = hodge->param;
3481 const cs_property_data_t *ptyd = hodge->pty_data;
3482
3483 /* Sanity check */
3484 assert(cb != NULL && hodgep != NULL);
3485 assert(hodgep->type == CS_HODGE_TYPE_EPFD);
3486 assert(hodgep->algo == CS_HODGE_ALGO_OCS2);
3487 assert(cs_eflag_test(cm->flag,
3488 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ |
3489 CS_FLAG_COMP_EV | CS_FLAG_COMP_SEF));
3490
3491 /* Initialize the local matrix related to this discrete Hodge operator */
3492 cs_sdm_t *hmat = hodge->matrix;
3493 cs_sdm_square_init(cm->n_ec, hmat);
3494
3495 /* Compute the upper right part of the local Hodge matrix */
3496 _compute_aniso_hepfd_ocs2_ur(3*hodgep->coef*hodgep->coef, ptyd->tensor,
3497 cm, cb, hmat);
3498
3499 /* Hodge operator leads to a symmetric matrix by construction */
3500 cs_sdm_symm_ur(hmat);
3501
3502 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3503 if (cm->c_id % CS_HODGE_MODULO == 0) {
3504 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3505 cs_sdm_dump(cm->c_id, NULL, NULL, hmat);
3506 }
3507 #endif
3508 }
3509
3510 /*----------------------------------------------------------------------------*/
3511 /*!
3512 * \brief Build a local Hodge operator for a given cell using VORONOI algo.
3513 * Hodge op. from primal faces to dual edges.
3514 * The discrete Hodge operator is stored in hodge->matrix
3515 * This function is related to cell-based schemes
3516 *
3517 * \param[in] cm pointer to a cs_cell_mesh_t structure
3518 * \param[in, out] hodge pointer to a cs_hodge_t structure
3519 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3520 */
3521 /*----------------------------------------------------------------------------*/
3522
3523 void
cs_hodge_fped_voro_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3524 cs_hodge_fped_voro_get(const cs_cell_mesh_t *cm,
3525 cs_hodge_t *hodge,
3526 cs_cell_builder_t *cb)
3527 {
3528 CS_UNUSED(cb);
3529
3530 const cs_property_data_t *ptyd = hodge->pty_data;
3531
3532 /* Sanity check */
3533 assert(ptyd != NULL);
3534 assert(hodge->param->type == CS_HODGE_TYPE_FPED);
3535 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
3536 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ));
3537
3538 /* Initialize the local matrix related to this discrete Hodge operator */
3539 cs_sdm_t *hmat = hodge->matrix;
3540 cs_sdm_square_init(cm->n_fc, hmat);
3541
3542 if (ptyd->is_iso) {
3543
3544 for (short int f = 0; f < cm->n_fc; f++)
3545 hmat->val[f*cm->n_fc+f] =
3546 ptyd->value*cm->face[f].meas/cm->dedge[f].meas;
3547
3548 }
3549 else {
3550
3551 cs_real_3_t mv;
3552 for (short int f = 0; f < cm->n_fc; f++) {
3553
3554 const cs_nvec3_t deq = cm->dedge[f];
3555 cs_math_33_3_product(ptyd->tensor, deq.unitv, mv);
3556 hmat->val[f*cm->n_fc+f] = deq.meas*_dp3(mv, deq.unitv)/cm->face[f].meas;
3557
3558 } /* Loop on cell faces */
3559
3560 } /* Isotropic or anisotropic */
3561
3562 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3563 if (cm->c_id % CS_HODGE_MODULO == 0) {
3564 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3565 cs_sdm_dump(cm->c_id, NULL, NULL, hodge->matrix);
3566 }
3567 #endif
3568 }
3569
3570 /*----------------------------------------------------------------------------*/
3571 /*!
3572 * \brief Build a local Hodge operator for a given cell using the COST algo.
3573 * Hodge op. from primal faces to dual edges.
3574 * The discrete Hodge operator is stored in hodge->matrix
3575 * This function is related to cell-based schemes
3576 *
3577 * \param[in] cm pointer to a cs_cell_mesh_t structure
3578 * \param[in, out] hodge pointer to a cs_hodge_t structure
3579 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3580 */
3581 /*----------------------------------------------------------------------------*/
3582
3583 void
cs_hodge_fped_cost_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3584 cs_hodge_fped_cost_get(const cs_cell_mesh_t *cm,
3585 cs_hodge_t *hodge,
3586 cs_cell_builder_t *cb)
3587 {
3588 const cs_hodge_param_t *hodgep = hodge->param;
3589 const cs_property_data_t *ptyd = hodge->pty_data;
3590
3591 /* Sanity check */
3592 assert(cb != NULL && hodgep != NULL && ptyd != NULL);
3593 assert(hodgep->type == CS_HODGE_TYPE_FPED);
3594 assert(hodgep->algo == CS_HODGE_ALGO_COST);
3595 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3596
3597 cs_real_3_t *pq = cb->vectors;
3598 cs_real_3_t *dq = cb->vectors + cm->n_fc;
3599
3600 /* Initialize the geometrical quantities related to this Hodge operator */
3601 _init_fb_geom_quant(cm, pq, dq);
3602
3603 /* Compute additional geometrical quantities:
3604 Initialize the local Hodge matrix with the consistency part which is
3605 constant over a cell.
3606 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3607 and discrete Hodge operator from DUAL->PRIMAL space */
3608 double *kappa = cb->values;
3609 double *alpha = cb->values + cm->n_fc;
3610
3611 /* Initialize the local matrix related to this discrete Hodge operator */
3612 cs_sdm_t *hmat = hodge->matrix;
3613 cs_sdm_square_init(cm->n_fc, hmat);
3614
3615 if (ptyd->is_unity)
3616 _compute_cost_quant_iso(cm->n_fc, 1/cm->vol_c, 1.0,
3617 (const cs_real_t (*)[3])pq,
3618 (const cs_real_t (*)[3])dq,
3619 alpha, kappa, hmat);
3620 else if (ptyd->is_iso)
3621 _compute_cost_quant_iso(cm->n_fc, 1/cm->vol_c, ptyd->value,
3622 (const cs_real_t (*)[3])pq,
3623 (const cs_real_t (*)[3])dq,
3624 alpha, kappa, hmat);
3625 else
3626 _compute_cost_quant(cm->n_fc, 1/cm->vol_c,
3627 (const cs_real_3_t *)ptyd->tensor,
3628 (const cs_real_t (*)[3])pq,
3629 (const cs_real_t (*)[3])dq,
3630 alpha, kappa, hmat);
3631
3632 _compute_hodge_cost(cm->n_fc, hodgep->coef*hodgep->coef, alpha, kappa,
3633 hmat->val);
3634
3635 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3636 _check_vector_hodge(cm->c_id,
3637 (const cs_real_t (*)[3])pq, (const cs_real_t (*)[3])dq,
3638 hodge, cb);
3639 #endif
3640 }
3641
3642 /*----------------------------------------------------------------------------*/
3643 /*!
3644 * \brief Build a local Hodge operator for a given cell using the Bubble algo.
3645 * Hodge op. from primal faces to dual edges.
3646 * The discrete Hodge operator is stored in hodge->matrix
3647 * This function is related to cell-based schemes
3648 *
3649 * \param[in] cm pointer to a cs_cell_mesh_t structure
3650 * \param[in, out] hodge pointer to a cs_hodge_t structure
3651 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3652 */
3653 /*----------------------------------------------------------------------------*/
3654
3655 void
cs_hodge_fped_bubble_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3656 cs_hodge_fped_bubble_get(const cs_cell_mesh_t *cm,
3657 cs_hodge_t *hodge,
3658 cs_cell_builder_t *cb)
3659 {
3660 const cs_hodge_param_t *hodgep = hodge->param;
3661 const cs_property_data_t *ptyd = hodge->pty_data;
3662
3663 /* Sanity check */
3664 assert(cb != NULL && hodgep != NULL && ptyd != NULL);
3665 assert(hodgep->type == CS_HODGE_TYPE_EDFP);
3666 assert(hodgep->algo == CS_HODGE_ALGO_BUBBLE);
3667 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3668
3669 /* Initialize the geometrical quantities related to this Hodge operator */
3670 cs_real_3_t *pq = cb->vectors;
3671 cs_real_3_t *dq = cb->vectors + cm->n_fc;
3672
3673 _init_fb_geom_quant(cm, pq, dq);
3674
3675 /* Compute additional geometrical quantities:
3676 Initialize the local Hodge matrix with the consistency part which is
3677 constant over a cell.
3678 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3679 and discrete Hodge operator from DUAL->PRIMAL space */
3680 cs_sdm_t *hmat = hodge->matrix;
3681 cs_sdm_square_init(cm->n_fc, hmat);
3682
3683 if (ptyd->is_iso || ptyd->is_unity)
3684 _compute_iso_bubble_hodge_ur(cm->n_fc,
3685 hodgep->coef,
3686 1./cm->vol_c,
3687 ptyd->value,
3688 (const cs_real_t (*)[3])dq,
3689 (const cs_real_t (*)[3])pq,
3690 cb, hmat);
3691 else
3692 _compute_aniso_bubble_hodge_ur(cm->n_fc,
3693 hodgep->coef,
3694 1./cm->vol_c,
3695 ptyd->tensor,
3696 (const cs_real_t (*)[3])dq,
3697 (const cs_real_t (*)[3])pq,
3698 cb, hmat);
3699
3700 /* Hodge operator is symmetric */
3701 cs_sdm_symm_ur(hmat);
3702
3703 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3704 _check_vector_hodge(cm->c_id,
3705 (const cs_real_t (*)[3])dq, (const cs_real_t (*)[3])pq,
3706 hodge, cb);
3707 #endif
3708 }
3709
3710 /*----------------------------------------------------------------------------*/
3711 /*!
3712 * \brief Build a local Hodge operator for a given cell using VORONOI algo.
3713 * Hodge op. from dual edges to primal faces.
3714 * The discrete Hodge operator is stored in hodge->matrix
3715 * This function is related to face-based schemes
3716 *
3717 * \param[in] cm pointer to a cs_cell_mesh_t structure
3718 * \param[in, out] hodge pointer to a cs_hodge_t structure
3719 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3720 */
3721 /*----------------------------------------------------------------------------*/
3722
3723 void
cs_hodge_edfp_voro_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3724 cs_hodge_edfp_voro_get(const cs_cell_mesh_t *cm,
3725 cs_hodge_t *hodge,
3726 cs_cell_builder_t *cb)
3727 {
3728 CS_UNUSED(cb);
3729
3730 const cs_property_data_t *ptyd = hodge->pty_data;
3731
3732 /* Sanity check */
3733 assert(ptyd != NULL);
3734 assert(hodge->param->type == CS_HODGE_TYPE_EDFP);
3735 assert(hodge->param->algo == CS_HODGE_ALGO_VORONOI);
3736 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ));
3737
3738 /* Initialize the local matrix related to this discrete Hodge operator */
3739 cs_sdm_t *hmat = hodge->matrix;
3740 cs_sdm_square_init(cm->n_fc, hmat);
3741
3742 if (ptyd->is_iso) {
3743
3744 for (short int f = 0; f < cm->n_fc; f++)
3745 hmat->val[f*cm->n_fc+f] =
3746 ptyd->value*cm->face[f].meas/cm->dedge[f].meas;
3747
3748 }
3749 else {
3750
3751 cs_real_3_t mv;
3752 for (short int f = 0; f < cm->n_fc; f++) {
3753
3754 const cs_quant_t pfq = cm->face[f];
3755 cs_math_33_3_product(ptyd->tensor, pfq.unitv, mv);
3756 hmat->val[f*cm->n_fc+f] = pfq.meas * _dp3(mv, pfq.unitv)/cm->edge[f].meas;
3757
3758 } /* Loop on cell faces */
3759
3760 } /* Isotropic or anistropic */
3761
3762 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
3763 if (cm->c_id % CS_HODGE_MODULO == 0) {
3764 cs_log_printf(CS_LOG_DEFAULT, " Hodge op. ");
3765 cs_sdm_dump(cm->c_id, NULL, NULL, hdg);
3766 }
3767 #endif
3768 }
3769
3770 /*----------------------------------------------------------------------------*/
3771 /*!
3772 * \brief Build a local Hodge operator for a given cell using the COST algo.
3773 * Hodge op. from dual edges to primal faces.
3774 * The discrete Hodge operator is stored in hodge->matrix
3775 * This function is related to face-based schemes
3776 *
3777 * \param[in] cm pointer to a cs_cell_mesh_t structure
3778 * \param[in, out] hodge pointer to a cs_hodge_t structure
3779 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3780 */
3781 /*----------------------------------------------------------------------------*/
3782
3783 void
cs_hodge_edfp_cost_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3784 cs_hodge_edfp_cost_get(const cs_cell_mesh_t *cm,
3785 cs_hodge_t *hodge,
3786 cs_cell_builder_t *cb)
3787 {
3788 const cs_hodge_param_t *hodgep = hodge->param;
3789 const cs_property_data_t *ptyd = hodge->pty_data;
3790
3791 /* Sanity check */
3792 assert(cb != NULL && hodgep != NULL && ptyd != NULL);
3793 assert(hodgep->type == CS_HODGE_TYPE_EDFP);
3794 assert(hodgep->algo == CS_HODGE_ALGO_COST);
3795 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3796
3797 /* Initialize the geometrical quantities related to this Hodge operator */
3798 cs_real_3_t *pq = cb->vectors;
3799 cs_real_3_t *dq = cb->vectors + cm->n_fc;
3800
3801 _init_fb_geom_quant(cm, pq, dq);
3802
3803 /* Compute additional geometrical quantities:
3804 Initialize the local Hodge matrix with the consistency part which is
3805 constant over a cell.
3806 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3807 and discrete Hodge operator from DUAL->PRIMAL space */
3808 double *kappa = cb->values;
3809 double *alpha = cb->values + cm->n_fc;
3810
3811 /* Initialize the local matrix related to this discrete Hodge operator */
3812 cs_sdm_t *hmat = hodge->matrix;
3813 cs_sdm_square_init(cm->n_fc, hmat);
3814
3815 if (ptyd->is_unity)
3816 _compute_cost_quant_iso(cm->n_fc, 1/cm->vol_c, 1.0,
3817 (const cs_real_t (*)[3])dq,
3818 (const cs_real_t (*)[3])pq,
3819 alpha, kappa, hmat);
3820 else if (ptyd->is_iso)
3821 _compute_cost_quant_iso(cm->n_fc, 1/cm->vol_c, ptyd->value,
3822 (const cs_real_t (*)[3])dq,
3823 (const cs_real_t (*)[3])pq,
3824 alpha, kappa, hmat);
3825 else
3826 _compute_cost_quant(cm->n_fc, 1/cm->vol_c,
3827 (const cs_real_3_t *)ptyd->tensor,
3828 (const cs_real_t (*)[3])dq,
3829 (const cs_real_t (*)[3])pq,
3830 alpha, kappa, hmat);
3831
3832 _compute_hodge_cost(cm->n_fc, hodgep->coef*hodgep->coef, alpha, kappa,
3833 hmat->val);
3834
3835 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3836 _check_vector_hodge(cm->c_id,
3837 (const cs_real_t (*)[3])dq, (const cs_real_t (*)[3])pq,
3838 hodge, cb);
3839 #endif
3840 }
3841
3842 /*----------------------------------------------------------------------------*/
3843 /*!
3844 * \brief Build a local Hodge operator for a given cell using the COST algo.
3845 * Hodge op. from dual edges to primal faces.
3846 * This function is related to face-based schemes
3847 *
3848 * \param[in] cm pointer to a cs_cell_mesh_t struct.
3849 * \param[in, out] hodge pointer to a cs_hodge_t structure
3850 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3851 */
3852 /*----------------------------------------------------------------------------*/
3853
3854 void
cs_hodge_edfp_cost_get_opt(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3855 cs_hodge_edfp_cost_get_opt(const cs_cell_mesh_t *cm,
3856 cs_hodge_t *hodge,
3857 cs_cell_builder_t *cb)
3858 {
3859 const cs_hodge_param_t *hodgep = hodge->param;
3860 const cs_property_data_t *ptyd = hodge->pty_data;
3861
3862 /* Sanity check */
3863 assert(cb != NULL && hodge != NULL && ptyd != NULL);
3864 assert(hodgep->type == CS_HODGE_TYPE_EDFP);
3865 assert(hodgep->algo == CS_HODGE_ALGO_COST);
3866 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3867
3868 /* Initialize the geometrical quantities related to this Hodge operator */
3869 cs_real_3_t *pq = cb->vectors;
3870 cs_real_3_t *dq = cb->vectors + cm->n_fc;
3871
3872 _init_fb_geom_quant(cm, pq, dq);
3873
3874 /* Compute additional geometrical quantities:
3875 Initialize the local Hodge matrix with the consistency part which is
3876 constant over a cell.
3877 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3878 and discrete Hodge operator from DUAL->PRIMAL space */
3879 cs_sdm_t *hmat = hodge->matrix;
3880 cs_sdm_square_init(cm->n_fc, hmat);
3881
3882 if (ptyd->is_iso || ptyd->is_unity)
3883 _compute_iso_hodge_ur(cm->n_fc,
3884 3*hodgep->coef*hodgep->coef,
3885 1./cm->vol_c,
3886 ptyd->value,
3887 (const cs_real_t (*)[3])dq,
3888 (const cs_real_t (*)[3])pq,
3889 cb, hmat);
3890 else
3891 _compute_aniso_hodge_ur(cm->n_fc,
3892 3*hodgep->coef*hodgep->coef,
3893 1./cm->vol_c,
3894 ptyd->tensor,
3895 (const cs_real_t (*)[3])dq,
3896 (const cs_real_t (*)[3])pq,
3897 cb, hmat);
3898
3899 /* Hodge operator is symmetric */
3900 cs_sdm_symm_ur(hmat);
3901
3902 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3903 _check_vector_hodge(cm->c_id,
3904 (const cs_real_t (*)[3])dq, (const cs_real_t (*)[3])pq,
3905 hodge, cb);
3906 #endif
3907 }
3908
3909 /*----------------------------------------------------------------------------*/
3910 /*!
3911 * \brief Build a local Hodge operator for a given cell using the Bubble algo.
3912 * Hodge op. from dual edges to primal faces.
3913 * The discrete Hodge operator is stored in hodge->matrix
3914 * This function is related to face-based schemes
3915 *
3916 * \param[in] cm pointer to a cs_cell_mesh_t structure
3917 * \param[in, out] hodge pointer to a cs_hodge_t structure
3918 * \param[in, out] cb pointer to a cs_cell_builder_t structure
3919 */
3920 /*----------------------------------------------------------------------------*/
3921
3922 void
cs_hodge_edfp_bubble_get(const cs_cell_mesh_t * cm,cs_hodge_t * hodge,cs_cell_builder_t * cb)3923 cs_hodge_edfp_bubble_get(const cs_cell_mesh_t *cm,
3924 cs_hodge_t *hodge,
3925 cs_cell_builder_t *cb)
3926 {
3927 const cs_hodge_param_t *hodgep = hodge->param;
3928 const cs_property_data_t *ptyd = hodge->pty_data;
3929
3930 /* Sanity check */
3931 assert(cb != NULL && hodgep != NULL && ptyd != NULL);
3932 assert(hodgep->type == CS_HODGE_TYPE_EDFP);
3933 assert(hodgep->algo == CS_HODGE_ALGO_BUBBLE);
3934 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
3935
3936 /* Initialize the geometrical quantities related to this Hodge operator */
3937 cs_real_3_t *pq = cb->vectors;
3938 cs_real_3_t *dq = cb->vectors + cm->n_fc;
3939
3940 _init_fb_geom_quant(cm, pq, dq);
3941
3942 /* Compute additional geometrical quantities:
3943 Initialize the local Hodge matrix with the consistency part which is
3944 constant over a cell.
3945 Switch arguments between discrete Hodge operator from PRIMAL->DUAL space
3946 and discrete Hodge operator from DUAL->PRIMAL space */
3947 cs_sdm_t *hmat = hodge->matrix;
3948 cs_sdm_square_init(cm->n_fc, hmat);
3949
3950 if (ptyd->is_iso || ptyd->is_unity)
3951 _compute_iso_bubble_hodge_ur(cm->n_fc,
3952 hodgep->coef,
3953 1./cm->vol_c,
3954 ptyd->value,
3955 (const cs_real_t (*)[3])dq,
3956 (const cs_real_t (*)[3])pq,
3957 cb, hmat);
3958 else
3959 _compute_aniso_bubble_hodge_ur(cm->n_fc,
3960 hodgep->coef,
3961 1./cm->vol_c,
3962 ptyd->tensor,
3963 (const cs_real_t (*)[3])dq,
3964 (const cs_real_t (*)[3])pq,
3965 cb, hmat);
3966
3967 /* Hodge operator is symmetric */
3968 cs_sdm_symm_ur(hmat);
3969
3970 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 0
3971 _check_vector_hodge(cm->c_id,
3972 (const cs_real_t (*)[3])dq, (const cs_real_t (*)[3])pq,
3973 hodge, cb);
3974 #endif
3975 }
3976
3977 /*----------------------------------------------------------------------------*/
3978 /*!
3979 * \brief Compute cellwise a discrete hodge operator and multiple it with
3980 * a vector
3981 *
3982 * \param[in] connect pointer to a cs_cdo_connect_t structure
3983 * \param[in] quant pointer to a cs_cdo_quantities_t structure
3984 * \param[in] hodgep cs_hodge_param_t structure
3985 * \param[in] pty pointer to a cs_property_t structure or NULL
3986 * \param[in] in_vals vector to multiply with the discrete Hodge op.
3987 * \param[in] t_eval time at which one performs the evaluation
3988 * \param[in, out] result array storing the resulting matrix-vector product
3989 */
3990 /*----------------------------------------------------------------------------*/
3991
3992 void
cs_hodge_matvec(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_hodge_param_t hodgep,const cs_property_t * pty,const cs_real_t in_vals[],cs_real_t t_eval,cs_real_t result[])3993 cs_hodge_matvec(const cs_cdo_connect_t *connect,
3994 const cs_cdo_quantities_t *quant,
3995 const cs_hodge_param_t hodgep,
3996 const cs_property_t *pty,
3997 const cs_real_t in_vals[],
3998 cs_real_t t_eval,
3999 cs_real_t result[])
4000 {
4001 if (in_vals == NULL)
4002 return;
4003
4004 const char *func_name = __func__;
4005
4006 if (result == NULL) {
4007 bft_error(__FILE__, __LINE__, 0,
4008 "%s: Resulting vector must be allocated", __func__);
4009 return; /* Avoid a warning */
4010 }
4011 assert(connect != NULL && quant != NULL); /* Sanity checks */
4012
4013 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
4014 shared(quant, connect, in_vals, t_eval, result, pty, func_name) \
4015 firstprivate(hodgep)
4016 {
4017 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
4018 int t_id = omp_get_thread_num();
4019 #else
4020 int t_id = 0;
4021 #endif
4022
4023 cs_eflag_t msh_flag = 0;
4024 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
4025 cs_cell_builder_t *cb = NULL;
4026 double *_in = NULL;
4027
4028 /* Each thread has its own pointer to a cs_hodge_t structure */
4029 cs_hodge_t *hodge = cs_hodge_create(connect, pty, &hodgep, true, false);
4030 cs_hodge_compute_t *compute = cs_hodge_get_func(func_name, hodgep);
4031 bool pty_uniform = cs_property_is_uniform(pty);
4032
4033 switch (hodgep.type) {
4034
4035 case CS_HODGE_TYPE_VPCD:
4036
4037 msh_flag |= CS_FLAG_COMP_PVQ;
4038 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOVB, connect);
4039 BFT_MALLOC(_in, connect->n_max_vbyc, double);
4040
4041 # pragma omp for
4042 for (cs_lnum_t i = 0; i < quant->n_vertices; i++) result[i] = 0;
4043
4044 switch (hodgep.algo) {
4045
4046 case CS_HODGE_ALGO_WBS:
4047 msh_flag |= CS_FLAG_COMP_PVQ |CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
4048 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_HFQ;
4049 break;
4050
4051 default: /* Minimal requirement for other Hodge algorithms */
4052 msh_flag |= CS_FLAG_COMP_PVQ;
4053 break;
4054
4055 }
4056 break;
4057
4058 case CS_HODGE_TYPE_EPFD:
4059
4060 msh_flag |= CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ;
4061 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOVB, connect);
4062 BFT_MALLOC(_in, connect->n_max_ebyc, double);
4063
4064 # pragma omp for
4065 for (cs_lnum_t i = 0; i < quant->n_edges; i++) result[i] = 0;
4066
4067 switch (hodgep.algo) {
4068
4069 case CS_HODGE_ALGO_OCS2:
4070 msh_flag |= CS_FLAG_COMP_EV | CS_FLAG_COMP_SEF;
4071 break;
4072 case CS_HODGE_ALGO_VORONOI:
4073 msh_flag |= CS_FLAG_COMP_SEF;
4074 break;
4075
4076 default: /* Minimal requirement for other Hodge algorithms */
4077 break;
4078
4079 }
4080 break;
4081
4082 case CS_HODGE_TYPE_EDFP:
4083
4084 msh_flag |= CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ;
4085 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOFB, connect);
4086 BFT_MALLOC(_in, connect->n_max_fbyc, double);
4087
4088 # pragma omp for
4089 for (cs_lnum_t i = 0; i < quant->n_faces; i++) result[i] = 0;
4090
4091 break;
4092
4093 case CS_HODGE_TYPE_FPED:
4094
4095 msh_flag |= CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ;
4096 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOFB, connect);
4097 BFT_MALLOC(_in, connect->n_max_fbyc, double);
4098
4099 # pragma omp for
4100 for (cs_lnum_t i = 0; i < quant->n_faces; i++) result[i] = 0;
4101
4102 break;
4103
4104 case CS_HODGE_TYPE_VC:
4105
4106 msh_flag |= CS_FLAG_COMP_PVQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_HFQ |
4107 CS_FLAG_COMP_DEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ;
4108 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOVCB, connect);
4109 BFT_MALLOC(_in, connect->n_max_vbyc + 1, double);
4110
4111 # pragma omp for
4112 for (cs_lnum_t i = 0; i < quant->n_vertices + quant->n_cells; i++)
4113 result[i] = 0;
4114
4115 break;
4116
4117 default:
4118 bft_error(__FILE__, __LINE__, 0,
4119 " %s: Invalid type of discrete Hodge operator", func_name);
4120 }
4121
4122 const cs_flag_t cell_flag = 0;
4123 if (pty_uniform) /* Get the value from the first cell */
4124 cs_hodge_set_property_value(0, t_eval, cell_flag, hodge);
4125
4126 # pragma omp for CS_CDO_OMP_SCHEDULE
4127 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
4128
4129 /* Set the local mesh structure for the current cell */
4130 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
4131
4132 /* Retrieve the value of the property inside the current cell */
4133 if (!pty_uniform)
4134 cs_hodge_set_property_value_cw(cm, t_eval, cell_flag, hodge);
4135
4136 /* Build the local discrete Hodge operator */
4137 compute(cm, hodge, cb);
4138
4139 /* Define a local vector to multiply for the current cell */
4140 switch (hodgep.type) {
4141
4142 case CS_HODGE_TYPE_VPCD:
4143 for (short int v = 0; v < cm->n_vc; v++)
4144 _in[v] = in_vals[cm->v_ids[v]];
4145
4146 /* Local matrix-vector operation */
4147 cs_sdm_square_matvec(hodge->matrix, _in, cb->values);
4148
4149 /* Assemble the resulting vector */
4150 for (short int v = 0; v < cm->n_vc; v++)
4151 # pragma omp atomic
4152 result[cm->v_ids[v]] += cb->values[v];
4153 break;
4154
4155 case CS_HODGE_TYPE_EPFD:
4156 for (short int e = 0; e < cm->n_ec; e++)
4157 _in[e] = in_vals[cm->e_ids[e]];
4158
4159 /* Local matrix-vector operation */
4160 cs_sdm_square_matvec(hodge->matrix, _in, cb->values);
4161
4162 /* Assemble the resulting vector */
4163 for (short int e = 0; e < cm->n_ec; e++)
4164 # pragma omp atomic
4165 result[cm->e_ids[e]] += cb->values[e];
4166 break;
4167
4168 case CS_HODGE_TYPE_FPED:
4169 for (short int f = 0; f < cm->n_fc; f++)
4170 _in[f] = in_vals[cm->f_ids[f]];
4171
4172 /* Local matrix-vector operation */
4173 cs_sdm_square_matvec(hodge->matrix, _in, cb->values);
4174
4175 /* Assemble the resulting vector */
4176 for (short int f = 0; f < cm->n_fc; f++)
4177 # pragma omp atomic
4178 result[cm->f_ids[f]] += cb->values[f];
4179 break;
4180
4181 default:
4182 bft_error(__FILE__, __LINE__, 0,
4183 " %s: Invalid type of discrete Hodge operator", func_name);
4184
4185 } /* Hodge type */
4186
4187 } /* Main loop on cells */
4188
4189 BFT_FREE(_in);
4190 cs_cell_builder_free(&cb);
4191 cs_hodge_free(&hodge);
4192
4193 } /* OpenMP Block */
4194
4195 }
4196
4197 /*----------------------------------------------------------------------------*/
4198 /*!
4199 * \brief Compute cellwise a discrete hodge operator in order to define
4200 * a circulation array from a flux array
4201 *
4202 * \param[in] connect pointer to a cs_cdo_connect_t structure
4203 * \param[in] quant pointer to a cs_cdo_quantities_t structure
4204 * \param[in] t_eval time at which one performs the evaluation
4205 * \param[in] hodgep cs_hodge_param_t structure
4206 * \param[in] pty pointer to a cs_property_t structure or NULL
4207 * \param[in] flux vector to multiply with the discrete Hodge op.
4208 * \param[in, out] circul array storing the resulting matrix-vector product
4209 */
4210 /*----------------------------------------------------------------------------*/
4211
4212 void
cs_hodge_circulation_from_flux(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,cs_real_t t_eval,const cs_hodge_param_t hodgep,const cs_property_t * pty,const cs_real_t flux[],cs_real_t circul[])4213 cs_hodge_circulation_from_flux(const cs_cdo_connect_t *connect,
4214 const cs_cdo_quantities_t *quant,
4215 cs_real_t t_eval,
4216 const cs_hodge_param_t hodgep,
4217 const cs_property_t *pty,
4218 const cs_real_t flux[],
4219 cs_real_t circul[])
4220 {
4221 if (flux == NULL)
4222 return;
4223
4224 const char *func_name = __func__;
4225
4226 if (circul == NULL) {
4227 bft_error(__FILE__, __LINE__, 0,
4228 "%s: Resulting vector must be allocated", __func__);
4229 return; /* Avoid a warning */
4230 }
4231 assert(connect != NULL && quant != NULL); /* Sanity checks */
4232
4233 # pragma omp parallel if (quant->n_cells > CS_THR_MIN) \
4234 shared(quant, connect, flux, t_eval, circul, pty, func_name) \
4235 firstprivate(hodgep)
4236 {
4237 #if defined(HAVE_OPENMP) /* Determine default number of OpenMP threads */
4238 int t_id = omp_get_thread_num();
4239 #else
4240 int t_id = 0;
4241 #endif
4242
4243 const cs_adjacency_t *c2f = connect->c2f;
4244 const cs_flag_t cell_flag = 0;
4245
4246 cs_eflag_t msh_flag = 0;
4247 cs_cell_mesh_t *cm = cs_cdo_local_get_cell_mesh(t_id);
4248 cs_cell_builder_t *cb = NULL;
4249 double *_fluxes = NULL;
4250
4251 /* Each thread has its own pointer to a cs_hodge_t structure */
4252 cs_hodge_t *hodge = cs_hodge_create(connect, pty, &hodgep, true, false);
4253 cs_hodge_compute_t *compute = cs_hodge_get_func(func_name, hodgep);
4254 bool pty_uniform = cs_property_is_uniform(pty);
4255
4256 switch (hodgep.type) {
4257
4258 /* Only this type of discrete Hodge operator makes sense for this
4259 operation */
4260 case CS_HODGE_TYPE_FPED:
4261
4262 msh_flag |= CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ;
4263 cb = _cell_builder_create(CS_SPACE_SCHEME_CDOFB, connect);
4264 BFT_MALLOC(_fluxes, connect->n_max_fbyc, double);
4265
4266 break;
4267
4268 default:
4269 bft_error(__FILE__, __LINE__, 0,
4270 " %s: Invalid type of discrete Hodge operator", func_name);
4271 }
4272
4273 if (pty_uniform) /* Get the value from the first cell */
4274 cs_hodge_set_property_value(0, t_eval, cell_flag, hodge);
4275
4276 # pragma omp for CS_CDO_OMP_SCHEDULE
4277 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
4278
4279 /* Set the local mesh structure for the current cell */
4280 cs_cell_mesh_build(c_id, msh_flag, connect, quant, cm);
4281
4282 /* Retrieve the value of the property inside the current cell */
4283 if (!pty_uniform)
4284 cs_hodge_set_property_value_cw(cm, t_eval, cell_flag, hodge);
4285
4286 /* Build the local discrete Hodge operator */
4287 compute(cm, hodge, cb);
4288
4289 /* Define a local vector to multiply for the current cell */
4290 for (short int f = 0; f < cm->n_fc; f++)
4291 _fluxes[f] = flux[cm->f_ids[f]];
4292 cs_real_t *_circ = circul + c2f->idx[c_id];
4293
4294 /* Local matrix-vector operation */
4295 cs_sdm_square_matvec(hodge->matrix, _fluxes, _circ);
4296
4297 } /* Main loop on cells */
4298
4299 BFT_FREE(_fluxes);
4300 cs_cell_builder_free(&cb);
4301 cs_hodge_free(&hodge);
4302
4303 } /* OpenMP Block */
4304
4305 }
4306
4307 /*----------------------------------------------------------------------------*/
4308 /*!
4309 * \brief Compute the hodge operator related to a face (i.e. a mass matrix
4310 * with unity property) using a Whitney Barycentric Subdivision (WBS)
4311 * algorithm
4312 *
4313 * \param[in] fm pointer to a cs_face_mesh_t structure
4314 * \param[in, out] hf pointer to a cs_sdm_t structure to define
4315 */
4316 /*----------------------------------------------------------------------------*/
4317
4318 void
cs_hodge_compute_wbs_surfacic(const cs_face_mesh_t * fm,cs_sdm_t * hf)4319 cs_hodge_compute_wbs_surfacic(const cs_face_mesh_t *fm,
4320 cs_sdm_t *hf)
4321 {
4322 assert(hf != NULL && fm != NULL);
4323
4324 /* Reset values */
4325 cs_sdm_square_init(fm->n_vf, hf);
4326
4327 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
4328
4329 double *hi = hf->val + vfi*hf->n_rows;
4330
4331 /* Default contribution */
4332 const double default_coef = 0.5 * fm->wvf[vfi] * fm->face.meas;
4333 for (short int vfj = 0; vfj < fm->n_vf; vfj++)
4334 hi[vfj] = default_coef * fm->wvf[vfj];
4335
4336 /* Specific diagonal contribution */
4337 hi[vfi] += 2 * default_coef * cs_math_1ov3;
4338
4339 } /* Loop on face vertices */
4340
4341 /* Specific extra-diag contribution */
4342 for (short int e = 0; e < fm->n_ef; e++) {
4343
4344 const short int v1 = fm->e2v_ids[2*e];
4345 const short int v2 = fm->e2v_ids[2*e+1];
4346 const double extra_val = cs_math_1ov12 * fm->tef[e];
4347
4348 hf->val[v1*hf->n_rows + v2] += extra_val;
4349 hf->val[v2*hf->n_rows + v1] += extra_val;
4350
4351 } /* Loop on face edges */
4352
4353 #if defined(DEBUG) && !defined(NDEBUG) && CS_HODGE_DBG > 1
4354 if (fm->c_id % CS_HODGE_MODULO == 0) {
4355 cs_log_printf(CS_LOG_DEFAULT, " Surfacic Hodge op. ");
4356 cs_sdm_dump(fm->f_id, NULL, NULL, hf);
4357 }
4358 #endif
4359 }
4360
4361 /*----------------------------------------------------------------------------*/
4362
4363 #undef _dp3
4364
4365 END_C_DECLS
4366