ot_fonctions.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_fonctions.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2025
//####//  Derniere modification par francois
//####//  mer 14 mai 2025 17:54:47 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "ot_fonctions.h"
 
 
double2 moment (double2 theta,std::vector<OT_VECTEUR_4DD>& ctrpts,int plan)
{
 
 
    int nb_pts=ctrpts.size();
    double ordre=4.;
    double2 zer=0.;
    double2 sum=0.;
    for (int i=0;i<nb_pts;i++)
    {
        //thi[i]=thi[i]*PI/180;
        //double2 nori=ctrpts[i].norme();
        double2 x=ctrpts[i].x()^2;
        double2 y=ctrpts[i].y()^2;
        double2 z=ctrpts[i].z()^2;
        double2 nori;
 
        // double2 nori=ctrpts[i].norme();
        double2 cosi=0.;
        double2 sini=0.;
        if (plan==12)
        {
            nori=(x+y)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_x()/nori;
                sini=ctrpts[i].get_y()/nori;
            }
        }
        if (plan==23)
        {
            nori=(y+z)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_y()/nori;
                sini=ctrpts[i].get_z()/nori;
            }
 
        }
        if (plan==13)
        {
            nori=(x+z)^0.5;
            if (nori!=zer) {
                cosi=ctrpts[i].get_x()/nori;
                sini=ctrpts[i].get_z()/nori;
            }
        }
        double2 costh=cos(theta);
        double2 sinth=sin(theta);
        double2 d_theta=sini*costh-cosi*sinth;
        //sum=sum+ pow(li[i],4)*pow(d_theta,4);
        double2 val1=nori^ordre;
        double2 val2=d_theta^ordre;
        double2 val=val1*val2;
        sum=sum+ val;
    }
    return sum;
}
 
 
double2  Max(double2 a,double2 b)
{
    return((a<b)? b : a);
}
 
double2  Min(double2 a,double2 b)
{
    return ((a>b)? b : a);
}
 
double2 Sign(double2 a,double2 b)
{
    double2 zer=0.;
    double2 c=-1.;
    if (b>=zer)
        return a;
    else
        return c*a;
}
 
void swap(double2& a,double2& b)
{
    double2 tmp;
    tmp = a;
    a = b;
    b = tmp;
}
void Shift(double2&a,double2&b,double2&c,double2 d)
{
    a=b;
    b=c;
    c=d;
}
 
