fct_taille_face.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/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       fct_taille_face.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:52 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
#include <math.h>
#include "fct_taille_face.h"
#include "ot_mathematique.h"
 
 
 
 
 
FCT_TAILLE_FACE::FCT_TAILLE_FACE(double eps,double dist,MG_FACE* mgface,LISTE_FRONTIERE *segment_frontiere):FCT_TAILLE_METRIQUE(),epsilon(eps),distance_maximale(dist),face(mgface),frontiere(segment_frontiere)
{
    double periode_u=mgface->get_surface()->get_periode_u();
    double periode_v=mgface->get_surface()->get_periode_v();
    decalage=new OT_DECALAGE_PARAMETRE(periode_u,periode_v);
}
 
 
FCT_TAILLE_FACE::~FCT_TAILLE_FACE()
{
    delete decalage;
}
 
void FCT_TAILLE_FACE::evaluer(double *param,double *resultat)
{
    double M1,M2,M3;
    face->get_M(param,M1,M2,M3);
    double facteur_eps;
    if (!(OPERATEUR::egal(M1+2.*M2+M3,0.,0.000001)))
        facteur_eps=sqrt(4.5*epsilon/(M1+2.*M2+M3));
    else
        facteur_eps=distance_maximale*2.;
    double facteur_dist=distance_maximale;
    double h=std::min(facteur_eps,facteur_dist);
    resultat[0]=1./h/h;
    resultat[1]=0.;
    resultat[2]=0.;
    resultat[3]=0.;
    resultat[4]=1./h/h;
    resultat[5]=0.;
    resultat[6]=0.;
    resultat[7]=0.;
    resultat[8]=1./h/h;
}
 
