ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/vct_point.cpp
Revision: 79
Committed: Thu Apr 10 15:14:15 2008 UTC (17 years, 1 month ago) by souaissa
Original Path: magic/lib/geometrie/geometrie/src/vct_point.cpp
File size: 8754 byte(s)
Log Message:

File Contents

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