geometrie/src/vct_surface.cpp

#include "gestionversion.h"

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


#pragma hdrstop
#include<math.h>
#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;

for(int pt=0;pt< nb_points;pt++)
{
           V1[0]=nurbs_params.get( indx_premier_ptctr+  4 * pt);
           V1[1]=nurbs_params.get( indx_premier_ptctr+  4 * pt + 1);
           V1[2]=nurbs_params.get( indx_premier_ptctr+  4 * pt + 2);
           V1[3]=nurbs_params.get( indx_premier_ptctr+  4 * pt + 3);

           lst_points.insert(lst_points.end(),V1);
}


  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);


           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);


           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);


           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);


           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_points_controle(void)
{
return lst_points;
}

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

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


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


OT_TENSEUR  VCT_SURFACE::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_SURFACE::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_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;
   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_SURFACE:: 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_SURFACE::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_SURFACE::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_SURFACE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}


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

ost<<"========================================"<<endl;
ost<<"% SURFACE:    "<<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;
 }


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