CAD4FE_mailleur2d_outil.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       CAD4FE_mailleur2d_outil.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:56 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
#include <vector>
using namespace std;
#include "gestionversion.h"
#include <math.h>
#include "CAD4FE_mailleur2d.h"
#include "ot_mathematique.h"    
#include "ot_boite_3d.h"
#include "ot_boite_2d.h"
#include "tpl_map_entite.h"
#include "CAD4FE_m3d_MCTriangle.h"
 
#include "CAD4FE_Intersection_MCSegment_MCSegment.h"
#include "CAD4FE_MCNodePolyline.h"
#include "CAD4FE_MCSegment_Middle.h"
#include "ot_algorithme_geometrique.h"
#include "CAD4FE_Intersection_Plane_PolySurface.h"    
#include "CAD4FE_PolySurface.h"
#include "CAD4FE_Geometric_Tools.h"
#include "CAD4FE_ShortestPathByClosestPointOnEdge.h"
#include "CAD4FE_InventorText_MG_MAILLAGE.h"
#include "CAD4FE_InventorText_MCSegment.h"
#include "CAD4FE_FaceBoundaryPoint.h"
#include "tpl_set.h"
#include <fstream>
 
 
using namespace CAD4FE;
 
 
int MAILLEUR2D::noeud_est_dans_triangle(MCNode* noeud,MG_TRIANGLE *triangle)
{
return noeud_est_dans_triangle(noeud,(MCNode*)triangle->get_noeud1(),(MCNode*)triangle->get_noeud2(),(MCNode*)triangle->get_noeud3());
}
 
/*
int MAILLEUR2D::noeud_est_dans_triangle(MG_NOEUD* noeud,MG_NOEUD *noeud1,MG_NOEUD *noeud2,MG_NOEUD *noeud3)
{
double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
 
double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
double u3=decalage->decalage_parametre_u(noeud3->get_u(),du);
double v3=decalage->decalage_parametre_v(noeud3->get_v(),dv);
double u=decalage->decalage_parametre_u(noeud->get_u(),du);
double v=decalage->decalage_parametre_v(noeud->get_v(),dv);
double eps=0.0001;
double delta=(u2-u1)*(v3-v1)-(v2-v1)*(u3-u1);
double precision=std::max(fabs(u1),fabs(v1));
precision=std::max(precision,fabs(u2));
precision=std::max(precision,fabs(v2));
precision=std::max(precision,fabs(u3));
precision=std::max(precision,fabs(v3));
precision=std::max(precision,fabs(u));
precision=std::max(precision,fabs(v));
if (OPERATEUR::egal(delta,0.0,precision*eps))  return(0);
double xsi=1.0/delta*((v3-v1)*(u-u1)-(u3-u1)*(v-v1));
double eta=1.0/delta*((u2-u1)*(v-v1)-(v2-v1)*(u-u1));
if (!((eta > eps) && (xsi > eps) && ((eta+xsi) < 1.0-eps)))  return(0);
return (1);
}      */
 
 
int MAILLEUR2D::noeud_est_dans_triangle(MCNode* noeud,MCNode *noeud1,MCNode *noeud2,MCNode *noeud3)
{
    // (1) la projection normale du noeud sur le triangle existe
 
    // the 3 triangle's vertices
    OT_VECTEUR_3D TV0(noeud1->get_coord());
    OT_VECTEUR_3D TV1(noeud2->get_coord());
    OT_VECTEUR_3D TV2(noeud3->get_coord());
 
    int res_projectionNoeudTriangle = 0;
    double pt_projectionNoeudTriangle[3];
    double signedDist_projectionNoeudTriangle;
    res_projectionNoeudTriangle = OT_ALGORITHME_GEOMETRIQUE::project_PointTriangle(noeud->get_coord(),TV0,TV1,TV2,pt_projectionNoeudTriangle,&signedDist_projectionNoeudTriangle);
 
    if (res_projectionNoeudTriangle == 0 // disjoint
        )
        return 0;
 
    OT_VECTEUR_3D TV0TV1=TV1-TV0;
    OT_VECTEUR_3D TV1TV2=TV2-TV1;  
    OT_VECTEUR_3D TV2TV0=TV0-TV2;
    double x = (TV0TV1&TV1TV2).get_longueur();
    x  /= (TV0TV1.get_longueur()+TV1TV2.get_longueur()+TV2TV0.get_longueur());
 
    if (fabs(signedDist_projectionNoeudTriangle) > x)
        return 0;
 
    return 1;
}
 