double  FCT_TAILLE_FACE::evaluer_facteur_distance_maximale(double *param)
{
    double du=decalage->calcul_decalage_parametre_u(param[0]);
    double dv=decalage->calcul_decalage_parametre_v(param[1]);
    double u=decalage->decalage_parametre_u(param[0],du);
    double v=decalage->decalage_parametre_v(param[1],dv);
    double E,F,G;
    face->get_EFG(param,E,F,G);
    TPL_MAP_ENTITE<MG_SEGMENT*> liste_trouvee;
    double dist_reference=4.*sqrt(E*G-F*F)*distance_maximale;
    frontiere->rechercher(u-du,v-dv,dist_reference,liste_trouvee);
    FRONTIERE_INFLUENTE frontiere_recherchee;
    int nb_liste_trouvee=liste_trouvee.get_nb();
    for (int i=0;i<nb_liste_trouvee;i++)
    {
        MG_SEGMENT* segment=liste_trouvee.get(i);
        MG_NOEUD* noeud1=segment->get_noeud1();
        MG_NOEUD* noeud2=segment->get_noeud2();
        double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
        double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
        double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
        double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
        double px=(u2-u1)*(u-u1)+(v2-v1)*(v-v1);
        double py=(u1-u2)*(u-u2)+(v1-v2)*(v-v2);
        double eps=0.0001*(u1-u2)*(u1-u2)+0.0001*(v1-v2)*(v1-v2);
        double distance;
        if ((px>(-eps)) && (py>(-eps)))
        {
            double det=(u1-u2)*(u1-u2)+(v1-v2)*(v1-v2);
            double ui=(u1*(v2-v1)*(v2-v1)+(u1-u2)*(v2-v1)*(v1-v)+u*(u1-u2)*(u1-u2))/det;
            double vi=(v1*(u1-u2)*(u1-u2)+(u1-u2)*(v1-v2)*(u-u1)+v*(v1-v2)*(v1-v2))/det;
            distance=calcule_distance(u,v,ui,vi,du,dv);
        }
        else distance=0.5*(calcule_distance(u,v,u1,v1,du,dv)+calcule_distance(u,v,u2,v2,du,dv));
        std::pair<const double,MG_SEGMENT*> tmp(distance,segment);
        frontiere_recherchee.insert(tmp);
    }
    int nombre=0;
    double numo=0.,deno=0.;
    double dismin;
    for (FRONTIERE_INFLUENTE::iterator j=frontiere_recherchee.begin();j!=frontiere_recherchee.end();j++)
    {
        double dens;
        if (((*j).first<dist_reference) && (nombre<4))
        {
            nombre++;
            MG_SEGMENT* segment=(*j).second;
            MG_NOEUD* noeud1=segment->get_noeud1();
            MG_NOEUD* noeud2=segment->get_noeud2();
            double dens1,dens2;
            for (int k1=0;k1<noeud1->get_lien_segment()->get_nb();k1++)
            {
                if (k1==0) dens1=noeud1->get_lien_segment()->get(k1)->get_longueur();
                else dens1=std::min(dens1,noeud1->get_lien_segment()->get(k1)->get_longueur());
            }
            for (int k2=0;k2<noeud2->get_lien_segment()->get_nb();k2++)
            {
                if (k2==0) dens2=noeud2->get_lien_segment()->get(k2)->get_longueur();
                else dens1=std::min(dens2,noeud2->get_lien_segment()->get(k2)->get_longueur());
            }
            dens=0.5*(dens1+dens2);
            if (nombre==1)
            {
                dismin=(*j).first;
            }
 
            numo=numo+(4*distance_maximale-(*j).first)*dens;
            deno=deno+4*distance_maximale-(*j).first;
        }
    }
    if (nombre!=0)
    {
        numo=numo+dismin*distance_maximale;
        deno=deno+dismin;
    }
    else return distance_maximale;
    return (numo/deno);
}
 
 
double FCT_TAILLE_FACE::calcule_distance(double u1,double v1,double u2,double v2,double du,double dv)
{
    double longueur_calculee=0.;
 
    double dt=0.03125;
    double ti;
    double tii=0;
    double ui=u1;
    double vi=v1;
    double uii;
    double vii;
    double dui;
    double dvi;
 
    while (tii<1.)
    {
        double param_ti[2]={ui-du,vi-dv};
        double Ei,Fi,Gi;
        face->get_EFG(param_ti,Ei,Fi,Gi);
        ti=tii;
        tii=ti+dt;
        uii=u1+tii*(u2-u1);
        vii=v1+tii*(v2-v1);
        dui=uii-ui;
        dvi=vii-vi;
        double t=0.7886751345*ti+0.2113248654*tii;
        double u=u1+t*(u2-u1);
        double v=v1+t*(v2-v1);
        double param_integration1[2]={u-du,v-dv};
        double E,F,G;
        face->get_EFG(param_integration1,E,F,G);
        longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(E*dui*dui+2*F*dui*dvi+G*dvi*dvi);
        t=0.7886751345*tii+0.2113248654*ti;
        u=u1+t*(u2-u1);
        v=v1+t*(v2-v1);
        double param_integration2[2]={u-du,v-dv};
        face->get_EFG(param_integration2,E,F,G);
        longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(E*dui*dui+2*F*dui*dvi+G*dvi*dvi);
        ui=uii;
        vi=vii;
    }
    return longueur_calculee;
}
 
void  FCT_TAILLE_FACE::deriver(double *param,double *resultat,int num_param)
{
    double uv[2];
    uv[0]=param[0];
    uv[1]=param[1];
    uv[num_param]=uv[num_param]+1e-6;
    double taille[9],taille_eps[9];
    evaluer(param,taille);
    evaluer(uv,taille_eps);
    resultat[0]=(taille_eps[0]-taille[0])/1e-6;
    resultat[1]=(taille_eps[1]-taille[1])/1e-6;
    resultat[2]=(taille_eps[2]-taille[2])/1e-6;
    resultat[3]=(taille_eps[3]-taille[3])/1e-6;
    resultat[4]=(taille_eps[4]-taille[4])/1e-6;
    resultat[5]=(taille_eps[5]-taille[5])/1e-6;
    resultat[6]=(taille_eps[6]-taille[6])/1e-6;
    resultat[7]=(taille_eps[7]-taille[7])/1e-6;
    resultat[8]=(taille_eps[8]-taille[8])/1e-6;
}
 
