microstructure/src/mstruct_outils.cpp

#include "mstruct_outils.h"
#include "mg_cg_groupe_forme.h"
#include "ot_geometrie.h"
#include "mg_maillage.h"
#include "fem_maillage.h"
#include "mg_cg_forme.h"
#include "mg_cg_forme_multi_volume.h"
#include "mg_cg_forme_volume.h"
#include "ot_tenseur.h"
#include "mg_cg_info.h"
#include "fem_maillage_quadratique_outils.h"
#include <math.h>

using namespace MICROSTRUCTURE;

int OUTILS::statistiques_champ_volumique(FEM_SOLUTION* fem_sol,
                                         MG_CG_GROUPE_FORME* groupe_forme,
                                         BOITE_3D* boite_analyse,
                                         double* moyenne,
                                         double* ecart_type,
                                         double* min,
                                         double* max,
                                         OT_HISTOGRAMME** tab_histogramme)
{
  int nb_champ = fem_sol->get_nb_champ();
  double volume_total=0.0;
  for(int i=0;i<nb_champ;i++)
  {
    moyenne[i]=0.0;
    ecart_type[i]=0.0;
    min[i]=numeric_limits< double >::max();
    max[i]=numeric_limits< double >::min();
  }
  std::vector<double*> vector_champ;
  std::vector<double> vector_volume;
  LISTE_FEM_ELEMENT3::iterator it;
  FEM_MAILLAGE* fem_mail = fem_sol->get_maillage();
  long k=0;
  for (FEM_ELEMENT3*  ele=fem_mail->get_premier_element3(it);ele!=NULL;ele=fem_mail->get_suivant_element3(it))
  {
    bool ok=false;
    if(groupe_forme!=NULL)
    {
      TPL_MAP_ENTITE<MG_VOLUME*> tpl_map_volume = groupe_forme->get_tpl_map_volume();
      TPL_MAP_ENTITE<MG_VOLUME*>::ITERATEUR it_volume;
      MG_VOLUME* mg_vol=tpl_map_volume.get_premier(it_volume);
      while(ok==false && mg_vol!=NULL)
      {
        if(ele->get_lien_topologie()==mg_vol) ok=true;
        else mg_vol=tpl_map_volume.get_suivant(it_volume);
      }
    } else ok=true;
    if(boite_analyse!=NULL && ok==true)
    {
      BOITE_3D b3d = ele->get_boite_3D();
      if(!boite_analyse->contient(b3d.get_xcentre(),b3d.get_ycentre(),b3d.get_zcentre())) ok=false;
    }
    if(ok==true)
    {
      double *champ=new double[nb_champ];
      double volume;
      for(int i=0;i<nb_champ;i++)
      {
        double unite=fem_mail->get_mg_geometrie()->get_valeur_unite();
        double val=0.;
        volume=0;
        int nbgauss=fem_sol->get_nb_gauss(k);
        for (int r=0;r<nbgauss;r++)
        {
          int degre=ele->get_degre_gauss(nbgauss);
          double w;
          double uvw[3];
          ele->get_pt_gauss(degre,r,w,uvw);
          double jac[9];
          int l,c;
          double det=ele->get_jacobien(jac,uvw,l,c,unite);
          double valeur=fem_sol->lire(k,i,0,r)*fabs(det);
          val=val+w*valeur;
          volume=volume+w*fabs(det);
        }
        champ[i] = val;
      }
      vector_champ.push_back(champ);
      vector_volume.push_back(volume);
      volume_total+=volume;
    }
    k++;
  }  
  long nb_echantillions=vector_champ.size();
  std::vector<double*>::iterator it_vector_champ;
  std::vector<double>::iterator it_vector_volume;
  for(it_vector_champ=vector_champ.begin(),it_vector_volume=vector_volume.begin();it_vector_champ!=vector_champ.end();it_vector_champ++,it_vector_volume++)
  {
    double *champ = *it_vector_champ;
    double volume = *it_vector_volume;
/*    for(int i=0;i<nb_champ;i++)
    {
      double champ_volumique = champ[i]/volume;
      moyenne[i]+=champ_volumique;
      if(champ_volumique<min[i]) min[i]=champ_volumique;
      if(champ_volumique>max[i]) max[i]=champ_volumique;
      tab_histogramme[i]->ajouter_valeur(champ[i],volume/volume_total);
    }*/
    for(int i=0;i<nb_champ;i++)
    {
      moyenne[i]+=champ[i];
      if(champ[i]<min[i]) min[i]=champ[i];
      if(champ[i]>max[i]) max[i]=champ[i];
      tab_histogramme[i]->ajouter_valeur(champ[i],volume/volume_total);
    }
  }
  for(int i=0;i<nb_champ;i++)
  {
//     moyenne[i]=moyenne[i]/nb_echantillions;
    moyenne[i]=moyenne[i]/volume_total;
  }
  for(it_vector_champ=vector_champ.begin(),it_vector_volume=vector_volume.begin();it_vector_champ!=vector_champ.end();it_vector_champ++,it_vector_volume++)
  {
    double *champ = *it_vector_champ;
    double volume = *it_vector_volume;
    for(int i=0;i<nb_champ;i++)
    {
//       double champ_volumique = champ[i]/volume;
//       ecart_type[i]+=(champ_volumique-moyenne[i])*(champ_volumique-moyenne[i]);
      ecart_type[i]+=(champ[i]-moyenne[i])*(champ[i]-moyenne[i]);
    }
  }
  for(it_vector_champ=vector_champ.begin();it_vector_champ!=vector_champ.end();it_vector_champ++) delete [] *it_vector_champ;
  for(int i=0;i<nb_champ;i++)
  {
    ecart_type[i]=ecart_type[i]*(1.0/(nb_echantillions-1.0));
    ecart_type[i]=sqrt(ecart_type[i]);
  }
  return OK;  
}

