CAD4FE_Geometric_Tools.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_Geometric_Tools.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:56 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
#pragma hdrstop
 
#include "gestionversion.h"
#include <mg_face.h>
#include <ot_algorithme_geometrique.h>
#include <CAD4FE_MCFace.h>             
#include <CAD4FE_MCEdge.h>              
#include <CAD4FE_MCVertex.h>       
#include <CAD4FE_PolySurface.h>   
#include <CAD4FE_PolyCurve.h>
 
 
#include "CAD4FE_Geometric_Tools.h"
 
 
#pragma package(smart_init)
 
using namespace CAD4FE;
 
bool CAD4FE::GeometricTools::IsInteriorDir(MG_FACE * __face, MG_ARETE * __edge, double __t, double __direction[3])
{
    double tEdge = __t;
    OT_VECTEUR_3D coedgeDerT, faceEdgeInteriorDir, vecTpCoedgeDerT, vecTpFaceEdgeInteriorDir;
    OT_VECTEUR_3D direction(__direction);
 
    double uv[2], xyz[3];
    __edge->evaluer(tEdge,xyz);
    __face->inverser(uv, xyz);
    OT_MATRICE_3D tangentPlaneFrame = TangentPlaneFrame( __face, uv );
    OT_MATRICE_3D tangentPlaneTranform = tangentPlaneFrame.inverse();
    OT_VECTEUR_3D vecTpDirection = tangentPlaneTranform*direction;
    vecTpDirection[2]=0; // project the direction in the tangent plane
    vecTpDirection.norme();
 
    MG_COARETE * coedge1=0;
    double coedge1T;
    int nb_coedge=__edge->get_nb_mg_coarete();
    for (int it_coe=0;it_coe<nb_coedge;it_coe++)
    {
        MG_COARETE *coedge =  __edge->get_mg_coarete(it_coe);
        MG_FACE * adjacent_face = coedge->get_boucle()->get_mg_face();
 
        if (adjacent_face == __face)
        {
              coedge1 = coedge; coedge1T = __t;
              break;
        }
    }
 
    double angle_coedge1_testdirection, cs_angle_coedge1_testdirection, sn_angle_coedge1_testdirection;
    OT_VECTEUR_3D coedge1DerT, coedge2DerT, vecTpCoedge1DerT, vecTpCoedge2DerT;
    coedge1->get_arete()->deriver(coedge1T, coedge1DerT);
    coedge1DerT *= coedge1->get_orientation();           
    vecTpCoedge1DerT = tangentPlaneTranform*coedge1DerT;
    vecTpCoedge1DerT.norme();
 
    cs_angle_coedge1_testdirection = vecTpCoedge1DerT*vecTpDirection;
    sn_angle_coedge1_testdirection = -vecTpCoedge1DerT[1]*vecTpDirection[0]+vecTpCoedge1DerT[0]*vecTpDirection[1];
    if (cs_angle_coedge1_testdirection>1) cs_angle_coedge1_testdirection=1;
    if (cs_angle_coedge1_testdirection < -1) cs_angle_coedge1_testdirection=-1;
    angle_coedge1_testdirection = acos(cs_angle_coedge1_testdirection);
    if (sn_angle_coedge1_testdirection < 0 ) angle_coedge1_testdirection *= -1;
 
    if ( angle_coedge1_testdirection < 0 )
        return false;
    else
        return true;
}
 
bool CAD4FE::GeometricTools::MCFace_MCTopoDir_IsInterior(MCFace * __mcFace, MG_ELEMENT_TOPOLOGIQUE * __mcTopo, double __xyz[3], double __direction[3])
{
    bool result = false;
    switch (__mcTopo->get_dimension())
    {
        case 1:
        {
            MCEdge * mcEdge = (MCEdge *) __mcTopo;
            double s;
            mcEdge->inverser(s,__xyz);
            result = MCFace_MCEdgeDir_IsInterior(__mcFace, mcEdge, s, __direction);
            break;
        }
        case 0:
        {
            MCVertex * mcVertex = (MCVertex*) __mcTopo;
            result = MCFace_MCVertexDir_IsInterior(__mcFace, mcVertex, __direction);
            break;
        }
        default:
        {
            result = false;
            break;
        }
    }
    return result;
}
 
