criaqoperation/src/criaqoperators.cpp

#include <stdio.h>
#include "mg_file.h"
#include "criaqoperators.h"
#include "tpl_grille.h"
#include "ot_mathematique.h"
#include <math.h>
#include "tpl_octree.h"
#include "fem_maillage_outils.h"
#include "mg_maillage_outils.h"
#include <sstream>
#include <string.h>
#include <stdlib.h> 

CRIAQOPERATORS::CRIAQOPERATORS():MAILLEUR(false)
{
}
CRIAQOPERATORS::~CRIAQOPERATORS()
{
}

void CRIAQOPERATORS::bumparea_calcul(char* magicfilename,int nummai,int mgsolid,double tolerance,char* areamagicfilename,int defnodid)
{
  MG_FILE gest(magicfilename);
  MG_MAILLAGE* mai;
  if (nummai==0) mai=gest.get_mg_maillage(nummai); else mai=gest.get_mg_maillageid(nummai);
  
  MG_SOLUTION* mgsol=gest.get_mg_solutionid(mgsolid);
  mgsol->active_solution(0);
  MG_NOEUD* defnd;
  LISTE_MG_NOEUD::iterator itnd;
  for(MG_NOEUD* nd=mai->get_premier_noeud(itnd);nd!=NULL;nd=mai->get_suivant_noeud(itnd))
  {
    if (nd->get_id()==defnodid)  
      defnd=nd;
  }
  //cout<<defnd->get_solution()<<endl;
  double max_dist=defnd->get_solution();
  TPL_MAP_ENTITE<MG_TRIANGLE*> areatris;
  for (int i=0;i<defnd->get_lien_triangle()->get_nb();i++)
  {
    MG_TRIANGLE* trii=defnd->get_lien_triangle()->get(i);
    double avdis=(trii->get_noeud1()->get_solution()+trii->get_noeud2()->get_solution()+trii->get_noeud3()->get_solution())/3.;
      if (avdis>tolerance)
        areatris.ajouter(trii);
  }
 int stoplvl=0; 
 while (areatris.get_nb()>stoplvl)
 {
   stoplvl=areatris.get_nb();
  TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR itareatris;
  for(MG_TRIANGLE* trii=areatris.get_premier(itareatris);trii!=NULL;trii=areatris.get_suivant(itareatris))
  {
    for(int i=0;i<trii->get_segment1()->get_lien_triangle()->get_nb();i++)
    {
      MG_TRIANGLE* trineib=trii->get_segment1()->get_lien_triangle()->get(i);
      double avdis=(trineib->get_noeud1()->get_solution()+trineib->get_noeud2()->get_solution()+trineib->get_noeud3()->get_solution())/3.;
      if (avdis>tolerance)
        areatris.ajouter(trineib);
    }
    for(int i=0;i<trii->get_segment2()->get_lien_triangle()->get_nb();i++)
    {
      MG_TRIANGLE* trineib=trii->get_segment2()->get_lien_triangle()->get(i);
      double avdis=(trineib->get_noeud1()->get_solution()+trineib->get_noeud2()->get_solution()+trineib->get_noeud3()->get_solution())/3.;
      if (avdis>tolerance)
        areatris.ajouter(trineib);
    }
    for(int i=0;i<trii->get_segment3()->get_lien_triangle()->get_nb();i++)
    {
      MG_TRIANGLE* trineib=trii->get_segment3()->get_lien_triangle()->get(i);
      double avdis=(trineib->get_noeud1()->get_solution()+trineib->get_noeud2()->get_solution()+trineib->get_noeud3()->get_solution())/3.;
      if (avdis>tolerance)
        areatris.ajouter(trineib);
    }
  }
 } 
  double area=0.;
  multimap<double,MG_NOEUD*,std::greater<double> > maxamplitude;
  TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR itareatris;
  for(MG_TRIANGLE* tri=areatris.get_premier(itareatris);tri!=NULL;tri=areatris.get_suivant(itareatris))
  {
    maxamplitude.insert(pair<double,MG_NOEUD*> (tri->get_noeud1()->get_solution(),tri->get_noeud1()));
    maxamplitude.insert(pair<double,MG_NOEUD*> (tri->get_noeud2()->get_solution(),tri->get_noeud2()));
    maxamplitude.insert(pair<double,MG_NOEUD*> (tri->get_noeud3()->get_solution(),tri->get_noeud3()));
    OT_VECTEUR_3D v1(tri->get_noeud1()->get_coord(),tri->get_noeud2()->get_coord());
    OT_VECTEUR_3D v2(tri->get_noeud1()->get_coord(),tri->get_noeud3()->get_coord());
    OT_VECTEUR_3D area1=v1&v2;
    double area2=0.5*area1.get_longueur();
    area=area+area2;
  }
  cout<<"area= "<<area<<endl;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itmaxamplitude=maxamplitude.begin();
  cout<<"maxamplitude= "<< (*itmaxamplitude).first<<endl;
  
  MG_MAILLAGE* mai1;
  if (nummai==0) mai1=gest.get_mg_maillage(nummai); else mai1=gest.get_mg_maillageid(nummai);
  TPL_MAP_ENTITE<MG_TRIANGLE*> deltris;
  LISTE_MG_TRIANGLE::iterator ittri;
  for (MG_TRIANGLE* tri=mai1->get_premier_triangle(ittri);tri!=NULL;tri=mai1->get_suivant_triangle(ittri))
  {
  //int triexist=0;
    for(MG_TRIANGLE* triarea=areatris.get_premier(itareatris);triarea!=NULL;triarea=areatris.get_suivant(itareatris))
    {
      if (triarea==tri)
        tri->change_origine(IMPOSE);
        //triexist++;
    }
    //if (triexist<1)
      //deltris.ajouter(tri);
  }
 /* for(MG_TRIANGLE* tridel=deltris.get_premier(itareatris);tridel!=NULL;tridel=deltris.get_suivant(itareatris))
  {
    tridel->change_origine(IMPOSE);
    //mai1->supprimer_mg_triangleid(tridel->get_id());
  }*/
     
  gest.ajouter_mg_maillage(mai1);
  gest.enregistrer(areamagicfilename);
}
void CRIAQOPERATORS::samplepoints_compare(char* samplepoints1,char* samplepoints2,char* sortedsmplpnts)
{
  std::vector<double*> lstsamplepoints1;
  FILE *in1=fopen(samplepoints1,"rt");
  if (in1==NULL) affiche((char*) "file is not available");
  while(!feof(in1))
  {
    char chaine[500];
    char* res=fgets(chaine,500,in1);
    if(res!=NULL)
    {
      double x,y,z;
      double q1,q2,q3;
      double* pcoordbc=new double[6];
      int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
      if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
      else if(nb==6)
      {
        pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
        lstsamplepoints1.push_back(pcoordbc);
      }
    }
  }
  fclose(in1);
  std::vector<double*> lstsamplepoints2;
  FILE *in2=fopen(samplepoints2,"rt");
  if (in2==NULL) affiche((char*) "file is not available");
  while(!feof(in2))
  {
    char chaine[500];
    char* res=fgets(chaine,500,in2);
    if(res!=NULL)
    {
      double x,y,z;
      double q1,q2,q3;
      double* pcoordbc=new double[6];
      int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
      if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
      else if(nb==6)
      {
        pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
        lstsamplepoints2.push_back(pcoordbc);
      }
    }
  }
  fclose(in2);
  cout<<lstsamplepoints1.size()<<endl;
  cout<<lstsamplepoints2.size()<<endl;
  
  std::vector<double*> lstsortedsmplpnts;
  for(std::vector<double*>::iterator itlstsamplepoints1=lstsamplepoints1.begin();itlstsamplepoints1!=lstsamplepoints1.end();itlstsamplepoints1++)
  {
    double* pcoordbc=new double[7];
    int countsmpcomp=0;
    double* xyzsmp1=(*itlstsamplepoints1);
    for(std::vector<double*>::iterator itlstsamplepoints2=lstsamplepoints2.begin();itlstsamplepoints2!=lstsamplepoints2.end();itlstsamplepoints2++)
    {
      double* xyzsmp2=(*itlstsamplepoints2);
      if(xyzsmp1[0]==xyzsmp2[0] && xyzsmp1[1]==xyzsmp2[1] && xyzsmp1[2]==xyzsmp2[2] && xyzsmp1[3]==xyzsmp2[3] && xyzsmp1[4]==xyzsmp2[4] && xyzsmp1[5]==xyzsmp2[5])
        countsmpcomp++;
    }
    if(countsmpcomp==0)
    {
      pcoordbc[0]=0; pcoordbc[1]=xyzsmp1[0]; pcoordbc[2]=xyzsmp1[1]; pcoordbc[3]=xyzsmp1[2]; pcoordbc[4]=xyzsmp1[3]; pcoordbc[5]=xyzsmp1[4]; pcoordbc[6]=xyzsmp1[5];
      lstsortedsmplpnts.push_back(pcoordbc);
    }
    else if (countsmpcomp>0)
    {
      pcoordbc[0]=1; pcoordbc[1]=xyzsmp1[0]; pcoordbc[2]=xyzsmp1[1]; pcoordbc[3]=xyzsmp1[2]; pcoordbc[4]=xyzsmp1[3]; pcoordbc[5]=xyzsmp1[4]; pcoordbc[6]=xyzsmp1[5];
      lstsortedsmplpnts.push_back(pcoordbc);
    }
      
    
  }
  cout<<lstsortedsmplpnts.size()<<endl;
  
  
 FILE* insmp=fopen(sortedsmplpnts,"wt");
    for(std::vector<double*>::iterator itlstsortedsmplpnts=lstsortedsmplpnts.begin();itlstsortedsmplpnts!=lstsortedsmplpnts.end();itlstsortedsmplpnts++)
    {
      double* xyzsmpsort=(*itlstsortedsmplpnts);
      //cout<<xyzsmpsort[0]<<" "<<xyzsmpsort[1]<<endl;
      fprintf(insmp,"%le %le %le %le %le %le %le \n",xyzsmpsort[0],xyzsmpsort[1],xyzsmpsort[2],xyzsmpsort[3],xyzsmpsort[4],xyzsmpsort[5],xyzsmpsort[6]);
    }
    fclose(insmp);

}

void CRIAQOPERATORS::meshdistance_compare(char* referencemagicfilename,int refmeshno,char* comparemagicfilename,int compmeshno,char* outputcomparefile,char* solfilename)
{
  MG_FILE refgest(referencemagicfilename);
  MG_MAILLAGE* refmgmai;
  if (refmeshno==0) refmgmai=refgest.get_mg_maillage(refmeshno); else refmgmai=refgest.get_mg_maillageid(refmeshno);
  MG_FILE compgest(comparemagicfilename);
  MG_MAILLAGE* compmgmai;
  if (compmeshno==0) compmgmai=compgest.get_mg_maillage(compmeshno); else compmgmai=compgest.get_mg_maillageid(compmeshno);

  //octree initialization 
 double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
LISTE_MG_NOEUD::iterator it;
for (MG_NOEUD* no=refmgmai->get_premier_noeud(it);no!=NULL;no=refmgmai->get_suivant_noeud(it))
  {
  if (no->get_x()<xmin) xmin=no->get_x();  
  if (no->get_y()<ymin) ymin=no->get_y();  
  if (no->get_z()<zmin) zmin=no->get_z();  
  if (no->get_x()>xmax) xmax=no->get_x();  
  if (no->get_y()>ymax) ymax=no->get_y();  
  if (no->get_z()>zmax) zmax=no->get_z();
  }
for (MG_NOEUD* no=compmgmai->get_premier_noeud(it);no!=NULL;no=compmgmai->get_suivant_noeud(it))
  {
  if (no->get_x()<xmin) xmin=no->get_x();  
  if (no->get_y()<ymin) ymin=no->get_y();  
  if (no->get_z()<zmin) zmin=no->get_z();  
  if (no->get_x()>xmax) xmax=no->get_x();  
  if (no->get_y()>ymax) ymax=no->get_y();  
  if (no->get_z()>zmax) zmax=no->get_z();
  lstnoeud.ajouter(no);
  }  
 OT_VECTEUR_3D vecmin(xmin,ymin,zmin);OT_VECTEUR_3D vecmmax(xmax,ymax,zmax);
 OT_VECTEUR_3D vec(vecmmax,vecmin);
 
 //double search_radius=0.001*(vec.get_longueur());
 double search_radius=3.;
 OT_VECTEUR_3D min(xmin,ymin,zmin); OT_VECTEUR_3D max(xmax,ymax,zmax); OT_VECTEUR_3D average=(min+max)/2.; OT_VECTEUR_3D lengthvec(min,max);
 double length=sqrt(lengthvec*lengthvec); double bxr=1.1;
 xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
 xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);

TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octree;
octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
for (MG_NOEUD* no=compmgmai->get_premier_noeud(it);no!=NULL;no=compmgmai->get_suivant_noeud(it))
  octree.inserer(no);

