geometrie/src/vct_courbe.cpp

#include "gestionversion.h"

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


#pragma hdrstop

#include "vct_courbe.h"
#include "mg_courbe.h"
#include "mg_arete.h"
#include <iomanip.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_COURBE::VCT_COURBE(MG_COURBE* elemgeo):VCT_ELEMENT_GEOMETRIQUE(elemgeo)
{
int indx_premier_ptctr;
OT_VECTEUR_4DD V1,V2,V;
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);
double2 ZERO=0.;
OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);


int nb_u=nurbs_params.get(3);

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


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

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

           V2[0]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1));
           V2[1]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 1);
           V2[2]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 2);
           V2[3]=nurbs_params.get(indx_premier_ptctr+  4 * (s + 1) + 3);

           lst_points.insert(lst_points.end(),V2);

           V=V2-V1;

           double2 norm_au_carre=(V[0]*V[0])+(V[1]*V[1])+(V[2]*V[2])+(V[3]*V[3]);
           double2 norm=norm_au_carre^0.5;

           if(V==VCT_NUL)
           lst_vecteurs.insert(lst_vecteurs.end(),VCT_NUL);

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


        }
}


VCT_COURBE::VCT_COURBE(VCT_COURBE& mdd):VCT_ELEMENT_GEOMETRIQUE(mdd.elem_geo)
{
lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;
}


VCT_COURBE::~VCT_COURBE()
{
}


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


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