bool CAD4FE::GeometricTools::MCFace_MCEdgeDir_IsInterior(MCFace * __mcFace, MCEdge * __mcEdge, double __s, double __direction[3])
{
    PolySurface * polysurface = __mcFace->GetPolySurface();
    PolyCurve * polycurve = __mcEdge->GetPolyCurve();
    MG_ARETE * edge;
    double tEdge, dtEdge;
    polycurve->Parameter_SToRefEdgeT(__s, &edge, &tEdge, &dtEdge, false);
 
    int nb_coedge=edge->get_nb_mg_coarete();
    for (int it_coe=0;it_coe<nb_coedge;it_coe++)
    {
        MG_FACE * adjacentFace = edge->get_mg_coarete(it_coe)->get_boucle()->get_mg_face();
        if ( polysurface->Contains(adjacentFace) )
        {
            if (IsInteriorDir(adjacentFace, edge, tEdge, __direction))
                return true;
        }
    }
 
    return false;
}
 
bool CAD4FE::GeometricTools::MCFace_MCVertexDir_IsInterior(MCFace * __mcFace, MCVertex * __mcVertex, double __direction[3])
{
    MG_SOMMET * vertex = __mcVertex->GetRefVertex();
    PolySurface * polysurface = __mcFace->GetPolySurface();
 
    std::set<MG_FACE*> adjacentFaces; MG_SOMMET_GetAdjacent_MG_FACE(vertex, adjacentFaces);
    for (std::set<MG_FACE*>::iterator itF = adjacentFaces.begin();
        itF != adjacentFaces.end(); itF++)
    {
        MG_FACE * adjacentFace = *itF;
        if ( polysurface->Contains(adjacentFace) )
        {
            if (IsInteriorDir(adjacentFace, vertex, __direction))
                return true;
        }
    }
 
    return false;
}
 
    
bool CAD4FE::GeometricTools::IsInteriorDir(MG_FACE * __face, MG_SOMMET * __vertex, double __direction[3])
{
    OT_VECTEUR_3D coedgeDerT, faceEdgeInteriorDir, vecTpCoedgeDerT, vecTpFaceEdgeInteriorDir;
    OT_VECTEUR_3D direction(__direction);
 
    double xyz[3], uv[2];
    __vertex->get_point()->evaluer(xyz);
    __face->inverser(uv, xyz);
    OT_MATRICE_3D tangentPlaneFrame = TangentPlaneFrame( __face, uv );
    OT_MATRICE_3D tangentPlaneTranform = tangentPlaneFrame.inverse();
    OT_VECTEUR_3D vecTpDirection = tangentPlaneTranform*direction;
    vecTpDirection[2]=0; // project the direction in the tangent plane
    vecTpDirection.norme();
 
    MG_COARETE * coedge1=0, * coedge2=0;
    double coedge1T, coedge2T;
    for (int it_cov = 0; it_cov < __vertex->get_nb_mg_cosommet(); it_cov++)
    {
        MG_COSOMMET * covertex = __vertex->get_mg_cosommet(it_cov);
        MG_ARETE * edge = covertex->get_arete();
        unsigned nb_adjacent_face = edge->get_nb_mg_coarete();
        for (unsigned i=0; i<nb_adjacent_face; i++)
        {
            MG_COARETE *coedge = edge->get_mg_coarete(i);
            MG_FACE * adjacent_face = coedge->get_boucle()->get_mg_face();
 
            if (adjacent_face == __face)
            {
                    MG_SOMMET *v1 = edge->get_cosommet1()->get_sommet();
                    MG_SOMMET *v2 = edge->get_cosommet2()->get_sommet();  
                    if (coedge->get_orientation() == -1 )
                        std::swap(v1,v2);
                    {                                                                
                        if (__vertex==v1) {coedge2 = coedge; coedge2->get_arete()->inverser(coedge2T,xyz); }
                        if (__vertex==v2) {coedge1 = coedge; coedge1->get_arete()->inverser(coedge1T,xyz); }
                    }
            }
        }
    }
 
    double angle_coedge1_coedge2, cs_angle_coedge1_coedge2, sn_angle_coedge1_coedge2;
    double angle_coedge1_testdirection, cs_angle_coedge1_testdirection, sn_angle_coedge1_testdirection;
    OT_VECTEUR_3D coedge1DerT, coedge2DerT, vecTpCoedge1DerT, vecTpCoedge2DerT;
    coedge1->get_arete()->deriver(coedge1T, coedge1DerT);
    coedge1DerT *= coedge1->get_orientation();           
    vecTpCoedge1DerT = tangentPlaneTranform*coedge1DerT;
    vecTpCoedge1DerT.norme();
    coedge2->get_arete()->deriver(coedge2T, coedge2DerT);
    coedge2DerT *= coedge2->get_orientation();
    vecTpCoedge2DerT = tangentPlaneTranform*coedge2DerT;
    vecTpCoedge2DerT.norme();
 
    cs_angle_coedge1_coedge2 = vecTpCoedge1DerT*vecTpCoedge2DerT;
    if (cs_angle_coedge1_coedge2>1) cs_angle_coedge1_coedge2=1;
    if (cs_angle_coedge1_coedge2<-1) cs_angle_coedge1_coedge2=-1;
    sn_angle_coedge1_coedge2 = -vecTpCoedge1DerT[1]*vecTpCoedge2DerT[0]+vecTpCoedge1DerT[0]*vecTpCoedge2DerT[1];
    angle_coedge1_coedge2 = acos(cs_angle_coedge1_coedge2);
    if (sn_angle_coedge1_coedge2<0) angle_coedge1_coedge2 *= -1;
    cs_angle_coedge1_testdirection = vecTpCoedge1DerT*vecTpDirection;
    sn_angle_coedge1_testdirection = -vecTpCoedge1DerT[1]*vecTpDirection[0]+vecTpCoedge1DerT[0]*vecTpDirection[1];
    if (cs_angle_coedge1_testdirection>1) cs_angle_coedge1_testdirection=1;
    if (cs_angle_coedge1_testdirection<-1) cs_angle_coedge1_testdirection=-1;
    angle_coedge1_testdirection = acos(cs_angle_coedge1_testdirection);
    if (sn_angle_coedge1_testdirection< 0) angle_coedge1_testdirection *= -1;
 
    if ( angle_coedge1_testdirection > M_PI || angle_coedge1_testdirection < angle_coedge1_coedge2)
        return false;
    else
        return true;
}
 