MG_SOLUTION* mgsol=new MG_SOLUTION(refmgmai,1,solfilename,1,(char*)"Normaldistance",MAGIC::ENTITE_SOLUTION::ENTITE_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
refgest.ajouter_mg_solution(mgsol);
mgsol->change_legende(0,"Normaldistance");
int i=0;
for (MG_NOEUD* nop=refmgmai->get_premier_noeud(it);nop!=NULL;nop=refmgmai->get_suivant_noeud(it))
{
  TPL_MAP_ENTITE<MG_NOEUD*> lstneiclose;
    int neindnb=lstneiclose.get_nb();
    while(neindnb==0) 
    {
      octree.rechercher(nop->get_x(),nop->get_y(),nop->get_z(),search_radius,lstneiclose);
      neindnb=lstneiclose.get_nb();
    }
    map<double,MG_NOEUD*,std::less<double> > lstdist;
    TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itnop;
    for(MG_NOEUD* no=lstneiclose.get_premier(itnop);no!=NULL;no=lstneiclose.get_suivant(itnop))
    {
       OT_VECTEUR_3D w_n_tot(0.,0.,0.);
       double norme_w_n_tot=0.;
       OT_VECTEUR_3D w2_n_tot(0.,0.,0.);
       double norme_w2_n_tot=0.;
       //int nbtri=nop->get_lien_triangle()->get_nb();//put sol on cad model
       int nbtri=no->get_lien_triangle()->get_nb();
    for (int i=0;i<nbtri;i++)
    {
        //MG_TRIANGLE* tri=(MG_TRIANGLE*)nop->get_lien_triangle()->get(i);//put sol on cad model
        MG_TRIANGLE* tri=(MG_TRIANGLE*)no->get_lien_triangle()->get(i);
        double *xyzn1=tri->get_noeud1()->get_coord();
        double *xyzn2=tri->get_noeud2()->get_coord();
        double *xyzn3=tri->get_noeud3()->get_coord();
        
        OT_VECTEUR_3D vec1(xyzn1,xyzn3);
        OT_VECTEUR_3D vec2(xyzn1,xyzn2);
        OT_VECTEUR_3D n_tri=vec1&vec2;  // Normale au triangle

        double w=n_tri.get_longueur();  // 2 fois l'aire du triangle (pour éviter une division par 2)
        n_tri.norme();
        OT_VECTEUR_3D w_n=w*n_tri;
        double norme_w_n=w_n.get_longueur();
        w_n_tot=w_n_tot+w_n;
        norme_w_n_tot=norme_w_n_tot+norme_w_n;
    }
    
    OT_VECTEUR_3D n=w_n_tot/norme_w_n_tot;      // Normale à la surface au noeud
    n.norme();
    OT_VECTEUR_3D distvec(nop->get_coord(),no->get_coord());  
    double normaldistance=distvec*n; 
    double normaldistanceabs=fabs(normaldistance);
    //cout<<"normaldistanceabs= "<<normaldistanceabs<<endl;
      lstdist.insert(pair<double,MG_NOEUD*>(normaldistanceabs,nop));
    }
    map<double,MG_NOEUD*,std::less<double> >::iterator itlstdist=lstdist.begin();
    double distnsclos=(*itlstdist).first;
    MG_NOEUD* nodcls=(*itlstdist).second;
    mgsol->ecrire(distnsclos,i,0,0);
    i++;
}
refgest.enregistrer(outputcomparefile);
}
void CRIAQOPERATORS::deformed_correspond_mgmaiadd(char* magicfilename,int meshno,int numsol1,int numsol2,int numsol3,int numchamp1,int numchamp2,int numchamp3,double coef,char* correspondfilename)
{
MG_FILE gest(magicfilename);
        FEM_MAILLAGE* femmai;
        if (meshno==0)    femmai=gest.get_fem_maillage(meshno); else femmai=gest.get_fem_maillageid(meshno);
        FEM_SOLUTION* sol1=gest.get_fem_solutionid(numsol1);
        FEM_SOLUTION* sol2=gest.get_fem_solutionid(numsol2);
        FEM_SOLUTION* sol3=gest.get_fem_solutionid(numsol3);
        femmai->calcul_deforme(sol1,numchamp1,sol2,numchamp2,sol3,numchamp3);
        if (coef>1e-16)
          {
          cout<<"   Création du maillage deforme avec correspondence text file"<<endl;
MG_MAILLAGE* defmgmai=new MG_MAILLAGE(gest.get_mg_geometrie(0));
gest.ajouter_mg_maillage(defmgmai);
{
if (femmai->get_degre()==1)
  {
  //gest=femmai->get_mg_maillage()->get_gestionnaire();
  //gest->ajouter_mg_maillage(this);
  std::vector<unsigned long*> lstcores;//sasan
LISTE_FEM_NOEUD::iterator it;
  for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
    {
     double x=no->get_x(coef); 
     double y=no->get_y(coef); 
     double z=no->get_z(coef);
     MG_NOEUD *newno=new MG_NOEUD(no->get_lien_topologie(),x,y,z,DEFORME); 
     defmgmai->ajouter_mg_noeud(newno);
     no->change_numero(newno->get_id());
     unsigned long* correscad_defcad=new unsigned long[2];//sasan
     correscad_defcad[0]=no->get_mg_element_maillage()->get_id();//sasan
     correscad_defcad[1]=newno->get_id();//sasan
     lstcores.push_back(correscad_defcad);//sasan
    }
    //sasan
    FILE* in1=fopen(correspondfilename,"wt");
    fprintf(in1,"MG_MAILLAGE No. : deformed_MAILLAGE No.\n");
    fprintf(in1,"%lu %lu\n",femmai->get_mg_maillage()->get_id(),defmgmai->get_id());
    for(std::vector<unsigned long*>::iterator itlstpinsert=lstcores.begin();itlstpinsert!=lstcores.end();itlstpinsert++)
      {
        unsigned long* xyzpins1;
        xyzpins1=(*itlstpinsert);
        //cout<<" xyzpins1[0]= "<<xyzpins1[0]<<" xyzpins1[1]= "<<xyzpins1[1]<<endl;
        fprintf(in1,"%lu %lu\n",xyzpins1[0],xyzpins1[1]);
      }
  fclose(in1);
//sasan till
  LISTE_FEM_ELEMENT1::iterator it1;
  for (FEM_ELEMENT1 *ele=femmai->get_premier_element1(it1);ele!=NULL;ele=femmai->get_suivant_element1(it1))
      {
      MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
      MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(ele->get_lien_topologie(),no1,no2,DEFORME);
      defmgmai->ajouter_mg_segment(seg);
      }
  LISTE_FEM_ELEMENT2::iterator it2;
  for (FEM_ELEMENT2 *ele=femmai->get_premier_element2(it2);ele!=NULL;ele=femmai->get_suivant_element2(it2))
      {
      if (ele->get_nb_fem_noeud()==3)
        {
        MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
        MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
        MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
        defmgmai->ajouter_mg_triangle(ele->get_lien_topologie(),no1,no2,no3,DEFORME);
        }
      if (ele->get_nb_fem_noeud()==4)
        {
        MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
        MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
        MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
        MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
        defmgmai->ajouter_mg_quadrangle(ele->get_lien_topologie(),no1,no2,no3,no4,DEFORME);
        }
      }
  LISTE_FEM_ELEMENT3::iterator it3;
  for (FEM_ELEMENT3 *ele=femmai->get_premier_element3(it3);ele!=NULL;ele=femmai->get_suivant_element3(it3))
      {
      if (ele->get_nb_fem_noeud()==4)
        {
        MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
        MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
        MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
        MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
        defmgmai->ajouter_mg_tetra(ele->get_lien_topologie(),no1,no2,no3,no4,DEFORME);
        }
      if (ele->get_nb_fem_noeud()==8)
        {
        MG_NOEUD* no1=defmgmai->get_mg_noeudid(ele->get_fem_noeud(0)->get_numero());
        MG_NOEUD* no2=defmgmai->get_mg_noeudid(ele->get_fem_noeud(1)->get_numero());
        MG_NOEUD* no3=defmgmai->get_mg_noeudid(ele->get_fem_noeud(2)->get_numero());
        MG_NOEUD* no4=defmgmai->get_mg_noeudid(ele->get_fem_noeud(3)->get_numero());
        MG_NOEUD* no5=defmgmai->get_mg_noeudid(ele->get_fem_noeud(4)->get_numero());
        MG_NOEUD* no6=defmgmai->get_mg_noeudid(ele->get_fem_noeud(5)->get_numero());
        MG_NOEUD* no7=defmgmai->get_mg_noeudid(ele->get_fem_noeud(6)->get_numero());
        MG_NOEUD* no8=defmgmai->get_mg_noeudid(ele->get_fem_noeud(7)->get_numero());
        defmgmai->ajouter_mg_hexa(ele->get_lien_topologie(),no1,no2,no3,no4,no5,no6,no7,no8,DEFORME);
        }
      }
  
  }
}
          }
        cout<<"   Enregistrement"<<endl;
        gest.enregistrer(magicfilename);
}
void CRIAQOPERATORS::curvtur_crit_remvsmp(char* magicfilename,char* correspondfilename,char* gnifinspointfile,char* out_gnifinspointfile_removondefect,double search_radius,double curvdif_coef,int gmsh_affiche,double relet_search_rad,int relet_curvature)
{
  MG_FILE gest(magicfilename);
  map<unsigned long, unsigned long> coresspondndidlst;
  FILE* in1=fopen(correspondfilename,"rt");
  if (in1==NULL) cout<<"correspondfile is not available"<<endl;
  char mess[500];
  char *res=fgets(mess,500,in1);
  if(feof(in1))
    cout<<"correspondfile has not been well written"<<endl;
  res=fgets(mess,500,in1);
  if(feof(in1))
    cout<<"correspondfile has not been well written"<<endl;
  int cadmeshno;
  int defcadmeshno;
  int nb1=sscanf(mess,"%d %d",&cadmeshno,&defcadmeshno);
while (!feof(in1))
{
  char chaine[500];
  res=fgets(chaine,500,in1);
  if(res!=NULL)
  {
  unsigned long cadndid;
  unsigned long defcadndid;
  int nb=sscanf(chaine,"%lu %lu",&cadndid,&defcadndid);
  if (nb!=-1 && nb!=2) cout<<"Wrong file format"<<endl;
  else if (nb==2)
  {
    coresspondndidlst.insert(pair<unsigned long, unsigned long> (defcadndid,cadndid) );
  }
  }
}
fclose(in1);

cout<<"cadmeshno= "<<cadmeshno<<" defcadmeshno= "<<defcadmeshno<<endl;
  MG_MAILLAGE* cadmgmai;
  if (cadmeshno==0) cadmgmai=gest.get_mg_maillage(cadmeshno); else cadmgmai=gest.get_mg_maillageid(cadmeshno);
  MG_MAILLAGE* defcadmgmai;
  if (defcadmeshno==0) defcadmgmai=gest.get_mg_maillage(defcadmeshno); else defcadmgmai=gest.get_mg_maillageid(defcadmeshno);

cout<<"octree initialization "<<endl;
//CAD octree initialization
LISTE_MG_NOEUD::iterator itcadndcurv;
double xmincad=1e308,ymincad=1e308,zmincad=1e308,xmaxcad=-1e308,ymaxcad=-1e308,zmaxcad=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstcadnoeud;
for(MG_NOEUD* cadndcurv=cadmgmai->get_premier_noeud(itcadndcurv);cadndcurv!=NULL;cadndcurv=cadmgmai->get_suivant_noeud(itcadndcurv))
  {
  if (cadndcurv->get_x()<xmincad) xmincad=cadndcurv->get_x();  
  if (cadndcurv->get_y()<ymincad) ymincad=cadndcurv->get_y();  
  if (cadndcurv->get_z()<zmincad) zmincad=cadndcurv->get_z();  
  if (cadndcurv->get_x()>xmaxcad) xmaxcad=cadndcurv->get_x();  
  if (cadndcurv->get_y()>ymaxcad) ymaxcad=cadndcurv->get_y();  
  if (cadndcurv->get_z()>zmaxcad) zmaxcad=cadndcurv->get_z();
  lstcadnoeud.ajouter(cadndcurv);
  }
 OT_VECTEUR_3D vecmincad(xmincad,ymincad,zmincad);
 OT_VECTEUR_3D vecmmaxcad(xmaxcad,ymaxcad,zmaxcad);
 OT_VECTEUR_3D vecad(vecmmaxcad,vecmincad);
 OT_VECTEUR_3D mincad(xmincad,ymincad,zmincad); OT_VECTEUR_3D maxcad(xmaxcad,ymaxcad,zmaxcad); OT_VECTEUR_3D averagecad=(mincad+maxcad)/2.; OT_VECTEUR_3D lengthvecad(mincad,maxcad);
 double lengthcad=sqrt(lengthvecad*lengthvecad); double bxrcad=1.1;
 xmincad=averagecad.get_x()-(lengthcad*bxrcad);ymincad=averagecad.get_y()-(lengthcad*bxrcad);zmincad=averagecad.get_z()-(lengthcad*bxrcad);
 xmaxcad=averagecad.get_x()+(lengthcad*bxrcad);ymaxcad=averagecad.get_y()+(lengthcad*bxrcad);zmaxcad=averagecad.get_z()+(lengthcad*bxrcad);
 TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreecad;
 octreecad.initialiser(&lstcadnoeud,1,xmincad,ymincad,zmincad,xmaxcad,ymaxcad,zmaxcad);
 for(MG_NOEUD* cadndcurv=cadmgmai->get_premier_noeud(itcadndcurv);cadndcurv!=NULL;cadndcurv=cadmgmai->get_suivant_noeud(itcadndcurv))
  octreecad.inserer(cadndcurv);
//deformed_CAD octree initialization  
LISTE_MG_NOEUD::iterator itdefcadndcurv; 
double xmindefcad=1e308,ymindefcad=1e308,zmindefcad=1e308,xmaxdefcad=-1e308,ymaxdefcad=-1e308,zmaxdefcad=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstdefcadnoeud;
for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
  {
  if (defcadndcurv->get_x()<xmindefcad) xmindefcad=defcadndcurv->get_x();  
  if (defcadndcurv->get_y()<ymindefcad) ymindefcad=defcadndcurv->get_y();  
  if (defcadndcurv->get_z()<zmindefcad) zmindefcad=defcadndcurv->get_z();  
  if (defcadndcurv->get_x()>xmaxdefcad) xmaxdefcad=defcadndcurv->get_x();  
  if (defcadndcurv->get_y()>ymaxdefcad) ymaxdefcad=defcadndcurv->get_y();  
  if (defcadndcurv->get_z()>zmaxdefcad) zmaxdefcad=defcadndcurv->get_z();
  lstdefcadnoeud.ajouter(defcadndcurv);
  }
 OT_VECTEUR_3D vecmindefcad(xmindefcad,ymindefcad,zmindefcad);OT_VECTEUR_3D vecmmaxdefcad(xmaxdefcad,ymaxdefcad,zmaxdefcad);
 OT_VECTEUR_3D vecdefcad(vecmmaxdefcad,vecmindefcad);
 OT_VECTEUR_3D mindefcad(xmindefcad,ymindefcad,zmindefcad); OT_VECTEUR_3D maxdefcad(xmaxdefcad,ymaxdefcad,zmaxdefcad); OT_VECTEUR_3D averagedefcad=(mindefcad+maxdefcad)/2.; OT_VECTEUR_3D lengthvecdefcad(mindefcad,maxdefcad);
 double lengthdefcad=sqrt(lengthvecdefcad*lengthvecdefcad); double bxrdefcad=1.1;
 xmindefcad=averagedefcad.get_x()-(lengthdefcad*bxrdefcad);ymindefcad=averagedefcad.get_y()-(lengthdefcad*bxrdefcad);zmindefcad=averagedefcad.get_z()-(lengthdefcad*bxrdefcad);
 xmaxdefcad=averagedefcad.get_x()+(lengthdefcad*bxrdefcad);ymaxdefcad=averagedefcad.get_y()+(lengthdefcad*bxrdefcad);zmaxdefcad=averagedefcad.get_z()+(lengthdefcad*bxrdefcad);
 TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreedefcad;
 octreedefcad.initialiser(&lstdefcadnoeud,1,xmindefcad,ymindefcad,zmindefcad,xmaxdefcad,ymaxdefcad,zmaxdefcad);
 for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
  octreedefcad.inserer(defcadndcurv); 

  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_defcadmgmai;
  map<MG_NOEUD*,double > lstkmin_cad;
  map<MG_NOEUD*,double > lstkmax_cad;
  map<MG_NOEUD*,double > lstkgau_cad;
  map<MG_NOEUD*,double > lstkmin_defcad;
  map<MG_NOEUD*,double > lstkmax_defcad;
  map<MG_NOEUD*,double > lstkgau_defcad;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_defcadmgmai;
     std::map<std::string,double> liste_kmin;
     std::map<std::string,double> liste_kmax;
     std::map<std::string,double> liste_gau;
     std::map<std::string,double> liste_courburemax1;
     std::map<std::string,double> liste_courburemax2;
     std::map<std::string,double> liste_courburemax3;
     std::map<std::string,double> liste_courburemax4;
     std::map<std::string,double> liste_courburemax5;  
     int opensurafece=1;

cout<<"calcul courbure cadmgmai"<<endl;
// calcul courbure cadmgmai ********************************
 for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv)) 
     {
       double kmin1=0.;
       double kmax1=0.;
       double kmin2=0.;
       double kmax2=0.;
       MG_MAILLAGE_OUTILS *cou;
       cou->calcul_courbure_face_arete_sommet(cadmgnd,kmin1,kmax1,kmin2,kmax2,liste_kmin,liste_kmax,liste_gau,liste_courburemax1,liste_courburemax2,liste_courburemax3,liste_courburemax4,liste_courburemax5,opensurafece);
     }
 for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv)) 
   {
     
     unsigned long idface=cadmgnd->get_lien_topologie()->get_id();
     unsigned long idnoeud=cadmgnd->get_id();
     char mess[500];
     sprintf(mess,"%lu_%lu",idnoeud,idface);
     for(std::map<std::string,double>::iterator itkmin=liste_kmin.begin();itkmin!=liste_kmin.end();itkmin++)
     {
       if(mess==(*itkmin).first)
       {
         lstkmin_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,cadmgnd) );
         lstkmin_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmin).second));
       }
      }
     for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
     {
       if(mess==(*itkmax).first)
       {
         lstkmax_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,cadmgnd) );
         lstkmax_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkmax).second));
       }
    }
      for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
     {
       if(mess==(*itkgau).first)
       {
         lstkgau_cadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,cadmgnd) );
         lstkgau_cad.insert(pair<MG_NOEUD*,double>(cadmgnd,(*itkgau).second));
       }
    }
   }
//calc reletive curvature CAD
map<MG_NOEUD*,double >::iterator itlstkmin_cad;
map<MG_NOEUD*,double >::iterator itlstkmax_cad;
map<MG_NOEUD*,double >::iterator itlstkgau_cad;
 for (MG_NOEUD* cadmgnd=cadmgmai->get_premier_noeud(itcadndcurv);cadmgnd!=NULL;cadmgnd=cadmgmai->get_suivant_noeud(itcadndcurv)) 
{
  itlstkmin_cad=lstkmin_cad.find(cadmgnd);
  double relet_kmincad=(*itlstkmin_cad).second;
  itlstkmax_cad=lstkmax_cad.find(cadmgnd);
  double relet_kmaxcad=(*itlstkmax_cad).second;
  itlstkgau_cad=lstkgau_cad.find(cadmgnd);
  double relet_kgaucad=(*itlstkgau_cad).second;
  TPL_MAP_ENTITE<MG_NOEUD*> lstnei_relcad;
  int relcad_neindnb=lstnei_relcad.get_nb();
  double relet_search_rad1=relet_search_rad;
  while(relcad_neindnb<1) 
  {
    octreecad.rechercher(cadmgnd->get_x(),cadmgnd->get_y(),cadmgnd->get_z(),relet_search_rad1,lstnei_relcad);
    relcad_neindnb=lstnei_relcad.get_nb();
    relet_search_rad1=relet_search_rad1+0.1*relet_search_rad;
  }
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
  for(MG_NOEUD* cadndnei=lstnei_relcad.get_premier(itndnei);cadndnei!=NULL;cadndnei=lstnei_relcad.get_suivant(itndnei))
  {
    itlstkmin_cad=lstkmin_cad.find(cadndnei);
    relet_kmincad=relet_kmincad+(*itlstkmin_cad).second;
    itlstkmax_cad=lstkmax_cad.find(cadndnei);
    relet_kmaxcad=relet_kmaxcad+(*itlstkmax_cad).second;
    itlstkgau_cad=lstkgau_cad.find(cadndnei);
    relet_kgaucad=relet_kgaucad+(*itlstkgau_cad).second;
  }
  relet_kmincad=relet_kmincad/lstnei_relcad.get_nb();
  relet_kmaxcad=relet_kmaxcad/lstnei_relcad.get_nb();
  relet_kgaucad=relet_kgaucad/lstnei_relcad.get_nb();
  lstkmin_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmincad,cadmgnd));
  lstkmax_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmaxcad,cadmgnd));
  lstkgau_relt_cadmgmai.insert(pair<double,MG_NOEUD*> (relet_kgaucad,cadmgnd));
}
// affichage reletive sur GMSH ******************************   
   if (gmsh_affiche>1)
    {
      cout<<"affiche: calcul_courbure_cadmgmai_relet "<<endl;
      std::ostringstream oss;
      oss << "calcul_courbure_cadmgmai_curvdif_relet_coef:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol1);
     sol1->change_legende(0,"Kmin_relet");
     sol1->change_legende(1,"Kmax_relet");
     sol1->change_legende(2,"K_gau_relet");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
     {
       multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     double valeur0;
     double valeur1;
     double valeur2;
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no1==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no1==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no1==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,0);
     sol1->ecrire(valeur1,compte,1,0);
     sol1->ecrire(valeur2,compte,2,0);
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no2==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no2==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no2==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,1);
     sol1->ecrire(valeur1,compte,1,1);
     sol1->ecrire(valeur2,compte,2,1);
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no3==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no3==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no3==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,2);
     sol1->ecrire(valeur1,compte,1,2);
     sol1->ecrire(valeur2,compte,2,2);
     compte++;  
     }
    } 
   // affichage sur GMSH ******************************   
   if (gmsh_affiche>1)
    {
      cout<<"affiche: calcul_courbure_cadmgmai "<<endl;
      std::ostringstream oss;
      oss << "calcul_courbure_cadmgmai_curvdif_coef:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol1);
     sol1->change_legende(0,"Kmin");
     sol1->change_legende(1,"Kmax");
     sol1->change_legende(2,"K_gau");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
     {
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     char mess[500];
     sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
     double valeur0=liste_kmin[mess];
     double valeur1=liste_kmax[mess];
     double valeur2=liste_gau[mess];
     sol1->ecrire(valeur0,compte,0,0);
     sol1->ecrire(valeur1,compte,1,0);
     sol1->ecrire(valeur2,compte,2,0);
     sprintf(mess,"%lu_%lu",no2->get_id(),tri->get_lien_topologie()->get_id());
      valeur0=liste_kmin[mess];
      valeur1=liste_kmax[mess];
      valeur2=liste_gau[mess];
     sol1->ecrire(valeur0,compte,0,1);
     sol1->ecrire(valeur1,compte,1,1);
     sol1->ecrire(valeur2,compte,2,1);
     sprintf(mess,"%lu_%lu",no3->get_id(),tri->get_lien_topologie()->get_id());
      valeur0=liste_kmin[mess];
      valeur1=liste_kmax[mess];
      valeur2=liste_gau[mess];
     sol1->ecrire(valeur0,compte,0,2);
     sol1->ecrire(valeur1,compte,1,2);
     sol1->ecrire(valeur2,compte,2,2);
     compte++;  
     }
    }
   liste_kmin.clear();liste_kmax.clear();liste_gau.clear();liste_courburemax1.clear();liste_courburemax2.clear();liste_courburemax3.clear();liste_courburemax4.clear();liste_courburemax5.clear();

   cout<<"calcul courbure defcadmgmai"<<endl;
   // calcul courbure defcadmgmai ********************************