int MAILLEUR2D::insere_segment(MCFace *__mcFace,MCSegment **nv_segment,MCNode* noeud1,MCNode* noeud2,int type_validation, double __tolerance)
{
OT_VECTEUR_3D noeud1_xyz(noeud1->get_coord());
OT_VECTEUR_3D noeud2_xyz(noeud2->get_coord());
OT_VECTEUR_3D milieu_segment_xyz = .5*noeud1_xyz + .5*noeud2_xyz;
double longueur = (noeud2_xyz - noeud1_xyz).get_longueur();
 
MCSegment * segment_candidat = new MCSegment(__mcFace,noeud1,noeud2,0);
/*BOITE_3D boite_candidat = segment_candidat->get_boite_3D();
{
    double offset=.5*segment_candidat->get_longueur();
    double xmin=boite_candidat.get_xmin()-offset;
    double ymin=boite_candidat.get_ymin()-offset;
    double zmin=boite_candidat.get_zmin()-offset;
    double xmax=boite_candidat.get_xmax()+offset;
    double ymax=boite_candidat.get_ymax()+offset;
    double zmax=boite_candidat.get_zmax()+offset;
    boite_candidat.reinit(xmin,ymin,zmin,xmax,ymax,zmax);
}  */
 
 
if (type_validation==TOUS_SEGMENT )
    {
    TPL_MAP_ENTITE<MCSegment*> liste_trouvee;   
    TPL_MAP_ENTITE<MCSegment*>::ITERATEUR it_seg;
    grille_de_segment->rechercher(milieu_segment_xyz(0),milieu_segment_xyz(1),milieu_segment_xyz(2),longueur,liste_trouvee);
    for (MCSegment* segment_trouve=liste_trouvee.get_premier(it_seg); segment_trouve; segment_trouve=liste_trouvee.get_suivant(it_seg))
       {
       if ( ((noeud1==segment_trouve->get_noeud1()) && (noeud2==segment_trouve->get_noeud2())) || ((noeud1==segment_trouve->get_noeud2()) && (noeud2==segment_trouve->get_noeud1())))
        {
            delete segment_candidat;
            return 0;
        }
        
        if (   (noeud1==segment_trouve->get_noeud1()&&noeud2!=segment_trouve->get_noeud2())
            || (noeud1==segment_trouve->get_noeud2()&&noeud2!=segment_trouve->get_noeud1())
            || (noeud2==segment_trouve->get_noeud2()&&noeud1!=segment_trouve->get_noeud1())
            || (noeud2==segment_trouve->get_noeud1()&&noeud1!=segment_trouve->get_noeud2()) )
            continue;
 
       /* BOITE_3D boite_trouve = segment_trouve->get_boite_3D();
       if (boite_candidat.intersecte(boite_trouve) == 0)
            continue;*/
 
       double uvIntersectionFace[2];
       MG_FACE * intersectionFace;
 
       // des intersections : voir bugs et livre Maillages George&Frey page 226 la remarque 6.11
       if ( Intersection_MCSegment_MCSegment_tolerance(segment_trouve, segment_candidat, uvIntersectionFace, &intersectionFace, __tolerance*segment_trouve->get_longueur()) )
       {
            delete segment_candidat;
            return 0;
       }
 
       }
     }
if (type_validation==TOUS_FRONT)
     {
     TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
     TPL_MAP_ENTITE<MG_FRONT_2D*>::ITERATEUR it_ft;
     grille_de_front->rechercher(milieu_segment_xyz(0),milieu_segment_xyz(1),milieu_segment_xyz(2),longueur,liste_trouvee);
     for (MG_FRONT_2D* ft=liste_trouvee.get_premier(it_ft);ft;ft=liste_trouvee.get_suivant(it_ft))
        {
        MCSegment* segment_trouve=ft->get_segment();
        if ( ((noeud1==segment_trouve->get_noeud1()) && (noeud2==segment_trouve->get_noeud2())) || ((noeud1==segment_trouve->get_noeud2()) && (noeud2==segment_trouve->get_noeud1())))
        {
            delete segment_candidat;
            return 0;
        }
 
        if ( (noeud1==segment_trouve->get_noeud1()&&noeud2!=segment_trouve->get_noeud2())
            || (noeud1==segment_trouve->get_noeud2()&&noeud2!=segment_trouve->get_noeud1())
            || (noeud2==segment_trouve->get_noeud2()&&noeud1!=segment_trouve->get_noeud1())
            || (noeud2==segment_trouve->get_noeud1()&&noeud1!=segment_trouve->get_noeud2()) )
            continue;
                               
       /* BOITE_3D boite_trouve = segment_trouve->get_boite_3D();
       if (boite_candidat.intersecte(boite_trouve) == 0)
            continue;*/
            
        double uvIntersectionFace[2];
        MG_FACE * intersectionFace;
       // des intersections : voir bugs et livre Maillages George&Frey page 226 la remarque 6.11
        if ( Intersection_MCSegment_MCSegment_tolerance(segment_trouve, segment_candidat, uvIntersectionFace, &intersectionFace, __tolerance*segment_trouve->get_longueur()) )
        {
            delete segment_candidat;
            return 0;
        }
 
        }
      }
 
mg_maillage->ajouter_mg_segment(segment_candidat);
grille_de_segment->inserer(segment_candidat);
*nv_segment=segment_candidat;
return(1);
}
 
 
void MAILLEUR2D::supprime_segment(MCSegment* __segment)
{
grille_de_segment->supprimer(__segment);
mg_maillage->supprimer_mg_segmentid(__segment->get_id());
}
 
int MAILLEUR2D::genere_noeud(MCFace* __mcFace, MCNode *__interiorIsolatedNode, MG_FRONT_2D **front_rencontre,MCNode **noeud_solution)
{
MCNode * iiNode = __interiorIsolatedNode;
 
double angle = M_PI * 44/180;
int nbRefFaces;
OT_VECTEUR_3D normale_face;
__mcFace->calcul_normale_unitaire(iiNode->GetRefFaceMapping(), normale_face, &nbRefFaces);
OT_MATRICE_3D planeFrame = GeometricTools::GetPlaneFrame(normale_face);
OT_VECTEUR_3D direction_placement = cos(angle)*planeFrame.get_vecteur1()+sin(angle)*planeFrame.get_vecteur2();
OT_VECTEUR_3D normale_plan_directeur = direction_placement & normale_face;
 
double taille;        
double density_tensor[9];
metrique->evaluer(iiNode->get_coord(), density_tensor);
taille = pow(density_tensor[0],-0.5);
double longueur_desiree=0.8660254037844386*taille;
 
Intersection_Plane_PolySurface ot_trajectoire_placement(iiNode, normale_plan_directeur, __mcFace, NULL);
ot_trajectoire_placement.MakeOffset(iiNode, direction_placement, longueur_desiree);
MCNode * noeud_decale = ot_trajectoire_placement.GetEndPoint();
OT_VECTEUR_3D xdecale(noeud_decale->get_coord());
 
TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
double longueur_recherche=1.1547*longueur_desiree;
OT_VECTEUR_3D normale_face_front=normale_face;
 
grille_de_front->rechercher(xdecale.get_x(),xdecale.get_y(),xdecale.get_z(),longueur_recherche,liste_trouvee);
TPL_MAP_ENTITE < MG_FRONT_2D * >::ITERATEUR it_liste_trouvee;
 
std::multimap <double, MCNode *> map_noeud_trouve;
std::map <MCNode *, MG_FRONT_2D*> map_front_trouve;
 for ( MG_FRONT_2D* front_courant=liste_trouvee.get_premier(it_liste_trouvee);
       front_courant;
       front_courant=liste_trouvee.get_suivant(it_liste_trouvee)
       )
 {                                
 
 MCNode* noeud_front1=(MCNode*)front_courant->get_noeud1();
 MCNode* noeud_front2=front_courant->get_noeud2();
 OT_VECTEUR_3D xnoeudfront1(noeud_front1->get_coord());
 OT_VECTEUR_3D xnoeudfront2(noeud_front2->get_coord());
 
    double distanceEuclidMPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xnoeudfront1, iiNode->get_coord());
    double distanceEuclidPcPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xnoeudfront1, xdecale);
    if ( distanceEuclidPcPf <= 1.1547*longueur_desiree && distanceEuclidMPf <= 1.5*longueur_desiree)
    if (map_front_trouve.find(noeud_front1) == map_front_trouve.end())
     {
     refresh();
     map_noeud_trouve.insert(std::make_pair(distanceEuclidPcPf,noeud_front1));
     map_front_trouve[noeud_front1]=front_courant;
     }
 }
 
