geometrie/src/vct_point.cpp

#include "gestionversion.h"
//---------------------------------------------------------------------------


#pragma hdrstop

#include "vct_point.h"
#include "mg_point.h"
//---------------------------------------------------------------------------
#include <iomanip.h>
#pragma package(smart_init)


VCT_POINT::VCT_POINT(MG_POINT* elemgeo):VCT_ELEMENT_GEOMETRIQUE(elemgeo)
{
int indx_premier_ptctr;

TPL_LISTE_ENTITE<double> nurbs_params;
elem_geo->get_param_NURBS(indx_premier_ptctr,nurbs_params);

OT_VECTEUR_4DD V;
           V[0]=nurbs_params.get(indx_premier_ptctr+  0);
           V[1]=nurbs_params.get(indx_premier_ptctr+  1);
           V[2]=nurbs_params.get(indx_premier_ptctr+  2);
           V[3]=nurbs_params.get(indx_premier_ptctr+  3);
lst_points.insert(lst_points.end(),V);
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
lst_vecteurs.insert(lst_vecteurs.end(), VCT_NUL);

nb_points=1;
}


VCT_POINT::VCT_POINT(VCT_POINT& mdd):VCT_ELEMENT_GEOMETRIQUE(mdd.elem_geo)
{
lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;
}


VCT_POINT::~ VCT_POINT()
{

}


std::vector<OT_VECTEUR_4DD> &VCT_POINT::get_vecteurs()
{
return lst_vecteurs;
}


OT_TENSEUR VCT_POINT:: calcule_tenseur_metrique()
{
OT_TENSEUR tns(lst_vecteurs);
return tns;
}

OT_VECTEUR_4DD VCT_POINT::calcule_barycentre()
{

OT_VECTEUR_4DD barycentre(0,0,0,0);

for(int i=0;i<nb_points;i++)
 {
 barycentre+=lst_vecteurs[i];
 }

barycentre*=1./nb_points;
return barycentre;
}


int VCT_POINT::get_nb_points()
 {
 return nb_points  ;
 }


OT_TENSEUR  VCT_POINT::calcule_covariance(void)
{

 std::vector<OT_VECTEUR_4DD> Points_Centre_au_barycentre(nb_points);
 OT_VECTEUR_4DD barycentre=calcule_barycentre();

   for( int i=0;i<nb_points;i++)
    {
    Points_Centre_au_barycentre[i]=lst_points[i]-barycentre;
    }

 OT_TENSEUR COVARIANCE(4);

     for(int r=0;r<4;r++) {
       for(int s=0;s<4;s++) {
             for(int i=0;i<nb_points;i++)  {
                 OT_VECTEUR_4DD v1=Points_Centre_au_barycentre[i];
                 for(int j=0;j<nb_points;j++) {
                    OT_VECTEUR_4DD v2=Points_Centre_au_barycentre[j];
                    COVARIANCE(r,s)= COVARIANCE(r,s)+v1[i]*v2[j];
                    }
               }
               COVARIANCE(r,s)=1./nb_points*COVARIANCE(r,s);
         }
      }
 return COVARIANCE;
}


void  VCT_POINT::calcule_axes_dinertie(OT_VECTEUR_4DD& vp1,OT_VECTEUR_4DD& vp2,OT_VECTEUR_4DD& vp3,OT_VECTEUR_4DD& vp4)
{
 int n=4;
 int nrot;
 OT_TENSEUR v(4);
 OT_TENSEUR cov=calcule_covariance();
 OT_VECTEUR_4DD d;

 double2 zro=0.0;

 if(cov(0,0)==zro&&cov(1,0)==zro&&cov(2,0)==zro&&cov(3,0)==zro)
   {
   cov(0,0)=1.0;
   cov(1,0)=zro;
   cov(2,0)=zro;
   cov(3,0)=zro;
   }
 if(cov(0,1)==zro&&cov(1,1)==zro&&cov(2,1)==zro&&cov(3,1)==zro)
   {
   cov(0,1)=zro;
   cov(1,1)=1.0;
   cov(2,1)=zro;
   cov(3,1)=zro;
   }
  if(cov(0,2)==zro&&cov(1,2)==zro&&cov(2,2)==zro&&cov(3,2)==zro)
   {
   cov(0,2)=zro;
   cov(1,2)=zro;
   cov(2,2)=1.0;
   cov(3,2)=zro;
   }
  if(cov(0,3)==zro&&cov(1,3)==zro&&cov(2,3)==zro&&cov(3,3)==zro)
   {
   cov(0,3)=zro;
   cov(1,3)=zro;
   cov(2,3)=zro;
   cov(3,3)=1.0;
   }


cov.get_orthogonalisation(cov,d,v,n,nrot);

}


OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_barycentre()  //repere globale
{
 OT_TENSEUR Inertie(4);
 OT_TENSEUR cov=calcule_covariance();
 double2 nb_pts=nb_points;
 double2 MOINS_UN=-1.0;
 Inertie=cov*MOINS_UN;
 Inertie=Inertie*nb_pts;

 Inertie(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
 Inertie(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
 Inertie(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
 Inertie(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
 return Inertie;
}


OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
{

OT_VECTEUR_4DD BARYCENTRE=this->calcule_barycentre();
OT_TENSEUR TENS=this->calcule_tenseur_inertie_au_barycentre();

OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;

double2 a=ABC[0] ;
double2 b=ABC[1] ;
double2 c=ABC[2] ;
double2 d=ABC[3] ;

TENS(0,0)=TENS(0,0)+b*b+c*c+d*d;
TENS(1,1)=TENS(1,1)+a*a+c*c+d*d;
TENS(2,2)=TENS(2,2)+a*a+b*b+d*d;
TENS(3,3)=TENS(3,3)+a*a+b*b+d*d;

TENS(0,1)=TENS(0,1)-nb_points*a*b;
TENS(0,2)=TENS(0,2)-nb_points*a*c;
TENS(0,3)=TENS(0,3)-nb_points*a*d;

TENS(1,0)=TENS(1,0)-nb_points*b*a;
TENS(1,2)=TENS(1,2)-nb_points*b*c;
TENS(1,3)=TENS(1,3)-nb_points*b*d;

TENS(2,0)=TENS(2,0)-nb_points*c*a;
TENS(2,1)=TENS(2,1)-nb_points*c*b;
TENS(2,3)=TENS(2,3)-nb_points*c*d;

TENS(3,0)=TENS(3,0)-nb_points*d*a;
TENS(3,1)=TENS(3,1)-nb_points*d*b;
TENS(3,2)=TENS(3,2)-nb_points*d*c;

return TENS;
}


OT_TENSEUR VCT_POINT::calcule_tenseur_inertie_base_locale()
{
// inertie calcul� dans le repere globale
 OT_VECTEUR_4DD v1,v2,v3,v4;
 OT_TENSEUR INERTIE;
 OT_TENSEUR I_GLOBALE;

 this->calcule_axes_dinertie(v1,v2,v3,v4);
 INERTIE=this->calcule_tenseur_inertie_au_barycentre();

 I_GLOBALE= INERTIE;

 OT_TENSEUR P(4);
 OT_TENSEUR I=I_GLOBALE;
 OT_TENSEUR P_TRSPOSE;
 OT_TENSEUR I_LOCALE;

 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= v1[i];
    if (j==1) P(i,j)= v2[i];
    if (j==2) P(i,j)= v3[i];
    if (j==3) P(i,j)= v4[i];
    }
  }

 P_TRSPOSE=P.transpose();

 I_LOCALE=P_TRSPOSE*I;
 I_LOCALE=I_LOCALE*P;
  return I_LOCALE;

}


void VCT_POINT::calcule_tenseur_inertie_base_entite( OT_VECTEUR_4DD& g1g2,OT_TENSEUR& tens2, VCT& vct_f)
{
OT_VECTEUR_4DD v1,v2,v3,v4,w1,w2,w3,w4,G1,G2,G1G2;
OT_TENSEUR P(4);
OT_TENSEUR pt;
OT_TENSEUR IV;
OT_TENSEUR Q;
OT_TENSEUR qt;
OT_TENSEUR R;
OT_TENSEUR Rt;
OT_TENSEUR IW;

this->calcule_axes_dinertie(v1,v2,v3,v4);
vct_f.calcule_axes_dinertie(w1,w2,w3,w4);


 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= v1[i];
    if (j==1) P(i,j)= v2[i];
    if (j==2) P(i,j)= v3[i];
    if (j==3) P(i,j)= v4[i];
    }
  }

 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= w1[i];
    if (j==1) P(i,j)= w2[i];
    if (j==2) P(i,j)= w3[i];
    if (j==3) P(i,j)= w4[i];
    }
  }

