mailleur/src/remailleur.cpp

//---------------------------------------------------------------------------

#include "gestionversion.h"
#include <math.h>
#include "remailleur.h"
#include "mg_gestionnaire.h"
#include "mailleur3d.h"
#include "mailleur0d.h"
#include "mailleur1d.h"
#include "mailleur2d.h"
#include "ot_mathematique.h"
#include "mg_geometrie_outils.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)


REMAILLEUR::REMAILLEUR(class MG_GESTIONNAIRE* g1,MG_GESTIONNAIRE* g2,MG_GEOMETRIE* geo1,MG_GEOMETRIE* geo2,class FCT_TAILLE* fct_taille1,FCT_TAILLE* fct_taille2,MG_MAILLAGE* mori,MG_MAILLAGE* mmodi,VCT_COMPARAISON_RESULTAT& cmp):MAILLEUR(false),gestorigine(g1),gestmodifie(g2),carteori(fct_taille1),cartemod(fct_taille2),rescmp(cmp),maiorigine(mori),maimodifie(mmodi),geoorigine(geo1),geomodifie(geo2)
{
    int nbvolume=geomodifie->get_nb_mg_volume();
    lsttrifront=new TPL_LISTE_ENTITE<MG_TRIANGLE*>[nbvolume];
}

REMAILLEUR::~REMAILLEUR()
{
//TPL_MAP_ENTITE<class MG_SEGMENT_FRONTIERE*>::ITERATEUR it;
//for (MG_SEGMENT_FRONTIERE* seg=lstsegfront.get_premier(it);seg!=NULL;seg=lstsegfront.get_suivant(it) ) delete seg;
    delete octree;
    delete [] lsttrifront;
}


void REMAILLEUR::maille(MG_GROUPE_TOPOLOGIQUE* mggt)
{
    maille(20);
}

void REMAILLEUR::maille(int etape)
{
//etape 0 initialisation
    initialise_octree();
//etape 1 Destruction des entites autour des entites disparues
    if (etape<1) return;
    affiche((char*)"   Destruction autour des disparitions");
    TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*> lstdetruire;
    int nb=rescmp.get_nb_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_DISPARUE);
    for (int i=0;i<nb;i++)
    {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_DISPARUE,i);
        int dim=ele->get_dimension();
        if (dim==0)
        {
            MG_SOMMET* som=(MG_SOMMET*)ele;
            int nblien=som->get_lien_maillage()->get_nb();
            for (int j=0;j<nblien;j++)
            {
                MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
                if (maiorigine->get_mg_noeudid(no->get_id())!=NULL) lstdetruire.ajouter(no);
            }
        }
        if (dim==1)
        {
            MG_ARETE* are=(MG_ARETE*)ele;
            int nblien=are->get_lien_maillage()->get_nb();
            for (int j=0;j<nblien;j++)
            {
                MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
                if (maiorigine->get_mg_segmentid(seg->get_id())!=NULL)
                {
                    lstdetruire.ajouter(seg->get_noeud1());
                    lstdetruire.ajouter(seg->get_noeud2());
                }
            }
        }
        if (dim==2)
        {
            MG_FACE* face=(MG_FACE*)ele;
            int nblien=face->get_lien_maillage()->get_nb();
            for (int j=0;j<nblien;j++)
            {
                MG_TRIANGLE* tri=(MG_TRIANGLE*)face->get_lien_maillage()->get(j);
                if (maiorigine->get_mg_triangleid(tri->get_id())!=NULL)
                {
                    lstdetruire.ajouter(tri->get_noeud1());
                    lstdetruire.ajouter(tri->get_noeud2());
                    lstdetruire.ajouter(tri->get_noeud3());
                }
            }
        }
    }
    TPL_MAP_ENTITE<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
    for (MG_NOEUD* no=(MG_NOEUD*)lstdetruire.get_premier(it);no!=NULL;no=(MG_NOEUD*)lstdetruire.get_suivant(it))
    {
        double ereelle=0.;
        int nbseg=no->get_lien_segment()->get_nb();
        for (int i=0;i<nbseg;i++)
            ereelle=ereelle+no->get_lien_segment()->get(i)->get_longueur();
        ereelle=ereelle/nbseg;
        detruit_noeud(no,2*ereelle);
    }
//etape 2 : transfert des noeuds de sommets;
    if (etape<2) return;
    affiche((char*)"   Transfert sommet");
    nb=rescmp.get_nb_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE);
    for (int i=0;i<nb;i++)
    {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE,i);
        int dim=ele->get_dimension();
        if (dim==0)
        {
            MG_SOMMET* som=(MG_SOMMET*)ele;
            int nblien=som->get_lien_maillage()->get_nb();
            for (int j=0;j<nblien;j++)
            {
                MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
                if ((maiorigine->get_mg_noeudid(no->get_id())!=NULL) && (no->get_nouveau_numero()==CONSERVE))
                {
                    MG_SOMMET* nvsom=geomodifie->get_mg_sommetid(rescmp.get_liste_correspondance_modifie(som->get_id()));
                    transfert_noeud(no,nvsom);
                }
            }

        }
    }
// etape 3 maillage 0d
    if (etape<3) return;
    affiche((char*)"   MAILLAGE 0D");
    MAILLEUR0D m0d(maimodifie,geomodifie);
    m0d.adapte();
// etape 4 destruction autour des nouveaux sommets
    if (etape<4) return;
    affiche((char*)"   Destruction autour des nouveaux sommets");
    lstdetruire.vide();
    nb=rescmp.get_nb_liste_topologie(VCT_COMPARAISON_RESULTAT::MODIFIE_APPARUE);
    for (int i=0;i<nb;i++)
    {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(VCT_COMPARAISON_RESULTAT::MODIFIE_APPARUE,i);
        int dim=ele->get_dimension();
        if (dim==0)
        {
            MG_SOMMET* som=(MG_SOMMET*)ele;
            int nblien=som->get_lien_maillage()->get_nb();
            for (int j=0;j<nblien;j++)
            {
                MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
                if (maimodifie->get_mg_noeudid(no->get_id())!=NULL) lstdetruire.ajouter(no);
            }
        }
    }
    for (MG_ELEMENT_MAILLAGE* ele=lstdetruire.get_premier(it);ele!=NULL;ele=lstdetruire.get_suivant(it))
    {
        detruit_noeud(ele,-1,1);
    }