for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv)) 
     {
       double kmin1=0.;
       double kmax1=0.;
       double kmin2=0.;
       double kmax2=0.;
       MG_MAILLAGE_OUTILS *cou;
       cou->calcul_courbure_face_arete_sommet(defcadmgnd,kmin1,kmax1,kmin2,kmax2,liste_kmin,liste_kmax,liste_gau,liste_courburemax1,liste_courburemax2,liste_courburemax3,liste_courburemax4,liste_courburemax5,opensurafece);
     }
 for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv))  
   {
     unsigned long idface=defcadmgnd->get_lien_topologie()->get_id();
     unsigned long idnoeud=defcadmgnd->get_id();
     char mess[500];
     sprintf(mess,"%lu_%lu",idnoeud,idface);
     for(std::map<std::string,double>::iterator itkmin=liste_kmin.begin();itkmin!=liste_kmin.end();itkmin++)
     {
       if(mess==(*itkmin).first)
       {
         lstkmin_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmin).second,defcadmgnd) );
         lstkmin_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmin).second));
       }
      }
     for(std::map<std::string,double>::iterator itkmax=liste_kmax.begin();itkmax!=liste_kmax.end();itkmax++)
     {
       if(mess==(*itkmax).first)
       {
         lstkmax_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkmax).second,defcadmgnd) );
         lstkmax_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkmax).second));
       }
    }
     for(std::map<std::string,double>::iterator itkgau=liste_gau.begin();itkgau!=liste_gau.end();itkgau++)
     {
       if(mess==(*itkgau).first)
       {
         lstkgau_defcadmgmai.insert(pair<double,MG_NOEUD*>((*itkgau).second,defcadmgnd) );
         lstkgau_defcad.insert(pair<MG_NOEUD*,double>(defcadmgnd,(*itkgau).second));
       }
    }
   }
//calc reletive curvature defCAD
map<MG_NOEUD*,double >::iterator itlstkmin_defcad;
map<MG_NOEUD*,double >::iterator itlstkmax_defcad;
map<MG_NOEUD*,double >::iterator itlstkgau_defcad;
 for (MG_NOEUD* defcadmgnd=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadmgnd!=NULL;defcadmgnd=defcadmgmai->get_suivant_noeud(itdefcadndcurv)) 
{
  itlstkmin_defcad=lstkmin_defcad.find(defcadmgnd);
  double relet_kmindefcad=(*itlstkmin_defcad).second;
  itlstkmax_defcad=lstkmax_defcad.find(defcadmgnd);
  double relet_kmaxdefcad=(*itlstkmax_defcad).second;
  itlstkgau_defcad=lstkgau_cad.find(defcadmgnd);
  double relet_kgaudefcad=(*itlstkgau_defcad).second;
  TPL_MAP_ENTITE<MG_NOEUD*> lstnei_reldefcad;
  int reldefcad_neindnb=lstnei_reldefcad.get_nb();
  double relet_search_rad1=relet_search_rad;
  while(reldefcad_neindnb<1) 
  {
    octreecad.rechercher(defcadmgnd->get_x(),defcadmgnd->get_y(),defcadmgnd->get_z(),relet_search_rad1,lstnei_reldefcad);
    reldefcad_neindnb=lstnei_reldefcad.get_nb();
    relet_search_rad1=relet_search_rad1+0.1*relet_search_rad;
  }
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
  for(MG_NOEUD* defcadndnei=lstnei_reldefcad.get_premier(itndnei);defcadndnei!=NULL;defcadndnei=lstnei_reldefcad.get_suivant(itndnei))
  {
    itlstkmin_cad=lstkmin_cad.find(defcadndnei);
    relet_kmindefcad=relet_kmindefcad+(*itlstkmin_defcad).second;
    itlstkmax_cad=lstkmax_cad.find(defcadndnei);
    relet_kmaxdefcad=relet_kmaxdefcad+(*itlstkmax_defcad).second;
    itlstkgau_cad=lstkgau_cad.find(defcadndnei);
    relet_kgaudefcad=relet_kgaudefcad+(*itlstkgau_defcad).second;
  }
  relet_kmindefcad=relet_kmindefcad/lstnei_reldefcad.get_nb();
  relet_kmaxdefcad=relet_kmaxdefcad/lstnei_reldefcad.get_nb();
  relet_kgaudefcad=relet_kgaudefcad/lstnei_reldefcad.get_nb();
  lstkmin_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmindefcad,defcadmgnd));
  lstkmax_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kmaxdefcad,defcadmgnd));
  lstkgau_relt_defcadmgmai.insert(pair<double,MG_NOEUD*> (relet_kgaudefcad,defcadmgnd));
}
// affichage reletive defcad sur GMSH ******************************   
   if (gmsh_affiche>1)
    {
      cout<<"affiche: calcul_courbure_cadmgmai_relet "<<endl;
      std::ostringstream oss;
      oss << "calcul_courbure_cadmgmai_curvdif_relet_coef:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol1=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol1);
     sol1->change_legende(0,"Kmin_relet");
     sol1->change_legende(1,"Kmax_relet");
     sol1->change_legende(2,"K_gau_relet");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
     {
       multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     double valeur0;
     double valeur1;
     double valeur2;
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no1==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no1==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no1==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,0);
     sol1->ecrire(valeur1,compte,1,0);
     sol1->ecrire(valeur2,compte,2,0);
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no2==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no2==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no2==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,1);
     sol1->ecrire(valeur1,compte,1,1);
     sol1->ecrire(valeur2,compte,2,1);
     for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
     {
       if (no3==(*itlstkmin_relt_cadmgmai).second)
         valeur0=(*itlstkmin_relt_cadmgmai).first;
    }
    for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
     {
       if (no3==(*itlstkmax_relt_cadmgmai).second)
         valeur1=(*itlstkmax_relt_cadmgmai).first;
    }
    for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
     {
       if (no3==(*itlstkgau_relt_cadmgmai).second)
         valeur2=(*itlstkgau_relt_cadmgmai).first;
    }
     sol1->ecrire(valeur0,compte,0,2);
     sol1->ecrire(valeur1,compte,1,2);
     sol1->ecrire(valeur2,compte,2,2);
     compte++;  
     }
    } 
    // affichage sur GMSH ******************************   
    if (gmsh_affiche>1)
    {
      cout<<"affiche: calcul_courbure_defcadmgmai "<<endl;
       std::ostringstream oss;
      oss << "calcul_courbure_defcadmgmai_curvdif_coef:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol2=new MG_SOLUTION(defcadmgmai,3,cstr,1,(char*)"discrete_curboure_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol2);
     sol2->change_legende(0,"Kmin");
     sol2->change_legende(1,"Kmax");
     sol2->change_legende(2,"K_gau");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=defcadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=defcadmgmai->get_suivant_triangle(ittri))
     {
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     char mess[500];
     sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
     double valeur0=liste_kmin[mess];
     double valeur1=liste_kmax[mess];
     double valeur2=liste_gau[mess];
     sol2->ecrire(valeur0,compte,0,0);
     sol2->ecrire(valeur1,compte,1,0);
     sol2->ecrire(valeur2,compte,2,0);
     sprintf(mess,"%lu_%lu",no2->get_id(),tri->get_lien_topologie()->get_id());
      valeur0=liste_kmin[mess];
      valeur1=liste_kmax[mess];
      valeur2=liste_gau[mess];
     sol2->ecrire(valeur0,compte,0,1);
     sol2->ecrire(valeur1,compte,1,1);
     sol2->ecrire(valeur2,compte,2,1);
     sprintf(mess,"%lu_%lu",no3->get_id(),tri->get_lien_topologie()->get_id());
      valeur0=liste_kmin[mess];
      valeur1=liste_kmax[mess];
      valeur2=liste_gau[mess];
     sol2->ecrire(valeur0,compte,0,2);
     sol2->ecrire(valeur1,compte,1,2);
     sol2->ecrire(valeur2,compte,2,2);
     compte++;  
     }
  
    }
   liste_kmin.clear();liste_kmax.clear();liste_gau.clear();liste_courburemax1.clear();liste_courburemax2.clear();liste_courburemax3.clear();liste_courburemax4.clear();liste_courburemax5.clear();

   
   ///make correspondence between cadmgmai and defcadmgmai nodes
  map<unsigned long, unsigned long>::iterator itcoresspondndidlst;
  for(MG_NOEUD* defcadndcurv=defcadmgmai->get_premier_noeud(itdefcadndcurv);defcadndcurv!=NULL;defcadndcurv=defcadmgmai->get_suivant_noeud(itdefcadndcurv))
  {
    itcoresspondndidlst=coresspondndidlst.find(defcadndcurv->get_id());
    defcadndcurv->change_nouveau_numero((*itcoresspondndidlst).second);
  }
 
cout<<"solution calculation"<<endl;  
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_diff;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_diff;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_diff;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_defcadmgmai;
cout<<"lstkmin_cadmgmai.size()= "<<lstkmin_cadmgmai.size()<<endl;
cout<<"lstkmin_defcadmgmai.size()= "<<lstkmin_defcadmgmai.size()<<endl;
cout<<"kmin diff. calc. started"<<endl;

for(itlstkmin_cadmgmai=lstkmin_cadmgmai.begin();itlstkmin_cadmgmai!=lstkmin_cadmgmai.end();itlstkmin_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkmin_cadmgmai).second;
  double kmindcad;
  double kmindcad_def;
  kmindcad=(*itlstkmin_cadmgmai).first;
  for(itlstkmin_defcadmgmai=lstkmin_defcadmgmai.begin();itlstkmin_defcadmgmai!=lstkmin_defcadmgmai.end();itlstkmin_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkmin_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
      kmindcad_def=(*itlstkmin_defcadmgmai).first;
  }
  
  double kmin_diff=kmindcad_def-kmindcad; 
  lstkmin_diff.insert(pair<double,MG_NOEUD*>(kmin_diff,cadndcurv) );
}
cout<<"kmax diff. calc. started"<<endl;
for(itlstkmax_cadmgmai=lstkmax_cadmgmai.begin();itlstkmax_cadmgmai!=lstkmax_cadmgmai.end();itlstkmax_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkmax_cadmgmai).second;
  double kmaxndcad;
  double kmaxndcad_def;
  kmaxndcad=(*itlstkmax_cadmgmai).first;
  for(itlstkmax_defcadmgmai=lstkmax_defcadmgmai.begin();itlstkmax_defcadmgmai!=lstkmax_defcadmgmai.end();itlstkmax_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkmax_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
    {
      kmaxndcad_def=(*itlstkmax_defcadmgmai).first;
    }
  }
  double kmax_diff=kmaxndcad_def-kmaxndcad;
  lstkmax_diff.insert(pair<double,MG_NOEUD*>(kmax_diff,cadndcurv) );
}
cout<<"kgau diff. calc. started"<<endl;
for(itlstkgau_cadmgmai=lstkgau_cadmgmai.begin();itlstkgau_cadmgmai!=lstkgau_cadmgmai.end();itlstkgau_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkgau_cadmgmai).second;
  double kgaundcad;
  double kgaundcad_def;
  kgaundcad=(*itlstkgau_cadmgmai).first;
  for(itlstkgau_defcadmgmai=lstkgau_defcadmgmai.begin();itlstkgau_defcadmgmai!=lstkgau_defcadmgmai.end();itlstkgau_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkgau_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
    {
      kgaundcad_def=(*itlstkgau_defcadmgmai).first;
    }
  }
  double kgau_diff=kgaundcad_def-kgaundcad;
  lstkgau_diff.insert(pair<double,MG_NOEUD*>(kgau_diff,cadndcurv) );
}
// mean diff
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmin_relt_diff;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkmax_relt_diff;
  multimap<double,MG_NOEUD*,std::greater<double> > lstkgau_relt_diff;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_cadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_defcadmgmai;
  multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_defcadmgmai;
cout<<"kmin _relt diff. calc. started"<<endl;
  for(itlstkmin_relt_cadmgmai=lstkmin_relt_cadmgmai.begin();itlstkmin_relt_cadmgmai!=lstkmin_relt_cadmgmai.end();itlstkmin_relt_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkmin_relt_cadmgmai).second;
  double kmindcad;
  double kmindcad_relt_def;
  kmindcad=(*itlstkmin_relt_cadmgmai).first;
  for(itlstkmin_relt_defcadmgmai=lstkmin_relt_defcadmgmai.begin();itlstkmin_relt_defcadmgmai!=lstkmin_relt_defcadmgmai.end();itlstkmin_relt_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkmin_relt_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
      kmindcad_relt_def=(*itlstkmin_relt_defcadmgmai).first;
  }
  
  double kmin_relt_diff=kmindcad_relt_def-kmindcad; 
  lstkmin_relt_diff.insert(pair<double,MG_NOEUD*>(kmin_relt_diff,cadndcurv) );
}
cout<<"kmax _relt diff. calc. started"<<endl;
for(itlstkmax_relt_cadmgmai=lstkmax_relt_cadmgmai.begin();itlstkmax_relt_cadmgmai!=lstkmax_relt_cadmgmai.end();itlstkmax_relt_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkmax_relt_cadmgmai).second;
  double kmaxndcad;
  double kmaxndcad_relt_def;
  kmaxndcad=(*itlstkmax_relt_cadmgmai).first;
  for(itlstkmax_relt_defcadmgmai=lstkmax_relt_defcadmgmai.begin();itlstkmax_relt_defcadmgmai!=lstkmax_relt_defcadmgmai.end();itlstkmax_relt_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkmax_relt_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
    {
      kmaxndcad_relt_def=(*itlstkmax_relt_defcadmgmai).first;
    }
  }
  double kmax_relt_diff=kmaxndcad_relt_def-kmaxndcad;
  lstkmax_relt_diff.insert(pair<double,MG_NOEUD*>(kmax_relt_diff,cadndcurv) );
}
cout<<"kgau_relt diff. calc. started"<<endl;
for(itlstkgau_relt_cadmgmai=lstkgau_relt_cadmgmai.begin();itlstkgau_relt_cadmgmai!=lstkgau_relt_cadmgmai.end();itlstkgau_relt_cadmgmai++)
{
  MG_NOEUD* cadndcurv=(*itlstkgau_relt_cadmgmai).second;
  double kgaundcad;
  double kgaundcad_relt_def;
  kgaundcad=(*itlstkgau_relt_cadmgmai).first;
  for(itlstkgau_relt_defcadmgmai=lstkgau_relt_defcadmgmai.begin();itlstkgau_relt_defcadmgmai!=lstkgau_relt_defcadmgmai.end();itlstkgau_relt_defcadmgmai++)
  {
    MG_NOEUD* defcadndcurv=(*itlstkgau_relt_defcadmgmai).second;
    if(defcadndcurv->get_nouveau_numero()==cadndcurv->get_id())
    {
      kgaundcad_relt_def=(*itlstkgau_relt_defcadmgmai).first;
    }
  }
  double kgau_relt_diff=kgaundcad_relt_def-kgaundcad;
  lstkgau_relt_diff.insert(pair<double,MG_NOEUD*>(kgau_relt_diff,cadndcurv) );
}
// affichage  reletive diff sur GMSH ******************************
if (gmsh_affiche>0)
    {
      cout<<"affiche: calcul_relet_difference_courbure "<<endl;
      std::ostringstream oss;
      oss << "KMaxMinGaudiff_curvdifcoef_relet:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol4=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"Kdiff_relt_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol4);
     sol4->change_legende(0,"Kmin_diff");
     sol4->change_legende(1,"Kmax_diff");
     sol4->change_legende(2,"Kgau_diff");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
     {
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diff;
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diff;
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_relt_diff;
     double valeurKmin1, valeurKmin2, valeurKmin3;
     double valeurKmax1, valeurKmax2, valeurKmax3;
     double valeurKgau1, valeurKgau2, valeurKgau3;
     for (itlstkmin_relt_diff=lstkmin_relt_diff.begin();itlstkmin_relt_diff!=lstkmin_relt_diff.end();itlstkmin_relt_diff++)
     {
       if (no1==(*itlstkmin_relt_diff).second)
         valeurKmin1=(*itlstkmin_relt_diff).first;
       if (no2==(*itlstkmin_relt_diff).second)
         valeurKmin2=(*itlstkmin_relt_diff).first;
       if (no3==(*itlstkmin_relt_diff).second)
         valeurKmin3=(*itlstkmin_relt_diff).first;
    }
    for (itlstkmax_relt_diff=lstkmax_relt_diff.begin();itlstkmax_relt_diff!=lstkmax_relt_diff.end();itlstkmax_relt_diff++)
    {
      if (no1==(*itlstkmax_relt_diff).second)
         valeurKmax1=(*itlstkmax_relt_diff).first;
      if (no2==(*itlstkmax_relt_diff).second)
         valeurKmax2=(*itlstkmax_relt_diff).first;
      if (no3==(*itlstkmax_relt_diff).second)
         valeurKmax3=(*itlstkmax_relt_diff).first;
    }
   for (itlstkgau_relt_diff=lstkgau_relt_diff.begin();itlstkgau_relt_diff!=lstkgau_relt_diff.end();itlstkgau_relt_diff++)
    {
      if (no1==(*itlstkgau_relt_diff).second)
         valeurKgau1=(*itlstkgau_relt_diff).first;
      if (no2==(*itlstkgau_relt_diff).second)
         valeurKgau2=(*itlstkgau_relt_diff).first;
      if (no3==(*itlstkgau_relt_diff).second)
         valeurKgau3=(*itlstkgau_relt_diff).first;
    }
     sol4->ecrire(valeurKmin1,compte,0,0);
     sol4->ecrire(valeurKmax1,compte,1,0);
     sol4->ecrire(valeurKgau1,compte,2,0);

     sol4->ecrire(valeurKmin2,compte,0,1);
     sol4->ecrire(valeurKmax2,compte,1,1);
     sol4->ecrire(valeurKgau2,compte,2,1);

     sol4->ecrire(valeurKmin3,compte,0,2);
     sol4->ecrire(valeurKmax3,compte,1,2);
     sol4->ecrire(valeurKgau3,compte,2,2);
     compte++;  
     }
    }