void GeometricTools::Surface_MoveParamInDomain(MG_FACE * __face, double * __uv, double * __xyz)
{
bool bEval = false;
double uvmax[2];
double uvmin[2];
uvmax[0]=__face->get_surface()->get_umax();
uvmin[0]=__face->get_surface()->get_umin();
uvmax[1]=__face->get_surface()->get_vmax();
uvmin[1]=__face->get_surface()->get_vmin();
double uvperiod[2];
uvperiod[0]=__face->get_surface()->get_periode_u();
uvperiod[1]=__face->get_surface()->get_periode_v();
for (int i=0; i<2; i++)
{
    if (uvperiod[i] != 0.0)
    {
        if (__uv[i] > uvmax[i])
            __uv[i] -= uvmax[i]-uvmin[i];
        else
            if (__uv[i] < uvmin[i])
                __uv[i] += uvmax[i]-uvmin[i];
    }
    else
    {
        if (__uv[i] > uvmax[i])
           { __uv[i] = uvmax[i]; bEval = true; }
        else
            if (__uv[i] < uvmin[i])
              {  __uv[i] = uvmin[i]; bEval = true; }
    }
}
if (bEval)
    __face->evaluer(__uv,__xyz);
}
 
/*MCNode * GeometricTools::FindNextPointInDirection(MG_FACE * __face, std::set<MCNode*> & __candidateEndPt, double __uv[2], double __tpDir[3], double __minDist, MCNode * __endPt )
{
    OT_DECALAGE_PARAMETRE ot_decalage(__face->get_surface()->get_periode_u(), __face->get_surface()->get_periode_v());
    std::set <MCNode *>::iterator it_LCI;
    MCNode * xn = 0;
 
    double decal_u=ot_decalage.calcul_decalage_parametre_u(__uv[0]);
    double decal_v=ot_decalage.calcul_decalage_parametre_v(__uv[1]);
    OT_VECTEUR_3D decal_uv, decal_uvc;
    decal_uv[0]=ot_decalage.decalage_parametre_u(__uv[0],decal_u);
    decal_uv[1]=ot_decalage.decalage_parametre_v(__uv[1],decal_v);
 
    double Du = __tpDir[0];
    double Dv = __tpDir[1];
    OT_VECTEUR_3D tpDir(__tpDir[0],__tpDir[1],0);
 
    double xyz[3];
    __face->evaluer(__uv,xyz);
 
    double minT = 1E99;
    for (it_LCI = __candidateEndPt.begin();
        it_LCI != __candidateEndPt.end();
        it_LCI++)
    {
        // projection of it_LCI on line passing by x[i] -- x[i]+Dt2d
        MCNode * xc = *it_LCI;
        double * xc_uv = xc->UV(__face);  
 
        decal_uvc[0]=ot_decalage.decalage_parametre_u(xc_uv[0],decal_u);
        decal_uvc[1]=ot_decalage.decalage_parametre_v(xc_uv[1],decal_v);
 
        double c1 = (decal_uvc[0]-decal_uv[0])*Du + (decal_uvc[1]-decal_uv[1])*Dv;
        double c2 = Du*Du + Dv*Dv;
        double xct = c1/c2;
 
                                         
        if ( xct > 0 && fabs(xct) < fabs(minT) )
        {                  
          double decal_proj[2];
          for (int i=0; i<2; i++) decal_proj[i]=decal_uv[i]+(__tpDir[i])*xct;
          double xc_proj_uv[2], xc_proj_xyz[3];
          xc_proj_uv[0] = ot_decalage.decalage_parametre_u(decal_proj[0],-decal_u);
          xc_proj_uv[1] = ot_decalage.decalage_parametre_v(decal_proj[1],-decal_v);
          Surface_MoveParamInDomain(__face,xc_proj_uv);
          __face->evaluer(xc_proj_uv, xc_proj_xyz);
          double distXcProjXc = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xc->get_coord(), xc_proj_xyz);
          double distX0ProjXc = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xc_proj_xyz, xyz );
          double distX0Xc = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xc->get_coord(), xyz);
          if ( distXcProjXc < __minDist )
          {
            if ( xc->GetRefEdgeMapping().size() > 0)
            {
                OT_VECTEUR_3D xc_vecWDir;
                OT_MATRICE_3D xc_matTpFrame = TangentPlaneFrame(__face, xc_uv);
                xc_vecWDir = xc_matTpFrame*tpDir;
 
                if ( !IsInteriorDir ( __face, xc, xc_vecWDir ) || xc == __endPt )
                {
                        minT = xct;
                        xn = xc;
                }
            }
            else
            {
                minT = xct;
                xn = xc;
            }
          }
        }
    }
 
    return xn;
}
*/
 