// etape 5 transferts des segments d arete
    if (etape<5) return;
    /*affiche("   Transfert arete");
    nb=rescmp.get_nb_liste_topologie(ORIGINE_CONSERVEE);
    for (int i=0;i<nb;i++)
        {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_CONSERVEE,i);
        int dim=ele->get_dimension();
        if (dim==1)
                {
                MG_ARETE* are=(MG_ARETE*)ele;
                int nblien=are->get_lien_maillage()->get_nb();
                for (int j=0;j<nblien;j++)
                        {
                        MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
                        if ((maiorigine->get_mg_segmentid(seg->get_id())!=NULL))
                                {
                                MG_NOEUD* n1=seg->get_noeud1();
                                MG_NOEUD* n2=seg->get_noeud2();
                                if (n1->get_nouveau_numero()==CONSERVE)
                                        if (n2->get_nouveau_numero()==CONSERVE)
                                                {
                                                MG_ARETE* nvare=geomodifie->get_mg_areteid(rescmp.get_liste_correspondance(are->get_id()));
                                                MG_NOEUD* nv1=get_liste_correspondance(n1);
                                                MG_NOEUD* nv2=get_liste_correspondance(n2);
                                                if (nv1==NULL) nv1=transfert_noeud(n1,nvare);
                                                if (nv2==NULL) nv2=transfert_noeud(n2,nvare);
                                                MG_SEGMENT *seg=new MG_SEGMENT(nvare,nv1,nv2,IMPOSE);
                                                maimodifie->ajouter_mg_segment(seg);
                                                }
                                }
                        }

                }
        }
    */
    affiche((char*)"   Creation des blocs de maille 1D");
    int nbarete=geomodifie->get_nb_mg_arete();
    nb=rescmp.get_nb_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE);
    for (int i=0;i<nb;i++)
    {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE,i);
        int dim=ele->get_dimension();
        if (dim==1)
        {
            std::multimap<unsigned long,MG_SEGMENT*,std::less <unsigned long> >  lstseg;
            MG_ARETE* arete=(MG_ARETE*)ele;
            TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
            for (MG_ELEMENT_MAILLAGE* ele=arete->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=arete->get_lien_maillage()->get_suivant(ittet))
            {
                MG_SEGMENT* seg=(MG_SEGMENT*)ele;
                if ((maiorigine->get_mg_segmentid(seg->get_id())!=NULL))
                {
                    MG_NOEUD* n1=seg->get_noeud1();
                    MG_NOEUD* n2=seg->get_noeud2();
                    if (n1->get_nouveau_numero()==CONSERVE)
                        if (n2->get_nouveau_numero()==CONSERVE)
                        {
                            seg->change_nouveau_numero(NONTRAITE);
                            std::pair<unsigned long,MG_SEGMENT*> tmp(n1->get_id(),seg);
                            lstseg.insert(tmp);
                            std::pair<unsigned long,MG_SEGMENT*> tmp2(n2->get_id(),seg);
                            lstseg.insert(tmp2);
                        }
                }
            }

            cree_bloc_maille(lstseg,arete);
        }
    }
// etape 6 Maillage des aretes
    if (etape<6) return;
    affiche((char*)"   MAILLAGE 1D");
    int nbare=geomodifie->get_nb_mg_arete();
    for (int i=0;i<nbare;i++)
    {
        MG_ARETE* arete=geomodifie->get_mg_arete(i);
        TPL_LISTE_ENTITE<double> param;
        int typeface=arete->get_courbe()->get_type_geometrique(param);
        std::multimap<unsigned long,BLOC_MAILLE_1D*,std::less <unsigned long> >::iterator it=lstb1d.begin();
        while (it!=lstb1d.end())
        {
            BLOC_MAILLE_1D *b1d=(*it).second;
            if ((b1d->etat==NONATTACHE) && (b1d->type==typeface))
            {
                MG_NOEUD* no1=b1d->lst[0]->get_noeud1();
                MG_NOEUD* no2=b1d->lst[0]->get_noeud2();
                double xyz1[3],xyz2[3],xyz[3];
                transfert_coord(no1->get_coord(),xyz1);
                transfert_coord(no2->get_coord(),xyz2);
                double t1,t2;
                arete->inverser(t1,xyz1);
                arete->inverser(t2,xyz2);
                int sensidentique=1;
                OT_VECTEUR_3D vec(xyz1,xyz2);
                vec.norme();
                double dxyz1[3];
                arete->deriver(t1,dxyz1);
                OT_VECTEUR_3D dvec(dxyz1);
                dvec.norme();
                double ps=dvec*vec;
                if (ps <0)  sensidentique=0;
                if (sensidentique==0)
                {
                    double t=t1;
                    t1=t2;
                    t2=t;
                    double xyz[3];
                    xyz[0]=xyz1[0];
                    xyz[1]=xyz1[1];
                    xyz[2]=xyz1[2];
                    xyz1[0]=xyz2[0];
                    xyz1[1]=xyz2[1];
                    xyz1[2]=xyz2[2];
                    xyz2[0]=xyz[0];
                    xyz2[1]=xyz[1];
                    xyz2[2]=xyz[2];
                }
                if (arete->get_courbe()->est_periodique()==1)
                {
                    //if (t1<arete->get_tmin()) t1=t1+arete->get_courbe()->get_periode();
                    if (t2<t1+1e-6) t2=t2+arete->get_courbe()->get_periode();
                }
                double xyz1tmp[3],xyz2tmp[3];
                arete->evaluer(t1,xyz1tmp);
                arete->evaluer(t2,xyz2tmp);
                OT_VECTEUR_3D vec1(xyz1,xyz1tmp);
                OT_VECTEUR_3D vec2(xyz2,xyz2tmp);
                double metrique[9];
                cartemod->evaluer(xyz1,metrique);
                double eps=(1./sqrt(metrique[0])*1e-6);
                double t=0.5*(t1+t2);
                if ((vec1.get_longueur()<eps)&&(vec2.get_longueur()<eps))
                {
                    if (arete->get_courbe()->est_periodique()==1)
                        if (t<arete->get_tmin()) t=t+arete->get_courbe()->get_periode();
                    if ((t>=arete->get_tmin()) && (t<=arete->get_tmax()))
                    {
                        int nbseg=b1d->lst.size();
                        b1d->etat=ATTACHE;
                        for (int i=0;i<nbseg;i++)
                        {
                            MG_SEGMENT *seg=b1d->lst[i];
                            MG_NOEUD* no1=seg->get_noeud1();
                            MG_NOEUD* no2=seg->get_noeud2();
                            double xyz1[3],xy2[3];
                            transfert_coord(no1->get_coord(),xyz1);
                            transfert_coord(no2->get_coord(),xyz2);
                            MG_NOEUD* nv1=get_liste_correspondance(no1);
                            MG_NOEUD* nv2=get_liste_correspondance(no2);
                            if (nv1==NULL) nv1=transfert_noeud(no1,arete);
                            if (nv2==NULL) nv2=transfert_noeud(no2,arete);
                            MG_SEGMENT *nvseg;
                            if (sensidentique==1)
                                nvseg=new MG_SEGMENT(arete,nv1,nv2,MAGIC::ORIGINE::IMPOSE);
                            else
                                nvseg=new MG_SEGMENT(arete,nv2,nv1,MAGIC::ORIGINE::IMPOSE);
                            maimodifie->ajouter_mg_segment(nvseg);
                        }
                    }
                }
            }
            it++;
        }
        MAILLEUR1D m1d(maimodifie,geomodifie,cartemod,arete);
        m1d.adapte();
    }
//cree_quadtree(maimodifie,quad,&lstsegfront);
// etape 7 Destruction autour des nouveaux segments
    if (etape<7) return;
    affiche((char*)"   Destruction autour des nouvelles aretes");
    lstdetruire.vide();
    /*nb=rescmp.get_nb_liste_topologie(MODIFIE_APPARUE);
    for (int i=0;i<nb;i++)
        {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(MODIFIE_APPARUE,i);
        int dim=ele->get_dimension();
        if (dim==1)
                {
                MG_ARETE* are=(MG_ARETE*)ele;
                int nblien=are->get_lien_maillage()->get_nb();
                for (int j=0;j<nblien;j++)
                        {
                        MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
                        if (maimodifie->get_mg_segmentid(seg->get_id())!=NULL)
                                {
                                lstnoeuddetruire.ajouter(seg->get_noeud1());
                                lstnoeuddetruire.ajouter(seg->get_noeud2());
                                }
                        }
                }
        }*/
    LISTE_MG_SEGMENT::iterator itnvseg;
    for (MG_SEGMENT* seg=maimodifie->get_premier_segment(itnvseg);seg!=NULL;seg=maimodifie->get_suivant_segment(itnvseg))
    {
        if (seg->get_origine()==MAGIC::ORIGINE::MAILLEUR_AUTO)
        {
            lstdetruire.ajouter(seg);
        }
    }
    for (MG_ELEMENT_MAILLAGE* ele=lstdetruire.get_premier(it);ele!=NULL;ele=lstdetruire.get_suivant(it))
    {
        detruit_noeud(ele,-1,2);
    }