// affichage sur GMSH ******************************   
if (gmsh_affiche>1)
    {
      cout<<"affiche: calcul_difference_courbure "<<endl;
       std::ostringstream oss;
      oss << "KMaxMinGaudiff_curvdifcoef:" <<curvdif_coef ;
      std::string namefic = oss.str();
      char *cstr = new char[namefic.length() + 1];
      strcpy(cstr, namefic.c_str());
     MG_SOLUTION* sol3=new MG_SOLUTION(cadmgmai,3,cstr,1,(char*)"Kdiff_solution",MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT2_NOEUD,MAGIC::TYPE_SOLUTION::SCALAIRE);
     gest.ajouter_mg_solution(sol3);
     sol3->change_legende(0,"Kmin_diff");
     sol3->change_legende(1,"Kmax_diff");
     sol3->change_legende(2,"Kgau_diff");
     LISTE_MG_TRIANGLE::iterator ittri;
     int compte=0;
     for (MG_TRIANGLE* tri=cadmgmai->get_premier_triangle(ittri);tri!=NULL;tri=cadmgmai->get_suivant_triangle(ittri))
     {
     MG_NOEUD *no1=tri->get_noeud1();
     MG_NOEUD *no2=tri->get_noeud2();
     MG_NOEUD *no3=tri->get_noeud3();
     //char mess[500];
     //sprintf(mess,"%lu_%lu",no1->get_id(),tri->get_lien_topologie()->get_id());
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_diff;
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_diff;
     multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkgau_diff;
     double valeurKmin1, valeurKmin2, valeurKmin3;
     double valeurKmax1, valeurKmax2, valeurKmax3;
     double valeurKgau1, valeurKgau2, valeurKgau3;
     for (itlstkmin_diff=lstkmin_diff.begin();itlstkmin_diff!=lstkmin_diff.end();itlstkmin_diff++)
     {
       if (no1==(*itlstkmin_diff).second)
         valeurKmin1=(*itlstkmin_diff).first;
       if (no2==(*itlstkmin_diff).second)
         valeurKmin2=(*itlstkmin_diff).first;
       if (no3==(*itlstkmin_diff).second)
         valeurKmin3=(*itlstkmin_diff).first;
    }
    for (itlstkmax_diff=lstkmax_diff.begin();itlstkmax_diff!=lstkmax_diff.end();itlstkmax_diff++)
    {
      if (no1==(*itlstkmax_diff).second)
         valeurKmax1=(*itlstkmax_diff).first;
      if (no2==(*itlstkmax_diff).second)
         valeurKmax2=(*itlstkmax_diff).first;
      if (no3==(*itlstkmax_diff).second)
         valeurKmax3=(*itlstkmax_diff).first;
    }
   for (itlstkgau_diff=lstkgau_diff.begin();itlstkgau_diff!=lstkgau_diff.end();itlstkgau_diff++)
    {
      if (no1==(*itlstkgau_diff).second)
         valeurKgau1=(*itlstkgau_diff).first;
      if (no2==(*itlstkgau_diff).second)
         valeurKgau2=(*itlstkgau_diff).first;
      if (no3==(*itlstkgau_diff).second)
         valeurKgau3=(*itlstkgau_diff).first;
    }
     sol3->ecrire(valeurKmin1,compte,0,0);
     sol3->ecrire(valeurKmax1,compte,1,0);
     sol3->ecrire(valeurKgau1,compte,2,0);

     sol3->ecrire(valeurKmin2,compte,0,1);
     sol3->ecrire(valeurKmax2,compte,1,1);
     sol3->ecrire(valeurKgau2,compte,2,1);

     sol3->ecrire(valeurKmin3,compte,0,2);
     sol3->ecrire(valeurKmax3,compte,1,2);
     sol3->ecrire(valeurKgau3,compte,2,2);
     compte++;  
     }
  
    }

     char chaine[1000];
   char* p;
   p=strrchr(magicfilename,'.');
   strncpy(chaine,magicfilename,p-magicfilename);
   std::string namfic=chaine;
   namfic=namfic  + "courbcal.magic";
   gest.enregistrer(namfic.c_str());
    cout<<"read insert points in lstpinsert"<<endl;  
//read insert points in lstpinsert
  std::vector<double*> lstpinsert;
double prop=0.001;
FILE *in=fopen(gnifinspointfile,"rt");
if (in==NULL) cout<<"file is not available"<<endl;
while (!feof(in))
{
  char chaine[500];
  char* res=fgets(chaine,500,in);
  if(res!=NULL)
  {
  double x,y,z;
  double q1,q2,q3;
  int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
  q1=prop*q1;
  q2=prop*q2;
  q3=prop*q3;
  if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
  else if (nb==6)
  {
  double* pcoordbc=new double[6];
  pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
  lstpinsert.push_back(pcoordbc);
  }
  }
}
fclose(in);
if (relet_curvature<1)
{
//criterion
multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmin_diff;
multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmax_diff;
multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_diff;
multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_diff;
ritlstkmin_diff=lstkmin_diff.rbegin();
double kmin_diff_min=(*ritlstkmin_diff).first;
itlstkmin_diff=lstkmin_diff.begin();
double kmin_diff_max=(*itlstkmin_diff).first;
ritlstkmax_diff=lstkmax_diff.rbegin();
double kmax_diff_min=(*ritlstkmax_diff).first;
itlstkmax_diff=lstkmax_diff.begin();
double kmax_diff_max=(*itlstkmax_diff).first;
cout<<"kmin_diff_min= "<<kmin_diff_min<<" kmin_diff_max= "<<kmin_diff_max<<endl;
cout<<"kmax_diff_min= "<<kmax_diff_min<<" kmax_diff_max= "<<kmax_diff_max<<endl;
  
  double crikmin_min=curvdif_coef;//sasan*kmin_diff_min;
  double crikmin_max=curvdif_coef;//sasan*kmin_diff_max;
  double crikmax_min=curvdif_coef;//sasan*kmax_diff_min;
  double crikmax_max=curvdif_coef;//sasan*kmax_diff_max;
cout<<"crikmin_min= "<<crikmin_min<<" crikmin_max= "<<crikmin_max<<" crikmax_min= "<<crikmax_min<<" crikmax_max= "<<crikmax_max<<endl;
cout<<"lstkmin_diff.size()= "<<lstkmin_diff.size()<<endl;
cout<<"lstkmax_diff.size()= "<<lstkmax_diff.size()<<endl;

cout<<"for Kmin"<<endl;
int Kmin_positive_nb=0;
double Kmin_positive_mean=0.;
double Kmin_positive_stdev=0.;
int Kmin_negative_nb=0;
double Kmin_negative_mean=0.;
double Kmin_negative_stdev=0.;
for (itlstkmin_diff=lstkmin_diff.begin();itlstkmin_diff!=lstkmin_diff.end();itlstkmin_diff++)
{
  if(((*itlstkmin_diff).first)>=0.)
  {
    Kmin_positive_nb++;
    Kmin_positive_mean += (*itlstkmin_diff).first;
  }
  if(((*itlstkmin_diff).first)<=0.)
  {
    Kmin_negative_nb++;
    Kmin_negative_mean += (*itlstkmin_diff).first;
  }
}
Kmin_positive_mean /= (1.*Kmin_positive_nb);
Kmin_negative_mean /= (1.*Kmin_negative_nb);
cout<<"Kmin_positive_mean= "<<Kmin_positive_mean<<" Kmin_negative_mean= "<<Kmin_negative_mean<<endl;

TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
for (itlstkmin_diff=lstkmin_diff.begin();itlstkmin_diff!=lstkmin_diff.end();itlstkmin_diff++)
{
  {
  if((crikmin_max<(*itlstkmin_diff).first) || (crikmin_min>(*itlstkmin_diff).first))
  {
    MG_NOEUD* piefect=(*itlstkmin_diff).second;
    int neindnb=lstneipdefect.get_nb();
    double search_radius1=search_radius;
    while(neindnb==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
      neindnb=lstneipdefect.get_nb();
      search_radius1=search_radius1+0.1*search_radius;
    }
    TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
    for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
    {
      for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins;
        xyzpins=(*itlstpinsert);
        if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
          lstpinsert.erase(itlstpinsert);
      }
    }
    lstneipdefect.vide();
  }
  }
}
cout<<"for Kmax"<<endl;
int Kmax_positive_nb=0;
double Kmax_positive_mean=0.;
double Kmax_positive_stdev=0.;
int Kmax_negative_nb=0;
double Kmax_negative_mean=0.;
double Kmax_negative_stdev=0.;
for (itlstkmax_diff=lstkmax_diff.begin();itlstkmax_diff!=lstkmax_diff.end();itlstkmax_diff++)
{
  if(((*itlstkmax_diff).first)>=0.)
  {
    Kmax_positive_nb++;
    Kmax_positive_mean += (*itlstkmax_diff).first;
  }
  if(((*itlstkmax_diff).first)<=0.)
  {
    Kmax_negative_nb++;
    Kmax_negative_mean += (*itlstkmax_diff).first;
  }
}
Kmax_positive_mean /= (1.*Kmax_positive_nb);
Kmax_negative_mean /= (1.*Kmax_negative_nb);
cout<<"Kmax_positive_mean= "<<Kmax_positive_mean<<" Kmax_negative_mean= "<<Kmax_negative_mean<<endl;

for (itlstkmax_diff=lstkmax_diff.begin();itlstkmax_diff!=lstkmax_diff.end();itlstkmax_diff++)
{
  if((crikmax_max<(*itlstkmax_diff).first) || (crikmax_min>(*itlstkmax_diff).first))
  {
    MG_NOEUD* piefect=(*itlstkmax_diff).second;
    TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
    int neindnb=lstneipdefect.get_nb();
    double search_radius1=search_radius;
    while(neindnb==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
      neindnb=lstneipdefect.get_nb();
      search_radius1=search_radius1+0.1*search_radius;
    }
    TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
    for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
    {
      for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins;
        xyzpins=(*itlstpinsert);
         if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
           lstpinsert.erase(itlstpinsert);
      }
    }
  
  }
}
}
if (relet_curvature>0)
{
//criterion
multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmin_relt_diff=lstkmin_relt_diff.rbegin();
multimap<double,MG_NOEUD*,std::greater<double> >::reverse_iterator ritlstkmax_relt_diff=lstkmax_relt_diff.rbegin(); 
multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diff=lstkmin_relt_diff.begin();
multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diff=lstkmax_relt_diff.begin();
double kmin_relt_diff_min=(*ritlstkmin_relt_diff).first;
double kmin_relt_diff_max=(*itlstkmin_relt_diff).first;
double kmax_relt_diff_min=(*ritlstkmax_relt_diff).first;
double kmax_relt_diff_max=(*itlstkmax_relt_diff).first;
cout<<"kmin_relt_diff_min= "<<kmin_relt_diff_min<<" kmin_relt_diff_max= "<<kmin_relt_diff_max<<endl;
cout<<"kmax_relt_diff_min= "<<kmax_relt_diff_min<<" kmax_relt_diff_max= "<<kmax_relt_diff_max<<endl;

  double crikmin_min;
  double crikmin_max;
  double crikmax_min;
  double crikmax_max;
  
cout<<"crikmin_min= "<<crikmin_min<<" crikmin_max= "<<crikmin_max<<" crikmax_min= "<<crikmax_min<<" crikmax_max= "<<crikmax_max<<endl;
cout<<"lstkmin_diff.size()= "<<lstkmin_diff.size()<<endl;
cout<<"lstkmax_diff.size()= "<<lstkmax_diff.size()<<endl;

cout<<"for Kmin_relet_diff"<<endl;
cout<<"lstkmin_relt_diff.size()= "<<lstkmin_relt_diff.size()<<endl;
int Kmin_positive_nb=0;
double Kmin_positive_mean=0.;
double Kmin_positive_stdev=0.;
int Kmin_negative_nb=0;
double Kmin_negative_mean=0.;
double Kmin_negative_stdev=0.;
int kmin_nb=0;
double kmin_mean=0.;
for (itlstkmin_relt_diff=lstkmin_relt_diff.begin();itlstkmin_relt_diff!=lstkmin_relt_diff.end();itlstkmin_relt_diff++)
{
  kmin_nb++;
  kmin_mean += (*itlstkmin_relt_diff).first;
  if(((*itlstkmin_relt_diff).first)>=0.)
  {
    Kmin_positive_nb++;
    Kmin_positive_mean += (*itlstkmin_relt_diff).first;
  }
  if(((*itlstkmin_relt_diff).first)<=0.)
  {
    Kmin_negative_nb++;
    Kmin_negative_mean += (*itlstkmin_relt_diff).first;
  }
}
Kmin_positive_mean /= (1.*Kmin_positive_nb);
Kmin_negative_mean /= (1.*Kmin_negative_nb);
kmin_mean /= (1.*kmin_nb);
cout<<"Kmin_positive_mean= "<<Kmin_positive_mean<<" Kmin_negative_mean= "<<Kmin_negative_mean<<" Kmin_mean= "<<kmin_mean<<endl;
crikmin_min=curvdif_coef*Kmin_negative_mean;
crikmin_max=curvdif_coef*Kmin_positive_mean;


for (itlstkmin_relt_diff=lstkmin_relt_diff.begin();itlstkmin_relt_diff!=lstkmin_relt_diff.end();itlstkmin_relt_diff++)
{
  if((crikmin_max<(*itlstkmin_relt_diff).first) || (crikmin_min>(*itlstkmin_relt_diff).first))
  {
    MG_NOEUD* piefect=(*itlstkmin_relt_diff).second;
    TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
    /*int neindnb=lstneipdefect.get_nb();
    double search_radius1=1.;
    while(neindnb==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
      neindnb=lstneipdefect.get_nb();
      search_radius1=search_radius1+0.1*1.;
    }*/
    for (int itt=0;itt<piefect->get_lien_segment()->get_nb();itt++)
    {
      MG_SEGMENT* segpiefect=piefect->get_lien_segment()->get(itt);
      if (segpiefect->get_noeud1()!=piefect)
        lstneipdefect.ajouter(segpiefect->get_noeud1());
      else if (segpiefect->get_noeud2()!=piefect)
        lstneipdefect.ajouter(segpiefect->get_noeud2());
    }
    
    
    TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
    int locnd_nb=0;
    double loc_curv=0.;
    for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
    {
      double curv=0.;
      if(ndnei!=piefect)
      {
        for (multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diffnei=lstkmin_relt_diff.begin();itlstkmin_relt_diffnei!=lstkmin_relt_diff.end();itlstkmin_relt_diffnei++)
        {
          if ((*itlstkmin_relt_diffnei).second==ndnei)
          {
            locnd_nb++;
            curv += (*itlstkmin_relt_diffnei).first;
          }
        }
      }
      curv /= locnd_nb;
      loc_curv=curv;
    }

int pass_kmin=0;    
if ((*itlstkmin_relt_diff).first>=0)
{
  if((*itlstkmin_relt_diff).first<=4.*loc_curv)
    pass_kmin++;
}
else if ((*itlstkmin_relt_diff).first<=0)
{
  if((*itlstkmin_relt_diff).first>=4.*loc_curv)
    pass_kmin++;
}  
if (pass_kmin>0)//(loc_curv>(0.2*(*itlstkmin_relt_diff).first) || loc_curv<(0.2*(*itlstkmin_relt_diff).first))    
{
  cout<<"loc_curv= "<<loc_curv<<" (*itlstkmin_relt_diff).first= "<<(*itlstkmin_relt_diff).first<<endl;
  TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect2;
    int neindnb2=lstneipdefect2.get_nb();
    double search_radius2=search_radius;
    while(neindnb2==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius2,lstneipdefect2);
      neindnb2=lstneipdefect2.get_nb();
      search_radius2=search_radius2+0.1*search_radius;
    }
    for(MG_NOEUD* ndnei=lstneipdefect2.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect2.get_suivant(itndnei))
    {
      double curv=0.;
      for (multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmin_relt_diffnei=lstkmin_relt_diff.begin();itlstkmin_relt_diffnei!=lstkmin_relt_diff.end();itlstkmin_relt_diffnei++)
        {
          if ((*itlstkmin_relt_diffnei).second==ndnei)
          {
            curv =(*itlstkmin_relt_diffnei).first;
          }
        }
        if(crikmin_max<curv || crikmin_min>curv)
      {     
      for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins;
        xyzpins=(*itlstpinsert);
        if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
        {
          lstpinsert.erase(itlstpinsert);
          itlstpinsert=lstpinsert.begin();
        }
      }  
      }
    }
    lstneipdefect.vide();
}  
  }
}
cout<<"for Kmax_relet_diff"<<endl;
int Kmax_positive_nb=0;
double Kmax_positive_mean=0.;
double Kmax_positive_stdev=0.;
int Kmax_negative_nb=0;
double Kmax_negative_mean=0.;
double Kmax_negative_stdev=0.;
int Kmax_nb=0;
double Kmax_mean=0.;
for (itlstkmax_relt_diff=lstkmax_relt_diff.begin();itlstkmax_relt_diff!=lstkmax_relt_diff.end();itlstkmax_relt_diff++)
{
  Kmax_nb++;
  Kmax_mean += (*itlstkmax_relt_diff).first;
  if(((*itlstkmax_relt_diff).first)>=0.)
  {
    Kmax_positive_nb++;
    Kmax_positive_mean += (*itlstkmax_relt_diff).first;
  }
  if(((*itlstkmax_relt_diff).first)<=0.)
  {
    Kmax_negative_nb++;
    Kmax_negative_mean += (*itlstkmax_relt_diff).first;
  }
}
Kmax_positive_mean /= (1.*Kmax_positive_nb);
Kmax_negative_mean /= (1.*Kmax_negative_nb);
Kmax_mean /= (1.*Kmax_nb);
cout<<"Kmax_positive_mean= "<<Kmax_positive_mean<<" Kmax_negative_mean= "<<Kmax_negative_mean<<" Kmax_mean= "<<Kmax_mean<<endl;
crikmax_min=curvdif_coef*Kmax_negative_mean;
crikmax_max=curvdif_coef*Kmax_positive_mean;