OT_MATRICE_3D GeometricTools::TangentPlaneFrame(MG_FACE * __sigma, double uv[2])
{
    OT_MATRICE_3D M;
    OT_VECTEUR_3D DerU, DerV, Normal;
    __sigma->deriver(uv,DerU,DerV);
    Normal = DerU & DerV;
 
    M.change_vecteur1(DerU);
    M.change_vecteur2(DerV);
    M.change_vecteur3(Normal);
 
    return M;
}
 
OT_MATRICE_3D GeometricTools::TangentPlaneFrame_SecondDersTransform(MG_FACE * __sigma, double uv[2], OT_MATRICE_3D & M, OT_MATRICE_3D & C, OT_VECTEUR_3D & xyz )
{
    OT_VECTEUR_3D DerU, DerV, Normal;
    OT_VECTEUR_3D DerUU, DerVV, DerUV;
    __sigma->deriver_seconde(uv, DerUU, DerUV, DerVV, xyz, DerU, DerV);
    Normal = DerU & DerV;
 
    M.change_vecteur1(DerU);
    M.change_vecteur2(DerV);
    M.change_vecteur3(Normal);
 
    C.change_vecteur1(.5*DerUU);
    C.change_vecteur2(.5*DerVV);
    C.change_vecteur3(DerUV);
 
    return M;
}
 