// etape 8 Creation des blocs de mailles 2D
    if (etape<8) return;
    affiche((char*)"   Creation des blocs de maille 2D");
    nb=rescmp.get_nb_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE);
    for (int i=0;i<nb;i++)
    {
        MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(VCT_COMPARAISON_RESULTAT::ORIGINE_CONSERVEE,i);
        int dim=ele->get_dimension();
        if (dim==2)
        {
            std::multimap<unsigned long,MG_TRIANGLE*,std::less <unsigned long> >  lsttri;
            MG_FACE* face=(MG_FACE*)ele;
            TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
            int compteur=0,compteur2=0;
            for (MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=face->get_lien_maillage()->get_suivant(ittet))
            {
                MG_TRIANGLE* tri=(MG_TRIANGLE*)ele;
                if ((maiorigine->get_mg_triangleid(tri->get_id())!=NULL))
                {
                    MG_NOEUD* n1=tri->get_noeud1();
                    MG_NOEUD* n2=tri->get_noeud2();
                    MG_NOEUD* n3=tri->get_noeud3();
                    compteur2++;
                    if (n1->get_nouveau_numero()==CONSERVE)
                        if (n2->get_nouveau_numero()==CONSERVE)
                            if (n3->get_nouveau_numero()==CONSERVE)
                            {
                                tri->change_nouveau_numero(NONTRAITE);
                                std::pair<unsigned long,MG_TRIANGLE*> tmp(n1->get_id(),tri);
                                lsttri.insert(tmp);
                                std::pair<unsigned long,MG_TRIANGLE*> tmp2(n2->get_id(),tri);
                                lsttri.insert(tmp2);
                                std::pair<unsigned long,MG_TRIANGLE*> tmp3(n3->get_id(),tri);
                                lsttri.insert(tmp3);
                                compteur++;
                            }
                }
            }

            cree_bloc_maille(lsttri,face);
        }
    }
// etape 9 MAILLAGE 2D
    if (etape<9) return;
    affiche((char*)"   MAILLAGE 2D");
    nbfacemod=geomodifie->get_nb_mg_face();
    for (int i=0;i<nbfacemod;i++)
    {
        MG_GEOMETRIE_OUTILS ot;
        TPL_LISTE_ENTITE<MG_TRIANGLE*> lsttricontraint;
        MG_FACE* face=geomodifie->get_mg_face(i);
        TPL_LISTE_ENTITE<double> param;
        int typeface=face->get_surface()->get_type_geometrique(param);
        std::multimap<unsigned long,BLOC_MAILLE_2D*,std::less <unsigned long> >::iterator it=lstb2d.begin();
        while (it!=lstb2d.end())
        {
            BLOC_MAILLE_2D *b2d=(*it).second;
            if ((b2d->etat==NONATTACHE) && (b2d->type==typeface))
            {
                MG_NOEUD* no1=b2d->lst[0]->get_noeud1();
                MG_NOEUD* no2=b2d->lst[0]->get_noeud2();
                MG_NOEUD* no3=b2d->lst[0]->get_noeud3();
                double uv[2],uv1[2],uv2[2],uv3[2];
                double xyz1[3],xyz2[3],xyz3[3];
                transfert_coord(no1->get_coord(),xyz1);
                transfert_coord(no2->get_coord(),xyz2);
                transfert_coord(no3->get_coord(),xyz3);
                face->inverser(uv1,xyz1);
                face->inverser(uv2,xyz2);
                face->inverser(uv3,xyz3);
                if (face->get_surface()->est_periodique_u()==1)
                      {
                      double eps=1e-10*face->get_surface()->get_periode_u();
                      if (uv1[0]<-eps) uv1[0]=uv1[0]+face->get_surface()->get_periode_u();
                      if (uv1[0]>face->get_surface()->get_periode_u()-eps) uv1[0]=uv1[0]-face->get_surface()->get_periode_u();
                      if (uv2[0]<-eps) uv2[0]=uv2[0]+face->get_surface()->get_periode_u();
                      if (uv2[0]>face->get_surface()->get_periode_u()-eps) uv2[0]=uv2[0]-face->get_surface()->get_periode_u();
                      if (uv3[0]<-eps) uv3[0]=uv3[0]+face->get_surface()->get_periode_u();
                      if (uv3[0]>face->get_surface()->get_periode_u()-eps) uv3[0]=uv3[0]-face->get_surface()->get_periode_u();
                      }   
                if (face->get_surface()->est_periodique_v()==1)
                      {
                      double eps=1e-10*face->get_surface()->get_periode_v();
                      if (uv1[1]<-eps) uv1[1]=uv1[1]+face->get_surface()->get_periode_v();
                      if (uv1[1]>face->get_surface()->get_periode_v()-eps) uv1[1]=uv1[1]-face->get_surface()->get_periode_v();
                      if (uv2[1]<-eps) uv2[1]=uv2[1]+face->get_surface()->get_periode_v();
                      if (uv2[1]>face->get_surface()->get_periode_v()-eps) uv2[1]=uv2[1]-face->get_surface()->get_periode_v();
                      if (uv3[1]<-eps) uv3[1]=uv3[1]+face->get_surface()->get_periode_v();
                      if (uv3[1]>face->get_surface()->get_periode_v()-eps) uv3[1]=uv3[1]-face->get_surface()->get_periode_v();
                      }   
            
                double xyz1tmp[3];
                face->evaluer(uv1,xyz1tmp);
                double xyz2tmp[3];
                face->evaluer(uv2,xyz2tmp);
                double xyz3tmp[3];
                face->evaluer(uv3,xyz3tmp);
                OT_VECTEUR_3D vec1(xyz1,xyz1tmp);
                OT_VECTEUR_3D vec2(xyz2,xyz2tmp);
                OT_VECTEUR_3D vec3(xyz3,xyz3tmp);
                double metrique[9];
                cartemod->evaluer(xyz1,metrique);
                double eps=(1./sqrt(metrique[0])*1e-6);
                if (vec1.get_longueur()<eps)
                    if (vec2.get_longueur()<eps)
                        if (vec3.get_longueur()<eps)
                        {
                            OT_DECALAGE_PARAMETRE deca(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
                            double du=deca.calcul_decalage_parametre_u(uv1[0]);
                            double dv=deca.calcul_decalage_parametre_v(uv1[1]);
                            uv[0]=0.333333333333333333*(deca.decalage_parametre_u(uv1[0],du)+deca.decalage_parametre_u(uv2[0],du)+deca.decalage_parametre_u(uv3[0],du));
                            uv[1]=0.333333333333333333*(deca.decalage_parametre_v(uv1[1],dv)+deca.decalage_parametre_v(uv2[1],dv)+deca.decalage_parametre_v(uv3[1],dv));
                            uv[0]=deca.decalage_parametre_u(uv[0],-du);
                            uv[1]=deca.decalage_parametre_v(uv[1],-dv);
                            double dist=ot.calcule_distance_contour_face_uv(uv,face);
                            if (dist>0.)
                            {
                                int sensidentique=1;
                                OT_VECTEUR_3D vec12(xyz1,xyz2);
                                OT_VECTEUR_3D vec13(xyz1,xyz3);
                                double normal[3];
                                face->calcul_normale_unitaire(uv1,normal);
                                OT_VECTEUR_3D nor=vec13 & vec12;
                                nor.norme();
                                if ((nor*normal)<0.) sensidentique=0;
                                int nbtri=b2d->lst.size();
                                b2d->etat=ATTACHE;
                                //cout << "               attachment du bloc " << b2d->face->get_id();
                                for (int i=0;i<nbtri;i++)
                                {
                                    MG_TRIANGLE *tri=b2d->lst[i];
                                    MG_NOEUD* no1=tri->get_noeud1();
                                    MG_NOEUD* no2=tri->get_noeud2();
                                    MG_NOEUD* no3=tri->get_noeud3();
                                    double xyz1[3],xy2[3],xyz3[3];
                                    transfert_coord(no1->get_coord(),xyz1);
                                    transfert_coord(no2->get_coord(),xyz2);
                                    transfert_coord(no3->get_coord(),xyz3);
                                    MG_NOEUD* nvno1=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],MAGIC::ORIGINE::IMPOSE);
                                    MG_NOEUD* nvno2=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],MAGIC::ORIGINE::IMPOSE);
                                    MG_NOEUD* nvno3=new MG_NOEUD(NULL,xyz3[0],xyz3[1],xyz3[2],MAGIC::ORIGINE::IMPOSE);
                                    MG_TRIANGLE* nvtri;
                                    if (sensidentique==1)
                                        nvtri=new MG_TRIANGLE(NULL,nvno1,nvno2,nvno3,NULL,NULL,NULL,MAGIC::ORIGINE::IMPOSE);
                                    else
                                        nvtri=new MG_TRIANGLE(NULL,nvno2,nvno1,nvno3,NULL,NULL,NULL,MAGIC::ORIGINE::IMPOSE);
                                    lsttricontraint.ajouter(nvtri);
                                    nvno1->change_nouveau_numero(no1->get_id());
                                    nvno2->change_nouveau_numero(no2->get_id());
                                    nvno3->change_nouveau_numero(no3->get_id());
                                }
                            }
                        }
            }
            it++;
        }
        MAILLEUR2D m2d(maimodifie,geomodifie,cartemod,face);
        char mess[50];
        sprintf(mess,"      Face %i identificateur %lu",i,face->get_id());
        affiche(mess);
        m2d.maille(face,NULL,&(lsttricontraint));
        int nbtricon=lsttricontraint.get_nb();
        for (int i=0;i<nbtricon;i++)
        {
            MG_TRIANGLE* tri=lsttricontraint.get(0);
            MG_NOEUD* noeud1=tri->get_noeud1();
            MG_NOEUD* noeud2=tri->get_noeud2();
            MG_NOEUD* noeud3=tri->get_noeud3();
            lsttricontraint.supprimer(tri);
            delete tri;
            delete noeud1;
            delete noeud2;
            delete noeud3;
        }
    }
    std::multimap<unsigned long,BLOC_MAILLE_2D*,std::less <unsigned long> >::iterator itfin2=lstb2d.begin();
    while (itfin2!=lstb2d.end())
    {
        delete (*itfin2).second;
        itfin2++;
    }
    cree_liste_frontiere(maimodifie);
