ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/vct_surface.cpp
Revision: 72
Committed: Mon Mar 31 23:21:51 2008 UTC (17 years, 1 month ago) by francois
Original Path: magic/lib/geometrie/geometrie/src/vct_surface.cpp
File size: 11903 byte(s)
Log Message: 
optimisation vectorisation

File Contents

# User Rev Content
1 souaissa 66 #include "gestionversion.h"
2    
3     //---------------------------------------------------------------------------
4    
5    
6     #pragma hdrstop
7    
8     #include "vct_surface.h"
9     #include "mg_surface.h"
10     #include <iomanip.h>
11     //---------------------------------------------------------------------------
12    
13     #pragma package(smart_init)
14    
15    
16     VCT_SURFACE::VCT_SURFACE(MG_SURFACE* elemgeo):VCT_ELEMENT_GEOMETRIQUE(elemgeo)
17     {
18 souaissa 69 int indx_premier_ptctr;
19     TPL_LISTE_ENTITE<double> nurbs_params;
20 souaissa 66 elem_geo->get_param_NURBS(indx_premier_ptctr,nurbs_params);
21     nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);
22 souaissa 69 int nb_u=nurbs_params.get(3);
23     int nb_v=nurbs_params.get(4);
24 souaissa 66 double2 ZERO=0.;
25 souaissa 69 OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
26 souaissa 66
27 souaissa 69 OT_VECTEUR_4DD V1,V2,V;
28 souaissa 66
29    
30     for(int r=0;r<nb_v;r++)
31     {
32     for(int s=0;s<nb_u-1;s++)
33     {
34 souaissa 69
35     V1[0]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s);
36     V1[1]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 1);
37     V1[2]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 2);
38     V1[3]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 3);
39 souaissa 71
40 souaissa 69 lst_points.insert(lst_points.end(),V1);
41 souaissa 71
42 souaissa 69 V2[0]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1));
43     V2[1]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1) + 1);
44     V2[2]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1) + 2);
45     V2[3]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1) + 3);
46 souaissa 71
47 souaissa 69 lst_points.insert(lst_points.end(),V2);
48 souaissa 66
49 souaissa 69 V=V2-V1;
50 souaissa 66
51 souaissa 69 double2 norm_au_carre=(V[0]*V[0])+(V[1]*V[1])+(V[2]*V[2])+(V[3]*V[3]);
52     double2 norm=norm_au_carre^0.5;
53 souaissa 66
54 souaissa 69 if(V==VCT_NUL)
55     lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
56    
57 souaissa 66 if(norm!=ZERO){
58 souaissa 69 V=1./norm*V;
59     lst_vecteurs.insert(lst_vecteurs.end(),V);
60 souaissa 66 }
61    
62     }
63     }
64    
65     for(int r=0;r<nb_v-1;r++)
66     {
67     for(int s=0;s<nb_u;s++)
68     {
69 souaissa 69 V1[0]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s);
70     V1[1]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 1);
71     V1[2]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 2);
72     V1[3]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 3);
73 souaissa 71
74 souaissa 69 lst_points.insert(lst_points.end(),V1);
75 souaissa 71
76 souaissa 69 V2[0]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 * nb_u + 4 * s);
77     V2[1]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 * nb_u + 4 * s + 1);
78     V2[2]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 * nb_u + 4 * s + 2);
79     V2[3]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 * nb_u + 4 * s + 3);
80 souaissa 71
81 souaissa 69 lst_points.insert(lst_points.end(),V2);
82     V=V2-V1;
83 souaissa 66
84 souaissa 69 double2 norm_au_carre=(V[0]*V[0])+(V[1]*V[1])+(V[2]*V[2])+(V[3]*V[3]);
85     double2 norm=norm_au_carre^0.5;
86 souaissa 66
87 souaissa 69 if(V==VCT_NUL)
88     lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
89 souaissa 66
90 souaissa 69 if (norm!=ZERO){
91     V=1./norm*V;
92     lst_vecteurs.insert(lst_vecteurs.end(),V);
93 souaissa 66
94     }
95     }
96     }
97    
98 souaissa 69
99    
100 souaissa 66 }
101    
102    
103    
104 souaissa 69
105    
106    
107     VCT_SURFACE::VCT_SURFACE(VCT_SURFACE& mdd):VCT_ELEMENT_GEOMETRIQUE(mdd.elem_geo)
108 souaissa 66 {
109 souaissa 69 lst_vecteurs=mdd.lst_vecteurs;
110     lst_points=mdd.lst_points;
111     nb_points=mdd.nb_points ;
112     }
113 souaissa 66
114    
115 souaissa 69 VCT_SURFACE::~VCT_SURFACE()
116     {
117     }
118    
119     std::vector<OT_VECTEUR_4DD> &VCT_SURFACE::get_vecteurs()
120     {
121     return lst_vecteurs;
122     }
123    
124 francois 72 std::vector<OT_VECTEUR_4DD>& VCT_SURFACE::get_points_controle(void)
125 souaissa 71 {
126     return lst_points;
127     }
128 souaissa 69
129 souaissa 71
130 souaissa 69 OT_TENSEUR VCT_SURFACE:: calcule_tenseur_metrique()
131     {
132     OT_TENSEUR tns(lst_vecteurs);
133     return tns;
134     }
135    
136     OT_VECTEUR_4DD VCT_SURFACE::calcule_barycentre()
137     {
138     OT_VECTEUR_4DD barycentre(0,0,0,0);
139    
140 souaissa 66 for(int i=0;i<nb_points;i++)
141     {
142 souaissa 71 barycentre+=lst_points[i];
143 souaissa 66 }
144    
145 souaissa 69 barycentre*=1./nb_points;
146     return barycentre;
147 souaissa 66 }
148    
149    
150    
151 souaissa 69 int VCT_SURFACE::get_nb_points()
152     {
153     return nb_points ;
154     }
155    
156    
157     OT_TENSEUR VCT_SURFACE::calcule_covariance(void)
158 souaissa 66 {
159 souaissa 71 OT_TENSEUR COVARIANCE(4);
160 souaissa 69 std::vector<OT_VECTEUR_4DD> Points_Centre_au_barycentre(nb_points);
161     OT_VECTEUR_4DD barycentre=calcule_barycentre();
162 souaissa 66
163 souaissa 69 for( int i=0;i<nb_points;i++)
164     {
165     Points_Centre_au_barycentre[i]=lst_points[i]-barycentre;
166     }
167    
168 souaissa 71 OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;
169 souaissa 69
170 souaissa 71 for( int i=0;i<nb_points;i++)
171     {
172     OT_VECTEUR_4DD pt=Points_Centre_au_barycentre[i];
173     PT_CENTRE(i,0)= pt[0] ;
174     PT_CENTRE(i,1)= pt[1] ;
175     PT_CENTRE(i,2)= pt[2] ;
176     PT_CENTRE(i,3)= pt[3] ;
177     }
178     PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();
179    
180     COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
181     double2 COEF= 1./nb_points;
182     COVARIANCE=COVARIANCE*COEF;
183    
184 souaissa 69 return COVARIANCE;
185 souaissa 66 }
186    
187 souaissa 69
188 souaissa 71 void VCT_SURFACE::calcule_axes_dinertie(OT_VECTEUR_4DD& D,OT_TENSEUR& V)
189 souaissa 66 {
190 souaissa 69 int n=4;
191     int nrot;
192    
193 souaissa 71 OT_TENSEUR COV=calcule_covariance();
194    
195    
196 souaissa 69 double2 zro=0.0;
197    
198 souaissa 71 if(COV(0,0)==zro&&COV(1,0)==zro&&COV(2,0)==zro&&COV(3,0)==zro)
199 souaissa 69 {
200 souaissa 71 COV(0,0)=1.0;
201     COV(1,0)=zro;
202     COV(2,0)=zro;
203     COV(3,0)=zro;
204 souaissa 69 }
205 souaissa 71 if(COV(0,1)==zro&&COV(1,1)==zro&&COV(2,1)==zro&&COV(3,1)==zro)
206 souaissa 69 {
207 souaissa 71 COV(0,1)=zro;
208     COV(1,1)=1.0;
209     COV(2,1)=zro;
210     COV(3,1)=zro;
211 souaissa 69 }
212 souaissa 71 if(COV(0,2)==zro&&COV(1,2)==zro&&COV(2,2)==zro&&COV(3,2)==zro)
213 souaissa 69 {
214 souaissa 71 COV(0,2)=zro;
215     COV(1,2)=zro;
216     COV(2,2)=1.0;
217     COV(3,2)=zro;
218 souaissa 69 }
219 souaissa 71 if(COV(0,3)==zro&&COV(1,3)==zro&&COV(2,3)==zro&&COV(3,3)==zro)
220 souaissa 69 {
221 souaissa 71 COV(0,3)=zro;
222     COV(1,3)=zro;
223     COV(2,3)=zro;
224     COV(3,3)=1.0;
225 souaissa 69 }
226    
227    
228 souaissa 71 COV.get_orthogonalisation(COV,D,V,n,nrot);
229 souaissa 69
230 souaissa 66 }
231    
232    
233 souaissa 69 OT_TENSEUR VCT_SURFACE:: calcule_tenseur_inertie_au_barycentre() //repere globale
234     {
235 souaissa 71 OT_TENSEUR INERTIE(4);
236 souaissa 69 OT_TENSEUR cov=calcule_covariance();
237     double2 nb_pts=nb_points;
238     double2 MOINS_UN=-1.0;
239 souaissa 71 INERTIE=cov*MOINS_UN;
240     INERTIE=INERTIE*nb_pts;
241 souaissa 69
242 souaissa 71 INERTIE(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
243     INERTIE(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
244     INERTIE(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
245     INERTIE(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
246    
247     return INERTIE;
248 souaissa 69 }
249    
250    
251     OT_TENSEUR VCT_SURFACE:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
252     {
253    
254 souaissa 71 OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre();
255     OT_TENSEUR TENS=calcule_tenseur_inertie_au_barycentre();
256 souaissa 69
257     OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;
258    
259     double2 a=ABC[0] ;
260     double2 b=ABC[1] ;
261     double2 c=ABC[2] ;
262     double2 d=ABC[3] ;
263    
264     TENS(0,0)=TENS(0,0)+b*b+c*c+d*d;
265     TENS(1,1)=TENS(1,1)+a*a+c*c+d*d;
266     TENS(2,2)=TENS(2,2)+a*a+b*b+d*d;
267 souaissa 71 TENS(3,3)=TENS(3,3)+a*a+b*b+c*c;
268 souaissa 69
269     TENS(0,1)=TENS(0,1)-nb_points*a*b;
270     TENS(0,2)=TENS(0,2)-nb_points*a*c;
271     TENS(0,3)=TENS(0,3)-nb_points*a*d;
272    
273     TENS(1,0)=TENS(1,0)-nb_points*b*a;
274     TENS(1,2)=TENS(1,2)-nb_points*b*c;
275     TENS(1,3)=TENS(1,3)-nb_points*b*d;
276    
277     TENS(2,0)=TENS(2,0)-nb_points*c*a;
278     TENS(2,1)=TENS(2,1)-nb_points*c*b;
279     TENS(2,3)=TENS(2,3)-nb_points*c*d;
280    
281     TENS(3,0)=TENS(3,0)-nb_points*d*a;
282     TENS(3,1)=TENS(3,1)-nb_points*d*b;
283     TENS(3,2)=TENS(3,2)-nb_points*d*c;
284    
285     return TENS;
286     }
287    
288    
289     OT_TENSEUR VCT_SURFACE::calcule_tenseur_inertie_base_locale()
290     {
291 souaissa 71
292 souaissa 69 OT_VECTEUR_4DD v1,v2,v3,v4;
293     OT_TENSEUR INERTIE;
294     OT_TENSEUR I_GLOBALE;
295 souaissa 71 OT_VECTEUR_4DD D;
296     OT_TENSEUR V(4);
297     this->calcule_axes_dinertie(D,V);
298 souaissa 69 INERTIE=this->calcule_tenseur_inertie_au_barycentre();
299    
300     I_GLOBALE= INERTIE;
301    
302     OT_TENSEUR P(4);
303     OT_TENSEUR I=I_GLOBALE;
304     OT_TENSEUR P_TRSPOSE;
305     OT_TENSEUR I_LOCALE;
306 souaissa 71 P=V;
307 souaissa 69 P_TRSPOSE=P.transpose();
308    
309     I_LOCALE=P_TRSPOSE*I;
310     I_LOCALE=I_LOCALE*P;
311    
312 souaissa 71 return I_LOCALE;
313 souaissa 69 }
314    
315    
316    
317    
318     void VCT_SURFACE::calcule_tenseur_inertie_base_entite( OT_VECTEUR_4DD& g1g2,OT_TENSEUR& tens2, VCT& vct_f)
319     {
320 souaissa 71 OT_VECTEUR_4DD Dv,Dw,G1,G2,G1G2;
321 souaissa 69 OT_TENSEUR P(4);
322     OT_TENSEUR pt;
323     OT_TENSEUR IV;
324 souaissa 71 OT_TENSEUR Q(4);
325 souaissa 69 OT_TENSEUR qt;
326     OT_TENSEUR R;
327     OT_TENSEUR Rt;
328     OT_TENSEUR IW;
329 souaissa 71 OT_TENSEUR MASSE_CONCENTRE(4);
330     OT_TENSEUR V(4),W(4);
331     this->calcule_axes_dinertie(Dv,V);
332     vct_f.calcule_axes_dinertie(Dw,W);
333 souaissa 69
334 souaissa 71 P=V;
335     Q=W;
336 souaissa 69
337 souaissa 71 IW=vct_f.calcule_tenseur_inertie_base_locale();
338 souaissa 69
339     pt=P.transpose();
340     qt=Q.transpose();
341     R=pt*Q;
342     Rt=qt*P;
343     IV=R*IW;
344     IV=IV*Rt;
345 souaissa 66
346    
347 souaissa 69 G1=this->calcule_barycentre();
348     G2=vct_f.calcule_barycentre();
349    
350     G1G2=G2-G1;
351    
352     OT_VECTEUR_4DD G1G2_LOC; //Calcul du vecteur G1G2 dans la base locale;
353     for(int i=0;i<4;i++)
354     {
355     for(int j=0;j<4;j++)
356     {
357     G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
358 souaissa 66 }
359 souaissa 69 }
360 souaissa 66
361 souaissa 69 MASSE_CONCENTRE(0,0)=nb_points*(G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
362     MASSE_CONCENTRE(0,1)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[1]);
363     MASSE_CONCENTRE(0,2)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);
364     MASSE_CONCENTRE(0,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);
365    
366     MASSE_CONCENTRE(1,0)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[0]);
367     MASSE_CONCENTRE(1,1)= nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
368     MASSE_CONCENTRE(1,2)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[2]);
369     MASSE_CONCENTRE(1,3)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);
370    
371     MASSE_CONCENTRE(2,0)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[0]);
372     MASSE_CONCENTRE(2,1)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[1]);
373     MASSE_CONCENTRE(2,2)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[3]*G1G2_LOC[3]);
374     MASSE_CONCENTRE(2,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);
375    
376     MASSE_CONCENTRE(3,0)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);
377     MASSE_CONCENTRE(3,1)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);
378     MASSE_CONCENTRE(3,2)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[3]);
379     MASSE_CONCENTRE(3,3)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]);
380    
381     IV=IV+MASSE_CONCENTRE;
382 souaissa 66 }
383    
384 souaissa 69
385     ostream& operator <<(ostream& os, VCT_SURFACE& vct_f)
386     {
387     vct_f.enregistrer(os) ;
388     return os;
389     }
390    
391 souaissa 71
392    
393    
394     void VCT_SURFACE::enregistrer(std::ostream& ost)
395     {
396    
397     ost<<"========================================"<<endl;
398     ost<<"% SURFACE: "<<endl;
399     ost<<"========================================"<<endl;
400     ost<<"POINTS_DE_CONTROLS: "<<endl;
401     ost<<"========================================"<<endl;
402    
403     for (unsigned int i=0;i< lst_points.size();i++)
404     {
405     OT_VECTEUR_4DD v= lst_points[i]; // Rmq: la precision nèest pas affocher, il faut rajouter l'affichage
406     ost<< v<<endl; // de la precision dans la classe doubleprecision
407     }
408    
409     ost<<endl;
410     ost<<"VECTORISATION: "<<endl;
411     ost<<"========================================"<<endl;
412     for (unsigned int i=0;i< lst_vecteurs.size();i++)
413     {
414     OT_VECTEUR_4DD v= lst_vecteurs[i];
415     ost<< v<<endl;
416     }
417    
418     ost<<endl<<endl;
419     ost<<"BARYCENTRE: "<<endl;
420     ost<<"========================================"<<endl;
421     OT_VECTEUR_4DD BARY= calcule_barycentre();
422     ost<<BARY<<endl;
423     ost<<endl<<endl;
424     ost<<"TENSEUR_METRIQUE: "<<endl;
425     ost<<"========================================"<<endl;
426     OT_TENSEUR TNS_MT= calcule_tenseur_metrique() ;
427     ost<< TNS_MT<<endl;
428     ost<<endl<<endl;
429     ost<<"TENSEUR_DE_COVARIANCE: "<<endl;
430     ost<<"========================================"<<endl;
431     OT_TENSEUR TNS_CV= calcule_covariance();
432     ost<< TNS_CV<<endl;
433     ost<<endl<<endl;
434     ost<<"AXES_D'INERTIE: "<<endl;
435     ost<<"========================================"<<endl;
436     OT_VECTEUR_4DD D;
437     OT_TENSEUR V(4);
438     calcule_axes_dinertie(D,V);
439     ost<< V<<endl;
440     ost<<endl<<endl;
441     ost<<"INERTIE_CALCULÉE_AU_BARYCENTRE: "<<endl;
442     ost<<"========================================"<<endl;
443     OT_TENSEUR I_BARY=calcule_tenseur_inertie_au_barycentre();
444     ost<<I_BARY<<endl;
445    
446     ost<<endl<<endl;
447     ost<<"INERTIE_TRANSPORTÉE_EN_UN_POINT: "<<endl;
448     ost<<"========================================"<<endl;
449     OT_TENSEUR I_TRANSP=calcule_tenseur_inertie_au_pt(BARY);
450     ost<<I_TRANSP<<endl;
451     ost<<endl<<endl;
452    
453     ost<<"INERTIE_CALCULÉ_DANS_LA_BASE_LOCALE: "<<endl;
454     ost<<"========================================"<<endl;
455     OT_TENSEUR I_BASE=calcule_tenseur_inertie_base_locale();
456     ost<<I_BASE<<endl;
457     ost<<endl<<endl;
458     }
459    
460