int nb_front_essaye=0;
for (std::multimap<double,MCNode*>::iterator itMCNode = map_noeud_trouve.begin();
    itMCNode != map_noeud_trouve.end();
    itMCNode++ )
{
   double distance_reference = itMCNode->first;
   MCNode * noeud = itMCNode->second; 
   MG_FRONT_2D * front_reference = map_front_trouve[noeud];
 
   if (distance_reference > longueur_desiree*1.1547) break;
   nb_front_essaye++;
 
    double distanceEuclidMPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(noeud->get_coord(), iiNode->get_coord());
    double distanceEuclidPcPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(noeud->get_coord(), noeud_decale->get_coord());
 
    if ( distanceEuclidPcPf > 1.1547*longueur_desiree || distanceEuclidMPf > 1.5*longueur_desiree)
        continue;
 
    ShortestPathByClosestPointOnEdge operateur_dist_MPf(__mcFace,noeud,iiNode,distanceEuclidMPf);
    double distanceMPf = operateur_dist_MPf.Find();
    if ( distanceMPf > 1.5*longueur_desiree )
        continue;
 
    ShortestPathByClosestPointOnEdge operateur_dist_PcPf(__mcFace,noeud,iiNode,distanceEuclidPcPf);
    double distancePcPf = operateur_dist_PcPf.Find();
    if ( distancePcPf > 1.1547*longueur_desiree )
        continue;
 
    (*front_rencontre)=front_reference;
    (*noeud_solution)=noeud;
    
    if (debug) // Begin Log
    {
        std::ostringstream out;
        out << "genere_noeud: found front segment " << front_reference->get_segment()->get_id() << " and connected isolated interior node " << __interiorIsolatedNode->get_id() << " to node " << noeud->get_id();
        WriteLog(out.str());
    } // End log
 
    return FRONT_RENCONTRE;
}
 
mg_maillage->ajouter_mg_noeud(noeud_decale);
if (debug) // Begin Log
{
    std::ostringstream out;
    out << "genere_noeud: Node " << noeud_decale->get_id() << " created from interior isolated node " << __interiorIsolatedNode->get_id();
    WriteLog(out.str()); 
} // End log
noeud_decale->incrementer();  
noeud_decale->Creator = MCNode::mailleur_2d;
(*noeud_solution)=noeud_decale;
(*front_rencontre)=NULL;
return NOEUD_CREE;
}
 
int MAILLEUR2D::genere_noeud(MCFace* __mcFace,MG_FRONT_2D* front,MG_FRONT_2D **front_rencontre,MCNode **noeud_rencontre)
{
MCNode* noeud1=front->get_noeud1();
MCNode* noeud2=front->get_noeud2();
 
MCNode * milieu_front=MCSegment_Middle(front->get_segment());
OT_VECTEUR_3D x1(noeud1->get_coord());
OT_VECTEUR_3D x2(noeud2->get_coord());
OT_VECTEUR_3D normale_plan_directeur = x2-x1;
OT_VECTEUR_3D normale_face(0,0,0);
int nbRefFaces;
__mcFace->calcul_normale_unitaire(milieu_front->GetRefFaceMapping(), normale_face, &nbRefFaces);
/*if (nbRefFaces == 0)
{
    double minDist=1E307;
    OT_VECTEUR_3D segMid=.5*(x1+x2);
    double minT=1E300;
    for (double t=0;t<=1;t+=.1)
    {                         
        delete milieu_front;
        milieu_front = new MCNode;
        front->get_segment()->evaluer_geo(t,milieu_front);
        int nbRefFacesTmp = 0;
        for (MCNode::FMapIterator itF = milieu_front->GetRefFaceMapping().begin();
            itF != milieu_front->GetRefFaceMapping().end();
            itF++)
        {
            MG_FACE * f = itF->first; 
            if (__mcFace->GetPolySurface()->Contains(f))
                nbRefFacesTmp++;
        }
        if (nbRefFacesTmp == 0)
        {
            continue;
        }
        OT_VECTEUR_3D xTmp(milieu_front->get_coord());
        double dist = (segMid-xTmp).get_longueur();
        if (dist<minDist)
        {
            minT=t;
            minDist=dist;
        }
    }          
    delete milieu_front;
    if (minDist < 1E9)
    {
        front->get_segment()->evaluer_geo(minT,milieu_front);
    }
}     */
OT_VECTEUR_3D direction_placement=normale_face&normale_plan_directeur;
 
double longueur_segment=front->get_segment()->get_longueur();
double taille;        
double density_tensor[9];
metrique->evaluer(milieu_front->get_coord(), density_tensor);
taille = pow(density_tensor[0],-0.5);
double longueur_desiree=0.8660254037844386*(MAILLEUR2D::priorite_metrique*taille+longueur_segment*(1-MAILLEUR2D::priorite_metrique));
longueur_desiree = longueur_desiree/front->get_ifail();
 
Intersection_Plane_PolySurface ot_trajectoire_placement(milieu_front,normale_plan_directeur,__mcFace,NULL);
ot_trajectoire_placement.MakeOffset(milieu_front, direction_placement, longueur_desiree);
MCNode * noeud_decale = ot_trajectoire_placement.GetEndPoint();
OT_VECTEUR_3D xdecale(noeud_decale->get_coord());
 
TPL_MAP_ENTITE<MG_FRONT_2D*> liste_trouvee;
double longueur_recherche=1.1547*longueur_desiree;
 
OT_VECTEUR_3D normale_face_front=normale_face;
 
grille_de_front->rechercher(xdecale.get_x(),xdecale.get_y(),xdecale.get_z(),longueur_recherche,liste_trouvee);
 TPL_MAP_ENTITE < MG_FRONT_2D * >::ITERATEUR it_liste_trouvee;
 
std::multimap <double, MCNode *> map_noeud_trouve;
std::map <MCNode *, MG_FRONT_2D*> map_front_trouve;
std::map <MCNode *, double> map_angle_front_trouve;
 for ( MG_FRONT_2D* front_courant=liste_trouvee.get_premier(it_liste_trouvee);
       front_courant;
       front_courant=liste_trouvee.get_suivant(it_liste_trouvee)
       )
 {                                
 
 MCNode* noeud_front1=(MCNode*)front_courant->get_noeud1();
 MCNode* noeud_front2=front_courant->get_noeud2();
 OT_VECTEUR_3D xnoeudfront1(noeud_front1->get_coord());
 OT_VECTEUR_3D xnoeudfront2(noeud_front2->get_coord());
 
    double distanceEuclidMPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xnoeudfront1, milieu_front->get_coord());
    double distanceEuclidPcPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xnoeudfront1, xdecale);
    if ( distanceEuclidPcPf <= 1.1547*longueur_desiree && distanceEuclidMPf <= 1.5*longueur_desiree)
    if (noeud_front1->get_id()!=noeud1->get_id())
    if (noeud_front1->get_id()!=noeud2->get_id())
    if ( !(noeud_front1->get_id()==noeud2->get_id() && noeud_front2->get_id()==noeud1->get_id()) ) 
    if (map_front_trouve.find(noeud_front1) != map_front_trouve.end())
     {
     refresh();
     /*{int nbRefFaces;
     __mcFace->calcul_normale_unitaire(noeud_front1->GetRefFaceMapping(), normale_face_front, &nbRefFaces);}*/
     noeud_front1->NormalMCFace(__mcFace,normale_face_front);
     //OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
     OT_VECTEUR_3D vecteur_baseu = x1 - xnoeudfront1;
     //OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
     OT_VECTEUR_3D vecteur_front =  xnoeudfront2 - xnoeudfront1;
     //     OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
     OT_VECTEUR_3D vecteur_basev = normale_face_front & vecteur_baseu;
     vecteur_baseu.norme();
     vecteur_basev.norme();
     vecteur_front.norme();
     double cosangle=vecteur_baseu*vecteur_front;
     double sinangle=vecteur_basev*vecteur_front;
     sinangle=-sinangle;
     if (cosangle>1.) cosangle=1.;
     if (cosangle<-1.) cosangle=(-1.);
     double angle=acos(cosangle);
     if (sinangle<(-0.0001)) angle=(-angle);
     if (angle<0.) angle=angle+2.*M_PI;
     if (angle<map_angle_front_trouve[noeud_front1])
         {
         map_angle_front_trouve[noeud_front1]=angle;  
         map_front_trouve[noeud_front1]=front_courant;
         }
     }
  else
     {
     refresh();
     /*{int nbRefFaces;
     __mcFace->calcul_normale_unitaire(noeud_front1->GetRefFaceMapping(), normale_face_front, &nbRefFaces);}*/
     noeud_front1->NormalMCFace(__mcFace,normale_face_front);
     // OT_VECTEUR_3D vecteur_baseu(u1-unoeudfront1,v1-vnoeudfront1,0.);
     OT_VECTEUR_3D vecteur_baseu = x1 - xnoeudfront1;
     // OT_VECTEUR_3D vecteur_front(unoeudfront2-unoeudfront1,vnoeudfront2-vnoeudfront1,0.);
     OT_VECTEUR_3D vecteur_front = xnoeudfront2 - xnoeudfront1;
     // OT_VECTEUR_3D vecteur_basev=w&vecteur_baseu;
     OT_VECTEUR_3D vecteur_basev = normale_face_front & vecteur_baseu;
     vecteur_baseu.norme();
     vecteur_basev.norme();
     vecteur_front.norme();
     double cosangle=vecteur_baseu*vecteur_front;
     double sinangle=vecteur_basev*vecteur_front;
     sinangle=-sinangle;
     if (cosangle>1.) cosangle=1.;
     if (cosangle<-1.) cosangle=(-1.);
     double angle=acos(cosangle);
     if (sinangle<(-0.0001)) angle=(-angle);
     if (angle<0.) angle=angle+2.*M_PI;
     map_noeud_trouve.insert(std::make_pair(distanceEuclidMPf,noeud_front1));
     map_front_trouve[noeud_front1]=front_courant;
     map_angle_front_trouve[noeud_front1]=angle;
     }
 }
 
 
 
 