bool GeometricTools::GetIntersectingPlaneDirection(MG_FACE * __face, double __uv[2], double __planeNormal[3], double __directionUV[2])
{
 
    OT_MATRICE_3D matTpFrame, tangentPlaneTransform;
    OT_VECTEUR_3D vecWDir, vecTpDir, vecWOrigin, planeRootPoint;
 
    __face->evaluer(__uv, planeRootPoint);
 
    matTpFrame = TangentPlaneFrame(__face, __uv);
    tangentPlaneTransform = matTpFrame.inverse();
    Intr3D_Surface_Plane_Get_LocalTangentPlane_Direction(__face, __planeNormal, planeRootPoint, __uv, matTpFrame, vecWDir);
    vecTpDir=tangentPlaneTransform*vecWDir;
 
    for (int i=0; i<2; i++)
        __directionUV[i] = vecTpDir[i];
 
    return true;
}
 
int GeometricTools::Intr3D_Surface_Plane_Get_LocalTangentPlane_Direction(
MG_FACE * __sigma,
double __N1[3], double __P1[3], double __uv[2], OT_MATRICE_3D & __tangentPlaneFrame, OT_VECTEUR_3D & __intrLineDirWorld)
{
    // tangent plane origin = surface point at UV
    double afTpOrigin[3];
    __sigma->evaluer(__uv, afTpOrigin);
 
    // Compute the plane #1 / tangent plane
    // intersection line
    OT_VECTEUR_3D afTpNormal = __tangentPlaneFrame.get_vecteur3();
 
    int intrLineResult = OT_ALGORITHME_GEOMETRIQUE::Intr3D_Plane_Plane(afTpNormal, afTpOrigin, __N1, __P1, __intrLineDirWorld);
 
    return  intrLineResult;
}
 
OT_MATRICE_3D GeometricTools::GetPlaneFrame(OT_VECTEUR_3D & __planeNormal)
{
OT_VECTEUR_3D dir1, dir2;
OT_VECTEUR_3D normal(__planeNormal);
normal.norme();
if (!(OPERATEUR::egal(normal.get_x(),0.,0.000001)))
                {
                dir1[0]=normal.get_y();
                dir1[1]=(-normal.get_x());
                dir1[2]=0.;
                }
        else if (!(OPERATEUR::egal(normal.get_y(),0.,0.000001)))
                {
                dir1[0]=0.;
                dir1[1]=(-normal.get_z());
                dir1[2]=normal.get_y();
                }
        else
                {
                dir1[0]=normal.get_z();
                dir1[1]=0.;
                dir1[2]=0.;
                }
dir1.norme();
dir2=normal&dir1;
return OT_MATRICE_3D(dir1, dir2, normal);
}
 
double GeometricTools::AngleInPlane(double *__axisX, double *__direction)
{
    double angle;
    double directionNorm = sqrt(__direction[0]*__direction[0]+__direction[1]*__direction[1]);
    double axisXNorm = sqrt(__axisX[0]*__axisX[0]+__axisX[1]*__axisX[1]);
    //double axisYNorm = sqrt(__axisY[0]*__axisY[0]+__axisY[1]*__axisY[1]);
    double cs = (+__direction[0]*__axisX[0] + __direction[1]*__axisX[1]) / (axisXNorm*directionNorm);
    double sn = (-__direction[0]*__axisX[1] + __direction[1]*__axisX[0]);// / (axisYNorm*directionNorm);
    if (cs > 1) cs=1;
    if (cs < -1) cs=-1;
    if (sn < 0) return -acos(cs);
    else return acos(cs);
}
 
 
bool GeometricTools::MG_FACE_Contains_MG_SOMMET(MG_FACE * __face, MG_SOMMET * v)
{
    for (int itE = 0; itE < v->get_nb_mg_cosommet(); itE++)
    {
        MG_ARETE * edge = v->get_mg_cosommet(itE)->get_arete();
        for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
        {
            MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
            if (face == __face)
                return true;
        }
    }
    return false;
}
 
bool GeometricTools::MG_FACE_Contains_MG_ARETE(MG_FACE * __face, MG_ARETE * e)
{
    int nb_coedge=e->get_nb_mg_coarete();
    for (int it_coe=0;it_coe<nb_coedge;it_coe++)
        if ( e->get_mg_coarete(it_coe)->get_boucle()->get_mg_face() == __face )
            return true;
 
    return false;
}
 
bool GeometricTools::MG_ARETE_Contains_MG_SOMMET(MG_ARETE * __edge, MG_SOMMET * __vertex)
{
    MG_SOMMET * v1 = __edge->get_cosommet1()->get_sommet();
    if (v1 == __vertex) return true;
    MG_SOMMET * v2 = __edge->get_cosommet2()->get_sommet();    
    if (v2 == __vertex) return true;
    return false;
}
 
