outil/src/ot_tenseur.cpp

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

#pragma hdrstop

#include "ot_tenseur.h"

#include <iomanip>
//#include <conio.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
//---------------------------------------------------------------------------

#pragma package(smart_init)


OT_TENSEUR::OT_TENSEUR()
{
 Dim1=0;
 Dim2=0;
 data=0;
}

OT_TENSEUR::OT_TENSEUR(int n)
{
 Dim1=n;
 Dim2=n;
  double2 zero(0.);
 data =new double2[n*n];
     for(int i=0;i<Dim1*Dim2;i++) data[i]=zero;
 }

OT_TENSEUR::OT_TENSEUR(int n,int m)
{
 Dim1=n;
 Dim2=m;
 double2 zero(0.);
 data =new double2[n*m];
      for(int i=0;i<Dim1*Dim2;i++) data[i]=zero;
}


OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
{

 Dim1=v.size();
 Dim2=Dim1;
 double2 zero=0.;
 double2 prec=1e-6;
 prec.set_dx(1e-8);
 double2 c=0.5;
 OT_VECTEUR_4DD VCT_NUL(1e-7,1e-7,1e-7,1e-7), w;
 int Dim= Dim1*Dim2  ;
 data=new double2[Dim];
    for(int i=0;i<Dim1;i++)
     {
      OT_VECTEUR_4DD z1=v[i];
     for (int j=0;j<Dim2;j++)
       {
       OT_VECTEUR_4DD z2=v[j];
       double2 prscal=z1*z2;
       w=z1^z2;
       double2 wnorm=w.norme();
         if(wnorm<prec)
          {
            if(prscal<zero && prscal.get_fabs()>c) prscal=  -1.;
            if(prscal>zero && prscal.get_fabs()>c) prscal=   1.;
            if(prscal<zero && prscal.get_fabs()<c) prscal=   0.;
            if(prscal>zero && prscal.get_fabs()<c) prscal=   0.;
          }
       double2 val= prscal.get_fabs();
       if(val.get_x()<=val.get_dx())  {
       prscal=0.;
       prscal.set_dx(val.get_dx());
       }
       data[i*Dim2+j]=prscal;
       }
     }


}

OT_TENSEUR::OT_TENSEUR(vector<OT_VECTEUR_4DD>&v1,vector<OT_VECTEUR_4DD>&v2)
{

 Dim1=v1.size();
 Dim2=v2.size();
 int Dim= Dim1*Dim2  ;
 data=new double2[Dim];
    for(int i=0;i<Dim1;i++)
     {
      OT_VECTEUR_4DD z1=v1[i];
     for (int j=0;j<Dim2;j++)
       {
       OT_VECTEUR_4DD z2=v2[j];
       double2 prscal=z1*z2;
       double2 val= prscal.get_fabs();
       if(val.get_x()<=val.get_dx())  {
       prscal=0.;
       prscal.set_dx(val.get_dx());
       }
       data[i*Dim2+j]=prscal;
       }
     }


}


OT_TENSEUR::~OT_TENSEUR()
{
Dim1=0;
Dim2=0;

if (data!=NULL) delete [] data;
}


OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
{
Dim1=ot.Dim1;
Dim2=ot.Dim2;
data=new double2[Dim1*Dim2];

    for(int i=0;i<Dim1;i++)
     {
     for (int j=0;j<Dim2;j++)
       {
       data[i*Dim2+j]=ot(i,j);
       }
     }
}


OT_TENSEUR& OT_TENSEUR:: operator =(const OT_TENSEUR& ot)
{

 if(this==&ot) return *this ;

 if(data)  delete [] data;

  Dim1=ot.Dim1;
  Dim2=ot.Dim2;
  data=new double2[Dim1*Dim2];

   for(int i=0;i<Dim1;i++) {
     for (int j=0;j<Dim2;j++) {
        data[i*Dim2+j]=ot.data[i*Dim2+j];
      }
     }


 return *this;
}

int operator==( OT_TENSEUR& t1, OT_TENSEUR&t2)
{
if(t1.Get_NbRows()!=t2.Get_NbRows()||t1.Get_NbCols()!=t2.Get_NbCols()) return 0;
else {
  for(int i=0;i<t1.Get_NbRows();i++)
    {
     for(int j=0;j<t1.Get_NbCols();j++)
       {
       double2 val1=t1(i,j);
       double2 val2=t2(i,j);
      if(val1==val2) continue;
      else return 0;
       }

    }
     }
return 1;
}

void OT_TENSEUR::identite(void)
{
double2 un(1.0);
for(int i=0;i<Dim1;i++) 
        data[i*Dim2+i]=un;
}


double2& OT_TENSEUR:: operator()(int i,int j)
{
assert(i<Dim1&&j<Dim2);
double2 prciseion=0.0;
if(data[i*Dim2+j].get_fabs()<=prciseion)
{
data[i*Dim2+j]=0.;

}
return data[i*Dim2+j];
}


double2& OT_TENSEUR:: operator()(int i,int j)const
{
assert(i<Dim1&&j<Dim2);
double2 prciseion=0.0;
if(data[i*Dim2+j].get_fabs()<=prciseion)
{
data[i*Dim2+j]=0.;

}
return data[i*Dim2+j];
}


int OT_TENSEUR:: Get_NbRows()const
{
return Dim1;
}