int nb_front_essaye=0;
for (std::multimap<double,MCNode*>::iterator itMCNode = map_noeud_trouve.begin();
    itMCNode != map_noeud_trouve.end();
    itMCNode++ )
{
   double distance_reference = itMCNode->first;
   MCNode * noeud = itMCNode->second; 
   MG_FRONT_2D * front_reference = map_front_trouve[noeud];
 
   if (distance_reference > 1.1547*longueur_desiree) break;
   nb_front_essaye++;
 
    double distanceEuclidMPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(noeud->get_coord(), milieu_front->get_coord());
    double distanceEuclidPcPf = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(noeud->get_coord(), noeud_decale->get_coord());
 
    if ( distanceEuclidPcPf > 1.1547*longueur_desiree || distanceEuclidMPf > 1.5*longueur_desiree)
        continue;
 
    std::vector< MCNode * > shortestPathNodes;
    std::vector< MG_ELEMENT_TOPOLOGIQUE * > shortestPathTopo;
 
    ShortestPathByClosestPointOnEdge operateur_dist_MPf(__mcFace,noeud,milieu_front,distanceEuclidMPf);
    double distanceMPf = operateur_dist_MPf.Find(&shortestPathNodes, &shortestPathTopo);
    if ( distanceMPf > 2 * distanceEuclidMPf )
        continue;
 
    // verify that the shortest path M-Pf does not intersect an interior MC Edge
    bool test_operateur_dist_MPf_Intersect_MCEdge = false;
    {
        for (int k=1; k+1<shortestPathNodes.size(); k++)
        {
            MCNode * mcNode_tmp = shortestPathNodes[k];
            if (mcNode_tmp->get_lien_topologie()->get_dimension() == 1)
            {
                test_operateur_dist_MPf_Intersect_MCEdge = true;
                break;
            }
        }
    }
    if (test_operateur_dist_MPf_Intersect_MCEdge)
        continue;
 
    ShortestPathByClosestPointOnEdge operateur_dist_PcPf(__mcFace,noeud,noeud_decale,distanceEuclidPcPf);
    double distancePcPf = operateur_dist_PcPf.Find();
    if ( distancePcPf > 1.1547*longueur_desiree )
        continue;
 
    if ( triangles_sont_dans_meme_sens( __mcFace, noeud1, noeud, noeud_decale, noeud2 ) == 0)
        continue;
 
    {
        MG_SEGMENT * segment1 = mg_maillage->get_mg_segment(noeud1->get_id(),noeud->get_id());
        if ( segment1 && front_reference->get_segment() != segment1 && segment1->get_lien_triangle()->get_nb() == 2 )
            continue;
        MG_SEGMENT * segment2 = mg_maillage->get_mg_segment(noeud2->get_id(),noeud->get_id());
        if ( segment2 && front_reference->get_segment() != segment2 && segment2->get_lien_triangle()->get_nb() == 2 )
            continue;
    }
 
    (*front_rencontre)=front_reference;
    (*noeud_rencontre)=noeud;
    if (debug) // Begin Log
    {
        std::ostringstream out;
        out << "genere_noeud: found front [ " << front_reference->get_noeud1()->get_id() << " - " << front_reference->get_noeud2()->get_id() << " ] and from candidate front [ " << front->get_segment()->get_noeud1()->get_id() << " - " << front->get_segment()->get_noeud2()->get_id() << " ] " << " to segment [ " << front->get_segment()->get_noeud1()->get_id() << " - " << front->get_segment()->get_noeud2()->get_id() << " ] ";
        WriteLog(out.str());
    } // End log
    return FRONT_RENCONTRE;
}
mg_maillage->ajouter_mg_noeud(noeud_decale);
noeud_decale->incrementer();
noeud_decale->Creator = MCNode::mailleur_2d;
if (debug) // Begin Log
{
    std::ostringstream out;
    out << "genere_noeud: created node " << noeud_decale->get_id() << " from candidate front [ " << front->get_noeud1()->get_id() << " - " << front->get_noeud2()->get_id() << " ] ";
    WriteLog(out.str());
} // End log
(*noeud_rencontre)=noeud_decale;
(*front_rencontre)=NULL;
return NOEUD_CREE;
}
 