pt=P.transpose();
qt=Q.transpose();
R=pt*Q;
Rt=qt*P;
IW=vct_f.calcule_tenseur_inertie_base_locale();

IV=R*IW;
IV=IV*Rt;


G1=this->calcule_barycentre();
G2=vct_f.calcule_barycentre();

G1G2=G2-G1;

OT_VECTEUR_4DD G1G2_LOC;  //Calcul du vecteur G1G2 dans la base locale;
  for(int i=0;i<4;i++)
  {
    for(int j=0;j<4;j++)
    {
    G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
    }
  }


OT_TENSEUR MASSE_CONCENTRE(4);

MASSE_CONCENTRE(0,0)=nb_points*(G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(0,1)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[1]);
MASSE_CONCENTRE(0,2)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);
MASSE_CONCENTRE(0,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);

MASSE_CONCENTRE(1,0)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[0]);
MASSE_CONCENTRE(1,1)=   nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(1,2)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[2]);
MASSE_CONCENTRE(1,3)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);

MASSE_CONCENTRE(2,0)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[0]);
MASSE_CONCENTRE(2,1)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[1]);
MASSE_CONCENTRE(2,2)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(2,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);

MASSE_CONCENTRE(3,0)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,1)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,2)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,3)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]);

IV=IV+MASSE_CONCENTRE;
}


void VCT_POINT::enregistrer(ostream& o)
{ /*
  int compt=0;
  o<<"%%"<<((MG_SOMMET*)elem_topo)->get_id()<<"="<<"(POINTS_CONTROL_SOMMET,"<< endl;
  o<<"%%"<<((MG_ARETE*)elem_topo)->get_id()<<"="<<"(POLYGONE_DE_CONTROL_SOMMET,"<< endl;
  compt=0;
  for(unsigned int i=0;i<LISTE_TOPO_VECTEUR_SOMMET.size();i++)
    {
      o<<" %  "<<setw(15)<<LISTE_TOPO_VECTEUR_SOMMET[i];//<< " % , ";
       compt++;
       if(compt==4) {
       o<<endl;
       compt=0;}
    } 
     o<<")" <<endl;
     */
}

