geometrie/src/vct_face.cpp

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


#pragma hdrstop

#include "vct_face.h"
#include "mg_face.h"
#include "mg_surface.h"
#include "vct_outils.h"
#include "mg_arete.h"
#include "ot_mathematique.h"
#include <math.h>
#include <iomanip>
#include "sld_fonction.h"
#include "constantegeo.h"
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_FACE::VCT_FACE(MG_FACE* face):VCT_ELEMENT_TOPOLOGIQUE(face)
{
int nb_boucle=face->get_nb_mg_boucle();
for(int j=0;j<nb_boucle;j++)
        {       
        MG_BOUCLE* boucle = face->get_mg_boucle(j);
        int nbarete = boucle->get_nb_mg_coarete();
        for (int w =0; w<nbarete;w++)
                {
                        MG_COARETE* coarete = boucle->get_mg_coarete(w);
                        MG_ARETE* arete = coarete->get_arete();
                        int sens=coarete->get_orientation();
                        std::vector<OT_VECTEUR_4DD>& list_points=arete->get_vectorisation().get_points_controle();
                        for(int i=0;i< list_points.size();i++)
                                lst_points.insert(lst_points.end(),list_points[i]);
                        std::vector<OT_VECTEUR_4DD>& list_vect=arete->get_vectorisation().get_vecteurs() ;
                        for(int i=0;i< list_vect.size();i++)
                                lst_vecteurs.insert(lst_vecteurs.end(),(list_vect[i]*sens));
                        
                }
        }
            
nb_points=lst_points.size();

/*
int indx;
nb_points=0;
TPL_LISTE_ENTITE<double> param_surf;
TPL_LISTE_ENTITE<double> param_courb;
MG_SURFACE*surf=face->get_surface();
int type=0;
if(surf->get_type_geometrique(param_surf)==MGCo_PLAN)
{
 type=1;
}
if(surf->get_type_geometrique(param_surf)==MGCo_CYLINDRE)
{
 type=2;
}
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
double2 ZERO=0.;

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

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

                        for(int i=0;i< nb_points_arete;i++)
                          {

                          OT_VECTEUR_4DD POINT= list_points[i];
                          if(type==2)
                           {
                           if(courb->get_type_geometrique(param_courb)==MGCo_CIRCLE)
                             if (i==0|| i==3)
                              {
                               lst_points_axes.insert(lst_points_axes.end(),POINT);
                              }
                           }

                          lst_points.insert(lst_points.end(),POINT);
                         // if (w==0&&i==0&&nb_points_arete==2)
                          if(!type||type==1)
                          lst_points_axes.insert(lst_points_axes.end(),POINT);

                          }

                          std::vector<OT_VECTEUR_4DD>& list_vect=arete->get_vectorisation().get_vecteurs() ;
                          if(!type||type==1)
                            for(unsigned int i=0;i<list_vect.size();i++)
                             {
                             OT_VECTEUR_4DD v=list_vect[i]  ;
                             if (sens<0) v=-1*v;
                             lst_vecteurs.insert(lst_vecteurs.end(),v);
                             lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),v);
                             }
                        nb_points+=nb_points_arete;
                      }

               }
   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);
                           }
                 }
     }
*/
}


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


VCT_FACE::~ VCT_FACE()
{

}


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


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


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

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


OT_TENSEUR  VCT_FACE::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;
            for (int j=0;j<4;j++)
              PT_CENTRE(i,j)= ptcent[j];

    }
   PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();


   COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
   double2 COEF= 1./nb_points;
   COVARIANCE=COVARIANCE*COEF;

   for(int i=0;i<4;i++)
      for(int j=0;j<4;j++)
        {
          if(fabs(COVARIANCE(i,j).get_x())<=1e-8)
          COVARIANCE(i,j).set_x(0.);

        }


   return  COVARIANCE;
}


