geometrie/src/vct_surface.cpp

#include "gestionversion.h"

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


#pragma hdrstop
#include<math.h>
#include "vct_surface.h"
#include "mg_surface.h"
#include <iomanip>
#include "vct_outils.h"
#include "constantegeo.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;
TPL_LISTE_ENTITE<double> param_surf;
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);

for(int pt=0;pt< nb_points;pt++)
        {
        OT_VECTEUR_4DD V1;
        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 j=0;j<nb_v-1;j++)
        for (int i=0;i<nb_u-1;i++)
                {
                OT_VECTEUR_4DD P1=lst_points[j*nb_u+i];
                OT_VECTEUR_4DD P2=lst_points[j*nb_u+i+1];
                OT_VECTEUR_4DD P3=lst_points[(j+1)*nb_u+i];
                OT_VECTEUR_4DD V1=P2-P1;
                //V1.vecteur_norme();
                OT_VECTEUR_4DD V2=P3-P1;
                //V2.vecteur_norme();
                lst_vecteurs.insert(lst_vecteurs.end(),V1);
                lst_vecteurs.insert(lst_vecteurs.end(),V2);
                }
for (int i=0;i<nb_u-1;i++)
                {
                OT_VECTEUR_4DD P1=lst_points[(nb_v-1)*nb_u+i];
                OT_VECTEUR_4DD P2=lst_points[(nb_v-1)*nb_u+i+1];
                OT_VECTEUR_4DD V1=P2-P1;
                //V1.vecteur_norme();
                lst_vecteurs.insert(lst_vecteurs.end(),V1);
                }
for (int i=0;i<nb_v-1;i++)
                {
                OT_VECTEUR_4DD P1=lst_points[i*nb_u+nb_u-1];
                OT_VECTEUR_4DD P2=lst_points[(i+1)*nb_u+nb_u-1];
                OT_VECTEUR_4DD V1=P2-P1;
                //V1.vecteur_norme();
                lst_vecteurs.insert(lst_vecteurs.end(),V1);
                }