int OT_TENSEUR:: Get_NbCols()const
{
return Dim2;
}


int OT_TENSEUR:: est_til_equivalent( OT_TENSEUR& tns)
 {

if (Dim1!=tns.Get_NbRows()||Dim2!=tns.Get_NbCols()) return 0;
vector<int> check;

 for(int i=0;i<tns.Get_NbRows();i++)
  {

  vector<double2> lig1;
  for(int ix1=0;ix1<Dim2;ix1++) lig1.insert(lig1.end(),data[i*Dim2+ix1]);

  for(int j=0;j<Dim2;j++)
    {
    int p=-1;
    if(!existe(check,j,p))
     {
    vector<int>trouve;
    vector<double2>valeur;
    vector<double2> lig2;
    for(int ix2=0;ix2<Dim2;ix2++) lig2.insert(lig2.end(),tns(j,ix2));
          for(int r=0;r<Dim2;r++)
           {
            double2 val=lig1[r];
            val=val.get_fabs();

            int indx1=-1;

            for(int s=0;s<Dim2;s++)
             {  int pos1;
                double2 vall=lig2[s];
                vall=vall.get_fabs();

                if (val==vall&&!existe(trouve,s,indx1))
                  {
                  trouve.insert(trouve.end(),s);
                  valeur.insert(valeur.end(),val);
                  break;
                  }

             }

           }
    int t_trouve=trouve.size();
    if (t_trouve==Dim2)
     {
     check.insert(check.end(),j);
     break;

     }
      
    }
    }
continue;
 }
int t_chek=check.size();
if (t_chek==Dim1)return 1;
else return 0;

 }


template <class T> int OT_TENSEUR:: existe(const vector<T>& vct,T val,int& idx)
{

  if(val<0)   val=-val;
  for(unsigned int i=0;i<vct.size();i++)
    {
    double x=vct[i];
    if(x<0)   val=-val;
    if(x==val){
       idx=i; return 1;  }
    }
 idx=-1;
return 0;
}


 OT_TENSEUR operator+( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
 {
 assert (ot1.Get_NbRows()==ot2.Get_NbRows()&&ot1.Get_NbCols()==ot2.Get_NbCols());
 OT_TENSEUR res(ot1.Get_NbRows(),ot1.Get_NbCols());
    for(int i=0;i<ot1.Get_NbRows();i++)
       {
        for(int j=0;j<ot2.Get_NbCols();j++)
           {
           res(i,j)=ot1(i,j)+ot2(i,j);
           }
        }

 return res;
 }


vector<double2> operator*( OT_TENSEUR& ot1, vector<double2>& v)
 {
 assert (ot1.Get_NbCols()==v.size());
 vector<double2> res(ot1.Get_NbCols());
    for(int i=0;i<ot1.Get_NbRows();i++)
       {
       res[i]=double2(0.0);
        for(int j=0;j<ot1.Get_NbCols();j++)
           {
           res[i]=res[i]+ot1(i,j)+v[j];
           }
        }

 return res;


 }

OT_TENSEUR operator*( double2& coef, OT_TENSEUR& ot)
{
OT_TENSEUR T=ot;
for(int i=0;i<ot.Get_NbRows();i++)
 for(int j=0;j<ot.Get_NbCols();j++)

  T(i,j)=coef*T(i,j);

return T;

}

OT_TENSEUR operator*(OT_TENSEUR& ot , double2& coef)
{
return coef*ot;
}

OT_TENSEUR OT_TENSEUR::inverse_homogene(void)
{
OT_TENSEUR TMP(4,4);
if (Dim1!=4) return TMP;
if (Dim2!=4) return TMP;
double2 zero(0.);
double2 un(1.);
OT_VECTEUR_4DD x1((*this)(0,0),(*this)(1,0),(*this)(2,0),zero);
OT_VECTEUR_4DD x2((*this)(0,1),(*this)(1,1),(*this)(2,1),zero);
OT_VECTEUR_4DD x3((*this)(0,2),(*this)(1,2),(*this)(2,2),zero);
OT_VECTEUR_4DD x4((*this)(0,3)/(*this)(3,3),(*this)(1,3)/(*this)(3,3),(*this)(2,3)/(*this)(3,3),un);
TMP(0,0)=(*this)(0,0);TMP(1,0)=(*this)(0,1);TMP(2,0)=(*this)(0,2);TMP(3,0)=zero;
TMP(0,1)=(*this)(1,0);TMP(1,1)=(*this)(1,1);TMP(2,1)=(*this)(1,2);TMP(3,1)=zero;
TMP(0,2)=(*this)(2,0);TMP(1,2)=(*this)(2,1);TMP(2,2)=(*this)(2,2);TMP(3,2)=zero;
TMP(0,3)=zero-x4*x1;
TMP(1,3)=zero-x4*x2;
TMP(2,3)=zero-x4*x3;
TMP(3,3)=un;
return TMP;
}


OT_TENSEUR OT_TENSEUR:: transpose()
{
OT_TENSEUR TMP(Dim2,Dim1);
 for(int i=0;i<Dim1;i++) {
    for(int j=0;j<Dim2;j++)
    {
TMP(j,i)=(*this)(i,j);  // ot(i,j);
    }
  }
return TMP;
}


OT_TENSEUR operator*( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
{
int dim1=ot1.Get_NbRows();
int dim3= ot1.Get_NbCols();
int dim2=ot2.Get_NbCols();
OT_TENSEUR TMP (dim1,dim2);

   for(int i=0;i<dim1;i++)
     {
     for(int j=0;j<dim2;j++)
       {
          for(int k=0;k<dim3;k++)
           TMP(i,j)=TMP(i,j)+ot1(i,k)*ot2(k,j);
       }
     }
return TMP;
}

OT_VECTEUR_4DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_4DD& ot2)
{
OT_VECTEUR_4DD TMP(0.,0.,0.,0.);
if (ot1.Get_NbRows()!=4) return TMP;
if (ot1.Get_NbCols()!=4) return TMP;
for(int i=0;i<4;i++)
     {
     TMP[i]=TMP[i]+ot1(i,0)*ot2[0]+ot1(i,1)*ot2[1]+ot1(i,2)*ot2[2]+ot1(i,3)*ot2[3];
     }
return TMP;
}

OT_VECTEUR_3DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_3DD& ot2)
{
OT_VECTEUR_3DD TMP(0.,0.,0.);
if (ot1.Get_NbRows()!=3) return TMP;
if (ot1.Get_NbCols()!=3) return TMP;
for(int i=0;i<3;i++)
     {
     TMP[i]=TMP[i]+ot1(i,0)*ot2[0]+ot1(i,1)*ot2[1]+ot1(i,2)*ot2[2];
     }
return TMP;
}