for (itlstkmax_relt_diff=lstkmax_relt_diff.begin();itlstkmax_relt_diff!=lstkmax_relt_diff.end();itlstkmax_relt_diff++)
{
  if((crikmax_max<(*itlstkmax_relt_diff).first) || (crikmax_min>(*itlstkmax_relt_diff).first))
  {
    MG_NOEUD* piefect=(*itlstkmax_relt_diff).second;
    TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect;
    /*int neindnb=lstneipdefect.get_nb();
    double search_radius1=1.;
    while(neindnb==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius1,lstneipdefect);
      neindnb=lstneipdefect.get_nb();
      search_radius1=search_radius1+0.1*1.;
    }*/
    for (int itt=0;itt<piefect->get_lien_segment()->get_nb();itt++)
    {
      MG_SEGMENT* segpiefect=piefect->get_lien_segment()->get(itt);
      if (segpiefect->get_noeud1()!=piefect)
        lstneipdefect.ajouter(segpiefect->get_noeud1());
      else if (segpiefect->get_noeud2()!=piefect)
        lstneipdefect.ajouter(segpiefect->get_noeud2());
    }
       
    
    TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itndnei;
    int locnd_nb=0;
    double loc_curv=0.;
    for(MG_NOEUD* ndnei=lstneipdefect.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect.get_suivant(itndnei))
    {
      double curv=0.;
      if(ndnei!=piefect)
      {
        for (multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diffnei=lstkmax_relt_diff.begin();itlstkmax_relt_diffnei!=lstkmax_relt_diff.end();itlstkmax_relt_diffnei++)
        {
          if ((*itlstkmax_relt_diffnei).second==ndnei)
          {
            locnd_nb++;
            curv += (*itlstkmax_relt_diffnei).first;
          }
        }
      }
      curv /= locnd_nb;
      loc_curv=curv;
    }

int pass_kmax=0;    
if ((*itlstkmax_relt_diff).first>=0)
{
  if((*itlstkmax_relt_diff).first<=4.*loc_curv)
    pass_kmax++;
}
else if ((*itlstkmax_relt_diff).first<=0)
{
  if((*itlstkmax_relt_diff).first>=4.*loc_curv)
    pass_kmax++;
}  
    
if (pass_kmax>0)   
{
  cout<<"loc_curv= "<<loc_curv<<" (*itlstkmax_relt_diff).first= "<<(*itlstkmax_relt_diff).first<<endl;
  TPL_MAP_ENTITE<MG_NOEUD*> lstneipdefect2;
    int neindnb2=lstneipdefect2.get_nb();
    double search_radius2=search_radius;
    while(neindnb2==0) 
    {
      octreecad.rechercher(piefect->get_x(),piefect->get_y(),piefect->get_z(),search_radius2,lstneipdefect2);
      neindnb2=lstneipdefect2.get_nb();
      search_radius2=search_radius2+0.1*search_radius;
    }
    for(MG_NOEUD* ndnei=lstneipdefect2.get_premier(itndnei);ndnei!=NULL;ndnei=lstneipdefect2.get_suivant(itndnei))
    {
      double curv=0.;
      for (multimap<double,MG_NOEUD*,std::greater<double> >::iterator itlstkmax_relt_diffnei=lstkmax_relt_diff.begin();itlstkmax_relt_diffnei!=lstkmax_relt_diff.end();itlstkmax_relt_diffnei++)
      {
        if ((*itlstkmax_relt_diffnei).second==ndnei)
        {
          curv=(*itlstkmax_relt_diffnei).first;
        }
      }
      if(crikmax_max<curv || crikmax_min>curv)
      {
      for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins;
        xyzpins=(*itlstpinsert);
         if(ndnei->get_x()==xyzpins[0] && ndnei->get_y()==xyzpins[1] && ndnei->get_z()==xyzpins[2])
         {
          lstpinsert.erase(itlstpinsert);
          itlstpinsert=lstpinsert.begin();
        }
      }
      }
    }
}  
  }
}
}
FILE* in2=fopen(out_gnifinspointfile_removondefect,"wt");
int countfuck=0;
for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins1;
        xyzpins1=(*itlstpinsert);
        fprintf(in2,"%le %le %le %le %le %le\n",xyzpins1[0],xyzpins1[1],xyzpins1[2],xyzpins1[3]*1000.,xyzpins1[4]*1000.,xyzpins1[5]*1000.);
      }
  fclose(in2);
}
void CRIAQOPERATORS::vm_crit_remvsmp(int femmeshno,char* magicfilename,int femsolid,int femsubsolno,char* gnifinspointfile,char* out_gnifinspointfile_removondefect,double search_radius,double vm_coef)
{
  MG_FILE gest(magicfilename);
  FEM_MAILLAGE* femmai;
  if (femmeshno==0) femmai=gest.get_fem_maillage(femmeshno); else femmai=gest.get_fem_maillageid(femmeshno);
  FEM_SOLUTION* femsol=gest.get_fem_solutionid(femsolid);
  femsol->active_solution(femsubsolno);

std::map<double,FEM_NOEUD*,std::greater<double> > mpsol;
LISTE_FEM_NOEUD::iterator itfemnds;
  for(FEM_NOEUD* nod=femmai->get_premier_noeud(itfemnds);nod!=NULL;nod=femmai->get_suivant_noeud(itfemnds))
  {
    //if(((MG_NOEUD*) nod->get_mg_element_maillage())->get_origine()==IMPOSE) // only calculate VM_mean for inserted nodes its good for models with noise but not good for models without noise
      mpsol.insert(pair<double,FEM_NOEUD*> (nod->get_solution(),nod));

  }
  std::map<double,FEM_NOEUD*,std::greater<double> >::iterator itmpsol;
  double vm_sum=0.;
  for(itmpsol=mpsol.begin();itmpsol!=mpsol.end();itmpsol++)
  {
    double vm=(*itmpsol).first;
    vm_sum=vm_sum+vm;
  }
  double vm_mean=vm_sum/(mpsol.size());
  double vm_sigma=0.;
  for(itmpsol=mpsol.begin();itmpsol!=mpsol.end();itmpsol++)
  {
    double vm=(*itmpsol).first;
    vm_sigma=vm_sigma+((vm-vm_mean)*(vm-vm_mean));
  }
  vm_sigma=vm_sigma/(mpsol.size());
  double vm_stdev=sqrt(vm_sigma);
itmpsol=mpsol.begin();

cout<<"VM_maximum= "<<(*itmpsol).first<<" vm_stdev= "<<vm_stdev<<" vm_mean= "<<vm_mean<<" VM_max_overSD= "<<(vm_mean/(vm_mean+vm_stdev)) <<" VM_difmax-vmsdoverVM___= "<<((*itmpsol).first-vm_stdev)/(((*itmpsol).first)/vm_mean)<<" VM___= "<<(((*itmpsol).first)/vm_mean)<<" 3*vm_stdev= "<<3*vm_stdev<<endl;
  double crit=vm_coef*(vm_mean+(vm_stdev*(vm_mean/(vm_mean+vm_stdev))));//(vm_mean/(vm_mean+vm_stdev))*vm_stdev;
  cout<<"critreion value= "<<crit<<endl;
std::vector<double*> lstpinsert;
double prop=0.001;
FILE *in=fopen(gnifinspointfile,"rt");
if (in==NULL) cout<<"file is not available"<<endl;
while (!feof(in))
{
  char chaine[500];
  char* res=fgets(chaine,500,in);
  if(res!=NULL)
  {
  double x,y,z;
  double q1,q2,q3;
  int nb=sscanf(chaine,"%le %le %le %le %le %le",&x,&y,&z,&q1,&q2,&q3);
  q1=prop*q1;
  q2=prop*q2;
  q3=prop*q3;
  if (nb!=-1 && nb!=6) cout<<"Wrong file format"<<endl;
  else if (nb==6)
  {
  double* pcoordbc=new double[6];
  pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z; pcoordbc[3]=q1; pcoordbc[4]=q2; pcoordbc[5]=q3;
  lstpinsert.push_back(pcoordbc);
  }
  }
}
fclose(in);
  
  
//octree initialization  
double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
TPL_MAP_ENTITE<FEM_NOEUD*> lstnoeud;
LISTE_FEM_NOEUD::iterator it;
for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
  {
  if (no->get_x()<xmin) xmin=no->get_x();  
  if (no->get_y()<ymin) ymin=no->get_y();  
  if (no->get_z()<zmin) zmin=no->get_z();  
  if (no->get_x()>xmax) xmax=no->get_x();  
  if (no->get_y()>ymax) ymax=no->get_y();  
  if (no->get_z()>zmax) zmax=no->get_z();
  lstnoeud.ajouter(no);
  }
 OT_VECTEUR_3D vecmin(xmin,ymin,zmin);OT_VECTEUR_3D vecmmax(xmax,ymax,zmax);
 OT_VECTEUR_3D vec(vecmmax,vecmin);
 
 OT_VECTEUR_3D min(xmin,ymin,zmin); OT_VECTEUR_3D max(xmax,ymax,zmax); OT_VECTEUR_3D average=(min+max)/2.; OT_VECTEUR_3D lengthvec(min,max);
 double length=sqrt(lengthvec*lengthvec); double bxr=1.1;
 xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
 xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);

TPL_OCTREE<FEM_NOEUD*,FEM_NOEUD*> octree;
octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
for (FEM_NOEUD* no=femmai->get_premier_noeud(it);no!=NULL;no=femmai->get_suivant_noeud(it))
  octree.inserer(no);
  
for (itmpsol=mpsol.begin();itmpsol!=mpsol.end();itmpsol++)
{
  if(crit<(*itmpsol).first)
  {
    FEM_NOEUD* pdefect=(*itmpsol).second;
    TPL_MAP_ENTITE<FEM_NOEUD*> lstneipdefect;
    int neindnb=lstneipdefect.get_nb();
    double search_radius1=search_radius;
    while(neindnb==0) 
    {
      octree.rechercher(pdefect->get_x(),pdefect->get_y(),pdefect->get_z(),search_radius1,lstneipdefect);
      neindnb=lstneipdefect.get_nb();
      search_radius1=search_radius1+0.1*search_radius;
    }
    TPL_MAP_ENTITE<FEM_NOEUD*>::ITERATEUR itfendnei;
    for(FEM_NOEUD* fendnei=lstneipdefect.get_premier(itfendnei);fendnei!=NULL;fendnei=lstneipdefect.get_suivant(itfendnei))
    {
//              std::map<double*,double*,std::less<double*> >::iterator itlstpinsert=lstpinsert.begin();//itlstpinsert!=lstpinsert.end();itlstpinsert++
      for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins;
        xyzpins=(*itlstpinsert);
         if(fendnei->get_x()==xyzpins[0] && fendnei->get_y()==xyzpins[1] && fendnei->get_z()==xyzpins[2])
         {
           lstpinsert.erase(itlstpinsert);
           itlstpinsert=lstpinsert.begin();
        }
      }
    }
  
  }
   
}

FILE* in1=fopen(out_gnifinspointfile_removondefect,"wt");
for(std::vector<double*>::iterator itlstpinsert=lstpinsert.begin();itlstpinsert!=lstpinsert.end();itlstpinsert++)
      {
        double* xyzpins1;
        xyzpins1=(*itlstpinsert);
        fprintf(in1,"%le %le %le %le %le %le\n",xyzpins1[0],xyzpins1[1],xyzpins1[2],xyzpins1[3]*1000.,xyzpins1[4]*1000.,xyzpins1[5]*1000.);
      }
  fclose(in1);
}
void CRIAQOPERATORS::surfmaker_e1(char* outputfilename,int meshno,char* magicfilename)
{
  MG_FILE gest(magicfilename);
  MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  //MG_FILE gest((char*)"longflxprtwithpockets_3d_3mm.magic");
  //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
  //MG_MAILLAGE* mai=gest.get_mg_maillageid(5674824);
  MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  gestnew->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  gestnew->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
    MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
    MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
    
    if(facevert->get_id()==76
     || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==13 || aretver->get_id()==110||aretver->get_id()==98
     || aretver->get_id()==81 || aretver->get_id()==83 || aretver->get_id()==10|| aretver->get_id()==12|| aretver->get_id()==90|| aretver->get_id()==97
     /* facevert->get_id()==133 
    || aretver->get_id()==355 || aretver->get_id()==362|| aretver->get_id()==374|| aretver->get_id()==298 || aretver->get_id()==291||aretver->get_id()==338
    || aretver->get_id()==326 || aretver->get_id()==333|| aretver->get_id()==319|| aretver->get_id()==312 || aretver->get_id()==305||aretver->get_id()==169
    || aretver->get_id()==162 || aretver->get_id()==141|| aretver->get_id()==174|| aretver->get_id()==148 || aretver->get_id()==155||aretver->get_id()==217
    || aretver->get_id()==212 || aretver->get_id()==191|| aretver->get_id()==184|| aretver->get_id()==64 || aretver->get_id()==280||aretver->get_id()==32
    || aretver->get_id()==387 || aretver->get_id()==198|| aretver->get_id()==205|| aretver->get_id()==269 || aretver->get_id()==262||aretver->get_id()==234
    || aretver->get_id()==227 || aretver->get_id()==255|| aretver->get_id()==248|| aretver->get_id()==241 || aretver->get_id()==380||aretver->get_id()==96
    || aretver->get_id()==369 || aretver->get_id()==344|| aretver->get_id()==128|| aretver->get_id()==274
    || sometver->get_id()==368|| sometver->get_id()==361|| sometver->get_id()==106|| sometver->get_id()==122|| sometver->get_id()==74|| sometver->get_id()==90
    || sometver->get_id()==233|| sometver->get_id()==240|| sometver->get_id()==247|| sometver->get_id()==254|| sometver->get_id()==261|| sometver->get_id()==268 
    || sometver->get_id()==226|| sometver->get_id()==224|| sometver->get_id()==190|| sometver->get_id()==197|| sometver->get_id()==204|| sometver->get_id()==211
    || sometver->get_id()==183|| sometver->get_id()==181|| sometver->get_id()==154|| sometver->get_id()==147|| sometver->get_id()==140|| sometver->get_id()==161
    || sometver->get_id()==168|| sometver->get_id()==138|| sometver->get_id()==318|| sometver->get_id()==311|| sometver->get_id()==304|| sometver->get_id()==297  
    || sometver->get_id()==325|| sometver->get_id()==332|| sometver->get_id()==290|| sometver->get_id()==288|| sometver->get_id()==42|| sometver->get_id()==58  
    || sometver->get_id()==352|| sometver->get_id()==354|| sometver->get_id()==10|| sometver->get_id()==26*/
      
    )
    {
      double x=vertices->get_x(); 
      double y=vertices->get_y(); 
      double z=vertices->get_z();
      //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
      MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
      mainew->ajouter_mg_noeud(newno);
      vertices->change_nouveau_numero(newno->get_id()); // what does this do?
    }
  }
  
  LISTE_MG_SEGMENT::iterator itSeg;
for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
      {
        //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
        //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
        //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
        MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
        MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
   if (faceseg->get_id()==76
     || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==13 || aretver->get_id()==110||aretver->get_id()==98
     //faceseg->get_id()==133 
    //|| aretver->get_id()==355 || aretver->get_id()==362|| aretver->get_id()==374|| aretver->get_id()==298 || aretver->get_id()==291||aretver->get_id()==338
    //|| aretver->get_id()==326 || aretver->get_id()==333|| aretver->get_id()==319|| aretver->get_id()==312 || aretver->get_id()==305||aretver->get_id()==169
    //|| aretver->get_id()==162 || aretver->get_id()==141|| aretver->get_id()==174|| aretver->get_id()==148 || aretver->get_id()==155||aretver->get_id()==217
    //|| aretver->get_id()==212 || aretver->get_id()==191|| aretver->get_id()==184|| aretver->get_id()==64 || aretver->get_id()==280||aretver->get_id()==32
    //|| aretver->get_id()==387 || aretver->get_id()==198|| aretver->get_id()==205|| aretver->get_id()==269 || aretver->get_id()==262||aretver->get_id()==234
   // || aretver->get_id()==227 || aretver->get_id()==255|| aretver->get_id()==248|| aretver->get_id()==241 || aretver->get_id()==380||aretver->get_id()==96
   // || aretver->get_id()==369 || aretver->get_id()==344|| aretver->get_id()==128|| aretver->get_id()==274
     
  )
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
  }
      }
      
LISTE_MG_TRIANGLE::iterator itTri;
for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
    {
      MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   if (facetri->get_id()==76)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
    MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
    MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
    mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
  }
    }
    
//gestnew->enregistrer("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump.magic");
gestnew->enregistrer(outputfilename);
//////////////////

////////////////////////
    
    
  /*
  mainew->supprimer_mg_triangleid(106625);
  mainew->supprimer_mg_triangleid(106662);
  mainew->supprimer_mg_triangleid(106626);
  mainew->supprimer_mg_triangleid(106701);
  mainew->supprimer_mg_triangleid(106702);
  mainew->supprimer_mg_triangleid(106639);*/
  //cout<<"output mesh no.= 495557 "<<endl;
  //gest.enregistrer("2dsurfmesh_cad.magic"); //CAD model output
   //scanned part output
  
  
}
void CRIAQOPERATORS::surfmaker_e2(char* outputfilename,int meshno,char* magicfilename)
{
  // surfmeshgenerator a partially scanned (deformed CAD model) for GNIF examples for simple exam with hole (badisplacenmentvector + smplwithpocket)
  MG_FILE gest(magicfilename);
  MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  //MG_FILE gest((char*)"2d_smplwithole_23mm.magic");
  //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
  //MG_MAILLAGE* mai=gest.get_mg_maillageid(1197);
  MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  gestnew->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  gestnew->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
    MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
    MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
    if(facevert->get_id()==76 
    || aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==98 || aretver->get_id()==13||aretver->get_id()==110
    || sometver->get_id()==10|| sometver->get_id()==12|| sometver->get_id()==81|| sometver->get_id()==83|| sometver->get_id()==90|| sometver->get_id()==97)
    {
      double x=vertices->get_x(); 
      double y=vertices->get_y(); 
      double z=vertices->get_z();
      //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
      MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
      mainew->ajouter_mg_noeud(newno);
      vertices->change_nouveau_numero(newno->get_id()); // what does this do?
    }
  }
  
  LISTE_MG_SEGMENT::iterator itSeg;
for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
      {
        //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
        //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
        //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
        MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
        //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
        MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
        //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
   if (faceseg->get_id()==76 ||
     aretver->get_id()==84 || aretver->get_id()==103|| aretver->get_id()==91|| aretver->get_id()==98 || aretver->get_id()==13|| aretver->get_id()==110)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
  }
      }
      