// etape 10 Destruction autour des nouveaux triangles
    if (etape<10) return;
    affiche((char*)"   Destruction autour des nouvelles faces");
    lstdetruire.vide();
    LISTE_MG_TRIANGLE::iterator itnvtri;
    for (MG_TRIANGLE* tri=maimodifie->get_premier_triangle(itnvtri);tri!=NULL;tri=maimodifie->get_suivant_triangle(itnvtri))
    {
        if (tri->get_origine()==MAGIC::ORIGINE::MAILLEUR_AUTO)
        {
            lstdetruire.ajouter(tri);
        }
    }
    for (MG_ELEMENT_MAILLAGE* ele=lstdetruire.get_premier(it);ele!=NULL;ele=lstdetruire.get_suivant(it))
    {
        detruit_noeud(ele,-1,3);
    }
// etape 11 Creation des blocs de mailles 3D
    if (etape<11) return;
    affiche((char*)"   Creation des blocs de maille 3D");
    int nbvolume=geoorigine->get_nb_mg_volume();
    for (int i=0;i<nbvolume;i++)
    {
        std::multimap<unsigned long,MG_TETRA*,std::less <unsigned long> >  lsttet;
        MG_VOLUME* vol=geoorigine->get_mg_volume(i);
        int nblien=vol->get_lien_maillage()->get_nb();
        TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
        for (MG_ELEMENT_MAILLAGE* ele=vol->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=vol->get_lien_maillage()->get_suivant(ittet))
        {
            MG_TETRA* tet=(MG_TETRA*)ele;
            if ((maiorigine->get_mg_tetraid(tet->get_id())!=NULL))
            {
                MG_NOEUD* n1=tet->get_noeud1();
                MG_NOEUD* n2=tet->get_noeud2();
                MG_NOEUD* n3=tet->get_noeud3();
                MG_NOEUD* n4=tet->get_noeud4();
                if (n1->get_nouveau_numero()==CONSERVE)
                    if (n2->get_nouveau_numero()==CONSERVE)
                        if (n3->get_nouveau_numero()==CONSERVE)
                            if (n4->get_nouveau_numero()==CONSERVE)
                            {
                                tet->change_nouveau_numero(NONTRAITE);
                                std::pair<unsigned long,MG_TETRA*> tmp(n1->get_id(),tet);
                                lsttet.insert(tmp);
                                std::pair<unsigned long,MG_TETRA*> tmp2(n2->get_id(),tet);
                                lsttet.insert(tmp2);
                                std::pair<unsigned long,MG_TETRA*> tmp3(n3->get_id(),tet);
                                lsttet.insert(tmp3);
                                std::pair<unsigned long,MG_TETRA*> tmp4(n4->get_id(),tet);
                                lsttet.insert(tmp4);
                            }
            }

        }
        cree_bloc_maille(lsttet,vol);
    }
// etape 12 Maillage 3D
    if (etape<12) return;
    affiche((char*)"   MAILLEUR 3D");
    nbvolume=geomodifie->get_nb_mg_volume();
    for (int i=0;i<nbvolume;i++)
    {
        TPL_LISTE_ENTITE<MG_TETRA*> lsttetcontraint;
        MG_VOLUME* vol=geomodifie->get_mg_volume(i);
        std::multimap<unsigned long,BLOC_MAILLE_3D*,std::less <unsigned long> >::iterator it=lstb3d.begin();
        while (it!=lstb3d.end())
        {
            BLOC_MAILLE_3D *b3d=(*it).second;
            if (b3d->etat==NONATTACHE)
            {
                MG_NOEUD* no1=b3d->lst[0]->get_noeud1();
                MG_NOEUD* no2=b3d->lst[0]->get_noeud2();
                MG_NOEUD* no3=b3d->lst[0]->get_noeud3();
                MG_NOEUD* no4=b3d->lst[0]->get_noeud4();
                double xyz1[3],xyz2[3],xyz3[3],xyz4[3],xyz[3];
                transfert_coord(no1->get_coord(),xyz1);
                transfert_coord(no2->get_coord(),xyz2);
                transfert_coord(no3->get_coord(),xyz3);
                transfert_coord(no4->get_coord(),xyz4);
                xyz[0]=0.25*(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0]);
                xyz[1]=0.25*(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1]);
                xyz[2]=0.25*(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2]);
                int res=point_appartient_volume(xyz,i);
                if (res==1)
                {
                    int nbtet=b3d->lst.size();
                    b3d->etat=ATTACHE;
                    //cout << "bloc recupere " << b3d->get_id() << endl;
                    for (int j=0;j<nbtet;j++)
                    {
                        MG_TETRA *tet=b3d->lst[j];
                        MG_NOEUD* no1=tet->get_noeud1();
                        MG_NOEUD* no2=tet->get_noeud2();
                        MG_NOEUD* no3=tet->get_noeud3();
                        MG_NOEUD* no4=tet->get_noeud4();
                        double xyz1[3],xy2[3],xyz3[3],xyz4[3];
                        transfert_coord(no1->get_coord(),xyz1);
                        transfert_coord(no2->get_coord(),xyz2);
                        transfert_coord(no3->get_coord(),xyz3);
                        transfert_coord(no4->get_coord(),xyz4);
                        MG_NOEUD* nvno1=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],MAGIC::ORIGINE::IMPOSE);
                        MG_NOEUD* nvno2=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],MAGIC::ORIGINE::IMPOSE);
                        MG_NOEUD* nvno3=new MG_NOEUD(NULL,xyz3[0],xyz3[1],xyz3[2],MAGIC::ORIGINE::IMPOSE);
                        MG_NOEUD* nvno4=new MG_NOEUD(NULL,xyz4[0],xyz4[1],xyz4[2],MAGIC::ORIGINE::IMPOSE);
                        MG_TETRA* nvtet=new MG_TETRA(NULL,nvno1,nvno2,nvno3,nvno4,NULL,NULL,NULL,NULL,MAGIC::ORIGINE::IMPOSE);
                        lsttetcontraint.ajouter(nvtet);
                        nvno1->change_nouveau_numero(no1->get_id());
                        nvno2->change_nouveau_numero(no2->get_id());
                        nvno3->change_nouveau_numero(no3->get_id());
                        nvno4->change_nouveau_numero(no4->get_id());
                        nvtet->change_nouveau_numero(tet->get_id());
                    }
                }
            }
            it++;
        }
        MAILLEUR3D m3d(maimodifie,geomodifie,cartemod,vol);
        char mess[50];
        sprintf(mess,"      Volume %i identificateur %lu",i,vol->get_id());
        affiche(mess);
        m3d.active_affichage(affiche2);
        m3d.maille(vol,NULL,&(lsttetcontraint));
        int nbtetcon=lsttetcontraint.get_nb();
        for (int i=0;i<nbtetcon;i++)
        {
            MG_TETRA* tet=lsttetcontraint.get(0);
            MG_NOEUD* noeud1=tet->get_noeud1();
            MG_NOEUD* noeud2=tet->get_noeud2();
            MG_NOEUD* noeud3=tet->get_noeud3();
            MG_NOEUD* noeud4=tet->get_noeud4();
            lsttetcontraint.supprimer(tet);
            delete tet;
            delete noeud1;
            delete noeud2;
            delete noeud3;
            delete noeud4;
        }
    }
    std::multimap<unsigned long,BLOC_MAILLE_3D*,std::less <unsigned long> >::iterator itfin3=lstb3d.begin();
    while (itfin3!=lstb3d.end())
    {
        delete (*itfin3).second;
        itfin3++;
    }
