geometrie/src/vct_courbe.cpp

#include "gestionversion.h"

//---------------------------------------------------------------------------


#pragma hdrstop

#include "vct_courbe.h"
#include "mg_courbe.h"
#include "mg_arete.h"
#include <iomanip.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_COURBE::VCT_COURBE(MG_COURBE* elemgeo):VCT_ELEMENT_GEOMETRIQUE(elemgeo)
{
int indx_premier_ptctr;
OT_VECTEUR_4DD V1,V2,V;
TPL_LISTE_ENTITE<double> nurbs_params;
elem_geo->get_param_NURBS(indx_premier_ptctr,nurbs_params);
nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);
double2 ZERO=0.;
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);


int nb_u=nurbs_params.get(3);


    for(int s=0;s<nb_u-1;s++)
       {
           V1[0]=nurbs_params.get(indx_premier_ptctr+  4 * s);
           V1[1]=nurbs_params.get(indx_premier_ptctr+  4 * s + 1);
           V1[2]=nurbs_params.get(indx_premier_ptctr+  4 * s + 2);
           V1[3]=nurbs_params.get(indx_premier_ptctr+  4 * s + 3);
           if (s==0)
           lst_points.insert(lst_points.end(),V1);

           V2[0]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1));
           V2[1]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 1);
           V2[2]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 2);
           V2[3]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 3);

           lst_points.insert(lst_points.end(),V2);

           V=V2-V1;

           double2 norm_au_carre=(V[0]*V[0])+(V[1]*V[1])+(V[2]*V[2])+(V[3]*V[3]);
           double2 norm=norm_au_carre^0.5;

           if(V==VCT_NUL)
           lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);

           if(norm!=ZERO){
           V=1./norm*V;
           lst_vecteurs.insert(lst_vecteurs.end(),V);
           }


        }
}


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


VCT_COURBE::~VCT_COURBE()
{
}


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

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


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

OT_VECTEUR_4DD VCT_COURBE::calcule_barycentre()
{

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

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

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


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


OT_TENSEUR  VCT_COURBE::calcule_covariance(void)
{
 OT_TENSEUR COVARIANCE(4);
 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 PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;

    for( int i=0;i<nb_points;i++)
    {
    OT_VECTEUR_4DD pt=Points_Centre_au_barycentre[i];
    PT_CENTRE(i,0)= pt[0] ;
    PT_CENTRE(i,1)= pt[1] ;
    PT_CENTRE(i,2)= pt[2] ;
    PT_CENTRE(i,3)= pt[3] ;
    }
   PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();

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

 return COVARIANCE;
}


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

 OT_TENSEUR COV=calcule_covariance();


 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_COURBE:: 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_COURBE:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
{

OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre();
OT_TENSEUR TENS=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+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_COURBE::calcule_tenseur_inertie_base_locale()
{

 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(D,V);
 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;
 P=V;
 P_TRSPOSE=P.transpose();

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

 return I_LOCALE;

}


OT_TENSEUR VCT_COURBE::calcule_tenseur_inertie_base_entite(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(Dv,V);
vct_f.calcule_axes_dinertie(Dw,W);

P=V;
Q=W;

IW=vct_f.calcule_tenseur_inertie_base_locale();

pt=P.transpose();
qt=Q.transpose();
R=pt*Q;
Rt=qt*P;
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];
    }
  }

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_COURBE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}

//==========================================================================
//Visualisation
//==========================================================================


void VCT_COURBE::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();
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();
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(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();
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(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();
ost<<I_BASE<<endl;
ost<<endl<<endl;


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


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