int OUTILS::statistiques_tenseur_orientation(MG_CG_GROUPE_FORME* groupe_forme,
                                             BOITE_3D* boite_analyse,
                                             double* moyenne,
                                             double* ecart_type,
                                             double* min,
                                             double* max,
                                             bool ponderation_volume)
{
  double volume_total=0;
  std::vector<double> vector_volume;
  for(int i=0;i<6;i++)
  {
    moyenne[i]=0.0;
    ecart_type[i]=0.0;
    min[i]=numeric_limits< double >::max();
    max[i]=numeric_limits< double >::min();
  }
  std::vector<double*> vector_tenseur;
  std::map<long,MG_CG_FORME*>::iterator it_forme;
  for(MG_CG_FORME* forme=groupe_forme->get_premiere_mgcg_forme(it_forme);forme!=NULL;forme=groupe_forme->get_suivante_mgcg_forme(it_forme))
  {
    BOITE_3D boite_forme = forme->get_boite_3D();
    bool ok=true;
    if(boite_analyse!=NULL)
    {
      if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymin(),boite_forme.get_zmin())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymin(),boite_forme.get_zmin())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymax(),boite_forme.get_zmin())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymax(),boite_forme.get_zmin())) ok=true;    
      else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymin(),boite_forme.get_zmax())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymin(),boite_forme.get_zmax())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmin(),boite_forme.get_ymax(),boite_forme.get_zmax())) ok=true;
      else if(boite_analyse->contient(boite_forme.get_xmax(),boite_forme.get_ymax(),boite_forme.get_zmax())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymin(),boite_analyse->get_zmin())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymin(),boite_analyse->get_zmin())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymax(),boite_analyse->get_zmin())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymax(),boite_analyse->get_zmin())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymin(),boite_analyse->get_zmax())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymin(),boite_analyse->get_zmax())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmin(),boite_analyse->get_ymax(),boite_analyse->get_zmax())) ok=true;
      else if(boite_forme.contient(boite_analyse->get_xmax(),boite_analyse->get_ymax(),boite_analyse->get_zmax())) ok=true;
      else ok=false;
    }
    if(ok)
    {
      double axe_forme[3];
      MG_CG_INFO_VCT_DOUBLE* info_axe = (MG_CG_INFO_VCT_DOUBLE*)forme->get_mgcg_info("AXE");
      std::vector<double> vct_double = info_axe->get_vct_valeur();
      axe_forme[0]=vct_double[0];
      axe_forme[1]=vct_double[1];
      axe_forme[2]=vct_double[2];
      double *tenseur = new double[6];
      tenseur[0] = axe_forme[0]*axe_forme[0];
      tenseur[1] = axe_forme[1]*axe_forme[1];
      tenseur[2] = axe_forme[2]*axe_forme[2];
      tenseur[3] = axe_forme[0]*axe_forme[1];
      tenseur[4] = axe_forme[0]*axe_forme[2];
      tenseur[5] = axe_forme[1]*axe_forme[2];
      if(ponderation_volume)
      {
        if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
        {
          MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
          double volume=OT_GEOMETRIE::get_volume(forme_volume->get_mg_volume(),1.0e-06);
          vector_volume.push_back(volume);
          volume_total+=volume;
        }
        else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
        {
          MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
          std::map<long,MG_VOLUME*>::iterator it_volume;
          double volume=0;
          for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
          {
            volume+=OT_GEOMETRIE::get_volume(vol,1.0e-06);
          }
          vector_volume.push_back(volume);
          volume_total+=volume;
        }
      }
      vector_tenseur.push_back(tenseur);
    }
  }
  long nb_forme=vector_tenseur.size();
  std::vector<double*>::iterator it_tenseur;
  std::vector<double>::iterator it_volume=vector_volume.begin();
  for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
  {
    double* tenseur = *it_tenseur;
    for(int i=0;i<6;i++)
    {
      if(!ponderation_volume)
      {
        moyenne[i]+=tenseur[i];
        if(tenseur[i]<min[i]) min[i]=tenseur[i];
        if(tenseur[i]>max[i]) max[i]=tenseur[i];
      }
      else
      {
        double volume=*it_volume;
        moyenne[i]+=tenseur[i]*volume/volume_total;
        if(tenseur[i]<min[i]) min[i]=tenseur[i];
        if(tenseur[i]>max[i]) max[i]=tenseur[i];
      }
    }
  }
  if(!ponderation_volume)
  for(int i=0;i<6;i++)
  {
    moyenne[i]=moyenne[i]/nb_forme;
  }
  for(it_tenseur=vector_tenseur.begin();it_tenseur!=vector_tenseur.end();it_tenseur++)
  {
    double* tenseur = *it_tenseur;
    for(int i=0;i<6;i++)
    {
      ecart_type[i]+=(tenseur[i]-moyenne[i])*(tenseur[i]-moyenne[i]);
    }
    delete [] *it_tenseur;
  }
  for(int i=0;i<6;i++)
  {
    ecart_type[i]=ecart_type[i]*(1.0/(nb_forme-1.0));
    ecart_type[i]=sqrt(ecart_type[i]);
  }
  return OK;
}