/*
double MAILLEUR2D::calcule_longueur_segment_metrique(MG_FACE* mgface,MG_SEGMENT* mgsegment)
{
MG_NOEUD* noeud1=mgsegment->get_noeud1();
MG_NOEUD* noeud2=mgsegment->get_noeud2();
 
double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
double u1=decalage->decalage_parametre_u(noeud1->get_u(),du);
double v1=decalage->decalage_parametre_v(noeud1->get_v(),dv);
double u2=decalage->decalage_parametre_u(noeud2->get_u(),du);
double v2=decalage->decalage_parametre_v(noeud2->get_v(),dv);
return calcule_distance_metrique(mgface,u1,v1,u2,v2,du,dv);
}           */
 
 
double MAILLEUR2D::calcule_distance_metrique(MG_FACE* mgface,double u1,double v1,double u2,double v2,double du,double dv)
{
double longueur_calculee=0.;
double param_debut[2]={u1-du,v1-dv};
double E,F,G;
mgface->get_EFG(param_debut,E,F,G);
 
double dt=0.03125;
double ti;
double tii=0;
 
while (tii<1.)
        {
        refresh();
        ti=tii;
        tii=ti+dt;
        if (tii>1.) tii=1.;
        double t=0.7886751345*ti+0.2113248654*tii;
        double u=u1+t*(u2-u1);
        double v=v1+t*(v2-v1);
        double ut=u2-u1;
        double vt=v2-v1;
        double param_integration1[2]={u-du,v-dv};
        double tenseur_metrique[9];
        double xyzdu[3];
        double xyzdv[3];
        if (creation_metrique) metrique->evaluer(param_integration1,tenseur_metrique);
        else
                {
                double xyz[3];
                mgface->evaluer(param_integration1,xyz);
                metrique->evaluer(xyz,tenseur_metrique);
                }
        mgface->deriver(param_integration1,xyzdu,xyzdv);
        double t3 = xyzdu[0]*ut+xyzdv[0]*vt;
        double t7 = xyzdu[1]*ut+xyzdv[1]*vt;
        double t11 = xyzdu[2]*ut+xyzdv[2]*vt;
        double t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
        double t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
        double t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
        double facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
        longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
        t=0.7886751345*tii+0.2113248654*ti;
        u=u1+t*(u2-u1);
        v=v1+t*(v2-v1);
        double param_integration2[2]={u-du,v-dv};
        if (creation_metrique) metrique->evaluer(param_integration2,tenseur_metrique);
        else
                {
                double xyz[3];
                mgface->evaluer(param_integration2,xyz);
                metrique->evaluer(xyz,tenseur_metrique);
                }
        mgface->deriver(param_integration2,xyzdu,xyzdv);
        t3 = xyzdu[0]*ut+xyzdv[0]*vt;
        t7 = xyzdu[1]*ut+xyzdv[1]*vt;
        t11 = xyzdu[2]*ut+xyzdv[2]*vt;
        t13 = t3*tenseur_metrique[0]+t7*tenseur_metrique[3]+t11*tenseur_metrique[6];
        t18 = t3*tenseur_metrique[1]+t7*tenseur_metrique[4]+t11*tenseur_metrique[7];
        t23 = t3*tenseur_metrique[2]+t7*tenseur_metrique[5]+t11*tenseur_metrique[8];
        facteur = (t13*xyzdu[0]+t18*xyzdu[1]+t23*xyzdu[2])*ut+(t13*xyzdv[0]+t18*xyzdv[1]+t23*xyzdv[2])*vt;
        longueur_calculee=longueur_calculee+0.5*(tii-ti)*sqrt(facteur);
        }
return longueur_calculee;
}
 
 
MCTriangle* MAILLEUR2D::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo, MCNode *mgnoeud1, MCNode *mgnoeud2, MCNode *mgnoeud3)
{
int i=0;
MCSegment* mgsegment[3];
MCNode *mgnoeud[3]={mgnoeud1,mgnoeud2,mgnoeud3};
for (i=0;i<3;i++)
{
    mgsegment[i]=(MCSegment*)mg_maillage->get_mg_segment(mgnoeud[i]->get_id(),mgnoeud[(i+1)%3]->get_id());
    if (mgsegment[i]==NULL)
    {
        printf("Pb : le segment n'existe pas pour inserer le triangle!!!\n");
        mgsegment[i] = new MCSegment(topo,mgnoeud[i],mgnoeud[(i+1)%3]);
        mg_maillage->ajouter_mg_segment(mgsegment[i]);
    }
}
M3D_MCTriangle* mtriangle=new M3D_MCTriangle(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment[0],mgsegment[1],mgsegment[2]);
mg_maillage->ajouter_mg_triangle(mtriangle);
_nbTriangles[topo]++;
double qual=OPERATEUR::qualite_triangle(mgnoeud1->get_coord(),mgnoeud2->get_coord(),mgnoeud3->get_coord());
mtriangle->change_qualite(qual);
std::pair<const double,M3D_MCTriangle*> tmp(mtriangle->get_qualite(),mtriangle);
lst_tri_qual.insert(tmp);
return mtriangle;
}
 
 
int MAILLEUR2D::triangles_sont_dans_meme_sens(MCFace* face,MCNode* noeud1,MCNode *noeud2,MCNode *noeud2bis,MCNode *noeud3)
{
double* xyz1=noeud1->get_coord();
double* xyz2=noeud2->get_coord();  
double* xyz2bis=noeud2bis->get_coord();
double* xyz3=noeud3->get_coord();
OT_VECTEUR_3D n1n3(xyz1,xyz3);
OT_VECTEUR_3D n1n2(xyz1,xyz2);
OT_VECTEUR_3D n1n2bis(xyz1,xyz2bis);
 
OT_VECTEUR_3D normal1=n1n3&n1n2;
OT_VECTEUR_3D normal1bis=n1n3&n1n2bis;
 
if (normal1 * normal1bis < 0) return 0;
return 1;
}
 