LISTE_MG_TRIANGLE::iterator itTri;
for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
    {
      MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   if (facetri->get_id()==76)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
    MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
    MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
    mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
  }
    }

  //gestnew->enregistrer("cadsurfmsh_smplwithole_23mm.magic");
  gestnew->enregistrer(outputfilename);
}
void CRIAQOPERATORS::surfmaker_e2t(char* outputfilename,int meshno,char* magicfilename)
{
  // surfmeshgenerator a partially scanned (deformed CAD model) for GNIF examples for simple exam with hole (badisplacenmentvector + smplwithpocket)
  MG_FILE gest(magicfilename);
  MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  //MG_FILE gest((char*)"2d_smplwithole_23mm.magic");
  //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
  //MG_MAILLAGE* mai=gest.get_mg_maillageid(1197);
  MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  gestnew->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  gestnew->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
    MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
    MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
    if(facevert->get_id()==37
      || aretver->get_id()==20 || aretver->get_id()==108|| aretver->get_id()==113|| aretver->get_id()==101 || aretver->get_id()==59||aretver->get_id()==64
      || aretver->get_id()==45 || aretver->get_id()==52|| aretver->get_id()==93|| aretver->get_id()==74 || aretver->get_id()==81||aretver->get_id()==88
      || sometver->get_id()==10|| sometver->get_id()==19|| sometver->get_id()==51|| sometver->get_id()==58|| sometver->get_id()==42|| sometver->get_id()==44
      || sometver->get_id()==87|| sometver->get_id()==73|| sometver->get_id()==80|| sometver->get_id()==71|| sometver->get_id()==100|| sometver->get_id()==107
//for e2t
    //|| aretver->get_id()==20 || aretver->get_id()==141|| aretver->get_id()==129|| aretver->get_id()==136 || aretver->get_id()==66||aretver->get_id()==59
    //|| aretver->get_id()==73 || aretver->get_id()==78|| aretver->get_id()==45|| aretver->get_id()==52 || aretver->get_id()==95||aretver->get_id()==102
    //|| aretver->get_id()==88 || aretver->get_id()==121|| aretver->get_id()==116|| aretver->get_id()==109 
    //|| sometver->get_id()==10|| sometver->get_id()==19|| sometver->get_id()==65|| sometver->get_id()==58|| sometver->get_id()==51|| sometver->get_id()==42
    //|| sometver->get_id()==44|| sometver->get_id()==72|| sometver->get_id()==115|| sometver->get_id()==87|| sometver->get_id()==108|| sometver->get_id()==101
    //|| sometver->get_id()==94|| sometver->get_id()==85|| sometver->get_id()==128|| sometver->get_id()==135
    )
    {
      double x=vertices->get_x(); 
      double y=vertices->get_y(); 
      double z=vertices->get_z();
      //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
      MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
      mainew->ajouter_mg_noeud(newno);
      vertices->change_nouveau_numero(newno->get_id()); // what does this do?
    }
  }
  
  LISTE_MG_SEGMENT::iterator itSeg;
for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
      {
        //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
        //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
        //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
        MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
        //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
        MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
        //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
   if (faceseg->get_id()==37
     || aretver->get_id()==20 || aretver->get_id()==108|| aretver->get_id()==113|| aretver->get_id()==101 || aretver->get_id()==59||aretver->get_id()==64
      || aretver->get_id()==45 || aretver->get_id()==52|| aretver->get_id()==93|| aretver->get_id()==74 || aretver->get_id()==81||aretver->get_id()==88
    //|| aretver->get_id()==20 || aretver->get_id()==141|| aretver->get_id()==129|| aretver->get_id()==136 || aretver->get_id()==66||aretver->get_id()==59
    //|| aretver->get_id()==73 || aretver->get_id()==78|| aretver->get_id()==45|| aretver->get_id()==52 || aretver->get_id()==95||aretver->get_id()==102
    //|| aretver->get_id()==88 || aretver->get_id()==121|| aretver->get_id()==116|| aretver->get_id()==109
     
  )
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
  }
      }
      
LISTE_MG_TRIANGLE::iterator itTri;
for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
    {
      MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   if (facetri->get_id()==37)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
    MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
    MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
    mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
  }
    }

  //gestnew->enregistrer("cadsurfmsh_smplwithole_23mm.magic");
  gestnew->enregistrer(outputfilename);
}
void CRIAQOPERATORS::surfmaker_e3(char* outputfilename,int meshno,char* magicfilename)
{
  //surfmesh maker (smplwithpocket+smpl withpocket)
  MG_FILE gest(magicfilename);
  MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
//MG_FILE gest((char*)"smplwithpocket_3d_13mm.magic");
  //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
  //MG_MAILLAGE* mai=gest.get_mg_maillageid(2982821);
  MG_GESTIONNAIRE* gestnew=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  gestnew->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  gestnew->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
    MG_ARETE* aretver=(MG_ARETE*) (vertices->get_lien_topologie());
    MG_SOMMET* sometver=(MG_SOMMET*) (vertices->get_lien_topologie());
    if(facevert->get_id()==5
    || aretver->get_id()==32 || aretver->get_id()==61|| aretver->get_id()==42|| aretver->get_id()==56 || aretver->get_id()==49||aretver->get_id()==27
    || aretver->get_id()==85 || aretver->get_id()==90|| aretver->get_id()==78|| aretver->get_id()==71 || aretver->get_id()==20||aretver->get_id()==13
    || sometver->get_id()==10|| sometver->get_id()==26|| sometver->get_id()==41|| sometver->get_id()==55|| sometver->get_id()==39|| sometver->get_id()==48
    || sometver->get_id()==84|| sometver->get_id()==77|| sometver->get_id()==70|| sometver->get_id()==68|| sometver->get_id()==12|| sometver->get_id()==19)
    {
      double x=vertices->get_x(); 
      double y=vertices->get_y(); 
      double z=vertices->get_z();
      //MG_NOEUD *newno=new MG_NOEUD(NULL,x,y,z,vertices->get_origine());
      MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
      mainew->ajouter_mg_noeud(newno);
      vertices->change_nouveau_numero(newno->get_id()); // what does this do?
    }
  }
  
  LISTE_MG_SEGMENT::iterator itSeg;
for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
      {
        //%275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010);
        //NOEUD linked with a topology (get_lien_topologie => $10 ) -> and %10=SOMMET(1,$9,0); %9=POINT_OCC(1);
        //%390=SEGMENT($13,$275,$277,1010); => %13=ARETE(3,$8,$14,$15,1,0); %275=NOEUD($10,93.94800850483597,-250.0000000000000,0.000000000000000,1010); %277=NOEUD($13,89.09047601136933,-249.4243468871926,0.000000000000000,1010);
        MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
        //cout<<"faceseg->get_id()= "<<faceseg->get_id()<<endl;
        MG_ARETE* aretver=(MG_ARETE*) (segment->get_lien_topologie());
        //cout<<"aretver->get_id()= "<<aretver->get_id()<<endl;
   if (faceseg->get_id()==5
    || aretver->get_id()==32 || aretver->get_id()==61|| aretver->get_id()==42|| aretver->get_id()==56 || aretver->get_id()==49||aretver->get_id()==27
    || aretver->get_id()==85 || aretver->get_id()==90|| aretver->get_id()==78|| aretver->get_id()==71 || aretver->get_id()==20||aretver->get_id()==13)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      //MG_SEGMENT *seg=new MG_SEGMENT(NULL,no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
  }
      }
      
