mg_geodesic.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/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mg_geodesic.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:56 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
#include "gestionversion.h"
#include "mg_geodesic.h"
#include "mg_fast_marching2D.h"
#include "mg_gestionnaire.h"
#include "fem_maillage.h"
#include "fem_triangle3.h"
#include "fem_solution.h"
#include "tpl_map_entite.h"
#include "tpl_liste_unique_entite.h"
#include "ot_mathematique.h"
#include <math.h>
#include "pars_argument.h"
#include "parse.h"
#include "fem_maillage_outils.h"
#include "ot_boite_3d.h"
 
 
namespace MAGIC
{
#define INFINI 1e300
 
  
MG_GEODESIC::MG_GEODESIC(char* chaine,MG_GESTIONNAIRE *g,FEM_MAILLAGE *f,int nbpt):MAGIC_AFFICHE(),fem(f),gest(g)
{
    char chaine2[500];
    sprintf(chaine2,"%s0_fm.sol",chaine);
    sol=gest->get_fem_solution(fem,nbpt+2,chaine2,(char*)"Geodesic");
    if (sol==NULL)
          {
          sol=new FEM_SOLUTION(fem,nbpt+2,chaine2,1,"Geodesic");
          for (int i=0;i<nbpt;i++)
        {
        char numlevel[500];
       sprintf(numlevel,"L%d",i+1);
        sol->change_legende(i,numlevel);
        }
          sol->change_legende(nbpt,(char*)"Li+Lj");
          sol->change_legende(nbpt+1,(char*)"Li-Lj");
          gest->ajouter_fem_solution(sol);
          }
/*    sprintf(chaine2,"%s1_fm.sol",chaine);
    solgrad=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
    if (solgrad==NULL)
          {
          solgrad=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::VECTEUR);
          solgrad->change_legende(0,(char*)"LA'");
          solgrad->change_legende(1,(char*)"LB'");
          gest->ajouter_fem_solution(solgrad);
          }
    sprintf(chaine2,"%s2_fm.sol",chaine);
    solgradnoeud=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
    if (solgradnoeud==NULL)
          {
          solgradnoeud=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::VECTEUR);
          solgradnoeud->change_legende(0,(char*)"LA'_ext");
          solgradnoeud->change_legende(1,(char*)"LB'_ext");
          gest->ajouter_fem_solution(solgradnoeud);
          }
    sprintf(chaine2,"%s3_fm.sol",chaine);
    solgradgrad=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
    if (solgradgrad==NULL)
          {
          solgradgrad=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::TENSEUR);
          solgradgrad->change_legende(0,(char*)"LA''");
          solgradgrad->change_legende(1,(char*)"LB''");
          gest->ajouter_fem_solution(solgradgrad);
          }
    sprintf(chaine2,"%s4_fm.sol",chaine);
    courburemax=gest->get_fem_solution(fem,2,chaine2,(char*)"Geodesic");
    if (courburemax==NULL)
          {
          courburemax=new FEM_SOLUTION(fem,2,chaine2,1,"Geodesic",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2,MAGIC::TYPE_SOLUTION::SCALAIRE);
          courburemax->change_legende(0,(char*)"CAmax");
          courburemax->change_legende(1,(char*)"CBmax");
          gest->ajouter_fem_solution(courburemax);
          }*/
    
}
 
 
 
MG_GEODESIC::MG_GEODESIC(MG_GEODESIC &mdd):sol(mdd.sol),fem(mdd.fem)
{
}
 
MG_GEODESIC::~MG_GEODESIC()
{
}
 
 
 
 
bool MG_GEODESIC::init_noeud(char* entite)
{
PARS_ARGUMENT param[100];
PARSE parse;
parse.decode(entite,"&",param);
int nb=param[0].argument.size();
MG_GEOMETRIE* geo=fem->get_mg_geometrie();
for (int i=0;i<nb;i++)
    {
    int num=atoi(param[0].argument[i].c_str());
    MG_ELEMENT_TOPOLOGIQUE *ele=NULL;
    ele=geo->get_mg_sommetid(num);
    if (ele!=NULL)
        {
        int nbele=ele->get_lien_fem_maillage()->get_nb();
        for (int j=0;j<nbele;j++)
            {
            FEM_NOEUD* no=(FEM_NOEUD*)ele->get_lien_fem_maillage()->get(j);
            if (fem->get_fem_noeudid(no->get_id())==no) lst.push_back(no);
            }
        }
    
      }
initdouble=false;
int nbnoeudinit=lst.size();
if (nbnoeudinit!=2) return false;
return true; 
}
bool MG_GEODESIC::init_coord(char* entite)
{
PARS_ARGUMENT param[100];
PARSE parse;
parse.decode(entite,"@,@,@,@,@,@",param);
lstdouble.push_back(atof(param[0].argument[0].c_str()));
lstdouble.push_back(atof(param[1].argument[0].c_str()));
lstdouble.push_back(atof(param[2].argument[0].c_str()));
lstdouble.push_back(atof(param[3].argument[0].c_str()));
lstdouble.push_back(atof(param[4].argument[0].c_str()));
lstdouble.push_back(atof(param[5].argument[0].c_str()));
initdouble=true;
return true;
}
 
 
 
 
double MG_GEODESIC::calcul_fm(int num,double xdep,double ydep,double zdep,std::vector<double> &lstpoints)
{
MG_FAST_MARCHING2D fm(sol,num);
if (affichageactif) fm.active_affichage(affiche2);
fm.init();
if (initdouble==false) 
  {
  std::vector<FEM_NOEUD*> lst2;
  lst2.push_back(lst[num]);     
  fm.init_noeud(lst2);
  }
else
  fm.init_coord(xdep,ydep,zdep);
fm.propage(1.0);
fm.finalise();
sol->active_solution(num);
double valA;
if (initdouble==false) valA=lst[(num+1)%2]->get_solution();
else 
    {
    int nbval=lstpoints.size()/3;
    for (int i=0;i<nbval;i++)
        {
        valA=fm.get_valeur(lstpoints[3*i+0],lstpoints[3*i+1],lstpoints[3*i+2]);
        lstpoints.push_back(valA);
        }
    }
return valA;    
}
 
void MG_GEODESIC::calcul_reseau(std::vector<double> lstpoint,double *matdist)
{
initdouble=true;
int nb=lstpoint.size()/3;
for (int i=0;i<nb;i++)
    for (int j=0;j<nb;j++)
        matdist[i*nb+j]=0;
for (int i=0;i<nb;i++)
    {
     std::vector<double> lst=lstpoint;
     lst.erase(lst.begin()+3*i,lst.begin()+3*i+3);
     calcul_fm(i,lstpoint[3*i],lstpoint[3*i+1],lstpoint[3*i+2],lst);
    int ii=0;
     for (int j=0;j<nb;j++)
        if (j!=i) 
            {
            double val=lst[3*nb-3+ii];
            matdist[i*nb+j]=val;
            //matdist[j*nb+i]=matdist[j*nb+i]+lst[3*nb-3+ii];
            ii++;
            }
    }
for (int i=0;i<nb;i++)
    for (int j=i+1;j<nb;j++)
        if (i!=j)
            {
            double vala,valb;
            calcul_milieux(i,j,nb,nb+1,vala,valb);  
            matdist[i*nb+j]=matdist[i*nb+j]+matdist[j*nb+i]-vala-valb;
            matdist[j*nb+i]=matdist[i*nb+j];
            }
}
 
double MG_GEODESIC::calcul(double *tab)
{       
std::vector<double> lst;
double xdep,ydep,zdep;
if (initdouble==true) 
    {
    lst.push_back(lstdouble[3]);
    lst.push_back(lstdouble[4]),
    lst.push_back(lstdouble[5]);
    xdep=lstdouble[0];
    ydep=lstdouble[1];
    zdep=lstdouble[2];
    }
double valA=calcul_fm(0,xdep,ydep,zdep,lst);
if (initdouble==true) 
    {
    lst.clear();
    lst.push_back(lstdouble[0]);
    lst.push_back(lstdouble[1]),
    lst.push_back(lstdouble[2]);
    xdep=lstdouble[3];
    ydep=lstdouble[4];
    zdep=lstdouble[5];
    }
double valB=calcul_fm(1,xdep,ydep,zdep,lst);
double vala,valb;
calcul_milieux(0,1,2,3,vala,valb);
double longueur=valA+valB-vala-valb;
if (tab!=NULL)
    {
    tab[0]=valA;  
    tab[1]=valB;  
    tab[2]=longueur;  
    }
return longueur;
}
 
void MG_GEODESIC::calcul_milieux(int num1,int num2,int num3, int num4,double &vala,double &valb)
{
FEM_MAILLAGE_OUTILS ot;
ot.operation_champs_solution(sol,num1,sol,num2,sol,num3,MAGIC::OPERATION_FEM_SOLUTION::ADDITION);
ot.operation_champs_solution(sol,num1,sol,num2,sol,num4,MAGIC::OPERATION_FEM_SOLUTION::SOUSTRACTION);
sol->active_solution(num4);
 
int ligne=0;
LISTE_FEM_NOEUD::iterator itn;
for (FEM_NOEUD* no=fem->get_premier_noeud(itn);no!=NULL;no=fem->get_suivant_noeud(itn))
  {
  double val=sol->lire(ligne,num3);
  no->change_solution(val,1);
  val=sol->lire(ligne,num1);
  no->change_solution(val,2);
  val=sol->lire(ligne,num2);
  no->change_solution(val,3);
  ligne++;
  }
double limit=0.;
double minApB=INFINI;
LISTE_FEM_ELEMENT2::iterator it;
double xm,ym,zm;
for (FEM_ELEMENT2 *ele=fem->get_premier_element2(it);ele!=NULL;ele=fem->get_suivant_element2(it))
    {
    double val1=ele->get_fem_noeud(0)->get_solution();  
    double val2=ele->get_fem_noeud(1)->get_solution();
    double val3=ele->get_fem_noeud(2)->get_solution();  
    if ((val1-limit)*(val2-limit)<1e-12)
            {
            double t=(limit-val1)/(val2-val1);
            double val1apb=ele->get_fem_noeud(0)->get_solution(1);
            double val2apb=ele->get_fem_noeud(1)->get_solution(1);
            double valApB=val1apb+t*(val2apb-val1apb);
            if (valApB<minApB) 
            {
            minApB=valApB;
            double val1a=ele->get_fem_noeud(0)->get_solution(2);
            double val2a=ele->get_fem_noeud(1)->get_solution(2);
            vala=val1a+t*(val2a-val1a);
                double val1b=ele->get_fem_noeud(0)->get_solution(3);
            double val2b=ele->get_fem_noeud(1)->get_solution(3);
            valb=val1b+t*(val2b-val1b);
            xm=ele->get_fem_noeud(0)->get_x()+t*(ele->get_fem_noeud(1)->get_x()-ele->get_fem_noeud(0)->get_x());
            ym=ele->get_fem_noeud(0)->get_y()+t*(ele->get_fem_noeud(1)->get_y()-ele->get_fem_noeud(0)->get_y());
            zm=ele->get_fem_noeud(0)->get_z()+t*(ele->get_fem_noeud(1)->get_z()-ele->get_fem_noeud(0)->get_z());
            }
            }
    if ((val2-limit)*(val3-limit)<1e-12)
            {
            double t=(limit-val2)/(val3-val2);
            double val2apb=ele->get_fem_noeud(1)->get_solution(1);
            double val3apb=ele->get_fem_noeud(2)->get_solution(1);
            double valApB=val2apb+t*(val3apb-val2apb);
            if (valApB<minApB) 
            {
            minApB=valApB;
            double val2a=ele->get_fem_noeud(1)->get_solution(2);
            double val3a=ele->get_fem_noeud(2)->get_solution(2);
            vala=val2a+t*(val3a-val2a);
                double val2b=ele->get_fem_noeud(1)->get_solution(3);
            double val3b=ele->get_fem_noeud(2)->get_solution(3);
            valb=val2b+t*(val3b-val2b);
            xm=ele->get_fem_noeud(1)->get_x()+t*(ele->get_fem_noeud(2)->get_x()-ele->get_fem_noeud(1)->get_x());
            ym=ele->get_fem_noeud(1)->get_y()+t*(ele->get_fem_noeud(2)->get_y()-ele->get_fem_noeud(1)->get_y());
            zm=ele->get_fem_noeud(1)->get_z()+t*(ele->get_fem_noeud(2)->get_z()-ele->get_fem_noeud(1)->get_z());
            }
            }
     if ((val1-limit)*(val3-limit)<1e-12)
            {
            double t=(limit-val1)/(val3-val1);
            double val1apb=ele->get_fem_noeud(0)->get_solution(1);
            double val3apb=ele->get_fem_noeud(2)->get_solution(1);
            double valApB=val1apb+t*(val3apb-val1apb);
            if (valApB<minApB) 
            {
            minApB=valApB;
            double val1a=ele->get_fem_noeud(0)->get_solution(2);
            double val3a=ele->get_fem_noeud(2)->get_solution(2);
            vala=val1a+t*(val3a-val1a);
                double val1b=ele->get_fem_noeud(0)->get_solution(3);
            double val3b=ele->get_fem_noeud(2)->get_solution(3);
            valb=val1b+t*(val3b-val1b);
            xm=ele->get_fem_noeud(0)->get_x()+t*(ele->get_fem_noeud(2)->get_x()-ele->get_fem_noeud(0)->get_x());
            ym=ele->get_fem_noeud(0)->get_y()+t*(ele->get_fem_noeud(2)->get_y()-ele->get_fem_noeud(0)->get_y());
            zm=ele->get_fem_noeud(0)->get_z()+t*(ele->get_fem_noeud(2)->get_z()-ele->get_fem_noeud(0)->get_z());
                        }
            } 
 
    }
char mess[500];
sprintf(mess,"    Le point milieu x=%lf, y=%lf, z=%lf la demi longueur de la LS1=%f la demi longueur de la LS2=%f ",xm,ym,zm,vala,valb) ;
affiche(mess);    
}
    
};