CAD4FE/src/CAD4FE_OptimizeEdgeSwap.cpp

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  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_OptimizeEdgeSwap.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:56 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
#include <set>

#pragma hdrstop
                    
#include "gestionversion.h"
#include <math.h>
#include "CAD4FE_OptimizeEdgeSwap.h"
#include "ot_mathematique.h"
#include "ot_boite_2d.h"
#include "tpl_map_entite.h"
#include "CAD4FE_m3d_MCTriangle.h"
#include "mg_maillage.h"   

#include "CAD4FE_Intersection_MCSegment_MCSegment.h"
#include "CAD4FE_MCSegment_Middle.h"                
#include "CAD4FE_MCSegment.h"                    
#include "CAD4FE_MCNode.h"
#include "ot_algorithme_geometrique.h"   
#include "CAD4FE_Geometric_Tools.h"  


#pragma package(smart_init)

using namespace CAD4FE;

OptimizeEdgeSwap::OptimizeEdgeSwap(MG_MAILLAGE* __mesh)
{
    _mesh=__mesh;
    maxDiedralAngle = 45*M_PI/180;
}

MG_SEGMENT* OptimizeEdgeSwap::GetSegmentAfter(MG_TRIANGLE *__triangle, MG_SEGMENT *__segment)
{
    if (__segment == __triangle->get_segment1())
        return __triangle->get_segment2();
    if (__segment == __triangle->get_segment2())
        return __triangle->get_segment3();
    if (__segment == __triangle->get_segment3())
        return __triangle->get_segment1();
    return NULL;
}

MG_SEGMENT* OptimizeEdgeSwap::GetSegmentBefore(MG_TRIANGLE *__triangle, MG_SEGMENT *__segment)
{                                    
    if (__segment == __triangle->get_segment1())
        return __triangle->get_segment3();
    if (__segment == __triangle->get_segment2())
        return __triangle->get_segment1();
    if (__segment == __triangle->get_segment3())
        return __triangle->get_segment2();
    return NULL;
}
                    
int OptimizeEdgeSwap::CheckOverlapTriangle(MG_NOEUD * __tri[2][3])
{
    for (int i=0;i<2;i++)
    for (int j=0;j<3;j++)
    {
        MG_NOEUD ** t = __tri[(i+1)%2];
        MG_NOEUD * v = __tri[i][j];
        
        if (t[0]==v||t[1]==v||t[2]==v) continue;

        double I[3],IT;
        int res = OT_ALGORITHME_GEOMETRIQUE::project_PointTriangle(v->get_coord(),t[0]->get_coord(),t[1]->get_coord(),t[2]->get_coord(),I,&IT);
        if (res != 0)
            return res;
    }
    return 0;
}

double OptimizeEdgeSwap::SwapScore(MG_SEGMENT *__segment)
{
    return Swap(__segment,false);
}

double OptimizeEdgeSwap::SwapSegment(MG_SEGMENT *__segment)
{
    MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
    MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
    no[2] = GetStartNode(triangle[0],__segment);
    no[1] = GetNodeBefore(triangle[0],no[2]);
    no[0] = GetNodeBefore(triangle[0],no[1]);
    no[3] = GetNodeBefore(triangle[1],no[0]);
    MG_SEGMENT * seg[4];
    seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id());  // AC
    seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id());  // CB
    seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id());  // BD
    seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id());  // DA

    double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
    double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
    MCSegment * swapSegment = new MCSegment(triangle[0]->get_lien_topologie(), (MCNode*)no[1], (MCNode*)no[3] );
    if (swapSegment->GetRefFaceMapping().size() == 0)
    {
        delete swapSegment;
        return 0;
    }
    _mesh->ajouter_mg_segment(swapSegment);
    M3D_MCTriangle * newTriangle[2];
    newTriangle[0] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode*)no[1], (MCNode*)no[2], (MCNode*)no[3], (MCSegment*)seg[1], (MCSegment*)seg[2], swapSegment);
    newTriangle[0]->change_qualite(crit21);
    newTriangle[1] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode*)no[3], (MCNode*)no[0], (MCNode*)no[1], (MCSegment*)seg[3], (MCSegment*)seg[0], swapSegment);
    newTriangle[1]->change_qualite(crit22);
    for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
    for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
    _mesh->supprimer_mg_segmentid(__segment->get_id());

    return 1;
}

