geometrie/src/vct_surface.cpp

#include "gestionversion.h"

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


#pragma hdrstop

#include "vct_surface.h"
#include "mg_surface.h"
#include <iomanip.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_SURFACE::VCT_SURFACE(MG_SURFACE* 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);
nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);
int nb_u=nurbs_params.get(3);
int nb_v=nurbs_params.get(4);
double2 ZERO=0.;
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);

OT_VECTEUR_4DD V1,V2,V;

double2 X1,Y1,Z1,w1,X2,Y2,Z2,w2,vx,vy,vz,w;

  for(int r=0;r<nb_v;r++)
    {
     for(int s=0;s<nb_u-1;s++)
       {
   
           V1[0]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s);
           V1[1]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 1);
           V1[2]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 2);
           V1[3]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * s + 3);
           lst_points.insert(lst_points.end(),V1);
           V2[0]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1));
           V2[1]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1) + 1);
           V2[2]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 4 * (s + 1) + 2);
           V2[3]=nurbs_params.get( indx_premier_ptctr+ r *4 * nb_u + 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);
           }

        }
       }

     for(int r=0;r<nb_v-1;r++)
       {
      for(int s=0;s<nb_u;s++)
       {
           V1[0]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s);
           V1[1]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 1);
           V1[2]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 2);
           V1[3]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 3);
           lst_points.insert(lst_points.end(),V1);
           V2[0]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 *  nb_u + 4 * s);
           V2[1]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 *  nb_u + 4 * s + 1);
           V2[2]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 *  nb_u + 4 * s + 2);
           V2[3]=nurbs_params.get(indx_premier_ptctr+ (r+1) *4 *  nb_u + 4 * s + 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_SURFACE::VCT_SURFACE(VCT_SURFACE& mdd):VCT_ELEMENT_GEOMETRIQUE(mdd.elem_geo)
{
lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;
}


VCT_SURFACE::~VCT_SURFACE()
{
}

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


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

OT_VECTEUR_4DD VCT_SURFACE::calcule_barycentre()
{

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

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

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


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


OT_TENSEUR  VCT_SURFACE::calcule_covariance(void)
{

 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 COVARIANCE(4);

     for(int r=0;r<4;r++) {
       for(int s=0;s<4;s++) {
             for(int i=0;i<nb_points;i++)  {
                 OT_VECTEUR_4DD v1=Points_Centre_au_barycentre[i];
                 for(int j=0;j<nb_points;j++) {
                    OT_VECTEUR_4DD v2=Points_Centre_au_barycentre[j];
                    COVARIANCE(r,s)= COVARIANCE(r,s)+v1[i]*v2[j];
                    }
               }
               COVARIANCE(r,s)=1./nb_points*COVARIANCE(r,s);
         }
      }
 return COVARIANCE;
}


void  VCT_SURFACE::calcule_axes_dinertie(OT_VECTEUR_4DD& vp1,OT_VECTEUR_4DD& vp2,OT_VECTEUR_4DD& vp3,OT_VECTEUR_4DD& vp4)
{
 int n=4;
 int nrot;
 OT_TENSEUR v(4);
 OT_TENSEUR cov=calcule_covariance();
 OT_VECTEUR_4DD d;

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

OT_VECTEUR_4DD BARYCENTRE=this->calcule_barycentre();
OT_TENSEUR TENS=this->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+d*d;

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_SURFACE::calcule_tenseur_inertie_base_locale()
{
// inertie calcul� dans le repere globale
 OT_VECTEUR_4DD v1,v2,v3,v4;
 OT_TENSEUR INERTIE;
 OT_TENSEUR I_GLOBALE;

 this->calcule_axes_dinertie(v1,v2,v3,v4);
 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;

 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= v1[i];
    if (j==1) P(i,j)= v2[i];
    if (j==2) P(i,j)= v3[i];
    if (j==3) P(i,j)= v4[i];
    }
  }

 P_TRSPOSE=P.transpose();

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

}


void VCT_SURFACE::calcule_tenseur_inertie_base_entite( OT_VECTEUR_4DD& g1g2,OT_TENSEUR& tens2, VCT& vct_f)
{
OT_VECTEUR_4DD v1,v2,v3,v4,w1,w2,w3,w4,G1,G2,G1G2;
OT_TENSEUR P(4);
OT_TENSEUR pt;
OT_TENSEUR IV;
OT_TENSEUR Q;
OT_TENSEUR qt;
OT_TENSEUR R;
OT_TENSEUR Rt;
OT_TENSEUR IW;

this->calcule_axes_dinertie(v1,v2,v3,v4);
vct_f.calcule_axes_dinertie(w1,w2,w3,w4);


 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= v1[i];
    if (j==1) P(i,j)= v2[i];
    if (j==2) P(i,j)= v3[i];
    if (j==3) P(i,j)= v4[i];
    }
  }

 for(int i=0;i<4;i++) {
    for(int j=0;j<4;j++)
    {
    if (j==0) P(i,j)= w1[i];
    if (j==1) P(i,j)= w2[i];
    if (j==2) P(i,j)= w3[i];
    if (j==3) P(i,j)= w4[i];
    }
  }

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

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


OT_TENSEUR MASSE_CONCENTRE(4);

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

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

 }

ostream& operator <<(ostream& os, VCT_SURFACE& vct_f)
{
vct_f.enregistrer(os) ;
return os;
}

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