std::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns){
os<<"------------------------------------------------------------------------------------"<<endl;
os<<"#  TENSEUR:   DIM1=  "<< tns.Get_NbRows()<<"    , DIM2:  "<<tns.Get_NbCols()<<"  #"<<endl;
os<<"------------------------------------------------------------------------------------"<<endl;
  for(int i=0;i<tns.Get_NbRows();i++) {
     for(int j=0;j<tns.Get_NbCols();j++) {
         os<<setw(18)<<tns(i,j);
                       }
            os<<endl;
            }
os<<"------------------------------------------------------------------------------------"<<endl;            
            return os;
     }


void OT_TENSEUR::get_orthogonalisation(OT_TENSEUR& a,OT_VECTEUR_4DD& d,OT_TENSEUR& v,int n,int&nrot)
{

        static const int ITMAX=50;
        double2 zro=0.0;
        int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
        double2 tresh,theta,tau,t,sm,s,h,g,c;
        double2 *b=new double2[n];                 // OT_VECTEUR_3D b,z;
        double2 *z=new double2[n];

        for(int ip=0;ip<n;ip++)
                {
                        for(int iq=0;iq<n;iq++)
                              //v[ip*n+iq]=0.0;
                              v(ip,iq)=0.0;
                        v(ip,ip)=1.0;
                        b[(unsigned int)ip]=d[(unsigned int)ip]=a(ip,ip);
                        z[(unsigned int)ip]=0.0;
                }
        nrot=0;
        for(int i=1;i<=ITMAX;i++)
                {
                sm=0.0;
                for(int ip=0;ip<n1;ip++)
                        for(int iq=ip+1;iq<n;iq++)
                               sm =sm+(a(ip,iq)).get_fabs();
                if(sm==zro)
                        return;
                if(i<4)
                        tresh=0.2*sm*n2;
                else
                        tresh=0.0;
                for(int ip=0;ip<n1;ip++)
                        for(int iq=ip+1;iq<n;iq++)
                                {
                                g=100*(a(ip,iq)).get_fabs();
                                if(i>4&&((d[(unsigned int)ip]+g).get_fabs()==d[(unsigned int)ip].get_fabs())&&((d[(unsigned int)iq])+g).get_fabs()==d[(unsigned int)iq].get_fabs())
                                        a(ip,iq)=0.0;
                                else if((a(ip,iq)).get_fabs()>tresh)
                                        {
                                                h=d[(unsigned int)iq]-d[(unsigned int)ip];
                                                if(((h).get_fabs()+g)==(h).get_fabs())
                                                        t=a(ip,iq)/h;
                                                else
                                                        {
                                                                theta=0.5*h/a(ip,iq);
                                                                t=1.0/((theta).get_fabs()+sqrt(1.0+theta*theta));
                                                                if(theta<zro)t=zro-t;
                                                        }
                                                c=1.0/sqrt(1+t*t);
                                                s=t*c;
                                                tau=s/(1.0+c);
                                                h=t*a(ip,iq);
                                                z[(unsigned int)ip]=z[(unsigned int)ip]-h;
                                                z[(unsigned int)iq]=z[(unsigned int)iq]+h;
                                                d[(unsigned int)ip]=d[(unsigned int)ip]-h;
                                                d[(unsigned int)iq]=d[(unsigned int)iq]+h;
                                                a(ip,iq)=0.0;
                                                ip1=ip-1;
                                                iq1=iq-1;
                                                for(int j=0;j<=ip1;j++)
                                                        {
                                                                g=a(j,ip);
                                                                h=a(j,iq);
                                                                a(j,ip)=g-s*(h+g*tau);
                                                                a(j,iq)=h+s*(g-h*tau);
                                                        }
                                                for(int j=ip+1;j<=iq1;j++)
                                                        {
                                                                g=a(ip,j);
                                                                h=a(j,iq);
                                                                a(ip,j)=g-s*(h+g*tau);
                                                                a(j,iq)=h+s*(g-h*tau);
                                                        }
                                                for(int j=iq+1;j<n;j++)
                                                        {
                                                                g=a(ip,j);
                                                                h=a(iq,j);
                                                                a(ip,j)=g-s*(h+g*tau);
                                                                a(iq,j)=h+s*(g-h*tau);
                                                        }
                                                for(int j=0;j<n;j++)
                                                        {
                                                                g=v(j,ip);
                                                                h=v(j,iq);
                                                                v(j,ip)=g-s*(h+g*tau);
                                                                v(j,iq)=h+s*(g-h*tau);
                                                        }
                                                ++nrot;
                                        }
                                }
                for(int ip=0;ip<n;ip++)
                        {
                                b[(unsigned int)ip]=b[(unsigned int)ip]+z[(unsigned int)ip];
                                d[(unsigned int)ip]=b[(unsigned int)ip];
                                z[(unsigned int)ip]=0.0;
                        }
                }


delete [] b;
delete [] z;

}