int MAILLEUR2D::triangle_est_dans_bon_sens(MCFace* face,MCNode* noeud1,MCNode *noeud2,MCNode *noeud3)
{
double* xyz1=noeud1->get_coord();
double* xyz2=noeud2->get_coord();
double* xyz3=noeud3->get_coord();
OT_VECTEUR_3D n1n3(xyz1,xyz3);
OT_VECTEUR_3D n1n2(xyz1,xyz2);
 
OT_VECTEUR_3D normal1(0,0,0);
noeud1->NormalMCFace(face, normal1);
OT_VECTEUR_3D nmat=normal1&n1n3;
nmat.norme();
n1n2.norme();
if (nmat*n1n2<0) return 0;
return 1;
}
 
int MAILLEUR2D::bouge_point(MCFace* __mcFace, MCNode *mg_noeud, double& crit, MCNode **__result)
{
if (mg_noeud->Creator != MCNode::mailleur_2d)
    return 0;
 
    //double du=decalage->calcul_decalage_parametre_u(mg_noeud->get_u());
OT_VECTEUR_3D xyz1(mg_noeud->get_coord());
 
OT_VECTEUR_3D xyzopt(0,0,0);
double qual_dep=1.;
int nb_tri=mg_noeud->get_lien_triangle()->get_nb();
for (int i=0;i<nb_tri;i++)
{
        M3D_MCTriangle* tri=(M3D_MCTriangle*)mg_noeud->get_lien_triangle()->get(i);
        qual_dep=std::min(qual_dep,tri->get_qualite());
        MG_NOEUD* no1=tri->get_noeud1();
        MG_NOEUD* no2=tri->get_noeud2();
        MG_NOEUD* no3=tri->get_noeud3();
        MG_NOEUD *noeud1=NULL,*noeud2=NULL;
        if (no1==mg_noeud)
                {
                noeud1=no2;
                noeud2=no3;
                }
        else if (no2==mg_noeud)
                {
                noeud1=no1;
                noeud2=no3;
                }
        else if (no3==mg_noeud)
                {
                noeud1=no1;
                noeud2=no2;
                }
        OT_VECTEUR_3D noeud1_xyz(noeud1->get_coord());
        OT_VECTEUR_3D noeud2_xyz(noeud2->get_coord());
        xyzopt += noeud1_xyz;               
        xyzopt += noeud2_xyz;
        //double u2=noeud1->get_u()+du;
        //double v2=noeud1->get_v()+dv;
        //uopt=uopt+u2;
        //vopt=vopt+v2;
        //u2=noeud2->get_u()+du;
        //v2=noeud2->get_v()+dv;
        //uopt=uopt+u2;
        //vopt=vopt+v2;
        }
xyzopt /= 2*nb_tri;
 
double ddeb=(xyzopt-xyz1).get_longueur();
double alpha=0.5;
double d=alpha*ddeb;
double qual=qual_dep;
double testu[8]={1.,-1.,0.,0.,0.707106781,-0.707106781,-0.707106781,0.707106781};
double testv[8]={0.,0.,1.,-1.,0.707106781,0.707106781,-0.707106781,-0.707106781};
MCNode *xi=mg_noeud;     
MCNode *xii;
MCNode *x=NULL;
while (alpha>0.09)
        {
        double qualcoq=0.;        
        OT_VECTEUR_3D xyznormal;
        int nbRefFaces;
        __mcFace->calcul_normale_unitaire(xi->GetRefFaceMapping(),xyznormal,&nbRefFaces);
        if (xyznormal.get_longueur() == 0) return 0;
        OT_MATRICE_3D planeFrame=GeometricTools::GetPlaneFrame(xyznormal);
        for (int nb_essai=0;nb_essai<8;nb_essai++)
        {
              //double uv[2];
                OT_VECTEUR_3D uv(0,0,0);
                //double xyz[3];
                OT_VECTEUR_3D xyz(0,0,0);
                //uv[0]=uu+d*testu[nb_essai]-du;
                //uv[1]=vv+d*testv[nb_essai]-dv;
                uv[0]=testu[nb_essai];
                uv[1]=testv[nb_essai];
                //int valide1=mgface->valide_parametre_u(uv[0]);
                //int valide2=mgface->valide_parametre_v(uv[1]);
                //if (!((valide1) && (valide2) )) break;
                //mgface->evaluer(uv,xyz);
                OT_VECTEUR_3D direction_placement=planeFrame*uv;
                OT_VECTEUR_3D normale_plan_directeur=direction_placement&xyznormal;
                double longueur_desiree=d;
                Intersection_Plane_PolySurface ot_trajectoire_placement(xi,normale_plan_directeur,__mcFace,NULL);
                ot_trajectoire_placement.MakeOffset(xi, direction_placement, longueur_desiree);
                xii = ot_trajectoire_placement.GetEndPoint();
                xyz = xii->get_coord();
 
                double qual1=1.;
                for (int i=0;i<nb_tri;i++)
                        {
                        M3D_MCTriangle* tri=(M3D_MCTriangle*)mg_noeud->get_lien_triangle()->get(i);
                        MG_NOEUD* no1=tri->get_noeud1();
                        MG_NOEUD* no2=tri->get_noeud2();
                        MG_NOEUD* no3=tri->get_noeud3();
                        double *xyz1=no1->get_coord();
                        double *xyz2=no2->get_coord();
                        double *xyz3=no3->get_coord();
                        if (no1==mg_noeud) xyz1=xyz;
                        if (no2==mg_noeud) xyz2=xyz;
                        if (no3==mg_noeud) xyz3=xyz;
                        double qualtmp=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3);
                        OT_VECTEUR_3D n1n3(xyz1,xyz3);
                        OT_VECTEUR_3D n1n2(xyz1,xyz2);
                        //double xyznormal[3];
                        //mgface->calcul_normale_unitaire(uv,xyznormal);
                        OT_VECTEUR_3D normalface(xyznormal);
                        OT_VECTEUR_3D normal=normalface&n1n3;
                        normal.norme();
                        n1n2.norme();
                        if (normal*n1n2<0.0001) qualtmp=0.;
                        if (qualtmp<qual1) qual1=qualtmp;
                        }
                if (qual1>qualcoq)
                        {
                        qualcoq=qual1;
                        crit=qualcoq;
                        
                        if ( x && xi != x )
                            x->decrementer();
 
                        if ( x && xi != x && x->get_nb_reference() == 0 )
                            { delete x; x=0; }
 
                        x=xii;
                        x->incrementer();
                        }
                }
 
        if (qualcoq<qual+0.0001)
        {
            alpha=alpha-0.1;
        }
        else
        {
                qual=qualcoq;
                if (xi != mg_noeud) xi->decrementer();
                if (xi && xi != mg_noeud && xi->get_nb_reference() == 0 ) {delete xi; xi=0;}
                xi=x;
        }
        d=ddeb*alpha;
 
        }