void min_Brak(double2 &ax,double2 &bx,double2 &cx,double2 &fa,double2 &fb,double2 &fc,std::vector<OT_VECTEUR_4DD>& ctrpts,int plan)
{
    static const double2 GOLD=1.618034;
    static const double2 GLIMIT=100.0;
    static const double2 TINY=1.0e-20;
    double2 ulim,u,r,q,fu;
    double2(*func)(double2,std::vector<OT_VECTEUR_4DD>&,int)=&moment;
    double2 dtheta=0.0001;
    ax=0.;
    bx=ax+dtheta;
    cx=bx+dtheta;
    double2 fva=(*func)(ax,ctrpts,plan);
    double2 fvb=(*func)(bx,ctrpts,plan);
    double2 fvc=(*func)(cx,ctrpts,plan);
 
    while (fvb> fva)
    {
        ax=bx;
        bx=ax+dtheta;
        fva=(*func)(ax,ctrpts,plan);
        fvb=(*func)(bx,ctrpts,plan);
    }
 
    if (fvb<fva)
    {
        cx=bx+dtheta;
        while (fvb> fvc)
        {
            bx=cx;
            cx=bx+dtheta;
            fvb=(*func)(bx,ctrpts,plan);
            fvc=(*func)(cx,ctrpts,plan);
        }
    }
 
    fa=(*func)(ax,ctrpts,plan);
    fb=(*func)(bx,ctrpts,plan);
    if (fb>fa)
    {
        swap(ax,bx);
        swap(fb,fa);
    }
    cx=bx+GOLD*(bx-ax);
    fc=(*func)(cx,ctrpts,plan);
    double2 d=2.;
    double2 zer=0.;
    while (fb>fc)
    {
        r=(bx-ax)*(fb-fc);
        q=(bx-cx)*(fb-fa);
        u=bx-((bx-cx)*q-(bx-ax)*r)/(d*Sign(Max((q-r).get_fabs(),TINY),q-r));//(d*Sign(Max(abs(q-r),TINY),q-r));
        ulim=bx+GLIMIT*(cx-bx);
        if ((bx-u)*(u-cx)>zer)
        {
            fu=(*func)(u,ctrpts,plan);
            if (fu<fc)
            {
                ax=bx;
                bx=u;
                fa=fb;
                fb=fu;
                return;
            }
            else if (fu>fb)
            {
                cx=u;
                fc=fu;
                return;
            }
            u=cx+GOLD*(cx-bx);
            fu=(*func)(u,ctrpts,plan);
        }
        else if ((cx-u)*(u-ulim)>zer)
        {
            fu=(*func)(u,ctrpts,plan);
            if (fu<fc)
            {
                bx=cx;
                cx=u;
                u=cx+GOLD*(cx-bx);
                Shift(fb,fc,fu,(*func)(u,ctrpts,plan));
            }
        }
        else if ((u-ulim)*(ulim-cx)>=zer)
        {
            u=ulim;
            fu=(*func)(u,ctrpts,plan);
        }
        else
        {
            u=cx+GOLD*(cx-bx);
            fu=(*func)(u,ctrpts,plan);
        }
        Shift(ax,bx,cx,u);
        Shift(fa,fb,fc,fu);
    }
}
 
 
 
 
std::vector<OT_VECTEUR_4DD> system_axes(std::vector<OT_VECTEUR_4DD>& ctrpts)
{
    std::vector<OT_VECTEUR_4DD> syst_axe;
    int nb_points=ctrpts.size();
    std::vector<OT_VECTEUR_4DD> ctrpts_proj;
    OT_VECTEUR_4DD axe1,axe2,axe3,axe4;
    double2 ax,bx,cx,fa,fb,fc;
    int plan=12;
    double2 z0=ctrpts[0].get_z() ;
    double2 y0=ctrpts[0].get_z() ;
    double2 x0=ctrpts[0].get_z() ;
    int cmpt1=0;
    int cmpt2=0;
    int cmpt3=0;
    for (int i=0;i<nb_points;i++)
    {
        if (ctrpts[i].get_x()==x0) cmpt1++;
        if (ctrpts[i].get_y()==y0) cmpt2++;
        if (ctrpts[i].get_z()==z0) cmpt3++;
    }
 
    if (cmpt1==nb_points) plan=23;
    if (cmpt2==nb_points) plan=13;
    if (cmpt3==nb_points) plan=12;
 
 
 
    min_Brak(ax,bx,cx,fa,fb,fc,ctrpts,plan);
    double theta= bx.get_x();
    double c=cos(theta);
    double s=sin(theta);
    axe1.change_x(c);
    axe1.change_y(s);
    axe1.change_z(0.);
    axe1.change_w(0.);
    double e1_proj[3];
    double e2_proj[3];
    double normal[3];
    normal[0] =(axe1.get_x()).get_x();
    normal[1] =(axe1.get_y()).get_x();
    normal[2] =(axe1.get_z()).get_x();
    double root[3]={0.,0.,0.};
    OT_ALGORITHME_GEOMETRIQUE algo;
    int pris=0;
    OT_MATRICE_3D P, P_1;
 
    for (int i=0;i<nb_points;i++)
    {
        double pnt[3];
        double proj_pnt[3];
        pnt[0] =(ctrpts[i].get_x()).get_x();
        pnt[1] =(ctrpts[i].get_y()).get_x();
        pnt[2] =(ctrpts[i].get_z()).get_x();
 
 
        algo.Proj3D_Point_Plan (normal, root, pnt, proj_pnt);
 
        OT_VECTEUR_4DD pt_proj;
 
 
        if ((proj_pnt[0]!=0.||proj_pnt[1]!=0.||proj_pnt[2]!=0)&&!pris)
        {
            e1_proj[0]=proj_pnt[0];
            e1_proj[1]=proj_pnt[1];
            e1_proj[2]=proj_pnt[2];
            double norm=sqrt(pow(proj_pnt[0],2)+pow(proj_pnt[1],2)+pow(proj_pnt[2],2));
 
            e1_proj[0]=e1_proj[0]/norm;
            e1_proj[1]=e1_proj[1]/norm;
            e1_proj[2]=e1_proj[2]/norm;
 
            e2_proj[0]=normal[1]*e1_proj[2]-normal[2]*e1_proj[1];
            e2_proj[1]=normal[2]*e1_proj[0]-normal[2]*e1_proj[2];
            e2_proj[2]=normal[0]*e1_proj[1]-normal[2]*e1_proj[0];
 
            OT_VECTEUR_3D V1(e1_proj[0],e2_proj[0],normal[0]);
            OT_VECTEUR_3D V2(e1_proj[1],e2_proj[1],normal[1]);
            OT_VECTEUR_3D V3(e1_proj[2],e2_proj[2],normal[2]);
 
 
 
            P.change_vecteur1(V1);
            P.change_vecteur2(V2);
            P.change_vecteur3(V3);
 
            P_1=P.inverse();
            pris=1;
        }
 
        OT_VECTEUR_3D V_PT(proj_pnt[0],proj_pnt[1],proj_pnt[2]);
 
        OT_VECTEUR_3D V=P_1*V_PT;
 
 
        pt_proj.change_x(V[0]);
        pt_proj.change_y(V[1]);
        pt_proj.change_z(V[2]);
 
        ctrpts_proj.insert(ctrpts_proj.end(),pt_proj);
    }
 
    min_Brak(ax,bx,cx,fa,fb,fc,ctrpts_proj,plan);
 
    theta= bx.get_x();
    c=cos(theta);
    s=sin(theta);
    OT_VECTEUR_3D Vp(c,s,0.);
    double   nnv_x=Vp*P.get_vecteur1();
    double   nnv_y=Vp*P.get_vecteur2();
    double   nnv_z=Vp*P.get_vecteur3();
 
    axe2.change_x(nnv_x);
    axe2.change_y(nnv_y);
    axe2.change_z(nnv_z);
 
    double2 val=axe1*axe2;
 
    axe3[0]=axe1[1]*axe2[2]-axe1[2]*axe2[1];
    axe3[1]=axe1[2]*axe2[0]-axe1[2]*axe2[2];
    axe3[2]=axe1[0]*axe2[1]-axe1[2]*axe2[0];
 
    axe3[3]=0.;
 
 
    syst_axe.insert( syst_axe.end(),axe1);
    syst_axe.insert( syst_axe.end(),axe2);
    syst_axe.insert( syst_axe.end(),axe3);
 
    return syst_axe;
}
 
 
OT_FONCTION_SYMBOLIQUE::OT_FONCTION_SYMBOLIQUE(char* f, char *var):formule(f)
{
parseur.SetExpr(formule.c_str());
char *p=var;
bool debut=true;
for (int i=0;i<strlen(var)+1;i++)
    {
    if (debut) 
        if ((var[i]!=' ') && (var[i]!=0))
            debut=false;
    if (!debut) 
        if ((var[i]==' ') || (var[i]==0))
            {
            char v[200];
            sscanf(p,"%s",v);
            variable.push_back(v);
            p=var+i;
            debut=true;
            }
        
    }
    nb_variable=variable.size();
}
 
 
 
