CAD4FE_OptimizeEdgeSwap.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_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);
}