// etape 13 Recuperation des correspondances
    if (etape<13) return;
    affiche((char*)"   Sauvegarde des correspondances");
    LISTE_MG_NOEUD::iterator itnoeud;
    for (MG_NOEUD* noeud=maimodifie->get_premier_noeud(itnoeud);noeud!=NULL;noeud=maimodifie->get_suivant_noeud(itnoeud))
    {
        if (noeud->get_origine()==MAGIC::ORIGINE::IMPOSE)
            rescmp.ajouter_correspondance(noeud->get_nouveau_numero(),noeud->get_id());
    }
    LISTE_MG_TETRA::iterator ittet;
    for (MG_TETRA* tet=maimodifie->get_premier_tetra(ittet);tet!=NULL;tet=maimodifie->get_suivant_tetra(ittet))
    {
        if (tet->get_origine()==MAGIC::ORIGINE::IMPOSE)
            rescmp.ajouter_correspondance(tet->get_nouveau_numero(),tet->get_id());
    }

}


void REMAILLEUR::initialise_octree(void)
{
    BOITE_3D boite=carteori->get_boite_3D();
    double xmin=boite.get_xmin();
    double ymin=boite.get_ymin();
    double zmin=boite.get_zmin();
    double xmax=boite.get_xmax();
    double ymax=boite.get_ymax();
    double zmax=boite.get_zmax();
    octree=new TPL_OCTREE_FCT<MG_NOEUD*,FCT_TAILLE >;
    octree->initialiser(*carteori,xmin,ymin,zmin,xmax,ymax,zmax);
    LISTE_MG_NOEUD::iterator it;
    for (MG_NOEUD* noeud=maiorigine->get_premier_noeud(it);noeud!=NULL;noeud=maiorigine->get_suivant_noeud(it))
    {
        noeud->change_nouveau_numero(CONSERVE);
        octree->inserer(noeud);
    }

}

void REMAILLEUR::detruit_noeud(MG_ELEMENT_MAILLAGE* elebase,double distance,int type)
{
    if (elebase->get_origine()!=MAGIC::ORIGINE::MAILLEUR_AUTO) return;
    double xyz[3];
    double distmin;
    OT_VECTEUR_3D direction;
    if (type==0)
    {
        MG_NOEUD* no=(MG_NOEUD*)elebase;
        OT_VECTEUR_3D vecnoeud(no->get_coord());
        //OT_VECTEUR_3D vecnoeud2=rescmp.change_coord_mod_ori(vecnoeud);
        OT_VECTEUR_3D vecnoeud2=vecnoeud;
        xyz[0]=vecnoeud2.get_x();
        xyz[1]=vecnoeud2.get_y();
        xyz[2]=vecnoeud2.get_z();
    }
    if (type==1)
    {
        MG_NOEUD* no=(MG_NOEUD*)elebase;
        OT_VECTEUR_3D vecnoeud(no->get_coord());
        OT_VECTEUR_3D vecnoeud2=rescmp.change_coord_mod_ori(vecnoeud);
        xyz[0]=vecnoeud2.get_x();
        xyz[1]=vecnoeud2.get_y();
        xyz[2]=vecnoeud2.get_z();
    }
    if (type==2)
    {
        MG_SEGMENT* seg=(MG_SEGMENT*)elebase;
        OT_VECTEUR_3D vecnoeud1(seg->get_noeud1()->get_coord());
        OT_VECTEUR_3D vecnoeud2(seg->get_noeud2()->get_coord());
        OT_VECTEUR_3D vecnoeud=0.5*(vecnoeud1+vecnoeud2);
        OT_VECTEUR_3D vecnoeud3=rescmp.change_coord_mod_ori(vecnoeud);
        xyz[0]=vecnoeud3.get_x();
        xyz[1]=vecnoeud3.get_y();
        xyz[2]=vecnoeud3.get_z();
        direction=vecnoeud2-vecnoeud1;
        distmin=direction.get_longueur()/2.;
        direction.norme();
    }
    if (type==3)
    {
        MG_TRIANGLE* tri=(MG_TRIANGLE*)elebase;
        OT_VECTEUR_3D vecnoeud1(tri->get_noeud1()->get_coord());
        OT_VECTEUR_3D vecnoeud2(tri->get_noeud2()->get_coord());
        OT_VECTEUR_3D vecnoeud3(tri->get_noeud3()->get_coord());
        OT_VECTEUR_3D vecnoeud=0.3333333333333333*(vecnoeud1+vecnoeud2+vecnoeud3);
        OT_VECTEUR_3D vecnoeud4=rescmp.change_coord_mod_ori(vecnoeud);
        xyz[0]=vecnoeud4.get_x();
        xyz[1]=vecnoeud4.get_y();
        xyz[2]=vecnoeud4.get_z();
        distmin=std::max((vecnoeud-vecnoeud1).get_longueur(),(vecnoeud-vecnoeud2).get_longueur());
        distmin=std::max((vecnoeud-vecnoeud3).get_longueur(),distmin);
        direction=(vecnoeud3-vecnoeud1)&(vecnoeud2-vecnoeud1);
        direction.norme();
    }
    TPL_MAP_ENTITE<MG_NOEUD*> lst;
    double Ea,Eb;
    if (type!=0)
    {
        double m[9];
        carteori->evaluer(xyz,m);
        Ea=1./sqrt(m[0]);
    }
    octree->rechercher(xyz[0],xyz[1],xyz[2],2*Ea,lst);
    double alpha=1.05;
    distmin=distmin*alpha;
    LISTE_MG_NOEUD::iterator it;
    for (MG_NOEUD* noeud=lst.get_premier(it);noeud!=NULL;noeud=lst.get_suivant(it))
    {
        double *xyztmp=noeud->get_coord();
        OT_VECTEUR_3D vec(xyz,xyztmp);
        if (type!=0)
        {
            double m[9];
            carteori->evaluer(xyztmp,m);
            Eb=1./sqrt(m[0]);
            distance=std::max(1.5*Ea,1.5*Eb);
            distance=std::max(distance,2.6*fabs(Ea-Eb));
        }
        if (type==2)
        {
            double cosinus=vec*direction/vec.get_longueur();
            double sinus=sqrt(1-cosinus*cosinus);
            double facteur;
            if (OPERATEUR::egal(cosinus,0.,0.000001)) facteur=distance;
            else
            {
                double delta=sinus*sinus+4.*distance*distance/distmin/distmin*cosinus*cosinus;
                facteur=(-sinus+sqrt(delta))/2./distance/cosinus/cosinus*distmin*distmin;
                distance=facteur;
            }
        }
        if (type==3)
        {
            double sinus=vec*direction/vec.get_longueur();
            if (sinus>-0.0000001)
            {
                double cosinus=sqrt(1-sinus*sinus);
                double facteur;
                if (OPERATEUR::egal(cosinus,0.,0.000001)) facteur=distance;
                else
                {
                    double delta=sinus*sinus+4*distance*distance/distmin/distmin*cosinus*cosinus;
                    facteur=(-sinus+sqrt(delta))/2./distance/cosinus/cosinus*distmin*distmin;
                    distance=facteur;
                }
            }

        }
        if (vec.get_longueur()<distance)
        {
            noeud->change_nouveau_numero(DETRUIT);
        }

    }

}


