geometrie/src/vct_arete.cpp

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


#pragma hdrstop
#include "vct_arete.h"
#include "vct_courbe.h"
#include "mg_courbe.h"
#include "mg_arete.h"
#include <iomanip.h>
#include<math.h>
#include "ot_mathematique.h"
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_ARETE::VCT_ARETE(MG_ARETE* arete):VCT_ELEMENT_TOPOLOGIQUE(arete)
{
double2 ZERO=0.;
int indx_premier_ptctr;
OT_VECTEUR_4DD V1,V2,V;

TPL_LISTE_ENTITE<double> nurbs_params;
arete->get_param_NURBS(indx_premier_ptctr,nurbs_params);
nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);

OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);

    for(int s=0;s<nb_points-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_ARETE::VCT_ARETE(VCT_ARETE& mdd):VCT_ELEMENT_TOPOLOGIQUE(mdd.elem_topo)
{

lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;

}


VCT_ARETE::~ VCT_ARETE()
{

}

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

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


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

OT_VECTEUR_4DD VCT_ARETE::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+=lst_points[i];
 barycentre+=pt;
 }

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


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


OT_TENSEUR  VCT_ARETE::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=lst_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_ARETE::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_ARETE:: 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_ARETE:: 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_ARETE::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_ARETE::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;

}

OT_VECTEUR_4DD VCT_ARETE::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;

}
ostream& operator <<(ostream& os, VCT_ARETE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}

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


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