/*double2 ZERO=0.;
double2 INF=1e6;
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
int type=0;
if(((MG_SURFACE*)elem_geo)->get_type_geometrique(param_surf)==MGCo_PLAN)
{
 type=1;
}
if(((MG_SURFACE*)elem_geo)->get_type_geometrique(param_surf)==MGCo_CYLINDRE)
{
 type=2;
}

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);
           //if (pt==0) lst_points.insert(lst_points.end(),V1);

           lst_points.insert(lst_points.end(),V1);
                          if(type==2)
                           {
                             if (pt==0|| pt==3||pt==7|| pt==10)
                              {
                               lst_points_axes.insert(lst_points_axes.end(),V1);
                              }
                           }
                          if(!type||type==1)
                          {
                          lst_points_axes.insert(lst_points_axes.end(),V1);

                          }
}


if(!type||type==1)
  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(!type||type==1)
               lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),VCT_NUL);
               }

           if(norm!=ZERO){
           V=1./norm*V;
           lst_vecteurs.insert(lst_vecteurs.end(),V);
           if(!type||type==1)
           lst_vecteurs_axes.insert(lst_vecteurs_axes.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(!type||type==1)
               lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),VCT_NUL);
               }

           if (norm!=ZERO){
           V=1./norm*V;
           lst_vecteurs.insert(lst_vecteurs.end(),V);
           if(!type||type==1)
           lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),V);

           }
           }
         }

   if(type==2)
     {
           vector<OT_VECTEUR_4DD> lst_pts;
           lst_pts.insert(lst_pts.end(),lst_points_axes[0] ) ;
           lst_pts.insert(lst_pts.end(),lst_points_axes[1] ) ;
           lst_pts.insert(lst_pts.end(),lst_points_axes[3] ) ;
           lst_pts.insert(lst_pts.end(),lst_points_axes[2] ) ;
              for(int r=0;r<lst_pts.size();r++)
                 {     OT_VECTEUR_4DD V;
                       if(r<lst_pts.size()-1)
                         V=lst_pts[r+1]-lst_pts[r];
                       if(r==lst_pts.size()-1)
                         V=lst_pts[0]-lst_pts[r];
                        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_axes.insert(lst_vecteurs_axes.end(),V);
                           }
                 }

                 OT_VECTEUR_4DD V=lst_vecteurs_axes[1];
                 lst_vecteurs_axes[1]=lst_vecteurs_axes[2];
                 lst_vecteurs_axes[2]=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_points_controle_pour_axes(void)
{
return lst_points_axes;
}
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_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_SURFACE::get_nb_points()
 {
 return nb_points  ;
 }


OT_TENSEUR  VCT_SURFACE::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=lst_points[i]-barycentre;
    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_SURFACE::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_SURFACE:: 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_SURFACE:: 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_SURFACE::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;
 return INERTIE;

}


OT_TENSEUR VCT_SURFACE::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;

}


ostream& operator <<(ostream& os, VCT_SURFACE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}


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

}


void VCT_SURFACE::enregistrer(std::ostream& ost)
 {
OT_TENSEUR TT(4);
for(int i=0;i<4;i++)
 TT(i,i)=1.0;
 T=TT;
 
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(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();
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();
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(T);
ost<<I_BASE<<endl;
ost<<endl<<endl;
 }


Revision:	222
Committed:	Thu Nov 19 19:14:16 2009 UTC (15 years, 5 months ago) by francois
Original Path:	*magic/lib/geometrie/geometrie/src/vct_surface.cpp*
File size:	17088 byte(s)
Log Message:	Divers corrections : vectorisation toxfem et extraction de peau.
#	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	francois	169	#include <iomanip>
11	souaissa	150	#include "vct_outils.h"
12			#include "constantegeo.h"
13	souaissa	66	//---------------------------------------------------------------------------
14
15			#pragma package(smart_init)
16
17
18			VCT_SURFACE::VCT_SURFACE(MG_SURFACE* elemgeo):VCT_ELEMENT_GEOMETRIQUE(elemgeo)
19			{
20	souaissa	69	int indx_premier_ptctr;
21			TPL_LISTE_ENTITE<double> nurbs_params;
22	souaissa	150	TPL_LISTE_ENTITE<double> param_surf;
23	souaissa	66	elem_geo->get_param_NURBS(indx_premier_ptctr,nurbs_params);
24			nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);
25	souaissa	69	int nb_u=nurbs_params.get(3);
26			int nb_v=nurbs_params.get(4);
27	francois	222
28			for(int pt=0;pt< nb_points;pt++)
29			{
30			OT_VECTEUR_4DD V1;
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			lst_points.insert(lst_points.end(),V1);
36			}
37
38			for (int j=0;j<nb_v-1;j++)
39			for (int i=0;i<nb_u-1;i++)
40			{
41			OT_VECTEUR_4DD P1=lst_points[j*nb_u+i];
42			OT_VECTEUR_4DD P2=lst_points[j*nb_u+i+1];
43			OT_VECTEUR_4DD P3=lst_points[(j+1)*nb_u+i];
44			OT_VECTEUR_4DD V1=P2-P1;
45			//V1.vecteur_norme();
46			OT_VECTEUR_4DD V2=P3-P1;
47			//V2.vecteur_norme();
48			lst_vecteurs.insert(lst_vecteurs.end(),V1);
49			lst_vecteurs.insert(lst_vecteurs.end(),V2);
50			}
51			for (int i=0;i<nb_u-1;i++)
52			{
53			OT_VECTEUR_4DD P1=lst_points[(nb_v-1)*nb_u+i];
54			OT_VECTEUR_4DD P2=lst_points[(nb_v-1)*nb_u+i+1];
55			OT_VECTEUR_4DD V1=P2-P1;
56			//V1.vecteur_norme();
57			lst_vecteurs.insert(lst_vecteurs.end(),V1);
58			}
59			for (int i=0;i<nb_v-1;i++)
60			{
61			OT_VECTEUR_4DD P1=lst_points[i*nb_u+nb_u-1];
62			OT_VECTEUR_4DD P2=lst_points[(i+1)*nb_u+nb_u-1];
63			OT_VECTEUR_4DD V1=P2-P1;
64			//V1.vecteur_norme();
65			lst_vecteurs.insert(lst_vecteurs.end(),V1);
66			}
67
68			/*double2 ZERO=0.;
69	souaissa	150	double2 INF=1e6;
70	souaissa	69	OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
71	souaissa	150	int type=0;
72			if(((MG_SURFACE*)elem_geo)->get_type_geometrique(param_surf)==MGCo_PLAN)
73			{
74			type=1;
75			}
76			if(((MG_SURFACE*)elem_geo)->get_type_geometrique(param_surf)==MGCo_CYLINDRE)
77			{
78			type=2;
79			}
80	souaissa	66
81	souaissa	69	OT_VECTEUR_4DD V1,V2,V;
82	souaissa	66
83	souaissa	79	for(int pt=0;pt< nb_points;pt++)
84			{
85			V1[0]=nurbs_params.get( indx_premier_ptctr+ 4 * pt);
86			V1[1]=nurbs_params.get( indx_premier_ptctr+ 4 * pt + 1);
87			V1[2]=nurbs_params.get( indx_premier_ptctr+ 4 * pt + 2);
88			V1[3]=nurbs_params.get( indx_premier_ptctr+ 4 * pt + 3);
89	souaissa	157	//if (pt==0) lst_points.insert(lst_points.end(),V1);
90	souaissa	66
91	souaissa	79	lst_points.insert(lst_points.end(),V1);
92	souaissa	150	if(type==2)
93			{
94			if (pt==0\|\| pt==3\|\|pt==7\|\| pt==10)
95			{
96			lst_points_axes.insert(lst_points_axes.end(),V1);
97			}
98			}
99			if(!type\|\|type==1)
100			{
101			lst_points_axes.insert(lst_points_axes.end(),V1);
102
103			}
104	souaissa	79	}
105
106
107
108	souaissa	150	if(!type\|\|type==1)
109	souaissa	66	for(int r=0;r<nb_v;r++)
110			{
111			for(int s=0;s<nb_u-1;s++)
112			{
113	souaissa	69
114			V1[0]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * s);
115			V1[1]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * s + 1);
116			V1[2]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * s + 2);
117			V1[3]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * s + 3);
118	souaissa	71
119
120	souaissa	79
121	souaissa	69	V2[0]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * (s + 1));
122			V2[1]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * (s + 1) + 1);
123			V2[2]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * (s + 1) + 2);
124			V2[3]=nurbs_params.get( indx_premier_ptctr+ r 4 nb_u + 4 * (s + 1) + 3);
125	souaissa	71
126	souaissa	66
127	souaissa	79
128	souaissa	69	V=V2-V1;
129	souaissa	66
130	souaissa	69	double2 norm_au_carre=(V[0]V[0])+(V[1]V[1])+(V[2]V[2])+(V[3]V[3]);
131			double2 norm=norm_au_carre^0.5;
132	souaissa	66
133	souaissa	150	if(V==VCT_NUL) {
134	souaissa	69	lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
135	souaissa	150	if(!type\|\|type==1)
136			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),VCT_NUL);
137			}
138	souaissa	69
139	souaissa	66	if(norm!=ZERO){
140	souaissa	69	V=1./norm*V;
141			lst_vecteurs.insert(lst_vecteurs.end(),V);
142	souaissa	150	if(!type\|\|type==1)
143			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),V);
144	souaissa	66	}
145
146			}
147			}
148
149			for(int r=0;r<nb_v-1;r++)
150			{
151			for(int s=0;s<nb_u;s++)
152			{
153	souaissa	69	V1[0]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s);
154			V1[1]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 1);
155			V1[2]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 2);
156			V1[3]=nurbs_params.get(indx_premier_ptctr+ r * 4 * nb_u + 4 * s + 3);
157	souaissa	71
158
159	souaissa	79
160	souaissa	69	V2[0]=nurbs_params.get(indx_premier_ptctr+ (r+1) 4 nb_u + 4 * s);
161			V2[1]=nurbs_params.get(indx_premier_ptctr+ (r+1) 4 nb_u + 4 * s + 1);
162			V2[2]=nurbs_params.get(indx_premier_ptctr+ (r+1) 4 nb_u + 4 * s + 2);
163			V2[3]=nurbs_params.get(indx_premier_ptctr+ (r+1) 4 nb_u + 4 * s + 3);
164	souaissa	71
165	souaissa	79
166	souaissa	69	V=V2-V1;
167	souaissa	66
168	souaissa	69	double2 norm_au_carre=(V[0]V[0])+(V[1]V[1])+(V[2]V[2])+(V[3]V[3]);
169			double2 norm=norm_au_carre^0.5;
170	souaissa	66
171	souaissa	150	if(V==VCT_NUL) {
172	souaissa	69	lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
173	souaissa	150	if(!type\|\|type==1)
174			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),VCT_NUL);
175			}
176	souaissa	66
177	souaissa	69	if (norm!=ZERO){
178			V=1./norm*V;
179			lst_vecteurs.insert(lst_vecteurs.end(),V);
180	souaissa	150	if(!type\|\|type==1)
181			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),V);
182	souaissa	66
183			}
184			}
185			}
186
187	souaissa	150	if(type==2)
188			{
189			vector<OT_VECTEUR_4DD> lst_pts;
190			lst_pts.insert(lst_pts.end(),lst_points_axes[0] ) ;
191			lst_pts.insert(lst_pts.end(),lst_points_axes[1] ) ;
192			lst_pts.insert(lst_pts.end(),lst_points_axes[3] ) ;
193			lst_pts.insert(lst_pts.end(),lst_points_axes[2] ) ;
194			for(int r=0;r<lst_pts.size();r++)
195			{ OT_VECTEUR_4DD V;
196			if(r<lst_pts.size()-1)
197			V=lst_pts[r+1]-lst_pts[r];
198			if(r==lst_pts.size()-1)
199			V=lst_pts[0]-lst_pts[r];
200			double2 norm_au_carre=(V[0]V[0])+(V[1]V[1])+(V[2]V[2])+(V[3]V[3]);
201			double2 norm=norm_au_carre^0.5;
202			if(V==VCT_NUL)
203			lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
204			if(norm!=ZERO)
205			{
206			V=1./norm*V;
207			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),V);
208			}
209			}
210	souaissa	69
211	souaissa	150	OT_VECTEUR_4DD V=lst_vecteurs_axes[1];
212			lst_vecteurs_axes[1]=lst_vecteurs_axes[2];
213			lst_vecteurs_axes[2]=V;
214	souaissa	69
215	souaissa	150	}
216
217	francois	222	*/
218	souaissa	66	}
219
220
221
222	souaissa	69
223
224
225			VCT_SURFACE::VCT_SURFACE(VCT_SURFACE& mdd):VCT_ELEMENT_GEOMETRIQUE(mdd.elem_geo)
226	souaissa	66	{
227	souaissa	69	lst_vecteurs=mdd.lst_vecteurs;
228			lst_points=mdd.lst_points;
229			nb_points=mdd.nb_points ;
230			}
231	souaissa	66
232
233	souaissa	69	VCT_SURFACE::~VCT_SURFACE()
234			{
235			}
236
237	souaissa	85	std::vector<OT_VECTEUR_4DD>& VCT_SURFACE::get_points_controle(void)
238	souaissa	69	{
239	souaissa	85	return lst_points;
240	souaissa	69	}
241	souaissa	150	std::vector<OT_VECTEUR_4DD>& VCT_SURFACE::get_points_controle_pour_axes(void)
242			{
243			return lst_points_axes;
244			}
245	souaissa	85	std::vector<OT_VECTEUR_4DD> & VCT_SURFACE::get_vecteurs()
246	souaissa	71	{
247	souaissa	85	return lst_vecteurs;
248	souaissa	71	}
249	souaissa	69
250	souaissa	71
251	souaissa	69	OT_TENSEUR VCT_SURFACE:: calcule_tenseur_metrique()
252			{
253			OT_TENSEUR tns(lst_vecteurs);
254			return tns;
255			}
256
257	souaissa	150	OT_VECTEUR_4DD VCT_SURFACE::calcule_barycentre(OT_TENSEUR& TT)
258	souaissa	69	{
259	souaissa	150	T=TT;
260	souaissa	69	OT_VECTEUR_4DD barycentre(0,0,0,0);
261
262	souaissa	66	for(int i=0;i<nb_points;i++)
263			{
264	souaissa	150	OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
265			//barycentre+=lst_points[i];
266			barycentre+= pt;
267	souaissa	66	}
268
269	souaissa	69	barycentre*=1./nb_points;
270	souaissa	150
271	souaissa	69	return barycentre;
272	souaissa	66	}
273
274
275
276	souaissa	69	int VCT_SURFACE::get_nb_points()
277			{
278			return nb_points ;
279			}
280
281
282	souaissa	150	OT_TENSEUR VCT_SURFACE::calcule_covariance(OT_TENSEUR& T)
283	souaissa	66	{
284	souaissa	71	OT_TENSEUR COVARIANCE(4);
285	souaissa	85
286	souaissa	150	OT_VECTEUR_4DD barycentre=calcule_barycentre(T);
287	souaissa	66
288	souaissa	69
289	souaissa	71	OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;
290	souaissa	69
291	souaissa	71	for( int i=0;i<nb_points;i++)
292			{
293	souaissa	150	//OT_VECTEUR_4DD pt=lst_points[i];
294			OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
295			//OT_VECTEUR_4DD ptcent=lst_points[i]-barycentre;
296			OT_VECTEUR_4DD ptcent=pt-barycentre;
297	souaissa	85	PT_CENTRE(i,0)= ptcent[0];
298			PT_CENTRE(i,1)= ptcent[1];
299			PT_CENTRE(i,2)= ptcent[2];
300			PT_CENTRE(i,3)= ptcent[3];
301	souaissa	71	}
302			PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();
303
304			COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
305			double2 COEF= 1./nb_points;
306			COVARIANCE=COVARIANCE*COEF;
307
308	souaissa	85	for(int i=0;i<4;i++)
309			for(int j=0;j<4;j++)
310			{
311			if(fabs(COVARIANCE(i,j).get_x())<=1e-8)
312			COVARIANCE(i,j).set_x(0.);
313
314			}
315	souaissa	69	return COVARIANCE;
316	souaissa	66	}
317
318	souaissa	69
319	souaissa	150	void VCT_SURFACE::calcule_axes_dinertie(OT_TENSEUR& T,OT_VECTEUR_4DD& D,OT_TENSEUR& V)
320	souaissa	66	{
321	souaissa	69	int n=4;
322			int nrot;
323
324	souaissa	150	OT_TENSEUR COV=calcule_covariance(T);
325	souaissa	71
326
327	souaissa	69	double2 zro=0.0;
328
329	souaissa	71	if(COV(0,0)==zro&&COV(1,0)==zro&&COV(2,0)==zro&&COV(3,0)==zro)
330	souaissa	69	{
331	souaissa	71	COV(0,0)=1.0;
332			COV(1,0)=zro;
333			COV(2,0)=zro;
334			COV(3,0)=zro;
335	souaissa	69	}
336	souaissa	71	if(COV(0,1)==zro&&COV(1,1)==zro&&COV(2,1)==zro&&COV(3,1)==zro)
337	souaissa	69	{
338	souaissa	71	COV(0,1)=zro;
339			COV(1,1)=1.0;
340			COV(2,1)=zro;
341			COV(3,1)=zro;
342	souaissa	69	}
343	souaissa	71	if(COV(0,2)==zro&&COV(1,2)==zro&&COV(2,2)==zro&&COV(3,2)==zro)
344	souaissa	69	{
345	souaissa	71	COV(0,2)=zro;
346			COV(1,2)=zro;
347			COV(2,2)=1.0;
348			COV(3,2)=zro;
349	souaissa	69	}
350	souaissa	71	if(COV(0,3)==zro&&COV(1,3)==zro&&COV(2,3)==zro&&COV(3,3)==zro)
351	souaissa	69	{
352	souaissa	71	COV(0,3)=zro;
353			COV(1,3)=zro;
354			COV(2,3)=zro;
355			COV(3,3)=1.0;
356	souaissa	69	}
357
358
359	souaissa	71	COV.get_orthogonalisation(COV,D,V,n,nrot);
360	souaissa	69
361	souaissa	66	}
362
363
364	souaissa	150	OT_TENSEUR VCT_SURFACE:: calcule_tenseur_inertie_au_barycentre(OT_TENSEUR& T) //repere globale
365	souaissa	69	{
366	souaissa	71	OT_TENSEUR INERTIE(4);
367	souaissa	150	OT_TENSEUR cov=calcule_covariance(T);
368	souaissa	69	double2 nb_pts=nb_points;
369			double2 MOINS_UN=-1.0;
370	souaissa	71	INERTIE=cov*MOINS_UN;
371			INERTIE=INERTIE*nb_pts;
372	souaissa	69
373	souaissa	71	INERTIE(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
374			INERTIE(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
375			INERTIE(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
376			INERTIE(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
377	souaissa	85	for(int i=0;i<4;i++)
378			for(int j=0;j<4;j++)
379			{
380			if(fabs(INERTIE(i,j).get_x())<=1e-12) INERTIE(i,j).set_x(0.);
381	souaissa	71
382	souaissa	85	}
383	souaissa	71	return INERTIE;
384	souaissa	150
385	souaissa	69	}
386
387
388	souaissa	150	OT_TENSEUR VCT_SURFACE:: calcule_tenseur_inertie_au_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& POINT)
389	souaissa	69	{
390
391	souaissa	150	OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre(T);
392			OT_TENSEUR TENS=calcule_tenseur_inertie_au_barycentre(T);
393	souaissa	69
394			OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;
395
396			double2 a=ABC[0] ;
397			double2 b=ABC[1] ;
398			double2 c=ABC[2] ;
399			double2 d=ABC[3] ;
400
401			TENS(0,0)=TENS(0,0)+bb+cc+d*d;
402			TENS(1,1)=TENS(1,1)+aa+cc+d*d;
403			TENS(2,2)=TENS(2,2)+aa+bb+d*d;
404	souaissa	71	TENS(3,3)=TENS(3,3)+aa+bb+c*c;
405	souaissa	69
406			TENS(0,1)=TENS(0,1)-nb_pointsab;
407			TENS(0,2)=TENS(0,2)-nb_pointsac;
408			TENS(0,3)=TENS(0,3)-nb_pointsad;
409
410			TENS(1,0)=TENS(1,0)-nb_pointsba;
411			TENS(1,2)=TENS(1,2)-nb_pointsbc;
412			TENS(1,3)=TENS(1,3)-nb_pointsbd;
413
414			TENS(2,0)=TENS(2,0)-nb_pointsca;
415			TENS(2,1)=TENS(2,1)-nb_pointscb;
416			TENS(2,3)=TENS(2,3)-nb_pointscd;
417
418			TENS(3,0)=TENS(3,0)-nb_pointsda;
419			TENS(3,1)=TENS(3,1)-nb_pointsdb;
420			TENS(3,2)=TENS(3,2)-nb_pointsdc;
421
422			return TENS;
423			}
424
425
426	souaissa	150	OT_TENSEUR VCT_SURFACE::calcule_tenseur_inertie_base_locale(OT_TENSEUR& T)
427	souaissa	69	{
428	souaissa	71
429	souaissa	69	OT_VECTEUR_4DD v1,v2,v3,v4;
430			OT_TENSEUR INERTIE;
431			OT_TENSEUR I_GLOBALE;
432	souaissa	71	OT_VECTEUR_4DD D;
433			OT_TENSEUR V(4);
434	souaissa	150	this->calcule_axes_dinertie(T,D,V);
435			INERTIE=this->calcule_tenseur_inertie_au_barycentre(T);
436	souaissa	69
437			I_GLOBALE= INERTIE;
438
439			OT_TENSEUR P(4);
440			OT_TENSEUR I=I_GLOBALE;
441			OT_TENSEUR P_TRSPOSE;
442			OT_TENSEUR I_LOCALE;
443	souaissa	71	P=V;
444	souaissa	69	P_TRSPOSE=P.transpose();
445
446			I_LOCALE=P_TRSPOSE*I;
447			I_LOCALE=I_LOCALE*P;
448	souaissa	85	for(int i=0;i<4;i++)
449			for(int j=0;j<4;j++)
450			{
451			if(fabs(I_LOCALE(i,j).get_x())<=1e-12) I_LOCALE(i,j).set_x(0.);
452			}
453			return I_LOCALE;
454	souaissa	150	return INERTIE;
455	souaissa	69
456			}
457
458
459
460
461	souaissa	150	OT_TENSEUR VCT_SURFACE::calcule_tenseur_inertie_base_entite(OT_TENSEUR& T,VCT& vct_f)
462	souaissa	69	{
463	souaissa	71	OT_VECTEUR_4DD Dv,Dw,G1,G2,G1G2;
464	souaissa	69	OT_TENSEUR P(4);
465			OT_TENSEUR pt;
466			OT_TENSEUR IV;
467	souaissa	71	OT_TENSEUR Q(4);
468	souaissa	69	OT_TENSEUR qt;
469			OT_TENSEUR R;
470			OT_TENSEUR Rt;
471			OT_TENSEUR IW;
472	souaissa	71	OT_TENSEUR MASSE_CONCENTRE(4);
473			OT_TENSEUR V(4),W(4);
474	souaissa	150	this->calcule_axes_dinertie(T,Dv,V);
475			vct_f.calcule_axes_dinertie(T,Dw,W);
476	souaissa	69
477	souaissa	71	P=V;
478			Q=W;
479	souaissa	69
480	souaissa	150	IW=vct_f.calcule_tenseur_inertie_base_locale(T);
481	souaissa	69
482			pt=P.transpose();
483			qt=Q.transpose();
484			R=pt*Q;
485			Rt=qt*P;
486			IV=R*IW;
487			IV=IV*Rt;
488	souaissa	66
489
490	souaissa	150	G1=this->calcule_barycentre(T);
491			G2=vct_f.calcule_barycentre(T);
492	souaissa	69
493			G1G2=G2-G1;
494
495			OT_VECTEUR_4DD G1G2_LOC; //Calcul du vecteur G1G2 dans la base locale;
496			for(int i=0;i<4;i++)
497			{
498			for(int j=0;j<4;j++)
499			{
500			G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
501	souaissa	66	}
502	souaissa	69	}
503	souaissa	66
504	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]);
505			MASSE_CONCENTRE(0,1)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[1]);
506			MASSE_CONCENTRE(0,2)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[2]);
507			MASSE_CONCENTRE(0,3)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
508
509			MASSE_CONCENTRE(1,0)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[0]);
510			MASSE_CONCENTRE(1,1)= nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[2]G1G2_LOC[2]+G1G2_LOC[3]G1G2_LOC[3]);
511			MASSE_CONCENTRE(1,2)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[2]);
512			MASSE_CONCENTRE(1,3)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
513
514			MASSE_CONCENTRE(2,0)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[0]);
515			MASSE_CONCENTRE(2,1)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[1]);
516			MASSE_CONCENTRE(2,2)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[3]G1G2_LOC[3]);
517	souaissa	79	MASSE_CONCENTRE(2,3)=-1nb_points(G1G2_LOC[3]*G1G2_LOC[2]);
518	souaissa	69
519			MASSE_CONCENTRE(3,0)=-1nb_points(G1G2_LOC[0]*G1G2_LOC[3]);
520			MASSE_CONCENTRE(3,1)=-1nb_points(G1G2_LOC[1]*G1G2_LOC[3]);
521			MASSE_CONCENTRE(3,2)=-1nb_points(G1G2_LOC[2]*G1G2_LOC[3]);
522			MASSE_CONCENTRE(3,3)=nb_points(G1G2_LOC[0]G1G2_LOC[0]+G1G2_LOC[1]G1G2_LOC[1]+G1G2_LOC[2]G1G2_LOC[2]);
523
524			IV=IV+MASSE_CONCENTRE;
525	souaissa	85
526
527	souaissa	79	return IV;
528	souaissa	85
529	souaissa	66	}
530
531	souaissa	69
532			ostream& operator <<(ostream& os, VCT_SURFACE& vct_f)
533			{
534	souaissa	85
535	souaissa	69	vct_f.enregistrer(os) ;
536			return os;
537			}
538
539	souaissa	71
540	souaissa	150	OT_VECTEUR_4DD VCT_SURFACE::get_nouveau_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& v)
541			{
542			OT_VECTEUR_4DD temp;
543	souaissa	71
544	souaissa	150	for( int i=0;i<4;i++)
545			{ temp[i]=0.;
546			for(int j=0;j<4;j++)
547			temp[i]= temp[i]+T(i,j)*v[j];
548			}
549			return temp;
550	souaissa	71
551	souaissa	150	}
552
553
554	souaissa	71	void VCT_SURFACE::enregistrer(std::ostream& ost)
555			{
556	souaissa	150	OT_TENSEUR TT(4);
557			for(int i=0;i<4;i++)
558			TT(i,i)=1.0;
559			T=TT;
560
561	souaissa	71	ost<<"========================================"<<endl;
562			ost<<"% SURFACE: "<<endl;
563			ost<<"========================================"<<endl;
564			ost<<"POINTS_DE_CONTROLS: "<<endl;
565			ost<<"========================================"<<endl;
566
567			for (unsigned int i=0;i< lst_points.size();i++)
568			{
569	francois	222	OT_VECTEUR_4DD v= lst_points[i]; // Rmq: la precision n�est pas affocher, il faut rajouter l'affichage
570	souaissa	71	ost<< v<<endl; // de la precision dans la classe doubleprecision
571			}
572
573	francois	222	ost<<endl;
574	souaissa	71	ost<<"VECTORISATION: "<<endl;
575			ost<<"========================================"<<endl;
576			for (unsigned int i=0;i< lst_vecteurs.size();i++)
577			{
578			OT_VECTEUR_4DD v= lst_vecteurs[i];
579			ost<< v<<endl;
580			}
581
582			ost<<endl<<endl;
583	francois	222	ost<<"BARYCENTRE: "<<endl;
584	souaissa	71	ost<<"========================================"<<endl;
585	souaissa	150	OT_VECTEUR_4DD BARY= calcule_barycentre(T);
586	souaissa	71	ost<<BARY<<endl;
587			ost<<endl<<endl;
588	francois	222
589	souaissa	71	ost<<"TENSEUR_METRIQUE: "<<endl;
590			ost<<"========================================"<<endl;
591			OT_TENSEUR TNS_MT= calcule_tenseur_metrique() ;
592			ost<< TNS_MT<<endl;
593			ost<<endl<<endl;
594	francois	222	/*
595	souaissa	71	ost<<"TENSEUR_DE_COVARIANCE: "<<endl;
596			ost<<"========================================"<<endl;
597			OT_TENSEUR TNS_CV= calcule_covariance();
598			ost<< TNS_CV<<endl;
599			ost<<endl<<endl;
600	souaissa	150	*/
601	souaissa	71	ost<<"AXES_D'INERTIE: "<<endl;
602			ost<<"========================================"<<endl;
603			OT_VECTEUR_4DD D;
604			OT_TENSEUR V(4);
605	souaissa	150	calcule_axes_dinertie(T,D,V);
606	souaissa	71	ost<< V<<endl;
607	souaissa	150	ost<<endl<<endl; /*
608	francois	222	ost<<"INERTIE_CALCUL�E_AU_BARYCENTRE: "<<endl;
609	souaissa	71	ost<<"========================================"<<endl;
610			OT_TENSEUR I_BARY=calcule_tenseur_inertie_au_barycentre();
611			ost<<I_BARY<<endl;
612	souaissa	150
613	souaissa	71	ost<<endl<<endl;
614	francois	222	ost<<"INERTIE_TRANSPORT�E_EN_UN_POINT: "<<endl;
615	souaissa	71	ost<<"========================================"<<endl;
616			OT_TENSEUR I_TRANSP=calcule_tenseur_inertie_au_pt(BARY);
617			ost<<I_TRANSP<<endl;
618	souaissa	150	ost<<endl<<endl; */
619	souaissa	71
620	francois	222	ost<<"INERTIE_CALCUL�_DANS_LA_BASE_LOCALE: "<<endl;
621	souaissa	71	ost<<"========================================"<<endl;
622	souaissa	150	OT_TENSEUR I_BASE=calcule_tenseur_inertie_base_locale(T);
623	souaissa	71	ost<<I_BASE<<endl;
624			ost<<endl<<endl;
625			}
626
627