if (qual>qual_dep)
        {
        //u=decalage->decalage_parametre_u(u,-du);
        //v=decalage->decalage_parametre_v(v,-dv);
        //double uv[2];
        //uv[0]=u;
        //uv[1]=v;
        //double xyz[3];
        //mgface->evaluer(uv,xyz);
        //x=xyz[0];
        //y=xyz[1];
        //z=xyz[2];
        *__result = x;
        x->Creator = MCNode::mailleur_2d_bouge_point;
        return 1;
        }
return 0;
}
 
int MAILLEUR2D::triangle_est_dans_bon_sens_front(MG_FRONT_2D * front,MCNode *noeud3)
{   
MCNode* noeud1=front->get_noeud1();
MCNode* noeud2=front->get_noeud2();
 
MCNode * milieu_front=MCSegment_Middle(front->get_segment());
OT_VECTEUR_3D x1(noeud1->get_coord());
OT_VECTEUR_3D x2(noeud2->get_coord());
OT_VECTEUR_3D x3(noeud3->get_coord());   
OT_VECTEUR_3D xm(milieu_front->get_coord());
OT_VECTEUR_3D normale_plan_directeur = x2-x1;
OT_VECTEUR_3D normale_face(0,0,0);
int nbRefFaces;
_mcFace->calcul_normale_unitaire(milieu_front->GetRefFaceMapping(), normale_face, &nbRefFaces);
OT_VECTEUR_3D direction_placement=normale_face&normale_plan_directeur;
OT_VECTEUR_3D direction_placement_euclidien = x3 - xm;
                 
double longueur_segment=front->get_segment()->get_longueur();
double taille;        
double density_tensor[9];
metrique->evaluer(milieu_front->get_coord(), density_tensor);
taille = pow(density_tensor[0],-0.5);
double longueur_desiree=0.8660254037844386*(MAILLEUR2D::priorite_metrique*taille+longueur_segment*(1-MAILLEUR2D::priorite_metrique));
longueur_desiree = longueur_desiree/front->get_ifail();
 
MCSegment * segment_M_N3 = new MCSegment(0, _mcFace, milieu_front, noeud3, 2 * longueur_desiree);
OT_VECTEUR_3D tangente_milieu, courbure_milieu;
 
std::vector <MG_ELEMENT_TOPOLOGIQUE*> topos = segment_M_N3->GetPolylineTopos();
std::vector <MG_ELEMENT_TOPOLOGIQUE*>::iterator topos_it;
MG_FACE * face1 = NULL;
for ( topos_it = topos.begin(); topos_it != topos.end(); topos_it ++)
{
    MG_ELEMENT_TOPOLOGIQUE* topo = *topos_it;
    if (topo->get_dimension() == 2)
    {
        face1 = (MG_FACE *) topo;
        break;
    }
}
 
int is_interior = 1;
if (face1)
{
    FaceBoundaryPoint p(milieu_front, face1);
    is_interior = p.TestInteriorDirection (direction_placement_euclidien);
}
 
delete segment_M_N3;
delete milieu_front;
 
if (is_interior == 1)
    return 1;
else
    return 0;
}
 