bool GeometricTools::PolySurface_Contains_RefVertex(PolySurface * __polysurface, MG_SOMMET * __refVertex)
{
    int nb_covertex = __refVertex->get_nb_mg_cosommet();
    for (int it_cov=0;it_cov<nb_covertex;it_cov++)
    {
        MG_COSOMMET * cov = __refVertex->get_mg_cosommet(it_cov);
        MG_ARETE * e = cov->get_arete();
        int nb_coedge = e->get_nb_mg_coarete();
        for (int it_coe=0;it_coe<nb_coedge;it_coe++)
        {
            MG_FACE * refFace = e->get_mg_coarete(it_coe)->get_boucle()->get_mg_face();
            if (__polysurface->Contains(refFace))
                return true;
        }
    }
    return false;
}
 
void GeometricTools::MG_SOMMET_GetAdjacent_MG_FACE(MG_SOMMET * __refVertex, std::set<MG_FACE*> & __adjacentfaces)
{
    int nb_covertex = __refVertex->get_nb_mg_cosommet();
    for (int it_cov=0;it_cov<nb_covertex;it_cov++)
    {
        MG_COSOMMET * cov = __refVertex->get_mg_cosommet(it_cov);
        MG_ARETE * e = cov->get_arete();
        int nb_coedge = e->get_nb_mg_coarete();
        for (int it_coe=0;it_coe<nb_coedge;it_coe++)
        {
            MG_FACE * refFace = e->get_mg_coarete(it_coe)->get_boucle()->get_mg_face();
            __adjacentfaces.insert(refFace);
        }
    }
}
 
double GeometricTools::FacePointCorrection(MG_FACE * __face, double __xyz[3], double __uv[2])
{
        OT_MATRICE_3D matTpFrame, tpTrans;
    OT_VECTEUR_3D d_uvn(1,1,1);
        int i=0;
        double error;
        OT_VECTEUR_3D xyz(__xyz);
        OT_VECTEUR_3D xyzEval,xyzEval2,dx2;
    __face->valide_parametre_u(__uv[0]);
    __face->valide_parametre_v(__uv[1]);
        __face->evaluer(__uv, xyzEval);
        OT_VECTEUR_3D dx=xyz-xyzEval;
 
 
        for (i=0;
                i<10 &&
                dx.get_longueur2() > 1E-28 &&
                (d_uvn[0]*d_uvn[0])+(d_uvn[1]*d_uvn[1]) > 1E-8*d_uvn[2]*d_uvn[2] // dx perpendicular to tangent plane ?
                ;
            i++)
        {
                matTpFrame=TangentPlaneFrame(__face,__uv);
                tpTrans=matTpFrame.inverse();
                d_uvn=tpTrans*dx;
                __uv[0]+=d_uvn[0];
                __uv[1]+=d_uvn[1];
                __face->valide_parametre_u(__uv[0]);
                __face->valide_parametre_v(__uv[1]);
                __face->evaluer(__uv, xyzEval);
                dx=xyz-xyzEval;
        }            
    error=dx.get_longueur();
 
    return error;
}
 