int  OT_TENSEUR::listes_equivalentes(std::vector<double2>& list1,std::vector<double2>& list2,std::vector<unsigned int>& indx)

{
 int compt=0;
 int pos=-1;
 int exist=0;
   for(int i=0;i<list1.size();i++)
      {
        double2 val1=list1[i];
        double prc=1e-3;
        double prcc;
           if(pos==-1)
             {
                for(int k=0;k< list2.size();k++)
                    {
                      double2 val2=list2[k];
                      prcc =val2.get_x()-val1.get_x();
                      if (fabs(val2.get_x())>val2.get_dx())
                      prcc =prcc/val2.get_x();
                           if( fabs(prcc)<prc){
                                indx.insert(indx.end(),k);
                                pos=k;
                                compt++;
                                break;
                                }
                    }
              }


  else  {
        for( int r=0;r<list2.size();r++)
           {
             double2 vval2=list2[r];
             for(int j=0;j<indx.size();j++)
                {
                  exist=0;
                  if (r==indx[j])
                     {
                      exist=1;
                      break;
                     }
                 }

              prcc =vval2.get_x()-val1.get_x();
              if(fabs(vval2.get_x())>vval2.get_dx())
              prcc =prcc/vval2.get_x();

              if(fabs(prcc)<prc&&!exist)   // vval2==val1
                 {
                  indx.insert(indx.end(),r);
                  pos=r;
                  compt++;
                  break;
                 }

            }
       }
   }

if (compt==list1.size())return 1;
else return 0;


}


