ot_tenseur.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  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-2025
//####//  Derniere modification par francois
//####//  mer 14 mai 2025 17:54:47 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
#pragma hdrstop
 
#include "ot_tenseur.h"
 
#include <iomanip>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>
 
#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 operator*( double& coef, OT_TENSEUR& ot)
{
    double2 coef2(coef);
    return coef2*ot;
 
}
 
OT_TENSEUR operator*(OT_TENSEUR& ot , double& coef)
{
    double2 coef2(coef);
    return coef2*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::inverse(void)
{
double2 zero(0.),eps(1e-14),un(1.);    
OT_TENSEUR tnul;
if (Dim1!=Dim2) return tnul;
int n=Dim1;
OT_TENSEUR tinv(n,n);
for (int i=0;i<n;i++)
    for (int j=0;j<n;j++)
        tinv(i,j)=(*this)(i,j);
double2 pd=un;
for (int l=0;l<n;l++)
    {
    double2 dd=zero;
    for (int k=0;k<n;k++)
        dd=dd+tinv(l,k)*tinv(l,k);
    dd=sqrt(dd);
    pd=pd*dd;
    }
    
std::vector<int> j(10000);
double2 detm=1.0;
for (int l=0;l<n;l++)
    j[l+20+n]=l;
for (int l=0;l<n;l++)
    {
    double2 cc=0.;
    int m=l;
    for (int k=l;k<n;k++)
        if (f2abs(cc)-f2abs(tinv(l,k))<zero)
            {
            m=k;
            cc=tinv(l,k);
                
            }
    if (l!=m) 
        {
        int k=j[m+20+n];
        j[m+20+n]=j[l+20+n];
        j[l+20+n]=k;
        for (int k=0;k<n;k++)
            {
            double2 s=tinv(k,l);
            tinv(k,l)=tinv(k,m);
            tinv(k,m)=s;
            }
        }
    tinv(l,l)=1.;
    detm=detm*cc;    
    for (int m=0;m<n;m++)
        tinv(l,m)=tinv(l,m)/cc;
    for (int m=0;m<n;m++)
        if (l!=m) 
            {   
            cc=tinv(m,l);
            if (f2abs(cc)>eps)
                {
                tinv(m,l)=zero;
                for (int k=0;k<n;k++)
                    tinv(m,k)=tinv(m,k)-cc*tinv(l,k);
                }
            }
        
    }
for (int l=0;l<n;l++)
    if (j[l+20+n]!=l)
    {
    int m=l;
    int ok;
    do 
        {
        ok=1;    
        m++;    
        if (j[m+20+n]!=l)
            if (n>m) ok=0;
    
        }
    while (ok==0);
    j[m+20+n]=j[l+20+n];
    for (int k=0;k<n;k++)
          {
          double2 cc=tinv(l,k);    
          tinv(l,k)=tinv(m,k);
          tinv(m,k)=cc;
          }
    j[l+20+n]=l;
        
        
    }
return tinv;
}
 
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;
}
 
double2 operator%( OT_TENSEUR& ot1, OT_TENSEUR& ot2)
{
    int dim11=ot1.Get_NbRows();
    int dim12= ot1.Get_NbCols();
    int dim21=ot2.Get_NbRows();
    int dim22=ot2.Get_NbCols();
    if (dim11!=dim21) return 0.;
    if (dim12!=dim22) return 0.;
    double2 res=0.;
 
    for (int i=0;i<dim11;i++)
        for (int j=0;j<dim12;j++)
            res=res+ot1(i,j)*ot2(i,j);
return res;
}
 
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;
 
 
}