geometrie/src/vct_sommet.cpp

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


#pragma hdrstop

#include "vct_sommet.h"
#include "mg_sommet.h"
#include "mg_point.h"
#include "vct_point.h"
#include <iomanip>
#include<math.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_SOMMET::VCT_SOMMET(MG_SOMMET* sommet):VCT_ELEMENT_TOPOLOGIQUE(sommet)
{

int indx_premier_ptctr;
MG_POINT* point=sommet->get_point();

TPL_LISTE_ENTITE<double> nurbs_params;
point->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_SOMMET::VCT_SOMMET(VCT_SOMMET& mdd):VCT_ELEMENT_TOPOLOGIQUE(mdd.elem_topo)
{
lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;
}


VCT_SOMMET::~ VCT_SOMMET()
{

}


std::vector<OT_VECTEUR_4DD>& VCT_SOMMET::get_points_controle(void)
{
return lst_points;
}

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


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

OT_VECTEUR_4DD VCT_SOMMET::calcule_barycentre(OT_TENSEUR& TT)
{
T=TT;
OT_VECTEUR_4DD barycentre(0,0,0,0);

for(int i=0;i<nb_points;i++)
 {
 OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
 barycentre+=pt;
 }

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


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


OT_TENSEUR  VCT_SOMMET::calcule_covariance(OT_TENSEUR& T)
{
 OT_TENSEUR COVARIANCE(4);

 OT_VECTEUR_4DD barycentre=calcule_barycentre(T);


   OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;

    for( int i=0;i<nb_points;i++)
    {
     OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
    OT_VECTEUR_4DD ptcent=pt-barycentre;
    PT_CENTRE(i,0)= ptcent[0];
    PT_CENTRE(i,1)= ptcent[1];
    PT_CENTRE(i,2)= ptcent[2];
    PT_CENTRE(i,3)= ptcent[3];
    }
   PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();

   COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
   double2 COEF= 1./nb_points;
   COVARIANCE=COVARIANCE*COEF;

   for(int i=0;i<4;i++)
      for(int j=0;j<4;j++)
        {
          if(fabs(COVARIANCE(i,j).get_x())<=1e-8)
          COVARIANCE(i,j).set_x(0.);

        }
 return COVARIANCE;
}


void  VCT_SOMMET::calcule_axes_dinertie(OT_TENSEUR& T,OT_VECTEUR_4DD& D,OT_TENSEUR& V)
{
 int n=4;
 int nrot;

 OT_TENSEUR COV=calcule_covariance(T);


 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_SOMMET:: calcule_tenseur_inertie_au_barycentre(OT_TENSEUR& T)  //repere globale
{
 OT_TENSEUR INERTIE(4);
 OT_TENSEUR cov=calcule_covariance(T);
 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;
   for(int i=0;i<4;i++)
      for(int j=0;j<4;j++)
        {
          if(fabs(INERTIE(i,j).get_x())<=1e-12) INERTIE(i,j).set_x(0.);

        }
 return INERTIE;
}


OT_TENSEUR VCT_SOMMET:: calcule_tenseur_inertie_au_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& POINT)
{

OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre(T);
OT_TENSEUR TENS=calcule_tenseur_inertie_au_barycentre(T);

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+c*c;

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_SOMMET::calcule_tenseur_inertie_base_locale(OT_TENSEUR& T)
{

 OT_VECTEUR_4DD v1,v2,v3,v4;
 OT_TENSEUR INERTIE;
 OT_TENSEUR I_GLOBALE;
 OT_VECTEUR_4DD D;
 OT_TENSEUR V(4);
 this->calcule_axes_dinertie(T,D,V);
 INERTIE=this->calcule_tenseur_inertie_au_barycentre(T);

 I_GLOBALE= INERTIE;

 OT_TENSEUR P(4);
 OT_TENSEUR I=I_GLOBALE;
 OT_TENSEUR P_TRSPOSE;
 OT_TENSEUR I_LOCALE;
 P=V;
 P_TRSPOSE=P.transpose();

 I_LOCALE=P_TRSPOSE*I;
 I_LOCALE=I_LOCALE*P;
   for(int i=0;i<4;i++)
      for(int j=0;j<4;j++)
        {
          if(fabs(I_LOCALE(i,j).get_x())<=1e-12) I_LOCALE(i,j).set_x(0.);
        }
 return I_LOCALE;

}


OT_TENSEUR VCT_SOMMET::calcule_tenseur_inertie_base_entite(OT_TENSEUR& T,VCT& vct_f)
{
OT_VECTEUR_4DD Dv,Dw,G1,G2,G1G2;
OT_TENSEUR P(4);
OT_TENSEUR pt;
OT_TENSEUR IV;
OT_TENSEUR Q(4);
OT_TENSEUR qt;
OT_TENSEUR R;
OT_TENSEUR Rt;
OT_TENSEUR IW;
OT_TENSEUR MASSE_CONCENTRE(4);
OT_TENSEUR V(4),W(4);
this->calcule_axes_dinertie(T,Dv,V);
vct_f.calcule_axes_dinertie(T,Dw,W);

P=V;
Q=W;

IW=vct_f.calcule_tenseur_inertie_base_locale(T);

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


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

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];
    }
  }

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[3]*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;


return IV;

}


