geometrie/src/vct_point.cpp

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


#pragma hdrstop
#include<math.h>
#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_points_controle(void)
{
return lst_points;
}

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_points[i];
 }

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


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


OT_TENSEUR  VCT_POINT::calcule_covariance(void)
{
 OT_TENSEUR COVARIANCE(4);

 OT_VECTEUR_4DD barycentre=calcule_barycentre();


   OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;

    for( int i=0;i<nb_points;i++)
    {
    OT_VECTEUR_4DD pt=lst_points[i];
    OT_VECTEUR_4DD ptcent=lst_points[i]-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_POINT::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_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;
   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_POINT:: 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_POINT::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;
   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_POINT::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_POINT& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}

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


void VCT_POINT::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:	85
Committed:	Fri May 2 14:24:42 2008 UTC (17 years ago) by souaissa
Original Path:	*magic/lib/geometrie/geometrie/src/vct_point.cpp*
File size:	10135 byte(s)
Log Message:
#	User	Rev	Content
1	souaissa	66	#include "gestionversion.h"
2			//---------------------------------------------------------------------------
3
4
5			#pragma hdrstop
6	souaissa	85	#include<math.h>
7	souaissa	66	#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	85	OT_VECTEUR_4DD barycentre=calcule_barycentre();
92	souaissa	66
93
94	souaissa	85	OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;
95	souaissa	66
96	souaissa	85	for( int i=0;i<nb_points;i++)
97			{
98			OT_VECTEUR_4DD pt=lst_points[i];
99			OT_VECTEUR_4DD ptcent=lst_points[i]-barycentre;
100			PT_CENTRE(i,0)= ptcent[0];
101			PT_CENTRE(i,1)= ptcent[1];
102			PT_CENTRE(i,2)= ptcent[2];
103			PT_CENTRE(i,3)= ptcent[3];
104			}
105			PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();
106
107			COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
108			double2 COEF= 1./nb_points;
109			COVARIANCE=COVARIANCE*COEF;
110
111			for(int i=0;i<4;i++)
112			for(int j=0;j<4;j++)
113			{
114			if(fabs(COVARIANCE(i,j).get_x())<=1e-8)
115			COVARIANCE(i,j).set_x(0.);
116
117			}
118			return COVARIANCE;
119	souaissa	69	}
120
121
122	souaissa	71	void VCT_POINT::calcule_axes_dinertie(OT_VECTEUR_4DD& D,OT_TENSEUR& V)
123	souaissa	69	{
124			int n=4;
125			int nrot;
126
127	souaissa	71	OT_TENSEUR COV=calcule_covariance();
128
129
130	souaissa	69	double2 zro=0.0;
131
132	souaissa	71	if(COV(0,0)==zro&&COV(1,0)==zro&&COV(2,0)==zro&&COV(3,0)==zro)
133	souaissa	69	{
134	souaissa	71	COV(0,0)=1.0;
135			COV(1,0)=zro;
136			COV(2,0)=zro;
137			COV(3,0)=zro;
138	souaissa	69	}
139	souaissa	71	if(COV(0,1)==zro&&COV(1,1)==zro&&COV(2,1)==zro&&COV(3,1)==zro)
140	souaissa	69	{
141	souaissa	71	COV(0,1)=zro;
142			COV(1,1)=1.0;
143			COV(2,1)=zro;
144			COV(3,1)=zro;
145	souaissa	69	}
146	souaissa	71	if(COV(0,2)==zro&&COV(1,2)==zro&&COV(2,2)==zro&&COV(3,2)==zro)
147	souaissa	69	{
148	souaissa	71	COV(0,2)=zro;
149			COV(1,2)=zro;
150			COV(2,2)=1.0;
151			COV(3,2)=zro;
152	souaissa	69	}
153	souaissa	71	if(COV(0,3)==zro&&COV(1,3)==zro&&COV(2,3)==zro&&COV(3,3)==zro)
154	souaissa	69	{
155	souaissa	71	COV(0,3)=zro;
156			COV(1,3)=zro;
157			COV(2,3)=zro;
158			COV(3,3)=1.0;
159	souaissa	69	}
160
161
162	souaissa	71	COV.get_orthogonalisation(COV,D,V,n,nrot);
163	souaissa	69
164			}
165
166
167	souaissa	85	OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_barycentre() //repere globale
168			{
169			OT_TENSEUR INERTIE(4);
170			OT_TENSEUR cov=calcule_covariance();
171			double2 nb_pts=nb_points;
172			double2 MOINS_UN=-1.0;
173			INERTIE=cov*MOINS_UN;
174			INERTIE=INERTIE*nb_pts;
175
176			INERTIE(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
177			INERTIE(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
178			INERTIE(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
179			INERTIE(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
180			for(int i=0;i<4;i++)
181			for(int j=0;j<4;j++)
182			{
183			if(fabs(INERTIE(i,j).get_x())<=1e-12) INERTIE(i,j).set_x(0.);
184
185			}
186			return INERTIE;
187			}
188
189
190	souaissa	69	OT_TENSEUR VCT_POINT:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
191			{
192
193	souaissa	71	OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre();
194			OT_TENSEUR TENS=calcule_tenseur_inertie_au_barycentre();
195	souaissa	69
196			OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;
197
198			double2 a=ABC[0] ;
199			double2 b=ABC[1] ;
200			double2 c=ABC[2] ;
201			double2 d=ABC[3] ;
202
203			TENS(0,0)=TENS(0,0)+bb+cc+d*d;
204			TENS(1,1)=TENS(1,1)+aa+cc+d*d;
205			TENS(2,2)=TENS(2,2)+aa+bb+d*d;
206	souaissa	71	TENS(3,3)=TENS(3,3)+aa+bb+c*c;
207	souaissa	69
208			TENS(0,1)=TENS(0,1)-nb_pointsab;
209			TENS(0,2)=TENS(0,2)-nb_pointsac;
210			TENS(0,3)=TENS(0,3)-nb_pointsad;
211
212			TENS(1,0)=TENS(1,0)-nb_pointsba;
213			TENS(1,2)=TENS(1,2)-nb_pointsbc;
214			TENS(1,3)=TENS(1,3)-nb_pointsbd;
215
216			TENS(2,0)=TENS(2,0)-nb_pointsca;
217			TENS(2,1)=TENS(2,1)-nb_pointscb;
218			TENS(2,3)=TENS(2,3)-nb_pointscd;
219
220			TENS(3,0)=TENS(3,0)-nb_pointsda;
221			TENS(3,1)=TENS(3,1)-nb_pointsdb;
222			TENS(3,2)=TENS(3,2)-nb_pointsdc;
223
224			return TENS;
225			}
226
227
228			OT_TENSEUR VCT_POINT::calcule_tenseur_inertie_base_locale()
229			{
230	souaissa	71
231	souaissa	69	OT_VECTEUR_4DD v1,v2,v3,v4;
232			OT_TENSEUR INERTIE;
233			OT_TENSEUR I_GLOBALE;
234	souaissa	71	OT_VECTEUR_4DD D;
235			OT_TENSEUR V(4);
236			this->calcule_axes_dinertie(D,V);
237	souaissa	69	INERTIE=this->calcule_tenseur_inertie_au_barycentre();
238
239			I_GLOBALE= INERTIE;
240
241			OT_TENSEUR P(4);
242			OT_TENSEUR I=I_GLOBALE;
243			OT_TENSEUR P_TRSPOSE;
244			OT_TENSEUR I_LOCALE;
245	souaissa	71	P=V;
246	souaissa	69	P_TRSPOSE=P.transpose();
247
248			I_LOCALE=P_TRSPOSE*I;
249			I_LOCALE=I_LOCALE*P;
250	souaissa	85	for(int i=0;i<4;i++)
251			for(int j=0;j<4;j++)
252			{
253			if(fabs(I_LOCALE(i,j).get_x())<=1e-12) I_LOCALE(i,j).set_x(0.);
254			}
255	souaissa	71	return I_LOCALE;
256
257	souaissa	66	}
258
259
260
261
262	souaissa	79	OT_TENSEUR VCT_POINT::calcule_tenseur_inertie_base_entite(VCT& vct_f)
263	souaissa	69	{
264	souaissa	71	OT_VECTEUR_4DD Dv,Dw,G1,G2,G1G2;
265	souaissa	69	OT_TENSEUR P(4);
266			OT_TENSEUR pt;
267			OT_TENSEUR IV;
268	souaissa	71	OT_TENSEUR Q(4);
269	souaissa	69	OT_TENSEUR qt;
270			OT_TENSEUR R;
271			OT_TENSEUR Rt;
272			OT_TENSEUR IW;
273	souaissa	71	OT_TENSEUR MASSE_CONCENTRE(4);
274			OT_TENSEUR V(4),W(4);
275			this->calcule_axes_dinertie(Dv,V);
276			vct_f.calcule_axes_dinertie(Dw,W);
277	souaissa	69
278	souaissa	71	P=V;
279			Q=W;
280	souaissa	69
281	souaissa	71	IW=vct_f.calcule_tenseur_inertie_base_locale();
282	souaissa	69
283			pt=P.transpose();
284			qt=Q.transpose();
285			R=pt*Q;
286			Rt=qt*P;
287			IV=R*IW;
288			IV=IV*Rt;
289
290
291			G1=this->calcule_barycentre();
292			G2=vct_f.calcule_barycentre();
293
294			G1G2=G2-G1;
295
296			OT_VECTEUR_4DD G1G2_LOC; //Calcul du vecteur G1G2 dans la base locale;
297			for(int i=0;i<4;i++)
298			{
299			for(int j=0;j<4;j++)
300			{
301			G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
302			}
303			}
304
305			MASSE_CONCENTRE(0,0)=nb_points(G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[2]G1G2_LOC[2]+G1G2_LOC[3]G1G2_LOC[3]);
306			MASSE_CONCENTRE(0,1)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[1]);
307			MASSE_CONCENTRE(0,2)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[2]);
308			MASSE_CONCENTRE(0,3)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
309
310			MASSE_CONCENTRE(1,0)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[0]);
311			MASSE_CONCENTRE(1,1)= nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[2]G1G2_LOC[2]+G1G2_LOC[3]G1G2_LOC[3]);
312			MASSE_CONCENTRE(1,2)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[2]);
313			MASSE_CONCENTRE(1,3)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
314
315			MASSE_CONCENTRE(2,0)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[0]);
316			MASSE_CONCENTRE(2,1)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[1]);
317			MASSE_CONCENTRE(2,2)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[3]G1G2_LOC[3]);
318	souaissa	79	MASSE_CONCENTRE(2,3)=-1nb_points(G1G2_LOC[3]*G1G2_LOC[2]);
319	souaissa	69
320			MASSE_CONCENTRE(3,0)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
321			MASSE_CONCENTRE(3,1)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
322			MASSE_CONCENTRE(3,2)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[3]);
323			MASSE_CONCENTRE(3,3)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[2]G1G2_LOC[2]);
324
325			IV=IV+MASSE_CONCENTRE;
326	souaissa	85
327
328	souaissa	79	return IV;
329	souaissa	85
330	souaissa	69	}
331
332
333	souaissa	71	ostream& operator <<(ostream& os, VCT_POINT& vct_f)
334			{
335
336			vct_f.enregistrer(os) ;
337			return os;
338	souaissa	69	}
339	souaissa	66
340	souaissa	71	//==========================================================================
341			//Visualisation
342			//==========================================================================
343
344
345
346
347			void VCT_POINT::enregistrer(std::ostream& ost)
348			{
349
350			ost<<"========================================"<<endl;
351			ost<<"% FACE: "<<endl;
352			ost<<"========================================"<<endl;
353			ost<<"POINTS_DE_CONTROLS: "<<endl;
354			ost<<"========================================"<<endl;
355
356			for (unsigned int i=0;i< lst_points.size();i++)
357			{
358			OT_VECTEUR_4DD v= lst_points[i]; // Rmq: la precision n�est pas affocher, il faut rajouter l'affichage
359			ost<< v<<endl; // de la precision dans la classe doubleprecision
360			}
361
362			ost<<endl;
363			ost<<"VECTORISATION: "<<endl;
364			ost<<"========================================"<<endl;
365			for (unsigned int i=0;i< lst_vecteurs.size();i++)
366			{
367			OT_VECTEUR_4DD v= lst_vecteurs[i];
368			ost<< v<<endl;
369			}
370
371			ost<<endl<<endl;
372			ost<<"BARYCENTRE: "<<endl;
373			ost<<"========================================"<<endl;
374			OT_VECTEUR_4DD BARY= calcule_barycentre();
375			ost<<BARY<<endl;
376			ost<<endl<<endl;
377			ost<<"TENSEUR_METRIQUE: "<<endl;
378			ost<<"========================================"<<endl;
379			OT_TENSEUR TNS_MT= calcule_tenseur_metrique() ;
380			ost<< TNS_MT<<endl;
381			ost<<endl<<endl;
382			ost<<"TENSEUR_DE_COVARIANCE: "<<endl;
383			ost<<"========================================"<<endl;
384			OT_TENSEUR TNS_CV= calcule_covariance();
385			ost<< TNS_CV<<endl;
386			ost<<endl<<endl;
387			ost<<"AXES_D'INERTIE: "<<endl;
388			ost<<"========================================"<<endl;
389			OT_VECTEUR_4DD D;
390			OT_TENSEUR V(4);
391			calcule_axes_dinertie(D,V);
392			ost<< V<<endl;
393			ost<<endl<<endl;
394			ost<<"INERTIE_CALCUL�E_AU_BARYCENTRE: "<<endl;
395			ost<<"========================================"<<endl;
396			OT_TENSEUR I_BARY=calcule_tenseur_inertie_au_barycentre();
397			ost<<I_BARY<<endl;
398
399			ost<<endl<<endl;
400			ost<<"INERTIE_TRANSPORT�E_EN_UN_POINT: "<<endl;
401			ost<<"========================================"<<endl;
402			OT_TENSEUR I_TRANSP=calcule_tenseur_inertie_au_pt(BARY);
403			ost<<I_TRANSP<<endl;
404			ost<<endl<<endl;
405
406			ost<<"INERTIE_CALCUL�_DANS_LA_BASE_LOCALE: "<<endl;
407			ost<<"========================================"<<endl;
408			OT_TENSEUR I_BASE=calcule_tenseur_inertie_base_locale();
409			ost<<I_BASE<<endl;
410			ost<<endl<<endl;
411
412
413			// void calcule_tenseur_inertie_base_entite(OT_VECTEUR_4DD& G1G2, OT_TENSEUR& tens,VCT& vct_f) ;
414
415			}