outil/src/ot_tenseur.cpp

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

#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=NULL;
}

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(std::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(std::vector<OT_VECTEUR_4DD>&v1,std::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;
    std::vector<int> check;

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

        std::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))
            {
                std::vector<int>trouve;
                std::vector<double2>valeur;
                std::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 std::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;
}


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