1 /*============================================================================
2 * Definitions des fonctions
3 * associees a la structure `ecs_table_t' decrivant un table
4 * et propres aux tables principaux de type "definition"
5 *============================================================================*/
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
12 This program is free software; you can redistribute it and/or modify it under
13 the terms of the GNU General Public License as published by the Free Software
14 Foundation; either version 2 of the License, or (at your option) any later
15 version.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 details.
21
22 You should have received a copy of the GNU General Public License along with
23 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24 Street, Fifth Floor, Boston, MA 02110-1301, USA.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29
30 /*============================================================================
31 * Visibilite
32 *============================================================================*/
33
34 /*----------------------------------------------------------------------------
35 * Fichiers `include' librairie standard C
36 *----------------------------------------------------------------------------*/
37
38 #include <assert.h>
39 #include <stdlib.h>
40 #include <string.h>
41
42
43 /*----------------------------------------------------------------------------
44 * Fichiers `include' visibles du paquetage global "Utilitaire"
45 *----------------------------------------------------------------------------*/
46
47 #include "ecs_elt_typ_liste.h"
48 #include "ecs_def.h"
49 #include "ecs_mem.h"
50 #include "ecs_tab.h"
51
52
53 /*----------------------------------------------------------------------------
54 * Fichiers `include' visibles des paquetages visibles
55 *----------------------------------------------------------------------------*/
56
57
58 /*----------------------------------------------------------------------------
59 * Fichiers `include' visibles du paquetage courant
60 *----------------------------------------------------------------------------*/
61
62 #include "ecs_table.h"
63
64
65 /*----------------------------------------------------------------------------
66 * Fichier `include' du paquetage courant associe au fichier courant
67 *----------------------------------------------------------------------------*/
68
69 #include "ecs_table_def.h"
70
71
72 /*----------------------------------------------------------------------------
73 * Fichiers `include' prives du paquetage courant
74 *----------------------------------------------------------------------------*/
75
76 #include "ecs_table_priv.h"
77
78
79 /*============================================================================
80 * Macros globales au fichier
81 *============================================================================*/
82
83 /* Longueur maximale du nom d'un type d'élément (+1 pour le `\0' !) */
84 #define ECS_LOC_LNG_MAX_NOM_TYP 11
85
86 #if !defined(FLT_MAX)
87 #define FLT_MAX HUGE_VAL
88 #endif
89
90 #define ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) ( \
91 prod_vect[0] = vect1[1] * vect2[2] - vect2[1] * vect1[2], \
92 prod_vect[1] = vect2[0] * vect1[2] - vect1[0] * vect2[2], \
93 prod_vect[2] = vect1[0] * vect2[1] - vect2[0] * vect1[1] )
94
95 #define ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) ( \
96 vect1[0] * vect2[0] + vect1[1] * vect2[1] + vect1[2] * vect2[2] )
97
98 #define ECS_LOC_MODULE(vect) \
99 sqrt(vect[0] * vect[0] + vect[1] * vect[1] + vect[2] * vect[2])
100
101 #define ECS_LOC_DETERMINANT(vect1, vect2, vect3) ( \
102 ((vect1[1] * vect2[2] - vect2[1] * vect1[2]) * vect3[0]) \
103 + ((vect2[0] * vect1[2] - vect1[0] * vect2[2]) * vect3[1]) \
104 + ((vect1[0] * vect2[1] - vect2[0] * vect1[1]) * vect3[2]) )
105
106 /*============================================================================
107 * Fonctions privees
108 *============================================================================*/
109
110 /*----------------------------------------------------------------------------
111 * Fonction qui met à jour la définition faces -> sommets en cas
112 * de fusion de sommets.
113 *----------------------------------------------------------------------------*/
114
115 static void
_table_def__maj_fac_som(ecs_table_t * table_def_fac,const ecs_tab_int_t * tab_som_old_new)116 _table_def__maj_fac_som(ecs_table_t *table_def_fac,
117 const ecs_tab_int_t *tab_som_old_new)
118 {
119 size_t cpt_val;
120 size_t nbr_fac;
121 size_t nbr_fac_mod;
122 size_t nbr_val_ini, nbr_val_fin;
123 size_t nbr_som_fac;
124
125 size_t num_som, num_som_prev;
126
127 size_t ind_fac;
128 size_t ind_fac_mod;
129 size_t ind_som;
130
131 size_t ipos_deb, ipos_fin;
132
133 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
134
135 /* Initialisations */
136 /* --------------- */
137
138 nbr_fac = table_def_fac->nbr;
139
140 nbr_val_ini = table_def_fac->pos[nbr_fac];
141
142 /* Mise à jour de la définition des faces */
143 /*----------------------------------------*/
144
145 for (ind_som = 0;
146 ind_som < table_def_fac->pos[table_def_fac->nbr] - 1;
147 ind_som++)
148
149 table_def_fac->val[ind_som]
150 = tab_som_old_new->val[table_def_fac->val[ind_som] - 1];
151
152 /* Suppression de sommets confondus de la définition des faces */
153 /*-------------------------------------------------------------*/
154
155 nbr_fac_mod = 0;
156
157 cpt_val = 0;
158 ipos_deb = 1;
159
160 for (ind_fac = 0; ind_fac < nbr_fac; ind_fac++) {
161
162 ind_fac_mod = 0;
163
164 ipos_fin = table_def_fac->pos[ind_fac + 1] - 1;
165
166 nbr_som_fac = ipos_fin - ipos_deb;
167
168 num_som_prev = table_def_fac->val[ipos_deb + nbr_som_fac - 1];
169
170 for (ind_som = 0; ind_som < nbr_som_fac; ind_som++) {
171
172 num_som = table_def_fac->val[ipos_deb + ind_som];
173
174 if (num_som != num_som_prev) {
175 num_som_prev = num_som;
176 table_def_fac->val[cpt_val++] = num_som;
177 }
178 else
179 ind_fac_mod = 1;
180 }
181
182 table_def_fac->pos[ind_fac + 1] = cpt_val + 1;
183
184 ipos_deb = ipos_fin;
185
186 nbr_fac_mod += ind_fac_mod;
187 }
188
189 nbr_val_fin = table_def_fac->pos[nbr_fac];
190
191 assert(nbr_val_fin <= nbr_val_ini);
192
193 if (nbr_val_fin != nbr_val_ini) {
194
195 ECS_REALLOC(table_def_fac->val, nbr_val_fin, ecs_int_t);
196 printf(_("\nMesh verification:\n\n"
197 " %d faces modified due to merged vertices.\n"),
198 (int)nbr_fac_mod);
199 }
200 }
201
202 /*----------------------------------------------------------------------------
203 * Fonction qui transforme un tableau d'équivalence en liste chaînée simple
204 *----------------------------------------------------------------------------*/
205
206 static void
_table_def__transf_equiv(size_t nbr_som,ecs_tab_int_t * tab_equiv_som)207 _table_def__transf_equiv(size_t nbr_som,
208 ecs_tab_int_t *tab_equiv_som)
209 {
210 size_t ind_som;
211 ecs_int_t ind_som_min;
212 ecs_int_t ind_som_max;
213 ecs_int_t ind_som_tmp;
214
215 ecs_tab_int_t tab_equiv_prec;
216
217 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
218
219 tab_equiv_prec.nbr = nbr_som;
220 ECS_MALLOC(tab_equiv_prec.val, tab_equiv_prec.nbr, ecs_int_t);
221
222 for (ind_som = 0; ind_som < nbr_som; ind_som++)
223 tab_equiv_prec.val[ind_som] = -1;
224
225 for (ind_som = 0; ind_som < nbr_som; ind_som++) {
226
227 if (tab_equiv_som->val[ind_som] != -1) {
228
229 ind_som_min = ind_som;
230 ind_som_max = tab_equiv_som->val[ind_som];
231
232 assert (ind_som_min < ind_som_max);
233
234 while ( ind_som_min != ind_som_max
235 && tab_equiv_prec.val[ind_som_max] != ind_som_min) {
236
237 /*
238 On parcourt la liste inverse correspondant à ind_som_max jusqu'au
239 point d'insertion (si pas de point precedent dans la chaine,
240 tab_equiv_prec.val[ind_som_max] = -1).
241 */
242
243 while (tab_equiv_prec.val[ind_som_max] > ind_som_min)
244 ind_som_max = tab_equiv_prec.val[ind_som_max];
245
246 ind_som_tmp = tab_equiv_prec.val[ind_som_max];
247
248 /*
249 Si l'on est en début de chaîne, on branche la liste inverse
250 correspondant à ind_som_min au début de celle correspondant
251 à ind_som_max. Sinon, on doit reboucler.
252 */
253
254 tab_equiv_prec.val[ind_som_max] = ind_som_min;
255
256 if (ind_som_tmp != -1) {
257 ind_som_max = ind_som_min;
258 ind_som_min = ind_som_tmp;
259 }
260
261 }
262 }
263 }
264
265 for (ind_som = 0; ind_som < nbr_som; ind_som++) {
266
267 if (tab_equiv_prec.val[ind_som] != -1)
268 tab_equiv_som->val[tab_equiv_prec.val[ind_som]] = ind_som;
269
270 }
271
272 tab_equiv_prec.nbr = 0;
273 ECS_FREE(tab_equiv_prec.val);
274 }
275
276 /*----------------------------------------------------------------------------
277 * Fonction qui fusionne des sommets équivalents
278 *
279 * Remarque : le tableau d'équivalence (fusion) des sommets est construit de
280 * manière à ce qu'à un sommet ne comportant pas d'équivalent
281 * où de plus petit indice parmi ses équivalents corresponde la
282 * valeur -1, alors qu'un un sommet possédant des équivalents de
283 * plus petit indice corresponde le plus grand indice parmi ces
284 * équivalents (ce qui constitue une sorte de liste chaînée).
285 *----------------------------------------------------------------------------*/
286
287 static ecs_tab_int_t
_table_def__fusion_som(size_t * n_vertices,ecs_coord_t ** coords,ecs_tab_int_t * tab_equiv_som)288 _table_def__fusion_som(size_t *n_vertices,
289 ecs_coord_t **coords,
290 ecs_tab_int_t *tab_equiv_som)
291 {
292 size_t ind_som, nbr_som_old, nbr_som_new, nbr_som_fus;
293 ecs_int_t ind_som_loc, ind_som_tmp, icoo;
294 ecs_coord_t som_poids;
295 double som_coord[3];
296
297 ecs_tab_int_t tab_som_old_new;
298
299 /* Initialization */
300 /* -------------- */
301
302 nbr_som_old = *n_vertices;
303
304 printf(_("\n Merging vertices:\n"));
305
306 /* Transform vertex equivalences array into simple linked list. */
307
308 _table_def__transf_equiv(nbr_som_old, tab_equiv_som);
309
310 printf(_(" Initial number of vertices : %10d\n"), (int)nbr_som_old);
311
312 tab_som_old_new.nbr = nbr_som_old;
313 ECS_MALLOC(tab_som_old_new.val, tab_som_old_new.nbr, ecs_int_t);
314
315 /* Initialize vertex renumbering array */
316
317 for (ind_som = 0; ind_som < nbr_som_old; ind_som++)
318 tab_som_old_new.val[ind_som] = 0;
319
320 /* Main loop on vertices */
321 /*-----------------------*/
322
323 nbr_som_new = 0;
324
325 for (ind_som = 0; ind_som < nbr_som_old; ind_som++) {
326
327 /* If the vertex has not been handled yet */
328
329 if (tab_som_old_new.val[ind_som] == 0) {
330
331 /* Initialize vertex */
332
333 som_poids = 0.0;
334
335 for (icoo = 0; icoo < 3; icoo++)
336 som_coord[icoo] = 0.0;
337
338 ind_som_loc = ind_som;
339
340 /* Mark equivalent vertices and compute their contribution */
341
342 do {
343
344 tab_som_old_new.val[ind_som_loc] = nbr_som_new + 1;
345
346 som_poids += 1.0;
347
348 for (icoo = 0; icoo < 3; icoo++)
349 som_coord[icoo] += (* coords)[3*ind_som_loc + icoo];
350
351 ind_som_loc = tab_equiv_som->val[ind_som_loc];
352
353 } while (ind_som_loc != -1);
354
355 /* Final coordinates */
356
357 for (icoo = 0; icoo < 3; icoo++)
358 (*coords)[3*nbr_som_new + icoo] = som_coord[icoo] / som_poids;
359
360 /* Do not forget to increment the number of vertices after merge */
361
362 nbr_som_new += 1;
363 }
364
365 }
366
367 /* End of main loop */
368 /* ---------------- */
369
370 /* We now know the number of vertices after merging */
371
372 *n_vertices = nbr_som_new;
373 ECS_REALLOC(*coords, nbr_som_new*3, ecs_coord_t);
374
375 printf(_(" Number of vertices after merging : %10d\n"), (int)nbr_som_new);
376
377 /* Mark vertices originating from a merge */
378 /*----------------------------------------*/
379
380 nbr_som_fus = 0;
381
382 /*
383 We will not need the vertex equivalence array any more; we thus modify
384 it so that only the first vertex of a linked equivalence list points
385 to its first equivalent, to use it as a marker.
386 */
387
388 for (ind_som = 0; ind_som < tab_equiv_som->nbr; ind_som++) {
389
390 if (tab_equiv_som->val[ind_som] != -1) {
391
392 nbr_som_fus += 1;
393
394 ind_som_loc = tab_equiv_som->val[ind_som];
395
396 while (ind_som_loc != -1) {
397
398 ind_som_tmp = ind_som_loc;
399 ind_som_loc = tab_equiv_som->val[ind_som_loc];
400
401 tab_equiv_som->val[ind_som_tmp] = -1;
402
403 }
404 }
405 }
406
407 printf(_(" Number of modified vertices : %10d\n"), (int)nbr_som_fus);
408
409 /* Return */
410
411 return tab_som_old_new;
412 }
413
414 /*----------------------------------------------------------------------------
415 * Fonction qui met à jour le tableau des fusions de sommets
416 *
417 * Remarque : le tableau d'équivalence (fusion) des sommets est construit de
418 * manière à ce qu'à un sommet ne comportant pas d'équivalent
419 * où de plus grand indice parmi ses équivalents corresponde la
420 * valeur -1, alors qu'un un sommet possédant des équivalents de
421 * plus grand indice corresponde le plus petit indice parmi ces
422 * équivalents (ce qui constitue une sorte de liste chaînée).
423 *----------------------------------------------------------------------------*/
424
425 static void
_table_def__maj_equiv_som(size_t ind_som_0,size_t ind_som_1,ecs_tab_int_t * tab_equiv_som)426 _table_def__maj_equiv_som(size_t ind_som_0,
427 size_t ind_som_1,
428 ecs_tab_int_t *tab_equiv_som)
429 {
430 ecs_int_t ind_som_max;
431 ecs_int_t ind_som_min;
432 ecs_int_t ind_som_tmp;
433
434 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
435
436 ind_som_min = ECS_MIN(ind_som_0, ind_som_1);
437 ind_som_max = ECS_MAX(ind_som_0, ind_som_1);
438
439 /* On fusionne deux listes chaînées */
440 /*----------------------------------*/
441
442 while ( ind_som_max != ind_som_min
443 && tab_equiv_som->val[ind_som_min] != ind_som_max) {
444
445 /*
446 On parcourt la liste chaînée correspondant à ind_som_min jusqu'au
447 point d'insertion (si pas de point suivant dans la chaine,
448 tab_equiv_som->val[ind_som_min] = -1).
449 */
450
451 while ( tab_equiv_som->val[ind_som_min] != -1
452 && tab_equiv_som->val[ind_som_min] < ind_som_max)
453 ind_som_min = tab_equiv_som->val[ind_som_min];
454
455 ind_som_tmp = tab_equiv_som->val[ind_som_min];
456
457 /*
458 Si l'on est en fin de chaîne, on branche la liste chaînée correspondant
459 à ind_som_max à la suite de celle correspondant à ind_som_min.
460 Sinon, on doit reboucler
461 */
462
463 tab_equiv_som->val[ind_som_min] = ind_som_max;
464
465 if (ind_som_tmp != -1) {
466 ind_som_min = ind_som_max;
467 ind_som_max = ind_som_tmp;
468 }
469 }
470 }
471
472 /*----------------------------------------------------------------------------
473 * Correction si nécessaire de l'orientation d'un quadrangle en connectivité
474 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
475 * ou que l'on demande une simple vérification (i.e. correc = false),
476 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
477 * l'orientation par une permutation de la connectivité locale.
478 *----------------------------------------------------------------------------*/
479
480 static ecs_int_t
_orient_quad(const ecs_coord_t coord[],ecs_int_t connect[],const bool correc)481 _orient_quad(const ecs_coord_t coord[],
482 ecs_int_t connect[],
483 const bool correc)
484 {
485 ecs_int_t isom_tmp;
486 ecs_int_t itri;
487 ecs_int_t nb_angle_obtus;
488 ecs_coord_t sgn;
489 ecs_coord_t normale[3];
490 ecs_coord_t vect1[3];
491 ecs_coord_t vect2[3];
492 ecs_coord_t prod_vect[3];
493
494 ecs_int_t passage = -1;
495
496 #define ECS_LOC_INIT_VECT(vect, i, j) ( \
497 vect[0] = coord[((connect[j-1] - 1) * 3) ] \
498 - coord[((connect[i-1] - 1) * 3) ], \
499 vect[1] = coord[((connect[j-1] - 1) * 3) + 1] \
500 - coord[((connect[i-1] - 1) * 3) + 1], \
501 vect[2] = coord[((connect[j-1] - 1) * 3) + 2] \
502 - coord[((connect[i-1] - 1) * 3) + 2] )
503
504 /* Calcul d'une direction normale approchée de la face
505 (peut être entrante ou sortante si la face est "croisée") */
506
507 ECS_LOC_INIT_VECT(vect1, 1, 2);
508 ECS_LOC_INIT_VECT(vect2, 1, 4);
509
510 ECS_LOC_PRODUIT_VECTORIEL(normale, vect1, vect2);
511
512 /* Boucle sur les renumérotations possibles */
513
514 do {
515
516 passage += 1;
517
518 nb_angle_obtus = 0;
519
520 /* Initialisation */
521
522 switch(passage) {
523
524 case 0:
525 break;
526
527 case 1:
528 if (correc == false)
529 return -1;
530 isom_tmp = connect[2];
531 connect[2] = connect[3];
532 connect[3] = isom_tmp;
533 break;
534
535 default:
536 return -1;
537
538 }
539
540 /* Boucle sur les coins */
541
542 /* On compte les angles obtus, qui devraient être au nombre de 2 sur
543 une face "croisée" (et 0 sur une face bien définie et convexe
544 1 sur une face bien définie non convexe, soit 3 en apparaence
545 sur une face bien définie convexe si la non-convexité se trouve
546 au niveau du sommet 1 et que l'on a donc calculé une normale
547 "à l'envers"). */
548
549 for (itri = 0; itri < 4; itri++) {
550
551 ECS_LOC_INIT_VECT(vect1, ((itri+2) % 4) + 1, ((itri+1) % 4) + 1);
552 ECS_LOC_INIT_VECT(vect2, ( itri % 4) + 1, ((itri+1) % 4) + 1);
553
554 ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2);
555
556 /* Angle obtus si produit mixte < 0, aigu sinon. */
557
558 sgn = ECS_LOC_PRODUIT_SCALAIRE(prod_vect, normale);
559
560 if (sgn < 0.0)
561 nb_angle_obtus += 1;
562
563 }
564
565 } while (nb_angle_obtus == 2);
566
567 return (ECS_MIN(passage, 1));
568
569 #undef ECS_LOC_INIT_VECT
570 }
571
572 /*----------------------------------------------------------------------------
573 * Correction si nécessaire de l'orientation d'un tétraèdre en connectivité
574 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation,
575 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
576 * l'orientation par une permutation de la connectivité locale.
577 *----------------------------------------------------------------------------*/
578
579 static ecs_int_t
_orient_tetra(const ecs_coord_t coord[],ecs_int_t connect[])580 _orient_tetra(const ecs_coord_t coord[],
581 ecs_int_t connect[])
582 {
583 ecs_int_t isom_tmp;
584 ecs_coord_t det;
585 ecs_coord_t vect12[3];
586 ecs_coord_t vect13[3];
587 ecs_coord_t vect14[3];
588
589 ecs_int_t passage = -1;
590
591 #define ECS_LOC_INIT_VECT(vect, i, j) ( \
592 vect[0] = coord[((connect[j-1] - 1) * 3) ] \
593 - coord[((connect[i-1] - 1) * 3) ], \
594 vect[1] = coord[((connect[j-1] - 1) * 3) + 1] \
595 - coord[((connect[i-1] - 1) * 3) + 1], \
596 vect[2] = coord[((connect[j-1] - 1) * 3) + 2] \
597 - coord[((connect[i-1] - 1) * 3) + 2] )
598
599 #define ECS_LOC_PERMUTE(i, j) ( \
600 isom_tmp = connect[i-1], \
601 connect[i-1] = connect[j-1], \
602 connect[j-1] = isom_tmp )
603
604 /* Boucle sur les renumérotations possibles */
605
606 do {
607
608 passage += 1;
609
610 /* Initialisation */
611
612 switch(passage) {
613
614 case 0:
615 break;
616
617 case 1:
618 ECS_LOC_PERMUTE(2, 3);
619 break;
620
621 default: /* Retour connectivité d'origine et sortie */
622 ECS_LOC_PERMUTE(2, 3);
623 return -1;
624
625 }
626
627 ECS_LOC_INIT_VECT(vect12, 1, 2);
628 ECS_LOC_INIT_VECT(vect13, 1, 3);
629 ECS_LOC_INIT_VECT(vect14, 1, 4);
630
631 det = ECS_LOC_DETERMINANT(vect12, vect13, vect14);
632
633 } while (det < 0.0);
634
635 return (ECS_MIN(passage, 1));
636
637 #undef ECS_LOC_INIT_VECT
638 #undef ECS_LOC_PERMUTE
639 }
640
641 /*----------------------------------------------------------------------------
642 * Correction de l'orientation issue de Foam2VTK d'un tétraèdre
643 * en connectivité nodale ; renvoie 1.
644 *----------------------------------------------------------------------------*/
645
646 static ecs_int_t
_orient_tetra_of(ecs_int_t connect[])647 _orient_tetra_of(ecs_int_t connect[])
648 {
649 ecs_int_t isom_tmp;
650
651 #define ECS_LOC_PERMUTE(i, j) ( \
652 isom_tmp = connect[i-1], \
653 connect[i-1] = connect[j-1], \
654 connect[j-1] = isom_tmp )
655
656 ECS_LOC_PERMUTE(2, 3);
657
658 return 1;
659
660 #undef ECS_LOC_PERMUTE
661 }
662
663 /*----------------------------------------------------------------------------
664 * Correction si nécessaire de l'orientation d'une pyramide en connectivité
665 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
666 * ou que l'on demande une simple vérification (i.e. correc = false),
667 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
668 * l'orientation par une permutation de la connectivité locale.
669 *
670 * Tous les cas ne sont pas détectés ou traités : on suppose que le
671 * sommet est toujours en position 5, mais que la base peut être
672 * parcourue dans l'ordre 1 4 3 2, 1 2 4 3, ou 1 3 4 2 au lieu de 1 2 3 4,
673 * dans quel cas on la réordonne.
674 *----------------------------------------------------------------------------*/
675
676 static ecs_int_t
_orient_pyram(const ecs_coord_t coord[],ecs_int_t connect[],const bool correc)677 _orient_pyram(const ecs_coord_t coord[],
678 ecs_int_t connect[],
679 const bool correc)
680 {
681 ecs_int_t isom;
682 ecs_int_t isom_tmp;
683 ecs_int_t connect_tmp[8];
684 ecs_int_t ret_base;
685 ecs_coord_t det;
686 ecs_coord_t vect1[3];
687 ecs_coord_t vect2[3];
688 ecs_coord_t vect3[3];
689
690 ecs_int_t retval = 0;
691 ecs_int_t passage = -1;
692
693 #define ECS_LOC_INIT_VECT(vect, i, j) ( \
694 vect[0] = coord[((connect_tmp[j-1] - 1) * 3) ] \
695 - coord[((connect_tmp[i-1] - 1) * 3) ], \
696 vect[1] = coord[((connect_tmp[j-1] - 1) * 3) + 1] \
697 - coord[((connect_tmp[i-1] - 1) * 3) + 1], \
698 vect[2] = coord[((connect_tmp[j-1] - 1) * 3) + 2] \
699 - coord[((connect_tmp[i-1] - 1) * 3) + 2] )
700
701 #define ECS_LOC_PERMUTE(i, j) ( \
702 isom_tmp = connect[i-1], \
703 connect_tmp[i-1] = connect_tmp[j-1], \
704 connect_tmp[j-1] = isom_tmp )
705
706
707 for (isom = 0; isom < 5; isom++)
708 connect_tmp[isom] = connect[isom];
709
710 /* Vérification et correction éventuelle de l'orientation de la base */
711
712 ret_base = _orient_quad(coord,
713 connect_tmp,
714 correc);
715
716 retval = ECS_MAX(ret_base, retval);
717
718 if ((correc == false && ret_base != 0) || ret_base < 0)
719 return - 1;
720
721
722 /* Boucle sur les renumérotations possibles */
723
724 do {
725
726 passage += 1;
727
728 /* Initialisation */
729
730 switch(passage) {
731
732 case 0:
733 break;
734
735 case 1:
736 if (correc == false)
737 return -1;
738 else
739 retval = 1;
740 ECS_LOC_PERMUTE(2, 4);
741 break;
742
743 default: /* Retour connectivité d'origine et sortie */
744 ECS_LOC_PERMUTE(2, 4);
745 return -1;
746
747 }
748
749 ECS_LOC_INIT_VECT(vect1, 1, 2);
750 ECS_LOC_INIT_VECT(vect2, 1, 4);
751 ECS_LOC_INIT_VECT(vect3, 1, 5);
752
753 det = ECS_LOC_DETERMINANT(vect1, vect2, vect3);
754
755 } while (det < 0.0);
756
757 if (retval > 0) {
758 for (isom = 0; isom < 5; isom++)
759 connect[isom] = connect_tmp[isom];
760 }
761
762 return retval;
763
764 #undef ECS_LOC_INIT_VECT
765 #undef ECS_LOC_PERMUTE
766 }
767
768 /*----------------------------------------------------------------------------
769 * Correction de l'orientation issue de Foam2VTK d'une pyramide
770 * en connectivité nodale ; renvoie 1.
771 *----------------------------------------------------------------------------*/
772
773 static ecs_int_t
_orient_pyram_of(ecs_int_t connect[])774 _orient_pyram_of(ecs_int_t connect[])
775 {
776 ecs_int_t isom;
777 ecs_int_t isom_tmp;
778 ecs_int_t connect_tmp[8];
779
780 #define ECS_LOC_PERMUTE(i, j) ( \
781 isom_tmp = connect[i-1], \
782 connect_tmp[i-1] = connect_tmp[j-1], \
783 connect_tmp[j-1] = isom_tmp )
784
785 for (isom = 0; isom < 5; isom++)
786 connect_tmp[isom] = connect[isom];
787
788 /* Boucle sur les renumérotations possibles */
789
790 ECS_LOC_PERMUTE(2, 4);
791
792 for (isom = 0; isom < 5; isom++)
793 connect[isom] = connect_tmp[isom];
794
795 return 1;
796
797 #undef ECS_LOC_INIT_VECT
798 #undef ECS_LOC_PERMUTE
799 }
800
801 /*----------------------------------------------------------------------------
802 * Correction si nécessaire de l'orientation d'un prisme en connectivité
803 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
804 * ou que l'on demande une simple vérification (i.e. correc = false),
805 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
806 * l'orientation par une permutation de la connectivité locale.
807 *
808 * Tous les cas ne sont pas détectés ou traités : on suppose que les
809 * bases peuvent être parcourues dans l'ordre 1 3 2 et 4 6 5 au lieu
810 * de 1 2 3 et 4 5 6, dans quel cas on les réordonne.
811 *----------------------------------------------------------------------------*/
812
813 static ecs_int_t
_orient_prism(const ecs_coord_t coord[],ecs_int_t connect[],const bool correc)814 _orient_prism(const ecs_coord_t coord[],
815 ecs_int_t connect[],
816 const bool correc)
817 {
818 ecs_int_t idim;
819 ecs_int_t isom_tmp;
820 ecs_coord_t pscal;
821 ecs_coord_t vect1[3];
822 ecs_coord_t vect2[3];
823 ecs_coord_t vect3[3];
824
825 ecs_int_t passage = -1;
826
827 #define ECS_LOC_INIT_VECT(vect, i, j) ( \
828 vect[0] = coord[((connect[j-1] - 1) * 3) ] \
829 - coord[((connect[i-1] - 1) * 3) ], \
830 vect[1] = coord[((connect[j-1] - 1) * 3) + 1] \
831 - coord[((connect[i-1] - 1) * 3) + 1], \
832 vect[2] = coord[((connect[j-1] - 1) * 3) + 2] \
833 - coord[((connect[i-1] - 1) * 3) + 2] )
834
835 #define ECS_LOC_PERMUTE(i, j) ( \
836 isom_tmp = connect[i-1], \
837 connect[i-1] = connect[j-1], \
838 connect[j-1] = isom_tmp )
839
840 /* Boucle sur les renumérotations possibles */
841
842 do {
843
844 passage += 1;
845
846 /* Initialisation */
847
848 switch(passage) {
849
850 case 0:
851 break;
852
853 case 1:
854 if (correc == false)
855 return -1;
856 ECS_LOC_PERMUTE(2, 3);
857 ECS_LOC_PERMUTE(5, 6);
858 break;
859
860 default: /* Retour connectivité d'origine et sortie */
861 ECS_LOC_PERMUTE(2, 3);
862 ECS_LOC_PERMUTE(5, 6);
863 return -1;
864
865 }
866
867 ECS_LOC_INIT_VECT(vect2, 1, 2);
868 ECS_LOC_INIT_VECT(vect3, 1, 3);
869
870 ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3);
871
872 for (idim = 0; idim < 3; idim++)
873 vect2[idim] = ( coord[((connect[4-1] - 1) * 3) + idim]
874 + coord[((connect[5-1] - 1) * 3) + idim]
875 + coord[((connect[6-1] - 1) * 3) + idim]
876 - coord[((connect[1-1] - 1) * 3) + idim]
877 - coord[((connect[2-1] - 1) * 3) + idim]
878 - coord[((connect[3-1] - 1) * 3) + idim]);
879
880 pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2);
881
882 } while (pscal < 0.0);
883
884 return (ECS_MIN(passage, 1));
885
886 #undef ECS_LOC_INIT_VECT
887 #undef ECS_LOC_PERMUTE
888 }
889
890 /*----------------------------------------------------------------------------
891 * Correction si nécessaire de l'orientation d'un hexaèdre en connectivité
892 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
893 * ou que l'on demande une simple vérification (i.e. correc = false),
894 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
895 * l'orientation par une permutation de la connectivité locale.
896 *
897 * Tous les cas ne sont pas détectés ou traités : on suppose que les
898 * bases peuvent être parcourues dans l'ordre soit 1 4 3 2 et 5 8 7 6,
899 * soit 1 2 4 3 et 5 6 8 7, soit 1 3 4 2 et 5 7 8 6 au lieu
900 * de 1 2 3 4 et 5 6 7 8, dans quel cas on les réordonne.
901 *----------------------------------------------------------------------------*/
902
903 static ecs_int_t
_orient_hexa(const ecs_coord_t coord[],ecs_int_t connect[],const bool correc)904 _orient_hexa(const ecs_coord_t coord[],
905 ecs_int_t connect[],
906 const bool correc)
907 {
908 ecs_int_t idim;
909 ecs_int_t isom;
910 ecs_int_t isom_tmp;
911 ecs_int_t connect_tmp[8];
912 ecs_int_t ret_base_1;
913 ecs_int_t ret_base_2;
914 ecs_coord_t pscal;
915 ecs_coord_t vect1[3];
916 ecs_coord_t vect2[3];
917 ecs_coord_t vect3[3];
918
919 ecs_int_t retval = 0;
920 ecs_int_t passage = -1;
921
922 #define ECS_LOC_INIT_VECT(vect, i, j) ( \
923 vect[0] = coord[((connect_tmp[j-1] - 1) * 3) ] \
924 - coord[((connect_tmp[i-1] - 1) * 3) ], \
925 vect[1] = coord[((connect_tmp[j-1] - 1) * 3) + 1] \
926 - coord[((connect_tmp[i-1] - 1) * 3) + 1], \
927 vect[2] = coord[((connect_tmp[j-1] - 1) * 3) + 2] \
928 - coord[((connect_tmp[i-1] - 1) * 3) + 2] )
929
930 #define ECS_LOC_PERMUTE(i, j) ( \
931 isom_tmp = connect[i-1], \
932 connect_tmp[i-1] = connect_tmp[j-1], \
933 connect_tmp[j-1] = isom_tmp )
934
935
936 for (isom = 0; isom < 8; isom++)
937 connect_tmp[isom] = connect[isom];
938
939 /* Vérification et correction éventuelle de l'orientation des bases */
940
941 ret_base_1 = _orient_quad(coord,
942 connect_tmp,
943 correc);
944
945 if ((correc == false && ret_base_1 != 0) || ret_base_1 < 0)
946 return - 1;
947
948 else if (ret_base_1 > 0) {
949 ret_base_2 = _orient_quad(coord,
950 connect_tmp + 4,
951 correc);
952 if (ret_base_2 != ret_base_1)
953 return - 1;
954 else
955 retval = 1;
956 }
957
958 /* Boucle sur les renumérotations possibles */
959
960 do {
961
962 passage += 1;
963
964 /* Initialisation */
965
966 switch(passage) {
967
968 case 0:
969 break;
970
971 case 1:
972 if (correc == false)
973 return -1;
974 else
975 retval = 1;
976 ECS_LOC_PERMUTE(2, 4);
977 ECS_LOC_PERMUTE(6, 8);
978 break;
979
980 default:
981 return -1;
982 break;
983
984 }
985
986 ECS_LOC_INIT_VECT(vect2, 1, 2);
987 ECS_LOC_INIT_VECT(vect3, 1, 4);
988
989 ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3);
990
991 for (idim = 0; idim < 3; idim++)
992 vect2[idim] = ( coord[((connect_tmp[5-1] - 1) * 3) + idim]
993 + coord[((connect_tmp[6-1] - 1) * 3) + idim]
994 + coord[((connect_tmp[7-1] - 1) * 3) + idim]
995 + coord[((connect_tmp[8-1] - 1) * 3) + idim]
996 - coord[((connect_tmp[1-1] - 1) * 3) + idim]
997 - coord[((connect_tmp[2-1] - 1) * 3) + idim]
998 - coord[((connect_tmp[3-1] - 1) * 3) + idim]
999 - coord[((connect_tmp[4-1] - 1) * 3) + idim]) * 0.25;
1000
1001 pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2);
1002
1003 } while (pscal < 0.0);
1004
1005 if (retval > 0) {
1006 for (isom = 0; isom < 8; isom++)
1007 connect[isom] = connect_tmp[isom];
1008 }
1009
1010 /* Vérification et correction éventuelle de l'orientation des cotés */
1011
1012 connect_tmp[1-1] = connect[2-1];
1013 connect_tmp[2-1] = connect[3-1];
1014 connect_tmp[3-1] = connect[7-1];
1015 connect_tmp[4-1] = connect[6-1];
1016
1017 ret_base_1 = _orient_quad(coord,
1018 connect_tmp,
1019 correc);
1020
1021 if (ret_base_1 != 0)
1022 return - 1;
1023
1024 connect_tmp[1-1] = connect[1-1];
1025 connect_tmp[2-1] = connect[2-1];
1026 connect_tmp[3-1] = connect[6-1];
1027 connect_tmp[4-1] = connect[5-1];
1028
1029 ret_base_1 = _orient_quad(coord,
1030 connect_tmp,
1031 correc);
1032
1033 if (ret_base_1 != 0)
1034 return - 1;
1035
1036 return retval;
1037
1038 #undef ECS_LOC_INIT_VECT
1039 #undef ECS_LOC_PERMUTE
1040 }
1041
1042 /*----------------------------------------------------------------------------
1043 * Correction de l'orientation issue de Foam2VTK d'un hexaedre
1044 * en connectivité nodale ; renvoie 1.
1045 *----------------------------------------------------------------------------*/
1046
1047 static ecs_int_t
_orient_hexa_of(ecs_int_t connect[])1048 _orient_hexa_of(ecs_int_t connect[])
1049 {
1050 ecs_int_t isom;
1051 ecs_int_t isom_tmp;
1052 ecs_int_t connect_tmp[8];
1053
1054 #define ECS_LOC_PERMUTE(i, j) ( \
1055 isom_tmp = connect[i-1], \
1056 connect_tmp[i-1] = connect_tmp[j-1], \
1057 connect_tmp[j-1] = isom_tmp )
1058
1059 for (isom = 0; isom < 8; isom++)
1060 connect_tmp[isom] = connect[isom];
1061
1062 /* Renumbering */
1063
1064 ECS_LOC_PERMUTE(2, 4);
1065 ECS_LOC_PERMUTE(6, 8);
1066
1067 for (isom = 0; isom < 8; isom++)
1068 connect[isom] = connect_tmp[isom];
1069
1070 return 1;
1071
1072 #undef ECS_LOC_PERMUTE
1073 }
1074
1075 /*----------------------------------------------------------------------------
1076 * Compute contribution of a polygonal face to a cell's volume,
1077 * using Stoke's theorem
1078 *
1079 * Pi+1
1080 * *---------* B : barycentre of the polygon
1081 * / . . \
1082 * / . . \ Pi : vertices of the polygon
1083 * / . . \
1084 * / . . Ti \ Ti : triangle
1085 * *.........B.........* Pi
1086 * Pn-1 \ . . /
1087 * \ . . /
1088 * \ . . /
1089 * \ . T0 . /
1090 * *---------*
1091 * P0
1092 *----------------------------------------------------------------------------*/
1093
1094 static double
_compute_face_vol_contrib(ecs_int_t n_face_vertices,const ecs_coord_t vtx_coord[],const ecs_int_t face_vtx_lst[])1095 _compute_face_vol_contrib(ecs_int_t n_face_vertices,
1096 const ecs_coord_t vtx_coord[],
1097 const ecs_int_t face_vtx_lst[])
1098 {
1099 ecs_int_t i, tri_id, vtx_id;
1100 ecs_coord_t tri_center[3], tri_norm[3];
1101 ecs_coord_t vect1[3], vect2[3];
1102
1103 ecs_coord_t face_barycentre[3] = {0., 0., 0.};
1104
1105 double inv3 = 1./3.;
1106
1107 double vol_contrib = 0.;
1108
1109 /* Compute barycentre (B) coordinates for the polygon (P) */
1110
1111 for (vtx_id = 0; vtx_id < n_face_vertices; vtx_id++) {
1112 size_t coord_id = face_vtx_lst[vtx_id] - 1;
1113 for (i = 0; i < 3; i++)
1114 face_barycentre[i] += vtx_coord[coord_id*3 + i];
1115 }
1116
1117 for (i = 0; i < 3; i++)
1118 face_barycentre[i] /= n_face_vertices;
1119
1120 /* Loop on triangles of the face */
1121
1122 for (tri_id = 0; tri_id < n_face_vertices; tri_id++) {
1123
1124 size_t coord_id_1 = face_vtx_lst[tri_id % n_face_vertices] - 1;
1125 size_t coord_id_2 = face_vtx_lst[(tri_id + 1)% n_face_vertices] - 1;
1126
1127 for (i = 0; i < 3; i++) {
1128 vect1[i] = vtx_coord[coord_id_1*3 + i] - face_barycentre[i];
1129 vect2[i] = vtx_coord[coord_id_2*3 + i] - face_barycentre[i];
1130 }
1131
1132 ECS_LOC_PRODUIT_VECTORIEL(tri_norm, vect1, vect2);
1133
1134 for (i = 0; i < 3; i++)
1135 tri_norm[i] *= 0.5;
1136
1137 /* Computation of the center of the triangle Ti */
1138
1139 for (i = 0; i < 3; i++)
1140 tri_center[i] = ( face_barycentre[i]
1141 + vtx_coord[coord_id_1*3 + i]
1142 + vtx_coord[coord_id_2*3 + i]) * inv3;
1143
1144 /* Contribution to cell volume using Stoke's formula */
1145
1146 for (i = 0; i < 3; i++)
1147 vol_contrib += (tri_norm[i] * tri_center[i]);
1148
1149 } /* End of loop on triangles of the face */
1150
1151 return vol_contrib;
1152 }
1153
1154 /*----------------------------------------------------------------------------
1155 * Correction si nécessaire de l'orientation d'un polyedre en connectivité
1156 * nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
1157 * ou que l'on demande une simple vérification (i.e. correc = false),
1158 * 0 si l'orientation initiale est bonne, et 1 en cas de correction de
1159 * l'orientation par une permutation de la connectivité locale.
1160 *
1161 * On suppose que toutes les faces polygonales sont bien définies, mais
1162 * que certaines faces peuvent être définies avec une normale intérieure
1163 * à la cellule et non pas extérieure. On suppose que la non-convexité
1164 * du polyèdre est limitée.
1165 *----------------------------------------------------------------------------*/
1166
1167 static ecs_int_t
_orient_polyhedron(const ecs_coord_t coord[],ecs_int_t connect[],ecs_int_t size,bool correc,ecs_tab_int_t * face_index,ecs_tab_int_t * face_marker,ecs_tab_int_t * edges)1168 _orient_polyhedron(const ecs_coord_t coord[],
1169 ecs_int_t connect[],
1170 ecs_int_t size,
1171 bool correc,
1172 ecs_tab_int_t *face_index,
1173 ecs_tab_int_t *face_marker,
1174 ecs_tab_int_t *edges)
1175 {
1176 size_t i, j, face_id, edge_id;
1177
1178 double cell_vol = 0.;
1179 size_t n_unmarked_faces = 0;
1180 size_t n_faces = 0;
1181 size_t n_edges = 0;
1182 ecs_int_t retval = 0;
1183
1184 /* Ensure working arrays are large enough */
1185
1186 /* Mark faces */
1187
1188 {
1189 ecs_int_t premier_som = -1;
1190
1191 for (i = 0; i < (size_t)size; i++) {
1192
1193 if (face_index->nbr < n_faces + 1) {
1194 face_index->nbr = ECS_MAX(n_faces + 1, face_index->nbr*2);
1195 ECS_REALLOC(face_index->val, face_index->nbr, ecs_int_t);
1196 }
1197
1198 if (premier_som == -1) {
1199 face_index->val[n_faces] = i;
1200 premier_som = connect[i];
1201 }
1202 else if (connect[i] == premier_som) {
1203 n_faces += 1;
1204 premier_som = -1;
1205 }
1206 }
1207
1208 if (face_index->nbr < n_faces + 1) {
1209 face_index->nbr = ECS_MAX(n_faces + 1, face_index->nbr*2);
1210 ECS_REALLOC(face_index->val, face_index->nbr, ecs_int_t);
1211 }
1212
1213 face_index->val[n_faces] = size;
1214
1215 /* face_marker: 0 initially, 1 if oriented, -1 if inverted */
1216
1217 if (face_marker->nbr < n_faces) {
1218 face_marker->nbr = ECS_MAX(n_faces, face_marker->nbr*2);
1219 ECS_REALLOC(face_marker->val, face_marker->nbr, ecs_int_t);
1220 }
1221
1222 for (face_id = 0; face_id < n_faces; face_id++)
1223 face_marker->val[face_id] = 0;
1224 }
1225
1226 if (n_faces == 0)
1227 return 0;
1228
1229 /* Extract edges and build edge-face connectivity */
1230 /*------------------------------------------------*/
1231
1232 for (face_id = 0; face_id < n_faces; face_id++) {
1233
1234 size_t start_id = face_index->val[face_id];
1235 size_t end_id = face_index->val[face_id + 1] - 1;
1236
1237 /* Loop on unassociated edges */
1238
1239 for (i = start_id; i < end_id; i++) {
1240
1241 ecs_int_t v_id_1 = connect[i], v_id_2 = connect[i + 1];
1242 ecs_int_t is_match = 0;
1243
1244 /* Compare with other edges */
1245
1246 for (edge_id = 0; edge_id < n_edges; edge_id++) {
1247
1248 ecs_int_t v_id_1_cmp = edges->val[edge_id*4];
1249 ecs_int_t v_id_2_cmp = edges->val[edge_id*4 + 1];
1250
1251 if ((v_id_1 == v_id_2_cmp) && (v_id_2 == v_id_1_cmp))
1252 is_match = 1;
1253
1254 else if ((v_id_1 == v_id_1_cmp) && (v_id_2 == v_id_2_cmp))
1255 is_match = -1;
1256
1257 /* Each edge should be shared by exactly 2 faces */
1258
1259 if (is_match != 0) {
1260 if (edges->val[edge_id*4 + 3] == 0) {
1261 edges->val[edge_id*4 + 3] = is_match * (face_id + 1);
1262 break;
1263 }
1264 else
1265 return -1;
1266 }
1267 }
1268
1269 /* If edge is unmatched, add it */
1270
1271 if (is_match == 0) {
1272
1273 if (edges->nbr < (n_edges + 1)*4) {
1274 edges->nbr = ECS_MAX(edges->nbr*2, (n_edges + 1)*4);
1275 ECS_REALLOC(edges->val, edges->nbr, ecs_int_t);
1276 }
1277
1278 edges->val[n_edges*4] = v_id_1;
1279 edges->val[n_edges*4 + 1] = v_id_2;
1280 edges->val[n_edges*4 + 2] = face_id + 1;
1281 edges->val[n_edges*4 + 3] = 0;
1282
1283 n_edges += 1;
1284 }
1285
1286 }
1287
1288 }
1289
1290 /* Check if each edge is associated with 2 faces
1291 (i.e., that cell is closed) */
1292
1293 for (edge_id = 0; edge_id < n_edges; edge_id++) {
1294 if (edges->val[edge_id*4 + 3] == 0)
1295 return -1;
1296 }
1297
1298 /* Now check if face orientations are compatible */
1299 /*-----------------------------------------------*/
1300
1301 face_marker->val[0] = 1; /* First face defines reference */
1302 n_unmarked_faces = n_faces - 1;
1303
1304 while (n_unmarked_faces > 0) {
1305
1306 for (edge_id = 0; edge_id < n_edges; edge_id++) {
1307
1308 ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
1309 ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
1310 ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;
1311
1312 if ( face_marker->val[face_id_1] == 0
1313 && face_marker->val[face_id_2] != 0) {
1314 if (face_num_2 > 0)
1315 face_marker->val[face_id_1] = face_marker->val[face_id_2];
1316 else
1317 face_marker->val[face_id_1] = - face_marker->val[face_id_2];
1318 n_unmarked_faces -= 1;
1319 }
1320
1321 else if ( face_marker->val[face_id_1] != 0
1322 && face_marker->val[face_id_2] == 0) {
1323 if (face_num_2 > 0)
1324 face_marker->val[face_id_2] = face_marker->val[face_id_1];
1325 else
1326 face_marker->val[face_id_2] = - face_marker->val[face_id_1];
1327 n_unmarked_faces -= 1;
1328 }
1329 }
1330
1331 }
1332
1333 for (edge_id = 0; edge_id < n_edges; edge_id++) {
1334
1335 ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
1336 ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
1337 ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;
1338
1339 ecs_int_t marker_product
1340 = face_marker->val[face_id_1]*face_marker->val[face_id_2];
1341
1342 if ( (face_num_2 < 0 && marker_product > 0)
1343 || (face_num_2 > 0 && marker_product < 0))
1344 return -1;
1345
1346 }
1347
1348 /* At this stage, topology is correct; check for inside-out cell */
1349 /*---------------------------------------------------------------*/
1350
1351 for (face_id = 0; face_id < n_faces; face_id++) {
1352
1353 size_t start_id = face_index->val[face_id];
1354 size_t end_id = face_index->val[face_id + 1] - 1;
1355
1356 double vol_contrib
1357 = _compute_face_vol_contrib(end_id - start_id,
1358 coord,
1359 connect + start_id);
1360 cell_vol += vol_contrib * (double)(face_marker->val[face_id]);
1361 }
1362
1363 /* Invert orientation if cell is inside_out */
1364
1365 if (cell_vol < 0.) {
1366 for (face_id = 0; face_id < n_faces; face_id++)
1367 face_marker->val[face_id] *= -1;
1368 }
1369
1370 /* Now correct connectivity if required */
1371 /*--------------------------------------*/
1372
1373 if (correc == true) {
1374
1375 for (face_id = 0; face_id < n_faces; face_id++) {
1376
1377 if (face_marker->val[face_id] == -1) {
1378
1379 for (i = face_index->val[face_id] + 1,
1380 j = face_index->val[face_id + 1] - 2;
1381 i < j;
1382 i++, j--) {
1383 ecs_int_t connect_tmp = connect[i];
1384 connect[i] = connect[j];
1385 connect[j] = connect_tmp;
1386 }
1387
1388 retval = 1;
1389
1390 }
1391 }
1392 }
1393
1394 else {
1395
1396 for (face_id = 0; face_id < n_faces; face_id++) {
1397 if (face_marker->val[face_id] == -1)
1398 retval = -1;
1399 }
1400 }
1401
1402 return retval;
1403 }
1404
1405 /*============================================================================
1406 * Fonctions publiques
1407 *============================================================================*/
1408
1409 /*----------------------------------------------------------------------------
1410 * Fonction qui réalise le tri des types géométriques
1411 * La fonction affiche le nombre d'éléments par type géométrique
1412 *----------------------------------------------------------------------------*/
1413
1414 ecs_tab_int_t
ecs_table_def__trie_typ(ecs_table_t * this_table_def,int dim_elt)1415 ecs_table_def__trie_typ(ecs_table_t *this_table_def,
1416 int dim_elt)
1417 {
1418 size_t ielt;
1419 size_t nbr_elt;
1420 size_t nbr_som;
1421
1422 ecs_tab_int_t tab_typ_geo_ord;
1423 ecs_tab_int_t tab_typ_geo;
1424 ecs_tab_int_t vect_renum;
1425
1426 ecs_elt_typ_t typ_geo;
1427
1428 const ecs_elt_typ_t *typ_geo_base;
1429
1430 const int lng_imp = ECS_LNG_AFF_STR - ECS_LOC_LNG_MAX_NOM_TYP;
1431
1432 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1433
1434 assert(this_table_def != NULL);
1435
1436 vect_renum.nbr = 0 ;
1437 vect_renum.val = NULL;
1438
1439 nbr_elt = this_table_def->nbr;
1440
1441 if (nbr_elt == 0)
1442 return vect_renum;
1443
1444 /* Détermination de base du type d'élément "classique" */
1445
1446 assert(dim_elt >=2 && dim_elt <= 3);
1447
1448 typ_geo_base = ecs_glob_typ_elt[dim_elt - 2];
1449
1450 /* Si tous les éléments sont de même type, rien à faire */
1451 /*------------------------------------------------------*/
1452
1453 if (this_table_def->pos == NULL) {
1454
1455 if (this_table_def->pas < 9)
1456 typ_geo = typ_geo_base[this_table_def->pas];
1457 else if (dim_elt == 2)
1458 typ_geo = ECS_ELT_TYP_FAC_POLY;
1459 else if (dim_elt == 3)
1460 typ_geo = ECS_ELT_TYP_CEL_POLY;
1461 else
1462 typ_geo = ECS_ELT_TYP_NUL;
1463
1464 printf(" %-*s%*.*s : %*lu\n",
1465 lng_imp, _("Number of elements"),
1466 ECS_LOC_LNG_MAX_NOM_TYP, ECS_LOC_LNG_MAX_NOM_TYP,
1467 ecs_fic_elt_typ_liste_c[typ_geo].nom,
1468 ECS_LNG_AFF_ENT, (unsigned long)nbr_elt);
1469
1470 }
1471
1472 /* Si tous les éléments ne sont pas de même type */
1473 /*-----------------------------------------------*/
1474
1475 else {
1476
1477 /* Construction d'un tableau temporaire de type géométrique */
1478
1479 tab_typ_geo.nbr = nbr_elt;
1480 ECS_MALLOC(tab_typ_geo.val, tab_typ_geo.nbr, ecs_int_t);
1481
1482 for (ielt = 0; ielt < nbr_elt; ielt++) {
1483
1484 nbr_som = this_table_def->pos[ielt + 1]
1485 - this_table_def->pos[ielt];
1486
1487 if (nbr_som < 9)
1488 tab_typ_geo.val[ielt] = typ_geo_base[nbr_som];
1489
1490 else if (dim_elt == 2)
1491 tab_typ_geo.val[ielt] = ECS_ELT_TYP_FAC_POLY;
1492
1493 else if (dim_elt == 3)
1494 tab_typ_geo.val[ielt] = ECS_ELT_TYP_CEL_POLY;
1495
1496 else
1497 tab_typ_geo.val[ielt] = ECS_ELT_TYP_NUL;
1498
1499 }
1500
1501 /* On regarde si les types géométriques ne sont pas deja ordonnés */
1502
1503 ielt = 1;
1504 while (ielt < nbr_elt &&
1505 tab_typ_geo.val[ielt] >= tab_typ_geo.val[ielt - 1] )
1506 ielt++;
1507
1508 if (ielt < nbr_elt) {
1509
1510 /* Les types géométriques ne sont pas ordonnés */
1511 /* On ordonne les types géométriques */
1512
1513 vect_renum.nbr = nbr_elt;
1514 ECS_MALLOC(vect_renum.val, nbr_elt, ecs_int_t);
1515
1516
1517 tab_typ_geo_ord = ecs_tab_int__trie_et_renvoie(tab_typ_geo,
1518 vect_renum);
1519
1520 /*
1521 `vect_renum' prend pour indice les indices nouveaux, et ses valeurs
1522 contiennent les indices anciens correspondants
1523 On inverse le contenu de `vect_renum' :
1524 à chaque indice ancien, `vect_renum' donne la valeur du nouvel indice
1525 */
1526
1527 ecs_tab_int__inverse(&vect_renum);
1528
1529 ECS_FREE(tab_typ_geo.val);
1530
1531 tab_typ_geo.val = tab_typ_geo_ord.val;
1532
1533 tab_typ_geo_ord.nbr = 0;
1534 tab_typ_geo_ord.val = NULL;
1535
1536 }
1537
1538 /* Message d'information sur la composition des éléments du maillage */
1539
1540 {
1541 ecs_int_t val_typ_ref;
1542 size_t nbr_val_typ_geo;
1543
1544 size_t cpt_ielt = 0;
1545
1546 while (cpt_ielt < nbr_elt) {
1547
1548 val_typ_ref = tab_typ_geo.val[cpt_ielt];
1549
1550 /* On compte le nombre d'éléments ayant le même type géométrique */
1551
1552 nbr_val_typ_geo = 0;
1553
1554 for (ielt = cpt_ielt;
1555 ielt < nbr_elt && tab_typ_geo.val[ielt] == val_typ_ref; ielt++)
1556 nbr_val_typ_geo++;
1557
1558 printf(" %-*s%*.*s : %*lu\n",
1559 lng_imp, _("Number of elements"),
1560 ECS_LOC_LNG_MAX_NOM_TYP, ECS_LOC_LNG_MAX_NOM_TYP,
1561 ecs_fic_elt_typ_liste_c[val_typ_ref].nom,
1562 ECS_LNG_AFF_ENT, (unsigned long)nbr_val_typ_geo);
1563
1564 cpt_ielt += nbr_val_typ_geo;
1565
1566 }
1567 }
1568
1569 /* Le tableau de type géométrique n'est plus nécessaire */
1570
1571 tab_typ_geo.nbr = 0;
1572 ECS_FREE(tab_typ_geo.val);
1573
1574 }
1575
1576 /* Renvoi du vecteur de renumérotation */
1577 /*-------------------------------------*/
1578
1579 return vect_renum;
1580 }
1581
1582 /*----------------------------------------------------------------------------
1583 * Fonction qui construit
1584 * les définitions des faces par décomposition des tables des cellules
1585 *----------------------------------------------------------------------------*/
1586
1587 void
ecs_table_def__decompose_cel(ecs_table_t ** table_def_fac,ecs_table_t * table_def_cel)1588 ecs_table_def__decompose_cel(ecs_table_t **table_def_fac,
1589 ecs_table_t *table_def_cel)
1590 {
1591 size_t nbr_cel;
1592 size_t nbr_def;
1593 size_t nbr_fac;
1594 size_t nbr_fac_old;
1595 size_t nbr_val_cel;
1596 size_t nbr_val_fac;
1597 size_t nbr_val_fac_old;
1598 size_t ind_pos_cel;
1599 size_t ind_pos_loc;
1600 size_t ind_pos_sui;
1601 size_t nbr_pos_loc;
1602 ecs_int_t num_def;
1603 ecs_int_t marqueur_fin;
1604
1605 size_t isom;
1606 size_t icel;
1607 int idef;
1608 size_t cpt_fac;
1609 int ifac;
1610
1611 ecs_int_t typ_geo_cel; /* Type géométrique élément en cours */
1612
1613 ecs_int_t typ_geo_base[9] = {ECS_ELT_TYP_NUL,
1614 ECS_ELT_TYP_NUL,
1615 ECS_ELT_TYP_NUL,
1616 ECS_ELT_TYP_NUL,
1617 ECS_ELT_TYP_CEL_TETRA,
1618 ECS_ELT_TYP_CEL_PYRAM,
1619 ECS_ELT_TYP_CEL_PRISM,
1620 ECS_ELT_TYP_NUL,
1621 ECS_ELT_TYP_CEL_HEXA};
1622
1623 ecs_size_t *def_cel_fac_pos = NULL;
1624 ecs_int_t *def_cel_fac_val = NULL;
1625
1626 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1627
1628 /*=================*/
1629 /* Initialisations */
1630 /*=================*/
1631
1632 nbr_cel = table_def_cel->nbr;
1633
1634 ecs_table__regle_en_pos(table_def_cel);
1635
1636 /* Construction, pour les cellules, des tables principaux */
1637 /*--------------------------------------------------------*/
1638
1639 /* Boucle de comptage pour l'allocation des sous-éléments */
1640 /*--------------------------------------------------------*/
1641
1642 if (*table_def_fac != NULL) {
1643 nbr_fac_old = (*table_def_fac)->nbr;
1644 nbr_val_fac_old = ecs_table__ret_val_nbr(*table_def_fac);
1645 }
1646 else {
1647 nbr_fac_old = 0;
1648 nbr_val_fac_old = 0;
1649 }
1650
1651 nbr_fac = 0;
1652 nbr_val_fac = 0;
1653 nbr_def = 0;
1654
1655 for (icel = 0; icel < nbr_cel; icel++ ) {
1656
1657 ind_pos_cel = table_def_cel->pos[icel] - 1;
1658
1659 nbr_val_cel = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1660
1661 /* Traitement des cellules "classiques" */
1662 /*--------------------------------------*/
1663
1664 if (nbr_val_cel < 9) {
1665
1666 typ_geo_cel = typ_geo_base[nbr_val_cel];
1667
1668 /* Boucle sur les sous-éléments définissant la cellulle */
1669
1670 for (ifac = 0;
1671 ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
1672 ifac++ ) {
1673
1674 const ecs_sous_elt_t * sous_elt
1675 = &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]);
1676
1677 nbr_fac++;
1678
1679 for (idef = 0;
1680 idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som;
1681 idef++);
1682
1683 nbr_val_fac += idef;
1684
1685 }
1686
1687 }
1688
1689 /* Traitement des éléments de type "polyèdre" */
1690 /*--------------------------------------------*/
1691
1692 else {
1693
1694 ind_pos_sui = table_def_cel->pos[icel + 1] - 1;
1695 nbr_pos_loc = ind_pos_sui - ind_pos_cel;
1696
1697 /* Convention : définition nodale cellule->sommets avec numéros de
1698 premiers sommets répétés en fin de liste pour marquer la fin
1699 de chaque face */
1700
1701 marqueur_fin = -1;
1702
1703 for (isom = ind_pos_cel; isom < ind_pos_sui; isom++) {
1704
1705 if (table_def_cel->val[isom] != marqueur_fin) {
1706 nbr_val_fac += 1;
1707 if (marqueur_fin == -1)
1708 marqueur_fin = table_def_cel->val[isom];
1709 }
1710 else {
1711 marqueur_fin = -1;
1712 nbr_fac += 1;
1713 }
1714
1715 }
1716
1717 }
1718
1719 } /* Fin de la boucle de comptage sur les cellules */
1720
1721 /* Allocation et initialisation pour les faces */
1722 /* des tableaux associés aux définitions */
1723 /*----------------------------------------------*/
1724
1725 nbr_fac += nbr_fac_old;
1726 nbr_val_fac += nbr_val_fac_old;
1727
1728 if (*table_def_fac != NULL) {
1729 ecs_table__regle_en_pos(*table_def_fac);
1730 (*table_def_fac)->nbr = nbr_fac;
1731 ECS_REALLOC((*table_def_fac)->pos, nbr_fac + 1, ecs_size_t);
1732 ECS_REALLOC((*table_def_fac)->val, nbr_val_fac, ecs_int_t);
1733 }
1734 else {
1735 *table_def_fac = ecs_table__alloue(nbr_fac,
1736 nbr_val_fac);
1737 (*table_def_fac)->pos[0] = 1;
1738 }
1739
1740 ECS_MALLOC(def_cel_fac_pos, nbr_cel + 1, ecs_size_t);
1741 ECS_MALLOC(def_cel_fac_val, nbr_fac - nbr_fac_old, ecs_int_t);
1742
1743 def_cel_fac_pos[0] = 1;
1744
1745 /*=======================================*/
1746 /* Boucle sur les cellules à transformer */
1747 /*=======================================*/
1748
1749 cpt_fac = nbr_fac_old;
1750 nbr_def = nbr_val_fac_old;
1751
1752 for (icel = 0; icel < nbr_cel; icel++ ) {
1753
1754 ind_pos_cel = table_def_cel->pos[icel] - 1;
1755
1756 nbr_val_cel = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1757
1758 /*--------------------------------------*/
1759 /* Traitement des éléments "classiques" */
1760 /*--------------------------------------*/
1761
1762 if (nbr_val_cel < 9) {
1763
1764 typ_geo_cel = typ_geo_base[nbr_val_cel];
1765
1766 /* Boucle sur les faces définissant la cellulle */
1767 /*==============================================*/
1768
1769 for (ifac = 0;
1770 ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
1771 ifac++ ) {
1772
1773 const ecs_sous_elt_t * sous_elt
1774 = &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]);
1775
1776 /* Boucle sur les sommets définissant la face */
1777 /*--------------------------------------------*/
1778
1779 for (idef = 0;
1780 idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som;
1781 idef++) {
1782
1783 /* Définition de la face en fonction des sommets */
1784
1785 num_def = sous_elt->som[idef];
1786
1787 (*table_def_fac)->val[nbr_def++]
1788 = table_def_cel->val[ind_pos_cel + num_def - 1];
1789
1790 }
1791
1792 /* Position de la face dans sa définition en fonction des sommets */
1793 /*----------------------------------------------------------------*/
1794
1795 (*table_def_fac)->pos[cpt_fac + 1]
1796 = (*table_def_fac)->pos[cpt_fac] + idef;
1797
1798 /* Détermination de la cellule en fonction des faces */
1799 /*---------------------------------------------------*/
1800
1801 def_cel_fac_val[cpt_fac - nbr_fac_old] = cpt_fac + 1;
1802
1803 cpt_fac++;
1804
1805 } /* Fin de la boucle sur les faces d'une cellule */
1806
1807 }
1808
1809
1810 /*--------------------------------------------*/
1811 /* Traitement des éléments de type "polyèdre" */
1812 /*--------------------------------------------*/
1813
1814 else {
1815
1816 /* Boucle sur les faces définissant le polyèdre */
1817 /*==============================================*/
1818
1819 ifac = 0;
1820 marqueur_fin = -1;
1821
1822 nbr_pos_loc = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1823
1824 for (ind_pos_loc = 0; ind_pos_loc < nbr_pos_loc; ind_pos_loc++) {
1825
1826 /* Définition de la face en fonction des sommets */
1827
1828 if (table_def_cel->val[ind_pos_cel + ind_pos_loc] != marqueur_fin) {
1829
1830 (*table_def_fac)->val[nbr_def++]
1831 = table_def_cel->val[ind_pos_cel + ind_pos_loc];
1832
1833 if (marqueur_fin == -1)
1834 marqueur_fin = table_def_cel->val[ind_pos_cel + ind_pos_loc];
1835
1836 }
1837
1838 /* Position de la face dans sa définition en fonction des sommets */
1839
1840 else {
1841
1842 (*table_def_fac)->pos[cpt_fac + 1] = nbr_def + 1;
1843
1844 marqueur_fin = -1;
1845
1846 def_cel_fac_val[cpt_fac - nbr_fac_old] = cpt_fac + 1;
1847
1848 /* Incrémentation du nombre de faces */
1849
1850 ifac++;
1851 cpt_fac++;
1852
1853 }
1854
1855 } /* Fin de la boucle sur les faces du polyèdre */
1856
1857 }
1858
1859 /* A ce point, on a : ifac
1860 = ecs_fic_elt_typ_liste_c[table_typ_geo_cel->val[icel]].nbr_sous_elt
1861 pour des éléments classiques, et ifac est égal au nombre de faces pour
1862 des éléments polyédriques. */
1863
1864 def_cel_fac_pos[icel + 1] = def_cel_fac_pos[icel] + ifac;
1865
1866 } /* Fin de la boucle sur les éléments */
1867
1868 /* Mise à jour de la table */
1869 /*-------------------------*/
1870
1871 ECS_FREE(table_def_cel->pos);
1872 ECS_FREE(table_def_cel->val);
1873
1874 table_def_cel->pos = def_cel_fac_pos;
1875 table_def_cel->val = def_cel_fac_val;
1876
1877 /* ecs_table__pos_en_regle(*table_def_fac); */
1878 ecs_table__pos_en_regle(table_def_cel);
1879 }
1880
1881 /*----------------------------------------------------------------------------
1882 * Fonction qui realise la fusion des definitions des elements
1883 *----------------------------------------------------------------------------*/
1884
1885 ecs_tab_int_t
ecs_table_def__fusionne(ecs_table_t * table_def,size_t * nbr_elt_cpct,ecs_tab_int_t * signe_elt)1886 ecs_table_def__fusionne(ecs_table_t *table_def,
1887 size_t *nbr_elt_cpct,
1888 ecs_tab_int_t *signe_elt)
1889 {
1890 size_t cpt_sup_fin;
1891 size_t ind_cmp;
1892 size_t ind_inf;
1893 size_t ind_loc_sup;
1894 size_t ind_loc_cmp;
1895 size_t ind_pos;
1896 size_t ind_pos_loc;
1897 size_t ind_sup;
1898 size_t nbr_sup_ini;
1899 size_t nbr_inf;
1900 ecs_int_t num_inf;
1901 ecs_int_t num_inf_loc;
1902 ecs_int_t num_inf_min;
1903 ecs_int_t num_inf_min_cmp;
1904 size_t pos_cmp;
1905 size_t pos_cpt;
1906 size_t pos_sup;
1907 int sgn;
1908
1909 size_t ind_pos_sup[3];
1910 size_t ind_pos_cmp[3];
1911
1912 ecs_tab_int_t cpt_ref_inf;
1913
1914 ecs_size_t *pos_recherche = NULL;
1915 ecs_int_t *val_recherche = NULL;
1916
1917 ecs_tab_int_t tab_transf; /* Tableau de transformation */
1918
1919 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1920
1921 ecs_table__regle_en_pos(table_def);
1922
1923 /* Fusion des definitions des elements sur la liste initiale */
1924 /*------------------------------------------------------------*/
1925
1926 /* Initialisations avec précautions pour cas vide */
1927
1928 tab_transf.nbr = 0;
1929 tab_transf.val = NULL;
1930
1931 if (table_def == NULL)
1932 return tab_transf;
1933
1934 if (table_def == NULL)
1935 return tab_transf;
1936
1937 nbr_sup_ini = table_def->nbr;
1938 cpt_sup_fin = 0;
1939
1940 if (nbr_sup_ini < 1)
1941 return tab_transf;
1942
1943 if (signe_elt != NULL) {
1944 signe_elt->nbr = nbr_sup_ini;
1945 ECS_MALLOC(signe_elt->val, nbr_sup_ini, ecs_int_t);
1946 }
1947
1948 /* Comptage du nombre de sous-entités */
1949
1950 nbr_inf = 0;
1951
1952 for (ind_sup = 0;
1953 ind_sup < table_def->pos[nbr_sup_ini] - 1;
1954 ind_sup++) {
1955 num_inf = table_def->val[ind_sup];
1956 if ((size_t)(ECS_ABS(num_inf)) > nbr_inf)
1957 nbr_inf = ECS_ABS(num_inf);
1958 }
1959
1960 /* Construction du tableau de recherche des sous-entités */
1961 /*-------------------------------------------------------*/
1962
1963 ECS_MALLOC(pos_recherche, nbr_inf + 1, ecs_size_t);
1964 ECS_MALLOC(val_recherche, table_def->nbr, ecs_int_t);
1965
1966 /*
1967 Comptage pour chaque élément de l'entité inférieure du nombre de fois où il
1968 apparaît comme élement d'indice minimal d'un élément de l'entité supérieure
1969 */
1970
1971 cpt_ref_inf = ecs_tab_int__cree_init(nbr_inf, 0);
1972
1973 for (ind_sup = 0; ind_sup < table_def->nbr; ind_sup++) {
1974
1975 ind_pos_loc = table_def->pos[ind_sup] - 1;
1976
1977 num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
1978 while (++ind_pos_loc < table_def->pos[ind_sup + 1] -1) {
1979 num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
1980 if (num_inf_loc < num_inf_min)
1981 num_inf_min = num_inf_loc;
1982 }
1983
1984 assert(num_inf_min > 0 && (size_t)num_inf_min <= nbr_inf);
1985
1986 cpt_ref_inf.val[num_inf_min - 1] += 1;
1987 }
1988
1989
1990 /* Construction du vecteur */
1991
1992 pos_recherche[0] = 1;
1993
1994 for (ind_inf = 0; ind_inf < nbr_inf; ind_inf++) {
1995
1996 pos_recherche[ind_inf + 1]
1997 = pos_recherche[ind_inf] + cpt_ref_inf.val[ind_inf];
1998
1999 cpt_ref_inf.val[ind_inf] = 0;
2000
2001 }
2002
2003
2004 for (ind_sup = 0; ind_sup < table_def->nbr; ind_sup++) {
2005
2006 ind_pos_loc = table_def->pos[ind_sup] - 1;
2007
2008 num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
2009 while (ind_pos_loc < table_def->pos[ind_sup + 1] -1) {
2010 num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
2011 ind_pos_loc++;
2012 if (num_inf_loc < num_inf_min)
2013 num_inf_min = num_inf_loc;
2014 }
2015
2016 ind_inf = num_inf_min - 1;
2017
2018 ind_pos = pos_recherche[ind_inf] - 1
2019 + cpt_ref_inf.val[ind_inf];
2020
2021 cpt_ref_inf.val[ind_inf] += 1;
2022
2023 val_recherche[ind_pos] = ind_sup + 1;
2024
2025 }
2026
2027
2028 /* Libération du tableau auxiliaire */
2029
2030 cpt_ref_inf.nbr = 0;
2031 ECS_FREE(cpt_ref_inf.val);
2032
2033 /* Allocation du tableau de transformation */
2034
2035 tab_transf = ecs_tab_int__cree(nbr_sup_ini);
2036
2037 for (ind_sup = 0; ind_sup < nbr_sup_ini; ind_sup++)
2038 tab_transf.val[ind_sup] = ind_sup;
2039
2040
2041 /* Boucle principale de recherche sur les éléments supérieurs */
2042 /*------------------------------------------------------------*/
2043
2044 /*
2045 Le premier élément ne peut pas être fusionné avec un élément précédent,
2046 la boucle commence donc au deuxième
2047 */
2048
2049 tab_transf.val[0] = 0;
2050
2051 if (signe_elt != NULL)
2052 signe_elt->val[0] = 1;
2053
2054 cpt_sup_fin = 1;
2055
2056 for (ind_sup = 1; ind_sup < table_def->nbr; ind_sup++) {
2057
2058 /* Recherche élement entité inférieure de plus petit numéro référencé */
2059
2060 ind_pos_sup[0] = table_def->pos[ind_sup ] - 1; /* début */
2061 ind_pos_sup[1] = table_def->pos[ind_sup + 1] - 1; /* fin */
2062 ind_pos_sup[2] = table_def->pos[ind_sup ] - 1; /* plus petit */
2063
2064 ind_pos_loc = ind_pos_sup[0];
2065
2066 num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
2067 while (++ind_pos_loc < ind_pos_sup[1]) {
2068 num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
2069 if (num_inf_loc < num_inf_min) {
2070 num_inf_min = num_inf_loc;
2071 ind_pos_sup[2] = ind_pos_loc;
2072 }
2073 }
2074
2075 /*
2076 On cherche des éléments de l'entité courante de plus petit numéro que
2077 l'entité courante ayant même plus petit élément de l'entité inférieure
2078 (recherche de candidats pour la fusion)
2079 */
2080
2081 ind_inf = num_inf_min - 1;
2082 sgn = 0;
2083
2084 for (pos_cmp = pos_recherche[ind_inf] - 1;
2085 pos_cmp < pos_recherche[ind_inf + 1] - 1;
2086 pos_cmp++) {
2087
2088 ind_cmp = val_recherche[pos_cmp] - 1;
2089
2090 /* Repérage point de départ pour comparaison */
2091
2092 if (ind_cmp < ind_sup) {
2093
2094 ind_pos_cmp[0] = table_def->pos[ind_cmp ] - 1; /* début */
2095 ind_pos_cmp[1] = table_def->pos[ind_cmp + 1] - 1; /* fin */
2096 ind_pos_cmp[2] = table_def->pos[ind_cmp ] - 1; /* plus petit */
2097
2098 assert(ind_pos_cmp[1] > ind_pos_cmp[0]);
2099
2100 ind_pos_cmp[1] = table_def->pos[ind_cmp + 1] - 1; /* fin */
2101
2102 ind_pos_loc = ind_pos_cmp[0];
2103
2104 num_inf_min_cmp = ECS_ABS(table_def->val[ind_pos_loc]);
2105 while ((++ind_pos_loc) < ind_pos_cmp[1]) {
2106 num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
2107 if (num_inf_loc < num_inf_min_cmp) {
2108 num_inf_min_cmp = num_inf_loc;
2109 ind_pos_cmp[2] = ind_pos_loc;
2110 }
2111 }
2112
2113 /* Comparaison des définitions */
2114
2115 for (sgn = 1; sgn > -2; sgn -= 2) {
2116
2117 ind_loc_sup = ind_pos_sup[2];
2118 ind_loc_cmp = ind_pos_cmp[2];
2119
2120 do {
2121
2122 ind_loc_sup++;
2123 if (ind_loc_sup == ind_pos_sup[1])
2124 ind_loc_sup = ind_pos_sup[0];
2125
2126 ind_loc_cmp += sgn;
2127 if (ind_loc_cmp == ind_pos_cmp[1])
2128 ind_loc_cmp = ind_pos_cmp[0];
2129 else if ( ind_loc_cmp < ind_pos_cmp[0]
2130 || ind_loc_cmp > ind_pos_cmp[1])
2131 ind_loc_cmp = ind_pos_cmp[1] - 1;
2132
2133 } while ( ( ECS_ABS(table_def->val[ind_loc_sup])
2134 == ECS_ABS(table_def->val[ind_loc_cmp]))
2135 && ind_loc_sup != ind_pos_sup[2]
2136 && ind_loc_cmp != ind_pos_cmp[2]);
2137
2138 if ( ind_loc_sup == ind_pos_sup[2]
2139 && ind_loc_cmp == ind_pos_cmp[2])
2140 break; /* Sortie boucle sur signe parcours (1, -1, -3) */
2141
2142 }
2143
2144 /*
2145 Si sgn = 1, les entités sont confondues, de même sens;
2146 Si sgn = -1, elles sont confondues, de sens inverse;
2147 Sinon, elles ne sont pas confondues
2148 */
2149
2150 if (sgn == 1 || sgn == -1)
2151 break; /* Sortie boucle sur pos_cmp pour recherche de candidats */
2152
2153 } /* Fin de la comparaison référence/candidat (if ind_cmp < ind_sup) */
2154
2155 } /* Fin boucle sur pos_cmp recherche de candidats */
2156
2157
2158 /* Si on a trouvé une entité à fusionner */
2159
2160 if (sgn == 1 || sgn == -1) {
2161
2162 assert(pos_cmp < pos_recherche[ind_inf + 1] - 1);
2163
2164 if (sgn == 1 && (ind_pos_cmp[1] - ind_pos_cmp[0] == 2)) {
2165
2166 /*
2167 Si l'on a deux entités inférieures par entité supérieure,
2168 la permutation cyclique peut trouver une égalité quel
2169 que soit le signe; on le corrige si nécessaire
2170 */
2171
2172 if ( ECS_ABS(table_def->val[ind_pos_sup[0]])
2173 != ECS_ABS(table_def->val[ind_pos_cmp[0]]))
2174 sgn = -1;
2175
2176 }
2177
2178 tab_transf.val[ind_sup] = tab_transf.val[ind_cmp];
2179
2180 if (signe_elt != NULL)
2181 signe_elt->val[ind_sup] = sgn;
2182
2183 }
2184 else {
2185
2186 tab_transf.val[ind_sup] = cpt_sup_fin++;
2187
2188 if (signe_elt != NULL)
2189 signe_elt->val[ind_sup] = 1;
2190
2191 }
2192
2193
2194 } /* Fin boucle sur les éléments de l'entité supérieure (que l'on fusionne) */
2195
2196 /* Libération du tableau de recherche */
2197
2198 ECS_FREE(pos_recherche);
2199 ECS_FREE(val_recherche);
2200
2201 /* Compactage de la définition */
2202 /*-----------------------------*/
2203
2204 /*
2205 Le premier élément ne peut pas être fusionné avec un élément précédent,
2206 la boucle commence donc au deuxième
2207 */
2208
2209 cpt_sup_fin = 1;
2210
2211 for (ind_sup = 1; ind_sup < table_def->nbr; ind_sup++) {
2212
2213 if (tab_transf.val[ind_sup] > (ecs_int_t)(cpt_sup_fin - 1)) {
2214
2215 /* Recopie définition sur partie compactée du tableau */
2216
2217 for (pos_cpt = table_def->pos[cpt_sup_fin],
2218 pos_sup = table_def->pos[ind_sup];
2219 pos_sup < table_def->pos[ind_sup + 1];
2220 pos_cpt++, pos_sup++)
2221 table_def->val[pos_cpt - 1] = table_def->val[pos_sup - 1];
2222
2223 cpt_sup_fin += 1;
2224
2225 table_def->pos[cpt_sup_fin] = pos_cpt;
2226
2227 }
2228 }
2229
2230 /* Redimensionnement de l'entité compactée */
2231
2232 table_def->nbr = cpt_sup_fin;
2233 ECS_REALLOC(table_def->pos, cpt_sup_fin + 1, ecs_size_t);
2234 ECS_REALLOC(table_def->val,
2235 table_def->pos[cpt_sup_fin] - 1,
2236 ecs_int_t);
2237
2238 ecs_table__pos_en_regle(table_def);
2239
2240 /* Affectation de la nouvelle liste compactee */
2241 /*--------------------------------------------*/
2242
2243 *nbr_elt_cpct = table_def->nbr;
2244
2245 /* Renvoi du tableau de transformation */
2246 /*-------------------------------------*/
2247
2248 return tab_transf;
2249 }
2250
2251 /*----------------------------------------------------------------------------
2252 * Fonction qui construit la liste des cellules attachées à une liste
2253 * de faces fournie en argument.
2254 *----------------------------------------------------------------------------*/
2255
2256 ecs_tab_int_t
ecs_table_def__liste_cel_fac(const size_t nbr_fac,ecs_table_t * table_def_cel,const ecs_tab_int_t liste_fac)2257 ecs_table_def__liste_cel_fac(const size_t nbr_fac,
2258 ecs_table_t *table_def_cel,
2259 const ecs_tab_int_t liste_fac)
2260 {
2261 size_t cpt_cel;
2262 size_t icel;
2263 size_t ifac;
2264 size_t iloc;
2265 size_t ipos;
2266 size_t nbr_cel;
2267 size_t nbr_fac_cel;
2268 ecs_int_t *indic_cel;
2269 ecs_int_t *indic_fac;
2270
2271 ecs_tab_int_t liste_cel;
2272
2273 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2274
2275 assert(table_def_cel != NULL);
2276
2277 ecs_table__regle_en_pos(table_def_cel);
2278
2279 /* Initialisations */
2280
2281 cpt_cel = 0;
2282
2283 nbr_cel = table_def_cel->nbr;
2284
2285 /* Indicateurs */
2286
2287 ECS_MALLOC(indic_cel, nbr_cel, ecs_int_t);
2288 ECS_MALLOC(indic_fac, nbr_fac, ecs_int_t);
2289
2290 for (icel = 0; icel < nbr_cel; icel++)
2291 indic_cel[icel] = 0;
2292
2293 for (ifac = 0; ifac < nbr_fac; ifac++)
2294 indic_fac[ifac] = 0;
2295
2296 for (ifac = 0; ifac < liste_fac.nbr; ifac++)
2297 indic_fac[liste_fac.val[ifac]] = 1;
2298
2299 /* Première boucle sur les cellules : marquage */
2300
2301 for (icel = 0; icel < nbr_cel; icel++) {
2302
2303 nbr_fac_cel
2304 = table_def_cel->pos[icel + 1]
2305 - table_def_cel->pos[icel];
2306
2307 ipos = table_def_cel->pos[icel] - 1;
2308
2309 for (iloc = 0; iloc < nbr_fac_cel; iloc++) {
2310
2311 ifac = ECS_ABS(table_def_cel->val[ipos]) - 1;
2312
2313 ipos++;
2314
2315 if (indic_fac[ifac] == 1 && indic_cel[icel] == 0) {
2316 indic_cel[icel] = 1;
2317 cpt_cel++;
2318 }
2319
2320 }
2321 }
2322
2323 ECS_FREE(indic_fac);
2324
2325 /* Seconde boucle sur les cellules : remplissage */
2326
2327 liste_cel.nbr = cpt_cel;
2328 ECS_MALLOC(liste_cel.val, liste_cel.nbr, ecs_int_t);
2329
2330 cpt_cel = 0;
2331
2332 for (icel = 0; icel < nbr_cel; icel++) {
2333
2334 if (indic_cel[icel] == 1)
2335 liste_cel.val[cpt_cel++] = icel;
2336
2337 }
2338
2339 ECS_FREE(indic_cel);
2340
2341 ecs_table__libere_pos(table_def_cel);
2342
2343 return liste_cel;
2344 }
2345
2346 /*----------------------------------------------------------------------------
2347 * Fonction qui remplace les références à des éléments
2348 * en des références à d'autres éléments liés aux premiers
2349 * par un tableau de renumérotation qui peut être signé.
2350 *----------------------------------------------------------------------------*/
2351
2352 void
ecs_table_def__remplace_ref(ecs_table_t * table_def,ecs_tab_int_t * tab_old_new)2353 ecs_table_def__remplace_ref(ecs_table_t *table_def,
2354 ecs_tab_int_t *tab_old_new)
2355 {
2356 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2357
2358 if (table_def != NULL && tab_old_new != NULL) {
2359
2360 ecs_int_t num_old, num_new;
2361 ecs_size_t ielt;
2362 ecs_size_t ival, ival_deb, ival_fin;
2363 ecs_size_t cpt_val = 0;
2364 size_t nbr_elt = table_def->nbr;
2365
2366 ecs_table__regle_en_pos(table_def);
2367
2368 ival_deb = 0;
2369
2370 for (ielt = 0; ielt < nbr_elt; ielt++) {
2371
2372 ival_fin = table_def->pos[ielt+1] - 1;
2373
2374 for (ival = ival_deb; ival < ival_fin; ival++) {
2375 num_new = 0;
2376 num_old = table_def->val[ival];
2377 if (num_old > 0)
2378 num_new = tab_old_new->val[num_old -1];
2379 else
2380 num_new = -tab_old_new->val[-num_old -1];
2381 if (num_new != 0)
2382 table_def->val[cpt_val++] = num_new;
2383 }
2384
2385 ival_deb = table_def->pos[ielt+1] - 1;
2386 table_def->pos[ielt+1] = cpt_val + 1;
2387 }
2388
2389 ECS_REALLOC(table_def->val, cpt_val, ecs_int_t);
2390
2391 ecs_table__pos_en_regle(table_def);
2392 }
2393 }
2394
2395 /*----------------------------------------------------------------------------
2396 * Fonction qui construit un tableau de booleens conforme a une liste
2397 * de sous-elements
2398 * Un sous-element est a `true'
2399 * s'il intervient dans la definition des elements
2400 *----------------------------------------------------------------------------*/
2401
2402 void
ecs_table_def__cree_masque(bool sselt_select[],ecs_table_t * table_def_elt)2403 ecs_table_def__cree_masque(bool sselt_select[],
2404 ecs_table_t *table_def_elt)
2405 {
2406 size_t nbr_elt;
2407 size_t num_sselt;
2408 size_t pos_inf;
2409 size_t pos_sup;
2410
2411 size_t ielt;
2412 size_t ipos;
2413
2414 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2415
2416 ecs_table__regle_en_pos(table_def_elt);
2417
2418 nbr_elt = table_def_elt->nbr;
2419
2420 for (ielt = 0; ielt < nbr_elt; ielt++) {
2421
2422 pos_inf = table_def_elt->pos[ielt ] - 1;
2423 pos_sup = table_def_elt->pos[ielt + 1] - 1;
2424
2425 for (ipos = pos_inf; ipos < pos_sup; ipos++) {
2426
2427 num_sselt = ECS_ABS(table_def_elt->val[ipos]) - 1;
2428
2429 /* If the sub-element has not yet been accounted for, mark it */
2430
2431 if (sselt_select[num_sselt] == false)
2432 sselt_select[num_sselt] = true;
2433
2434 }
2435 }
2436
2437 ecs_table__libere_pos(table_def_elt);
2438 }
2439
2440 /*----------------------------------------------------------------------------
2441 * Suppression des sommets ne participant pas à la connectivité
2442 * et mise à jour de la connectivité.
2443 *----------------------------------------------------------------------------*/
2444
2445 void
ecs_table_def__nettoie_nodal(size_t * n_vertices,ecs_coord_t ** vtx_coords,ecs_table_t * table_def_fac,ecs_table_t * table_def_cel)2446 ecs_table_def__nettoie_nodal(size_t *n_vertices,
2447 ecs_coord_t **vtx_coords,
2448 ecs_table_t *table_def_fac,
2449 ecs_table_t *table_def_cel)
2450 {
2451 size_t cpt_som;
2452 size_t iloc;
2453 size_t ipos_new;
2454 size_t ipos_old;
2455 size_t isom;
2456 size_t ival;
2457 size_t nbr_som;
2458 size_t nbr_val;
2459
2460 ecs_tab_int_t tab_som_old_new;
2461
2462 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2463
2464 tab_som_old_new.nbr = 0;
2465 tab_som_old_new.val = NULL;
2466
2467 if (vtx_coords == NULL)
2468 return;
2469
2470 /* Initialisation du vecteur de renumérotation comme marqueur */
2471
2472 nbr_som = *n_vertices;
2473
2474 tab_som_old_new = ecs_tab_int__cree_init(nbr_som, 0);
2475
2476 /* faces */
2477
2478 if (table_def_fac != NULL) {
2479
2480 nbr_val = ecs_table__ret_val_nbr(table_def_fac);
2481
2482 for (ival = 0; ival < nbr_val; ival++)
2483 tab_som_old_new.val[table_def_fac->val[ival] - 1] = 1;
2484 }
2485
2486 /* cellules */
2487
2488 if (table_def_cel != NULL) {
2489
2490 nbr_val = ecs_table__ret_val_nbr(table_def_cel);
2491
2492 for (ival = 0; ival < nbr_val; ival++)
2493 tab_som_old_new.val[table_def_cel->val[ival] - 1] = 1;
2494 }
2495
2496 /* Préparation de la renumérotation */
2497
2498 cpt_som = 0;
2499
2500 for (isom = 0; isom < nbr_som; isom++) {
2501 if (tab_som_old_new.val[isom] != 0)
2502 cpt_som++;
2503 }
2504
2505 /* Si tous les sommets sont référencés, rien à faire */
2506
2507 if (cpt_som == nbr_som) {
2508 tab_som_old_new.nbr = 0;
2509 ECS_FREE(tab_som_old_new.val);
2510 return;
2511 }
2512
2513 /* Compactage du tableau des sommets
2514 et transformation du marquer en renumérotation */
2515
2516 cpt_som = 0;
2517
2518 for (isom = 0; isom < nbr_som; isom++) {
2519
2520 if (tab_som_old_new.val[isom] != 0) {
2521
2522 ipos_new = 3 * (cpt_som);
2523 ipos_old = 3 * isom;
2524
2525 for (iloc = 0; iloc < 3; iloc++)
2526 (*vtx_coords)[ipos_new++] = (*vtx_coords)[ipos_old++];
2527
2528 tab_som_old_new.val[isom] = cpt_som + 1;
2529
2530 cpt_som++;
2531
2532 }
2533 }
2534
2535 *n_vertices = cpt_som;
2536 ECS_REALLOC(*vtx_coords, cpt_som*3, ecs_coord_t);
2537
2538 /* Mise à jour de la connectivité nodale */
2539
2540 if (table_def_fac != NULL)
2541 ecs_table_def__remplace_ref(table_def_fac,
2542 &tab_som_old_new);
2543
2544 if (table_def_cel != NULL)
2545 ecs_table_def__remplace_ref(table_def_cel,
2546 &tab_som_old_new);
2547
2548 tab_som_old_new.nbr = 0;
2549 ECS_FREE(tab_som_old_new.val);
2550 }
2551
2552 /*----------------------------------------------------------------------------
2553 * Correction si nécessaire de l'orientation des éléments en connectivité
2554 * nodale. L'argument liste_cel_err est optionnel.
2555 *----------------------------------------------------------------------------*/
2556
2557 void
ecs_table_def__orient_nodal(ecs_coord_t * vtx_coords,ecs_table_t * table_def_fac,ecs_table_t * table_def_cel,ecs_tab_int_t * liste_cel_err,bool correc_orient)2558 ecs_table_def__orient_nodal(ecs_coord_t *vtx_coords,
2559 ecs_table_t *table_def_fac,
2560 ecs_table_t *table_def_cel,
2561 ecs_tab_int_t *liste_cel_err,
2562 bool correc_orient)
2563 {
2564 size_t ipos_cel;
2565 size_t ipos_fac;
2566 size_t icel;
2567 size_t ifac;
2568 size_t nbr_som_loc;
2569 size_t nbr_fac;
2570 size_t nbr_cel;
2571 ecs_int_t typ_elt;
2572 ecs_int_t ret_orient;
2573
2574 ecs_int_t cpt_orient_correc[ECS_ELT_TYP_FIN];
2575 ecs_int_t cpt_orient_erreur[ECS_ELT_TYP_FIN];
2576
2577 ecs_int_t typ_cel_base[9] = {ECS_ELT_TYP_NUL,
2578 ECS_ELT_TYP_NUL,
2579 ECS_ELT_TYP_NUL,
2580 ECS_ELT_TYP_NUL,
2581 ECS_ELT_TYP_CEL_TETRA,
2582 ECS_ELT_TYP_CEL_PYRAM,
2583 ECS_ELT_TYP_CEL_PRISM,
2584 ECS_ELT_TYP_NUL,
2585 ECS_ELT_TYP_CEL_HEXA};
2586
2587 ecs_int_t cpt_cel_erreur = 0;
2588
2589 int foam2vtk = 0; /* tentative adaptation to meshes transformed with
2590 foam2vtk; solves most issues but a few remain */
2591 if (getenv("CS_PREPROCESS_FOAM_2_VTK_SOURCE") != NULL)
2592 foam2vtk = 1;
2593
2594 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2595
2596 if (vtx_coords == NULL)
2597 return;
2598
2599 if (table_def_fac != NULL)
2600 ecs_table__regle_en_pos(table_def_fac);
2601
2602 if (table_def_cel != NULL)
2603 ecs_table__regle_en_pos(table_def_cel);
2604
2605 assert(vtx_coords != NULL);
2606
2607 for (typ_elt = 0; typ_elt < ECS_ELT_TYP_FIN; typ_elt++) {
2608 cpt_orient_correc[typ_elt] = 0;
2609 cpt_orient_erreur[typ_elt] = 0;
2610 }
2611
2612 printf(_("\n Element orientation check.\n\n"));
2613
2614 /* faces */
2615
2616 if (table_def_fac != NULL) {
2617
2618 nbr_fac = table_def_fac->nbr;
2619
2620 for (ifac = 0; ifac < nbr_fac; ifac++) {
2621
2622 ipos_fac = table_def_fac->pos[ifac ] - 1;
2623 nbr_som_loc = table_def_fac->pos[ifac + 1] - 1 - ipos_fac;
2624
2625 if (nbr_som_loc == 4) {
2626
2627 ret_orient
2628 = _orient_quad(vtx_coords,
2629 &(table_def_fac->val[ipos_fac]),
2630 correc_orient);
2631
2632 if (ret_orient < 0)
2633 cpt_orient_erreur[ECS_ELT_TYP_FAC_QUAD] += 1;
2634 else if (ret_orient > 0)
2635 cpt_orient_correc[ECS_ELT_TYP_FAC_QUAD] += 1;
2636
2637 }
2638 }
2639 }
2640
2641 /* cellules */
2642
2643 if (table_def_cel != NULL) {
2644
2645 ecs_tab_int_t face_index, face_marker, edges;
2646
2647 face_index.nbr = 0;
2648 face_index.val = NULL;
2649
2650 face_marker.nbr = 0;
2651 face_marker.val = NULL;
2652
2653 edges.nbr = 0;
2654 edges.val = NULL;
2655
2656 nbr_cel = table_def_cel->nbr;
2657
2658 for (icel = 0; icel < nbr_cel; icel++) {
2659
2660 ipos_cel = table_def_cel->pos[icel] - 1;
2661
2662 nbr_som_loc = table_def_cel->pos[icel + 1] - 1 - ipos_cel;
2663
2664 if (nbr_som_loc < 9)
2665 typ_elt = typ_cel_base[nbr_som_loc];
2666 else
2667 typ_elt = ECS_ELT_TYP_CEL_POLY;
2668
2669 switch (typ_elt) {
2670
2671 case ECS_ELT_TYP_CEL_TETRA:
2672
2673 if (foam2vtk)
2674 ret_orient = _orient_tetra_of(&(table_def_cel->val[ipos_cel]));
2675 else
2676 ret_orient = _orient_tetra(vtx_coords,
2677 &(table_def_cel->val[ipos_cel]));
2678 break;
2679
2680 case ECS_ELT_TYP_CEL_PYRAM:
2681
2682 if (foam2vtk)
2683 ret_orient
2684 = _orient_pyram_of(&(table_def_cel->val[ipos_cel]));
2685 else
2686 ret_orient
2687 = _orient_pyram(vtx_coords,
2688 &(table_def_cel->val[ipos_cel]),
2689 correc_orient);
2690 break;
2691
2692 case ECS_ELT_TYP_CEL_PRISM:
2693
2694 if (foam2vtk)
2695 ret_orient = 0;
2696 else
2697 ret_orient
2698 = _orient_prism(vtx_coords,
2699 &(table_def_cel->val[ipos_cel]),
2700 correc_orient);
2701 break;
2702
2703 case ECS_ELT_TYP_CEL_HEXA:
2704
2705 if (foam2vtk)
2706 ret_orient
2707 = _orient_hexa_of(&(table_def_cel->val[ipos_cel]));
2708 else
2709 ret_orient
2710 = _orient_hexa(vtx_coords,
2711 &(table_def_cel->val[ipos_cel]),
2712 correc_orient);
2713 break;
2714
2715 default: /* ECS_ELT_TYP_CEL_POLY */
2716
2717 ret_orient = _orient_polyhedron(vtx_coords,
2718 &(table_def_cel->val[ipos_cel]),
2719 ( table_def_cel->pos[icel + 1]
2720 - table_def_cel->pos[icel]),
2721 correc_orient,
2722 &face_index,
2723 &face_marker,
2724 &edges);
2725
2726 };
2727
2728 if (ret_orient < 0) {
2729 cpt_orient_erreur[typ_elt] += 1;
2730 if (liste_cel_err != NULL) {
2731 if (liste_cel_err->nbr == 0) {
2732 liste_cel_err->nbr = nbr_cel;
2733 ECS_MALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t);
2734 }
2735 liste_cel_err->val[cpt_cel_erreur] = icel;
2736 }
2737 cpt_cel_erreur += 1;
2738 }
2739
2740 else if (ret_orient > 0) {
2741 cpt_orient_correc[typ_elt] += 1;
2742 }
2743 }
2744
2745 if (face_index.nbr > 0)
2746 ECS_FREE(face_index.val);
2747
2748 if (face_marker.nbr > 0)
2749 ECS_FREE(face_marker.val);
2750
2751 if (edges.nbr > 0)
2752 ECS_FREE(edges.val);
2753 }
2754
2755 /* Impression de messages d'avertissement */
2756
2757 for (typ_elt = 0; typ_elt < ECS_ELT_TYP_FIN; typ_elt++) {
2758 if (cpt_orient_correc[typ_elt] > 0) {
2759 ecs_warn();
2760 printf(_("%d elements of type %s had to be re-oriented\n"),
2761 (int)(cpt_orient_correc[typ_elt]),
2762 _(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2763 }
2764 if (cpt_orient_erreur[typ_elt] > 0) {
2765 if (correc_orient == true) {
2766 ecs_warn();
2767 printf(_("%d elements of type %s were impossible to re-orient\n"),
2768 (int)(cpt_orient_erreur[typ_elt]),
2769 _(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2770 }
2771 else {
2772 ecs_warn();
2773 printf(_("%d elements of type %s are mis-oriented or highly warped\n"),
2774 (int)(cpt_orient_erreur[typ_elt]),
2775 _(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2776 }
2777 }
2778 }
2779
2780 /* Redimensionnement de la liste des cellules toujours mal orientées */
2781
2782 if (liste_cel_err != NULL) {
2783 if (liste_cel_err->nbr > 0) {
2784 liste_cel_err->nbr = cpt_cel_erreur;
2785 ECS_REALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t);
2786 }
2787 }
2788
2789 if (table_def_fac != NULL)
2790 ecs_table__libere_pos(table_def_fac);
2791
2792 if (table_def_cel != NULL)
2793 ecs_table__libere_pos(table_def_cel);
2794 }
2795
2796 /*----------------------------------------------------------------------------
2797 * Fusion des sommets confondus d'après la longueur des arêtes des faces.
2798 * La connectivité des faces est mise à jour.
2799 *----------------------------------------------------------------------------*/
2800
2801 void
ecs_table_def__nettoie_som_fac(size_t * n_vertices,ecs_coord_t ** vtx_coords,ecs_table_t * table_def_fac)2802 ecs_table_def__nettoie_som_fac(size_t *n_vertices,
2803 ecs_coord_t **vtx_coords,
2804 ecs_table_t *table_def_fac)
2805 {
2806 size_t cpt_fusion;
2807
2808 size_t icoo;
2809 size_t ifac;
2810
2811 size_t ipos_deb;
2812 size_t ipos_0;
2813 size_t ipos_1;
2814 size_t isom;
2815 size_t isom_0;
2816 size_t isom_1;
2817
2818 size_t nbr_fac;
2819 size_t nbr_som;
2820 size_t nbr_som_fac;
2821
2822 double delta_coo;
2823 double lng_are_min;
2824 double lng_are_2;
2825 double lng_min_2;
2826
2827 ecs_tab_int_t tab_equiv_som;
2828
2829 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2830
2831 if (vtx_coords == NULL || table_def_fac == NULL)
2832 return;
2833
2834 ecs_table__regle_en_pos(table_def_fac);
2835
2836 nbr_fac = table_def_fac->nbr;
2837 nbr_som = *n_vertices;
2838
2839 lng_are_min = 1.0e-15;
2840 if (getenv("CS_PREPROCESS_MIN_EDGE_LEN") != NULL) {
2841 lng_are_min = atof(getenv("CS_PREPROCESS_MIN_EDGE_LEN"));
2842 }
2843 lng_min_2 = lng_are_min * ECS_ABS(lng_are_min);
2844
2845 /* Marquage des sommets */
2846 /*----------------------*/
2847
2848 tab_equiv_som = ecs_tab_int__cree_init(nbr_som, -1);
2849
2850 cpt_fusion = 0;
2851
2852 /* Boucle sur les faces pour détecter les sommets confondus */
2853 /*----------------------------------------------------------*/
2854
2855 for (ifac = 0; ifac < nbr_fac; ifac++) {
2856
2857 ipos_deb = table_def_fac->pos[ifac] - 1;
2858 nbr_som_fac = (table_def_fac->pos[ifac + 1] - 1) - ipos_deb;
2859
2860 for (isom = 0; isom < nbr_som_fac; isom++) {
2861
2862 isom_0 = table_def_fac->val[ipos_deb + isom] - 1;
2863 isom_1 = table_def_fac->val[ ipos_deb + ((isom + 1)%nbr_som_fac)] -1;
2864 ipos_0 = isom_0 * 3;
2865 ipos_1 = isom_1 * 3;
2866
2867 lng_are_2 = 0.0;
2868
2869 for (icoo = 0; icoo < 3; icoo++) {
2870 delta_coo = (*vtx_coords)[ipos_1 + icoo]
2871 - (*vtx_coords)[ipos_0 + icoo];
2872 lng_are_2 += delta_coo * delta_coo;
2873 }
2874
2875 /* Sommets confondus détectés */
2876
2877 if (lng_are_2 < lng_min_2) {
2878 _table_def__maj_equiv_som(isom_0, isom_1, &tab_equiv_som);
2879 cpt_fusion += 1;
2880 }
2881
2882 }
2883 }
2884
2885 /* Fusion si nécessaire */
2886
2887 if (cpt_fusion > 0) {
2888
2889 ecs_tab_int_t tab_som_old_new;
2890
2891 printf(_("\nMesh verification:\n\n"
2892 " %d vertices belong to edges of length less than %g\n"
2893 " and will be merged; (this tolerance may be modified\n"
2894 " using the CS_PREPROCESS_MIN_EDGE_LEN environment variable).\n"),
2895 (int)cpt_fusion, lng_are_min);
2896
2897 tab_som_old_new = _table_def__fusion_som(n_vertices,
2898 vtx_coords,
2899 &tab_equiv_som);
2900
2901 _table_def__maj_fac_som(table_def_fac, &tab_som_old_new);
2902
2903 ECS_FREE(tab_som_old_new.val);
2904
2905 printf("\n");
2906 }
2907
2908 /* Libération mémoire et retour */
2909
2910 ECS_FREE(tab_equiv_som.val);
2911
2912 ecs_table__pos_en_regle(table_def_fac);
2913 }
2914
2915 /*----------------------------------------------------------------------------
2916 * Fonction qui supprime les éventuelles faces dégénérées
2917 *----------------------------------------------------------------------------*/
2918
2919 ecs_tab_int_t
ecs_table_def__nettoie_fac(ecs_table_t * table_def_fac)2920 ecs_table_def__nettoie_fac(ecs_table_t *table_def_fac)
2921 {
2922 ecs_int_t nbr_fac_old;
2923 ecs_int_t nbr_fac_new;
2924
2925 ecs_int_t ind_fac;
2926
2927 ecs_int_t cpt_fac;
2928 ecs_int_t cpt_val;
2929 ecs_int_t ind_val;
2930 ecs_int_t ind_val_deb;
2931 ecs_int_t ind_val_fin;
2932
2933 ecs_tab_int_t tab_fac_old_new;
2934
2935 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2936
2937 tab_fac_old_new.nbr = 0;
2938 tab_fac_old_new.val = NULL;
2939
2940 if (table_def_fac == NULL)
2941 return tab_fac_old_new;
2942
2943 ecs_table__regle_en_pos(table_def_fac);
2944
2945 /* Initialisations */
2946
2947 nbr_fac_old = table_def_fac->nbr;
2948
2949 /* Boucle sur les faces (renumérotation) */
2950 /* ------------------------------------- */
2951
2952 tab_fac_old_new.nbr = nbr_fac_old;
2953 ECS_MALLOC(tab_fac_old_new.val, tab_fac_old_new.nbr, ecs_int_t);
2954
2955 cpt_fac = 0;
2956 cpt_val = 0;
2957
2958 for (ind_fac = 0; ind_fac < nbr_fac_old; ind_fac++) {
2959
2960 ind_val_deb = table_def_fac->pos[ind_fac ] - 1;
2961 ind_val_fin = table_def_fac->pos[ind_fac + 1] - 1;
2962
2963 if (ind_val_fin - ind_val_deb > 2) {
2964
2965 for (ind_val = ind_val_deb; ind_val < ind_val_fin; ind_val++)
2966 table_def_fac->val[cpt_val++] = table_def_fac->val[ind_val];
2967
2968 tab_fac_old_new.val[ind_fac] = 1 + cpt_fac++;
2969
2970 table_def_fac->pos[cpt_fac] = cpt_val + 1;
2971
2972 }
2973 else
2974 tab_fac_old_new.val[ind_fac] = 0;
2975
2976 }
2977
2978 nbr_fac_new = cpt_fac;
2979
2980 /* Si on n'a pas de faces dégénérées, on peut sortir */
2981
2982 if (nbr_fac_new == nbr_fac_old) {
2983
2984 tab_fac_old_new.nbr = 0;
2985 ECS_FREE(tab_fac_old_new.val);
2986
2987 ecs_table__pos_en_regle(table_def_fac);
2988
2989 return tab_fac_old_new;
2990 }
2991
2992 printf(_("\nMesh verification:\n\n"
2993 " Removal of %d degenerate faces:\n"
2994 " Initial number of faces : %10d\n"
2995 " Number of faces after processing : %10d\n"),
2996 (int)(nbr_fac_old - nbr_fac_new), (int)nbr_fac_old, (int)nbr_fac_new);
2997
2998 /* On redimensionne le tableau des faces */
2999
3000 table_def_fac->nbr = nbr_fac_new;
3001 ECS_REALLOC(table_def_fac->pos, nbr_fac_new + 1, ecs_size_t);
3002 ECS_REALLOC(table_def_fac->val, cpt_val, ecs_int_t);
3003
3004 ecs_table__pos_en_regle(table_def_fac);
3005
3006 return tab_fac_old_new;
3007 }
3008
3009 /*----------------------------------------------------------------------------
3010 * Fonction qui renvoie un tableau associant un type à chaque face, sous
3011 * forme de masque : 0 pour face isolée, 1 ou 2 pour face de bord (1 si
3012 * cellule avec cette face normale sortante, 2 si cellule avec cette face
3013 * normale entrante), 1+2 = 3 pour face interne, et 4 ou plus pour tous
3014 * les autres cas, correspondant à une erreur de connectivité (+4 pour faces
3015 * voyant au moins deux cellules avec face normale sortante, +8 pour faces
3016 * voyant au moins deux cellules avec face normale entrante).
3017 *
3018 * Le type de chaque face pourra être modifié ultérieurement en fonction
3019 * des informations de périodicité.
3020 *----------------------------------------------------------------------------*/
3021
3022 ecs_tab_int_t
ecs_table_def__typ_fac_cel(ecs_table_t * table_def_cel,ecs_table_t * table_def_fac)3023 ecs_table_def__typ_fac_cel(ecs_table_t *table_def_cel,
3024 ecs_table_t *table_def_fac)
3025 {
3026 size_t nbr_cel;
3027 size_t nbr_fac;
3028 ecs_int_t num_fac;
3029 size_t icel;
3030 size_t ifac;
3031 size_t ipos;
3032
3033 ecs_tab_int_t typ_fac;
3034
3035 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
3036
3037 typ_fac.nbr = 0;
3038 typ_fac.val = NULL;
3039
3040 if (table_def_cel == NULL)
3041 return typ_fac;
3042
3043 /* Initialisations */
3044
3045 ecs_table__regle_en_pos(table_def_cel);
3046
3047 nbr_cel = table_def_cel->nbr;
3048 nbr_fac = table_def_fac->nbr;
3049
3050 /* Type de face selon le nombre de cellules voisines : 0 pour face isolée,
3051 1 ou 2 pour face de bord, 3 pour faces internes, et 4 pour tous les
3052 autres cas (faces voyant au moins deux cellules sur un même côté) */
3053
3054 typ_fac.nbr = nbr_fac;
3055
3056 ECS_MALLOC(typ_fac.val, nbr_fac, ecs_int_t);
3057
3058 for (ifac = 0; ifac < nbr_fac; ifac++)
3059 typ_fac.val[ifac] = 0;
3060
3061 /* Boucle sur les cellules : marquage */
3062 /*------------------------------------*/
3063
3064 for (icel = 0; icel < nbr_cel; icel++ ) {
3065
3066 for (ipos = table_def_cel->pos[icel] - 1;
3067 ipos < table_def_cel->pos[icel+1] - 1;
3068 ipos++) {
3069
3070 num_fac = table_def_cel->val[ipos];
3071
3072 ifac = ECS_ABS(num_fac) - 1;
3073
3074 if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 0)
3075 typ_fac.val[ifac] += 1;
3076
3077 else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 0)
3078 typ_fac.val[ifac] += 2;
3079
3080 else {
3081 if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 1)
3082 typ_fac.val[ifac] = typ_fac.val[ifac] | 4;
3083 else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 2)
3084 typ_fac.val[ifac] = typ_fac.val[ifac] | 8;
3085 }
3086
3087 }
3088
3089 }
3090
3091 ecs_table__libere_pos(table_def_cel);
3092
3093 return typ_fac;
3094 }
3095
3096 /*----------------------------------------------------------------------------
3097 * Fonction qui renvoie un tableau associant un type à chaque face les
3098 * numéros des cellules définies par cette face (normale sortante,
3099 * puis normale entrante). On affecte une valeur 0 lorsqu'il n'y a pas de
3100 * cellule correspondante directe (la périodicité n'est donc pas prise en
3101 * compte à ce niveau).
3102 *
3103 * On suppose que la cohérence du maillage a déjà été véridifiée et
3104 * qu'aucune face n'appartient à plus d'une cellule par côté.
3105 *----------------------------------------------------------------------------*/
3106
3107 ecs_tab_int_t
ecs_table_def__fac_cel(ecs_table_t * table_def_cel,ecs_table_t * table_def_fac)3108 ecs_table_def__fac_cel(ecs_table_t *table_def_cel,
3109 ecs_table_t *table_def_fac)
3110 {
3111 size_t nbr_cel;
3112 size_t nbr_fac;
3113 ecs_int_t num_fac;
3114 size_t icel;
3115 size_t ifac;
3116 size_t ipos;
3117
3118 ecs_tab_int_t fac_cel;
3119
3120 /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
3121
3122 fac_cel.nbr = 0;
3123 fac_cel.val = NULL;
3124
3125 if (table_def_cel == NULL)
3126 return fac_cel;
3127
3128 /* Initialisations */
3129
3130 ecs_table__regle_en_pos(table_def_cel);
3131
3132 nbr_cel = table_def_cel->nbr;
3133 nbr_fac = table_def_fac->nbr;
3134
3135 /* Allocation et mise à zéro des connectivités */
3136
3137 fac_cel.nbr = nbr_fac * 2;
3138
3139 ECS_MALLOC(fac_cel.val, fac_cel.nbr, ecs_int_t);
3140
3141 for (ipos = 0; ipos < fac_cel.nbr; ipos++)
3142 fac_cel.val[ipos] = 0;
3143
3144 /* Boucle sur les cellules : marquage */
3145 /*------------------------------------*/
3146
3147 for (icel = 0; icel < nbr_cel; icel++ ) {
3148
3149 for (ipos = table_def_cel->pos[icel] - 1;
3150 ipos < table_def_cel->pos[icel+1] - 1;
3151 ipos++) {
3152
3153 num_fac = table_def_cel->val[ipos];
3154
3155 ifac = ECS_ABS(num_fac) - 1;
3156
3157 if (num_fac > 0) {
3158 assert(fac_cel.val[ifac*2] == 0);
3159 fac_cel.val[ifac*2] = icel + 1;
3160 }
3161 else {
3162 assert(fac_cel.val[ifac*2 + 1] == 0);
3163 fac_cel.val[ifac*2 + 1] = icel + 1;
3164 }
3165
3166 }
3167
3168 }
3169
3170 ecs_table__libere_pos(table_def_cel);
3171
3172 return fac_cel;
3173 }
3174
3175 /*----------------------------------------------------------------------------*/
3176
3177