int OUTILS::statistiques_CAO(BOITE_3D boite_VES,
                             TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
                             long int& nb_forme,
                             long int& nb_volume,
                             double& volume_forme_moyenne,
                             double& volume_forme_ecart_type,
                             double& volume_forme_min,
                             double& volume_forme_max,
                             OT_HISTOGRAMME& volume_forme_histogramme,
                             double& volume_total,
                             double& fraction_volumique,
                             double precision)
{
  double volume_ves = boite_VES.get_volume(); 
  nb_forme=0;
  nb_volume=0;
  volume_forme_moyenne=0;
  volume_forme_ecart_type=0;
  volume_forme_min=numeric_limits< double >::max();
  volume_forme_max=numeric_limits< double >::min();
  volume_total=0;  
  std::vector<double> vector_volume;
  
  TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
  for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
  {
    TPL_MAP_ENTITE<MG_VOLUME*> tpl_map_volume;
    if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
    {
      MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
      tpl_map_volume.ajouter(forme_volume->get_mg_volume());
    }
    else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
    {
      MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
      std::map<long,MG_VOLUME*>::iterator it_volume;
      for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
      {
        tpl_map_volume.ajouter(vol);
      }      
    }
    double volume_forme=0;
    nb_volume += tpl_map_volume.get_nb();
    TPL_MAP_ENTITE<MG_VOLUME*>::ITERATEUR it_mg_volume;
    for(MG_VOLUME* vol=tpl_map_volume.get_premier(it_mg_volume);vol!=NULL;vol=tpl_map_volume.get_suivant(it_mg_volume))
    {
      volume_forme += OT_GEOMETRIE::get_volume(vol,precision);
    }
    vector_volume.push_back(volume_forme);   
  }
  std::vector<double>::iterator it_data;
  for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
  {
    double volume = *it_data;
    volume_total+=volume;
    if(volume<volume_forme_min) volume_forme_min=volume;
    if(volume>volume_forme_max) volume_forme_max=volume;
  }
  for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
  {
    double volume = *it_data;
    volume_forme_histogramme.ajouter_valeur(volume,volume/volume_total);
  }
  nb_forme=vector_volume.size();   
  volume_forme_moyenne=volume_total/nb_forme;
  if(nb_forme>1)
  {
  for(it_data=vector_volume.begin();it_data!=vector_volume.end();it_data++)
  {
    double volume = *it_data;
    volume_forme_ecart_type+=(volume-volume_forme_moyenne)*(volume-volume_forme_moyenne);
  }
  volume_forme_ecart_type=volume_forme_ecart_type*(1.0/(nb_forme-1.0));
  volume_forme_ecart_type=sqrt(volume_forme_ecart_type);
  }
  else volume_forme_ecart_type = nan("");
  fraction_volumique = volume_total/volume_ves;  
}

