outil/src/ot_tenseur.cpp

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       ot_tenseur.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:54:00 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------

#pragma hdrstop

#include "ot_tenseur.h"

#include <iomanip>
#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;
}

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

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