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
*============================================================================*/
8
This file is part of Code_Saturne, a general-purpose CFD tool.
10
Copyright (C) 1998-2011 EDF S.A.
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
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
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.
27
/*----------------------------------------------------------------------------*/
30
/*============================================================================
32
*============================================================================*/
34
/*----------------------------------------------------------------------------
35
* Fichiers `include' librairie standard C
36
*----------------------------------------------------------------------------*/
43
/*----------------------------------------------------------------------------
44
* Fichiers `include' visibles du paquetage global "Utilitaire"
45
*----------------------------------------------------------------------------*/
47
#include "ecs_elt_typ_liste.h"
53
/*----------------------------------------------------------------------------
54
* Fichiers `include' visibles des paquetages visibles
55
*----------------------------------------------------------------------------*/
58
/*----------------------------------------------------------------------------
59
* Fichiers `include' visibles du paquetage courant
60
*----------------------------------------------------------------------------*/
62
#include "ecs_table.h"
65
/*----------------------------------------------------------------------------
66
* Fichier `include' du paquetage courant associe au fichier courant
67
*----------------------------------------------------------------------------*/
69
#include "ecs_table_def.h"
72
/*----------------------------------------------------------------------------
73
* Fichiers `include' prives du paquetage courant
74
*----------------------------------------------------------------------------*/
76
#include "ecs_table_priv.h"
79
/*============================================================================
80
* Macros globales au fichier
81
*============================================================================*/
83
/* Longueur maximale du nom d'un type d'élément (+1 pour le `\0' !) */
84
#define ECS_LOC_LNG_MAX_NOM_TYP 11
87
#define FLT_MAX HUGE_VAL
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] )
95
#define ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) ( \
96
vect1[0] * vect2[0] + vect1[1] * vect2[1] + vect1[2] * vect2[2] )
98
#define ECS_LOC_MODULE(vect) \
99
sqrt(vect[0] * vect[0] + vect[1] * vect[1] + vect[2] * vect[2])
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]) )
106
/*============================================================================
108
*============================================================================*/
110
/*----------------------------------------------------------------------------
111
* Fonction qui met à jour la définition faces -> sommets en cas
112
* de fusion de sommets.
113
*----------------------------------------------------------------------------*/
116
_table_def__maj_fac_som(ecs_table_t *table_def_fac,
117
const ecs_tab_int_t *tab_som_old_new)
122
size_t nbr_val_ini, nbr_val_fin;
125
size_t num_som, num_som_prev;
131
size_t ipos_deb, ipos_fin;
133
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
135
/* Initialisations */
136
/* --------------- */
138
nbr_fac = table_def_fac->nbr;
140
nbr_val_ini = table_def_fac->pos[nbr_fac];
142
/* Mise à jour de la définition des faces */
143
/*----------------------------------------*/
146
ind_som < table_def_fac->pos[table_def_fac->nbr] - 1;
149
table_def_fac->val[ind_som]
150
= tab_som_old_new->val[table_def_fac->val[ind_som] - 1];
152
/* Suppression de sommets confondus de la définition des faces */
153
/*-------------------------------------------------------------*/
160
for (ind_fac = 0; ind_fac < nbr_fac; ind_fac++) {
164
ipos_fin = table_def_fac->pos[ind_fac + 1] - 1;
166
nbr_som_fac = ipos_fin - ipos_deb;
168
num_som_prev = table_def_fac->val[ipos_deb + nbr_som_fac - 1];
170
for (ind_som = 0; ind_som < nbr_som_fac; ind_som++) {
172
num_som = table_def_fac->val[ipos_deb + ind_som];
174
if (num_som != num_som_prev) {
175
num_som_prev = num_som;
176
table_def_fac->val[cpt_val++] = num_som;
182
table_def_fac->pos[ind_fac + 1] = cpt_val + 1;
186
nbr_fac_mod += ind_fac_mod;
189
nbr_val_fin = table_def_fac->pos[nbr_fac];
191
assert(nbr_val_fin <= nbr_val_ini);
193
if (nbr_val_fin != nbr_val_ini) {
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"),
202
/*----------------------------------------------------------------------------
203
* Fonction qui transforme un tableau d'équivalence en liste chaînée simple
204
*----------------------------------------------------------------------------*/
207
_table_def__transf_equiv(size_t nbr_som,
208
ecs_tab_int_t *tab_equiv_som)
211
ecs_int_t ind_som_min;
212
ecs_int_t ind_som_max;
213
ecs_int_t ind_som_tmp;
215
ecs_tab_int_t tab_equiv_prec;
217
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
219
tab_equiv_prec.nbr = nbr_som;
220
ECS_MALLOC(tab_equiv_prec.val, tab_equiv_prec.nbr, ecs_int_t);
222
for (ind_som = 0; ind_som < nbr_som; ind_som++)
223
tab_equiv_prec.val[ind_som] = -1;
225
for (ind_som = 0; ind_som < nbr_som; ind_som++) {
227
if (tab_equiv_som->val[ind_som] != -1) {
229
ind_som_min = ind_som;
230
ind_som_max = tab_equiv_som->val[ind_som];
232
assert (ind_som_min < ind_som_max);
234
while ( ind_som_min != ind_som_max
235
&& tab_equiv_prec.val[ind_som_max] != ind_som_min) {
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).
243
while (tab_equiv_prec.val[ind_som_max] > ind_som_min)
244
ind_som_max = tab_equiv_prec.val[ind_som_max];
246
ind_som_tmp = tab_equiv_prec.val[ind_som_max];
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.
254
tab_equiv_prec.val[ind_som_max] = ind_som_min;
256
if (ind_som_tmp != -1) {
257
ind_som_max = ind_som_min;
258
ind_som_min = ind_som_tmp;
265
for (ind_som = 0; ind_som < nbr_som; ind_som++) {
267
if (tab_equiv_prec.val[ind_som] != -1)
268
tab_equiv_som->val[tab_equiv_prec.val[ind_som]] = ind_som;
272
tab_equiv_prec.nbr = 0;
273
ECS_FREE(tab_equiv_prec.val);
276
/*----------------------------------------------------------------------------
277
* Fonction qui fusionne des sommets équivalents
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
*----------------------------------------------------------------------------*/
288
_table_def__fusion_som(size_t *n_vertices,
289
ecs_coord_t **coords,
290
ecs_tab_int_t *tab_equiv_som)
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;
297
ecs_tab_int_t tab_som_old_new;
302
nbr_som_old = *n_vertices;
304
printf(_("\n Merging vertices:\n"));
306
/* Transform vertex equivalences array into simple linked list. */
308
_table_def__transf_equiv(nbr_som_old, tab_equiv_som);
310
printf(_(" Initial number of vertices : %10d\n"), (int)nbr_som_old);
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);
315
/* Initialize vertex renumbering array */
317
for (ind_som = 0; ind_som < nbr_som_old; ind_som++)
318
tab_som_old_new.val[ind_som] = 0;
320
/* Main loop on vertices */
321
/*-----------------------*/
325
for (ind_som = 0; ind_som < nbr_som_old; ind_som++) {
327
/* If the vertex has not been handled yet */
329
if (tab_som_old_new.val[ind_som] == 0) {
331
/* Initialize vertex */
335
for (icoo = 0; icoo < 3; icoo++)
336
som_coord[icoo] = 0.0;
338
ind_som_loc = ind_som;
340
/* Mark equivalent vertices and compute their contribution */
344
tab_som_old_new.val[ind_som_loc] = nbr_som_new + 1;
348
for (icoo = 0; icoo < 3; icoo++)
349
som_coord[icoo] += (* coords)[3*ind_som_loc + icoo];
351
ind_som_loc = tab_equiv_som->val[ind_som_loc];
353
} while (ind_som_loc != -1);
355
/* Final coordinates */
357
for (icoo = 0; icoo < 3; icoo++)
358
(*coords)[3*nbr_som_new + icoo] = som_coord[icoo] / som_poids;
360
/* Do not forget to increment the number of vertices after merge */
367
/* End of main loop */
368
/* ---------------- */
370
/* We now know the number of vertices after merging */
372
*n_vertices = nbr_som_new;
373
ECS_REALLOC(*coords, nbr_som_new*3, ecs_coord_t);
375
printf(_(" Number of vertices after merging : %10d\n"), (int)nbr_som_new);
377
/* Mark vertices originating from a merge */
378
/*----------------------------------------*/
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.
388
for (ind_som = 0; ind_som < tab_equiv_som->nbr; ind_som++) {
390
if (tab_equiv_som->val[ind_som] != -1) {
394
ind_som_loc = tab_equiv_som->val[ind_som];
396
while (ind_som_loc != -1) {
398
ind_som_tmp = ind_som_loc;
399
ind_som_loc = tab_equiv_som->val[ind_som_loc];
401
tab_equiv_som->val[ind_som_tmp] = -1;
407
printf(_(" Number of modified vertices : %10d\n"), (int)nbr_som_fus);
411
return tab_som_old_new;
414
/*----------------------------------------------------------------------------
415
* Fonction qui met à jour le tableau des fusions de sommets
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
*----------------------------------------------------------------------------*/
426
_table_def__maj_equiv_som(size_t ind_som_0,
428
ecs_tab_int_t *tab_equiv_som)
430
ecs_int_t ind_som_max;
431
ecs_int_t ind_som_min;
432
ecs_int_t ind_som_tmp;
434
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
436
ind_som_min = ECS_MIN(ind_som_0, ind_som_1);
437
ind_som_max = ECS_MAX(ind_som_0, ind_som_1);
439
/* On fusionne deux listes chaînées */
440
/*----------------------------------*/
442
while ( ind_som_max != ind_som_min
443
&& tab_equiv_som->val[ind_som_min] != ind_som_max) {
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).
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];
455
ind_som_tmp = tab_equiv_som->val[ind_som_min];
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
463
tab_equiv_som->val[ind_som_min] = ind_som_max;
465
if (ind_som_tmp != -1) {
466
ind_som_min = ind_som_max;
467
ind_som_max = ind_som_tmp;
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
*----------------------------------------------------------------------------*/
481
_orient_quad(const ecs_coord_t coord[],
487
ecs_int_t nb_angle_obtus;
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];
494
ecs_int_t passage = -1;
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] )
504
/* Calcul d'une direction normale approchée de la face
505
(peut être entrante ou sortante si la face est "croisée") */
507
ECS_LOC_INIT_VECT(vect1, 1, 2);
508
ECS_LOC_INIT_VECT(vect2, 1, 4);
510
ECS_LOC_PRODUIT_VECTORIEL(normale, vect1, vect2);
512
/* Boucle sur les renumérotations possibles */
530
isom_tmp = connect[2];
531
connect[2] = connect[3];
532
connect[3] = isom_tmp;
540
/* Boucle sur les coins */
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
549
for (itri = 0; itri < 4; itri++) {
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);
554
ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2);
556
/* Angle obtus si produit mixte < 0, aigu sinon. */
558
sgn = ECS_LOC_PRODUIT_SCALAIRE(prod_vect, normale);
565
} while (nb_angle_obtus == 2);
567
return (ECS_MIN(passage, 1));
569
#undef ECS_LOC_INIT_VECT
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
*----------------------------------------------------------------------------*/
580
_orient_tetra(const ecs_coord_t coord[],
585
ecs_coord_t vect12[3];
586
ecs_coord_t vect13[3];
587
ecs_coord_t vect14[3];
589
ecs_int_t passage = -1;
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] )
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 )
604
/* Boucle sur les renumérotations possibles */
618
ECS_LOC_PERMUTE(2, 3);
621
default: /* Retour connectivité d'origine et sortie */
622
ECS_LOC_PERMUTE(2, 3);
627
ECS_LOC_INIT_VECT(vect12, 1, 2);
628
ECS_LOC_INIT_VECT(vect13, 1, 3);
629
ECS_LOC_INIT_VECT(vect14, 1, 4);
631
det = ECS_LOC_DETERMINANT(vect12, vect13, vect14);
635
return (ECS_MIN(passage, 1));
637
#undef ECS_LOC_INIT_VECT
638
#undef ECS_LOC_PERMUTE
641
/*----------------------------------------------------------------------------
642
* Correction si nécessaire de l'orientation d'une pyramide en connectivité
643
* nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
644
* ou que l'on demande une simple vérification (i.e. correc = false),
645
* 0 si l'orientation initiale est bonne, et 1 en cas de correction de
646
* l'orientation par une permutation de la connectivité locale.
648
* Tous les cas ne sont pas détectés ou traités : on suppose que le
649
* sommet est toujours en position 5, mais que la base peut être
650
* parcourue dans l'ordre 1 4 3 2, 1 2 4 3, ou 1 3 4 2 au lieu de 1 2 3 4,
651
* dans quel cas on la réordonne.
652
*----------------------------------------------------------------------------*/
655
_orient_pyram(const ecs_coord_t coord[],
661
ecs_int_t connect_tmp[8];
664
ecs_coord_t vect1[3];
665
ecs_coord_t vect2[3];
666
ecs_coord_t vect3[3];
668
ecs_int_t retval = 0;
669
ecs_int_t passage = -1;
671
#define ECS_LOC_INIT_VECT(vect, i, j) ( \
672
vect[0] = coord[((connect_tmp[j-1] - 1) * 3) ] \
673
- coord[((connect_tmp[i-1] - 1) * 3) ], \
674
vect[1] = coord[((connect_tmp[j-1] - 1) * 3) + 1] \
675
- coord[((connect_tmp[i-1] - 1) * 3) + 1], \
676
vect[2] = coord[((connect_tmp[j-1] - 1) * 3) + 2] \
677
- coord[((connect_tmp[i-1] - 1) * 3) + 2] )
679
#define ECS_LOC_PERMUTE(i, j) ( \
680
isom_tmp = connect[i-1], \
681
connect_tmp[i-1] = connect_tmp[j-1], \
682
connect_tmp[j-1] = isom_tmp )
685
for (isom = 0; isom < 5; isom++)
686
connect_tmp[isom] = connect[isom];
688
/* Vérification et correction éventuelle de l'orientation de la base */
690
ret_base = _orient_quad(coord,
694
retval = ECS_MAX(ret_base, retval);
696
if ((correc == false && ret_base != 0) || ret_base < 0)
700
/* Boucle sur les renumérotations possibles */
718
ECS_LOC_PERMUTE(2, 4);
721
default: /* Retour connectivité d'origine et sortie */
722
ECS_LOC_PERMUTE(2, 4);
727
ECS_LOC_INIT_VECT(vect1, 1, 2);
728
ECS_LOC_INIT_VECT(vect2, 1, 4);
729
ECS_LOC_INIT_VECT(vect3, 1, 5);
731
det = ECS_LOC_DETERMINANT(vect1, vect2, vect3);
736
for (isom = 0; isom < 5; isom++)
737
connect[isom] = connect_tmp[isom];
742
#undef ECS_LOC_INIT_VECT
743
#undef ECS_LOC_PERMUTE
746
/*----------------------------------------------------------------------------
747
* Correction si nécessaire de l'orientation d'un prisme en connectivité
748
* nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
749
* ou que l'on demande une simple vérification (i.e. correc = false),
750
* 0 si l'orientation initiale est bonne, et 1 en cas de correction de
751
* l'orientation par une permutation de la connectivité locale.
753
* Tous les cas ne sont pas détectés ou traités : on suppose que les
754
* bases peuvent être parcourues dans l'ordre 1 3 2 et 4 6 5 au lieu
755
* de 1 2 3 et 4 5 6, dans quel cas on les réordonne.
756
*----------------------------------------------------------------------------*/
759
_orient_prism(const ecs_coord_t coord[],
766
ecs_coord_t vect1[3];
767
ecs_coord_t vect2[3];
768
ecs_coord_t vect3[3];
770
ecs_int_t passage = -1;
772
#define ECS_LOC_INIT_VECT(vect, i, j) ( \
773
vect[0] = coord[((connect[j-1] - 1) * 3) ] \
774
- coord[((connect[i-1] - 1) * 3) ], \
775
vect[1] = coord[((connect[j-1] - 1) * 3) + 1] \
776
- coord[((connect[i-1] - 1) * 3) + 1], \
777
vect[2] = coord[((connect[j-1] - 1) * 3) + 2] \
778
- coord[((connect[i-1] - 1) * 3) + 2] )
780
#define ECS_LOC_PERMUTE(i, j) ( \
781
isom_tmp = connect[i-1], \
782
connect[i-1] = connect[j-1], \
783
connect[j-1] = isom_tmp )
785
/* Boucle sur les renumérotations possibles */
801
ECS_LOC_PERMUTE(2, 3);
802
ECS_LOC_PERMUTE(5, 6);
805
default: /* Retour connectivité d'origine et sortie */
806
ECS_LOC_PERMUTE(2, 3);
807
ECS_LOC_PERMUTE(5, 6);
812
ECS_LOC_INIT_VECT(vect2, 1, 2);
813
ECS_LOC_INIT_VECT(vect3, 1, 3);
815
ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3);
817
for (idim = 0; idim < 3; idim++)
818
vect2[idim] = ( coord[((connect[4-1] - 1) * 3) + idim]
819
+ coord[((connect[5-1] - 1) * 3) + idim]
820
+ coord[((connect[6-1] - 1) * 3) + idim]
821
- coord[((connect[1-1] - 1) * 3) + idim]
822
- coord[((connect[2-1] - 1) * 3) + idim]
823
- coord[((connect[3-1] - 1) * 3) + idim]);
825
pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2);
827
} while (pscal < 0.0);
829
return (ECS_MIN(passage, 1));
831
#undef ECS_LOC_INIT_VECT
832
#undef ECS_LOC_PERMUTE
835
/*----------------------------------------------------------------------------
836
* Correction si nécessaire de l'orientation d'un hexaèdre en connectivité
837
* nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
838
* ou que l'on demande une simple vérification (i.e. correc = false),
839
* 0 si l'orientation initiale est bonne, et 1 en cas de correction de
840
* l'orientation par une permutation de la connectivité locale.
842
* Tous les cas ne sont pas détectés ou traités : on suppose que les
843
* bases peuvent être parcourues dans l'ordre soit 1 4 3 2 et 5 8 7 6,
844
* soit 1 2 4 3 et 5 6 8 7, soit 1 3 4 2 et 5 7 8 6 au lieu
845
* de 1 2 3 4 et 5 6 7 8, dans quel cas on les réordonne.
846
*----------------------------------------------------------------------------*/
849
_orient_hexa(const ecs_coord_t coord[],
856
ecs_int_t connect_tmp[8];
857
ecs_int_t ret_base_1;
858
ecs_int_t ret_base_2;
860
ecs_coord_t vect1[3];
861
ecs_coord_t vect2[3];
862
ecs_coord_t vect3[3];
864
ecs_int_t retval = 0;
865
ecs_int_t passage = -1;
867
#define ECS_LOC_INIT_VECT(vect, i, j) ( \
868
vect[0] = coord[((connect_tmp[j-1] - 1) * 3) ] \
869
- coord[((connect_tmp[i-1] - 1) * 3) ], \
870
vect[1] = coord[((connect_tmp[j-1] - 1) * 3) + 1] \
871
- coord[((connect_tmp[i-1] - 1) * 3) + 1], \
872
vect[2] = coord[((connect_tmp[j-1] - 1) * 3) + 2] \
873
- coord[((connect_tmp[i-1] - 1) * 3) + 2] )
875
#define ECS_LOC_PERMUTE(i, j) ( \
876
isom_tmp = connect[i-1], \
877
connect_tmp[i-1] = connect_tmp[j-1], \
878
connect_tmp[j-1] = isom_tmp )
881
for (isom = 0; isom < 8; isom++)
882
connect_tmp[isom] = connect[isom];
884
/* Vérification et correction éventuelle de l'orientation des bases */
886
ret_base_1 = _orient_quad(coord,
890
if ((correc == false && ret_base_1 != 0) || ret_base_1 < 0)
893
else if (ret_base_1 > 0) {
894
ret_base_2 = _orient_quad(coord,
897
if (ret_base_2 != ret_base_1)
903
/* Boucle sur les renumérotations possibles */
921
ECS_LOC_PERMUTE(2, 4);
922
ECS_LOC_PERMUTE(6, 8);
931
ECS_LOC_INIT_VECT(vect2, 1, 2);
932
ECS_LOC_INIT_VECT(vect3, 1, 4);
934
ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3);
936
for (idim = 0; idim < 3; idim++)
937
vect2[idim] = ( coord[((connect_tmp[5-1] - 1) * 3) + idim]
938
+ coord[((connect_tmp[6-1] - 1) * 3) + idim]
939
+ coord[((connect_tmp[7-1] - 1) * 3) + idim]
940
+ coord[((connect_tmp[8-1] - 1) * 3) + idim]
941
- coord[((connect_tmp[1-1] - 1) * 3) + idim]
942
- coord[((connect_tmp[2-1] - 1) * 3) + idim]
943
- coord[((connect_tmp[3-1] - 1) * 3) + idim]
944
- coord[((connect_tmp[4-1] - 1) * 3) + idim]) * 0.25;
946
pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2);
948
} while (pscal < 0.0);
951
for (isom = 0; isom < 8; isom++)
952
connect[isom] = connect_tmp[isom];
955
/* Vérification et correction éventuelle de l'orientation des cotés */
957
connect_tmp[1-1] = connect[2-1];
958
connect_tmp[2-1] = connect[3-1];
959
connect_tmp[3-1] = connect[7-1];
960
connect_tmp[4-1] = connect[6-1];
962
ret_base_1 = _orient_quad(coord,
969
connect_tmp[1-1] = connect[1-1];
970
connect_tmp[2-1] = connect[2-1];
971
connect_tmp[3-1] = connect[6-1];
972
connect_tmp[4-1] = connect[5-1];
974
ret_base_1 = _orient_quad(coord,
983
#undef ECS_LOC_INIT_VECT
984
#undef ECS_LOC_PERMUTE
987
/*----------------------------------------------------------------------------
988
* Compute contribution of a polygonal face to a cell's volume,
989
* using Stoke's theorem
992
* *---------* B : barycentre of the polygon
994
* / . . \ Pi : vertices of the polygon
996
* / . . Ti \ Ti : triangle
997
* *.........B.........* Pi
1004
*----------------------------------------------------------------------------*/
1007
_compute_face_vol_contrib(ecs_int_t n_face_vertices,
1008
const ecs_coord_t vtx_coord[],
1009
const ecs_int_t face_vtx_lst[])
1011
ecs_int_t i, tri_id, vtx_id;
1012
ecs_coord_t tri_center[3], tri_norm[3];
1013
ecs_coord_t vect1[3], vect2[3];
1015
ecs_coord_t face_barycentre[3] = {0., 0., 0.};
1017
double inv3 = 1./3.;
1019
double vol_contrib = 0.;
1021
/* Compute barycentre (B) coordinates for the polygon (P) */
1023
for (vtx_id = 0; vtx_id < n_face_vertices; vtx_id++) {
1024
size_t coord_id = face_vtx_lst[vtx_id] - 1;
1025
for (i = 0; i < 3; i++)
1026
face_barycentre[i] += vtx_coord[coord_id*3 + i];
1029
for (i = 0; i < 3; i++)
1030
face_barycentre[i] /= n_face_vertices;
1032
/* Loop on triangles of the face */
1034
for (tri_id = 0; tri_id < n_face_vertices; tri_id++) {
1036
size_t coord_id_1 = face_vtx_lst[tri_id % n_face_vertices] - 1;
1037
size_t coord_id_2 = face_vtx_lst[(tri_id + 1)% n_face_vertices] - 1;
1039
for (i = 0; i < 3; i++) {
1040
vect1[i] = vtx_coord[coord_id_1*3 + i] - face_barycentre[i];
1041
vect2[i] = vtx_coord[coord_id_2*3 + i] - face_barycentre[i];
1044
ECS_LOC_PRODUIT_VECTORIEL(tri_norm, vect1, vect2);
1046
for (i = 0; i < 3; i++)
1049
/* Computation of the center of the triangle Ti */
1051
for (i = 0; i < 3; i++)
1052
tri_center[i] = ( face_barycentre[i]
1053
+ vtx_coord[coord_id_1*3 + i]
1054
+ vtx_coord[coord_id_2*3 + i]) * inv3;
1056
/* Contribution to cell volume using Stoke's formula */
1058
for (i = 0; i < 3; i++)
1059
vol_contrib += (tri_norm[i] * tri_center[i]);
1061
} /* End of loop on triangles of the face */
1066
/*----------------------------------------------------------------------------
1067
* Correction si nécessaire de l'orientation d'un polyedre en connectivité
1068
* nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
1069
* ou que l'on demande une simple vérification (i.e. correc = false),
1070
* 0 si l'orientation initiale est bonne, et 1 en cas de correction de
1071
* l'orientation par une permutation de la connectivité locale.
1073
* On suppose que toutes les faces polygonales sont bien définies, mais
1074
* que certaines faces peuvent être définies avec une normale intérieure
1075
* à la cellule et non pas extérieure. On suppose que la non-convexité
1076
* du polyèdre est limitée.
1077
*----------------------------------------------------------------------------*/
1080
_orient_polyhedron(const ecs_coord_t coord[],
1081
ecs_int_t connect[],
1084
ecs_tab_int_t *face_index,
1085
ecs_tab_int_t *face_marker,
1086
ecs_tab_int_t *edges)
1088
size_t i, j, face_id, edge_id;
1089
ecs_int_t ind_premier, ind_dernier;
1091
double cell_vol = 0.;
1092
size_t n_unmarked_faces = 0;
1095
ecs_int_t retval = 0;
1097
/* Ensure working arrays are large enough */
1102
ecs_int_t premier_som = -1;
1104
for (i = 0; i < (size_t)size; i++) {
1106
if (face_index->nbr < n_faces + 1) {
1107
face_index->nbr = ECS_MAX(n_faces + 1, face_index->nbr*2);
1108
ECS_REALLOC(face_index->val, face_index->nbr, ecs_int_t);
1111
if (premier_som == -1) {
1112
face_index->val[n_faces] = i;
1113
premier_som = connect[i];
1116
else if (connect[i] == premier_som) {
1123
if (face_index->nbr < n_faces + 1) {
1124
face_index->nbr = ECS_MAX(n_faces + 1, face_index->nbr*2);
1125
ECS_REALLOC(face_index->val, face_index->nbr, ecs_int_t);
1128
face_index->val[n_faces] = size;
1130
/* face_marker: 0 initially, 1 if oriented, -1 if inverted */
1132
if (face_marker->nbr < n_faces) {
1133
face_marker->nbr = ECS_MAX(n_faces, face_marker->nbr*2);
1134
ECS_REALLOC(face_marker->val, face_marker->nbr, ecs_int_t);
1137
for (face_id = 0; face_id < n_faces; face_id++)
1138
face_marker->val[face_id] = 0;
1144
/* Extract edges and build edge-face connectivity */
1145
/*------------------------------------------------*/
1147
for (face_id = 0; face_id < n_faces; face_id++) {
1149
size_t start_id = face_index->val[face_id];
1150
size_t end_id = face_index->val[face_id + 1] - 1;
1152
/* Loop on unassociated edges */
1154
for (i = start_id; i < end_id; i++) {
1156
ecs_int_t v_id_1 = connect[i], v_id_2 = connect[i + 1];
1157
ecs_int_t is_match = 0;
1159
/* Compare with other edges */
1161
for (edge_id = 0; edge_id < n_edges; edge_id++) {
1163
ecs_int_t v_id_1_cmp = edges->val[edge_id*4];
1164
ecs_int_t v_id_2_cmp = edges->val[edge_id*4 + 1];
1166
if ((v_id_1 == v_id_2_cmp) && (v_id_2 == v_id_1_cmp))
1169
else if ((v_id_1 == v_id_1_cmp) && (v_id_2 == v_id_2_cmp))
1172
/* Each edge should be shared by exactly 2 faces */
1174
if (is_match != 0) {
1175
if (edges->val[edge_id*4 + 3] == 0) {
1176
edges->val[edge_id*4 + 3] = is_match * (face_id + 1);
1184
/* If edge is unmatched, add it */
1186
if (is_match == 0) {
1188
if (edges->nbr < (n_edges + 1)*4) {
1189
edges->nbr = ECS_MAX(edges->nbr*2, (n_edges + 1)*4);
1190
ECS_REALLOC(edges->val, edges->nbr, ecs_int_t);
1193
edges->val[n_edges*4] = v_id_1;
1194
edges->val[n_edges*4 + 1] = v_id_2;
1195
edges->val[n_edges*4 + 2] = face_id + 1;
1196
edges->val[n_edges*4 + 3] = 0;
1205
/* Check if each edge is associated with 2 faces
1206
(i.e., that cell is closed) */
1208
for (edge_id = 0; edge_id < n_edges; edge_id++) {
1209
if (edges->val[edge_id*4 + 3] == 0)
1213
/* Now check if face orientations are compatible */
1214
/*-----------------------------------------------*/
1216
face_marker->val[0] = 1; /* First face defines reference */
1217
n_unmarked_faces = n_faces - 1;
1219
while (n_unmarked_faces > 0) {
1221
for (edge_id = 0; edge_id < n_edges; edge_id++) {
1223
ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
1224
ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
1225
ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;
1227
if ( face_marker->val[face_id_1] == 0
1228
&& face_marker->val[face_id_2] != 0) {
1230
face_marker->val[face_id_1] = face_marker->val[face_id_2];
1232
face_marker->val[face_id_1] = - face_marker->val[face_id_2];
1233
n_unmarked_faces -= 1;
1236
else if ( face_marker->val[face_id_1] != 0
1237
&& face_marker->val[face_id_2] == 0) {
1239
face_marker->val[face_id_2] = face_marker->val[face_id_1];
1241
face_marker->val[face_id_2] = - face_marker->val[face_id_1];
1242
n_unmarked_faces -= 1;
1248
for (edge_id = 0; edge_id < n_edges; edge_id++) {
1250
ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
1251
ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
1252
ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;
1254
ecs_int_t marker_product
1255
= face_marker->val[face_id_1]*face_marker->val[face_id_2];
1257
if ( (face_num_2 < 0 && marker_product > 0)
1258
|| (face_num_2 > 0 && marker_product < 0))
1263
/* At this stage, topology is correct; check for inside-out cell */
1264
/*---------------------------------------------------------------*/
1266
for (face_id = 0; face_id < n_faces; face_id++) {
1268
size_t start_id = face_index->val[face_id];
1269
size_t end_id = face_index->val[face_id + 1] - 1;
1272
= _compute_face_vol_contrib(end_id - start_id,
1274
connect + start_id);
1275
cell_vol += vol_contrib * (double)(face_marker->val[face_id]);
1278
/* Invert orientation if cell is inside_out */
1280
if (cell_vol < 0.) {
1281
for (face_id = 0; face_id < n_faces; face_id++)
1282
face_marker->val[face_id] *= -1;
1285
/* Now correct connectivity if required */
1286
/*--------------------------------------*/
1288
if (correc == true) {
1290
for (face_id = 0; face_id < n_faces; face_id++) {
1292
if (face_marker->val[face_id] == -1) {
1294
for (i = face_index->val[face_id] + 1,
1295
j = face_index->val[face_id + 1] - 2;
1298
ecs_int_t connect_tmp = connect[i];
1299
connect[i] = connect[j];
1300
connect[j] = connect_tmp;
1311
for (face_id = 0; face_id < n_faces; face_id++) {
1312
if (face_marker->val[face_id] == -1)
1320
/*============================================================================
1321
* Fonctions publiques
1322
*============================================================================*/
1324
/*----------------------------------------------------------------------------
1325
* Fonction qui réalise le tri des types géométriques
1326
* La fonction affiche le nombre d'éléments par type géométrique
1327
*----------------------------------------------------------------------------*/
1330
ecs_table_def__trie_typ(ecs_table_t *this_table_def,
1337
ecs_tab_int_t tab_typ_geo_ord;
1338
ecs_tab_int_t tab_typ_geo;
1339
ecs_tab_int_t vect_renum;
1341
ecs_elt_typ_t typ_geo;
1343
const ecs_elt_typ_t *typ_geo_base;
1345
const int lng_imp = ECS_LNG_AFF_STR - ECS_LOC_LNG_MAX_NOM_TYP;
1347
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1349
assert(this_table_def != NULL);
1351
vect_renum.nbr = 0 ;
1352
vect_renum.val = NULL;
1354
nbr_elt = this_table_def->nbr;
1359
/* Détermination de base du type d'élément "classique" */
1361
assert(dim_elt >=2 && dim_elt <= 3);
1363
typ_geo_base = ecs_glob_typ_elt[dim_elt - 2];
1365
/* Si tous les éléments sont de même type, rien à faire */
1366
/*------------------------------------------------------*/
1368
if (this_table_def->pos == NULL) {
1370
if (this_table_def->pas < 9)
1371
typ_geo = typ_geo_base[this_table_def->pas];
1372
else if (dim_elt == 2)
1373
typ_geo = ECS_ELT_TYP_FAC_POLY;
1374
else if (dim_elt == 3)
1375
typ_geo = ECS_ELT_TYP_CEL_POLY;
1377
typ_geo = ECS_ELT_TYP_NUL;
1379
printf(" %-*s%*.*s : %*lu\n",
1380
lng_imp, _("Number of elements"),
1381
ECS_LOC_LNG_MAX_NOM_TYP, ECS_LOC_LNG_MAX_NOM_TYP,
1382
ecs_fic_elt_typ_liste_c[typ_geo].nom,
1383
ECS_LNG_AFF_ENT, (unsigned long)nbr_elt);
1387
/* Si tous les éléments ne sont pas de même type */
1388
/*-----------------------------------------------*/
1392
/* Construction d'un tableau temporaire de type géométrique */
1394
tab_typ_geo.nbr = nbr_elt;
1395
ECS_MALLOC(tab_typ_geo.val, tab_typ_geo.nbr, ecs_int_t);
1397
for (ielt = 0; ielt < nbr_elt; ielt++) {
1399
nbr_som = this_table_def->pos[ielt + 1]
1400
- this_table_def->pos[ielt];
1403
tab_typ_geo.val[ielt] = typ_geo_base[nbr_som];
1405
else if (dim_elt == 2)
1406
tab_typ_geo.val[ielt] = ECS_ELT_TYP_FAC_POLY;
1408
else if (dim_elt == 3)
1409
tab_typ_geo.val[ielt] = ECS_ELT_TYP_CEL_POLY;
1412
tab_typ_geo.val[ielt] = ECS_ELT_TYP_NUL;
1416
/* On regarde si les types géométriques ne sont pas deja ordonnés */
1419
while (ielt < nbr_elt &&
1420
tab_typ_geo.val[ielt] >= tab_typ_geo.val[ielt - 1] )
1423
if (ielt < nbr_elt) {
1425
/* Les types géométriques ne sont pas ordonnés */
1426
/* On ordonne les types géométriques */
1428
vect_renum.nbr = nbr_elt;
1429
ECS_MALLOC(vect_renum.val, nbr_elt, ecs_int_t);
1432
tab_typ_geo_ord = ecs_tab_int__trie_et_renvoie(tab_typ_geo,
1436
`vect_renum' prend pour indice les indices nouveaux, et ses valeurs
1437
contiennent les indices anciens correspondants
1438
On inverse le contenu de `vect_renum' :
1439
à chaque indice ancien, `vect_renum' donne la valeur du nouvel indice
1442
ecs_tab_int__inverse(&vect_renum);
1444
ECS_FREE(tab_typ_geo.val);
1446
tab_typ_geo.val = tab_typ_geo_ord.val;
1448
tab_typ_geo_ord.nbr = 0;
1449
tab_typ_geo_ord.val = NULL;
1453
/* Message d'information sur la composition des éléments du maillage */
1456
ecs_int_t val_typ_ref;
1457
size_t nbr_val_typ_geo;
1459
size_t cpt_ielt = 0;
1461
while (cpt_ielt < nbr_elt) {
1463
val_typ_ref = tab_typ_geo.val[cpt_ielt];
1465
/* On compte le nombre d'éléments ayant le même type géométrique */
1467
nbr_val_typ_geo = 0;
1469
for (ielt = cpt_ielt;
1470
ielt < nbr_elt && tab_typ_geo.val[ielt] == val_typ_ref; ielt++)
1473
printf(" %-*s%*.*s : %*lu\n",
1474
lng_imp, _("Number of elements"),
1475
ECS_LOC_LNG_MAX_NOM_TYP, ECS_LOC_LNG_MAX_NOM_TYP,
1476
ecs_fic_elt_typ_liste_c[val_typ_ref].nom,
1477
ECS_LNG_AFF_ENT, (unsigned long)nbr_val_typ_geo);
1479
cpt_ielt += nbr_val_typ_geo;
1484
/* Le tableau de type géométrique n'est plus nécessaire */
1486
tab_typ_geo.nbr = 0;
1487
ECS_FREE(tab_typ_geo.val);
1491
/* Renvoi du vecteur de renumérotation */
1492
/*-------------------------------------*/
1497
/*----------------------------------------------------------------------------
1498
* Fonction qui construit
1499
* les définitions des faces par décomposition des tables des cellules
1500
*----------------------------------------------------------------------------*/
1503
ecs_table_def__decompose_cel(ecs_table_t **table_def_fac,
1504
ecs_table_t *table_def_cel)
1512
size_t nbr_val_fac_old;
1518
ecs_int_t marqueur_fin;
1526
ecs_int_t typ_geo_cel; /* Type géométrique élément en cours */
1528
ecs_int_t typ_geo_base[9] = {ECS_ELT_TYP_NUL,
1532
ECS_ELT_TYP_CEL_TETRA,
1533
ECS_ELT_TYP_CEL_PYRAM,
1534
ECS_ELT_TYP_CEL_PRISM,
1536
ECS_ELT_TYP_CEL_HEXA};
1538
ecs_size_t *def_cel_fac_pos = NULL;
1539
ecs_int_t *def_cel_fac_val = NULL;
1541
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1543
/*=================*/
1544
/* Initialisations */
1545
/*=================*/
1547
nbr_cel = table_def_cel->nbr;
1549
ecs_table__regle_en_pos(table_def_cel);
1551
/* Construction, pour les cellules, des tables principaux */
1552
/*--------------------------------------------------------*/
1554
/* Boucle de comptage pour l'allocation des sous-éléments */
1555
/*--------------------------------------------------------*/
1557
if (*table_def_fac != NULL) {
1558
nbr_fac_old = (*table_def_fac)->nbr;
1559
nbr_val_fac_old = ecs_table__ret_val_nbr(*table_def_fac);
1563
nbr_val_fac_old = 0;
1570
for (icel = 0; icel < nbr_cel; icel++ ) {
1572
ind_pos_cel = table_def_cel->pos[icel] - 1;
1574
nbr_val_cel = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1576
/* Traitement des cellules "classiques" */
1577
/*--------------------------------------*/
1579
if (nbr_val_cel < 9) {
1581
typ_geo_cel = typ_geo_base[nbr_val_cel];
1583
/* Boucle sur les sous-éléments définissant la cellulle */
1586
ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
1589
const ecs_sous_elt_t * sous_elt
1590
= &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]);
1595
idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som;
1598
nbr_val_fac += idef;
1604
/* Traitement des éléments de type "polyèdre" */
1605
/*--------------------------------------------*/
1609
ind_pos_sui = table_def_cel->pos[icel + 1] - 1;
1610
nbr_pos_loc = ind_pos_sui - ind_pos_cel;
1612
/* Convention : définition nodale cellule->sommets avec numéros de
1613
premiers sommets répétés en fin de liste pour marquer la fin
1618
for (isom = ind_pos_cel; isom < ind_pos_sui; isom++) {
1620
if (table_def_cel->val[isom] != marqueur_fin) {
1622
if (marqueur_fin == -1)
1623
marqueur_fin = table_def_cel->val[isom];
1634
} /* Fin de la boucle de comptage sur les cellules */
1636
/* Allocation et initialisation pour les faces */
1637
/* des tableaux associés aux définitions */
1638
/*----------------------------------------------*/
1640
nbr_fac += nbr_fac_old;
1641
nbr_val_fac += nbr_val_fac_old;
1643
if (*table_def_fac != NULL) {
1644
ecs_table__regle_en_pos(*table_def_fac);
1645
(*table_def_fac)->nbr = nbr_fac;
1646
ECS_REALLOC((*table_def_fac)->pos, nbr_fac + 1, ecs_size_t);
1647
ECS_REALLOC((*table_def_fac)->val, nbr_val_fac, ecs_int_t);
1650
*table_def_fac = ecs_table__alloue(nbr_fac,
1652
(*table_def_fac)->pos[0] = 1;
1655
ECS_MALLOC(def_cel_fac_pos, nbr_cel + 1, ecs_size_t);
1656
ECS_MALLOC(def_cel_fac_val, nbr_fac - nbr_fac_old, ecs_int_t);
1658
def_cel_fac_pos[0] = 1;
1660
/*=======================================*/
1661
/* Boucle sur les cellules à transformer */
1662
/*=======================================*/
1664
cpt_fac = nbr_fac_old;
1665
nbr_def = nbr_val_fac_old;
1667
for (icel = 0; icel < nbr_cel; icel++ ) {
1669
ind_pos_cel = table_def_cel->pos[icel] - 1;
1671
nbr_val_cel = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1673
/*--------------------------------------*/
1674
/* Traitement des éléments "classiques" */
1675
/*--------------------------------------*/
1677
if (nbr_val_cel < 9) {
1679
typ_geo_cel = typ_geo_base[nbr_val_cel];
1681
/* Boucle sur les faces définissant la cellulle */
1682
/*==============================================*/
1685
ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
1688
const ecs_sous_elt_t * sous_elt
1689
= &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]);
1691
/* Boucle sur les sommets définissant la face */
1692
/*--------------------------------------------*/
1695
idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som;
1698
/* Définition de la face en fonction des sommets */
1700
num_def = sous_elt->som[idef];
1702
(*table_def_fac)->val[nbr_def++]
1703
= table_def_cel->val[ind_pos_cel + num_def - 1];
1707
/* Position de la face dans sa définition en fonction des sommets */
1708
/*----------------------------------------------------------------*/
1710
(*table_def_fac)->pos[cpt_fac + 1]
1711
= (*table_def_fac)->pos[cpt_fac] + idef;
1713
/* Détermination de la cellule en fonction des faces */
1714
/*---------------------------------------------------*/
1716
def_cel_fac_val[cpt_fac - nbr_fac_old] = cpt_fac + 1;
1720
} /* Fin de la boucle sur les faces d'une cellule */
1725
/*--------------------------------------------*/
1726
/* Traitement des éléments de type "polyèdre" */
1727
/*--------------------------------------------*/
1731
/* Boucle sur les faces définissant le polyèdre */
1732
/*==============================================*/
1737
nbr_pos_loc = table_def_cel->pos[icel + 1] - 1 - ind_pos_cel;
1739
for (ind_pos_loc = 0; ind_pos_loc < nbr_pos_loc; ind_pos_loc++) {
1741
/* Définition de la face en fonction des sommets */
1743
if (table_def_cel->val[ind_pos_cel + ind_pos_loc] != marqueur_fin) {
1745
(*table_def_fac)->val[nbr_def++]
1746
= table_def_cel->val[ind_pos_cel + ind_pos_loc];
1748
if (marqueur_fin == -1)
1749
marqueur_fin = table_def_cel->val[ind_pos_cel + ind_pos_loc];
1753
/* Position de la face dans sa définition en fonction des sommets */
1757
(*table_def_fac)->pos[cpt_fac + 1] = nbr_def + 1;
1761
def_cel_fac_val[cpt_fac - nbr_fac_old] = cpt_fac + 1;
1763
/* Incrémentation du nombre de faces */
1770
} /* Fin de la boucle sur les faces du polyèdre */
1774
/* A ce point, on a : ifac
1775
= ecs_fic_elt_typ_liste_c[table_typ_geo_cel->val[icel]].nbr_sous_elt
1776
pour des éléments classiques, et ifac est égal au nombre de faces pour
1777
des éléments polyédriques. */
1779
def_cel_fac_pos[icel + 1] = def_cel_fac_pos[icel] + ifac;
1781
} /* Fin de la boucle sur les éléments */
1783
/* Mise à jour de la table */
1784
/*-------------------------*/
1786
ECS_FREE(table_def_cel->pos);
1787
ECS_FREE(table_def_cel->val);
1789
table_def_cel->pos = def_cel_fac_pos;
1790
table_def_cel->val = def_cel_fac_val;
1792
/* ecs_table__pos_en_regle(*table_def_fac); */
1793
ecs_table__pos_en_regle(table_def_cel);
1796
/*----------------------------------------------------------------------------
1797
* Fonction qui realise la fusion des definitions des elements
1798
*----------------------------------------------------------------------------*/
1801
ecs_table_def__fusionne(ecs_table_t *table_def,
1802
size_t *nbr_elt_cpct,
1803
ecs_tab_int_t *signe_elt)
1816
ecs_int_t num_inf_loc;
1817
ecs_int_t num_inf_min;
1818
ecs_int_t num_inf_min_cmp;
1824
size_t ind_pos_sup[3];
1825
size_t ind_pos_cmp[3];
1827
ecs_tab_int_t cpt_ref_inf;
1829
ecs_size_t *pos_recherche = NULL;
1830
ecs_int_t *val_recherche = NULL;
1832
ecs_tab_int_t tab_transf; /* Tableau de transformation */
1834
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
1836
ecs_table__regle_en_pos(table_def);
1838
/* Fusion des definitions des elements sur la liste initiale */
1839
/*------------------------------------------------------------*/
1841
/* Initialisations avec précautions pour cas vide */
1844
tab_transf.val = NULL;
1846
if (table_def == NULL)
1849
if (table_def == NULL)
1852
nbr_sup_ini = table_def->nbr;
1855
if (nbr_sup_ini < 1)
1858
if (signe_elt != NULL) {
1859
signe_elt->nbr = nbr_sup_ini;
1860
ECS_MALLOC(signe_elt->val, nbr_sup_ini, ecs_int_t);
1863
/* Comptage du nombre de sous-entités */
1868
ind_sup < table_def->pos[nbr_sup_ini] - 1;
1870
num_inf = table_def->val[ind_sup];
1871
if ((size_t)(ECS_ABS(num_inf)) > nbr_inf)
1872
nbr_inf = ECS_ABS(num_inf);
1875
/* Construction du tableau de recherche des sous-entités */
1876
/*-------------------------------------------------------*/
1878
ECS_MALLOC(pos_recherche, nbr_inf + 1, ecs_size_t);
1879
ECS_MALLOC(val_recherche, table_def->nbr, ecs_int_t);
1882
Comptage pour chaque élément de l'entité inférieure du nombre de fois où il
1883
apparaît comme élement d'indice minimal d'un élément de l'entité supérieure
1886
cpt_ref_inf = ecs_tab_int__cree_init(nbr_inf, 0);
1888
for (ind_sup = 0; ind_sup < table_def->nbr; ind_sup++) {
1890
ind_pos_loc = table_def->pos[ind_sup] - 1;
1892
num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
1893
while (++ind_pos_loc < table_def->pos[ind_sup + 1] -1) {
1894
num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
1895
if (num_inf_loc < num_inf_min)
1896
num_inf_min = num_inf_loc;
1899
assert(num_inf_min > 0 && (size_t)num_inf_min <= nbr_inf);
1901
cpt_ref_inf.val[num_inf_min - 1] += 1;
1905
/* Construction du vecteur */
1907
pos_recherche[0] = 1;
1909
for (ind_inf = 0; ind_inf < nbr_inf; ind_inf++) {
1911
pos_recherche[ind_inf + 1]
1912
= pos_recherche[ind_inf] + cpt_ref_inf.val[ind_inf];
1914
cpt_ref_inf.val[ind_inf] = 0;
1919
for (ind_sup = 0; ind_sup < table_def->nbr; ind_sup++) {
1921
ind_pos_loc = table_def->pos[ind_sup] - 1;
1923
num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
1924
while (ind_pos_loc < table_def->pos[ind_sup + 1] -1) {
1925
num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
1927
if (num_inf_loc < num_inf_min)
1928
num_inf_min = num_inf_loc;
1931
ind_inf = num_inf_min - 1;
1933
ind_pos = pos_recherche[ind_inf] - 1
1934
+ cpt_ref_inf.val[ind_inf];
1936
cpt_ref_inf.val[ind_inf] += 1;
1938
val_recherche[ind_pos] = ind_sup + 1;
1943
/* Libération du tableau auxiliaire */
1945
cpt_ref_inf.nbr = 0;
1946
ECS_FREE(cpt_ref_inf.val);
1948
/* Allocation du tableau de transformation */
1950
tab_transf = ecs_tab_int__cree(nbr_sup_ini);
1952
for (ind_sup = 0; ind_sup < nbr_sup_ini; ind_sup++)
1953
tab_transf.val[ind_sup] = ind_sup;
1956
/* Boucle principale de recherche sur les éléments supérieurs */
1957
/*------------------------------------------------------------*/
1960
Le premier élément ne peut pas être fusionné avec un élément précédent,
1961
la boucle commence donc au deuxième
1964
tab_transf.val[0] = 0;
1966
if (signe_elt != NULL)
1967
signe_elt->val[0] = 1;
1971
for (ind_sup = 1; ind_sup < table_def->nbr; ind_sup++) {
1973
/* Recherche élement entité inférieure de plus petit numéro référencé */
1975
ind_pos_sup[0] = table_def->pos[ind_sup ] - 1; /* début */
1976
ind_pos_sup[1] = table_def->pos[ind_sup + 1] - 1; /* fin */
1977
ind_pos_sup[2] = table_def->pos[ind_sup ] - 1; /* plus petit */
1979
ind_pos_loc = ind_pos_sup[0];
1981
num_inf_min = ECS_ABS(table_def->val[ind_pos_loc]);
1982
while (++ind_pos_loc < ind_pos_sup[1]) {
1983
num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
1984
if (num_inf_loc < num_inf_min) {
1985
num_inf_min = num_inf_loc;
1986
ind_pos_sup[2] = ind_pos_loc;
1991
On cherche des éléments de l'entité courante de plus petit numéro que
1992
l'entité courante ayant même plus petit élément de l'entité inférieure
1993
(recherche de candidats pour la fusion)
1996
ind_inf = num_inf_min - 1;
1999
for (pos_cmp = pos_recherche[ind_inf] - 1;
2000
pos_cmp < pos_recherche[ind_inf + 1] - 1;
2003
ind_cmp = val_recherche[pos_cmp] - 1;
2005
/* Repérage point de départ pour comparaison */
2007
if (ind_cmp < ind_sup) {
2009
ind_pos_cmp[0] = table_def->pos[ind_cmp ] - 1; /* début */
2010
ind_pos_cmp[1] = table_def->pos[ind_cmp + 1] - 1; /* fin */
2011
ind_pos_cmp[2] = table_def->pos[ind_cmp ] - 1; /* plus petit */
2013
assert(ind_pos_cmp[1] > ind_pos_cmp[0]);
2015
ind_pos_cmp[1] = table_def->pos[ind_cmp + 1] - 1; /* fin */
2017
ind_pos_loc = ind_pos_cmp[0];
2019
num_inf_min_cmp = ECS_ABS(table_def->val[ind_pos_loc]);
2020
while ((++ind_pos_loc) < ind_pos_cmp[1]) {
2021
num_inf_loc = ECS_ABS(table_def->val[ind_pos_loc]);
2022
if (num_inf_loc < num_inf_min_cmp) {
2023
num_inf_min_cmp = num_inf_loc;
2024
ind_pos_cmp[2] = ind_pos_loc;
2028
/* Comparaison des définitions */
2030
for (sgn = 1; sgn > -2; sgn -= 2) {
2032
ind_loc_sup = ind_pos_sup[2];
2033
ind_loc_cmp = ind_pos_cmp[2];
2038
if (ind_loc_sup == ind_pos_sup[1])
2039
ind_loc_sup = ind_pos_sup[0];
2042
if (ind_loc_cmp == ind_pos_cmp[1])
2043
ind_loc_cmp = ind_pos_cmp[0];
2044
else if ( ind_loc_cmp < ind_pos_cmp[0]
2045
|| ind_loc_cmp > ind_pos_cmp[1])
2046
ind_loc_cmp = ind_pos_cmp[1] - 1;
2048
} while ( ( ECS_ABS(table_def->val[ind_loc_sup])
2049
== ECS_ABS(table_def->val[ind_loc_cmp]))
2050
&& ind_loc_sup != ind_pos_sup[2]
2051
&& ind_loc_cmp != ind_pos_cmp[2]);
2053
if ( ind_loc_sup == ind_pos_sup[2]
2054
&& ind_loc_cmp == ind_pos_cmp[2])
2055
break; /* Sortie boucle sur signe parcours (1, -1, -3) */
2060
Si sgn = 1, les entités sont confondues, de même sens;
2061
Si sgn = -1, elles sont confondues, de sens inverse;
2062
Sinon, elles ne sont pas confondues
2065
if (sgn == 1 || sgn == -1)
2066
break; /* Sortie boucle sur pos_cmp pour recherche de candidats */
2068
} /* Fin de la comparaison référence/candidat (if ind_cmp < ind_sup) */
2070
} /* Fin boucle sur pos_cmp recherche de candidats */
2073
/* Si on a trouvé une entité à fusionner */
2075
if (sgn == 1 || sgn == -1) {
2077
assert(pos_cmp < pos_recherche[ind_inf + 1] - 1);
2079
if (sgn == 1 && (ind_pos_cmp[1] - ind_pos_cmp[0] == 2)) {
2082
Si l'on a deux entités inférieures par entité supérieure,
2083
la permutation cyclique peut trouver une égalité quel
2084
que soit le signe; on le corrige si nécessaire
2087
if ( ECS_ABS(table_def->val[ind_pos_sup[0]])
2088
!= ECS_ABS(table_def->val[ind_pos_cmp[0]]))
2093
tab_transf.val[ind_sup] = tab_transf.val[ind_cmp];
2095
if (signe_elt != NULL)
2096
signe_elt->val[ind_sup] = sgn;
2101
tab_transf.val[ind_sup] = cpt_sup_fin++;
2103
if (signe_elt != NULL)
2104
signe_elt->val[ind_sup] = 1;
2109
} /* Fin boucle sur les éléments de l'entité supérieure (que l'on fusionne) */
2111
/* Libération du tableau de recherche */
2113
ECS_FREE(pos_recherche);
2114
ECS_FREE(val_recherche);
2116
/* Compactage de la définition */
2117
/*-----------------------------*/
2120
Le premier élément ne peut pas être fusionné avec un élément précédent,
2121
la boucle commence donc au deuxième
2126
for (ind_sup = 1; ind_sup < table_def->nbr; ind_sup++) {
2128
if (tab_transf.val[ind_sup] > (ecs_int_t)(cpt_sup_fin - 1)) {
2130
/* Recopie définition sur partie compactée du tableau */
2132
for (pos_cpt = table_def->pos[cpt_sup_fin],
2133
pos_sup = table_def->pos[ind_sup];
2134
pos_sup < table_def->pos[ind_sup + 1];
2135
pos_cpt++, pos_sup++)
2136
table_def->val[pos_cpt - 1] = table_def->val[pos_sup - 1];
2140
table_def->pos[cpt_sup_fin] = pos_cpt;
2145
/* Redimensionnement de l'entité compactée */
2147
table_def->nbr = cpt_sup_fin;
2148
ECS_REALLOC(table_def->pos, cpt_sup_fin + 1, ecs_size_t);
2149
ECS_REALLOC(table_def->val,
2150
table_def->pos[cpt_sup_fin] - 1,
2153
ecs_table__pos_en_regle(table_def);
2155
/* Affectation de la nouvelle liste compactee */
2156
/*--------------------------------------------*/
2158
*nbr_elt_cpct = table_def->nbr;
2160
/* Renvoi du tableau de transformation */
2161
/*-------------------------------------*/
2166
/*----------------------------------------------------------------------------
2167
* Fonction qui construit la liste des cellules attachées à une liste
2168
* de faces fournie en argument.
2169
*----------------------------------------------------------------------------*/
2172
ecs_table_def__liste_cel_fac(const size_t nbr_fac,
2173
ecs_table_t *table_def_cel,
2174
const ecs_tab_int_t liste_fac)
2183
ecs_int_t *indic_cel;
2184
ecs_int_t *indic_fac;
2186
ecs_tab_int_t liste_cel;
2188
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2190
assert(table_def_cel != NULL);
2192
ecs_table__regle_en_pos(table_def_cel);
2194
/* Initialisations */
2198
nbr_cel = table_def_cel->nbr;
2202
ECS_MALLOC(indic_cel, nbr_cel, ecs_int_t);
2203
ECS_MALLOC(indic_fac, nbr_fac, ecs_int_t);
2205
for (icel = 0; icel < nbr_cel; icel++)
2206
indic_cel[icel] = 0;
2208
for (ifac = 0; ifac < nbr_fac; ifac++)
2209
indic_fac[ifac] = 0;
2211
for (ifac = 0; ifac < liste_fac.nbr; ifac++)
2212
indic_fac[liste_fac.val[ifac]] = 1;
2214
/* Première boucle sur les cellules : marquage */
2216
for (icel = 0; icel < nbr_cel; icel++) {
2219
= table_def_cel->pos[icel + 1]
2220
- table_def_cel->pos[icel];
2222
ipos = table_def_cel->pos[icel] - 1;
2224
for (iloc = 0; iloc < nbr_fac_cel; iloc++) {
2226
ifac = ECS_ABS(table_def_cel->val[ipos]) - 1;
2230
if (indic_fac[ifac] == 1 && indic_cel[icel] == 0) {
2231
indic_cel[icel] = 1;
2238
ECS_FREE(indic_fac);
2240
/* Seconde boucle sur les cellules : remplissage */
2242
liste_cel.nbr = cpt_cel;
2243
ECS_MALLOC(liste_cel.val, liste_cel.nbr, ecs_int_t);
2247
for (icel = 0; icel < nbr_cel; icel++) {
2249
if (indic_cel[icel] == 1)
2250
liste_cel.val[cpt_cel++] = icel;
2254
ECS_FREE(indic_cel);
2256
ecs_table__libere_pos(table_def_cel);
2261
/*----------------------------------------------------------------------------
2262
* Fonction qui remplace les références à des éléments
2263
* en des références à d'autres éléments liés aux premiers
2264
* par un tableau de renumérotation qui peut être signé.
2265
*----------------------------------------------------------------------------*/
2268
ecs_table_def__remplace_ref(ecs_table_t *table_def,
2269
ecs_tab_int_t *tab_old_new)
2271
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2273
if (table_def != NULL && tab_old_new != NULL) {
2275
ecs_int_t num_old, num_new;
2277
ecs_size_t ival, ival_deb, ival_fin;
2278
ecs_size_t cpt_val = 0;
2279
size_t nbr_elt = table_def->nbr;
2281
ecs_table__regle_en_pos(table_def);
2285
for (ielt = 0; ielt < nbr_elt; ielt++) {
2287
ival_fin = table_def->pos[ielt+1] - 1;
2289
for (ival = ival_deb; ival < ival_fin; ival++) {
2291
num_old = table_def->val[ival];
2293
num_new = tab_old_new->val[num_old -1];
2295
num_new = -tab_old_new->val[-num_old -1];
2297
table_def->val[cpt_val++] = num_new;
2300
ival_deb = table_def->pos[ielt+1] - 1;
2301
table_def->pos[ielt+1] = cpt_val + 1;
2304
ECS_REALLOC(table_def->val, cpt_val, ecs_int_t);
2306
ecs_table__pos_en_regle(table_def);
2310
/*----------------------------------------------------------------------------
2311
* Fonction qui construit un tableau de booleens conforme a une liste
2313
* Un sous-element est a `true'
2314
* s'il intervient dans la definition des elements
2315
*----------------------------------------------------------------------------*/
2318
ecs_table_def__cree_masque(ecs_tab_bool_t bool_sselt_select,
2319
ecs_table_t *table_def_elt)
2329
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2331
ecs_table__regle_en_pos(table_def_elt);
2333
nbr_elt = table_def_elt->nbr;
2335
for (ielt = 0; ielt < nbr_elt; ielt++) {
2337
pos_inf = table_def_elt->pos[ielt ] - 1;
2338
pos_sup = table_def_elt->pos[ielt + 1] - 1;
2340
for (ipos = pos_inf; ipos < pos_sup; ipos++) {
2342
num_sselt = ECS_ABS(table_def_elt->val[ipos]) - 1;
2344
/* If the sub-element has not yet been accounted for, mark it */
2346
if (bool_sselt_select.val[num_sselt] == false)
2347
bool_sselt_select.val[num_sselt] = true;
2352
ecs_table__libere_pos(table_def_elt);
2355
/*----------------------------------------------------------------------------
2356
* Suppression des sommets ne participant pas à la connectivité
2357
* et mise à jour de la connectivité.
2358
*----------------------------------------------------------------------------*/
2361
ecs_table_def__nettoie_nodal(size_t *n_vertices,
2362
ecs_coord_t **vtx_coords,
2363
ecs_table_t *table_def_fac,
2364
ecs_table_t *table_def_cel)
2375
ecs_tab_int_t tab_som_old_new;
2377
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2379
tab_som_old_new.nbr = 0;
2380
tab_som_old_new.val = NULL;
2382
if (vtx_coords == NULL)
2385
/* Initialisation du vecteur de renumérotation comme marqueur */
2387
nbr_som = *n_vertices;
2389
tab_som_old_new = ecs_tab_int__cree_init(nbr_som, 0);
2393
if (table_def_fac != NULL) {
2395
nbr_val = ecs_table__ret_val_nbr(table_def_fac);
2397
for (ival = 0; ival < nbr_val; ival++)
2398
tab_som_old_new.val[table_def_fac->val[ival] - 1] = 1;
2403
if (table_def_cel != NULL) {
2405
nbr_val = ecs_table__ret_val_nbr(table_def_cel);
2407
for (ival = 0; ival < nbr_val; ival++)
2408
tab_som_old_new.val[table_def_cel->val[ival] - 1] = 1;
2411
/* Préparation de la renumérotation */
2415
for (isom = 0; isom < nbr_som; isom++) {
2416
if (tab_som_old_new.val[isom] != 0)
2420
/* Si tous les sommets sont référencés, rien à faire */
2422
if (cpt_som == nbr_som) {
2423
tab_som_old_new.nbr = 0;
2424
ECS_FREE(tab_som_old_new.val);
2428
/* Compactage du tableau des sommets
2429
et transformation du marquer en renumérotation */
2433
for (isom = 0; isom < nbr_som; isom++) {
2435
if (tab_som_old_new.val[isom] != 0) {
2437
ipos_new = 3 * (cpt_som);
2438
ipos_old = 3 * isom;
2440
for (iloc = 0; iloc < 3; iloc++)
2441
(*vtx_coords)[ipos_new++] = (*vtx_coords)[ipos_old++];
2443
tab_som_old_new.val[isom] = cpt_som + 1;
2450
*n_vertices = cpt_som;
2451
ECS_REALLOC(*vtx_coords, cpt_som*3, ecs_coord_t);
2453
/* Mise à jour de la connectivité nodale */
2455
if (table_def_fac != NULL)
2456
ecs_table_def__remplace_ref(table_def_fac,
2459
if (table_def_cel != NULL)
2460
ecs_table_def__remplace_ref(table_def_cel,
2463
tab_som_old_new.nbr = 0;
2464
ECS_FREE(tab_som_old_new.val);
2467
/*----------------------------------------------------------------------------
2468
* Correction si nécessaire de l'orientation des éléments en connectivité
2469
* nodale. L'argument liste_cel_err est optionnel.
2470
*----------------------------------------------------------------------------*/
2473
ecs_table_def__orient_nodal(ecs_coord_t *vtx_coords,
2474
ecs_table_t *table_def_fac,
2475
ecs_table_t *table_def_cel,
2476
ecs_tab_int_t *liste_cel_err,
2487
ecs_int_t ret_orient;
2489
ecs_int_t cpt_orient_correc[ECS_ELT_TYP_FIN];
2490
ecs_int_t cpt_orient_erreur[ECS_ELT_TYP_FIN];
2492
ecs_int_t typ_cel_base[9] = {ECS_ELT_TYP_NUL,
2496
ECS_ELT_TYP_CEL_TETRA,
2497
ECS_ELT_TYP_CEL_PYRAM,
2498
ECS_ELT_TYP_CEL_PRISM,
2500
ECS_ELT_TYP_CEL_HEXA};
2502
ecs_int_t cpt_cel_erreur = 0;
2503
ecs_int_t cpt_cel_correc = 0;
2505
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2507
if (vtx_coords == NULL)
2510
if (table_def_fac != NULL)
2511
ecs_table__regle_en_pos(table_def_fac);
2513
if (table_def_cel != NULL)
2514
ecs_table__regle_en_pos(table_def_cel);
2516
assert(vtx_coords != NULL);
2518
for (typ_elt = 0; typ_elt < ECS_ELT_TYP_FIN; typ_elt++) {
2519
cpt_orient_correc[typ_elt] = 0;
2520
cpt_orient_erreur[typ_elt] = 0;
2523
printf(_("\n Element orientation check.\n\n"));
2527
if (table_def_fac != NULL) {
2529
nbr_fac = table_def_fac->nbr;
2531
for (ifac = 0; ifac < nbr_fac; ifac++) {
2533
ipos_fac = table_def_fac->pos[ifac ] - 1;
2534
nbr_som_loc = table_def_fac->pos[ifac + 1] - 1 - ipos_fac;
2536
if (nbr_som_loc == 4) {
2539
= _orient_quad(vtx_coords,
2540
&(table_def_fac->val[ipos_fac]),
2544
cpt_orient_erreur[ECS_ELT_TYP_FAC_QUAD] += 1;
2545
else if (ret_orient > 0)
2546
cpt_orient_correc[ECS_ELT_TYP_FAC_QUAD] += 1;
2554
if (table_def_cel != NULL) {
2556
ecs_tab_int_t face_index, face_marker, edges;
2559
face_index.val = NULL;
2561
face_marker.nbr = 0;
2562
face_marker.val = NULL;
2567
nbr_cel = table_def_cel->nbr;
2569
for (icel = 0; icel < nbr_cel; icel++) {
2571
ipos_cel = table_def_cel->pos[icel] - 1;
2573
nbr_som_loc = table_def_cel->pos[icel + 1] - 1 - ipos_cel;
2575
if (nbr_som_loc < 9)
2576
typ_elt = typ_cel_base[nbr_som_loc];
2578
typ_elt = ECS_ELT_TYP_CEL_POLY;
2582
case ECS_ELT_TYP_CEL_TETRA:
2584
ret_orient = _orient_tetra(vtx_coords,
2585
&(table_def_cel->val[ipos_cel]));
2588
case ECS_ELT_TYP_CEL_PYRAM:
2591
= _orient_pyram(vtx_coords,
2592
&(table_def_cel->val[ipos_cel]),
2596
case ECS_ELT_TYP_CEL_PRISM:
2599
= _orient_prism(vtx_coords,
2600
&(table_def_cel->val[ipos_cel]),
2605
case ECS_ELT_TYP_CEL_HEXA:
2608
= _orient_hexa(vtx_coords,
2609
&(table_def_cel->val[ipos_cel]),
2614
default: /* ECS_ELT_TYP_CEL_POLY */
2616
ret_orient = _orient_polyhedron(vtx_coords,
2617
&(table_def_cel->val[ipos_cel]),
2618
( table_def_cel->pos[icel + 1]
2619
- table_def_cel->pos[icel]),
2627
if (ret_orient < 0) {
2628
cpt_orient_erreur[typ_elt] += 1;
2629
if (liste_cel_err != NULL) {
2630
if (liste_cel_err->nbr == 0) {
2631
liste_cel_err->nbr = nbr_cel;
2632
ECS_MALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t);
2634
liste_cel_err->val[cpt_cel_erreur] = icel;
2636
cpt_cel_erreur += 1;
2639
else if (ret_orient > 0) {
2640
cpt_orient_correc[typ_elt] += 1;
2641
cpt_cel_correc += 1;
2645
if (face_index.nbr > 0)
2646
ECS_FREE(face_index.val);
2648
if (face_marker.nbr > 0)
2649
ECS_FREE(face_marker.val);
2652
ECS_FREE(edges.val);
2655
/* Impression de messages d'avertissement */
2657
for (typ_elt = 0; typ_elt < ECS_ELT_TYP_FIN; typ_elt++) {
2658
if (cpt_orient_correc[typ_elt] > 0) {
2660
printf(_("%d elements of type %s had to be re-oriented\n"),
2661
(int)(cpt_orient_correc[typ_elt]),
2662
_(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2664
if (cpt_orient_erreur[typ_elt] > 0) {
2665
if (correc_orient == true) {
2667
printf(_("%d elements of type %s were impossible to re-orient\n"),
2668
(int)(cpt_orient_erreur[typ_elt]),
2669
_(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2673
printf(_("%d elements of type %s are mis-oriented or highly warped\n"),
2674
(int)(cpt_orient_erreur[typ_elt]),
2675
_(ecs_fic_elt_typ_liste_c[typ_elt].nom));
2680
/* Redimensionnement de la liste des cellules toujours mal orientées */
2682
if (liste_cel_err != NULL) {
2683
if (liste_cel_err->nbr > 0) {
2684
liste_cel_err->nbr = cpt_cel_erreur;
2685
ECS_REALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t);
2689
if (table_def_fac != NULL)
2690
ecs_table__libere_pos(table_def_fac);
2692
if (table_def_cel != NULL)
2693
ecs_table__libere_pos(table_def_cel);
2696
/*----------------------------------------------------------------------------
2697
* Fusion des sommets confondus d'après la longueur des arêtes des faces.
2698
* La connectivité des faces est mise à jour.
2699
*----------------------------------------------------------------------------*/
2702
ecs_table_def__nettoie_som_fac(size_t *n_vertices,
2703
ecs_coord_t **vtx_coords,
2704
ecs_table_t *table_def_fac)
2727
ecs_tab_int_t tab_equiv_som;
2729
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2731
if (vtx_coords == NULL || table_def_fac == NULL)
2734
ecs_table__regle_en_pos(table_def_fac);
2736
nbr_fac = table_def_fac->nbr;
2737
nbr_som = *n_vertices;
2739
lng_are_min = 1.0e-15;
2740
if (getenv("CS_PREPROCESS_MIN_EDGE_LEN") != NULL) {
2741
lng_are_min = atof(getenv("CS_PREPROCESS_MIN_EDGE_LEN"));
2743
lng_min_2 = lng_are_min * ECS_ABS(lng_are_min);
2745
/* Marquage des sommets */
2746
/*----------------------*/
2748
tab_equiv_som = ecs_tab_int__cree_init(nbr_som, -1);
2752
/* Boucle sur les faces pour détecter les sommets confondus */
2753
/*----------------------------------------------------------*/
2755
for (ifac = 0; ifac < nbr_fac; ifac++) {
2757
ipos_deb = table_def_fac->pos[ifac] - 1;
2758
nbr_som_fac = (table_def_fac->pos[ifac + 1] - 1) - ipos_deb;
2760
for (isom = 0; isom < nbr_som_fac; isom++) {
2762
isom_0 = table_def_fac->val[ipos_deb + isom] - 1;
2763
isom_1 = table_def_fac->val[ ipos_deb + ((isom + 1)%nbr_som_fac)] -1;
2764
ipos_0 = isom_0 * 3;
2765
ipos_1 = isom_1 * 3;
2769
for (icoo = 0; icoo < 3; icoo++) {
2770
delta_coo = (*vtx_coords)[ipos_1 + icoo]
2771
- (*vtx_coords)[ipos_0 + icoo];
2772
lng_are_2 += delta_coo * delta_coo;
2775
/* Sommets confondus détectés */
2777
if (lng_are_2 < lng_min_2) {
2778
_table_def__maj_equiv_som(isom_0, isom_1, &tab_equiv_som);
2785
/* Fusion si nécessaire */
2787
if (cpt_fusion > 0) {
2789
ecs_tab_int_t tab_som_old_new;
2791
printf(_("\nMesh verification:\n\n"
2792
" %d vertices belong to edges of length less than %g\n"
2793
" and will be merged; (this tolerance may be modified\n"
2794
" using the CS_PREPROCESS_MIN_EDGE_LEN environment variable).\n"),
2795
(int)cpt_fusion, lng_are_min);
2797
tab_som_old_new = _table_def__fusion_som(n_vertices,
2801
_table_def__maj_fac_som(table_def_fac, &tab_som_old_new);
2803
ECS_FREE(tab_som_old_new.val);
2808
/* Libération mémoire et retour */
2810
ECS_FREE(tab_equiv_som.val);
2812
ecs_table__pos_en_regle(table_def_fac);
2815
/*----------------------------------------------------------------------------
2816
* Fonction qui supprime les éventuelles faces dégénérées
2817
*----------------------------------------------------------------------------*/
2820
ecs_table_def__nettoie_fac(ecs_table_t *table_def_fac)
2822
ecs_int_t nbr_fac_old;
2823
ecs_int_t nbr_fac_new;
2830
ecs_int_t ind_val_deb;
2831
ecs_int_t ind_val_fin;
2833
ecs_tab_int_t tab_fac_old_new;
2835
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2837
tab_fac_old_new.nbr = 0;
2838
tab_fac_old_new.val = NULL;
2840
if (table_def_fac == NULL)
2841
return tab_fac_old_new;
2843
ecs_table__regle_en_pos(table_def_fac);
2845
/* Initialisations */
2847
nbr_fac_old = table_def_fac->nbr;
2849
/* Boucle sur les faces (renumérotation) */
2850
/* ------------------------------------- */
2852
tab_fac_old_new.nbr = nbr_fac_old;
2853
ECS_MALLOC(tab_fac_old_new.val, tab_fac_old_new.nbr, ecs_int_t);
2858
for (ind_fac = 0; ind_fac < nbr_fac_old; ind_fac++) {
2860
ind_val_deb = table_def_fac->pos[ind_fac ] - 1;
2861
ind_val_fin = table_def_fac->pos[ind_fac + 1] - 1;
2863
if (ind_val_fin - ind_val_deb > 2) {
2865
for (ind_val = ind_val_deb; ind_val < ind_val_fin; ind_val++)
2866
table_def_fac->val[cpt_val++] = table_def_fac->val[ind_val];
2868
tab_fac_old_new.val[ind_fac] = 1 + cpt_fac++;
2870
table_def_fac->pos[cpt_fac] = cpt_val + 1;
2874
tab_fac_old_new.val[ind_fac] = 0;
2878
nbr_fac_new = cpt_fac;
2880
/* Si on n'a pas de faces dégénérées, on peut sortir */
2882
if (nbr_fac_new == nbr_fac_old) {
2884
tab_fac_old_new.nbr = 0;
2885
ECS_FREE(tab_fac_old_new.val);
2887
ecs_table__pos_en_regle(table_def_fac);
2889
return tab_fac_old_new;
2892
printf(_("\nMesh verification:\n\n"
2893
" Removal of %d degenerate faces:\n"
2894
" Initial number of faces : %10d\n"
2895
" Number of faces after processing : %10d\n"),
2896
(int)(nbr_fac_old - nbr_fac_new), (int)nbr_fac_old, (int)nbr_fac_new);
2898
/* On redimensionne le tableau des faces */
2900
table_def_fac->nbr = nbr_fac_new;
2901
ECS_REALLOC(table_def_fac->pos, nbr_fac_new + 1, ecs_size_t);
2902
ECS_REALLOC(table_def_fac->val, cpt_val, ecs_int_t);
2904
ecs_table__pos_en_regle(table_def_fac);
2906
return tab_fac_old_new;
2909
/*----------------------------------------------------------------------------
2910
* Fonction qui renvoie un tableau associant un type à chaque face, sous
2911
* forme de masque : 0 pour face isolée, 1 ou 2 pour face de bord (1 si
2912
* cellule avec cette face normale sortante, 2 si cellule avec cette face
2913
* normale entrante), 1+2 = 3 pour face interne, et 4 ou plus pour tous
2914
* les autres cas, correspondant à une erreur de connectivité (+4 pour faces
2915
* voyant au moins deux cellules avec face normale sortante, +8 pour faces
2916
* voyant au moins deux cellules avec face normale entrante).
2918
* Le type de chaque face pourra être modifié ultérieurement en fonction
2919
* des informations de périodicité.
2920
*----------------------------------------------------------------------------*/
2923
ecs_table_def__typ_fac_cel(ecs_table_t *table_def_cel,
2924
ecs_table_t *table_def_fac)
2933
ecs_tab_int_t typ_fac;
2935
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
2940
if (table_def_cel == NULL)
2943
/* Initialisations */
2945
ecs_table__regle_en_pos(table_def_cel);
2947
nbr_cel = table_def_cel->nbr;
2948
nbr_fac = table_def_fac->nbr;
2950
/* Type de face selon le nombre de cellules voisines : 0 pour face isolée,
2951
1 ou 2 pour face de bord, 3 pour faces internes, et 4 pour tous les
2952
autres cas (faces voyant au moins deux cellules sur un même côté) */
2954
typ_fac.nbr = nbr_fac;
2956
ECS_MALLOC(typ_fac.val, nbr_fac, ecs_int_t);
2958
for (ifac = 0; ifac < nbr_fac; ifac++)
2959
typ_fac.val[ifac] = 0;
2961
/* Boucle sur les cellules : marquage */
2962
/*------------------------------------*/
2964
for (icel = 0; icel < nbr_cel; icel++ ) {
2966
for (ipos = table_def_cel->pos[icel] - 1;
2967
ipos < table_def_cel->pos[icel+1] - 1;
2970
num_fac = table_def_cel->val[ipos];
2972
ifac = ECS_ABS(num_fac) - 1;
2974
if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 0)
2975
typ_fac.val[ifac] += 1;
2977
else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 0)
2978
typ_fac.val[ifac] += 2;
2981
if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 1)
2982
typ_fac.val[ifac] = typ_fac.val[ifac] | 4;
2983
else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 1)
2984
typ_fac.val[ifac] = typ_fac.val[ifac] | 8;
2991
ecs_table__libere_pos(table_def_cel);
2996
/*----------------------------------------------------------------------------
2997
* Fonction qui renvoie un tableau associant un type à chaque face les
2998
* numéros des cellules définies par cette face (normale sortante,
2999
* puis normale entrante). On affecte une valeur 0 lorsqu'il n'y a pas de
3000
* cellule correspondante directe (la périodicité n'est donc pas prise en
3001
* compte à ce niveau).
3003
* On suppose que la cohérence du maillage a déjà été véridifiée et
3004
* qu'aucune face n'appartient à plus d'une cellule par côté.
3005
*----------------------------------------------------------------------------*/
3008
ecs_table_def__fac_cel(ecs_table_t *table_def_cel,
3009
ecs_table_t *table_def_fac)
3018
ecs_tab_int_t fac_cel;
3020
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
3025
if (table_def_cel == NULL)
3028
/* Initialisations */
3030
ecs_table__regle_en_pos(table_def_cel);
3032
nbr_cel = table_def_cel->nbr;
3033
nbr_fac = table_def_fac->nbr;
3035
/* Allocation et mise à zéro des connectivités */
3037
fac_cel.nbr = nbr_fac * 2;
3039
ECS_MALLOC(fac_cel.val, fac_cel.nbr, ecs_int_t);
3041
for (ipos = 0; ipos < fac_cel.nbr; ipos++)
3042
fac_cel.val[ipos] = 0;
3044
/* Boucle sur les cellules : marquage */
3045
/*------------------------------------*/
3047
for (icel = 0; icel < nbr_cel; icel++ ) {
3049
for (ipos = table_def_cel->pos[icel] - 1;
3050
ipos < table_def_cel->pos[icel+1] - 1;
3053
num_fac = table_def_cel->val[ipos];
3055
ifac = ECS_ABS(num_fac) - 1;
3058
assert(fac_cel.val[ifac*2] == 0);
3059
fac_cel.val[ifac*2] = icel + 1;
3062
assert(fac_cel.val[ifac*2 + 1] == 0);
3063
fac_cel.val[ifac*2 + 1] = icel + 1;
3070
ecs_table__libere_pos(table_def_cel);
3075
/*----------------------------------------------------------------------------*/
added
removed
preprocessor/base/ecs_table_def.c