ostream& operator <<(ostream& os, VCT_SOMMET& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}
//==========================================================================
//Visualisation
//==========================================================================

OT_VECTEUR_4DD VCT_SOMMET::get_nouveau_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& v)
{
OT_VECTEUR_4DD temp;

    for( int i=0;i<4;i++)
     { temp[i]=0.;
       for(int j=0;j<4;j++)
        temp[i]= temp[i]+T(i,j)*v[j];
     }
return temp;

}


void VCT_SOMMET::enregistrer(std::ostream& ost)
 {

ost<<"========================================"<<endl;
ost<<"% FACE:    "<<endl;
ost<<"========================================"<<endl;
ost<<"POINTS_DE_CONTROLS:    "<<endl;
ost<<"========================================"<<endl;

 for (unsigned int i=0;i< lst_points.size();i++)
   {
      OT_VECTEUR_4DD v= lst_points[i];   // Rmq: la precision n�est pas affocher, il faut rajouter l'affichage
        ost<< v<<endl;                   // de la precision dans la classe doubleprecision
   }

ost<<endl;
ost<<"VECTORISATION:    "<<endl;
ost<<"========================================"<<endl;
  for (unsigned int i=0;i< lst_vecteurs.size();i++)
   {
      OT_VECTEUR_4DD v= lst_vecteurs[i];
        ost<< v<<endl;
   }

ost<<endl<<endl;
ost<<"BARYCENTRE:    "<<endl;
ost<<"========================================"<<endl;
OT_VECTEUR_4DD BARY= calcule_barycentre(T);
ost<<BARY<<endl;
ost<<endl<<endl;
ost<<"TENSEUR_METRIQUE:    "<<endl;
ost<<"========================================"<<endl;
OT_TENSEUR TNS_MT= calcule_tenseur_metrique() ;
ost<< TNS_MT<<endl;
ost<<endl<<endl;
ost<<"TENSEUR_DE_COVARIANCE:    "<<endl;
ost<<"========================================"<<endl;
OT_TENSEUR TNS_CV=  calcule_covariance(T);
ost<< TNS_CV<<endl;
ost<<endl<<endl;
ost<<"AXES_D'INERTIE:    "<<endl;
ost<<"========================================"<<endl;
OT_VECTEUR_4DD D;
OT_TENSEUR V(4);
calcule_axes_dinertie(T,D,V);
ost<< V<<endl;
ost<<endl<<endl;
ost<<"INERTIE_CALCUL�E_AU_BARYCENTRE:    "<<endl;
ost<<"========================================"<<endl;
OT_TENSEUR I_BARY=calcule_tenseur_inertie_au_barycentre(T);
ost<<I_BARY<<endl;

ost<<endl<<endl;
ost<<"INERTIE_TRANSPORT�E_EN_UN_POINT:    "<<endl;
ost<<"========================================"<<endl;
OT_TENSEUR I_TRANSP=calcule_tenseur_inertie_au_pt(T,BARY);
ost<<I_TRANSP<<endl;
ost<<endl<<endl;

ost<<"INERTIE_CALCUL�_DANS_LA_BASE_LOCALE:    "<<endl;
ost<<"========================================"<<endl;
OT_TENSEUR I_BASE=calcule_tenseur_inertie_base_locale(T);
ost<<I_BASE<<endl;
ost<<endl<<endl;


// void calcule_tenseur_inertie_base_entite(OT_VECTEUR_4DD& G1G2, OT_TENSEUR& tens,VCT& vct_f) ;


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