geometrie/src/vct_arete.cpp

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


#pragma hdrstop
#include "vct_arete.h"
#include "vct_courbe.h"
#include "mg_courbe.h"
#include "mg_arete.h"
#include <iomanip.h>
#include "ot_mathematique.h"
//---------------------------------------------------------------------------

#pragma package(smart_init)


VCT_ARETE::VCT_ARETE(MG_ARETE* arete):VCT_ELEMENT_TOPOLOGIQUE(arete)
{
double2 ZERO=0.;
int indx_premier_ptctr;
OT_VECTEUR_4DD V1,V2,V;

TPL_LISTE_ENTITE<double> nurbs_params;
arete->get_param_NURBS(indx_premier_ptctr,nurbs_params);
nb_points=1/4.*(nurbs_params.get_nb()-indx_premier_ptctr);

OT_VECTEUR_4DD VCT_NUL(0.,0.,0.,0.);

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


    for(int s=0;s<nb_points-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_ARETE::VCT_ARETE(VCT_ARETE& mdd):VCT_ELEMENT_TOPOLOGIQUE(mdd.elem_topo)
{

lst_vecteurs=mdd.lst_vecteurs;
lst_points=mdd.lst_points;
nb_points=mdd.nb_points ;

}


VCT_ARETE::~ VCT_ARETE()
{

}

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


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

OT_VECTEUR_4DD VCT_ARETE::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_ARETE::get_nb_points()
 {
 return nb_points  ;
 }


OT_TENSEUR  VCT_ARETE::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_ARETE::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_ARETE:: 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_ARETE:: 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_ARETE::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_ARETE::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_ARETE::enregistrer(ostream& o)
{
// o<<"ARETE__"<<((MG_ARETE*)elem_topo)->get_id()<<"  :     ";
// o<<"__CTRS_PTS"<<endl;
//  int cmpt=0;
   /*
   for(int i=indx_premier_ptctr;i<nurbs_params.get_nb();i++)
    {
      //if(cmpt!=3)
      o<<setprecision(9)<<nurbs_params.get(i)<<"  ;  ";
      //cmpt++;
      //if(cmpt==4) cmpt=0;

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