void  VCT_FACE::calcule_axes_dinertie(OT_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_VECTEUR_4DD G1;
       OT_TENSEUR p(4),q(4),pT,qT;

       G1=calcule_barycentre(T);
       //OT_VECTEUR_4DD VNul(0.0001,0.0001,0.0001,0.0001);
       double2 z=0.0001;
       double2 cof=-1;
       int plan=-1;

       pT=p.transpose();

       vector<OT_VECTEUR_4DD>& list_points1=get_points_controle_pour_axes();
       vector<OT_VECTEUR_4DD>& list_vecteurs1=lst_vecteurs_axes;//get_vecteurs();
       OT_VECTEUR_4DD Normal1(0.,0.,0.,0.),N1;
       OT_VECTEUR_4DD axe11,axe12,axe13;
       int cmpt=0;

       OT_VECTEUR_4DD v11= list_vecteurs1[0];
       int nb_vecteur1=list_vecteurs1.size();
       for(unsigned int pi=0;pi<nb_vecteur1-1;pi++)
           {
               OT_VECTEUR_4DD v2= list_vecteurs1[pi+1];
               if(pi==0)
                  Normal1=v11^v2;
               if (pi>0)
                  N1=v11^v2;
               OT_VECTEUR_4DD colineaire= Normal1^N1;
               if(colineaire[0].get_fabs()<z&&colineaire[1].get_fabs()<z&&colineaire[2].get_fabs()<z&&colineaire[3].get_fabs()<z)
                  {
                    cmpt++;
                  }
           }
                vector<OT_VECTEUR_4DD> axes1;
                if (cmpt== list_vecteurs1.size()-1)
                  {
                  plan=1;
                  axe11=list_vecteurs1[0];
                  axe12=Normal1^axe11;
                  axe13=Normal1;
                  axes1.insert(axes1.end(),axe11);
                  axes1.insert(axes1.end(),axe12);
                  axes1.insert(axes1.end(),axe13);
                  }

         VCT_OUTILS VOUTIL(4);
         OT_MATRICE_3D  sys_axes1=VOUTIL.get_system_axes(G1,axes1,list_points1) ;
         for(int i=0;i<4;i++)
            {
             for(int j=0;j<4;j++)
                {
                if(i<3&&j<3)
                V(i,j)= sys_axes1(i,j);
                if(i<3 &&j>2)
                V(i,j)=0.;
                if(i==3&&j==3)  V(i,j)=1.;
                }
           }   */
 }


OT_TENSEUR VCT_FACE:: 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_FACE:: 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_FACE::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(4);
 OT_TENSEUR I_LOCALE;
 P=V;
 //---------------
 OT_MATRICE_3D p1;
    for(int t=0;t<3;t++)
       for(int tt=0;tt<3;tt++)
         p1(t,tt)=P(t,tt).get_x();
  p1=p1.inverse();
    for(int t=0;t<3;t++)
       for(int tt=0;tt<3;tt++)
         P_TRSPOSE(t,tt)=p1(t,tt);
  P_TRSPOSE(3,3)=1.;
 //---------------
 //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_FACE::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_FACE& vct_f)
{

vct_f.enregistrer(os) ;
return os;
}


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

}
//==========================================================================
//Visualisation
//==========================================================================


