poly_voro.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/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       poly_voro.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:54:00 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
#include "poly_voro.h"
#include "poly_point.h"
#include "poly_noeud.h"
#include "poly_face.h"
#include "poly_cellule.h"
 
 
#include "tpl_grille.h"
#include "ot_mathematique.h"
 
 
#include <algorithm>
 
#include "../voro++-0.4.6/src/voro++.hh"
 
Poly_Voro::Poly_Voro(std::vector<Poly_Point*> list_pnts,double dg): POLY_AFFICHE(),list_points(list_pnts),epsfusion(dg)
{
}
 
void Poly_Voro::construit(void)
{
  // Géométrie du conteneur
  double x_min= list_points[0]->get_x(), x_max= list_points[6]->get_x();
  double y_min= list_points[0]->get_y(), y_max= list_points[6]->get_y();
  double z_min= list_points[0]->get_z(), z_max= list_points[6]->get_z();
  double cvol = (x_max - x_min) * (y_max - y_min) * (x_max - x_min);
 
  // Nombre de blocs qui divise le container (blocs calculatoires)
  int n_x= 6, n_y= 6, n_z= 6;
 
  // Cré un container avec la géométrie données
  voro::container con(x_min, x_max, y_min, y_max, z_min, z_max, n_x, n_y, n_z, false, false, false, 8);
 
  // Ajoute les particules au container, on n'insère pas les 8 coins du cube
  double x,y,z;
  for(int i=8; i < list_points.size(); i++) {
      x = list_points[i]->get_x();
      y = list_points[i]->get_y();
      z = list_points[i]->get_z();
      con.put(i, x, y, z);
  }
  
  // Additionne les volumes et vérifi que ça concorde avec le volume du container
  double vvol = con.sum_cell_volumes();
  char message[2000];
  sprintf(message,"   Volume de l'echantillon = %lf",cvol);
  affiche(message);
  voro::voronoicell voro_cell;
  voro::c_loop_all loop(con);
  int nb=0;
  if(loop.start()) do nb++;while(loop.inc());
  sprintf(message,"   Nombre de cellules de Voronoi = %d",nb);
  affiche(message);
  sprintf(message,"   Volume des cellules de Voronoi = %lf",vvol);
  affiche(message);
  
  
  
  
  
  // loop a travers les cellules
  if(loop.start()) do if(con.compute_cell(voro_cell,loop)) 
  {
    Poly_Cellule* cell = new Poly_Cellule;
    add_cell(cell);
    loop.pos(x,y,z); // position globale de la particule en cours
    
    // récupération des noeuds
    std::vector< double > vert_xyz;
    voro_cell.vertices(x, y, z, vert_xyz); // récupère tous les vertex de la cellule dans le système global
    
    for(int i=0; i < vert_xyz.size(); i+=3)
    {
      static int numeronoeud=0;
      numeronoeud++;
      Poly_Noeud* nd = new Poly_Noeud(numeronoeud, vert_xyz[i], vert_xyz[i+1], vert_xyz[i+2]);
      cell->add_noeud(nd);
    }
    
    // récupération des faces
    std::vector< int > face_vert;
    voro_cell.face_vertices(face_vert); // liste de vertex pour chaque face
    
    for(int i=0; i < face_vert.size(); i++)
    {
      Poly_Face* face = new Poly_Face;
      int nbVert = face_vert[i]; // première valeur = le nombre de vertex de la face
      for(int j=i+1; j <  i+nbVert+1; j++)
      {
        face->add_noeud(face_vert[j]);
      }
      cell->add_face(face);
      i += nbVert;
    }
  } while (loop.inc());
}
 
Poly_Voro::~Poly_Voro()
{
    for (int i=0; i<get_nb_cell(); i++)
    {
        Poly_Cellule* c = list_cellules[i];
        delete c;
    }      
}
 