/*void REMAILLEUR::cree_quadtree(MG_MAILLAGE* mai,TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*> **quad,TPL_MAP_ENTITE<MG_SEGMENT_FRONTIERE*> *lstsegfront)
{
int nbface=mai->get_mg_geometrie()->get_nb_mg_face();
for (int i=0;i<nbface;i++)
        {
        ((quad[i]))=new TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*>;
        TPL_LISTE_ENTITE<MG_NOEUD*> lstn;
        double umax=-1e300;
        double vmax=-1e300;
        double umin=1e300;
        double vmin=1e300;
        MG_FACE* face=mai->get_mg_geometrie()->get_mg_face(i);
        int nbboucle=face->get_nb_mg_boucle();
        int nbsegfravant=lstsegfront->get_nb();
        for (int j=0;j<nbboucle;j++)
                {
                MG_BOUCLE* bou=face->get_mg_boucle(j);
                int nbarete=bou->get_nb_mg_coarete();
                for (int k=0;k<nbarete;k++)
                        {
                        MG_ARETE* are=bou->get_mg_coarete(k)->get_arete();
                        int nbseg=are->get_lien_maillage()->get_nb();
                                for (int l=0;l<nbseg;l++)
                                        {
                                        MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(l);
                                        if (mai->get_mg_segmentid(seg->get_id())==NULL) continue;
                                        double *xyz=seg->get_noeud1()->get_coord();
                                        double uv[2];
                                        face->inverser(uv,xyz);
                                        seg->get_noeud1()->change_u(uv[0]);
                                        seg->get_noeud1()->change_v(uv[1]);
                                        xyz=seg->get_noeud2()->get_coord();
                                        face->inverser(uv,xyz);
                                        seg->get_noeud2()->change_u(uv[0]);
                                        seg->get_noeud2()->change_v(uv[1]);
                                        lstn.ajouter(seg->get_noeud2());
                                        if (uv[0]<umin) umin=uv[0];
                                        if (uv[0]>umax) umax=uv[0];
                                        if (uv[1]<vmin) vmin=uv[1];
                                        if (uv[1]>vmax) vmax=uv[1];
                                        MG_SEGMENT_FRONTIERE* segfr=new MG_SEGMENT_FRONTIERE(seg);
                                        lstsegfront->ajouter(segfr);
                                        }
                        }
                }
        double periodeenu=face->get_surface()->get_periode_u();
        double periodeenv=face->get_surface()->get_periode_v();
        if (periodeenu!=0.0)
                {
                 umin=0.;
                 umax=periodeenu;
                }
        else
                {
                 double diff=umax-umin;
                 umin=umin-0.125*diff;
                 umax=umax+0.125*diff;
                 }
        if (periodeenv!=0.0)
                 {
                 vmin=0.;
                 vmax=periodeenv;
                 }
        else
                {
                double diff=vmax-vmin;
                vmin=vmin-0.125*diff;
                vmax=vmax+0.125*diff;
                }
        (quad[i])->initialiser(&lstn,1,umin,vmin,umax,vmax,face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
        int nbsegfr=lstsegfront->get_nb();
        for (int j=nbsegfravant;j<nbsegfr;j++)
                {
                double *uv=lstsegfront->get(j)->get_uv1();
                uv[0]=lstsegfront->get(j)->get_segment()->get_noeud1()->get_u();
                uv[1]=lstsegfront->get(j)->get_segment()->get_noeud1()->get_v();
                uv=lstsegfront->get(j)->get_uv2();
                uv[0]=lstsegfront->get(j)->get_segment()->get_noeud2()->get_u();
                uv[1]=lstsegfront->get(j)->get_segment()->get_noeud2()->get_v();
                (quad[i])->inserer( lstsegfront->get(j));
                }
        }

}
*/
void REMAILLEUR::cree_liste_frontiere(MG_MAILLAGE* mai)
{
    int nbvolume=mai->get_mg_geometrie()->get_nb_mg_volume();
    for (int i=0;i<nbvolume;i++)
    {
        MG_VOLUME* vol=mai->get_mg_geometrie()->get_mg_volume(i);
        int nbco=vol->get_nb_mg_coquille();
        for (int j=0;j<nbco;j++)
        {
            MG_COQUILLE* coq=vol->get_mg_coquille(j);
            int nbface=coq->get_nb_mg_coface();
            for (int k=0;k<nbface;k++)
            {
                MG_FACE* face=coq->get_mg_coface(k)->get_face();
                int nbseg=face->get_lien_maillage()->get_nb();
                for (int l=0;l<nbseg;l++)
                {
                    MG_TRIANGLE* tri=(MG_TRIANGLE*)face->get_lien_maillage()->get(l);
                    if (mai->get_mg_triangleid(tri->get_id())==NULL) continue;
                    lsttrifront[i].ajouter(tri);
                }
            }
        }


    }
}


void REMAILLEUR::ajouter_liste_correspondance_noeud(MG_NOEUD* n1,MG_NOEUD* n2)
{
    CORRESPONDANCENOEUD corr(n1,n2);
    std::pair<unsigned long,CORRESPONDANCENOEUD> tmp(n1->get_id(),corr);
    lstcorrnoeud.insert(tmp);
}

MG_NOEUD* REMAILLEUR::get_liste_correspondance(MG_NOEUD* n1)
{
    MG_NOEUD *n2=NULL;
    std::map<unsigned long,CORRESPONDANCENOEUD,std::less<unsigned long> >::iterator j=lstcorrnoeud.find(n1->get_id());
    if (j!=lstcorrnoeud.end())
        n2=(*j).second.noeudmod;
    return n2;
}


MG_NOEUD* REMAILLEUR::transfert_noeud(MG_NOEUD* no,MG_ELEMENT_TOPOLOGIQUE* ele)
{
    double *xyz=no->get_coord();
    OT_VECTEUR_3D pt(xyz);
    OT_VECTEUR_3D nvpt=rescmp.change_coord_ori_mod(pt);
    MG_NOEUD* noeud=new MG_NOEUD(ele,nvpt.get_x(),nvpt.get_y(),nvpt.get_z(),MAGIC::ORIGINE::IMPOSE);
    maimodifie->ajouter_mg_noeud(noeud);
    ajouter_liste_correspondance_noeud(no,noeud);
    noeud->change_nouveau_numero(no->get_id());
    return noeud;
}