/*int MAILLEUR2D::bouge_point(MG_FACE* mgface,MG_NOEUD* mg_noeud,double& crit,double &u,double& v,double& x,double& y, double& z)
{
if (mg_noeud->get_lien_topologie()!=mgface) return 0;
double du=decalage->calcul_decalage_parametre_u(mg_noeud->get_u());
double dv=decalage->calcul_decalage_parametre_v(mg_noeud->get_v());
double u1=decalage->decalage_parametre_u(mg_noeud->get_u(),du);
double v1=decalage->decalage_parametre_v(mg_noeud->get_v(),dv);
 
double uopt=0.;
double vopt=0.;
double qual_dep=1.;
int nb_tri=mg_noeud->get_lien_triangle()->get_nb();
for (int i=0;i<nb_tri;i++)
        {
        M3D_MCTriangle* tri=(M3D_MCTriangle*)mg_noeud->get_lien_triangle()->get(i);
        qual_dep=std::min(qual_dep,tri->get_qualite());
        MG_NOEUD* no1=tri->get_noeud1();
        MG_NOEUD* no2=tri->get_noeud2();
        MG_NOEUD* no3=tri->get_noeud3();
        MG_NOEUD *noeud1,*noeud2;
        if (no1==mg_noeud)
                {
                noeud1=no2;
                noeud2=no3;
                }
        if (no2==mg_noeud)
                {
                noeud1=no1;
                noeud2=no3;
                }
        if (no3==mg_noeud)
                {
                noeud1=no1;
                noeud2=no2;
                }
        double u2=noeud1->get_u()+du;
        double v2=noeud1->get_v()+dv;
        uopt=uopt+u2;
        vopt=vopt+v2;
        u2=noeud2->get_u()+du;
        v2=noeud2->get_v()+dv;
        uopt=uopt+u2;
        vopt=vopt+v2;
        }
uopt=uopt/2./nb_tri;
vopt=vopt/2./nb_tri;
 
double ddeb=sqrt((u1-uopt)*(u1-uopt)+(v1-vopt)*(v1-vopt));
double alpha=0.5;
double d=alpha*ddeb;
double qual=qual_dep;
double testu[8]={1.,-1.,0.,0.,0.707106781,-0.707106781,-0.707106781,0.707106781};
double testv[8]={0.,0.,1.,-1.,0.707106781,0.707106781,-0.707106781,-0.707106781};
double uu=u1;
double vv=v1;
while (alpha>0.09)
        {
        double qualcoq=0.;
        for (int nb_essai=0;nb_essai<8;nb_essai++)
                {
                double uv[2];
                double xyz[3];
                uv[0]=uu+d*testu[nb_essai]-du;
                uv[1]=vv+d*testv[nb_essai]-dv;
                int valide1=mgface->valide_parametre_u(uv[0]);
                int valide2=mgface->valide_parametre_v(uv[1]);
                if (!((valide1) && (valide2) )) break;
                mgface->evaluer(uv,xyz);
                double qual1=1.;
                for (int i=0;i<nb_tri;i++)
                        {
                        M3D_MCTriangle* tri=(M3D_MCTriangle*)mg_noeud->get_lien_triangle()->get(i);
                        MG_NOEUD* no1=tri->get_noeud1();
                        MG_NOEUD* no2=tri->get_noeud2();
                        MG_NOEUD* no3=tri->get_noeud3();
                        double *xyz1=no1->get_coord();
                        double *xyz2=no2->get_coord();
                        double *xyz3=no3->get_coord();
                        if (no1==mg_noeud) xyz1=xyz;
                        if (no2==mg_noeud) xyz2=xyz;
                        if (no3==mg_noeud) xyz3=xyz;
                        double qualtmp=OPERATEUR::qualite_triangle(xyz1,xyz2,xyz3);
                        OT_VECTEUR_3D n1n3(xyz1,xyz3);
                        OT_VECTEUR_3D n1n2(xyz1,xyz2);
                        double xyznormal[3];
                        mgface->calcul_normale_unitaire(uv,xyznormal);
                        OT_VECTEUR_3D normalface(xyznormal);
                        OT_VECTEUR_3D normal=normalface&n1n3;
                        normal.norme();
                        n1n2.norme();
                        if (normal*n1n2<0.0001) qualtmp=0.;
                        if (qualtmp<qual1) qual1=qualtmp;
                        }
                if (qual1>qualcoq)
                        {
                        qualcoq=qual1;
                        crit=qualcoq;
                        u=uv[0]+du;
                        v=uv[1]+dv;
                        }
                }
 
        if (qualcoq<qual+0.0001) alpha=alpha-0.1;
        else
                {
                qual=qualcoq;
                uu=u;
                vv=v;
                }
        d=ddeb*alpha;
 
        }
if (qual>qual_dep)
        {
        u=decalage->decalage_parametre_u(u,-du);
        v=decalage->decalage_parametre_v(v,-dv);
        double uv[2];
        uv[0]=u;
        uv[1]=v;
        double xyz[3];
        mgface->evaluer(uv,xyz);
        x=xyz[0];
        y=xyz[1];
        z=xyz[2];
        return 1;
        }
return 0;
}     */
 
 
 
 
 
 
 
 
 
 
 
 
/*
 
 
 
#define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
#define EGAL(x,y,eps) (float)fabs((double)(x-y))<eps
#define DETER(a,b,c,d) (a*d-b*c)
 
int MAILLEUR2D::intersection_segment_segment(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
{
 
double du=decalage->calcul_decalage_parametre_u(noeud1->get_u());
double dv=decalage->calcul_decalage_parametre_v(noeud1->get_v());
double ua=decalage->decalage_parametre_u(noeud1->get_u(),du);
double va=decalage->decalage_parametre_v(noeud1->get_v(),dv);
double ub=decalage->decalage_parametre_u(noeud2->get_u(),du);
double vb=decalage->decalage_parametre_v(noeud2->get_v(),dv);
double um=decalage->decalage_parametre_u(noeud3->get_u(),du);
double vm=decalage->decalage_parametre_v(noeud3->get_v(),dv);
double un=decalage->decalage_parametre_u(noeud4->get_u(),du);
double vn=decalage->decalage_parametre_v(noeud4->get_v(),dv);
double ab[3];
double nm[3];
double am[3];
ab[0]=ub-ua;
ab[1]=vb-va;
ab[2]=0.;
nm[0]=um-un;
nm[1]=vm-vn;
nm[2]=0.;
am[0]=um-ua;
am[1]=vm-va;
am[2]=0.;
int equation[4];
equation[0]=1; // etat de l'equation 0 
equation[1]=1;
equation[2]=1;
equation[3]=3; // cette variable comporte le bilan du nombre d'equation
double eps2=PSCA(ab,ab);
double eps=sqrt(eps2);
eps=eps*0.0001;
eps2=eps2*0.0001;
if ( (EGAL(ab[0],0,eps)) && (EGAL(nm[0],0,eps)) )
        if (EGAL(am[0],0,eps)) equation[0]=0; else return(0);
if ( (EGAL(ab[1],0,eps)) && (EGAL(nm[1],0,eps)) )
        if (EGAL(am[1],0,eps)) equation[1]=0; else return(0);
if ( (EGAL(ab[2],0,eps)) && (EGAL(nm[2],0,eps)) )
        if (EGAL(am[2],0,eps)) equation[2]=0; else return(0);
equation[3]=equation[0]+equation[1]+equation[2]; 
if (equation[3]==3)
        {
        double det=DETER(ab[0],nm[0],ab[1],nm[1]);
        if (fabs(det)>eps2)  
                {
                det=1/det;
                double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
                double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
                if ( (float)fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps2) return(0);
                return(examine_solution(sol1,sol2,1));
                }
        else
                {
                equation[0]=0;
                equation[3]=2;
                // on verifie la compatibilite des deux equations dont le det est nul
                double tmp;
                if (ab[0]!=0) tmp=ab[1]*am[0]/ab[0]; else tmp=nm[1]*am[0]/nm[0];
                if (!(EGAL(tmp,am[1],eps))) return(0);
                }
        }
if (equation[3]==2)
        {
        // on repere les equations qui existent
        int ne1;
        int ne2;
        if (equation[0]!=0) 
                {
                ne1=0;
                if (equation[1]!=0) ne2=1; else ne2=2;
                }
        else 
                {
                ne1=1;
                ne2=2;
                }
        
        double det=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
        if (fabs(det)>eps2)  
                {
                det=1/det;
                double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
                double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
                return(examine_solution(sol1,sol2,1));
                }
        else
                {
                equation[ne1]=0;
                equation[3]=1;
                // on verifie la compatibilite des deux equations dont le det est nul 
                double tmp;
                if (ab[ne1]!=0) tmp=ab[ne2]*am[ne1]/ab[ne1]; else tmp=nm[ne2]*am[ne1]/nm[ne1];
                if (!(EGAL(tmp,am[ne2],eps))) return(0);
                }
        
        }
if (equation[3]==1)
        {
        // on repere l' equation qui existe
        int ne1;
        if (equation[0]!=0) ne1=0; else
        if (equation[1]!=0) ne1=1; else ne1=2;
        double an[3];
        an[0]=un-ua;
        an[1]=vn-va;
        an[2]=0.;
        double tmp=1/ab[ne1];
        double sol1=am[ne1]*tmp;
        double sol2=an[ne1]*tmp;
        return(examine_solution(sol1,sol2,2));
        }
return(0);
}
 
 
 
int MAILLEUR2D::examine_solution(double sol1,double sol2,int type)
{
double epsilon=0.0001;
 
if (type==1)
        {
        if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
        if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
        if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
        if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
        if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
 
        }
if (type==2)
        {
        if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
        if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
        if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
        if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
        }
return 0;
}
#undef EGAL
#undef PSCA
#undef DETER
 
*/