void FCT_TAILLE_FACE::evaluer_decompose(double *metrique_depart,double *metrique_decompose)
{
    const double precision = 1e-13;
    double v1[2];
    double v2[2];
    double lambda1,lambda2;
    //     ( a11  a12 )     ( a   b )     ( m(1)  m(3) )
    //     (          )  =  (       )  =  (            )
    //     ( a21  a22 )     ( c   d )     ( m(2)  m(4) )
 
    // Puisque tenseur est un RM_TENSEUR2X2_SYM, b et c sont egaux.
    double a = metrique_depart[0];
    double b = metrique_depart[1];
    double c = metrique_depart[2];
    double d = metrique_depart[3];
 
    double coef1  = (a + d) * (a + d);
    double coef2  = 4.0 * (a*d - b*c);
    double discriminant = coef1 - coef2;
    // NOTE: une matrice 2v2[0] est dite positive quand discriminant est positif!
 
    //
    //     CAS 1
    //
    //     Si disc est resolument positif (plus que la valeur de troncature
    //     prec), ce qui sera presque toujours le cas, alors la matrice a
    //     deux valeurs propres reelles distinctes. A chacune de ces valeurs
    //     propres correspond un vecteur propre non-nul. Une valeur propre
    //     peut etre nulle (si a*d - b*c = 0, si det(m) = 0) mais dans ce
    //     cas-ci, on s'en fout, on obient toujours deux vecteurs propres
    //     distincts et non-nuls.
    //
    if ( discriminant > precision * ( fabs(coef1) + fabs(coef2) ) )
    {
        discriminant = sqrt ( discriminant );
        lambda1   = ( a + d - discriminant ) * 0.5;
        lambda2   = ( a + d + discriminant ) * 0.5;
        v1[0] = b;
        v1[1] = lambda1 - a;
        if ( fabs(v1[1]) <=  precision * ( fabs(lambda1) + fabs(a) ) )
        {
            v1[0] = lambda1 - d;
            v1[1] = c;
        }
        double normV1 = 1.0 / sqrt ( v1[0]*v1[0] + v1[1]*v1[1] );
        v1[0]   = v1[0] * normV1;
        v1[1]   = v1[1] * normV1;
        v2[0] = b;
        v2[1] = lambda2 - a;
        if ( fabs(v2[1]) <= precision * ( fabs(lambda2) + fabs(a) ) )
        {
            v2[0] = lambda2 - d;
            v2[1] = c;
        }
        double normV2 = 1.0 / sqrt ( v2[0]*v2[0] + v2[1]*v2[1] );
        v2[0]   = v2[0] * normV2;
        v2[1]   = v2[1] * normV2;
    }
    else if ( fabs(discriminant) <= precision*( fabs(coef1) + fabs(coef2) ) )
    {
        if ( fabs(b) + fabs(c) > precision * ( fabs(a) + fabs(d) ) )
        {
            lambda1 = ( a + d ) * 0.5;
            lambda2 = lambda1;
            if ( fabs(b) > fabs(c) )
            {
                v1[0] = 1;
                v1[1] = - ( a - d ) * v1[0] / ( 2.0 * b );
            }
            else
            {
                v1[1] = 1;
                v1[0] = ( a - d ) * v1[1] / ( 2.0 * c );
            }
            double norm   = 1.0 / sqrt ( v1[0]*v1[0] + v1[1]*v1[1] );
            v1[0]     = v1[0] * norm;
            v1[1]     = v1[1] * norm;
            v2[0]     = -v1[1];
            v2[1]     =  v1[0];
        }
        else
        {
            if ( a > d )
            {
                lambda1 = a;
                v1[0] = 1;
                v1[1] = 0;
                lambda2 = d;
                v2[0] = 0;
                v2[1] = 1;
            }
            else
            {
                lambda1 = d;
                v1[0] = 0;
                v1[1] = 1;
                lambda2 = a;
                v2[0] = 1;
                v2[1] = 0;
            }
        }
    }
    else
    {
        if ( fabs(b) > fabs(c) )
        {
            lambda1 = b;
            lambda2 = b;
        }
        else
        {
            lambda1 = c;
            lambda2 = c;
        }
        v1[0] = 1;
        v1[1] = 0;
        v2[0] = 0;
        v2[1] = 1;
    }
 
    metrique_decompose[0]=lambda1;
    metrique_decompose[1]=v1[0];
    metrique_decompose[2]=v1[1];
    metrique_decompose[3]=lambda2;
    metrique_decompose[4]=v2[0];
    metrique_decompose[5]=v2[1];
}
 
int FCT_TAILLE_FACE::valide_parametre(double *param)
{
    return 1;
}
 
 
double FCT_TAILLE_FACE::get_valeur_maximale(int num)
{
    return distance_maximale;
}