OT_VECTEUR_4DD VCT_COURBE::calcule_barycentre()
{

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

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

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


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


OT_TENSEUR  VCT_COURBE::calcule_covariance(void)
{

 std::vector<OT_VECTEUR_4DD> Points_Centre_au_barycentre(nb_points);
 OT_VECTEUR_4DD barycentre=calcule_barycentre();

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

 OT_TENSEUR COVARIANCE(4);

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


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

 double2 zro=0.0;

 if(cov(0,0)==zro&&cov(1,0)==zro&&cov(2,0)==zro&&cov(3,0)==zro)
   {
   cov(0,0)=1.0;
   cov(1,0)=zro;
   cov(2,0)=zro;
   cov(3,0)=zro;
   }
 if(cov(0,1)==zro&&cov(1,1)==zro&&cov(2,1)==zro&&cov(3,1)==zro)
   {
   cov(0,1)=zro;
   cov(1,1)=1.0;
   cov(2,1)=zro;
   cov(3,1)=zro;
   }
  if(cov(0,2)==zro&&cov(1,2)==zro&&cov(2,2)==zro&&cov(3,2)==zro)
   {
   cov(0,2)=zro;
   cov(1,2)=zro;
   cov(2,2)=1.0;
   cov(3,2)=zro;
   }
  if(cov(0,3)==zro&&cov(1,3)==zro&&cov(2,3)==zro&&cov(3,3)==zro)
   {
   cov(0,3)=zro;
   cov(1,3)=zro;
   cov(2,3)=zro;
   cov(3,3)=1.0;
   }


cov.get_orthogonalisation(cov,d,v,n,nrot);

}


OT_TENSEUR VCT_COURBE:: calcule_tenseur_inertie_au_barycentre()  //repere globale
{
 OT_TENSEUR Inertie(4);
 OT_TENSEUR cov=calcule_covariance();
 double2 nb_pts=nb_points;
 double2 MOINS_UN=-1.0;
 Inertie=cov*MOINS_UN;
 Inertie=Inertie*nb_pts;

 Inertie(0,0)=(cov(1,1)+cov(2,2)+cov(3,3))*nb_pts;
 Inertie(1,1)=(cov(0,0)+cov(2,2)+cov(3,3))*nb_pts;
 Inertie(2,2)=(cov(1,1)+cov(0,0)+cov(3,3))*nb_pts;
 Inertie(3,3)=(cov(1,1)+cov(2,2)+cov(0,0))*nb_pts;
 return Inertie;
}


OT_TENSEUR VCT_COURBE:: calcule_tenseur_inertie_au_pt(OT_VECTEUR_4DD& POINT)
{

OT_VECTEUR_4DD BARYCENTRE=this->calcule_barycentre();
OT_TENSEUR TENS=this->calcule_tenseur_inertie_au_barycentre();

OT_VECTEUR_4DD ABC=POINT-BARYCENTRE;

double2 a=ABC[0] ;
double2 b=ABC[1] ;
double2 c=ABC[2] ;
double2 d=ABC[3] ;

TENS(0,0)=TENS(0,0)+b*b+c*c+d*d;
TENS(1,1)=TENS(1,1)+a*a+c*c+d*d;
TENS(2,2)=TENS(2,2)+a*a+b*b+d*d;
TENS(3,3)=TENS(3,3)+a*a+b*b+d*d;

TENS(0,1)=TENS(0,1)-nb_points*a*b;
TENS(0,2)=TENS(0,2)-nb_points*a*c;
TENS(0,3)=TENS(0,3)-nb_points*a*d;

TENS(1,0)=TENS(1,0)-nb_points*b*a;
TENS(1,2)=TENS(1,2)-nb_points*b*c;
TENS(1,3)=TENS(1,3)-nb_points*b*d;

TENS(2,0)=TENS(2,0)-nb_points*c*a;
TENS(2,1)=TENS(2,1)-nb_points*c*b;
TENS(2,3)=TENS(2,3)-nb_points*c*d;

TENS(3,0)=TENS(3,0)-nb_points*d*a;
TENS(3,1)=TENS(3,1)-nb_points*d*b;
TENS(3,2)=TENS(3,2)-nb_points*d*c;

return TENS;
}


OT_TENSEUR VCT_COURBE::calcule_tenseur_inertie_base_locale()
{
// inertie calcul� dans le repere globale
 OT_VECTEUR_4DD v1,v2,v3,v4;
 OT_TENSEUR INERTIE;
 OT_TENSEUR I_GLOBALE;

 this->calcule_axes_dinertie(v1,v2,v3,v4);
 INERTIE=this->calcule_tenseur_inertie_au_barycentre();

 I_GLOBALE= INERTIE;

 OT_TENSEUR P(4);
 OT_TENSEUR I=I_GLOBALE;
 OT_TENSEUR P_TRSPOSE;
 OT_TENSEUR I_LOCALE;

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

 P_TRSPOSE=P.transpose();

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

}


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

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


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

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

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

IV=R*IW;
IV=IV*Rt;


G1=this->calcule_barycentre();
G2=vct_f.calcule_barycentre();

G1G2=G2-G1;

OT_VECTEUR_4DD G1G2_LOC;  //Calcul du vecteur G1G2 dans la base locale;
  for(int i=0;i<4;i++)
  {
    for(int j=0;j<4;j++)
    {
    G1G2_LOC[i]= G1G2_LOC[i]+pt(i,j)*G1G2[j];
    }
  }


OT_TENSEUR MASSE_CONCENTRE(4);

MASSE_CONCENTRE(0,0)=nb_points*(G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(0,1)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[1]);
MASSE_CONCENTRE(0,2)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);
MASSE_CONCENTRE(0,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);

MASSE_CONCENTRE(1,0)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[0]);
MASSE_CONCENTRE(1,1)=   nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[2]*G1G2_LOC[2]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(1,2)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[2]);
MASSE_CONCENTRE(1,3)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);

MASSE_CONCENTRE(2,0)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[0]);
MASSE_CONCENTRE(2,1)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[1]);
MASSE_CONCENTRE(2,2)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[3]*G1G2_LOC[3]);
MASSE_CONCENTRE(2,3)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[2]);

MASSE_CONCENTRE(3,0)=-1*nb_points*(G1G2_LOC[0]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,1)=-1*nb_points*(G1G2_LOC[1]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,2)=-1*nb_points*(G1G2_LOC[2]*G1G2_LOC[3]);
MASSE_CONCENTRE(3,3)=nb_points*(G1G2_LOC[0]*G1G2_LOC[0]+G1G2_LOC[1]*G1G2_LOC[1]+G1G2_LOC[2]*G1G2_LOC[2]);

IV=IV+MASSE_CONCENTRE;
}


void VCT_COURBE::enregistrer(std::ostream& o)
 {
   /*
  int compt=0;
   for(int i=indx_premier_ptctr;i<nurbs_params.get_nb();i++)
    {
      o<<setw(15)<<nurbs_params.get(i);
       compt++;
       if(compt==4) {
       o<<endl;
       compt=0;}
    }
  */
 }

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