geometrie/src/vct_face.cpp

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


#pragma hdrstop

#include "vct_face.h"
#include "mg_face.h"
#include "mg_surface.h"

#include "mg_arete.h"
#include "ot_mathematique.h"
#include <math.h>
#include <iomanip>
#include "sld_fonction.h"

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

#pragma package(smart_init)


VCT_FACE::VCT_FACE(MG_FACE* face):VCT_ELEMENT_TOPOLOGIQUE(face)
{

int indx;
nb_points=0;

           int nb_boucle=((MG_FACE*)elem_topo)->get_nb_mg_boucle();

            for(int j=0;j<nb_boucle;j++)
               {
                MG_BOUCLE* Boucle = ((MG_FACE*)elem_topo)->get_mg_boucle(j);
                int nbarete = Boucle->get_nb_mg_coarete();
                for (int w =0; w<nbarete;w++)
                        {
                        TPL_LISTE_ENTITE<double> lst_points_ctrls;
                        MG_COARETE* coArete = Boucle->get_mg_coarete(w);
                        MG_ARETE* arete = coArete->get_arete();
                        arete->get_param_NURBS(indx,lst_points_ctrls);
                        int nb_points_arete=arete->get_vectorisation().get_nb_points();
                        vector<OT_VECTEUR_4DD>& list_points=arete->get_vectorisation().get_points_controle();

                        for(int i=0;i< nb_points_arete;i++)
                          {
                          OT_VECTEUR_4DD POINT= list_points[i];
                          lst_points.insert(lst_points.end(),POINT);
                          }

                        std::vector<OT_VECTEUR_4DD>& list_vect=arete->get_vectorisation().get_vecteurs() ;

                        for(unsigned int i=0;i<list_vect.size();i++)
                            {
                             lst_vecteurs.insert(lst_vecteurs.end(),list_vect[i]);
                            }
                        nb_points+=nb_points_arete;
                      }

               }


}


VCT_FACE::VCT_FACE(VCT_FACE& mdd):VCT_ELEMENT_TOPOLOGIQUE(mdd.elem_topo)
{
lst_vecteurs=mdd.lst_vecteurs;
lst_points =mdd.lst_points;
nb_points=mdd.nb_points;
}


VCT_FACE::~ VCT_FACE()
{

}


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

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


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

OT_VECTEUR_4DD VCT_FACE::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_FACE::get_nb_points()
 {
 return nb_points  ;
 }


OT_TENSEUR  VCT_FACE::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_FACE::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_FACE:: 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_FACE:: 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_FACE::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_FACE::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_FACE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}

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


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

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;

}


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