LISTE_MG_TRIANGLE::iterator itTri;
for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
    {
      MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   if (facetri->get_id()==5)
   {
     MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
    MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
    MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
    mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
  }
    }

  //gestnew->enregistrer("scnsurf_smplwithpocket_3d_13mm_232mmdeform.magic");
  gestnew->enregistrer(outputfilename);
}
void CRIAQOPERATORS::gnifformatmaker(char* elementxtoutfile,char* nodtxtoutfile,int meshno,char* magicfilename)
{
  // gnifformatmaker
  MG_FILE gest(magicfilename);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
//MG_FILE gest((char*)"scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump.magic");
 // MG_MAILLAGE* mai=gest.get_mg_maillageid(2);
  //std::string fichiermaillage="CAD_output_elements.txt";
  //FILE* in=fopen(fichiermaillage.c_str(),"wt");
  //FILE* in1=fopen("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump_nodes.txt","wt");
  FILE* in1=fopen(nodtxtoutfile,"wt");
  //FILE* in1=fopen("scannedsurface_nodes.txt","wt"); 
  LISTE_MG_NOEUD::iterator itnod;
  unsigned long firstno=mai->get_premier_noeud(itnod)->get_id();
  firstno=firstno-1;
  cout<<firstno<<endl;
  int counter1=1;
  cout<<"mai->get_nb_mg_noeud()="<<mai->get_nb_mg_noeud()<<endl;
  for (MG_NOEUD* nod=mai->get_premier_noeud(itnod);nod!=NULL;nod=mai->get_suivant_noeud(itnod))
    {
      nod->change_nouveau_numero(counter1);
      //nod->change_id(counter1);
      //fprintf(in1,"%d node coordinate %le %le %le\n",nod->get_id(), nod->get_x(),nod->get_y(), nod->get_z() );
      //fprintf(in,"%d node coordinate %le %le %le\n",nod->get_id(), nod->get_x(),nod->get_y(), nod->get_z() );
      fprintf(in1,"%19d%32.9le%16.9le%8d\n",nod->get_nouveau_numero(), nod->get_x(),nod->get_y(),nod->get_nouveau_numero());
      fprintf(in1,"%7d%16.9le\n",nod->get_nouveau_numero(), nod->get_z());
      counter1++;
    }
  //FILE* in2=fopen("scnsurf_longflxprtwithpockets_3d_3mm_8mmdef_0bump_elements.txt","wt"); // CAD text output elements
  FILE* in2=fopen(elementxtoutfile,"wt");
  //FILE* in2=fopen("scannedsurface_elements.txt","wt"); //scanned text output elements
  int counter=1;
  LISTE_MG_TRIANGLE::iterator ittri;
  for (MG_TRIANGLE* tri=mai->get_premier_triangle(ittri);tri!=NULL;tri=mai->get_suivant_triangle(ittri))
    {
      fprintf(in2,"%10d        1%8d%8d%8d\n",counter,tri->get_noeud1()->get_nouveau_numero(),tri->get_noeud2()->get_nouveau_numero(),tri->get_noeud3()->get_nouveau_numero() );
      counter++;
    }
  fclose(in1);
  fclose(in2);
}
void CRIAQOPERATORS::rmovscnprtfromdefrmdcad(char* outputfilename,int meshno,char* magicfilename)
{
  MG_FILE gest(magicfilename);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  MG_GESTIONNAIRE* newgest=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  newgest->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  newgest->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mai->get_premier_noeud(itnod);vertices!=NULL;vertices=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
    //if(facevert->get_id()!=76 && topnd->get_dimension()!=3) //e1
    //if(facevert->get_id()!=5 && topnd->get_dimension()!=3) //e2
    {
      double x=vertices->get_x(); 
      double y=vertices->get_y(); 
      double z=vertices->get_z();
      MG_NOEUD *newno=new MG_NOEUD(vertices->get_lien_topologie(),x,y,z,vertices->get_origine());
      mainew->ajouter_mg_noeud(newno);
      vertices->change_nouveau_numero(newno->get_id()); // what does this do?
    }
  }
  
  LISTE_MG_SEGMENT::iterator itSeg;
  for (MG_SEGMENT * segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
  {
    MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
    MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
    //if (faceseg->get_id()!=76 && topseg->get_dimension()!=3)
     if(faceseg->get_id()!=5 && topseg->get_dimension()!=3) //e2 
    {
      MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
    }
  }
  
  LISTE_MG_TRIANGLE::iterator itTri;
  for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
  {
    MG_ELEMENT_TOPOLOGIQUE *toptri=triangle->get_lien_topologie();
    MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   // if (facetri->get_id()!=76 && toptri->get_dimension()!=3)
     if (facetri->get_id()!=5 && toptri->get_dimension()!=3) //e2
    {
      MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
      MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
      mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
    }
  }
newgest->enregistrer(outputfilename);
//gest.enregistrer("deformedcad_sansscannedpart.magic"); 
  
}
void CRIAQOPERATORS::msh2dmakerfrommsh3d(char* outputfilename,int meshno,char* magicfilename)
{
  //to make a 2D mesh out of 3D mesh :)
/*  
MG_MAILLAGE* mainew=mai->dupliquer(&gest);
LISTE_MG_TETRA::iterator ittet;
for(MG_TETRA* )
{
  
}

  cout<<"mainew->get_nb_mg_noeud()= "<<mainew->get_nb_mg_noeud()<<endl;
  cout<<"mainew->get_nb_mg_segment()= "<<mainew->get_nb_mg_segment()<<endl;
  cout<<"mainew->get_nb_mg_triangle()= "<<mainew->get_nb_mg_triangle()<<endl;
  cout<<"mainew->get_nb_mg_tetra()= "<<mainew->get_nb_mg_tetra()<<endl;
  
  LISTE_MG_TRIANGLE::iterator ittri;
  for (MG_TRIANGLE* triangle=mainew->get_premier_triangle(ittri);triangle!=NULL;triangle=mainew->get_suivant_triangle(ittri))
  {
    MG_ELEMENT_TOPOLOGIQUE *top=triangle->get_lien_topologie();
    int topdimension=top->get_dimension();
    if (topdimension==3)
      mainew->supprimer_mg_triangleid(triangle->get_id());
  }
  
  LISTE_MG_TETRA::iterator ittet;
  for (MG_TETRA* tetra=mainew->get_premier_tetra(ittet);tetra!=NULL;tetra=mainew->get_suivant_tetra(ittet))
    mainew->supprimer_mg_tetraid(tetra->get_id());
  //MG_FACE* face=geo->get_mg_faceid(5);
    

  LISTE_MG_SEGMENT::iterator itseg;
  for (MG_SEGMENT* segment=mainew->get_premier_segment(itseg);segment!=NULL;segment=mainew->get_suivant_segment(itseg))
  {
    MG_FACE* faceseg=(MG_FACE*) (segment->get_lien_topologie());
   if (faceseg->get_id()!=181)
     mainew->supprimer_mg_segmentid(segment->get_id());
   
    MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
    if(topseg->get_dimension()==3)
      cout<<"topseg->get_dimension()==3"<<endl;
  }
  
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* vertices=mainew->get_premier_noeud(itnod);vertices!=NULL;vertices=mainew->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=vertices->get_lien_topologie();
    MG_FACE* facevert=(MG_FACE*) (vertices->get_lien_topologie());
   //if (facevert->get_id()!=181)
    if (topnd->get_dimension()==0)
     mainew->supprimer_mg_segmentid(vertices->get_id());
  }
  
  for (MG_TRIANGLE* triangle=mainew->get_premier_triangle(ittri);triangle!=NULL;triangle=mainew->get_suivant_triangle(ittri))
  {
   MG_FACE* facetri=(MG_FACE*) (triangle->get_lien_topologie());
   if (facetri->get_id()!=181)
     mainew->supprimer_mg_triangleid(triangle->get_id());
  }
*/

  MG_FILE gest(magicfilename);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  MG_GESTIONNAIRE* newgest=new MG_GESTIONNAIRE();
  MG_GEOMETRIE* geonew=mai->get_mg_geometrie();
  newgest->ajouter_mg_geometrie(geonew);
  MG_MAILLAGE* mainew=new MG_MAILLAGE(geonew);
  newgest->ajouter_mg_maillage(mainew);
  
  LISTE_MG_NOEUD::iterator itNod;
  for (MG_NOEUD* nd=mai->get_premier_noeud(itNod);nd!=NULL;nd=mai->get_suivant_noeud(itNod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=nd->get_lien_topologie();
    if(topnd->get_dimension()<3)
    {
      double x=nd->get_x(); 
      double y=nd->get_y(); 
      double z=nd->get_z();
      MG_NOEUD *newno=new MG_NOEUD(nd->get_lien_topologie(),x,y,z,nd->get_origine());
      mainew->ajouter_mg_noeud(newno);
      nd->change_nouveau_numero(newno->get_id());
    }
  }
  LISTE_MG_SEGMENT::iterator itSeg;
  for (MG_SEGMENT* segment =mai->get_premier_segment(itSeg); segment; segment=mai->get_suivant_segment(itSeg))
  {
    MG_ELEMENT_TOPOLOGIQUE *topseg=segment->get_lien_topologie();
    if(topseg->get_dimension()<3)
    {
      MG_NOEUD* no1=mainew->get_mg_noeudid(segment->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(segment->get_noeud2()->get_nouveau_numero());
      MG_SEGMENT *seg=new MG_SEGMENT(segment->get_lien_topologie(),no1,no2,segment->get_origine());
      mainew->ajouter_mg_segment(seg);
    }
  }
  LISTE_MG_TRIANGLE::iterator itTri;
  for (MG_TRIANGLE * triangle =mai->get_premier_triangle(itTri); triangle; triangle=mai->get_suivant_triangle(itTri))      
  {
    MG_ELEMENT_TOPOLOGIQUE *toptri=triangle->get_lien_topologie();
    if(toptri->get_dimension()<3)
    {
      MG_NOEUD* no1=mainew->get_mg_noeudid(triangle->get_noeud1()->get_nouveau_numero());
      MG_NOEUD* no2=mainew->get_mg_noeudid(triangle->get_noeud2()->get_nouveau_numero());
      MG_NOEUD* no3=mainew->get_mg_noeudid(triangle->get_noeud3()->get_nouveau_numero());
      mainew->ajouter_mg_triangle(triangle->get_lien_topologie(),no1,no2,no3,triangle->get_origine());
    }
  }
  newgest->enregistrer(outputfilename);
}
void CRIAQOPERATORS::bumpmaker(char* outputfilename,int bumpndid, double bumptip,int meshno,char* magicfilename)
{
  MG_FILE gest(magicfilename);
  //MG_GEOMETRIE* geo=gest.get_mg_geometrie(0);
  MG_MAILLAGE* mai;
  if (meshno==0) mai=gest.get_mg_maillage(meshno); else mai=gest.get_mg_maillageid(meshno);
  //MG_FILE gest((char*)"scnsurf_smplwithpocket_3d_13mm_64mmdeform.magic");
  //MG_GEOMETRIE* geo=gest.get_mg_geometrieid(1);
  //MG_MAILLAGE* mai=gest.get_mg_maillageid(2);
  MG_NOEUD* bptipnd;
  TPL_MAP_ENTITE<MG_NOEUD*> usednd;
  LISTE_MG_NOEUD::iterator itnod;
  for (MG_NOEUD* nd=mai->get_premier_noeud(itnod);nd!=NULL;nd=mai->get_suivant_noeud(itnod))
  {
    MG_ELEMENT_TOPOLOGIQUE *topnd=nd->get_lien_topologie();
    if (topnd->get_dimension() <3 && nd->get_id()==bumpndid)
    {
      bptipnd=nd;
      double new_xyz[3];
      new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+bumptip;
      nd->change_coord(new_xyz);
      usednd.ajouter(nd);
    }
  }
  TPL_MAP_ENTITE<MG_NOEUD*> usedndnd;
  TPL_MAP_ENTITE<MG_NOEUD*> neignd;
  for(int ineignd=0;ineignd<bptipnd->get_lien_segment()->get_nb();ineignd++)
  {
    MG_SEGMENT* segnei=bptipnd->get_lien_segment()->get(ineignd);
    if(segnei->get_noeud1()!=bptipnd) neignd.ajouter(segnei->get_noeud1());
    else if(segnei->get_noeud2()!=bptipnd) neignd.ajouter(segnei->get_noeud2());
  }
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd;
  for (MG_NOEUD* nd=neignd.get_premier(itneignd);nd!=NULL;nd=neignd.get_suivant(itneignd))
  {
    double new_xyz[3];
    new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*3./4.);
    nd->change_coord(new_xyz);
    usednd.ajouter(nd);
    usedndnd.ajouter(nd);
  }
  
  TPL_MAP_ENTITE<MG_NOEUD*> neignd1;
  TPL_MAP_ENTITE<MG_NOEUD*> usedndnd1;
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusednd;
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusedndnd;
  for(MG_NOEUD* nd=usedndnd.get_premier(itusedndnd);nd!=NULL;nd=usedndnd.get_suivant(itusedndnd))
  {
    for(int ineignd=0;ineignd<nd->get_lien_segment()->get_nb();ineignd++)
    {
      MG_SEGMENT* segnei=nd->get_lien_segment()->get(ineignd);
      int sgnd1=0;
      int sgnd2=0;
      for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
      {
        if(segnei->get_noeud1()==nd1) sgnd1++;
      }
      for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
      {
        if(segnei->get_noeud2()==nd1) sgnd2++;
      }
      if(sgnd1==0) neignd1.ajouter(segnei->get_noeud1());
      else if(sgnd2==0) neignd1.ajouter(segnei->get_noeud2());
    }
  }
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd1;
  for (MG_NOEUD* nd=neignd1.get_premier(itneignd1);nd!=NULL;nd=neignd1.get_suivant(itneignd1))
  {
    double new_xyz[3];
    new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*2./4.);
    nd->change_coord(new_xyz);
    usednd.ajouter(nd);
    usedndnd1.ajouter(nd);
  }
  
  TPL_MAP_ENTITE<MG_NOEUD*> neignd2;
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itusedndnd1;
  for(MG_NOEUD* nd=usedndnd1.get_premier(itusedndnd1);nd!=NULL;nd=usedndnd1.get_suivant(itusedndnd1))
  {
    for(int ineignd=0;ineignd<nd->get_lien_segment()->get_nb();ineignd++)
    {
      MG_SEGMENT* segnei=nd->get_lien_segment()->get(ineignd);
      int sgnd1=0;
      int sgnd2=0;
      for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
      {
        if(segnei->get_noeud1()==nd1) sgnd1++;
      }
      for(MG_NOEUD* nd1=usednd.get_premier(itusednd);nd1!=NULL;nd1=usednd.get_suivant(itusednd))
      {
        if(segnei->get_noeud2()==nd1) sgnd2++;
      }
      if(sgnd1==0) neignd2.ajouter(segnei->get_noeud1());
      else if(sgnd2==0) neignd2.ajouter(segnei->get_noeud2());
    }
  }
  TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR itneignd2;
  for (MG_NOEUD* nd=neignd2.get_premier(itneignd2);nd!=NULL;nd=neignd2.get_suivant(itneignd2))
  {
    double new_xyz[3];
    new_xyz[0]=nd->get_x(); new_xyz[1]=nd->get_y(); new_xyz[2]=nd->get_z()+(bumptip*1./4.);
    nd->change_coord(new_xyz);
  }
  
  
  //gest.enregistrer("scnsurf_smplwithpocket_3d_13mm_64mmdeform_withbump5mm.magic");
  gest.enregistrer(outputfilename);
}
int CRIAQOPERATORS::import_triangulation_gnif(char* outputfilename,char* triangulationfile)
{
FILE* in=fopen(triangulationfile,"rt");
    if (in==NULL) return 0;


    MG_GESTIONNAIRE gesttmp;
    MG_MAILLAGE* mai=new MG_MAILLAGE(NULL);
    gesttmp.ajouter_mg_maillage(mai);
    char mess[500];

    double xmax=-1e308;
    double xmin=1e308;
    double ymax=-1e308;
    double ymin=1e308;
    double zmax=-1e308;
    double zmin=1e308;
    while (!feof(in))
    {
        double x,y,z;
        char nom[20];
        char *res=fgets(mess,500,in);
        if (feof(in)) break;
        if(res!=NULL)
        {
        sscanf(mess,"%le %le %le",&x,&y,&z);
        if (x<xmin) xmin=x;
        if (x>xmax) xmax=x;
        if (y<ymin) ymin=y;
        if (y>ymax) ymax=y;
        if (z<zmin) zmin=z;
        if (z>zmax) zmax=z;
        if (feof(in)) return 1;
        MG_NOEUD* noeud1=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
        noeud1->change_nouveau_numero(-1);
        res=fgets(mess,500,in);
        if (feof(in)) break;
        sscanf(mess,"%le %le %le",&x,&y,&z);
        if (x<xmin) xmin=x;
        if (x>xmax) xmax=x;
        if (y<ymin) ymin=y;
        if (y>ymax) ymax=y;
        if (z<zmin) zmin=z;
        if (z>zmax) zmax=z;
        if (feof(in)) return 1;
        res=fgets(mess,500,in);
        if (feof(in)) break;
        MG_NOEUD* noeud2=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
        noeud2->change_nouveau_numero(-1);
        sscanf(mess,"%le %le %le",&x,&y,&z);
        if (x<xmin) xmin=x;
        if (x>xmax) xmax=x;
        if (y<ymin) ymin=y;
        if (y>ymax) ymax=y;
        if (z<zmin) zmin=z;
        if (z>zmax) zmax=z;
        if (feof(in)) return 1;
        MG_NOEUD* noeud3=mai->ajouter_mg_noeud(NULL,x,y,z,TRIANGULATION);
        noeud3->change_nouveau_numero(-1);
        mai->ajouter_mg_triangle(NULL,noeud1,noeud2,noeud3,TRIANGULATION);
        }
    }
    fclose(in);
    //to merge the nodes of each triangulation which are close (the same)
    double eps=1e-4;
    BOITE_3D boite(xmin,ymin,zmin,xmax,ymax,zmax);
    double rayon=boite.get_rayon()*1.1;
    double x=boite.get_xcentre();
    double y=boite.get_ycentre();
    double z=boite.get_zcentre();
    boite.reinit(x-rayon,y-rayon,z-rayon,x+rayon,y+rayon,z+rayon);
    TPL_GRILLE<MG_NOEUD*> grille;
    int nb_noeud=mai->get_nb_mg_noeud();
    int pas=(int)(pow(nb_noeud,0.33333333333333333))+1;
    grille.initialiser(boite.get_xmin(),boite.get_ymin(),boite.get_zmin(),boite.get_xmax(),boite.get_ymax(),boite.get_zmax(),pas,pas,pas);
    std::vector<MG_NOEUD*> lst_noeud;
    LISTE_MG_NOEUD::iterator  it;
    int i=0;
    for (MG_NOEUD* noeud=mai->get_premier_noeud(it);noeud!=NULL;noeud=mai->get_suivant_noeud(it))
    {
        //MG_NOEUD* noeud=mai->get_mg_noeud(i);
        grille.inserer(noeud);
        noeud->change_nouveau_numero(i);
        lst_noeud.insert(lst_noeud.end(),noeud);
        i++;
    }

    int nb_cell=grille.get_nb_cellule();
    for (int i=0;i<nb_cell;i++)
    {
        TPL_CELLULE_GRILLE<MG_NOEUD*> *cell=grille.get_cellule(i);
        int nb_noeud_cell=cell->get_nb_entite();
        for (int j=0;j<nb_noeud_cell;j++)
        {
            MG_NOEUD* noeud1=cell->get_entite(j);
//                  if (noeud1->get_nouveau_numero()!=(-1)) continue;
//                  noeud1->change_nouveau_numero(noeud1->get_id());
            double *xyz1=noeud1->get_coord();
            int nb=noeud1->get_lien_segment()->get_nb();
            double longref=0.;
            for (int jj=0;jj<nb;jj++)
              longref=longref+noeud1->get_lien_segment()->get(jj)->get_longueur();
            longref=longref/nb;
            for (int k=j+1;k<nb_noeud_cell;k++)
            {
                MG_NOEUD* noeud2=cell->get_entite(k);
                double *xyz2=noeud2->get_coord();
                OT_VECTEUR_3D vec(xyz1,xyz2);
                double longueur=vec.get_longueur();
                if (longueur<eps*longref)
                {
                    lst_noeud[noeud2->get_nouveau_numero()]=noeud1;
                }
            }
        }
    }
    MG_GESTIONNAIRE gest;
    MG_MAILLAGE* mai2=new MG_MAILLAGE(NULL);
    gest.ajouter_mg_maillage(mai2);
    for (int i=0;i<nb_noeud;i++)
    {
        MG_NOEUD* noeud=lst_noeud[i];
        MG_NOEUD* nvnoeud;
        if (noeud->get_nouveau_numero()==i) nvnoeud=mai2->ajouter_mg_noeud(NULL,noeud->get_x(),noeud->get_y(),noeud->get_z(),TRIANGULATION);
        else nvnoeud=lst_noeud[noeud->get_nouveau_numero()];
        lst_noeud[i]=nvnoeud;
    }
    int nb_tri=mai->get_nb_mg_triangle();
    for (int i=0;i<nb_tri;i++)
    {
        MG_TRIANGLE* tri=mai->get_mg_triangle(i);
        MG_NOEUD* noeud1=tri->get_noeud1();
        MG_NOEUD* noeud2=tri->get_noeud2();
        MG_NOEUD* noeud3=tri->get_noeud3();
        mai2->ajouter_mg_triangle(NULL,lst_noeud[noeud1->get_nouveau_numero()],lst_noeud[noeud2->get_nouveau_numero()],lst_noeud[noeud3->get_nouveau_numero()],TRIANGULATION);
    }
    gest.enregistrer(outputfilename);
    return 1;
}
std::vector<double*> CRIAQOPERATORS::sp_project_onCAD(char* magicfilename,int nummai,char* samplepoints1)
{
  std::vector<double*> lstsmplpoints;
std::vector<double*> lstsmplpoints_afterprojectioncad; 
FILE *in=fopen(samplepoints1,"rt");
if (in==NULL) cout<<"file is not available"<<endl;//affiche((char*) "file is not available");
while(!feof(in))
{
  char chaine[500];
  char* res=fgets(chaine,500,in);
  if(res!=NULL)
  {
    double x,y,z;
    int nb=sscanf(chaine,"%le %le %le",&x,&y,&z);
    double* pcoordbc=new double[3];
    if (nb!=-1 && nb!=3)
      cout<<"Wrong file format"<<endl;
    if(nb==3)
    {
      pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z;
      lstsmplpoints.push_back(pcoordbc);
    }
  }
}
//cout<<"lstsmplpoints.size()= "<<lstsmplpoints.size()<<endl;
for(std::vector<double*>::iterator ilstsmplpoints=lstsmplpoints.begin();ilstsmplpoints!=lstsmplpoints.end();ilstsmplpoints++)
{
double* p=(*ilstsmplpoints);
///cout<<"p[0]= "<<p[0]<<" p[0]= "<<p[1]<<" p[0]= "<<p[2]<<endl;
MG_FILE gest(magicfilename);
MG_MAILLAGE* mai=gest.get_mg_maillageid(nummai);
int geonb=gest.get_nb_mg_geometrie();
int nogeo=0;
MG_GEOMETRIE* geo;
if(geonb>0)
  geo=gest.get_mg_geometrie(0);
else if(geonb<1)
{
  nogeo++;
  cout<<"no geometry is found !!!"<<endl;
}
  
double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
LISTE_MG_NOEUD::iterator it;
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
{
 if (no->get_x()<xmin) xmin=no->get_x();  
 if (no->get_y()<ymin) ymin=no->get_y();  
 if (no->get_z()<zmin) zmin=no->get_z();  
 if (no->get_x()>xmax) xmax=no->get_x();  
 if (no->get_y()>ymax) ymax=no->get_y();  
 if (no->get_z()>zmax) zmax=no->get_z();
 lstnoeud.ajouter(no); 
}
OT_VECTEUR_3D vecmin(xmin,ymin,zmin);
OT_VECTEUR_3D vecmax(xmax,ymax,zmax);
OT_VECTEUR_3D lengthvec(vecmin,vecmax);
double search_radius=0.0000001;
//double search_radius=0.001*(lengthvec.get_longueur());
//char bboxdiag[1000];
//sprintf(bboxdiag,"the bounding box diagonal is: %lf , and the neighborhood search radius is: %lf ",vec.get_longueur(),search_radius);
//affiche((char*) bboxdiag);
OT_VECTEUR_3D average=(vecmin+vecmax)/2.;
double length=sqrt(lengthvec*lengthvec);
double bxr=1.1;
xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);
TPL_OCTREE<MG_TRIANGLE*,MG_NOEUD*> octree;
octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
LISTE_MG_TRIANGLE::iterator it9;
for (MG_TRIANGLE* tri=mai->get_premier_triangle(it9);tri!=NULL;tri=mai->get_suivant_triangle(it9))
  octree.inserer(tri);
TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreends;
octreends.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
  octreends.inserer(no);

// determine the normal of geometry in the inserting point coordination and saving in "normale"
double normale[3];
std::multimap<double,MG_TRIANGLE*,std::less<double> > delinsqualtri; 
int insideface=0;
std::multimap<double,MG_FACE*,std::less<double> > lstfacepointdist;
int facenb=geo->get_nb_mg_face();
            LISTE_MG_FACE::iterator itf;
            double xyz[3];
              xyz[0]=p[0];
              xyz[1]=p[1];
              xyz[2]=p[2];
            for(MG_FACE* fac=geo->get_premier_face(itf);fac!=NULL;fac=geo->get_suivant_face(itf))
            {
              double uv[2];
              fac->inverser(uv,xyz);
              double xyz2[3];
              fac->evaluer(uv,xyz2);
              OT_VECTEUR_3D vechk(xyz,xyz2);
              lstfacepointdist.insert(pair<double,MG_FACE*> (vechk.get_longueur(),fac));
              //double dis;
              //MG_GEOMETRIE_OUTILS f1;
              //int ch=f1.calcule_distance_contour_face_xyz(xyz,fac,&dis);
              ////cout<<"xyz= "<<xyz[0]<<" xyz= "<<xyz[1]<<" xyz= "<<xyz[2]<<"  xyz2= "<<xyz2[0]<<" xyz2= "<<xyz2[1]<<" xyz2= "<<xyz2[2]<<endl;
              //cout<<"dis= "<<dis<<"   vechk.get_longueur()= "<<vechk.get_longueur()<<endl;
              //if (ch==1 && dis>=0.)
             // if (vechk.get_longueur()<5e-3)
            }
std::multimap<double,MG_FACE*,std::less<double> >::iterator itlstfacepointdist=lstfacepointdist.begin();
if (lstfacepointdist.size()>0)      
{
  //cout<<"lstfacepointdist.size()= "<<lstfacepointdist.size()<<endl;
  //cout<<"(*itlstfacepointdist).first= "<<(*itlstfacepointdist).first<<endl;
  MG_FACE* face=(*itlstfacepointdist).second;
  double uv[2];
  face->inverser(uv,xyz);
  face->calcul_normale_unitaire(uv,normale);
  MG_COFACE *coface=face->get_mg_coface(0);
  int signe=coface->get_orientation();
  
  TPL_MAP_ENTITE<MG_TRIANGLE*> lstneitri;
  int neitrinb=lstneitri.get_nb();
  while(neitrinb==0) 
{
  octree.rechercher(p[0],p[1],p[2],search_radius,lstneitri);
  neitrinb=lstneitri.get_nb();
  search_radius=search_radius+0.01*search_radius;
}
TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR itnormtst;
OT_VECTEUR_3D  normale_discrete(0.,0.,0.);
for (MG_TRIANGLE* tri1=lstneitri.get_premier(itnormtst);tri1!=NULL;tri1=lstneitri.get_suivant(itnormtst))
    {
      OT_VECTEUR_3D vec1(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
      OT_VECTEUR_3D vec2(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
      OT_VECTEUR_3D normtst=(vec1&vec2); normtst.norme();
      normale_discrete=normtst+normale_discrete;
    }
normale_discrete.norme();
double n_ps=normale_discrete*normale;
n_ps=n_ps*-1;
if (n_ps>1.) n_ps=1.;
if (n_ps<-1.) n_ps=-1.;
double n_angle=acos(n_ps);
//cout<<"n_angle= "<<n_angle<<endl;
if(n_angle<(0.5*3.14159265358979323846) || n_angle>-(0.5*3.14159265358979323846))             
 insideface++;
}

//cout<<"insideface= "<<insideface<<endl;
if(insideface<1)
  cout<<"insideface is not found !!!"<<endl;
if(insideface>0)
{           
TPL_MAP_ENTITE<MG_TRIANGLE*> lstneitri1;
int neitrinb=lstneitri1.get_nb();
double search_radius1=0.01;
while(neitrinb==0) 
{
  octree.rechercher(p[0],p[1],p[2],500.*search_radius1,lstneitri1);//sasan
  neitrinb=lstneitri1.get_nb();
  search_radius1=search_radius1+0.00001*search_radius1;
}
//cout<<"lstneitri1.get_nb()= "<<lstneitri1.get_nb()<<endl;
double pproj[3];
TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR inspi;
int tri_found=0;
int beshmar=0;


for (MG_TRIANGLE* insphertri=lstneitri1.get_premier(inspi);insphertri!=NULL;insphertri=lstneitri1.get_suivant(inspi))
{
  //cout<<insphertri->get_id()<<endl;
  //while (tri_found<1 || beshmar<lstneitri.get_nb())
{
  beshmar++;
    MG_NOEUD* inspnd1=insphertri->get_noeud1();
    MG_NOEUD* inspnd2=insphertri->get_noeud2();
    MG_NOEUD* inspnd3=insphertri->get_noeud3();
    OT_VECTEUR_3D tvec1(inspnd1->get_coord(),inspnd2->get_coord());
    OT_VECTEUR_3D tvec2(inspnd1->get_coord(),inspnd3->get_coord());
    OT_VECTEUR_3D tnorm=(tvec1&tvec2); tnorm.norme();
    OT_VECTEUR_3D norm=normale;

    double parallel=norm*tnorm;
    //cout<<"parallel= "<<parallel<<endl;
    if (parallel>0.7 || parallel<-0.7)
    {
      
      //cout<<"itis"<<endl;
    OT_VECTEUR_3D vvec(inspnd1->get_coord(),p);
  double t= (tnorm*vvec)/parallel;
  t=-t;
  pproj[0]=p[0]+t*norm.get_x();
  pproj[1]=p[1]+t*norm.get_y();
  pproj[2]=p[2]+t*norm.get_z();
//////////////////////////////////////////////////////////////
//check if the projected p is inside the triangle !!!!!  
  MG_MAILLAGE_OUTILS tribas;
  int tribasval=tribas.estdansletriangle(insphertri,pproj[0],pproj[1],pproj[2]);
  int intri=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::INTERIEUR);
  int insidtri=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::STRICTINTERIEUR);
  int onedge=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::SUR_ARETE);
  //printf("projection before is: %16.14lf %16.14lf %16.14lf \n",pproj[0],pproj[1],pproj[2]);
  //cout<<"intri= "<<intri<<endl;
if(intri==1)
{
  //cout<<"normale= "<<norm.get_x()<<" normale= "<<norm.get_y()<<" normale= "<<norm.get_z()<<endl;
  //printf("projection is: %16.14lf %16.14lf %16.14lf \n",pproj[0],pproj[1],pproj[2]);
  double* proj=new double[3];
  proj[0]=pproj[0]; proj[1]=pproj[1]; proj[2]=pproj[2];
  OT_VECTEUR_3D disprojp(p,pproj);
    cout<<"CAD distance proj-p= "<<disprojp.get_longueur()<<endl;
  lstsmplpoints_afterprojectioncad.push_back(proj);
  tri_found++;
  break;
  
}

}
}
}
}
}

 //for(std::vector<double*>::iterator itsprojcad=lstsmplpoints_afterprojectioncad.begin();itsprojcad!=lstsmplpoints_afterprojectioncad.end();itsprojcad++)
  {
    //double* xyz_cadproj=(*itsprojcad);
    //cout<<"xyz_cadproj[0]= "<<xyz_cadproj[0]<<" xyz_cadproj[1]= "<<xyz_cadproj[1]<<" xyz_cadproj[2]= "<<xyz_cadproj[2]<<endl;
    
  }