/*double OptimizeEdgeSwap::Swap(MG_SEGMENT *__segment, bool __execute)
{
    double score1;
    if (__segment->get_lien_triangle()->get_nb() != 2)
    {
        return -1;
    }
    if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
    {
        return -2;
    }
    if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
    {
        return -3;
    }

    MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
    MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
    no[2] = GetStartNode(triangle[0],__segment);
    no[1] = GetNodeBefore(triangle[0],no[2]);
    no[0] = GetNodeBefore(triangle[0],no[1]);
    no[3] = GetNodeBefore(triangle[1],no[0]);
    MG_SEGMENT * seg[4];
    seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id());  // AC
    seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id());  // CB
    seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id());  // BD
    seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id());  // DA

    double crit11,crit12,crit21,crit22;
    crit11 = OPERATEUR::qualite_triangle(no[0]->get_coord(),no[1]->get_coord(),no[2]->get_coord());
    crit12 = OPERATEUR::qualite_triangle(no[2]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
    crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
    crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
                                                
    double crit1=std::min(crit11,crit12);
    double crit2=std::min(crit21,crit22);

    score1 = (crit1 > .5) ? -1 : (crit2-crit1);
    if (score1 < 0) return score1;

    MG_NOEUD* tris2[2][3]={{no[3],no[1],no[2]},{no[1],no[3],no[0]}};
    int testOverlapTriangles = CheckOverlapTriangle(tris2);
    score1 = ( testOverlapTriangles != 0 ) ? -1 : score1;
    if (score1 < 0) return score1;

    double diedralAngle2 = GeometricTools::ComputeDiedralAngle(tris2);
    score1 = ( diedralAngle2 > maxDiedralAngle ) ? -1 : score1; 
    if (score1 < 0) return score1;

    if (__execute && score1 > 0)
    {
        double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
        double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
        MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
        MCSegment * swapSegment = new MCSegment(topo, (MCNode*)no[1], (MCNode*)no[3] );
        _mesh->ajouter_mg_segment(swapSegment);
        M3D_MCTriangle * newTriangle[2];
        newTriangle[0] = new M3D_MCTriangle(topo, (MCNode*)no[1], (MCNode*)no[2], (MCNode*)no[3], (MCSegment*)seg[1], (MCSegment*)seg[2], swapSegment);
        newTriangle[0]->change_qualite(crit21);
        newTriangle[1] = new M3D_MCTriangle(topo, (MCNode*)no[3], (MCNode*)no[0], (MCNode*)no[1], (MCSegment*)seg[3], (MCSegment*)seg[0], swapSegment);
        newTriangle[1]->change_qualite(crit22);
        for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
        for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
        _mesh->supprimer_mg_segmentid(__segment->get_id());
    }

    return score1;
}                 */