OT_FONCTION_SYMBOLIQUE::OT_FONCTION_SYMBOLIQUE(OT_FONCTION_SYMBOLIQUE& mdd):nb_variable(mdd.nb_variable),formule(mdd.formule)
{
}
 
 
 
OT_FONCTION_SYMBOLIQUE::~OT_FONCTION_SYMBOLIQUE()
{
}
 
int OT_FONCTION_SYMBOLIQUE::get_nb_variable()
{
return nb_variable;
}
 
std::string OT_FONCTION_SYMBOLIQUE::get_nom_variable(int num)
{
return variable[num];    
}
 
 
void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double f(double,double))
{
parseur.DefineFun(nom,f);
}
 
void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double f( double )  )
{
parseur.DefineFun(nom,f);
}
 
void OT_FONCTION_SYMBOLIQUE::def_fonction ( char* nom, double  f( double, double, double ))
{
parseur.DefineFun(nom,f);
}
 
double OT_FONCTION_SYMBOLIQUE::evalue ( double* var )
{
for (int i=0;i<nb_variable;i++)
    parseur.DefineVar(variable[i], var+i);
double result;
try 
    {
    result = parseur.Eval();    
    }
catch (mu::Parser::exception_type &e)
    {
    result=0.;
    }
return result;
}
 
double OT_FONCTION_SYMBOLIQUE::integre ( double a, double b, int nb_pas)
{
if (nb_variable!=1) return 0.;
    
double intg=0.;
for (int i=0;i<nb_pas;i++)
    {
    double ti=a+(b-a)/nb_pas*i;    
    double tii=a+(b-a)/nb_pas*(i+1.);
    double xsi=-1./sqrt(3);
    double t=(tii-ti)/2.*xsi+(tii+ti)/2.;
    double res=evalue(&t);
    intg=intg+(tii-ti)/2.*(1.*res);
    xsi=1./sqrt(3);    
    t=(tii-ti)/2.*xsi+(tii+ti)/2.;
    res=evalue(&t);
    intg=intg+(tii-ti)/2.*(1.*res);
    }
return intg;    
}