return lstsmplpoints_afterprojectioncad;
}
std::vector<double*> CRIAQOPERATORS::sp_project_onSCAN(char* magicfilename,int nummai,char* samplepoints1)
{
  std::vector<double*> lstsmplpoints;
  std::vector<double*> lstsmplpoints_afterprojection; 
FILE *in=fopen(samplepoints1,"rt");
if (in==NULL) cout<<"file is not available"<<endl;//affiche((char*) "file is not available");
while(!feof(in))
{
  char chaine[500];
  char* res=fgets(chaine,500,in);
  if(res!=NULL)
  {
    double x,y,z;
    int nb=sscanf(chaine,"%le %le %le",&x,&y,&z);
    double* pcoordbc=new double[3];
    if (nb!=-1 && nb!=3)
      cout<<"Wrong file format"<<endl;
    if(nb==3)
    {
      pcoordbc[0]=x; pcoordbc[1]=y; pcoordbc[2]=z;
      lstsmplpoints.push_back(pcoordbc);
    }
  }
}
cout<<"lstsmplpoints.size()= "<<lstsmplpoints.size()<<endl;
for(std::vector<double*>::iterator ilstsmplpoints=lstsmplpoints.begin();ilstsmplpoints!=lstsmplpoints.end();ilstsmplpoints++)
{
  ///cout<<"<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"<<endl;
double* p=(*ilstsmplpoints);
///cout<<"p[0]= "<<p[0]<<" p[0]= "<<p[1]<<" p[0]= "<<p[2]<<endl;
MG_FILE gest(magicfilename);
MG_MAILLAGE* mai=gest.get_mg_maillageid(nummai);
  
double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
LISTE_MG_NOEUD::iterator it;
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
{
 if (no->get_x()<xmin) xmin=no->get_x();  
 if (no->get_y()<ymin) ymin=no->get_y();  
 if (no->get_z()<zmin) zmin=no->get_z();  
 if (no->get_x()>xmax) xmax=no->get_x();  
 if (no->get_y()>ymax) ymax=no->get_y();  
 if (no->get_z()>zmax) zmax=no->get_z();
 lstnoeud.ajouter(no); 
}
OT_VECTEUR_3D vecmin(xmin,ymin,zmin);
OT_VECTEUR_3D vecmax(xmax,ymax,zmax);
OT_VECTEUR_3D lengthvec(vecmin,vecmax);
double search_radius=5.1;//0.0000001;
//double search_radius=0.001*(lengthvec.get_longueur());
//char bboxdiag[1000];
//sprintf(bboxdiag,"the bounding box diagonal is: %lf , and the neighborhood search radius is: %lf ",vec.get_longueur(),search_radius);
//affiche((char*) bboxdiag);
OT_VECTEUR_3D average=(vecmin+vecmax)/2.;
double length=sqrt(lengthvec*lengthvec);
double bxr=30.1; //this can be changed sasan
xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);
TPL_OCTREE<MG_TRIANGLE*,MG_NOEUD*> octree;
octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
LISTE_MG_TRIANGLE::iterator it9;
for (MG_TRIANGLE* tri=mai->get_premier_triangle(it9);tri!=NULL;tri=mai->get_suivant_triangle(it9))
  octree.inserer(tri);
TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreends;
octreends.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
  octreends.inserer(no);

// determine the normal of geometry in the inserting point coordination and saving in "normale"
double normale[3];
//////////////////////////////// if no geometry

  TPL_MAP_ENTITE<MG_TRIANGLE*> lstneitri;
  int neitrinb=lstneitri.get_nb();
  while(neitrinb==0) 
{
  octree.rechercher(p[0],p[1],p[2],search_radius,lstneitri);
  neitrinb=lstneitri.get_nb();
  search_radius=search_radius+0.00001*search_radius;
}
//cout<<"lstneitri.get_nb()= "<<lstneitri.get_nb()<<endl;
TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR itnormtst;
OT_VECTEUR_3D  normale_discrete(0.,0.,0.);
for (MG_TRIANGLE* tri1=lstneitri.get_premier(itnormtst);tri1!=NULL;tri1=lstneitri.get_suivant(itnormtst))
    {
      OT_VECTEUR_3D vec1(tri1->get_noeud1()->get_coord(),tri1->get_noeud2()->get_coord());
      OT_VECTEUR_3D vec2(tri1->get_noeud1()->get_coord(),tri1->get_noeud3()->get_coord());
      OT_VECTEUR_3D normtst=(vec1&vec2); normtst.norme();
      //cout<<"normtst.get_x()= "<<normtst.get_x()<<" normtst.get_y()= "<<normtst.get_y()<<" normtst.get_z()= "<<normtst.get_z()<<endl;
      normale_discrete=normtst+normale_discrete;
    }
normale_discrete.norme();
normale[0]=normale_discrete.get_x();    normale[1]=normale_discrete.get_y();    normale[2]=normale_discrete.get_z();
///cout<<"normale[0]= "<<normale[0]<<" normale[1]= "<<normale[1]<<" normale[2]= "<<normale[2]<<endl;
////////////////////////////////
    
TPL_MAP_ENTITE<MG_TRIANGLE*> lstneitri1;
int neitrinb1=lstneitri1.get_nb();
double search_radius1=2.;//sasan
while(neitrinb1==0) 
{
  octree.rechercher(p[0],p[1],p[2],search_radius1,lstneitri1);
  neitrinb1=lstneitri1.get_nb();
  search_radius1=search_radius1+0.00001*search_radius1;
}
//cout<<"lstneitri1.get_nb()= "<<lstneitri1.get_nb()<<endl;
double pproj[3];
TPL_MAP_ENTITE<MG_TRIANGLE*>::ITERATEUR inspi;
int tri_found=0;
int beshmar=0;


for (MG_TRIANGLE* insphertri=lstneitri1.get_premier(inspi);insphertri!=NULL;insphertri=lstneitri1.get_suivant(inspi))
{
  //while (tri_found<1 || beshmar<lstneitri.get_nb())

  beshmar++;
    MG_NOEUD* inspnd1=insphertri->get_noeud1();
    MG_NOEUD* inspnd2=insphertri->get_noeud2();
    MG_NOEUD* inspnd3=insphertri->get_noeud3();
    OT_VECTEUR_3D tvec1(inspnd1->get_coord(),inspnd2->get_coord());
    OT_VECTEUR_3D tvec2(inspnd1->get_coord(),inspnd3->get_coord());
    OT_VECTEUR_3D tnorm=(tvec1&tvec2); tnorm.norme();
    OT_VECTEUR_3D norm=normale;

    double parallel=norm*tnorm;
    //cout<<"parallel= "<<parallel<<endl;
    if (parallel>0.7 || parallel<-0.7)
    {
      //cout<<"itis"<<endl;
    OT_VECTEUR_3D vvec(inspnd1->get_coord(),p);
  double t= (tnorm*vvec)/parallel;
  t=-t;
  pproj[0]=p[0]+t*norm.get_x();
  pproj[1]=p[1]+t*norm.get_y();
  pproj[2]=p[2]+t*norm.get_z();
//////////////////////////////////////////////////////////////
//check if the projected p is inside the triangle !!!!!  
  MG_MAILLAGE_OUTILS tribas;
  int tribasval=tribas.estdansletriangle(insphertri,pproj[0],pproj[1],pproj[2]);
  int intri=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::INTERIEUR);
  int insidtri=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::STRICTINTERIEUR);
  int onedge=tribas.compare_etat_triangle(tribasval,MG_MAILLAGE_OUTILS::SUR_ARETE);
  //printf("projection before is: %16.14lf %16.14lf %16.14lf \n",pproj[0],pproj[1],pproj[2]);
  //cout<<"intri= "<<intri<<endl;
  if(intri==1)
  {
    //printf("projection is: %16.14lf %16.14lf %16.14lf \n",pproj[0],pproj[1],pproj[2]);
    OT_VECTEUR_3D disprojp(p,pproj);
    cout<<"SCAN distance proj-p= "<<disprojp.get_longueur()<<endl;
  double* proj=new double[3];
  proj[0]=pproj[0]; proj[1]=pproj[1]; proj[2]=pproj[2];  
  lstsmplpoints_afterprojection.push_back(proj);
  tri_found++;
  //break;
  }
   } 

}
}
return lstsmplpoints_afterprojection;
}
void CRIAQOPERATORS::sp_project_onCAD_SCAN(char* magicfilenamecad,int nummaicad,char* samplepointscad,char* magicfilenamescn,int nummaiscn,char* samplepointsscn,char* outputspbcfilename)
{
  std::vector<double*> sprojcad;
  sprojcad=sp_project_onCAD(magicfilenamecad,nummaicad,samplepointscad);
  std::vector<double*> sprojscn;
  sprojscn=sp_project_onSCAN(magicfilenamescn,nummaiscn,samplepointsscn);
  

  int sizesm;
  if(sprojcad.size()==sprojscn.size())
    sizesm=sprojcad.size();
  else
    cout<<"not sprojcad.size()==sprojscn.size() !!!!!!!!!!!"<<endl;
  cout<<"sprojcad.size()= "<<sprojcad.size()<<endl;
  cout<<"sprojscn.size()= "<<sprojscn.size()<<endl;
  FILE* insmp=fopen(outputspbcfilename,"wt");/*
  for(std::vector<double*>::iterator itsprojcad=sprojcad.begin();itsprojcad!=sprojcad.end();itsprojcad++)
  {
    double* xyz_cadproj=(*itsprojcad);
    cout<<"xyz_cadproj[0]= "<<xyz_cadproj[0]<<endl;
    
  }*/
  std::vector<double*>::iterator itsprojcad=sprojcad.begin();
  std::vector<double*>::iterator itsprojscn=sprojscn.begin();
  
  for(int iitt=0;iitt<sizesm;iitt++)
  {
    if(iitt>0)
    {
      itsprojcad++;
      itsprojscn++;
    }
    double* xyz_cadproj=(*itsprojcad);
    double* xyz_scnproj=(*itsprojscn);
    fprintf(insmp,"%le %le %le %le %le %le \n",xyz_cadproj[0],xyz_cadproj[1],xyz_cadproj[2],(xyz_scnproj[0]-xyz_cadproj[0]),(xyz_scnproj[1]-xyz_cadproj[1]),(xyz_scnproj[2]-xyz_cadproj[2]));
  }
  fclose(insmp);
  
  
}

void CRIAQOPERATORS::ajouter_nois_normdir(char* magicfilename,int nummai,char* noisefile,char* outputfilename)
{
  MG_FILE gest(magicfilename);
  MG_MAILLAGE* mai;
  if (nummai==0) mai=gest.get_mg_maillage(nummai); else mai=gest.get_mg_maillageid(nummai);
/*
double xmin=1e308,ymin=1e308,zmin=1e308,xmax=-1e308,ymax=-1e308,zmax=-1e308;
TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
LISTE_MG_NOEUD::iterator it;
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
{
 if (no->get_x()<xmin) xmin=no->get_x();  
 if (no->get_y()<ymin) ymin=no->get_y();  
 if (no->get_z()<zmin) zmin=no->get_z();  
 if (no->get_x()>xmax) xmax=no->get_x();  
 if (no->get_y()>ymax) ymax=no->get_y();  
 if (no->get_z()>zmax) zmax=no->get_z();
 lstnoeud.ajouter(no); 
}
OT_VECTEUR_3D vecmin(xmin,ymin,zmin);
OT_VECTEUR_3D vecmax(xmax,ymax,zmax);
OT_VECTEUR_3D lengthvec(vecmin,vecmax);
double search_radius=5.0;//0.0000001;
//double search_radius=0.001*(lengthvec.get_longueur());
//char bboxdiag[1000];
//sprintf(bboxdiag,"the bounding box diagonal is: %lf , and the neighborhood search radius is: %lf ",vec.get_longueur(),search_radius);
//affiche((char*) bboxdiag);
OT_VECTEUR_3D average=(vecmin+vecmax)/2.;
double length=sqrt(lengthvec*lengthvec);
double bxr=2.0; //this can be changed sasan
xmin=average.get_x()-(length*bxr);ymin=average.get_y()-(length*bxr);zmin=average.get_z()-(length*bxr);
xmax=average.get_x()+(length*bxr);ymax=average.get_y()+(length*bxr);zmax=average.get_z()+(length*bxr);
TPL_OCTREE<MG_TRIANGLE*,MG_NOEUD*> octree;
octree.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
LISTE_MG_TRIANGLE::iterator it9;
for (MG_TRIANGLE* tri=mai->get_premier_triangle(it9);tri!=NULL;tri=mai->get_suivant_triangle(it9))
  octree.inserer(tri);
TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> octreends;
octreends.initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
for (MG_NOEUD* no=mai->get_premier_noeud(it);no!=NULL;no=mai->get_suivant_noeud(it))
  octreends.inserer(no);
*/

TPL_LISTE_ENTITE<double>  randomnumbr;
FILE *in=fopen(noisefile,"rt");
if (in==NULL)
  cout<< "file is not available"<<endl; 
while (!feof(in))
{
 char chaine[500];
 char* res=fgets(chaine,500,in);
 if(res!=NULL)
 {
  double rndnbr;
  int nb=sscanf(chaine,"%le",&rndnbr);
  if (nb!=-1 && nb!=1)
    cout<<"Wrong file format"<<endl;
  else if (nb!=-1 && nb==1)
  {
    randomnumbr.ajouter(rndnbr);
  }
 }
}
fclose(in);
cout<<"randomnumbr.get_nb()= "<<randomnumbr.get_nb()<<endl;
cout<<"mai->get_nb_mg_noeud()= "<<mai->get_nb_mg_noeud()<<endl;
LISTE_MG_NOEUD::iterator itnd;
MG_NOEUD* noeu=mai->get_premier_noeud(itnd);
for(int i=0;i<randomnumbr.get_nb();i++)
{
  //cout<<"i= "<<i<<endl;
  double rndno=randomnumbr.get(i);
  //for(MG_NOEUD* noeu=mai->get_premier_noeud(itnd);noeu!=NULL;noeu=mai->get_suivant_noeud(itnd))
  {
    TPL_MAP_ENTITE<MG_TRIANGLE*> lstneitri;
    /*double search_rad1=search_radius;
    int neitrinb=lstneitri.get_nb();
    while(neitrinb==0) 
    {
      octree.rechercher(noeu->get_x(),noeu->get_y(),noeu->get_z(),search_rad1,lstneitri);
      neitrinb=lstneitri.get_nb();
      search_rad1=search_rad1+0.00001*search_rad1;
    }
    cout<<"lstneitri.get_nb()= "<<lstneitri.get_nb()<<endl;
    */
     OT_VECTEUR_3D w_n_tot(0.,0.,0.);
     double norme_w_n_tot=0.;
     OT_VECTEUR_3D w2_n_tot(0.,0.,0.);
     double norme_w2_n_tot=0.;
     int nbtri=noeu->get_lien_triangle()->get_nb();
    for (int ii=0;ii<nbtri;ii++)
    {
        MG_TRIANGLE* tri=(MG_TRIANGLE*)noeu->get_lien_triangle()->get(ii);      
        //tri->change_origine(IMPOSE);
        double *xyzn1=tri->get_noeud1()->get_coord();
        double *xyzn2=tri->get_noeud2()->get_coord();
        double *xyzn3=tri->get_noeud3()->get_coord();
        OT_VECTEUR_3D vec1(xyzn1,xyzn3);
        OT_VECTEUR_3D vec2(xyzn1,xyzn2);
        OT_VECTEUR_3D n_tri=vec1&vec2;  // Normale au triangle

        double w=n_tri.get_longueur();  // 2 fois l'aire du triangle (pour éviter une division par 2)
        n_tri.norme();
        OT_VECTEUR_3D w_n=w*n_tri;
        double norme_w_n=w_n.get_longueur();
        w_n_tot=w_n_tot+w_n;
        norme_w_n_tot=norme_w_n_tot+norme_w_n;
    }
    
    OT_VECTEUR_3D normale_discrete=w_n_tot/norme_w_n_tot;       // Normale à la surface au noeud
    normale_discrete.norme();
    
    double xyz_new[3]; 
    xyz_new[0]=noeu->get_x()+(rndno*normale_discrete.get_x());
    xyz_new[1]=noeu->get_y()+(rndno*normale_discrete.get_y());
    xyz_new[2]=noeu->get_z()+(rndno*normale_discrete.get_z());
    noeu->change_coord(xyz_new);
  }
  noeu=mai->get_suivant_noeud(itnd);
}
gest.enregistrer(outputfilename);
}


Revision:	761
Committed:	Wed Nov 18 21:28:15 2015 UTC (9 years, 8 months ago) by sattarpa
File size:	141284 byte(s)
Log Message:	modification sur "bumparea_calcul" et "ins_point_withbc".