ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/diamesh/src/m3d_move.cpp
Revision: 5
Committed: Tue Jun 12 20:26:34 2007 UTC (17 years, 11 months ago)
Original Path: magic/lib/diamesh/diamesh/src/m3d_move.cpp
File size: 7494 byte(s)
Log Message:

File Contents

# User Rev Content
1 5 /* BOUGE DE POINTS SAUCE P.L.GEORGE Article INRIA STRUCOME 1994 */
2     #include <stdio.h>
3     #include <math.h>
4     #include "m3d_struct.h"
5     #include "m3d_const.h"
6     #include "m3d_hotes.h"
7     #include "m3d_erreur.h"
8     #include "prototype.h"
9     extern GEST_MEM *gest ;
10     extern int debug ;
11     extern float crit_min ;
12     int m3d_move(int nump)
13     {
14     NOEUD *n1 ;
15     CONNEC *connec ;
16     TETRAEDRE *tab_tetra[1000] ;
17     int nb_ele, i, j, k, *numele ;
18     float tab_crit[1000], vec_opt[3], cmin ;
19     float tab_opt[3000] ;
20     int tabele[3000], sens ;
21     float gamma, vn[3], xi, yi, zi, vab[3], vac[3], vbc[3], pvec[3], hauteur, vap[3] ;
22     float dab, dac, dbc, perimetre, norme, coeff, alpha ;
23     float delta[3] ;
24     float vec_save[3], bmin, bmax, *coord, eps, crit_save ;
25     float calpha, calpha_moins_eps, vresu[3] ;
26    
27     /* coordonnees et connectivites */
28     coord = gest->coord ;
29     numele = gest->numele ;
30    
31     if (nump<gest->nb_init) return(VRAI) ; /* on ne bouge pas un noeud sur la peau */
32     ADRESSE(noeud,nump,n1)
33     if (n1 == NULL)
34     {
35     m3d_erreur(ERR_SYST) ;
36     return(FAUX) ;
37     }
38     if (n1->mark == KILLED) return(VRAI) ;
39    
40     /* sauvegarde des coordonnees du noeud */
41     vec_save[0] = coord[x(nump)] ;
42     vec_save[1] = coord[y(nump)] ;
43     vec_save[2] = coord[z(nump)] ;
44    
45     /* tetraedres connectes au noeud */
46     nb_ele = 0 ;
47     connec = n1->lis_con ;
48     while (connec != NULL)
49     {
50     if ((connec->tetra)->mark != KILLED)
51     {
52     tab_tetra[nb_ele] = connec->tetra ;
53     /* determination des noeuds de la face de la boule */
54     k = 0 ;
55     for (j=0;j<4;j++)
56     if (gest->numele[4*(connec->tetra->num)+ j] != n1->num)
57     {
58     tabele[3*nb_ele + k] = gest->numele[4*(connec->tetra->num)+ j] ;
59     k++ ;
60     }
61     if (k!=3) return(FAUX) ;
62     /* determination du point optimal de la face nb_ele */
63     xi = (coord[x(tabele[3*nb_ele])] + coord[x(tabele[3*nb_ele+1])] + coord[x(tabele[3*nb_ele+2])])/3. ;
64     yi = (coord[y(tabele[3*nb_ele])] + coord[y(tabele[3*nb_ele+1])] + coord[y(tabele[3*nb_ele+2])])/3. ;
65     zi = (coord[z(tabele[3*nb_ele])] + coord[z(tabele[3*nb_ele+1])] + coord[z(tabele[3*nb_ele+2])])/3. ;
66    
67     /* calcul de la normale a la face */
68     vab[0] = coord[x(tabele[3*nb_ele+1])] - coord[x(tabele[3*nb_ele])] ;
69     vab[1] = coord[y(tabele[3*nb_ele+1])] - coord[y(tabele[3*nb_ele])] ;
70     vab[2] = coord[z(tabele[3*nb_ele+1])] - coord[z(tabele[3*nb_ele])] ;
71     dab = NORME(vab) ;
72     vac[0] = coord[x(tabele[3*nb_ele+2])] - coord[x(tabele[3*nb_ele])] ;
73     vac[1] = coord[y(tabele[3*nb_ele+2])] - coord[y(tabele[3*nb_ele])] ;
74     vac[2] = coord[z(tabele[3*nb_ele+2])] - coord[z(tabele[3*nb_ele])] ;
75     dac = NORME(vac) ;
76     vbc[0] = coord[x(tabele[3*nb_ele+2])] - coord[x(tabele[3*nb_ele+1])] ;
77     vbc[1] = coord[y(tabele[3*nb_ele+2])] - coord[y(tabele[3*nb_ele+1])] ;
78     vbc[2] = coord[z(tabele[3*nb_ele+2])] - coord[z(tabele[3*nb_ele+1])] ;
79     dbc = NORME(vbc) ;
80    
81     /* calcul du critere 2 D */
82     pvec[0] = PVECX(vab,vac) ;
83     pvec[1] = PVECY(vab,vac) ;
84     pvec[2] = PVECZ(vab,vac) ;
85     perimetre = dab + dbc + dac ;
86     hauteur = (perimetre/3.) * 0.8 ;
87    
88     vn[0] = PVECX(vab,vac) ;
89     vn[1] = PVECY(vab,vac) ;
90     vn[2] = PVECZ(vab,vac) ;
91     norme = NORME(vn) ;
92     if (norme<EPSILON)
93     {
94     // if (debug) aff_text("normale nulle M3D_MOVE\n ") ;
95     m3d_erreur(ERR_SYST) ;
96     return(FAUX) ;
97     }
98     for (i=0;i<3;i++) vn[i] = vn[i]/norme ;
99     /* calcul du vecteur AP et orientation de la face */
100     vap[0] = coord[x(nump)] - coord[x(tabele[3*nb_ele])] ;
101     vap[1] = coord[y(nump)] - coord[y(tabele[3*nb_ele])] ;
102     vap[2] = coord[z(nump)] - coord[z(tabele[3*nb_ele])] ;
103     if (PROSCA(vap,vn) < EPSILON)
104     {
105     for (i=0;i<3;i++) vn[i] = -1. * vn[i] ;
106     i = tabele[3*nb_ele+1] ;
107     tabele[3*nb_ele+1] = tabele[3*nb_ele+2] ;
108     tabele[3*nb_ele+2] = i ;
109     }
110     /* point optimal */
111     tab_opt[x(nb_ele)] = xi + hauteur * vn[0] ;
112     tab_opt[y(nb_ele)] = yi + hauteur * vn[1] ;
113     tab_opt[z(nb_ele)] = zi + hauteur * vn[2] ;
114    
115     nb_ele ++ ;
116     }
117     connec = connec->suivant ;
118     }
119     if (nb_ele == 0) return(VRAI) ;
120     /* calcul de la qualite de chaque tetraedre */
121     cmin = 1. ;
122     for (i=0;i<nb_ele;i++)
123     {
124     /* calcul du critere */
125     tab_crit[i] = m3d_e_qual2(gest->coord,tabele[3*i],tabele[3*i+1],tabele[3*i+2],nump) ;
126     if (tab_crit[i]<EPSILON)
127     {
128     // if (debug) aff_text("critere nul M3D_MOVE\n ") ;
129     m3d_erreur(ERR_SYST) ;
130     return(FAUX) ;
131     }
132     if (tab_crit[i] < cmin) cmin = tab_crit[i] ;
133     }
134     /* calcul des coordonnees du point ideal */
135     gamma = 0. ;
136     for (i=0;i<nb_ele;i++)
137     {
138     gamma = gamma + 1./(tab_crit[i]*tab_crit[i]) ;/* somme des carres des inverses */
139     }
140     /* calcul du point optimal pour la boule */
141     gamma = 1./gamma ;
142     vec_opt[0] = 0. ;
143     vec_opt[1] = 0. ;
144     vec_opt[2] = 0. ;
145    
146     for (i=0;i<nb_ele;i++)
147     {
148     coeff = gamma/(tab_crit[i]*tab_crit[i]) ;
149     vec_opt[0] = vec_opt[0] + coeff * tab_opt[x(i)] ;
150     vec_opt[1] = vec_opt[1] + coeff * tab_opt[y(i)] ;
151     vec_opt[2] = vec_opt[2] + coeff * tab_opt[z(i)] ;
152     }
153    
154     delta[0] = vec_opt[0] - coord[x(nump)] ;
155     delta[1] = vec_opt[1] - coord[y(nump)] ;
156     delta[2] = vec_opt[2] - coord[z(nump)] ;
157     /* sauvegarde du critere mini */
158     crit_save = cmin ;
159    
160     /* initialisation des bornes du parametre */
161     bmin = 0. ;
162     bmax = 1. ;
163     sens = 1 ;/* augmnentation */
164    
165     vresu[0] = vec_save[0] ;
166     vresu[1] = vec_save[1] ;
167     vresu[2] = vec_save[2] ;
168    
169     for (i=0;i<5;i++)
170     {
171     alpha = 0.5 * (bmin + bmax) ;
172     coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
173     coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
174     coord[z(nump)] = vec_save[2] + alpha * delta[2] ;
175    
176     /* evaluation de la qualite pour la valeur alpha */
177     /* on evalue la qualite cmin pour alpha */
178     calpha = 1. ;
179     for (j=0;j<nb_ele;j++)
180     {
181     /* calcul du critere */
182     tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3*j],tabele[3*j+1],tabele[3*j+2],nump) ;
183     if (tab_crit[j] < calpha) calpha = tab_crit[j] ;
184     }
185     /* la valeur du critere pour alpha est superieure, il faut la stocker */
186     if (calpha > crit_save)
187     {
188     vresu[0] = coord[x(nump)] ;
189     vresu[1] = coord[y(nump)] ;
190     vresu[2] = coord[z(nump)] ;
191     crit_save = calpha ;
192     }
193    
194     /* on evalue la qualite c1 pour alpha - eps */
195     eps = (bmax-bmin)/50. ;
196     coord[x(nump)] = vec_save[0] + (alpha - eps) * delta[0] ;
197     coord[y(nump)] = vec_save[1] + (alpha - eps) * delta[1] ;
198     coord[z(nump)] = vec_save[2] + (alpha - eps) * delta[2] ;
199     calpha_moins_eps = 1. ;
200     for (j=0;j<nb_ele;j++)
201     {
202     /* calcul du critere */
203     tab_crit[j] = m3d_e_qual2(gest->coord,tabele[3*j],tabele[3*j+1],tabele[3*j+2],nump) ;
204     if (tab_crit[j] < calpha_moins_eps) calpha_moins_eps = tab_crit[j] ;
205     }
206    
207     coord[x(nump)] = vec_save[0] + alpha * delta[0] ;
208     coord[y(nump)] = vec_save[1] + alpha * delta[1] ;
209     coord[z(nump)] = vec_save[2] + alpha * delta[2] ;
210    
211     if (calpha_moins_eps < calpha)
212     {
213     bmin = alpha ;
214     /* bmax est inchange */
215     }
216     else
217     {
218     /* bmin est inchange */
219     bmax = alpha ;
220     }
221     }
222     coord[x(nump)] = vresu[0] ;
223     coord[y(nump)] = vresu[1] ;
224     coord[z(nump)] = vresu[2] ;
225    
226     return(VRAI) ;
227     }
228    
229    
230    
231    
232    
233    
234