void REMAILLEUR::transfert_coord(double *xyz,double *xyz2)
{
    OT_VECTEUR_3D pt(xyz);
    OT_VECTEUR_3D nvpt=rescmp.change_coord_ori_mod(pt);
    xyz2[0]=nvpt.get_x();
    xyz2[1]=nvpt.get_y();
    xyz2[2]=nvpt.get_z();
    return;
}

void REMAILLEUR::cree_bloc_maille(std::multimap<unsigned long,MG_SEGMENT*,std::less <unsigned long> >  &lst,MG_ARETE* arete)
{
    while (lst.size()>0)
    {
        std::multimap<unsigned long,MG_SEGMENT*,std::less <unsigned long> >::iterator it=lst.begin();
        MG_SEGMENT* seg=(*it).second;
        lst.erase(it);
        TPL_LISTE_ENTITE<double> param;
        BLOC_MAILLE_1D *b1d=new BLOC_MAILLE_1D(arete,arete->get_courbe()->get_type_geometrique(param));
        std::pair<unsigned long,BLOC_MAILLE_1D*> tmp(b1d->get_id(),b1d);
        lstb1d.insert(tmp);
        b1d->lst.push_back(seg);
        seg->change_nouveau_numero(TRAITE);
        char mess[255];
        sprintf(mess,"     Arete %lu , Bloc %d",arete->get_id(),(int)lstb1d.size());
        //affiche(mess);
        int num=0;
        while (num!=b1d->lst.size())
        {
            MG_SEGMENT* itseg=b1d->lst[num++];
            MG_SEGMENT* seg=trouve_segment(lst,itseg->get_noeud1());
            while (seg!=NULL)
            {
                if (seg->get_nouveau_numero()==NONTRAITE)
                {
                    b1d->lst.push_back(seg);
                    seg->change_nouveau_numero(TRAITE);
                }
                seg=trouve_segment(lst,itseg->get_noeud1());
            }
            seg=trouve_segment(lst,itseg->get_noeud2());
            while (seg!=NULL)
            {
                if (seg->get_nouveau_numero()==NONTRAITE)
                {
                    b1d->lst.push_back(seg);
                    seg->change_nouveau_numero(TRAITE);
                }
                seg=trouve_segment(lst,itseg->get_noeud2());
            }
        }

        sprintf(mess,"         Segment %d ",(int)b1d->lst.size());
        //affiche(mess);
    }
}

void REMAILLEUR::cree_bloc_maille(std::multimap<unsigned long,MG_TRIANGLE*,std::less <unsigned long> >  &lst,MG_FACE* face)
{
    while (lst.size()>0)
    {
        std::multimap<unsigned long,MG_TRIANGLE*,std::less <unsigned long> >::iterator it=lst.begin();
        MG_TRIANGLE* tri=(*it).second;
        lst.erase(it);
        TPL_LISTE_ENTITE<double> param;
        BLOC_MAILLE_2D *b2d=new BLOC_MAILLE_2D(geoorigine->get_mg_faceid(face->get_id()),face->get_surface()->get_type_geometrique(param));
        std::pair<unsigned long,BLOC_MAILLE_2D*> tmp(b2d->get_id(),b2d);
        lstb2d.insert(tmp);
        b2d->lst.push_back(tri);
        tri->change_nouveau_numero(TRAITE);
        char mess[255];
        sprintf(mess,"     Face %lu , Bloc %d",face->get_id(),(int)lstb2d.size());
        //affiche(mess);
        int num=0;
        while (num!=b2d->lst.size())
        {
            MG_TRIANGLE* ittri=b2d->lst[num++];
            MG_TRIANGLE* tri=trouve_triangle(lst,ittri->get_noeud1());
            while (tri!=NULL)
            {
                if (tri->get_nouveau_numero()==NONTRAITE)
                {
                    b2d->lst.push_back(tri);
                    tri->change_nouveau_numero(TRAITE);
                }
                tri=trouve_triangle(lst,ittri->get_noeud1());
            }
            tri=trouve_triangle(lst,ittri->get_noeud2());
            while (tri!=NULL)
            {
                if (tri->get_nouveau_numero()==NONTRAITE)
                {
                    b2d->lst.push_back(tri);
                    tri->change_nouveau_numero(TRAITE);
                }
                tri=trouve_triangle(lst,ittri->get_noeud2());
            }
            tri=trouve_triangle(lst,ittri->get_noeud3());
            while (tri!=NULL)
            {
                if (tri->get_nouveau_numero()==NONTRAITE)
                {
                    b2d->lst.push_back(tri);
                    tri->change_nouveau_numero(TRAITE);
                }
                tri=trouve_triangle(lst,ittri->get_noeud3());
            }
        }

        sprintf(mess,"         Triangle %d ",(int)b2d->lst.size());
        //affiche(mess);
    }
}

void REMAILLEUR::cree_bloc_maille(std::multimap<unsigned long,MG_TETRA*,std::less <unsigned long> >  &lst,MG_VOLUME* vol)
{
    while (lst.size()>0)
    {
        std::multimap<unsigned long,MG_TETRA*,std::less <unsigned long> >::iterator it=lst.begin();
        MG_TETRA* tet=(*it).second;
        lst.erase(it);
        BLOC_MAILLE_3D *b3d=new BLOC_MAILLE_3D(geoorigine->get_mg_volumeid(vol->get_id()));
        std::pair<unsigned long,BLOC_MAILLE_3D*> tmp(b3d->get_id(),b3d);
        lstb3d.insert(tmp);
        b3d->lst.push_back(tet);
        tet->change_nouveau_numero(TRAITE);
        char mess[255];
        sprintf(mess,"     Volume %lu , Bloc %d",vol->get_id(),(int)lstb3d.size());
        //affiche(mess);
        int num=0;
        while (num!=b3d->lst.size())
        {
            MG_TETRA* ittet=b3d->lst[num++];
            MG_TETRA* tet=trouve_tetra(lst,ittet->get_noeud1());
            while (tet!=NULL)
            {
                if (tet->get_nouveau_numero()==NONTRAITE)
                {
                    b3d->lst.push_back(tet);
                    tet->change_nouveau_numero(TRAITE);
                }
                tet=trouve_tetra(lst,ittet->get_noeud1());
            }
            tet=trouve_tetra(lst,ittet->get_noeud2());
            while (tet!=NULL)
            {
                if (tet->get_nouveau_numero()==NONTRAITE)
                {
                    b3d->lst.push_back(tet);
                    tet->change_nouveau_numero(TRAITE);
                }
                tet=trouve_tetra(lst,ittet->get_noeud2());
            }
            tet=trouve_tetra(lst,ittet->get_noeud3());
            while (tet!=NULL)
            {
                if (tet->get_nouveau_numero()==NONTRAITE)
                {
                    b3d->lst.push_back(tet);
                    tet->change_nouveau_numero(TRAITE);
                }
                tet=trouve_tetra(lst,ittet->get_noeud3());
            }
            tet=trouve_tetra(lst,ittet->get_noeud4());
            while (tet!=NULL)
            {
                if (tet->get_nouveau_numero()==NONTRAITE)
                {
                    b3d->lst.push_back(tet);
                    tet->change_nouveau_numero(TRAITE);
                }
                tet=trouve_tetra(lst,ittet->get_noeud4());
            }
        }

        sprintf(mess,"         Tetra %d ",(int)b3d->lst.size());
        //affiche(mess);
    }
}