void VCT_FACE::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<<"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<<"========================================"<<endl;
ost<<"POINTS_DE_CONTROLS_AXES:    "<<endl;
ost<<"========================================"<<endl;

 for (unsigned int i=0;i< lst_points_axes.size();i++)
   {
      OT_VECTEUR_4DD v= lst_points_axes[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;
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;


}


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_face.cpp*
File size:	18021 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			#pragma hdrstop
6
7			#include "vct_face.h"
8			#include "mg_face.h"
9			#include "mg_surface.h"
10	souaissa	150	#include "vct_outils.h"
11	souaissa	66	#include "mg_arete.h"
12			#include "ot_mathematique.h"
13	francois	102	#include <math.h>
14	souaissa	66	#include <iomanip>
15			#include "sld_fonction.h"
16	souaissa	150	#include "constantegeo.h"
17	souaissa	66	//---------------------------------------------------------------------------
18
19			#pragma package(smart_init)
20
21
22
23			VCT_FACE::VCT_FACE(MG_FACE* face):VCT_ELEMENT_TOPOLOGIQUE(face)
24			{
25	francois	222	int nb_boucle=face->get_nb_mg_boucle();
26			for(int j=0;j<nb_boucle;j++)
27			{
28			MG_BOUCLE* boucle = face->get_mg_boucle(j);
29			int nbarete = boucle->get_nb_mg_coarete();
30			for (int w =0; w<nbarete;w++)
31			{
32			MG_COARETE* coarete = boucle->get_mg_coarete(w);
33			MG_ARETE* arete = coarete->get_arete();
34			int sens=coarete->get_orientation();
35			std::vector<OT_VECTEUR_4DD>& list_points=arete->get_vectorisation().get_points_controle();
36			for(int i=0;i< list_points.size();i++)
37			lst_points.insert(lst_points.end(),list_points[i]);
38			std::vector<OT_VECTEUR_4DD>& list_vect=arete->get_vectorisation().get_vecteurs() ;
39			for(int i=0;i< list_vect.size();i++)
40			lst_vecteurs.insert(lst_vecteurs.end(),(list_vect[i]*sens));
41
42			}
43			}
44
45			nb_points=lst_points.size();
46	souaissa	66
47	francois	222	/*
48	souaissa	69	int indx;
49	souaissa	68	nb_points=0;
50	souaissa	150	TPL_LISTE_ENTITE<double> param_surf;
51			TPL_LISTE_ENTITE<double> param_courb;
52			MG_SURFACE*surf=face->get_surface();
53			int type=0;
54			if(surf->get_type_geometrique(param_surf)==MGCo_PLAN)
55			{
56			type=1;
57			}
58			if(surf->get_type_geometrique(param_surf)==MGCo_CYLINDRE)
59			{
60			type=2;
61			}
62			OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);
63			double2 ZERO=0.;
64	souaissa	66
65			int nb_boucle=((MG_FACE*)elem_topo)->get_nb_mg_boucle();
66	souaissa	69
67	souaissa	66	for(int j=0;j<nb_boucle;j++)
68			{
69			MG_BOUCLE* Boucle = ((MG_FACE*)elem_topo)->get_mg_boucle(j);
70			int nbarete = Boucle->get_nb_mg_coarete();
71			for (int w =0; w<nbarete;w++)
72			{
73	souaissa	150	double xyz[3];
74	souaissa	69	TPL_LISTE_ENTITE<double> lst_points_ctrls;
75	souaissa	66	MG_COARETE* coArete = Boucle->get_mg_coarete(w);
76	souaissa	69	MG_ARETE* arete = coArete->get_arete();
77	souaissa	150	int sens=coArete->get_orientation();
78	souaissa	69	arete->get_param_NURBS(indx,lst_points_ctrls);
79	souaissa	150	MG_COURBE*courb=arete->get_courbe();
80	souaissa	71	int nb_points_arete=arete->get_vectorisation().get_nb_points();
81	francois	72	vector<OT_VECTEUR_4DD>& list_points=arete->get_vectorisation().get_points_controle();
82
83	souaissa	71	for(int i=0;i< nb_points_arete;i++)
84			{
85	souaissa	150
86	souaissa	71	OT_VECTEUR_4DD POINT= list_points[i];
87	souaissa	150	if(type==2)
88			{
89			if(courb->get_type_geometrique(param_courb)==MGCo_CIRCLE)
90			if (i==0\|\| i==3)
91			{
92			lst_points_axes.insert(lst_points_axes.end(),POINT);
93			}
94			}
95
96	souaissa	71	lst_points.insert(lst_points.end(),POINT);
97	souaissa	150	// if (w==0&&i==0&&nb_points_arete==2)
98			if(!type\|\|type==1)
99			lst_points_axes.insert(lst_points_axes.end(),POINT);
100
101	souaissa	71	}
102	souaissa	67
103	souaissa	150	std::vector<OT_VECTEUR_4DD>& list_vect=arete->get_vectorisation().get_vecteurs() ;
104			if(!type\|\|type==1)
105			for(unsigned int i=0;i<list_vect.size();i++)
106			{
107			OT_VECTEUR_4DD v=list_vect[i] ;
108			if (sens<0) v=-1*v;
109			lst_vecteurs.insert(lst_vecteurs.end(),v);
110			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),v);
111			}
112	souaissa	71	nb_points+=nb_points_arete;
113	souaissa	66	}
114
115			}
116	souaissa	150	if(type==2)
117			{
118			vector<OT_VECTEUR_4DD> lst_pts;
119			lst_pts.insert(lst_pts.end(),lst_points_axes[0] ) ;
120			lst_pts.insert(lst_pts.end(),lst_points_axes[1] ) ;
121			lst_pts.insert(lst_pts.end(),lst_points_axes[3] ) ;
122			lst_pts.insert(lst_pts.end(),lst_points_axes[2] ) ;
123			for(int r=0;r<lst_pts.size();r++)
124			{ OT_VECTEUR_4DD V;
125			if(r<lst_pts.size()-1)
126			V=lst_pts[r+1]-lst_pts[r];
127			if(r==lst_pts.size()-1)
128			V=lst_pts[0]-lst_pts[r];
129			double2 norm_au_carre=(V[0]V[0])+(V[1]V[1])+(V[2]V[2])+(V[3]V[3]);
130			double2 norm=norm_au_carre^0.5;
131			if(V==VCT_NUL)
132			lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);
133			if(norm!=ZERO)
134			{
135			V=1./norm*V;
136			lst_vecteurs_axes.insert(lst_vecteurs_axes.end(),V);
137			}
138			}
139			}
140	francois	222	*/
141	souaissa	66	}
142
143
144	souaissa	69
145			VCT_FACE::VCT_FACE(VCT_FACE& mdd):VCT_ELEMENT_TOPOLOGIQUE(mdd.elem_topo)
146	souaissa	66	{
147	souaissa	69	lst_vecteurs=mdd.lst_vecteurs;
148			lst_points =mdd.lst_points;
149			nb_points=mdd.nb_points;
150	souaissa	66	}
151
152
153	souaissa	69	VCT_FACE::~ VCT_FACE()
154			{
155	souaissa	66
156	souaissa	69	}
157
158
159	francois	72	std::vector<OT_VECTEUR_4DD>& VCT_FACE::get_points_controle(void)
160	souaissa	71	{
161			return lst_points;
162			}
163	souaissa	150	std::vector<OT_VECTEUR_4DD>& VCT_FACE::get_points_controle_pour_axes(void)
164			{
165			return lst_points_axes;
166			}
167	souaissa	69
168	souaissa	150
169	souaissa	69	std::vector<OT_VECTEUR_4DD> & VCT_FACE::get_vecteurs()
170	souaissa	66	{
171	souaissa	69	return lst_vecteurs;
172			}
173	souaissa	66
174
175	souaissa	69	OT_TENSEUR VCT_FACE:: calcule_tenseur_metrique()
176			{
177			OT_TENSEUR tns(lst_vecteurs);
178			return tns;
179			}
180
181	souaissa	150	OT_VECTEUR_4DD VCT_FACE::calcule_barycentre(OT_TENSEUR& TT)
182	souaissa	69	{
183	souaissa	150	T=TT;
184	souaissa	69	OT_VECTEUR_4DD barycentre(0,0,0,0);
185
186			for(int i=0;i<nb_points;i++)
187	souaissa	66	{
188	souaissa	150	OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
189			//barycentre+=lst_points[i];
190			barycentre+=pt;
191	souaissa	66	}
192
193	souaissa	69	barycentre*=1./nb_points;
194	souaissa	150
195	souaissa	69	return barycentre;
196	souaissa	150
197	souaissa	66	}
198
199
200
201	souaissa	68	int VCT_FACE::get_nb_points()
202	souaissa	66	{
203	souaissa	68	return nb_points ;
204	souaissa	66	}
205
206
207	souaissa	150	OT_TENSEUR VCT_FACE::calcule_covariance(OT_TENSEUR& T)
208	souaissa	66	{
209	souaissa	150	OT_TENSEUR COVARIANCE(4);
210			OT_VECTEUR_4DD barycentre=calcule_barycentre( T);
211			OT_TENSEUR PT_CENTRE(nb_points,4),PT_CENTRE_TRANSPOSE;
212	souaissa	79
213	souaissa	66
214	souaissa	150	for( int i=0;i<nb_points;i++)
215			{
216			//OT_VECTEUR_4DD pt=lst_points[i];
217			OT_VECTEUR_4DD pt=get_nouveau_pt(T,lst_points[i]);
218			//OT_VECTEUR_4DD ptcent=lst_points[i]-barycentre;
219			OT_VECTEUR_4DD ptcent=pt-barycentre;
220			for (int j=0;j<4;j++)
221			PT_CENTRE(i,j)= ptcent[j];
222	souaissa	66
223	souaissa	71	}
224			PT_CENTRE_TRANSPOSE =PT_CENTRE.transpose();
225
226	souaissa	150
227	souaissa	71	COVARIANCE=PT_CENTRE_TRANSPOSE*PT_CENTRE;
228			double2 COEF= 1./nb_points;
229			COVARIANCE=COVARIANCE*COEF;
230
231	souaissa	79	for(int i=0;i<4;i++)
232			for(int j=0;j<4;j++)
233			{
234	souaissa	85	if(fabs(COVARIANCE(i,j).get_x())<=1e-8)
235	souaissa	79	COVARIANCE(i,j).set_x(0.);
236
237			}
238	souaissa	150
239
240			return COVARIANCE;
241	souaissa	66	}
242
243
244	souaissa	150	void VCT_FACE::calcule_axes_dinertie(OT_TENSEUR& T,OT_VECTEUR_4DD& D,OT_TENSEUR& V)
245	souaissa	66	{
246	souaissa	157
247	souaissa	69	int n=4;
248			int nrot;
249	souaissa	66
250	souaissa	150	OT_TENSEUR COV=calcule_covariance(T);
251	souaissa	71
252
253	souaissa	69	double2 zro=0.0;
254	souaissa	66
255	souaissa	71	if(COV(0,0)==zro&&COV(1,0)==zro&&COV(2,0)==zro&&COV(3,0)==zro)
256	souaissa	66	{
257	souaissa	71	COV(0,0)=1.0;
258			COV(1,0)=zro;
259			COV(2,0)=zro;
260			COV(3,0)=zro;
261	souaissa	66	}
262	souaissa	71	if(COV(0,1)==zro&&COV(1,1)==zro&&COV(2,1)==zro&&COV(3,1)==zro)
263	souaissa	66	{
264	souaissa	71	COV(0,1)=zro;
265			COV(1,1)=1.0;
266			COV(2,1)=zro;
267			COV(3,1)=zro;
268	souaissa	66	}
269	souaissa	71	if(COV(0,2)==zro&&COV(1,2)==zro&&COV(2,2)==zro&&COV(3,2)==zro)
270	souaissa	66	{
271	souaissa	71	COV(0,2)=zro;
272			COV(1,2)=zro;
273			COV(2,2)=1.0;
274			COV(3,2)=zro;
275	souaissa	66	}
276	souaissa	71	if(COV(0,3)==zro&&COV(1,3)==zro&&COV(2,3)==zro&&COV(3,3)==zro)
277	souaissa	69	{
278	souaissa	71	COV(0,3)=zro;
279			COV(1,3)=zro;
280			COV(2,3)=zro;
281			COV(3,3)=1.0;
282	souaissa	69	}
283	souaissa	66
284
285	souaissa	157	COV.get_orthogonalisation(COV,D,V,n,nrot); /*
286	souaissa	66
287	souaissa	157
288	souaissa	150	OT_VECTEUR_4DD G1;
289			OT_TENSEUR p(4),q(4),pT,qT;
290	souaissa	66
291	souaissa	150	G1=calcule_barycentre(T);
292			//OT_VECTEUR_4DD VNul(0.0001,0.0001,0.0001,0.0001);
293			double2 z=0.0001;
294			double2 cof=-1;
295			int plan=-1;
296	souaissa	66
297	souaissa	150	pT=p.transpose();
298
299			vector<OT_VECTEUR_4DD>& list_points1=get_points_controle_pour_axes();
300			vector<OT_VECTEUR_4DD>& list_vecteurs1=lst_vecteurs_axes;//get_vecteurs();
301			OT_VECTEUR_4DD Normal1(0.,0.,0.,0.),N1;
302			OT_VECTEUR_4DD axe11,axe12,axe13;
303			int cmpt=0;
304
305			OT_VECTEUR_4DD v11= list_vecteurs1[0];
306			int nb_vecteur1=list_vecteurs1.size();
307			for(unsigned int pi=0;pi<nb_vecteur1-1;pi++)
308			{
309			OT_VECTEUR_4DD v2= list_vecteurs1[pi+1];
310			if(pi==0)
311			Normal1=v11^v2;
312			if (pi>0)
313			N1=v11^v2;
314			OT_VECTEUR_4DD colineaire= Normal1^N1;
315			if(colineaire[0].get_fabs()<z&&colineaire[1].get_fabs()<z&&colineaire[2].get_fabs()<z&&colineaire[3].get_fabs()<z)
316			{
317			cmpt++;
318			}
319			}
320			vector<OT_VECTEUR_4DD> axes1;
321			if (cmpt== list_vecteurs1.size()-1)
322			{
323			plan=1;
324			axe11=list_vecteurs1[0];
325			axe12=Normal1^axe11;
326			axe13=Normal1;
327			axes1.insert(axes1.end(),axe11);
328			axes1.insert(axes1.end(),axe12);
329			axes1.insert(axes1.end(),axe13);
330			}
331
332			VCT_OUTILS VOUTIL(4);
333			OT_MATRICE_3D sys_axes1=VOUTIL.get_system_axes(G1,axes1,list_points1) ;
334			for(int i=0;i<4;i++)
335			{
336			for(int j=0;j<4;j++)
337			{
338			if(i<3&&j<3)
339			V(i,j)= sys_axes1(i,j);
340			if(i<3 &&j>2)
341			V(i,j)=0.;
342			if(i==3&&j==3) V(i,j)=1.;
343			}
344	souaissa	157	} */
345	souaissa	150	}
346
347
348
349			OT_TENSEUR VCT_FACE:: calcule_tenseur_inertie_au_barycentre(OT_TENSEUR& T) //repere globale
350	souaissa	66	{
351	souaissa	71	OT_TENSEUR INERTIE(4);
352	souaissa	150	OT_TENSEUR cov=calcule_covariance(T);
353	souaissa	69	double2 nb_pts=nb_points;
354			double2 MOINS_UN=-1.0;
355	souaissa	71	INERTIE=cov*MOINS_UN;
356			INERTIE=INERTIE*nb_pts;
357	souaissa	66
358	souaissa	71	INERTIE(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
359			INERTIE(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
360			INERTIE(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
361			INERTIE(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
362	souaissa	79	for(int i=0;i<4;i++)
363			for(int j=0;j<4;j++)
364			{
365			if(fabs(INERTIE(i,j).get_x())<=1e-12) INERTIE(i,j).set_x(0.);
366	souaissa	71
367	souaissa	79	}
368	souaissa	150
369
370	souaissa	71	return INERTIE;
371	souaissa	150
372
373	souaissa	66	}
374
375
376	souaissa	150	OT_TENSEUR VCT_FACE:: calcule_tenseur_inertie_au_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& POINT)
377	souaissa	66	{
378
379	souaissa	150	OT_VECTEUR_4DD BARYCENTRE=calcule_barycentre(T);
380			OT_TENSEUR TENS=calcule_tenseur_inertie_au_barycentre(T);
381	souaissa	66
382	souaissa	69	OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;
383	souaissa	66
384	souaissa	69	double2 a=ABC[0] ;
385			double2 b=ABC[1] ;
386			double2 c=ABC[2] ;
387			double2 d=ABC[3] ;
388	souaissa	66
389	souaissa	69	TENS(0,0)=TENS(0,0)+bb+cc+d*d;
390			TENS(1,1)=TENS(1,1)+aa+cc+d*d;
391			TENS(2,2)=TENS(2,2)+aa+bb+d*d;
392	souaissa	71	TENS(3,3)=TENS(3,3)+aa+bb+c*c;
393	souaissa	66
394	souaissa	69	TENS(0,1)=TENS(0,1)-nb_pointsab;
395			TENS(0,2)=TENS(0,2)-nb_pointsac;
396			TENS(0,3)=TENS(0,3)-nb_pointsad;
397	souaissa	66
398	souaissa	69	TENS(1,0)=TENS(1,0)-nb_pointsba;
399			TENS(1,2)=TENS(1,2)-nb_pointsbc;
400			TENS(1,3)=TENS(1,3)-nb_pointsbd;
401	souaissa	66
402	souaissa	69	TENS(2,0)=TENS(2,0)-nb_pointsca;
403			TENS(2,1)=TENS(2,1)-nb_pointscb;
404			TENS(2,3)=TENS(2,3)-nb_pointscd;
405	souaissa	66
406	souaissa	69	TENS(3,0)=TENS(3,0)-nb_pointsda;
407			TENS(3,1)=TENS(3,1)-nb_pointsdb;
408			TENS(3,2)=TENS(3,2)-nb_pointsdc;
409
410			return TENS;
411	souaissa	66	}
412
413
414	souaissa	150	OT_TENSEUR VCT_FACE::calcule_tenseur_inertie_base_locale(OT_TENSEUR& T)
415	souaissa	66	{
416	souaissa	71
417	souaissa	69	OT_VECTEUR_4DD v1,v2,v3,v4;
418			OT_TENSEUR INERTIE;
419			OT_TENSEUR I_GLOBALE;
420	souaissa	71	OT_VECTEUR_4DD D;
421			OT_TENSEUR V(4);
422	souaissa	150	this->calcule_axes_dinertie(T,D,V);
423			INERTIE=this->calcule_tenseur_inertie_au_barycentre(T);
424	souaissa	66
425	souaissa	69	I_GLOBALE= INERTIE;
426	souaissa	66
427	souaissa	69	OT_TENSEUR P(4);
428			OT_TENSEUR I=I_GLOBALE;
429	souaissa	150	OT_TENSEUR P_TRSPOSE(4);
430	souaissa	69	OT_TENSEUR I_LOCALE;
431	souaissa	71	P=V;
432	souaissa	150	//---------------
433			OT_MATRICE_3D p1;
434			for(int t=0;t<3;t++)
435			for(int tt=0;tt<3;tt++)
436			p1(t,tt)=P(t,tt).get_x();
437			p1=p1.inverse();
438			for(int t=0;t<3;t++)
439			for(int tt=0;tt<3;tt++)
440			P_TRSPOSE(t,tt)=p1(t,tt);
441			P_TRSPOSE(3,3)=1.;
442			//---------------
443			//P_TRSPOSE=P.transpose();
444	souaissa	66
445	souaissa	69	I_LOCALE=P_TRSPOSE*I;
446			I_LOCALE=I_LOCALE*P;
447	souaissa	79	for(int i=0;i<4;i++)
448			for(int j=0;j<4;j++)
449			{
450			if(fabs(I_LOCALE(i,j).get_x())<=1e-12) I_LOCALE(i,j).set_x(0.);
451			}
452	souaissa	71	return I_LOCALE;
453	souaissa	150	return INERTIE;
454	souaissa	71
455	souaissa	66	}
456
457
458
459
460	souaissa	150
461			OT_TENSEUR VCT_FACE::calcule_tenseur_inertie_base_entite(OT_TENSEUR& T,VCT& vct_f)
462	souaissa	66	{
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	66
477	souaissa	71	P=V;
478			Q=W;
479	souaissa	66
480	souaissa	150	IW=vct_f.calcule_tenseur_inertie_base_locale(T);
481	souaissa	69
482			pt=P.transpose();
483			qt=Q.transpose();
484	souaissa	66	R=pt*Q;
485			Rt=qt*P;
486			IV=R*IW;
487			IV=IV*Rt;
488
489
490	souaissa	150	G1=this->calcule_barycentre(T);
491			G2=vct_f.calcule_barycentre(T);
492	souaissa	66
493	souaissa	69	G1G2=G2-G1;
494	souaissa	66
495	souaissa	69	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	souaissa	66	{
500	souaissa	69	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	souaissa	66
509	souaissa	69	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	souaissa	66
514	souaissa	69	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	66
519	souaissa	69	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	souaissa	66
524	souaissa	69	IV=IV+MASSE_CONCENTRE;
525	souaissa	66
526	souaissa	79
527			return IV;
528
529	souaissa	69	}
530	souaissa	66
531
532	souaissa	150
533
534
535
536	souaissa	69	ostream& operator <<(ostream& os, VCT_FACE& vct_f)
537	souaissa	66	{
538	souaissa	71
539	souaissa	66	vct_f.enregistrer(os) ;
540			return os;
541			}
542
543	souaissa	150
544			OT_VECTEUR_4DD VCT_FACE::get_nouveau_pt(OT_TENSEUR& T,OT_VECTEUR_4DD& v)
545			{
546			OT_VECTEUR_4DD temp;
547
548			for( int i=0;i<4;i++)
549			{ temp[i]=0.;
550			for(int j=0;j<4;j++)
551			temp[i]= temp[i]+T(i,j)*v[j];
552			}
553			return temp;
554
555			}
556	souaissa	66	//==========================================================================
557			//Visualisation
558			//==========================================================================
559
560
561
562
563			void VCT_FACE::enregistrer(std::ostream& ost)
564			{
565	souaissa	150	OT_TENSEUR TT(4);
566			for(int i=0;i<4;i++)
567			TT(i,i)=1.0;
568			T=TT;
569	souaissa	71	ost<<"========================================"<<endl;
570			ost<<"POINTS_DE_CONTROLS: "<<endl;
571			ost<<"========================================"<<endl;
572	souaissa	66
573	souaissa	71	for (unsigned int i=0;i< lst_points.size();i++)
574			{
575	francois	222	OT_VECTEUR_4DD v= lst_points[i]; // Rmq: la precision n�est pas affocher, il faut rajouter l'affichage
576	souaissa	71	ost<< v<<endl; // de la precision dans la classe doubleprecision
577			}
578	souaissa	150	ost<<endl;
579	francois	222
580	souaissa	150	ost<<"========================================"<<endl;
581			ost<<"POINTS_DE_CONTROLS_AXES: "<<endl;
582			ost<<"========================================"<<endl;
583	souaissa	66
584	souaissa	150	for (unsigned int i=0;i< lst_points_axes.size();i++)
585			{
586	francois	222	OT_VECTEUR_4DD v= lst_points_axes[i]; // Rmq: la precision n�est pas affocher, il faut rajouter l'affichage
587	souaissa	150	ost<< v<<endl; // de la precision dans la classe doubleprecision
588			}
589			ost<<endl;
590	souaissa	71	ost<<"VECTORISATION: "<<endl;
591			ost<<"========================================"<<endl;
592			for (unsigned int i=0;i< lst_vecteurs.size();i++)
593			{
594			OT_VECTEUR_4DD v= lst_vecteurs[i];
595			ost<< v<<endl;
596	souaissa	79	}
597	souaissa	150	ost<<endl;
598			ost<<endl<<endl;
599	francois	222	ost<<"BARYCENTRE: "<<endl;
600	souaissa	71	ost<<"========================================"<<endl;
601	souaissa	150	OT_VECTEUR_4DD BARY= calcule_barycentre(T);
602	souaissa	71	ost<<BARY<<endl;
603	francois	222	ost<<endl<<endl;
604	souaissa	71	ost<<"TENSEUR_METRIQUE: "<<endl;
605			ost<<"========================================"<<endl;
606			OT_TENSEUR TNS_MT= calcule_tenseur_metrique() ;
607			ost<< TNS_MT<<endl;
608			ost<<endl<<endl;
609	francois	222	/*
610	souaissa	71	ost<<"TENSEUR_DE_COVARIANCE: "<<endl;
611			ost<<"========================================"<<endl;
612			OT_TENSEUR TNS_CV= calcule_covariance();
613			ost<< TNS_CV<<endl;
614	souaissa	150	ost<<endl<<endl; */
615	souaissa	71	ost<<"AXES_D'INERTIE: "<<endl;
616			ost<<"========================================"<<endl;
617			OT_VECTEUR_4DD D;
618			OT_TENSEUR V(4);
619	souaissa	150	calcule_axes_dinertie(T,D,V);
620
621			ost<<V<<endl;
622			ost<<endl<<endl; /*
623	francois	222	ost<<"INERTIE_CALCUL�E_AU_BARYCENTRE: "<<endl;
624	souaissa	71	ost<<"========================================"<<endl;
625			OT_TENSEUR I_BARY=calcule_tenseur_inertie_au_barycentre();
626			ost<<I_BARY<<endl;
627	souaissa	79	/*
628	souaissa	71	ost<<endl<<endl;
629	francois	222	ost<<"INERTIE_TRANSPORT�E_EN_UN_POINT: "<<endl;
630	souaissa	71	ost<<"========================================"<<endl;
631			OT_TENSEUR I_TRANSP=calcule_tenseur_inertie_au_pt(BARY);
632			ost<<I_TRANSP<<endl;
633	souaissa	150	ost<<endl<<endl; / */
634	souaissa	66
635	francois	222	ost<<"INERTIE_CALCUL�_DANS_LA_BASE_LOCALE: "<<endl;
636	souaissa	71	ost<<"========================================"<<endl;
637	souaissa	150	OT_TENSEUR I_BASE=calcule_tenseur_inertie_base_locale(T);
638	souaissa	71	ost<<I_BASE<<endl;
639	souaissa	66
640	souaissa	150
641	souaissa	66	}
642
643
644
645
646
647
648