double GeometricTools::Vector_UV_To_3D_Plane(MG_FACE *__face, double *__planeNormal, double * __planeRootPoint, double * __uv1, double * __uv12)
{
        OT_VECTEUR_3D P1, P2, vec_UV1_UV2, vec_P1_P2;
 
    // Compute the distance vector between a P(uv) and the plane
    __face->evaluer(__uv1, P1);
    OT_VECTEUR_3D vecP1P2;
    double distance3D = OT_ALGORITHME_GEOMETRIQUE::Closest_Point_to_Plane_3d(__planeNormal, __planeRootPoint, P1, P2);
    vec_P1_P2 = P2 - P1;
 
    // Compute tangent plane tranformation matrix
    OT_MATRICE_3D tangentPlaneFrame3D = GeometricTools::TangentPlaneFrame(__face, __uv1);
    OT_MATRICE_3D tangentPlaneTransform = tangentPlaneFrame3D.inverse();
 
    // compute the value of the vector in the tangent plane frame of the face
    vec_UV1_UV2 = tangentPlaneTransform * vec_P1_P2;
 
    for (int i=0;i<2;i++)
        __uv12[i] = vec_UV1_UV2[i];
 
    return distance3D;
}
 
 
double GeometricTools::MoveParamInPlane(MG_FACE * __face, double __planeNormal[3], double __planeRootPoint[3], double deltamax3D, int nbPointToPlaneResidueCorrection, double __uv[2], double __xyz[3])
{
    OT_MATRICE_3D tangentPlaneFrame3D, tangentPlaneTransform;
    OT_VECTEUR_3D uv1(__uv);
    OT_VECTEUR_3D uv2;
        OT_VECTEUR_3D P1(__xyz);
    OT_VECTEUR_3D P2, vec_UV1_UV2, vec_P1_P2;
    OT_VECTEUR_3D vecP1P2;
    int nbPlaneCorrections=0;
    double error;
 
    error=OT_ALGORITHME_GEOMETRIQUE::Closest_Point_to_Plane_3d(__planeNormal, __planeRootPoint, P1, P2);
    for(nbPlaneCorrections = 0;
        nbPlaneCorrections < nbPointToPlaneResidueCorrection && deltamax3D<error;
        nbPlaneCorrections ++)
    {
            tangentPlaneFrame3D = GeometricTools::TangentPlaneFrame(__face, uv1);
            tangentPlaneTransform = tangentPlaneFrame3D.inverse();
            vec_P1_P2 = P2 - P1;
            vec_UV1_UV2 = tangentPlaneTransform*vec_P1_P2;
            uv2 = uv1 + vec_UV1_UV2;
            uv1=uv2;
            Surface_MoveParamInDomain(__face, uv1, P1); 
            __face->evaluer(uv1,P1);
            error=OT_ALGORITHME_GEOMETRIQUE::Closest_Point_to_Plane_3d(__planeNormal, __planeRootPoint, P1, P2);
    }
 
    for (int i=0; i<2; i++)
        __uv[i] = uv1[i];
 
    for (int i=0; i<3; i++)
        __xyz[i] = P1[i];
 
    return error;
}
 
double GeometricTools::Segment2dCurvilinearLength(MG_FACE * __face, double __uv1[2], double __uv2[2], unsigned __nbNubdivisions)
{
    double E,F,G;                 
    double decal_uv2[2],decal_uv1[2],duv[2],uv[2],decal[2]={0,0};
    double period[2];period[0]=__face->get_surface()->get_periode_u();period[1]=__face->get_surface()->get_periode_v();
    OT_DECALAGE_PARAMETRE d(period[0],period[1]);
 
    decal[0]=d.calcul_decalage_parametre_u(__uv1[0]);
    decal[1]=d.calcul_decalage_parametre_v(__uv1[1]);
    decal_uv1[0]=d.decalage_parametre_u(__uv1[0],decal[0]);
    decal_uv1[1]=d.decalage_parametre_v(__uv1[1],decal[1]);
    decal_uv2[0]=d.decalage_parametre_u(__uv2[0],decal[0]);
    decal_uv2[1]=d.decalage_parametre_v(__uv2[1],decal[1]);
 
    duv[0]=1.0/__nbNubdivisions*(decal_uv2[0]-decal_uv1[0]);
    duv[1]=1.0/__nbNubdivisions*(decal_uv2[1]-decal_uv1[1]);
 
    double length=0;
    for (unsigned i=0; i<=__nbNubdivisions; i++)
    {
        uv[0]=d.decalage_parametre_u(decal_uv1[0]+i*duv[0],-decal[0]);
        uv[1]=d.decalage_parametre_v(decal_uv1[1]+i*duv[1],-decal[1]);
        __face->valide_parametre_u(uv[0]);
        __face->valide_parametre_v(uv[1]);
 
        __face->get_EFG(uv, E, F, G);
        if (i > 0)
        {
            double dist2=duv[0]*duv[0]*E+duv[1]*duv[1]*G+duv[1]*duv[0]*F;
            double dist=0;
            if ( dist2 > 1E-308)
                dist=sqrt(dist2);
            length += dist;
        }
    }
    
    return length;
}
                 
 
double GeometricTools::ComputeDiedralAngle(MG_NOEUD * __tri[2][3])
{
    OT_VECTEUR_3D normal[2];
 
    for (int i=0;i<2;i++)
    {
        OT_VECTEUR_3D s[2];
        
        for (int j=0;j<2;j++)
        for (int k=0;k<3;k++)
            s[j][k]=__tri[i][j+1]->get_coord()[k]-__tri[i][j]->get_coord()[k];
 
        normal[i]=s[1]&s[0];
 
        if (normal[i].get_longueur2() == 0) // degenerated triangle
            return 1E308;
 
        normal[i].norme();
    }
 
    double cs = normal[0]*normal[1];
    if (cs>1)cs=1;
    if (cs<-1)cs=-1;
    double angle = acos(cs);
 
    return angle;
}
 