MG_SEGMENT* REMAILLEUR::trouve_segment(std::multimap<unsigned long,MG_SEGMENT*,std::less <unsigned long> >  &lst,MG_NOEUD* no)
{
    MG_SEGMENT* seg=NULL;
    unsigned long id=no->get_id();
    std::multimap<unsigned long,MG_SEGMENT*,std::less <unsigned long> >::iterator it=lst.find(id);
    if (it!=lst.end())
    {
        seg=(*it).second;
        lst.erase(it);
    }
    return seg;
}

MG_TRIANGLE* REMAILLEUR::trouve_triangle(std::multimap<unsigned long,MG_TRIANGLE*,std::less <unsigned long> >  &lst,MG_NOEUD* no)
{
    MG_TRIANGLE* tri=NULL;
    unsigned long id=no->get_id();
    std::multimap<unsigned long,MG_TRIANGLE*,std::less <unsigned long> >::iterator it=lst.find(id);
    if (it!=lst.end())
    {
        tri=(*it).second;
        lst.erase(it);
    }
    return tri;
}

MG_TETRA* REMAILLEUR::trouve_tetra(std::multimap<unsigned long,MG_TETRA*,std::less <unsigned long> >  &lst,MG_NOEUD* no)
{
    MG_TETRA* tet=NULL;
    unsigned long id=no->get_id();
    std::multimap<unsigned long,MG_TETRA*,std::less <unsigned long> >::iterator it=lst.find(id);
    if (it!=lst.end())
    {
        tet=(*it).second;
        lst.erase(it);
    }
    return tet;
}

int REMAILLEUR::point_appartient_volume(double *xyz,int numvol)
{
    int ok=0;
    int i=0;
    std::multimap<double,double,std::less<double> > intersection;
    OT_VECTEUR_3D vecteur_dir;
    while (ok==0)

    {
        int ok2=0;
        while (ok2==0)
        {
            MG_TRIANGLE* tri=lsttrifront[numvol].get(i);
            MG_NOEUD* no1=tri->get_noeud1();
            MG_NOEUD* no2=tri->get_noeud2();
            MG_NOEUD* no3=tri->get_noeud3();
            double xb=(no1->get_x()+no2->get_x()+no3->get_x())/3.;
            double yb=(no1->get_y()+no2->get_y()+no3->get_y())/3.;
            double zb=(no1->get_z()+no2->get_z()+no3->get_z())/3.;
            OT_VECTEUR_3D directeur(xb-xyz[0],yb-xyz[1],zb-xyz[2]);
            directeur.norme();
            vecteur_dir=directeur;
            OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
            OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
            OT_VECTEUR_3D n=n1n3&n1n2;
            n.norme();
            double ps1=n*directeur;
            if (fabs(ps1)<1e-6) i++;
            else ok2=1;
        }
        ok2=0;
        intersection.clear();
        for (int j=0;j<lsttrifront[numvol].get_nb();j++)
        {
            MG_TRIANGLE* tri=lsttrifront[numvol].get(j);

            double t;
            int res=inter_droite_triangle(xyz,vecteur_dir,tri,&t);
            if (res==2) {
                ok2=1;
                break;
            }
            if (res==1)
            {
                std::pair<double,double> tmp(t,t);
                intersection.insert(tmp);
            }
        }
        if (ok2==0) ok=1;
        i++;

    }
    std::multimap<double,double,std::less <double> >::iterator it=intersection.begin();
    int nbinterneg=0;
    double tavant=1e308;
    while (it!=intersection.end())
    {
        if ((*it).second<0.)
        {
            double t=(*it).second;
            if (!(OPERATEUR::egal(t,tavant,1e-6)))
                nbinterneg++;
            tavant=t;
        }
        it++;
    }
    if (nbinterneg%2==1) return 1;
    return 0;
}


int REMAILLEUR::inter_droite_triangle(double *xyz,double *dir,MG_TRIANGLE* tri,double *t)
{
    int inter=0;
    MG_NOEUD *noeud1=tri->get_noeud1();
    MG_NOEUD *noeud2=tri->get_noeud2();
    MG_NOEUD *noeud3=tri->get_noeud3();
    OT_VECTEUR_3D g1(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
    OT_VECTEUR_3D g2(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
    OT_VECTEUR_3D g3(dir);
    OT_MATRICE_3D systeme(g1,g2,g3);
    double det=systeme.get_determinant();
    double eps=0.333333333333*(g1.diff()+g2.diff()+g3.diff());
    eps=eps*eps*0.0018;
    if (OPERATEUR::egal(det,0.0,eps)==true)
    {
        OT_VECTEUR_3D g3b(xyz[0]-noeud1->get_x(),xyz[1]-noeud1->get_y(),xyz[2]-noeud1->get_z());
        OT_MATRICE_3D systeme2(g1,g2,g3b);
        double det2=systeme2.get_determinant();
        double eps2=0.333333333333*(g1.diff()+g2.diff()+g3b.diff());
        eps2=18.*eps2*eps2*eps2*1e-6;
        if (OPERATEUR::egal(det2,0.0,eps)==true)    // cas 2D
        {
            return 2;
        }
        else return 0;
    }
    else
    {
        double x,y,z;
        x=1.0/det*(g2.get_y()*g3.get_z()-g2.get_z()*g3.get_y());
        y=1.0/det*(g3.get_x()*g2.get_z()-g2.get_x()*g3.get_z());
        z=1.0/det*(g2.get_x()*g3.get_y()-g2.get_y()*g3.get_x());
        OT_VECTEUR_3D g1b(x,y,z);
        x=1.0/det*(g3.get_y()*g1.get_z()-g1.get_y()*g3.get_z());
        y=1.0/det*(g1.get_x()*g3.get_z()-g3.get_x()*g1.get_z());
        z=1.0/det*(g3.get_x()*g1.get_y()-g1.get_x()*g3.get_y());
        OT_VECTEUR_3D g2b(x,y,z);
        x=1.0/det*(g1.get_y()*g2.get_z()-g1.get_z()*g2.get_y());
        y=1.0/det*(g2.get_x()*g1.get_z()-g1.get_x()*g2.get_z());
        z=1.0/det*(g1.get_x()*g2.get_y()-g1.get_y()*g2.get_x());
        OT_VECTEUR_3D g3b(x,y,z);
        OT_VECTEUR_3D n1n4(xyz[0]-noeud1->get_x(),xyz[1]-noeud1->get_y(),xyz[2]-noeud1->get_z());
        double alpha1=n1n4*g1b;
        double alpha2=n1n4*g2b;
        double alpha3=-(n1n4*g3b);
        double alpha4=1-alpha1-alpha2;
        double eps=0.000001;
        if ((alpha1>-eps) && (alpha1<1.+eps))
            if ((alpha2>-eps) && (alpha2<1.+eps))
                if ((alpha4>-eps) && (alpha4<1.+eps)) {
                    *t=alpha3;
                    return 1;
                }
        return 0;
    }

}
Revision:	951
Committed:	Fri Aug 10 15:17:17 2018 UTC (6 years, 9 months ago) by couturad
Original Path:	*magic/lib/mailleur_auto/src/remailleur.cpp*
File size:	54696 byte(s)
Log Message:	-> Ajout de Project Chrono (voir CMakeLists.txt). -> Ajout d'un générateur de microstructure basé sur la dynamique des corps rigides (MSTRUCT_GENERATEUR_DCR). -> Ajout d'un opérateur de décallage de la topologie (MG_CG_OP_TRANSF_DECALLAGE). -> Retrait de «using namespace std» (conflit avec namespace chrono) et modification des fichiers affectés. -> Modification de mailleur2d.cpp afin d'enregistrer un fichier MAGiC (void.magic) lorsque le nombre d'itération dépasse la valeur maximale.