double OptimizeEdgeSwap::Swap(MG_SEGMENT *__segment, bool __execute)
{
    double score1;
    if (__segment->get_lien_triangle()->get_nb() != 2)
    {
        return -1;
    }
    if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
    {
        return -2;
    }
    if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
    {
        return -3;
    }

    MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
    MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
    no[2] = GetStartNode(triangle[0],__segment);
    no[1] = GetNodeBefore(triangle[0],no[2]);
    no[0] = GetNodeBefore(triangle[0],no[1]); 
    no[3] = GetNodeBefore(triangle[1],no[0]);
    MG_SEGMENT * seg[4];
    seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id());  // AC
    seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id());  // CB
    seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id());  // BD
    seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id());  // DA

    MG_NOEUD* tris[2][2][3]={{{no[2],no[3],no[0]},{no[0],no[1],no[2]}}, // triangles of the current configuration
    {{no[3],no[1],no[2]},{no[1],no[3],no[0]}}}; // triangles of the candidate configuration (after segment swapping)


    int testOverlapTriangles = CheckOverlapTriangle(tris[1]);
    if (testOverlapTriangles) return -1; 
    
    // projection des triangles sur le plan tangent
    double crit[2][2];
    double area[2][2];
    for (int k=0;k<2;k++)
    for (int i=0;i<2;i++)
    {
        OT_VECTEUR_3D normal[4];
        MCFace * face = (MCFace *)triangle[i]->get_lien_topologie();
        for (int j=0;j<3;j++)
            ((MCNode*)tris[k][i][j])->NormalMCFace(face,normal[j]);
        normal[3] = .3333*(normal[0]+normal[1]+normal[2]);
        // frame (x,y,n)
        OT_MATRICE_3D tangentPlane = GeometricTools::GetPlaneFrame(normal[3]);
        OT_MATRICE_3D trf; tangentPlane.transpose(trf);
        OT_VECTEUR_3D tpProjNodes[3];
        for (int j=0;j<3;j++)
        {
            OT_VECTEUR_3D p(tris[k][i][j]->get_coord());
            tpProjNodes[j] = trf*p;
            // projection --> n component = 0
            tpProjNodes[j][2] = 0;
        }
        crit[k][i]=OPERATEUR::qualite_triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
        area[k][i]=OT_ALGORITHME_GEOMETRIQUE::Area_Triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
    }
    double crit1=std::min(crit[0][0],crit[0][1]);
    double crit2=std::min(crit[1][0],crit[1][1]);
                                                   
    // dist_P[i] = la distance entre la diagonale de la ieme configuration et la MCFace
    // voir page 110 de these_jfl.pdf
    double critGeomAccuracy[2];
    {
        MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
        MCSegment * AB = (MCSegment*)__segment;
        MCNode Pgeo;
        //    MCSegment * CD = new MCSegment(topo, (MCNode*)no[1], (MCNode*)no[3] );
        //    CD->evaluer_geo(((area[0][0]*0+area[0][1]*1)/(area[0][0]+area[0][1])),P[0]);
        {
            double wA, wB;
            wA=(AB->get_noeud1() == no[0])?0:1;
            if(AB->get_orientation() == -1) wA=(wA==0)?1:0;
            wB=(wA==0)?1:0; 
            AB->evaluer_geo(((area[1][0]*wA+area[1][1]*wB)/(area[1][0]+area[1][1])),&Pgeo);
        }

        double P_Euclid[2][3];
        for (int i=0;i<3;i++) {                               
            P_Euclid[1][i]=((area[0][0]*no[1]->get_coord()[i]+area[0][1]*no[3]->get_coord()[i])/(area[0][0]+area[0][1]));
            P_Euclid[0][i]=((area[1][0]*no[0]->get_coord()[i]+area[1][1]*no[2]->get_coord()[i])/(area[1][0]+area[1][1]));
        }                   
        double distP[2],segLength[2];
        for (int i=0;i<2;i++)
            distP[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(P_Euclid[i],Pgeo.get_coord());
        for (int i=0;i<2;i++)
            segLength[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(no[i]->get_coord(),no[i+2]->get_coord());
        double maxDistP=.5*std::max(segLength[0],segLength[1]);
        for (int i=0;i<2;i++)
            if ( distP[i] < maxDistP )
                critGeomAccuracy[i] = 1-distP[i]/maxDistP;
            else
                critGeomAccuracy[i] = 0;
    }

    crit1 += critGeomAccuracy[0];
    crit2 += critGeomAccuracy[1];

    if (__execute && crit2 > crit1 )
    {
        SwapSegment(__segment);
    }

    return (crit2>crit1)?crit2:-1;
}         


MG_NOEUD * OptimizeEdgeSwap::GetStartNode(MG_TRIANGLE *__t, MG_SEGMENT *__s)
{
    MG_NOEUD *TN[3];
    TN[0] = __t->get_noeud1();
    TN[1] = __t->get_noeud2();
    TN[2] = __t->get_noeud3();
    for (int i=0; i<3; i++)
    {
        if (__s->get_noeud1() == TN[i])
        {
            if (__s->get_noeud2() == TN[(i+1)%3])
                return __s->get_noeud1();
            else
                return __s->get_noeud2();
        }
    }
    return NULL;
}

MG_NOEUD * OptimizeEdgeSwap::GetNodeAfter(MG_TRIANGLE *__t, MG_NOEUD *__n)
{
    MG_NOEUD *TN[3];
    TN[0] = __t->get_noeud1();
    TN[1] = __t->get_noeud2();
    TN[2] = __t->get_noeud3();
    for (int i=0; i<3; i++)
    {
        if (TN[i]==__n)
            return TN[(i+1)%3];
    }
    return NULL;
}

MG_NOEUD * OptimizeEdgeSwap::GetNodeBefore(MG_TRIANGLE *__t, MG_NOEUD *__n)
{
    MG_NOEUD *TN[3];
    TN[0] = __t->get_noeud1();
    TN[1] = __t->get_noeud2();
    TN[2] = __t->get_noeud3();
    for (int i=0; i<3; i++)
    {
        if (TN[i]==__n)
            return TN[(i+2)%3];
    }
    return NULL;
}

MG_NOEUD * OptimizeEdgeSwap::GetOppositeNode(MG_TRIANGLE *__t, MG_SEGMENT *__s)
{             
    MG_NOEUD *TN[3];
    TN[0] = __t->get_noeud1();
    TN[1] = __t->get_noeud2();
    TN[2] = __t->get_noeud3();
    for (int i=0; i<3; i++)
        if (TN[i] != __s->get_noeud1() && TN[i] != __s->get_noeud2())
            return TN[i];
    return NULL;
}

void OptimizeEdgeSwap::RemoveVolumeMesh()
{
        LISTE_MG_SEGMENT::iterator itSeg;
        LISTE_MG_NOEUD::iterator itNo;   
        LISTE_MG_TRIANGLE::iterator itTri;
        _mesh->supprimer_tout_mg_tetra();
        for (MG_TRIANGLE * tri = _mesh->get_premier_triangle(itTri);tri; tri = _mesh->get_suivant_triangle(itTri))
        {                 
            if (tri->get_lien_topologie()->get_dimension()==3)
                _mesh->supprimer_mg_triangleid(tri->get_id());
        }     
        for (MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
        {
            if (seg->get_lien_topologie()->get_dimension()==3)
                _mesh->supprimer_mg_segmentid(seg->get_id());
        }
        for (MG_NOEUD * no = _mesh->get_premier_noeud(itNo);no; no = _mesh->get_suivant_noeud(itNo))
        {
            if (no->get_lien_topologie()->get_dimension()==3)
                _mesh->supprimer_mg_noeudid(no->get_id());
        }
}

int OptimizeEdgeSwap::Optimize(MG_FACE * __face)
{                         
        int nbSwap = 0;               
        LISTE_MG_SEGMENT::iterator itSeg;
        MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);   
        std::map<double, MG_SEGMENT*> lst;

        for (int j=0;j<3;j++)
        {
            lst.clear();
            int tmpnbSwap = nbSwap;
            for (seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
            {
                MG_ELEMENT_TOPOLOGIQUE * topo = seg->get_lien_topologie();
                if (__face != NULL && topo != __face)
                    continue;
                if (topo->get_dimension()!=2)
                    continue;
                lst[SwapScore(seg)]=seg;
            }

            for (std::map<double, MG_SEGMENT*>::reverse_iterator i=lst.rbegin();i!=lst.rend();i++)
            {
                double score1 = i->first;
                seg=i->second;

                if (score1<0)
                    break;

                double score2 = SwapSegment(seg);
                if (score2 > 0)
                {
                    nbSwap++;
                }
            }
            if (nbSwap-tmpnbSwap==0)
                break;
        }

        return nbSwap;
}

int OptimizeEdgeSwap::OptimizeAllFaces()
{
    return Optimize(0);
}

Revision:	1158
Committed:	Thu Jun 13 22:18:49 2024 UTC (11 months ago) by francois
File size:	17512 byte(s)
Log Message:	compatibilité Ubuntu 22.04 Suppression des refeences à Windows Ajout d'une banière
#	User	Rev	Content
1	francois	1158	//####//------------------------------------------------------------
2			//####//------------------------------------------------------------
3			//####// MAGiC
4			//####// Jean Christophe Cuilliere et Vincent FRANCOIS
5			//####// Departement de Genie Mecanique - UQTR
6			//####//------------------------------------------------------------
7			//####// MAGIC est un projet de recherche de l equipe ERICCA
8			//####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9			//####// http://www.uqtr.ca/ericca
10			//####// http://www.uqtr.ca/
11			//####//------------------------------------------------------------
12			//####//------------------------------------------------------------
13			//####//
14			//####// CAD4FE_OptimizeEdgeSwap.cpp
15			//####//
16			//####//------------------------------------------------------------
17			//####//------------------------------------------------------------
18			//####// COPYRIGHT 2000-2024
19			//####// jeu 13 jun 2024 11:58:56 EDT
20			//####//------------------------------------------------------------
21			//####//------------------------------------------------------------
22	foucault	27	#include <set>
23
24			#pragma hdrstop
25
26			#include "gestionversion.h"
27			#include <math.h>
28			#include "CAD4FE_OptimizeEdgeSwap.h"
29			#include "ot_mathematique.h"
30	foucault	569	#include "ot_boite_2d.h"
31	foucault	27	#include "tpl_map_entite.h"
32			#include "CAD4FE_m3d_MCTriangle.h"
33			#include "mg_maillage.h"
34
35			#include "CAD4FE_Intersection_MCSegment_MCSegment.h"
36			#include "CAD4FE_MCSegment_Middle.h"
37			#include "CAD4FE_MCSegment.h"
38			#include "CAD4FE_MCNode.h"
39	foucault	569	#include "ot_algorithme_geometrique.h"
40			#include "CAD4FE_Geometric_Tools.h"
41	foucault	27
42
43			#pragma package(smart_init)
44
45			using namespace CAD4FE;
46
47			OptimizeEdgeSwap::OptimizeEdgeSwap(MG_MAILLAGE* __mesh)
48			{
49			_mesh=__mesh;
50			maxDiedralAngle = 45*M_PI/180;
51			}
52
53			MG_SEGMENT* OptimizeEdgeSwap::GetSegmentAfter(MG_TRIANGLE __triangle, MG_SEGMENT __segment)
54			{
55			if (__segment == __triangle->get_segment1())
56			return __triangle->get_segment2();
57			if (__segment == __triangle->get_segment2())
58			return __triangle->get_segment3();
59			if (__segment == __triangle->get_segment3())
60			return __triangle->get_segment1();
61			return NULL;
62			}
63
64			MG_SEGMENT* OptimizeEdgeSwap::GetSegmentBefore(MG_TRIANGLE __triangle, MG_SEGMENT __segment)
65			{
66			if (__segment == __triangle->get_segment1())
67			return __triangle->get_segment3();
68			if (__segment == __triangle->get_segment2())
69			return __triangle->get_segment1();
70			if (__segment == __triangle->get_segment3())
71			return __triangle->get_segment2();
72			return NULL;
73			}
74
75			int OptimizeEdgeSwap::CheckOverlapTriangle(MG_NOEUD * __tri[2][3])
76			{
77			for (int i=0;i<2;i++)
78			for (int j=0;j<3;j++)
79			{
80			MG_NOEUD ** t = __tri[(i+1)%2];
81			MG_NOEUD * v = __tri[i][j];
82
83			if (t[0]==v\|\|t[1]==v\|\|t[2]==v) continue;
84
85			double I[3],IT;
86			int res = OT_ALGORITHME_GEOMETRIQUE::project_PointTriangle(v->get_coord(),t[0]->get_coord(),t[1]->get_coord(),t[2]->get_coord(),I,&IT);
87			if (res != 0)
88			return res;
89			}
90			return 0;
91			}
92
93			double OptimizeEdgeSwap::SwapScore(MG_SEGMENT *__segment)
94			{
95			return Swap(__segment,false);
96			}
97
98			double OptimizeEdgeSwap::SwapSegment(MG_SEGMENT *__segment)
99			{
100			MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
101			MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
102			no[2] = GetStartNode(triangle[0],__segment);
103			no[1] = GetNodeBefore(triangle[0],no[2]);
104			no[0] = GetNodeBefore(triangle[0],no[1]);
105			no[3] = GetNodeBefore(triangle[1],no[0]);
106			MG_SEGMENT * seg[4];
107			seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
108			seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
109			seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
110			seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
111
112			double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
113			double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
114			MCSegment * swapSegment = new MCSegment(triangle[0]->get_lien_topologie(), (MCNode)no[1], (MCNode)no[3] );
115			if (swapSegment->GetRefFaceMapping().size() == 0)
116			{
117			delete swapSegment;
118			return 0;
119			}
120			_mesh->ajouter_mg_segment(swapSegment);
121			M3D_MCTriangle * newTriangle[2];
122			newTriangle[0] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode)no[1], (MCNode)no[2], (MCNode)no[3], (MCSegment)seg[1], (MCSegment*)seg[2], swapSegment);
123			newTriangle[0]->change_qualite(crit21);
124			newTriangle[1] = new M3D_MCTriangle(triangle[0]->get_lien_topologie(), (MCNode)no[3], (MCNode)no[0], (MCNode)no[1], (MCSegment)seg[3], (MCSegment*)seg[0], swapSegment);
125			newTriangle[1]->change_qualite(crit22);
126			for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
127			for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
128			_mesh->supprimer_mg_segmentid(__segment->get_id());
129
130			return 1;
131			}
132
133			/double OptimizeEdgeSwap::Swap(MG_SEGMENT __segment, bool __execute)
134			{
135			double score1;
136			if (__segment->get_lien_triangle()->get_nb() != 2)
137			{
138			return -1;
139			}
140			if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
141			{
142			return -2;
143			}
144			if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
145			{
146			return -3;
147			}
148
149			MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
150			MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
151			no[2] = GetStartNode(triangle[0],__segment);
152			no[1] = GetNodeBefore(triangle[0],no[2]);
153			no[0] = GetNodeBefore(triangle[0],no[1]);
154			no[3] = GetNodeBefore(triangle[1],no[0]);
155			MG_SEGMENT * seg[4];
156			seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
157			seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
158			seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
159			seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
160
161			double crit11,crit12,crit21,crit22;
162			crit11 = OPERATEUR::qualite_triangle(no[0]->get_coord(),no[1]->get_coord(),no[2]->get_coord());
163			crit12 = OPERATEUR::qualite_triangle(no[2]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
164			crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
165			crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
166
167			double crit1=std::min(crit11,crit12);
168			double crit2=std::min(crit21,crit22);
169
170			score1 = (crit1 > .5) ? -1 : (crit2-crit1);
171			if (score1 < 0) return score1;
172
173			MG_NOEUD* tris2[2][3]={{no[3],no[1],no[2]},{no[1],no[3],no[0]}};
174			int testOverlapTriangles = CheckOverlapTriangle(tris2);
175			score1 = ( testOverlapTriangles != 0 ) ? -1 : score1;
176			if (score1 < 0) return score1;
177
178			double diedralAngle2 = GeometricTools::ComputeDiedralAngle(tris2);
179			score1 = ( diedralAngle2 > maxDiedralAngle ) ? -1 : score1;
180			if (score1 < 0) return score1;
181
182			if (__execute && score1 > 0)
183			{
184			double crit21 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[3]->get_coord(),no[0]->get_coord());
185			double crit22 = OPERATEUR::qualite_triangle(no[1]->get_coord(),no[2]->get_coord(),no[3]->get_coord());
186			MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
187			MCSegment * swapSegment = new MCSegment(topo, (MCNode)no[1], (MCNode)no[3] );
188			_mesh->ajouter_mg_segment(swapSegment);
189			M3D_MCTriangle * newTriangle[2];
190			newTriangle[0] = new M3D_MCTriangle(topo, (MCNode)no[1], (MCNode)no[2], (MCNode)no[3], (MCSegment)seg[1], (MCSegment*)seg[2], swapSegment);
191			newTriangle[0]->change_qualite(crit21);
192			newTriangle[1] = new M3D_MCTriangle(topo, (MCNode)no[3], (MCNode)no[0], (MCNode)no[1], (MCSegment)seg[3], (MCSegment*)seg[0], swapSegment);
193			newTriangle[1]->change_qualite(crit22);
194			for (int i=0; i<2; i++)_mesh->ajouter_mg_triangle(newTriangle[i]);
195			for (int i=0; i<2; i++)_mesh->supprimer_mg_triangleid(triangle[i]->get_id());
196			_mesh->supprimer_mg_segmentid(__segment->get_id());
197			}
198
199			return score1;
200			} */
201
202
203			double OptimizeEdgeSwap::Swap(MG_SEGMENT *__segment, bool __execute)
204			{
205			double score1;
206			if (__segment->get_lien_triangle()->get_nb() != 2)
207			{
208			return -1;
209			}
210			if (__segment->get_lien_triangle()->get(0)->get_lien_topologie() != __segment->get_lien_triangle()->get(1)->get_lien_topologie())
211			{
212			return -2;
213			}
214			if (__segment->get_lien_topologie()&&__segment->get_lien_topologie()->get_dimension()!=2)
215			{
216			return -3;
217			}
218
219			MG_TRIANGLE * triangle[2]={__segment->get_lien_triangle()->get(0),__segment->get_lien_triangle()->get(1)};
220			MG_NOEUD * no[4]; // no[]={A,C,B,D} where AB is the segment shared by both triangles AB, C and D nodes opposed to segment AB
221			no[2] = GetStartNode(triangle[0],__segment);
222			no[1] = GetNodeBefore(triangle[0],no[2]);
223			no[0] = GetNodeBefore(triangle[0],no[1]);
224			no[3] = GetNodeBefore(triangle[1],no[0]);
225			MG_SEGMENT * seg[4];
226			seg[0] = _mesh->get_mg_segment(no[0]->get_id(),no[1]->get_id()); // AC
227			seg[1] = _mesh->get_mg_segment(no[1]->get_id(),no[2]->get_id()); // CB
228			seg[2] = _mesh->get_mg_segment(no[2]->get_id(),no[3]->get_id()); // BD
229			seg[3] = _mesh->get_mg_segment(no[3]->get_id(),no[0]->get_id()); // DA
230
231			MG_NOEUD* tris[2][2][3]={{{no[2],no[3],no[0]},{no[0],no[1],no[2]}}, // triangles of the current configuration
232			{{no[3],no[1],no[2]},{no[1],no[3],no[0]}}}; // triangles of the candidate configuration (after segment swapping)
233
234
235			int testOverlapTriangles = CheckOverlapTriangle(tris[1]);
236			if (testOverlapTriangles) return -1;
237
238			// projection des triangles sur le plan tangent
239			double crit[2][2];
240			double area[2][2];
241			for (int k=0;k<2;k++)
242			for (int i=0;i<2;i++)
243			{
244			OT_VECTEUR_3D normal[4];
245			MCFace * face = (MCFace *)triangle[i]->get_lien_topologie();
246			for (int j=0;j<3;j++)
247			((MCNode*)tris[k][i][j])->NormalMCFace(face,normal[j]);
248			normal[3] = .3333*(normal[0]+normal[1]+normal[2]);
249			// frame (x,y,n)
250			OT_MATRICE_3D tangentPlane = GeometricTools::GetPlaneFrame(normal[3]);
251			OT_MATRICE_3D trf; tangentPlane.transpose(trf);
252			OT_VECTEUR_3D tpProjNodes[3];
253			for (int j=0;j<3;j++)
254			{
255			OT_VECTEUR_3D p(tris[k][i][j]->get_coord());
256			tpProjNodes[j] = trf*p;
257			// projection --> n component = 0
258			tpProjNodes[j][2] = 0;
259			}
260			crit[k][i]=OPERATEUR::qualite_triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
261			area[k][i]=OT_ALGORITHME_GEOMETRIQUE::Area_Triangle(tpProjNodes[0],tpProjNodes[1],tpProjNodes[2]);
262			}
263			double crit1=std::min(crit[0][0],crit[0][1]);
264			double crit2=std::min(crit[1][0],crit[1][1]);
265
266			// dist_P[i] = la distance entre la diagonale de la ieme configuration et la MCFace
267			// voir page 110 de these_jfl.pdf
268			double critGeomAccuracy[2];
269			{
270			MG_ELEMENT_TOPOLOGIQUE * topo = __segment->get_lien_triangle()->get(0)->get_lien_topologie() ;
271			MCSegment * AB = (MCSegment*)__segment;
272			MCNode Pgeo;
273			// MCSegment * CD = new MCSegment(topo, (MCNode)no[1], (MCNode)no[3] );
274			// CD->evaluer_geo(((area[0][0]0+area[0][1]1)/(area[0][0]+area[0][1])),P[0]);
275			{
276			double wA, wB;
277			wA=(AB->get_noeud1() == no[0])?0:1;
278			if(AB->get_orientation() == -1) wA=(wA==0)?1:0;
279			wB=(wA==0)?1:0;
280			AB->evaluer_geo(((area[1][0]wA+area[1][1]wB)/(area[1][0]+area[1][1])),&Pgeo);
281			}
282
283			double P_Euclid[2][3];
284			for (int i=0;i<3;i++) {
285			P_Euclid[1][i]=((area[0][0]no[1]->get_coord()[i]+area[0][1]no[3]->get_coord()[i])/(area[0][0]+area[0][1]));
286			P_Euclid[0][i]=((area[1][0]no[0]->get_coord()[i]+area[1][1]no[2]->get_coord()[i])/(area[1][0]+area[1][1]));
287			}
288			double distP[2],segLength[2];
289			for (int i=0;i<2;i++)
290			distP[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(P_Euclid[i],Pgeo.get_coord());
291			for (int i=0;i<2;i++)
292			segLength[i]=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(no[i]->get_coord(),no[i+2]->get_coord());
293			double maxDistP=.5*std::max(segLength[0],segLength[1]);
294			for (int i=0;i<2;i++)
295			if ( distP[i] < maxDistP )
296			critGeomAccuracy[i] = 1-distP[i]/maxDistP;
297			else
298			critGeomAccuracy[i] = 0;
299			}
300
301			crit1 += critGeomAccuracy[0];
302			crit2 += critGeomAccuracy[1];
303
304			if (__execute && crit2 > crit1 )
305			{
306			SwapSegment(__segment);
307			}
308
309			return (crit2>crit1)?crit2:-1;
310			}
311
312
313			MG_NOEUD * OptimizeEdgeSwap::GetStartNode(MG_TRIANGLE __t, MG_SEGMENT __s)
314			{
315			MG_NOEUD *TN[3];
316			TN[0] = __t->get_noeud1();
317			TN[1] = __t->get_noeud2();
318			TN[2] = __t->get_noeud3();
319			for (int i=0; i<3; i++)
320			{
321			if (__s->get_noeud1() == TN[i])
322			{
323			if (__s->get_noeud2() == TN[(i+1)%3])
324			return __s->get_noeud1();
325			else
326			return __s->get_noeud2();
327			}
328			}
329			return NULL;
330			}
331
332			MG_NOEUD * OptimizeEdgeSwap::GetNodeAfter(MG_TRIANGLE __t, MG_NOEUD __n)
333			{
334			MG_NOEUD *TN[3];
335			TN[0] = __t->get_noeud1();
336			TN[1] = __t->get_noeud2();
337			TN[2] = __t->get_noeud3();
338			for (int i=0; i<3; i++)
339			{
340			if (TN[i]==__n)
341			return TN[(i+1)%3];
342			}
343			return NULL;
344			}
345
346			MG_NOEUD * OptimizeEdgeSwap::GetNodeBefore(MG_TRIANGLE __t, MG_NOEUD __n)
347			{
348			MG_NOEUD *TN[3];
349			TN[0] = __t->get_noeud1();
350			TN[1] = __t->get_noeud2();
351			TN[2] = __t->get_noeud3();
352			for (int i=0; i<3; i++)
353			{
354			if (TN[i]==__n)
355			return TN[(i+2)%3];
356			}
357			return NULL;
358			}
359
360			MG_NOEUD * OptimizeEdgeSwap::GetOppositeNode(MG_TRIANGLE __t, MG_SEGMENT __s)
361			{
362			MG_NOEUD *TN[3];
363			TN[0] = __t->get_noeud1();
364			TN[1] = __t->get_noeud2();
365			TN[2] = __t->get_noeud3();
366			for (int i=0; i<3; i++)
367			if (TN[i] != __s->get_noeud1() && TN[i] != __s->get_noeud2())
368			return TN[i];
369			return NULL;
370			}
371
372			void OptimizeEdgeSwap::RemoveVolumeMesh()
373			{
374			LISTE_MG_SEGMENT::iterator itSeg;
375			LISTE_MG_NOEUD::iterator itNo;
376			LISTE_MG_TRIANGLE::iterator itTri;
377			_mesh->supprimer_tout_mg_tetra();
378			for (MG_TRIANGLE * tri = _mesh->get_premier_triangle(itTri);tri; tri = _mesh->get_suivant_triangle(itTri))
379			{
380			if (tri->get_lien_topologie()->get_dimension()==3)
381			_mesh->supprimer_mg_triangleid(tri->get_id());
382			}
383			for (MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
384			{
385			if (seg->get_lien_topologie()->get_dimension()==3)
386			_mesh->supprimer_mg_segmentid(seg->get_id());
387			}
388			for (MG_NOEUD * no = _mesh->get_premier_noeud(itNo);no; no = _mesh->get_suivant_noeud(itNo))
389			{
390			if (no->get_lien_topologie()->get_dimension()==3)
391			_mesh->supprimer_mg_noeudid(no->get_id());
392			}
393			}
394
395			int OptimizeEdgeSwap::Optimize(MG_FACE * __face)
396			{
397			int nbSwap = 0;
398			LISTE_MG_SEGMENT::iterator itSeg;
399			MG_SEGMENT * seg = _mesh->get_premier_segment(itSeg);
400			std::map<double, MG_SEGMENT*> lst;
401
402			for (int j=0;j<3;j++)
403			{
404			lst.clear();
405			int tmpnbSwap = nbSwap;
406			for (seg = _mesh->get_premier_segment(itSeg);seg; seg = _mesh->get_suivant_segment(itSeg))
407			{
408			MG_ELEMENT_TOPOLOGIQUE * topo = seg->get_lien_topologie();
409			if (__face != NULL && topo != __face)
410			continue;
411			if (topo->get_dimension()!=2)
412			continue;
413			lst[SwapScore(seg)]=seg;
414			}
415
416			for (std::map<double, MG_SEGMENT*>::reverse_iterator i=lst.rbegin();i!=lst.rend();i++)
417			{
418			double score1 = i->first;
419			seg=i->second;
420
421			if (score1<0)
422			break;
423
424			double score2 = SwapSegment(seg);
425			if (score2 > 0)
426			{
427			nbSwap++;
428			}
429			}
430			if (nbSwap-tmpnbSwap==0)
431			break;
432			}
433
434			return nbSwap;
435			}
436
437			int OptimizeEdgeSwap::OptimizeAllFaces()
438			{
439			return Optimize(0);
440			}
441