int GeometricTools::MG_TOPO_GetColor(MG_ELEMENT_TOPOLOGIQUE* ele, double & couleur)
{
couleur=0;
int nbccf=ele->get_nb_ccf();
for (int k=0;k<nbccf;k++)
    {
        char typeccf[10];
        ele->get_type_ccf(k,typeccf);
        if ((typeccf[0]=='C') && (typeccf[1]=='c'))
            {
                 couleur = ele->get_valeur_ccf(k);
                 return 1;
            }
    }
return 0;
}
 
int GeometricTools::MG_TOPO_GetColor(MG_ELEMENT_TOPOLOGIQUE* ele, long & rgba)
{
    double val;
    int result = MG_TOPO_GetColor(ele, val);
 
    // if there is no color information
    if (result == 0)
    {
        rgba = -1;
        return result;
    }
 
    // Codage couleur : 4 octets dans un double (8 octets)
    // -0-1-2-3-4-5-6-7-
    // |R|G|B|A|0|0|0|0|
    rgba = 0;
    unsigned char *p=(unsigned char*)&val;
    rgba += (*p);p++;
    rgba += ((*p)<<8);p++;
    rgba += ((*p)<<16);p++;
    rgba += ((*p)<<24);
    return result;
}
 
int GeometricTools::MG_TOPO_GetColor(MG_ELEMENT_TOPOLOGIQUE* ele, int & rgba)
{
    double val;
    int result = MG_TOPO_GetColor(ele, val);
 
    // Codage couleur : 4 octets dans un double (8 octets)
    // -0-1-2-3-4-5-6-7-
    // |R|G|B|A|0|0|0|0|
    unsigned char *p=(unsigned char*)&val;
    rgba = 0;
    rgba += (*p);p++;
    rgba += ((*p)<<8);p++;
    rgba += ((*p)<<16);p++;
    rgba += ((*p)<<24);   
    return result;
}
 
int GeometricTools::MG_TOPO_GetColor(MG_ELEMENT_TOPOLOGIQUE* ele, unsigned char rgba[4])
{
    double val;
    int result = MG_TOPO_GetColor(ele, val);
    ccf_double_to_uc4(val, rgba);   
    return result;
}
 
void GeometricTools::ccf_uc4_to_double(unsigned char rgba[4], double & val)
{
    val = 0;
    unsigned char *p=(unsigned char*)&val;
    for (int i=0; i<4; i++)
        p[i] = rgba[i];
}
 
 
void GeometricTools::ccf_double_to_uc4(double val, unsigned char rgba[4])
{
    // Codage couleur : 4 octets dans un double (8 octets)
    // -0-1-2-3-4-5-6-7-
    // |R|G|B|A|0|0|0|0|
    unsigned char *p=(unsigned char*)&val;  
    for (int i=0; i<4; i++)
        rgba[i] = p[i];
}
 
 
 
void GeometricTools::MG_TOPO_SetColor(MG_ELEMENT_TOPOLOGIQUE* ele, unsigned char rgba[4])
{
    double couleur = 0;
    ccf_uc4_to_double(rgba, couleur);
    MG_TOPO_SetColor(ele, couleur);
}
                 
void GeometricTools::MG_TOPO_SetColor(MG_ELEMENT_TOPOLOGIQUE* ele, double d)
{
    ele->ajouter_ccf((char*)"Cc", d);
}
 
void GeometricTools::MG_TOPO_SetColor(MG_ELEMENT_TOPOLOGIQUE* ele, unsigned char r, unsigned char g, unsigned char b, unsigned char a)
{
    unsigned char rgba[4] = {r,g,b,a};
    MG_TOPO_SetColor(ele, rgba);
}