Revision:	249
Committed:	Fri Jun 18 14:03:27 2010 UTC (14 years, 10 months ago) by francois
Original Path:	*magic/lib/outil/outil/src/ot_tenseur.cpp*
File size:	14195 byte(s)
Log Message:	Mise a jour toxfem au 18 juin 2010
#	User	Rev	Content
1	souaissa	58	//---------------------------------------------------------------------------
2	francois	17	#include "gestionversion.h"
3
4			#pragma hdrstop
5	souaissa	58
6	francois	17	#include "ot_tenseur.h"
7	souaissa	58
8			#include <iomanip>
9	francois	102	//#include <conio.h>
10	souaissa	58	#include <stdio.h>
11	souaissa	69	#include <assert.h>
12	souaissa	86	#include <math.h>
13	francois	17	//---------------------------------------------------------------------------
14	souaissa	58
15	francois	17	#pragma package(smart_init)
16
17
18
19	souaissa	69
20			OT_TENSEUR::OT_TENSEUR()
21	souaissa	58	{
22	souaissa	69	Dim1=0;
23			Dim2=0;
24	souaissa	70	data=0;
25	souaissa	58	}
26	souaissa	70
27	souaissa	69	OT_TENSEUR::OT_TENSEUR(int n)
28	souaissa	58	{
29	souaissa	69	Dim1=n;
30			Dim2=n;
31	souaissa	80	double2 zero(0.);
32	souaissa	70	data =new double2[n*n];
33	souaissa	80	for(int i=0;i<Dim1*Dim2;i++) data[i]=zero;
34	souaissa	70	}
35
36	souaissa	69	OT_TENSEUR::OT_TENSEUR(int n,int m)
37	francois	17	{
38	souaissa	69	Dim1=n;
39			Dim2=m;
40	souaissa	80	double2 zero(0.);
41	souaissa	70	data =new double2[n*m];
42	souaissa	80	for(int i=0;i<Dim1*Dim2;i++) data[i]=zero;
43	souaissa	58	}
44	souaissa	70
45
46
47	souaissa	69	OT_TENSEUR:: OT_TENSEUR(vector<OT_VECTEUR_4DD>& v)
48	francois	17	{
49	souaissa	70
50	souaissa	69	Dim1=v.size();
51			Dim2=Dim1;
52	souaissa	86	double2 zero=0.;
53	souaissa	93	double2 prec=1e-6;
54			prec.set_dx(1e-8);
55	souaissa	86	double2 c=0.5;
56	souaissa	93	OT_VECTEUR_4DD VCT_NUL(1e-7,1e-7,1e-7,1e-7), w;
57	souaissa	69	int Dim= Dim1*Dim2 ;
58	souaissa	70	data=new double2[Dim];
59			for(int i=0;i<Dim1;i++)
60			{
61			OT_VECTEUR_4DD z1=v[i];
62			for (int j=0;j<Dim2;j++)
63			{
64			OT_VECTEUR_4DD z2=v[j];
65			double2 prscal=z1*z2;
66	souaissa	86	w=z1^z2;
67	souaissa	93	double2 wnorm=w.norme();
68			if(wnorm<prec)
69	souaissa	86	{
70			if(prscal<zero && prscal.get_fabs()>c) prscal= -1.;
71			if(prscal>zero && prscal.get_fabs()>c) prscal= 1.;
72			if(prscal<zero && prscal.get_fabs()<c) prscal= 0.;
73			if(prscal>zero && prscal.get_fabs()<c) prscal= 0.;
74			}
75	souaissa	80	double2 val= prscal.get_fabs();
76			if(val.get_x()<=val.get_dx()) {
77			prscal=0.;
78			prscal.set_dx(val.get_dx());
79			}
80	souaissa	70	data[i*Dim2+j]=prscal;
81			}
82	souaissa	69	}
83	francois	17
84
85			}
86
87	souaissa	69	OT_TENSEUR::OT_TENSEUR(vector<OT_VECTEUR_4DD>&v1,vector<OT_VECTEUR_4DD>&v2)
88	souaissa	58	{
89	francois	17
90	souaissa	69	Dim1=v1.size();
91			Dim2=v2.size();
92			int Dim= Dim1*Dim2 ;
93	souaissa	70	data=new double2[Dim];
94			for(int i=0;i<Dim1;i++)
95			{
96			OT_VECTEUR_4DD z1=v1[i];
97			for (int j=0;j<Dim2;j++)
98			{
99			OT_VECTEUR_4DD z2=v2[j];
100			double2 prscal=z1*z2;
101	souaissa	80	double2 val= prscal.get_fabs();
102			if(val.get_x()<=val.get_dx()) {
103			prscal=0.;
104			prscal.set_dx(val.get_dx());
105			}
106	souaissa	70	data[i*Dim2+j]=prscal;
107			}
108	souaissa	69	}
109	francois	17
110	souaissa	69
111	souaissa	58	}
112	francois	17
113	souaissa	69
114			OT_TENSEUR::~OT_TENSEUR()
115	souaissa	58	{
116	souaissa	69	Dim1=0;
117			Dim2=0;
118	souaissa	70
119			if (data!=NULL) delete [] data;
120	souaissa	58	}
121	francois	17
122
123	souaissa	69	OT_TENSEUR::OT_TENSEUR(const OT_TENSEUR& ot)
124	francois	17	{
125	souaissa	69	Dim1=ot.Dim1;
126			Dim2=ot.Dim2;
127	souaissa	70	data=new double2[Dim1*Dim2];
128
129			for(int i=0;i<Dim1;i++)
130			{
131			for (int j=0;j<Dim2;j++)
132			{
133			data[i*Dim2+j]=ot(i,j);
134			}
135			}
136	francois	17	}
137
138	souaissa	58
139	souaissa	69	OT_TENSEUR& OT_TENSEUR:: operator =(const OT_TENSEUR& ot)
140	francois	17	{
141	souaissa	58
142	souaissa	69	if(this==&ot) return *this ;
143	souaissa	70
144			if(data) delete [] data;
145
146	souaissa	69	Dim1=ot.Dim1;
147			Dim2=ot.Dim2;
148	souaissa	70	data=new double2[Dim1*Dim2];
149
150			for(int i=0;i<Dim1;i++) {
151			for (int j=0;j<Dim2;j++) {
152			data[iDim2+j]=ot.data[iDim2+j];
153			}
154			}
155
156
157	souaissa	69	return *this;
158	francois	17	}
159
160	souaissa	69	int operator==( OT_TENSEUR& t1, OT_TENSEUR&t2)
161	francois	17	{
162	souaissa	69	if(t1.Get_NbRows()!=t2.Get_NbRows()\|\|t1.Get_NbCols()!=t2.Get_NbCols()) return 0;
163	souaissa	58	else {
164	souaissa	69	for(int i=0;i<t1.Get_NbRows();i++)
165	souaissa	58	{
166	souaissa	69	for(int j=0;j<t1.Get_NbCols();j++)
167	souaissa	58	{
168	souaissa	75	double2 val1=t1(i,j);
169			double2 val2=t2(i,j);
170			if(val1==val2) continue;
171			else return 0;
172	souaissa	58	}
173	francois	17
174	souaissa	58	}
175			}
176			return 1;
177			}
178	francois	17
179	francois	175	void OT_TENSEUR::identite(void)
180			{
181			double2 un(1.0);
182			for(int i=0;i<Dim1;i++)
183			data[i*Dim2+i]=un;
184			}
185	francois	17
186	souaissa	69
187	francois	175
188			double2& OT_TENSEUR:: operator()(int i,int j)
189	souaissa	58	{
190	souaissa	69	assert(i<Dim1&&j<Dim2);
191	souaissa	80	double2 prciseion=0.0;
192			if(data[i*Dim2+j].get_fabs()<=prciseion)
193			{
194			data[i*Dim2+j]=0.;
195	souaissa	69
196	souaissa	80	}
197	souaissa	70	return data[i*Dim2+j];
198	souaissa	58	}
199	francois	17
200
201	souaissa	70
202			double2& OT_TENSEUR:: operator()(int i,int j)const
203	souaissa	58	{
204	souaissa	69	assert(i<Dim1&&j<Dim2);
205	souaissa	80	double2 prciseion=0.0;
206			if(data[i*Dim2+j].get_fabs()<=prciseion)
207			{
208			data[i*Dim2+j]=0.;
209	souaissa	69
210	souaissa	80	}
211	souaissa	70	return data[i*Dim2+j];
212	souaissa	58	}
213	francois	17
214	souaissa	58
215	souaissa	70
216	souaissa	69	int OT_TENSEUR:: Get_NbRows()const
217	francois	17	{
218	souaissa	69	return Dim1;
219	francois	17	}
220
221	souaissa	69
222
223			int OT_TENSEUR:: Get_NbCols()const
224	souaissa	58	{
225	souaissa	69	return Dim2;
226	souaissa	58	}
227	francois	17
228
229	souaissa	58
230	souaissa	70	int OT_TENSEUR:: est_til_equivalent( OT_TENSEUR& tns)
231	souaissa	58	{
232	souaissa	86
233			if (Dim1!=tns.Get_NbRows()\|\|Dim2!=tns.Get_NbCols()) return 0;
234	souaissa	58	vector<int> check;
235
236	souaissa	69	for(int i=0;i<tns.Get_NbRows();i++)
237	souaissa	58	{
238
239	souaissa	69	vector<double2> lig1;
240			for(int ix1=0;ix1<Dim2;ix1++) lig1.insert(lig1.end(),data[i*Dim2+ix1]);
241
242			for(int j=0;j<Dim2;j++)
243	francois	17	{
244	souaissa	69	int p=-1;
245	souaissa	58	if(!existe(check,j,p))
246			{
247			vector<int>trouve;
248			vector<double2>valeur;
249			vector<double2> lig2;
250	souaissa	69	for(int ix2=0;ix2<Dim2;ix2++) lig2.insert(lig2.end(),tns(j,ix2));
251			for(int r=0;r<Dim2;r++)
252	souaissa	58	{
253	souaissa	69	double2 val=lig1[r];
254			val=val.get_fabs();
255	souaissa	80
256	souaissa	58	int indx1=-1;
257	francois	17
258	souaissa	69	for(int s=0;s<Dim2;s++)
259	souaissa	58	{ int pos1;
260	souaissa	69	double2 vall=lig2[s];
261			vall=vall.get_fabs();
262	souaissa	80
263			if (val==vall&&!existe(trouve,s,indx1))
264	souaissa	58	{
265			trouve.insert(trouve.end(),s);
266			valeur.insert(valeur.end(),val);
267			break;
268			}
269
270			}
271
272			}
273			int t_trouve=trouve.size();
274	souaissa	69	if (t_trouve==Dim2)
275	souaissa	58	{
276			check.insert(check.end(),j);
277			break;
278
279	francois	17	}
280	souaissa	58
281			}
282			}
283			continue;
284			}
285			int t_chek=check.size();
286	souaissa	69	if (t_chek==Dim1)return 1;
287	souaissa	58	else return 0;
288	francois	17
289	souaissa	58	}
290
291
292
293	souaissa	69	template <class T> int OT_TENSEUR:: existe(const vector<T>& vct,T val,int& idx)
294	souaissa	58	{
295	francois	17
296	souaissa	58	if(val<0) val=-val;
297	souaissa	69	for(unsigned int i=0;i<vct.size();i++)
298	souaissa	58	{
299			double x=vct[i];
300			if(x<0) val=-val;
301			if(x==val){
302			idx=i; return 1; }
303			}
304			idx=-1;
305			return 0;
306			}
307	souaissa	69
308
309			OT_TENSEUR operator+( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
310			{
311			assert (ot1.Get_NbRows()==ot2.Get_NbRows()&&ot1.Get_NbCols()==ot2.Get_NbCols());
312	souaissa	70	OT_TENSEUR res(ot1.Get_NbRows(),ot1.Get_NbCols());
313	souaissa	69	for(int i=0;i<ot1.Get_NbRows();i++)
314			{
315			for(int j=0;j<ot2.Get_NbCols();j++)
316			{
317			res(i,j)=ot1(i,j)+ot2(i,j);
318			}
319			}
320
321			return res;
322			}
323
324
325	souaissa	151	vector<double2> operator*( OT_TENSEUR& ot1, vector<double2>& v)
326			{
327			assert (ot1.Get_NbCols()==v.size());
328			vector<double2> res(ot1.Get_NbCols());
329			for(int i=0;i<ot1.Get_NbRows();i++)
330			{
331			res[i]=double2(0.0);
332			for(int j=0;j<ot1.Get_NbCols();j++)
333			{
334			res[i]=res[i]+ot1(i,j)+v[j];
335			}
336			}
337
338			return res;
339
340
341			}
342
343	souaissa	70	OT_TENSEUR operator*( double2& coef, OT_TENSEUR& ot)
344			{
345			OT_TENSEUR T=ot;
346			for(int i=0;i<ot.Get_NbRows();i++)
347			for(int j=0;j<ot.Get_NbCols();j++)
348	souaissa	69
349	souaissa	70	T(i,j)=coef*T(i,j);
350
351			return T;
352
353			}
354
355	francois	102	OT_TENSEUR operator*(OT_TENSEUR& ot , double2& coef)
356	souaissa	70	{
357			return coef*ot;
358	francois	102	}
359	francois	249
360			OT_TENSEUR OT_TENSEUR::inverse_homogene(void)
361			{
362			OT_TENSEUR TMP(4,4);
363			if (Dim1!=4) return TMP;
364			if (Dim2!=4) return TMP;
365			double2 zero(0.);
366			double2 un(1.);
367			OT_VECTEUR_4DD x1((this)(0,0),(this)(1,0),(*this)(2,0),zero);
368			OT_VECTEUR_4DD x2((this)(0,1),(this)(1,1),(*this)(2,1),zero);
369			OT_VECTEUR_4DD x3((this)(0,2),(this)(1,2),(*this)(2,2),zero);
370			OT_VECTEUR_4DD x4((this)(0,3)/(this)(3,3),(this)(1,3)/(this)(3,3),(this)(2,3)/(this)(3,3),un);
371			TMP(0,0)=(this)(0,0);TMP(1,0)=(this)(0,1);TMP(2,0)=(*this)(0,2);TMP(3,0)=zero;
372			TMP(0,1)=(this)(1,0);TMP(1,1)=(this)(1,1);TMP(2,1)=(*this)(1,2);TMP(3,1)=zero;
373			TMP(0,2)=(this)(2,0);TMP(1,2)=(this)(2,1);TMP(2,2)=(*this)(2,2);TMP(3,2)=zero;
374			TMP(0,3)=zero-x4*x1;
375			TMP(1,3)=zero-x4*x2;
376			TMP(2,3)=zero-x4*x3;
377			TMP(3,3)=un;
378			return TMP;
379			}
380
381
382
383	francois	102	OT_TENSEUR OT_TENSEUR:: transpose()
384			{
385			OT_TENSEUR TMP(Dim2,Dim1);
386			for(int i=0;i<Dim1;i++) {
387	souaissa	70	for(int j=0;j<Dim2;j++)
388			{
389			TMP(j,i)=(*this)(i,j); // ot(i,j);
390			}
391	francois	102	}
392			return TMP;
393	souaissa	70	}
394
395	francois	249
396	francois	102	OT_TENSEUR operator*( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
397			{
398			int dim1=ot1.Get_NbRows();
399			int dim3= ot1.Get_NbCols();
400			int dim2=ot2.Get_NbCols();
401			OT_TENSEUR TMP (dim1,dim2);
402
403			for(int i=0;i<dim1;i++)
404			{
405			for(int j=0;j<dim2;j++)
406			{
407			for(int k=0;k<dim3;k++)
408			TMP(i,j)=TMP(i,j)+ot1(i,k)*ot2(k,j);
409			}
410			}
411			return TMP;
412			}
413	francois	249
414			OT_VECTEUR_4DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_4DD& ot2)
415			{
416			OT_VECTEUR_4DD TMP(0.,0.,0.,0.);
417			if (ot1.Get_NbRows()!=4) return TMP;
418			if (ot1.Get_NbCols()!=4) return TMP;
419			for(int i=0;i<4;i++)
420			{
421			TMP[i]=TMP[i]+ot1(i,0)ot2[0]+ot1(i,1)ot2[1]+ot1(i,2)ot2[2]+ot1(i,3)ot2[3];
422			}
423			return TMP;
424			}
425
426			OT_VECTEUR_3DD operator*( OT_TENSEUR& ot1, OT_VECTEUR_3DD& ot2)
427			{
428			OT_VECTEUR_3DD TMP(0.,0.,0.);
429			if (ot1.Get_NbRows()!=3) return TMP;
430			if (ot1.Get_NbCols()!=3) return TMP;
431			for(int i=0;i<3;i++)
432			{
433			TMP[i]=TMP[i]+ot1(i,0)ot2[0]+ot1(i,1)ot2[1]+ot1(i,2)*ot2[2];
434			}
435			return TMP;
436			}
437	francois	102
438	souaissa	69	std::ostream& operator<<(std::ostream& os, OT_TENSEUR& tns){
439	souaissa	70	os<<"------------------------------------------------------------------------------------"<<endl;
440			os<<"# TENSEUR: DIM1= "<< tns.Get_NbRows()<<" , DIM2: "<<tns.Get_NbCols()<<" #"<<endl;
441			os<<"------------------------------------------------------------------------------------"<<endl;
442	souaissa	69	for(int i=0;i<tns.Get_NbRows();i++) {
443			for(int j=0;j<tns.Get_NbCols();j++) {
444			os<<setw(18)<<tns(i,j);
445			}
446			os<<endl;
447			}
448	souaissa	70	os<<"------------------------------------------------------------------------------------"<<endl;
449	souaissa	69	return os;
450			}
451
452
453
454			void OT_TENSEUR::get_orthogonalisation(OT_TENSEUR& a,OT_VECTEUR_4DD& d,OT_TENSEUR& v,int n,int&nrot)
455	francois	102	{
456
457			static const int ITMAX=50;
458			double2 zro=0.0;
459			int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
460			double2 tresh,theta,tau,t,sm,s,h,g,c;
461			double2 *b=new double2[n]; // OT_VECTEUR_3D b,z;
462			double2 *z=new double2[n];
463
464			for(int ip=0;ip<n;ip++)
465			{
466			for(int iq=0;iq<n;iq++)
467			//v[ip*n+iq]=0.0;
468			v(ip,iq)=0.0;
469			v(ip,ip)=1.0;
470			b[(unsigned int)ip]=d[(unsigned int)ip]=a(ip,ip);
471			z[(unsigned int)ip]=0.0;
472			}
473			nrot=0;
474			for(int i=1;i<=ITMAX;i++)
475			{
476			sm=0.0;
477			for(int ip=0;ip<n1;ip++)
478			for(int iq=ip+1;iq<n;iq++)
479			sm =sm+(a(ip,iq)).get_fabs();
480			if(sm==zro)
481			return;
482			if(i<4)
483			tresh=0.2smn2;
484			else
485			tresh=0.0;
486			for(int ip=0;ip<n1;ip++)
487			for(int iq=ip+1;iq<n;iq++)
488			{
489			g=100*(a(ip,iq)).get_fabs();
490			if(i>4&&((d[(unsigned int)ip]+g).get_fabs()==d[(unsigned int)ip].get_fabs())&&((d[(unsigned int)iq])+g).get_fabs()==d[(unsigned int)iq].get_fabs())
491			a(ip,iq)=0.0;
492			else if((a(ip,iq)).get_fabs()>tresh)
493			{
494			h=d[(unsigned int)iq]-d[(unsigned int)ip];
495			if(((h).get_fabs()+g)==(h).get_fabs())
496			t=a(ip,iq)/h;
497			else
498			{
499			theta=0.5*h/a(ip,iq);
500			t=1.0/((theta).get_fabs()+sqrt(1.0+theta*theta));
501			if(theta<zro)t=zro-t;
502			}
503			c=1.0/sqrt(1+t*t);
504			s=t*c;
505			tau=s/(1.0+c);
506			h=t*a(ip,iq);
507			z[(unsigned int)ip]=z[(unsigned int)ip]-h;
508			z[(unsigned int)iq]=z[(unsigned int)iq]+h;
509			d[(unsigned int)ip]=d[(unsigned int)ip]-h;
510			d[(unsigned int)iq]=d[(unsigned int)iq]+h;
511			a(ip,iq)=0.0;
512			ip1=ip-1;
513			iq1=iq-1;
514			for(int j=0;j<=ip1;j++)
515			{
516			g=a(j,ip);
517			h=a(j,iq);
518			a(j,ip)=g-s(h+gtau);
519			a(j,iq)=h+s(g-htau);
520			}
521			for(int j=ip+1;j<=iq1;j++)
522			{
523			g=a(ip,j);
524			h=a(j,iq);
525			a(ip,j)=g-s(h+gtau);
526			a(j,iq)=h+s(g-htau);
527			}
528			for(int j=iq+1;j<n;j++)
529			{
530			g=a(ip,j);
531			h=a(iq,j);
532			a(ip,j)=g-s(h+gtau);
533			a(iq,j)=h+s(g-htau);
534			}
535			for(int j=0;j<n;j++)
536			{
537			g=v(j,ip);
538			h=v(j,iq);
539			v(j,ip)=g-s(h+gtau);
540			v(j,iq)=h+s(g-htau);
541			}
542			++nrot;
543			}
544			}
545			for(int ip=0;ip<n;ip++)
546			{
547			b[(unsigned int)ip]=b[(unsigned int)ip]+z[(unsigned int)ip];
548			d[(unsigned int)ip]=b[(unsigned int)ip];
549			z[(unsigned int)ip]=0.0;
550			}
551			}
552
553
554			delete [] b;
555			delete [] z;
556
557	souaissa	58	}
558	francois	17
559	souaissa	69
560	souaissa	58
561	souaissa	80
562	francois	102	int OT_TENSEUR::listes_equivalentes(std::vector<double2>& list1,std::vector<double2>& list2,std::vector<unsigned int>& indx)
563	souaissa	80
564	francois	102	{
565	souaissa	151	int compt=0;
566			int pos=-1;
567			int exist=0;
568			for(int i=0;i<list1.size();i++)
569			{
570			double2 val1=list1[i];
571			double prc=1e-3;
572			double prcc;
573			if(pos==-1)
574			{
575			for(int k=0;k< list2.size();k++)
576			{
577			double2 val2=list2[k];
578			prcc =val2.get_x()-val1.get_x();
579	francois	163	if (fabs(val2.get_x())>val2.get_dx())
580			prcc =prcc/val2.get_x();
581	souaissa	151	if( fabs(prcc)<prc){
582			indx.insert(indx.end(),k);
583			pos=k;
584			compt++;
585			break;
586			}
587			}
588			}
589	souaissa	80
590
591	souaissa	151	else {
592			for( int r=0;r<list2.size();r++)
593			{
594			double2 vval2=list2[r];
595			for(int j=0;j<indx.size();j++)
596			{
597			exist=0;
598			if (r==indx[j])
599			{
600			exist=1;
601			break;
602			}
603			}
604	souaissa	80
605	souaissa	151	prcc =vval2.get_x()-val1.get_x();
606	francois	163	if(fabs(vval2.get_x())>vval2.get_dx())
607	souaissa	151	prcc =prcc/vval2.get_x();
608	souaissa	80
609	souaissa	151	if(fabs(prcc)<prc&&!exist) // vval2==val1
610			{
611			indx.insert(indx.end(),r);
612			pos=r;
613			compt++;
614			break;
615			}
616	souaissa	80
617	souaissa	151	}
618	francois	102	}
619			}
620
621			if (compt==list1.size())return 1;
622			else return 0;
623
624
625			}
626
627
628
629
630
631