void Poly_Voro::fusion_noeuds(void)
{
  TPL_MAP_ENTITE<Poly_Noeud*> lstnoeud;
  TPL_GRILLE<Poly_Noeud*>  grille;
  double x_min= list_points[0]->get_x(), x_max= list_points[6]->get_x();
  double y_min= list_points[0]->get_y(), y_max= list_points[6]->get_y();
  double z_min= list_points[0]->get_z(), z_max= list_points[6]->get_z();
  BOITE_3D boiteenglobante(x_min,y_min,z_min,x_max,y_max,z_max);
  boiteenglobante.change_grosseur(1.1);
  int nb=pow(list_points.size(),0.333333333333333333)+1;
  grille.initialiser(boiteenglobante.get_xmin(),boiteenglobante.get_ymin(),boiteenglobante.get_zmin(),boiteenglobante.get_xmax(),boiteenglobante.get_ymax(),boiteenglobante.get_zmax(),nb,nb,nb);
  
  double eps = epsfusion;
  
  
  for(int i=0; i < get_nb_cell(); i++)
  {
    Poly_Cellule* c = get_cell(i);
    for(int j=0; j < c->get_nb_noeud(); j++)
    {
      grille.inserer(c->get_noeud(j));
      lstnoeud.ajouter(c->get_noeud(j));
    }
  }
  
  
  TPL_MAP_ENTITE<Poly_Noeud*>::ITERATEUR it;
  
  for (Poly_Noeud* nd=lstnoeud.get_premier(it);nd!=NULL;nd=lstnoeud.get_suivant(it))
    if (nd->est_fusionne()==false)
        {
    double xyz[3];
    xyz[0]=nd->get_x();
        xyz[1]=nd->get_y();
        xyz[2]=nd->get_z();
        TPL_MAP_ENTITE<Poly_Noeud*> liste_entite_trouve;
        grille.rechercher(xyz[0],xyz[1],xyz[2],eps,liste_entite_trouve);
        TPL_MAP_ENTITE<Poly_Noeud*>::ITERATEUR it2;
        for (Poly_Noeud* nd2=liste_entite_trouve.get_premier(it2);nd2!=NULL;nd2=liste_entite_trouve.get_suivant(it2))
                if (nd2->est_fusionne()==false)
                  {
                  double xyz2[3];
          xyz2[0]=nd2->get_x();
                  xyz2[1]=nd2->get_y();
                  xyz2[2]=nd2->get_z();
                  OT_VECTEUR_3D vec(xyz,xyz2);
          if (vec.get_longueur()<eps)
                          {
                          nd2->set_x(xyz[0]);  
                          nd2->set_y(xyz[1]);  
                          nd2->set_z(xyz[2]);  
                          nd2->active_fusion();
                          nd2->change_maitre_fusion(nd);
              nd->active_fusion();
                          }
                    
                  }
        
        nd->change_maitre_fusion(nd);
        }
        
        
#ifdef TYTY
        
        
        for(int i=0; i < get_nb_cell(); i++)
        {
        Poly_Cellule* c = get_cell(i);
        for(int j=0; j < c->get_nb_face(); j++)
            {
           // std::cout << "face :" << j+1 << std::endl;
            Poly_Face *face=c->get_face(j);
            /*int numdep= face->get_noeud(0);
            Poly_Noeud *dep= c->get_noeud(numdep);
            int numfin= face->get_last_noeud();
            Poly_Noeud *fin= c->get_noeud(numfin);
            int fusionsuite=0;
            Poly_Noeud *lastfus;
            int numlastfus;
            int numnextfus;
            int numfus;
            Poly_Noeud *nextfus;
            Poly_Noeud *fus;
            std::vector<int> lstnewface;*/
            std::vector<int> aretes(c->get_nb_face());
            for (int i=0;i<face->get_nb_noeud();i++)
                {
                int num1= face->get_noeud(i);   
                Poly_Noeud* nd1=c->get_noeud(num1);
                Poly_Noeud* nd2;
                if (i+1!=face->get_nb_noeud()) 
                    {
                    int num2= face->get_noeud(i+1);   
                    nd2=c->get_noeud(num2);
                    }
                else 
                    {
                    int num2= face->get_noeud(0);   
                    nd2=c->get_noeud(num2);
                    }
                if (nd1->get_maitre_fusion()==nd2->get_maitre_fusion()) 
                    aretes[i]=2;
                else
                    {
                    double xyz1[3],xyz2[3];
                    xyz1[0]=nd1->get_x();
                    xyz1[1]=nd1->get_y();
                    xyz1[2]=nd1->get_z();
                    xyz2[0]=nd2->get_x();
                    xyz2[1]=nd2->get_y();
                    xyz2[2]=nd2->get_z();
                    OT_VECTEUR_3D vec(xyz1,xyz2);
                    if (vec.get_longueur()<epsfusion)
                        aretes[i]=0;    
                    else 
                    if (vec.get_longueur()<1e-10)
                        aretes[i]=2;
                    else
                        aretes[i]=1;    
                    }
                }
        int numlastfus;
        int numfus;
        int numnextfus;
        int fusionsuite=0;
        std::vector<int> lstnewface;
        if (aretes[0]==1) 
            numlastfus=1;
        else 
            {
            int num=face->get_nb_noeud()-1;
            while (aretes[num]==0) 
                num--;
            numlastfus=num+1;
            if (numlastfus>=face->get_nb_noeud()) numlastfus=0;
            fusionsuite=face->get_nb_noeud()-num;
            }
        for (int i=0;i<face->get_nb_noeud();i++)        
            {
            if (aretes[i]==0)
                {
                fusionsuite++; 
                numfus=i;
                }
            else
                {
                if (fusionsuite>1) 
                    {
                    numnextfus=i;
                    lstnewface.push_back(numlastfus);
                    lstnewface.push_back(numfus);
                    lstnewface.push_back(numnextfus);
                    }
                fusionsuite=0;
                numlastfus=i+1;
                }
                
            }
            
            
            
       
            
            for (int i=0;i<face->get_nb_noeud();i++)
                    if (aretes[i]==0)
                        face->supprimer_noeud(i+1);
 
            for (int i=0;i<lstnewface.size();i=i+3)
                {
                Poly_Face *face=new Poly_Face;
                face->add_noeud(lstnewface[i]);
                face->add_noeud(lstnewface[i+1]);
                face->add_noeud(lstnewface[i+2]+1);
                c->add_face(face);
                }            
            }
    
  }
#endif
}
 
 
 
void Poly_Voro::add_cell(Poly_Cellule* element)
{ 
  list_cellules.push_back(element); 
}
 
Poly_Point* Poly_Voro::get_point(int num)
{ 
  return list_points[num]; 
}
 
Poly_Cellule* Poly_Voro::get_cell(int num)
{ 
  return list_cellules[num]; 
}
 
int Poly_Voro::get_nb_cell(void)
{ 
  return list_cellules.size();
}