Revision:	69
Committed:	Thu Mar 27 13:20:26 2008 UTC (17 years, 1 month ago) by souaissa
Original Path:	*magic/lib/geometrie/geometrie/src/vct_point.cpp*
File size:	8405 byte(s)
Log Message:	Mise à jour des classes de la vectorisation et des calsses: ot_mathematique,ot_tenseur,ot_doubleprecision dans outil
#	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
49	souaissa	69	std::vector<OT_VECTEUR_4DD> &VCT_POINT::get_vecteurs()
50			{
51			return lst_vecteurs;
52	souaissa	66	}
53
54
55	souaissa	69	OT_TENSEUR VCT_POINT:: calcule_tenseur_metrique()
56	souaissa	66	{
57	souaissa	69	OT_TENSEUR tns(lst_vecteurs);
58			return tns;
59			}
60	souaissa	66
61	souaissa	69	OT_VECTEUR_4DD VCT_POINT::calcule_barycentre()
62			{
63	souaissa	66
64	souaissa	69	OT_VECTEUR_4DD barycentre(0,0,0,0);
65
66	souaissa	66	for(int i=0;i<nb_points;i++)
67			{
68	souaissa	69	barycentre+=lst_vecteurs[i];
69	souaissa	66	}
70
71	souaissa	69	barycentre*=1./nb_points;
72			return barycentre;
73	souaissa	66	}
74
75
76	souaissa	69
77			int VCT_POINT::get_nb_points()
78			{
79			return nb_points ;
80			}
81
82
83			OT_TENSEUR VCT_POINT::calcule_covariance(void)
84	souaissa	66	{
85
86	souaissa	69	std::vector<OT_VECTEUR_4DD> Points_Centre_au_barycentre(nb_points);
87			OT_VECTEUR_4DD barycentre=calcule_barycentre();
88	souaissa	66
89	souaissa	69	for( int i=0;i<nb_points;i++)
90			{
91			Points_Centre_au_barycentre[i]=lst_points[i]-barycentre;
92			}
93	souaissa	66
94	souaissa	69	OT_TENSEUR COVARIANCE(4);
95	souaissa	66
96	souaissa	69	for(int r=0;r<4;r++) {
97			for(int s=0;s<4;s++) {
98			for(int i=0;i<nb_points;i++) {
99			OT_VECTEUR_4DD v1=Points_Centre_au_barycentre[i];
100			for(int j=0;j<nb_points;j++) {
101			OT_VECTEUR_4DD v2=Points_Centre_au_barycentre[j];
102			COVARIANCE(r,s)= COVARIANCE(r,s)+v1[i]*v2[j];
103			}
104			}
105			COVARIANCE(r,s)=1./nb_points*COVARIANCE(r,s);
106			}
107			}
108			return COVARIANCE;
109			}
110
111
112			void VCT_POINT::calcule_axes_dinertie(OT_VECTEUR_4DD& vp1,OT_VECTEUR_4DD& vp2,OT_VECTEUR_4DD& vp3,OT_VECTEUR_4DD& vp4)
113			{
114			int n=4;
115			int nrot;
116			OT_TENSEUR v(4);
117			OT_TENSEUR cov=calcule_covariance();
118			OT_VECTEUR_4DD d;
119
120			double2 zro=0.0;
121
122			if(cov(0,0)==zro&&cov(1,0)==zro&&cov(2,0)==zro&&cov(3,0)==zro)
123			{
124			cov(0,0)=1.0;
125			cov(1,0)=zro;
126			cov(2,0)=zro;
127			cov(3,0)=zro;
128			}
129			if(cov(0,1)==zro&&cov(1,1)==zro&&cov(2,1)==zro&&cov(3,1)==zro)
130			{
131			cov(0,1)=zro;
132			cov(1,1)=1.0;
133			cov(2,1)=zro;
134			cov(3,1)=zro;
135			}
136			if(cov(0,2)==zro&&cov(1,2)==zro&&cov(2,2)==zro&&cov(3,2)==zro)
137			{
138			cov(0,2)=zro;
139			cov(1,2)=zro;
140			cov(2,2)=1.0;
141			cov(3,2)=zro;
142			}
143			if(cov(0,3)==zro&&cov(1,3)==zro&&cov(2,3)==zro&&cov(3,3)==zro)
144			{
145			cov(0,3)=zro;
146			cov(1,3)=zro;
147			cov(2,3)=zro;
148			cov(3,3)=1.0;
149			}
150
151
152			cov.get_orthogonalisation(cov,d,v,n,nrot);
153
154			}
155
156
157			OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_barycentre() //repere globale
158			{
159			OT_TENSEUR Inertie(4);
160			OT_TENSEUR cov=calcule_covariance();
161			double2 nb_pts=nb_points;
162			double2 MOINS_UN=-1.0;
163			Inertie=cov*MOINS_UN;
164			Inertie=Inertie*nb_pts;
165
166			Inertie(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
167			Inertie(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
168			Inertie(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
169			Inertie(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
170			return Inertie;
171			}
172
173
174			OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
175			{
176
177			OT_VECTEUR_4DD BARYCENTRE=this->calcule_barycentre();
178			OT_TENSEUR TENS=this->calcule_tenseur_inertie_au_barycentre();
179
180			OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;
181
182			double2 a=ABC[0] ;
183			double2 b=ABC[1] ;
184			double2 c=ABC[2] ;
185			double2 d=ABC[3] ;
186
187			TENS(0,0)=TENS(0,0)+bb+cc+d*d;
188			TENS(1,1)=TENS(1,1)+aa+cc+d*d;
189			TENS(2,2)=TENS(2,2)+aa+bb+d*d;
190			TENS(3,3)=TENS(3,3)+aa+bb+d*d;
191
192			TENS(0,1)=TENS(0,1)-nb_pointsab;
193			TENS(0,2)=TENS(0,2)-nb_pointsac;
194			TENS(0,3)=TENS(0,3)-nb_pointsad;
195
196			TENS(1,0)=TENS(1,0)-nb_pointsba;
197			TENS(1,2)=TENS(1,2)-nb_pointsbc;
198			TENS(1,3)=TENS(1,3)-nb_pointsbd;
199
200			TENS(2,0)=TENS(2,0)-nb_pointsca;
201			TENS(2,1)=TENS(2,1)-nb_pointscb;
202			TENS(2,3)=TENS(2,3)-nb_pointscd;
203
204			TENS(3,0)=TENS(3,0)-nb_pointsda;
205			TENS(3,1)=TENS(3,1)-nb_pointsdb;
206			TENS(3,2)=TENS(3,2)-nb_pointsdc;
207
208			return TENS;
209			}
210
211
212			OT_TENSEUR VCT_POINT::calcule_tenseur_inertie_base_locale()
213			{
214			// inertie calcul� dans le repere globale
215			OT_VECTEUR_4DD v1,v2,v3,v4;
216			OT_TENSEUR INERTIE;
217			OT_TENSEUR I_GLOBALE;
218
219			this->calcule_axes_dinertie(v1,v2,v3,v4);
220			INERTIE=this->calcule_tenseur_inertie_au_barycentre();
221
222			I_GLOBALE= INERTIE;
223
224			OT_TENSEUR P(4);
225			OT_TENSEUR I=I_GLOBALE;
226			OT_TENSEUR P_TRSPOSE;
227			OT_TENSEUR I_LOCALE;
228
229			for(int i=0;i<4;i++) {
230			for(int j=0;j<4;j++)
231			{
232			if (j==0) P(i,j)= v1[i];
233			if (j==1) P(i,j)= v2[i];
234			if (j==2) P(i,j)= v3[i];
235			if (j==3) P(i,j)= v4[i];
236			}
237	souaissa	66	}
238	souaissa	69
239			P_TRSPOSE=P.transpose();
240
241			I_LOCALE=P_TRSPOSE*I;
242			I_LOCALE=I_LOCALE*P;
243			return I_LOCALE;
244
245	souaissa	66	}
246
247
248
249
250	souaissa	69	void VCT_POINT::calcule_tenseur_inertie_base_entite( OT_VECTEUR_4DD& g1g2,OT_TENSEUR& tens2, VCT& vct_f)
251			{
252			OT_VECTEUR_4DD v1,v2,v3,v4,w1,w2,w3,w4,G1,G2,G1G2;
253			OT_TENSEUR P(4);
254			OT_TENSEUR pt;
255			OT_TENSEUR IV;
256			OT_TENSEUR Q;
257			OT_TENSEUR qt;
258			OT_TENSEUR R;
259			OT_TENSEUR Rt;
260			OT_TENSEUR IW;
261
262			this->calcule_axes_dinertie(v1,v2,v3,v4);
263			vct_f.calcule_axes_dinertie(w1,w2,w3,w4);
264
265
266			for(int i=0;i<4;i++) {
267			for(int j=0;j<4;j++)
268			{
269			if (j==0) P(i,j)= v1[i];
270			if (j==1) P(i,j)= v2[i];
271			if (j==2) P(i,j)= v3[i];
272			if (j==3) P(i,j)= v4[i];
273			}
274			}
275
276			for(int i=0;i<4;i++) {
277			for(int j=0;j<4;j++)
278			{
279			if (j==0) P(i,j)= w1[i];
280			if (j==1) P(i,j)= w2[i];
281			if (j==2) P(i,j)= w3[i];
282			if (j==3) P(i,j)= w4[i];
283			}
284			}
285
286			pt=P.transpose();
287			qt=Q.transpose();
288			R=pt*Q;
289			Rt=qt*P;
290			IW=vct_f.calcule_tenseur_inertie_base_locale();
291
292			IV=R*IW;
293			IV=IV*Rt;
294
295
296			G1=this->calcule_barycentre();
297			G2=vct_f.calcule_barycentre();
298
299			G1G2=G2-G1;
300
301			OT_VECTEUR_4DD G1G2_LOC; //Calcul du vecteur G1G2 dans la base locale;
302			for(int i=0;i<4;i++)
303			{
304			for(int j=0;j<4;j++)
305			{
306			G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
307			}
308			}
309
310
311
312			OT_TENSEUR MASSE_CONCENTRE(4);
313
314			MASSE_CONCENTRE(0,0)=nb_points(G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[2]G1G2_LOC[2]+G1G2_LOC[3]G1G2_LOC[3]);
315			MASSE_CONCENTRE(0,1)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[1]);
316			MASSE_CONCENTRE(0,2)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[2]);
317			MASSE_CONCENTRE(0,3)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
318
319			MASSE_CONCENTRE(1,0)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[0]);
320			MASSE_CONCENTRE(1,1)= nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[2]G1G2_LOC[2]+G1G2_LOC[3]G1G2_LOC[3]);
321			MASSE_CONCENTRE(1,2)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[2]);
322			MASSE_CONCENTRE(1,3)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
323
324			MASSE_CONCENTRE(2,0)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[0]);
325			MASSE_CONCENTRE(2,1)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[1]);
326			MASSE_CONCENTRE(2,2)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[3]G1G2_LOC[3]);
327			MASSE_CONCENTRE(2,3)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[2]);
328
329			MASSE_CONCENTRE(3,0)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
330			MASSE_CONCENTRE(3,1)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
331			MASSE_CONCENTRE(3,2)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[3]);
332			MASSE_CONCENTRE(3,3)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[2]G1G2_LOC[2]);
333
334			IV=IV+MASSE_CONCENTRE;
335			}
336
337
338			void VCT_POINT::enregistrer(ostream& o)
339			{ /*
340	souaissa	66	int compt=0;
341	souaissa	69	o<<"%%"<<((MG_SOMMET*)elem_topo)->get_id()<<"="<<"(POINTS_CONTROL_SOMMET,"<< endl;
342			o<<"%%"<<((MG_ARETE*)elem_topo)->get_id()<<"="<<"(POLYGONE_DE_CONTROL_SOMMET,"<< endl;
343			compt=0;
344			for(unsigned int i=0;i<LISTE_TOPO_VECTEUR_SOMMET.size();i++)
345	souaissa	66	{
346	souaissa	69	o<<" % "<<setw(15)<<LISTE_TOPO_VECTEUR_SOMMET[i];//<< " % , ";
347	souaissa	66	compt++;
348			if(compt==4) {
349			o<<endl;
350			compt=0;}
351	souaissa	69	}
352			o<<")" <<endl;
353			*/
354			}
355	souaissa	66