int OUTILS::statistiques_mg_maillage(BOITE_3D boite_VES,
                                     MG_MAILLAGE* mai,
                                     TPL_MAP_ENTITE<MG_CG_FORME*> &tpl_map_forme,
                                     long &nb_element_2D,long &nb_element_3D,
                                     double& qualite_moyenne_2D, double& qualite_ecart_type_2D, double& qualite_min_2D, double& qualite_max_2D, OT_HISTOGRAMME& histogramme_qualite_2D,
                                     double& qualite_moyenne_3D, double& qualite_ecart_type_3D, double& qualite_min_3D, double& qualite_max_3D, OT_HISTOGRAMME& histogramme_qualite_3D,
                                     double& taille_moyenne_2D, double& taille_ecart_type_2D, double& taille_min_2D, double& taille_max_2D, OT_HISTOGRAMME& histogramme_taille_2D,
                                     double& taille_moyenne_3D, double& taille_ecart_type_3D, double& taille_min_3D, double& taille_max_3D, OT_HISTOGRAMME& histogramme_taille_3D,
                                     double& volume, double& fraction_volumique)
{
  double volume_ves = boite_VES.get_volume();
  nb_element_2D=0;nb_element_3D=0;
  qualite_moyenne_2D=0;qualite_ecart_type_2D=0;qualite_min_2D=numeric_limits< double >::max();;qualite_max_2D=numeric_limits< double >::min();
  qualite_moyenne_3D=0;qualite_ecart_type_3D=0;qualite_min_3D=numeric_limits< double >::max();;qualite_max_3D=numeric_limits< double >::min();
  taille_moyenne_2D=0;taille_ecart_type_2D=0;taille_min_2D=numeric_limits< double >::max();;taille_max_2D=numeric_limits< double >::min();
  taille_moyenne_3D=0;taille_ecart_type_3D=0;taille_min_3D=numeric_limits< double >::max();;taille_max_3D=numeric_limits< double >::min();
  volume=0;fraction_volumique=0; 
  TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_3D;  
  TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> tpl_map_element_maill_2D;  
  TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
  for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
  {
    if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
    {
      MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
      TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_maillage();
      TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
      for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_3D.ajouter(ele);
      TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
      OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
      TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
      for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
      {
        tmp_maill = face->get_lien_maillage();
        for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_2D.ajouter(ele);
      }
      double volume_i;
      OT_VECTEUR_3D cdm;
      OT_MATRICE_3D inertieglobale;
      OT_MATRICE_3D inertiecdm;
      forme_volume->get_mg_volume()->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
      volume+=volume_i;
    }
    else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
    {
      MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
      std::map<long,MG_VOLUME*>::iterator it_volume;
      for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
      {
        TPL_SET<MG_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_maillage();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it_maill;
        for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_3D.ajouter(ele);
        TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
        OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
        TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
        for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
        {
          tmp_maill = face->get_lien_maillage();
          for(MG_ELEMENT_MAILLAGE* ele=tmp_maill->get_premier(it_maill);ele!=NULL;ele=tmp_maill->get_suivant(it_maill)) tpl_map_element_maill_2D.ajouter(ele);
        }
        double volume_i;
        OT_VECTEUR_3D cdm;
        OT_MATRICE_3D inertieglobale;
        OT_MATRICE_3D inertiecdm;
        vol->get_propriete_massique(mai,volume_i,cdm,inertieglobale,inertiecdm);
        volume+=volume_i;
      }      
    }
  }
  fraction_volumique = volume/volume_ves;
  std::vector<double> vector_qualite_3D;
  std::vector<double> vector_taille_3D;
  LISTE_MG_TETRA::iterator it_tetra;
  for(MG_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
  {
    if(tpl_map_element_maill_3D.existe(tet))
    {
      double qual=OPERATEUR::qualite_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
      vector_qualite_3D.push_back(qual);
      if(qual<qualite_min_3D) qualite_min_3D=qual;
      if(qual>qualite_max_3D) qualite_max_3D=qual;
      qualite_moyenne_3D+=qual;
      double taille=OPERATEUR::taille_tetra(tet->get_noeud1()->get_coord(),tet->get_noeud2()->get_coord(),tet->get_noeud3()->get_coord(),tet->get_noeud4()->get_coord());
      vector_taille_3D.push_back(taille);
      if(taille<taille_min_3D) taille_min_3D=taille;
      if(taille>taille_max_3D) taille_max_3D=taille;
      taille_moyenne_3D+=taille;
      nb_element_3D++;
    }
  }    
  qualite_moyenne_3D = qualite_moyenne_3D/nb_element_3D;
  taille_moyenne_3D = taille_moyenne_3D/nb_element_3D;
  std::vector<double>::iterator it_qual_3D;
  
  for(it_qual_3D=vector_qualite_3D.begin();it_qual_3D!=vector_qualite_3D.end();it_qual_3D++)
  {
    double qual = *it_qual_3D;
    qualite_ecart_type_3D+=(qual-qualite_moyenne_3D)*(qual-qualite_moyenne_3D);
    histogramme_qualite_3D.ajouter_valeur(qual,1.0/nb_element_3D);
  }
  qualite_ecart_type_3D=qualite_ecart_type_3D*(1.0/(nb_element_3D-1.0));
  qualite_ecart_type_3D=sqrt(qualite_ecart_type_3D);
  
  std::vector<double>::iterator it_taille_3D;
  for(it_taille_3D=vector_taille_3D.begin();it_taille_3D!=vector_taille_3D.end();it_taille_3D++)
  {
    double taille = *it_taille_3D;
    taille_ecart_type_3D=(taille-taille_moyenne_3D)*(taille-taille_moyenne_3D);
    histogramme_taille_3D.ajouter_valeur(taille,1.0/nb_element_3D);
  }
  taille_ecart_type_3D=taille_ecart_type_3D*(1.0/(nb_element_3D-1.0));
  taille_ecart_type_3D=sqrt(taille_ecart_type_3D);
  
  
  std::vector<double> vector_qualite_2D;
  std::vector<double> vector_taille_2D;
  LISTE_MG_TRIANGLE::iterator it_tri;
  for(MG_TRIANGLE* tri=mai->get_premier_triangle(it_tri);tri!=NULL;tri=mai->get_suivant_triangle(it_tri))
  {
    if(tpl_map_element_maill_2D.existe(tri))
    {
      double qual=OPERATEUR::qualite_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
      vector_qualite_2D.push_back(qual);
      if(qual<qualite_min_2D) qualite_min_2D=qual;
      if(qual>qualite_max_2D) qualite_max_2D=qual;
      qualite_moyenne_2D+=qual;
      double taille=OPERATEUR::taille_triangle(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord(),tri->get_noeud3()->get_coord());
      vector_taille_2D.push_back(taille);
      if(taille<taille_min_2D) taille_min_2D=taille;
      if(taille>taille_max_2D) taille_max_2D=taille;
      taille_moyenne_2D+=taille;
      nb_element_2D++;
    }
  }    
  qualite_moyenne_2D = qualite_moyenne_2D/nb_element_2D;
  taille_moyenne_2D = taille_moyenne_2D/nb_element_2D;
  std::vector<double>::iterator it_qual_2D;
  
  for(it_qual_2D=vector_qualite_2D.begin();it_qual_2D!=vector_qualite_2D.end();it_qual_2D++)
  {
    double qual = *it_qual_2D;
    qualite_ecart_type_2D+=(qual-qualite_moyenne_2D)*(qual-qualite_moyenne_2D);
    histogramme_qualite_2D.ajouter_valeur(qual,1.0/nb_element_2D);
  }
  qualite_ecart_type_2D=qualite_ecart_type_2D*(1.0/(nb_element_2D-1.0));
  qualite_ecart_type_2D=sqrt(qualite_ecart_type_2D);
  
  std::vector<double>::iterator it_taille_2D;
  for(it_taille_2D=vector_taille_2D.begin();it_taille_2D!=vector_taille_2D.end();it_taille_2D++)
  {
    double taille = *it_taille_2D;
    taille_ecart_type_2D=(taille-taille_moyenne_2D)*(taille-taille_moyenne_2D);
    histogramme_taille_2D.ajouter_valeur(taille,1.0/nb_element_2D);
  }
  taille_ecart_type_2D=taille_ecart_type_2D*(1.0/(nb_element_2D-1.0));
  taille_ecart_type_2D=sqrt(taille_ecart_type_2D);
  
}


int OUTILS::statistiques_fem_maillage(FEM_MAILLAGE* fem,
                                      BOITE_3D *boite_analyse,
                                      TPL_MAP_ENTITE< MG_CG_FORME* >& tpl_map_forme,
                                      long int& nb_element_2D, long int& nb_element_3D,
                                      double& jacobien_min_moyenne_2D, double& jacobien_min_ecart_type_2D, double& jacobien_min_min_2D, double& jacobien_min_max_2D, OT_HISTOGRAMME& histogramme_jacobien_min_2D,
                                      double& jacobien_max_moyenne_2D, double& jacobien_max_ecart_type_2D, double& jacobien_max_min_2D, double& jacobien_max_max_2D, OT_HISTOGRAMME& histogramme_jacobien_max_2D,
                                      double& jacobien_min_moyenne_3D, double& jacobien_min_ecart_type_3D, double& jacobien_min_min_3D, double& jacobien_min_max_3D, OT_HISTOGRAMME& histogramme_jacobien_min_3D,
                                      double& jacobien_max_moyenne_3D, double& jacobien_max_ecart_type_3D, double& jacobien_max_min_3D, double& jacobien_max_max_3D, OT_HISTOGRAMME& histogramme_jacobien_max_3D,
                                      double& distortion_moyenne_2D, double& distortion_ecart_type_2D, double& distortion_min_2D, double& distortion_max_2D, OT_HISTOGRAMME& histogramme_distortion_2D,
                                      double& distortion_moyenne_3D, double& distortion_ecart_type_3D, double& distortion_min_3D, double& distortion_max_3D, OT_HISTOGRAMME& histogramme_distortion_3D,
                                      double& volume,
                                      double& fraction_volumique)
{
  double volume_ves=0;
  nb_element_2D=0;nb_element_3D=0;
  jacobien_min_moyenne_2D=0.0; jacobien_min_ecart_type_2D=0.0; jacobien_min_min_2D=numeric_limits< double >::max(); jacobien_min_max_2D=numeric_limits< double >::min();
  jacobien_max_moyenne_2D=0.0; jacobien_max_ecart_type_2D=0.0; jacobien_max_min_2D=numeric_limits< double >::max(); jacobien_max_max_2D=numeric_limits< double >::min();
  jacobien_min_moyenne_3D=0.0; jacobien_min_ecart_type_3D=0.0; jacobien_min_min_3D=numeric_limits< double >::max(); jacobien_min_max_3D=numeric_limits< double >::min();
  jacobien_max_moyenne_3D=0.0; jacobien_max_ecart_type_3D=0.0; jacobien_max_min_3D=numeric_limits< double >::max(); jacobien_max_max_3D=numeric_limits< double >::min();
  distortion_moyenne_2D=0.0; distortion_ecart_type_2D=0.0; distortion_min_2D=numeric_limits< double >::max(); distortion_max_2D=numeric_limits< double >::min();
  distortion_moyenne_3D=0.0; distortion_ecart_type_3D=0.0; distortion_min_3D=numeric_limits< double >::max(); distortion_max_3D=numeric_limits< double >::min();
  volume=0;
  fraction_volumique=0;
  
  LISTE_FEM_ELEMENT3::iterator it_element3;
  for(FEM_ELEMENT3* ele3=fem->get_premier_element3(it_element3);ele3!=NULL;ele3=fem->get_suivant_element3(it_element3))
  {
    BOITE_3D b3d_ele = ele3->get_boite_3D();
    if(boite_analyse!=NULL)
          if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
    int nb_pt_gauss=ele3->get_nb_pt_gauss(ele3->get_degremax_fonction_interpolation());
    for (int r=0;r<nb_pt_gauss;r++)
    {
      int degre=ele3->get_degre_gauss(nb_pt_gauss);
      double w;
      double uvw[3];
      ele3->get_pt_gauss(degre,r,w,uvw);
      double jac[9];
      int l,c;
      double det=ele3->get_jacobien(jac,uvw,l,c);
      volume_ves+=w*det;
    }
  }  
  TPL_MAP_ENTITE<FEM_ELEMENT3*> tpl_map_element_maill_3D;  
  TPL_MAP_ENTITE<FEM_ELEMENT2*> tpl_map_element_maill_2D;  
  TPL_MAP_ENTITE<MG_CG_FORME*>::ITERATEUR it_forme;
  for(MG_CG_FORME* forme=tpl_map_forme.get_premier(it_forme);forme!=NULL;forme=tpl_map_forme.get_suivant(it_forme))
  {
    if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::VOLUME)
    {
      MG_CG_FORME_VOLUME* forme_volume = (MG_CG_FORME_VOLUME*)forme;
      TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = forme_volume->get_mg_volume()->get_lien_fem_maillage();
      long tmp_nb = tmp_maill->get_nb();
      for(long i=0;i<tmp_nb;i++)
      {
        FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
        BOITE_3D b3d_ele = ele3->get_boite_3D();
        if(boite_analyse!=NULL)
          if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
        tpl_map_element_maill_3D.ajouter(ele3);
      }
      TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
      OT_GEOMETRIE::get_map_mg_face_sous_jacent(forme_volume->get_mg_volume(),tpl_map_face);
      TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
      for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
      {
        tmp_maill = face->get_lien_fem_maillage();
        tmp_nb = tmp_maill->get_nb();   
        for(long i=0;i<tmp_nb;i++)
        {
          FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
          BOITE_3D b3d_ele = ele2->get_boite_3D();
          if(boite_analyse!=NULL)
            if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
          tpl_map_element_maill_2D.ajouter(ele2);
        }
      }
    }
    else if(forme->get_type_forme()==MG_CG_FORME::TYPE_FORME::MULTI_VOLUME)
    {
      MG_CG_FORME_MULTI_VOLUME* forme_multi_volume = (MG_CG_FORME_MULTI_VOLUME*)forme;
      std::map<long,MG_VOLUME*>::iterator it_volume;
      for(MG_VOLUME* vol=forme_multi_volume->get_premier_mg_volume(it_volume);vol!=NULL;vol=forme_multi_volume->get_suivant_mg_volume(it_volume))
      {
        TPL_LISTE_ENTITE<FEM_ELEMENT_MAILLAGE*> *tmp_maill = vol->get_lien_fem_maillage();
        long tmp_nb = tmp_maill->get_nb();
        for(long i=0;i<tmp_nb;i++)
        {
          FEM_ELEMENT3* ele3 = (FEM_ELEMENT3*)tmp_maill->get(i);
          BOITE_3D b3d_ele = ele3->get_boite_3D();
          if(boite_analyse!=NULL)
            if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
          tpl_map_element_maill_3D.ajouter(ele3);
        }
        TPL_MAP_ENTITE<MG_FACE*> tpl_map_face;
        OT_GEOMETRIE::get_map_mg_face_sous_jacent(vol,tpl_map_face);
        TPL_MAP_ENTITE<MG_FACE*>::ITERATEUR it_face;
        for(MG_FACE* face=tpl_map_face.get_premier(it_face);face!=NULL;face=tpl_map_face.get_suivant(it_face))
        {
          tmp_maill = face->get_lien_fem_maillage();
          tmp_nb = tmp_maill->get_nb();
          for(long i=0;i<tmp_nb;i++)
          {
            FEM_ELEMENT2* ele2 = (FEM_ELEMENT2*)tmp_maill->get(i);
            BOITE_3D b3d_ele = ele2->get_boite_3D();
            if(boite_analyse!=NULL)
              if(!boite_analyse->contient(b3d_ele.get_xcentre(),b3d_ele.get_ycentre(),b3d_ele.get_zcentre())) continue;
            tpl_map_element_maill_2D.ajouter(ele2);
          }
        }
      }      
    }
  }
  nb_element_3D=tpl_map_element_maill_3D.get_nb();
  nb_element_2D=tpl_map_element_maill_2D.get_nb();
  std::vector<double> vector_jacobien_min_2D;
  std::vector<double> vector_jacobien_max_2D;
  std::vector<double> vector_jacobien_min_3D;
  std::vector<double> vector_jacobien_max_3D;
  std::vector<double> vector_distortion_2D;
  std::vector<double> vector_distortion_3D;
  FEM_MAILLAGE_QUADRATIQUE_OUTILS ot;
  TPL_MAP_ENTITE<FEM_ELEMENT3*>::ITERATEUR it_ele3;
  for(FEM_ELEMENT3* ele3=tpl_map_element_maill_3D.get_premier(it_ele3);ele3!=NULL;ele3=tpl_map_element_maill_3D.get_suivant(it_ele3))
  {
    int nb_pt_gauss=ele3->get_nb_pt_gauss(ele3->get_degremax_fonction_interpolation());
    for (int r=0;r<nb_pt_gauss;r++)
    {
      int degre=ele3->get_degre_gauss(nb_pt_gauss);
      double w;
      double uvw[3];
      ele3->get_pt_gauss(degre,r,w,uvw);
      double jac[9];
      int l,c;
      double det=ele3->get_jacobien(jac,uvw,l,c);
      volume+=w*det;
    }
    double jmin=ot.get_jmin(ele3);
    double jmax=ot.get_jmin(ele3);
    double distortion = ot.get_distorsion2(ele3);
    vector_jacobien_min_3D.push_back(jmin);
    vector_jacobien_max_3D.push_back(jmax);
    vector_distortion_3D.push_back(distortion);
    jacobien_min_moyenne_3D+=jmin;
    jacobien_max_moyenne_3D+=jmax;
    distortion_moyenne_3D+=distortion;
  }
  fraction_volumique=volume/volume_ves;
  jacobien_min_moyenne_3D=jacobien_min_moyenne_3D/nb_element_3D;
  jacobien_max_moyenne_3D=jacobien_max_moyenne_3D/nb_element_3D;
  distortion_moyenne_3D=distortion_moyenne_3D/nb_element_3D;    
  std::vector<double>::iterator it_jac_min_3D=vector_jacobien_min_3D.begin();
  std::vector<double>::iterator it_jac_max_3D=vector_jacobien_max_3D.begin();
  std::vector<double>::iterator it_dist_3D=vector_distortion_3D.begin();
  for(long i=0;i<nb_element_3D;i++)
  {
    double jmin=*it_jac_min_3D;
    double jmax=*it_jac_max_3D;
    double dist=*it_dist_3D;
    if(jmin<jacobien_min_min_3D)jacobien_min_min_3D=jmin;
    if(jmin>jacobien_min_max_3D)jacobien_min_max_3D=jmin;    
    if(jmax<jacobien_max_min_3D)jacobien_max_min_3D=jmax;
    if(jmax>jacobien_max_max_3D)jacobien_max_max_3D=jmax;
    if(dist<distortion_min_3D)distortion_min_3D=dist;
    if(dist>distortion_max_3D)distortion_max_3D=dist;    
    jacobien_min_ecart_type_3D=(jmin-jacobien_min_moyenne_3D)*(jmin-jacobien_min_moyenne_3D);
    jacobien_max_ecart_type_3D=(jmax-jacobien_max_moyenne_3D)*(jmax-jacobien_max_moyenne_3D);
    distortion_ecart_type_3D=(dist-distortion_moyenne_3D)*(dist-distortion_moyenne_3D);
    histogramme_jacobien_min_3D.ajouter_valeur(jmin,1.0/nb_element_3D);
    histogramme_jacobien_max_3D.ajouter_valeur(jmax,1.0/nb_element_3D);
    histogramme_distortion_3D.ajouter_valeur(dist,1.0/nb_element_3D);
    it_jac_min_3D++;
    it_jac_max_3D++;
    it_dist_3D++;
  }
  jacobien_min_ecart_type_3D=sqrt(jacobien_min_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
  jacobien_max_ecart_type_3D=sqrt(jacobien_max_ecart_type_3D*(1.0/(nb_element_3D-1.0)));
  distortion_ecart_type_3D=sqrt(distortion_ecart_type_3D*(1.0/(nb_element_3D-1.0)));  
  
  TPL_MAP_ENTITE<FEM_ELEMENT2*>::ITERATEUR it_ele2;
  for(FEM_ELEMENT2* ele2=tpl_map_element_maill_2D.get_premier(it_ele2);ele2!=NULL;ele2=tpl_map_element_maill_2D.get_suivant(it_ele2))
  {
    double jmin=ot.get_jmin(ele2);
    double jmax=ot.get_jmin(ele2);
    double distortion = ot.get_distorsion2(ele2);
    vector_jacobien_min_2D.push_back(jmin);
    vector_jacobien_max_2D.push_back(jmax);
    vector_distortion_2D.push_back(distortion);
    jacobien_min_moyenne_2D+=jmin;
    jacobien_max_moyenne_2D+=jmax;
    distortion_moyenne_2D+=distortion;
  }
  jacobien_min_moyenne_2D=jacobien_min_moyenne_2D/nb_element_2D;
  jacobien_max_moyenne_2D=jacobien_max_moyenne_2D/nb_element_2D;
  distortion_moyenne_2D=distortion_moyenne_2D/nb_element_2D;    
  std::vector<double>::iterator it_jac_min_2D=vector_jacobien_min_2D.begin();
  std::vector<double>::iterator it_jac_max_2D=vector_jacobien_max_2D.begin();
  std::vector<double>::iterator it_dist_2D=vector_distortion_2D.begin();
  for(long i=0;i<nb_element_2D;i++)
  {
    double jmin=*it_jac_min_2D;
    double jmax=*it_jac_max_2D;
    double dist=*it_dist_2D;
    if(jmin<jacobien_min_min_2D)jacobien_min_min_2D=jmin;
    if(jmin>jacobien_min_max_2D)jacobien_min_max_2D=jmin;    
    if(jmax<jacobien_max_min_2D)jacobien_max_min_2D=jmax;
    if(jmax>jacobien_max_max_2D)jacobien_max_max_2D=jmax;
    if(dist<distortion_min_2D)distortion_min_2D=dist;
    if(dist>distortion_max_2D)distortion_max_2D=dist;    
    jacobien_min_ecart_type_2D=(jmin-jacobien_min_moyenne_2D)*(jmin-jacobien_min_moyenne_2D);
    jacobien_max_ecart_type_2D=(jmax-jacobien_max_moyenne_2D)*(jmax-jacobien_max_moyenne_2D);
    distortion_ecart_type_2D=(dist-distortion_moyenne_2D)*(dist-distortion_moyenne_2D);
    histogramme_jacobien_min_2D.ajouter_valeur(jmin,1.0/nb_element_2D);
    histogramme_jacobien_max_2D.ajouter_valeur(jmax,1.0/nb_element_2D);
    histogramme_distortion_2D.ajouter_valeur(dist,1.0/nb_element_2D);
    it_jac_min_2D++;
    it_jac_max_2D++;
    it_dist_2D++;
  }
  jacobien_min_ecart_type_2D=sqrt(jacobien_min_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
  jacobien_max_ecart_type_2D=sqrt(jacobien_max_ecart_type_2D*(1.0/(nb_element_2D-1.0)));
  distortion_ecart_type_2D=sqrt(distortion_ecart_type_2D*(1.0/(nb_element_2D-1.0)));  
}


Revision:	943
Committed:	Fri Jun 8 18:47:16 2018 UTC (6 years, 11 months ago) by couturad
File size:	32902 byte(s)
Log Message:	Correction bugs dans MICROSTRUCTURE