CAD4FE/src/CAD4FE_MCSegment.cpp

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

#include <fstream>

#pragma hdrstop

#include "gestionversion.h"
#include "CAD4FE_MCSegment.h"      
#include "CAD4FE_MCFace.h"         
#include "CAD4FE_MCEdge.h"

#include "ot_algorithme_geometrique.h" 
#include "ot_decalage_parametre.h"
#include "CAD4FE_MCNode.h"  
#include "CAD4FE_PolySurface.h"        
#include "CAD4FE_PolyCurve.h"
#include "CAD4FE_MCVertex.h"    
#include "CAD4FE_MCEdge.h"
#include "CAD4FE_Intersection_Plane_PolySurface.h"
#include "CAD4FE_ShortestPathByClosestPointOnEdge.h"     
#include "CAD4FE_MCNodePolyline.h"                     
#include "CAD4FE_Geometric_Tools.h"
#include "CAD4FE_InventorText_MG_MAILLAGE.h"
#include "ot_root_find.h"
#include "CAD4FE_mailleur2d.h"

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

#pragma package(smart_init)

using namespace CAD4FE;

//---------------------------------------------------------------------------
MCSegment::MCSegment (unsigned long num, MG_ELEMENT_TOPOLOGIQUE *topo, MCNode *mgnoeud1, MCNode *mgnoeud2, double longue)
: MG_SEGMENT(num,topo,mgnoeud1,mgnoeud2,longue),_polylineEvaluator(0),_sense(1), _saveFormat(0)
{
    _cstrMCFaceScheme = CstrMCFaceByPlaneIntr;
    ConstructGeometry();
}     
//---------------------------------------------------------------------------
MCSegment::MCSegment (MG_ELEMENT_TOPOLOGIQUE *topo, MCNode *mgnoeud1, MCNode *mgnoeud2, double longue, long __cstrMCFaceScheme)
: MG_SEGMENT(topo,mgnoeud1,mgnoeud2,longue),_polylineEvaluator(0),_sense(1),_cstrMCFaceScheme(__cstrMCFaceScheme), _saveFormat(0)
{
    ConstructGeometry();
}
//---------------------------------------------------------------------------
MCSegment::MCSegment( MCSegment & __s)
: MG_SEGMENT(0,__s.liaison_topologique,__s.get_noeud1(),__s.get_noeud2(),__s.get_longueur()), _saveFormat(0)
{
    unsigned i;
    for (FMapCIterator itF = __s.F.begin(); itF != __s.F.end(); itF++)
    {
        MG_FACE * face = itF->first;
        F.insert(std::make_pair(face,itF->second));
    }
    for (EMapCIterator itE = __s.E.begin(); itE != __s.E.end(); itE++)
    {
        MG_ARETE * edge = itE->first;
        E.insert(std::make_pair(edge,itE->second));
    }
    //V = __s.V;
    _polylineEvaluator=new MCNodePolyline(*__s._polylineEvaluator);
    _sense = __s._sense;
}                     
//---------------------------------------------------------------------------
MCSegment::FMap & MCSegment::GetRefFaceMapping()
{
    return F;
}
//---------------------------------------------------------------------------
MCSegment::EMap & MCSegment::GetRefEdgeMapping()
{
    return E;
}
//---------------------------------------------------------------------------
MCSegment::~MCSegment()
{
    unsigned i,j;       
    if (_polylineEvaluator) delete _polylineEvaluator;
}
//---------------------------------------------------------------------------
std::vector <MCNode*> MCSegment::GetPolylineNodes()
{
    return _polylineEvaluator->GetPolylineNodes();
}
//---------------------------------------------------------------------------
MCNode* MCSegment::GetPolylineNode(unsigned __index)
{
    return _polylineEvaluator->GetPolylineNode(__index);
}
//---------------------------------------------------------------------------
std::vector <MG_ELEMENT_TOPOLOGIQUE*> MCSegment::GetPolylineTopos()
{
    return _polylineEvaluator->GetPolylineTopos();
}
//---------------------------------------------------------------------------
MG_ELEMENT_TOPOLOGIQUE* MCSegment::GetPolylineTopo(unsigned __index)
{
    return _polylineEvaluator->GetPolylineTopo(__index);
}
//---------------------------------------------------------------------------
unsigned MCSegment::GetPolylineNodeCount()
{
    return _polylineEvaluator->GetPolylineNodes().size();
}
//---------------------------------------------------------------------------
double MCSegment::get_longueur_geo()
{
    double length = _polylineEvaluator->GetLength();
    return length;
}                                                                   
//---------------------------------------------------------------------------
void MCSegment::evaluer_geo(double __t, MCNode * __result, double * tangent, double * curvature)
{
    double t,L;
    if (_sense == -1)
        t = (1-__t);
    else
        t = __t;

    L=get_longueur_geo();
    * __result = _polylineEvaluator->Evaluate( t * L, tangent, curvature );
    if (tangent)
        for (int i=0; i<3; i++) tangent[i] *= L;
    if (_sense == -1)
    {
        if (tangent)
            for (int i=0; i<3; i++) tangent[i]*=-1;
    }
}
//---------------------------------------------------------------------------
void MCSegment::inverser_geo(double & __t, MCNode * __node)
{
    _polylineEvaluator->Inverse(__t, __node);
    double length = _polylineEvaluator->GetLength();
    __t = __t/length;

    if (_sense == -1)
        __t = (1-__t);
}                                      
//---------------------------------------------------------------------------
int MCSegment::get_orientation()
{
    return _sense;
}                                       
//---------------------------------------------------------------------------
void MCSegment::change_orientation(int __sense)
{
    _sense = __sense;
}                            
//---------------------------------------------------------------------------
void MCSegment::UpdateGeometry ()
{              
    calcule_longueur();
    F.clear();
    E.clear();
    //V.clear();
    delete _polylineEvaluator;
    _polylineEvaluator = 0;
    ConstructGeometry ();
}
//---------------------------------------------------------------------------
void MCSegment::ConstructGeometry ()
{
    MCNode * mcNode1 = (MCNode*)get_noeud1();
    MCNode * mcNode2 = (MCNode*)get_noeud2();
    CAD4FE::MCFace * mcFace = NULL;
    CAD4FE::MCEdge * mcEdge = NULL;
    if (get_lien_topologie()->get_dimension() == 2)
    {
        mcFace = (MCFace*)get_lien_topologie();
    }
    if (get_lien_topologie()->get_dimension() == 1)
    {
        mcEdge = ((CAD4FE::MCEdge*)get_lien_topologie());
    }
    if (mcEdge)
    {
        Construct_MCEdge(mcEdge);
        Construct_MergedVertices(mcNode1);
        Construct_MergedVertices(mcNode2);
    }
    else if (mcFace)
    {
        int res;
        if (_cstrMCFaceScheme == CstrMCFaceByPlaneIntr )
        {
            res = Construct_MCFace(mcFace,mcNode1,mcNode2);
            if (res != 1 )
            {
                res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
                if (res == 1) _cstrMCFaceScheme = CstrMCFaceByShortestPath;
            }
        }
        else if (_cstrMCFaceScheme == CstrMCFaceByPlaneIntrNormalByShortestPath)
        {
            res = Construct_MCFace(mcFace,mcNode1,mcNode2);
            if (res != 1)
            {                                
                res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
                if (res == 1)
                    _cstrMCFaceScheme = CstrMCFaceByShortestPath;
            }
        }
        else if (_cstrMCFaceScheme == CstrMCFaceByShortestPath )
        {
            res = Construct_MCFaceByShortestPath(mcFace,mcNode1,mcNode2);
            if (res != 1)
            {  
                res = Construct_MCFace(mcFace,mcNode1,mcNode2);
                if (res == 1)
                    _cstrMCFaceScheme = CstrMCFaceByPlaneIntr;
            }
        }
        else if (_cstrMCFaceScheme == CstrNoRefTopologyMapping)
        {
            res = Construct_NoTopology(mcNode1,mcNode2);
        }
    }

    _sense = 1;
    ConstructBoundaryBox();
}

void MCSegment::ConstructBoundaryBox()
{
    for (unsigned j=0; j<_polylineEvaluator->GetPolylineNodes().size();j++)
    {
        MCNode * n = _polylineEvaluator->GetPolylineNode(j);
        double * x = n->get_coord();
        if (j == 0)
        {
            for (int i=0;i<3;i++)
                _bbox[i] = _bbox[i+3] = x[i];
        }
        else
        {
            for (int i=0;i<3;i++)
            {
                if (_bbox[i]>x[i])
                    _bbox[i] = x[i];
                else if (_bbox[i+3]<x[i])
                    _bbox[i+3] = x[i];
            }
        }
    }
    for (double t=.25; t<=.75; t+=.25)
    {
        MCNode n = _polylineEvaluator->Evaluate(t);
        double * x = n.get_coord();

            for (int i=0;i<3;i++)
            {
                if (_bbox[i]>x[i])
                    _bbox[i] = x[i];
                else if (_bbox[i+3]<x[i])
                    _bbox[i+3] = x[i];
            }
    }
}

extern MG_MAILLAGE * debugMesh;

int MCSegment::Construct_MCFaceByShortestPath(MCFace * __mcFace, MCNode* __n1, MCNode* __n2)
{
    unsigned i,j;
    double dN1N2=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n1->get_coord(),__n2->get_coord());
    int shortestPathOptions = 0;
    ShortestPathByClosestPointOnEdge shortestPath(__mcFace,__n1,__n2,dN1N2, shortestPathOptions);
    ShortestPathByClosestPointOnEdge *shortestPathDoubleDist = NULL;
    std::vector<MCNode*> vecPts;
    std::vector<MG_ELEMENT_TOPOLOGIQUE*> vecTopo;
    double distance = shortestPath.Find(&vecPts,&vecTopo);
    
    if (distance > 100) // failure, second try with 2*dist(n1,n2)
    {
        unsigned maxNbTries = 4;
        unsigned maxDistFactor = 1;
        for (i=0; i<maxNbTries; i++)
        {
            maxDistFactor *= 2;
            vecPts.clear();vecTopo.clear();
            shortestPathDoubleDist = new ShortestPathByClosestPointOnEdge (__mcFace,__n1,__n2,maxDistFactor*dN1N2, shortestPathOptions);
            distance = shortestPathDoubleDist->Find(&vecPts,&vecTopo);
            if (distance > 100)
            {
                delete shortestPathDoubleDist;     
                vecPts.clear();vecTopo.clear();
                shortestPathOptions=ShortestPathByClosestPointOnEdge::NoCheckFaceIntr;
                shortestPathDoubleDist = new ShortestPathByClosestPointOnEdge (__mcFace,__n1,__n2,maxDistFactor*dN1N2, shortestPathOptions);
                distance = shortestPathDoubleDist->Find(&vecPts,&vecTopo);
                if (distance > 100)
                {
                    delete shortestPathDoubleDist;
                    shortestPathDoubleDist = NULL;
                }
            }
            else
            {
                break;
            }
        }
    }

    if (distance > 100)
    {
        printf("Warning : Construct_MCFaceByShortestPath failed !\n");
        distance = shortestPath.Find(&vecPts,&vecTopo);
        Construct_MCFace(__mcFace,__n1,__n2);
        return 0;
    }
    if (distance < .8*dN1N2)
    {                
        printf("Warning : %d < .8*%f in Construct_MCFaceByShortestPath !\n", distance, dN1N2);
    }

    if (_polylineEvaluator == NULL) _polylineEvaluator = new MCNodePolyline(liaison_topologique);

    for (i=0;i+1<vecPts.size();i++)
    {
        MCNode * xi=vecPts[i], *xii=vecPts[i+1];
        MG_ELEMENT_TOPOLOGIQUE * topo = vecTopo[i];

        if (topo->get_dimension() == 2)
        {
                MG_FACE * face = (MG_FACE *) topo;
                std::vector<MCNode*> Ci;
                OT_VECTEUR_3D normalXi,normalXii,normalAvg,dir,v,xi_xyz(xi->get_coord()),xii_xyz(xii->get_coord());
                xi->NormalMCFace(__mcFace, normalXi);
                xii->NormalMCFace(__mcFace, normalXii);
                normalAvg = .5*normalXi+.5*normalXii;
                dir=xii_xyz-xi_xyz;
                v=normalAvg&dir;
                double L=dir.get_longueur();
                Intersection_Plane_MG_FACE * intr = NULL;
                if (L > 1E-6 && ((shortestPathOptions&ShortestPathByClosestPointOnEdge::NoCheckFaceIntr)==0))
                {
                    intr = new Intersection_Plane_MG_FACE(xi,v,face,__mcFace);
                    intr->MakeSegment(xi,xii,L/50,L/3000,3,dir);
                    Ci = intr->Get()[0];
                    if (Ci.size()<2||Ci[Ci.size()-1]!=xii||Ci[0]!=xi)
                    {
                        delete intr;
                        intr = 0;
                        Ci.clear();
                        Ci.push_back(xi);
                        Ci.push_back(xii);
                    }
                }
                else
                {
                    Ci.push_back(xi);
                    Ci.push_back(xii);
                }
                F.insert( std::make_pair(face, Ci) );

                for (j=0; j<Ci.size(); j++)
                {
                    _polylineEvaluator->Add(Ci[j],face);
                }
                if (intr != NULL)
                    delete intr;
        }
        if (topo->get_dimension() == 1)
        {
            MG_ARETE * edge = (MG_ARETE *) topo;

            //double dist  = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xi->get_coord(), xii->get_coord());
            //double epsilon = dist*1E-2;
            //int insertedNodeCount = RefineInRefEdge(__mcFace, edge, xi, xii, epsilon, 1E300, refinedNodes, refinedEdges );

            VPts Ci;
            Ci.push_back(xi);
            Ci.push_back(xii);

            E.insert (std::make_pair(edge, Ci));
            for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
            {
                MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
                F.insert(std::make_pair(face, Ci));
            }

            for (j=0; j<Ci.size(); j++)
            {
                _polylineEvaluator->Add(Ci[j],edge);
            }
        }  
        if (topo->get_dimension() == 0)
        {
                
        }
    }

    if (shortestPathDoubleDist != NULL)
        delete shortestPathDoubleDist;

    return 1;
}

int MCSegment::Construct_NoTopology(MCNode * __n1, MCNode * __n2)
{
    _polylineEvaluator = new MCNodePolyline(liaison_topologique);
    _polylineEvaluator->Add(__n1,NULL);
    _polylineEvaluator->Add(__n2,NULL);
    return 1;
}

int MCSegment::Construct_MCFace(MCFace * __mcFace, MCNode* __n1, MCNode* __n2)
{
    unsigned i,j;
    PolySurface * polysurface = __mcFace->GetPolySurface();
                                                  
    OT_VECTEUR_3D normal1, normal2, xyz1, xyz2, v1, v2;    
    if ( _cstrMCFaceScheme != CstrMCFaceByPlaneIntrNormalByShortestPath )
    {    
            int nbRefFaces;
            polysurface->calcul_normale_unitaire(__n1->GetRefFaceMapping(),normal1,&nbRefFaces);
            polysurface->calcul_normale_unitaire(__n2->GetRefFaceMapping(),normal2,&nbRefFaces);
            v2 = .5*(normal1+normal2);
            if (v2.get_longueur2() < 1E-100)
            {
                printf("Warning: surface's normal at node 1 is opposite to surface's normal at node 2\n");
                printf("\t trying to evaluate surface's normal with the shortest path algorithm\n"); 
                _cstrMCFaceScheme = CstrMCFaceByPlaneIntrNormalByShortestPath;
            }
    }
    if (_cstrMCFaceScheme==CstrMCFaceByPlaneIntrNormalByShortestPath)
    {
        double maxDistBetween3DSegmentAndShortestPath=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n1->get_coord(),__n2->get_coord());
        ShortestPathByClosestPointOnEdge * shortestPath = new ShortestPathByClosestPointOnEdge(__mcFace,__n1,__n2,maxDistBetween3DSegmentAndShortestPath);
        std::vector<MCNode*> vecPts;
        std::vector<MG_ELEMENT_TOPOLOGIQUE*> vecTopo;
        double distance=shortestPath->Find(&vecPts, &vecTopo);
        if (distance > 100)
        {
//            delete shortestPath;
//            vecPts.clear(); vecTopo.clear();
//            shortestPath = new ShortestPathByClosestPointOnEdge(__mcFace,__n1,__n2,100);
            double *x1=__n1->get_coord(),*x2=__n2->get_coord();
            printf("%f %f %f %f %f %f ",x1[0],x1[1],x1[2],x2[0],x2[1],x2[2]);
//            double distance = shortestPath->Find(&vecPts, &vecTopo);
            if (distance > 100)
                printf("Error in average normal\n");

        }
        if (distance == 0.0)
        {            
            delete shortestPath;
            return 0;
        }
        // estimation of the face normal by the average of the normals
        // in the polyline
        OT_VECTEUR_3D avgNormal(0,0,0);
        MCNodePolyline polyline(liaison_topologique,vecPts,vecTopo);
        unsigned nb_points=10, real_nb_points=0;
        double s_max=polyline.GetLength();
        for (double s=0; real_nb_points < 1 || s < s_max; s+=.9999*(s_max/nb_points))
        {
            real_nb_points++;
            MCNode mcNode=polyline.Evaluate(s);
            OT_VECTEUR_3D tmpNormal;
            int nbRefFaces;
            polysurface->calcul_normale_unitaire(mcNode.GetRefFaceMapping(),tmpNormal,&nbRefFaces);
            avgNormal += tmpNormal;
        }
        avgNormal /= real_nb_points;
        v2=avgNormal;
        delete shortestPath;
    }

    xyz1=OT_VECTEUR_3D(__n1->get_coord());
    xyz2=OT_VECTEUR_3D(__n2->get_coord());
    v1=xyz2-xyz1;

    // Make a special treatment to manage nodes of merged vertices
    MCNode * polysurfaceNodes[2]={__n1,__n2};
    MCNode * n[2]={__n1,__n2};
    VPts MergedVertexCurveNodes[2];
    std::vector < MG_ELEMENT_TOPOLOGIQUE * > MergedVertexCurveTopo[2];
    for (int i=0;i<2;i++)
    if (n[i]->GetRefVertexMapping().size() > 1) // if node is in a merged vertex
    {
        MCVertex * mcVertex = (MCVertex*) n[i]->get_lien_topologie();
        MG_SOMMET * refVertex = (MG_SOMMET*) n[i]->get_lien_topologie_reference();
        MG_SOMMET * mergedVertex = NULL;
        if (GeometricTools::PolySurface_Contains_RefVertex(polysurface,refVertex) == false)
        {
            for (MCNode::VMapIterator itV = n[i]->GetRefVertexMapping().begin();
                itV != n[i]->GetRefVertexMapping().end();
                itV++)
                if ( GeometricTools::PolySurface_Contains_RefVertex(polysurface,*itV) )
                {
                    mergedVertex = *itV;
                    break;
                }

            double xyz[3];mergedVertex->get_point()->evaluer(xyz);
            double d=OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(xyz,n[i]->get_coord());
            MCNode * mergedVertexNode = new MCNode(mcVertex,mergedVertex,xyz[0],xyz[1],xyz[2]);
            polysurfaceNodes[i] = mergedVertexNode;

            MCFace * mcFacePath = NULL;
            std::set<MG_FACE*> mergedVertexMCFaces;
            GeometricTools::MG_SOMMET_GetAdjacent_MG_FACE(mcVertex,mergedVertexMCFaces);
            for (std::set<MG_FACE*>::iterator itF=mergedVertexMCFaces.begin();
                itF != mergedVertexMCFaces.end(); itF++)
            {
                MCFace* adjacentMCFace = (MCFace*)(*itF);
                if (GeometricTools::PolySurface_Contains_RefVertex(adjacentMCFace->GetPolySurface(),mergedVertex)
                && GeometricTools::PolySurface_Contains_RefVertex(adjacentMCFace->GetPolySurface(),refVertex)
                )
                {
                    mcFacePath = adjacentMCFace;
                    break;
                }
            }
            ShortestPathByClosestPointOnEdge shortestPath(mcFacePath,n[i],mergedVertexNode,d);
            double distance = shortestPath.Find(&(MergedVertexCurveNodes[i]),&(MergedVertexCurveTopo[i]));

            if (distance < 100)
                for (unsigned j=0; j<MergedVertexCurveNodes[i].size();j++)
                    MergedVertexCurveNodes[i][j]->incrementer();
        }
    }

    Intersection_Plane_PolySurface intrPolySurface (polysurfaceNodes[0], v2&v1, __mcFace, polysurfaceNodes[1]);
    MCNodePolyline * intersectionPlaneMCFace = intrPolySurface.MakeSegment(polysurfaceNodes[0], polysurfaceNodes[1], v1);

    if (intersectionPlaneMCFace != NULL)
    {
        _polylineEvaluator = intersectionPlaneMCFace;
        _polylineEvaluator->incrementer();

        Intersection_Plane_PolySurface::FMapIterator it;
        Intersection_Plane_MG_FACE * intr = NULL;
        for (intr = intrPolySurface.GetFirstCurve(it); intr; intr = intrPolySurface.GetNextCurve(it) )
        {
            MG_FACE * face = intr->GetFace();
            VPts nodes = intr->Get()[0];
            F.insert(std::make_pair(face, nodes));
        }

        for ( unsigned i=0;1+i<GetPolylineNodeCount();i++)
        {
            MG_ELEMENT_TOPOLOGIQUE * topo = GetPolylineTopo(i);
            if ( topo != NULL && topo->get_dimension()==1)
            {
                MCNode * xi = GetPolylineNode(i);
                MCNode * xii = GetPolylineNode(i+1);
                MG_ARETE * edge = (MG_ARETE *) topo;
                VPts Ci;
                Ci.push_back(xi);
                Ci.push_back(xii);

                E.insert (std::make_pair(edge, Ci));
                for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
                {
                    MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
                    F.insert(std::make_pair(face, Ci));
                }
            }
        }

        return 1;
    }    
    else
    {
        return 0;
    }
}

int MCSegment::RefineInRefEdge(MCEdge * __mcEdge, MG_ARETE * __refEdge, MCNode * __n1, MCNode * __n2, double __epsilon, double __lastDistanceToEdge, std::map < double , MCNode * > & __mapMCEdgeNodesBySParameter, std::map < double , MG_ARETE * > & __mapMCEdgeRefEdgesBySParameter, int __iRecursiveCalls)
{
    static const int RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT = 100;
    if (__iRecursiveCalls > RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT)
    {
        printf("Warning : RefineInRefEdge : The number of recursive calls has been reached (%d)\n",RefineInRefEdge_MAX_RECURSIVE_CALL_COUNT);
        return 0;
    }

    int insertedNodeCount = 0;
    MCNode * n[2] = { __n1, __n2};
    for (int i=0; i<2; i++)
    {
        MG_SOMMET * v = 0;
        if (n[i]->get_lien_topologie_reference()->get_dimension() == 0)
            v = (MG_SOMMET * ) n[i]->get_lien_topologie_reference();
        if (v)
        if (__mcEdge->GetPolyCurve()->Contains( v ) == false )
            return 0;
        if (n[i]->IsInEdge(__refEdge) == false)
            return 0;
        // nodes have merged ref vertices
        if (n[(i+1)%2]->IsInVertex(v))
            return 0;
    }                               
    __n1->incrementer();__n2->incrementer();
    double t1 = __n1->T(__refEdge);
    double t2 = __n2->T(__refEdge);
    if (fabs(t2-t1) < 1E-6*(__refEdge->get_tmax()-__refEdge->get_tmin()))
        return 0;

    double period = __refEdge->get_courbe()->get_periode();
    if (period != 0.0 && fabs(t2 - t1) > .5*period)
        if (t2<t1)
            t2 += period;
        else t1 += period;
    double t3 = t1 + (t2-t1)*.5;
    double s3;
    double xyz[3]; __refEdge->evaluer(t3,xyz);
    MCNode edgeNode (__mcEdge, __refEdge, t3, xyz);
    for (int i=0; i<__refEdge->get_nb_mg_coarete(); i++)
    {
        MG_FACE * face = __refEdge->get_mg_coarete(i)->get_boucle()->get_mg_face();

        MCNodePolyline faceLine(liaison_topologique);
        faceLine.Add(__n1,face);
        faceLine.Add(__n2,face);
        double L =  faceLine.GetLength();
        MCNode faceNode = faceLine.Evaluate(L*.5);
        double distanceToEdge = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(faceNode.get_coord(),edgeNode.get_coord());
        double ratioConvergence = distanceToEdge / __lastDistanceToEdge;
        double distN1N2 = OT_ALGORITHME_GEOMETRIQUE::VEC3_DISTANCE_VEC3(__n2->get_coord(),__n1->get_coord());
        bool testConvergence = ( fabs(ratioConvergence - 1) < .001 || distN1N2 < __epsilon );
        if ( (distanceToEdge > __epsilon  )
            && testConvergence == false )
        {
            insertedNodeCount++;
            __mcEdge->GetPolyCurve()->Parameter_RefEdgeTToS ( t3, __refEdge, &s3, false);
            if (__mcEdge->GetPolyCurve()->est_periodique()
            && __mcEdge->get_tmin() > s3+1E-5*(__mcEdge->GetPolyCurve()->get_periode())
               )
                s3  += __mcEdge->GetPolyCurve()->get_periode();
            MCNode * newMCNode = new MCNode ( edgeNode );
            __mapMCEdgeNodesBySParameter.insert(std::make_pair(s3,newMCNode));
            __mapMCEdgeRefEdgesBySParameter.insert(std::make_pair(s3, __refEdge));
            newMCNode->incrementer();
            insertedNodeCount += RefineInRefEdge(__mcEdge, __refEdge, __n1, newMCNode, __epsilon, distanceToEdge, __mapMCEdgeNodesBySParameter,__mapMCEdgeRefEdgesBySParameter, __iRecursiveCalls+1);
            insertedNodeCount += RefineInRefEdge(__mcEdge, __refEdge, newMCNode, __n2, __epsilon, distanceToEdge, __mapMCEdgeNodesBySParameter,__mapMCEdgeRefEdgesBySParameter, __iRecursiveCalls+1);
            newMCNode->decrementer();
            break;
        }
    }
    __n1->decrementer();__n2->decrementer();
    return insertedNodeCount;
}

int MCSegment::Construct_MCEdge(MCEdge * __mcEdge)
{
    double s[2];        
    double t,dt;
    int i,j;
    unsigned iEdge;
    MG_ARETE * edge;
    PolyCurve * polycurve = __mcEdge->GetPolyCurve();
    double period = polycurve->get_periode();
    double smin = __mcEdge->get_tmin();
    double smax = __mcEdge->get_tmax();

    VPts CurveNodes;
    std::vector < MG_ELEMENT_TOPOLOGIQUE * > CurveTopo;

    MCNode * n[2];
    n[0]=(MCNode *)get_noeud1();
    n[1]=(MCNode *)get_noeud2();

    std::map < double , MCNode * > nodes;
    std::map < double , MG_ARETE * > edges;

    // get the parameter of extremities n1 n2 of the segment
    for (i=0; i<2; i++)
    {
        if (n[i]->GetRefVertexMapping().size() >= 1) // MC node is on a ref vertex
        {
            MG_SOMMET * refVertex = (MG_SOMMET*)n[i]->get_lien_topologie_reference();
            if (polycurve->ContainsRefVertex(refVertex) == false)
            {
                for (std::set<MG_SOMMET*>::iterator itV = n[i]->GetRefVertexMapping().begin();
                    itV !=  n[i]->GetRefVertexMapping().end();
                    itV++)
                {
                    MG_SOMMET * v = *itV;
                    if (polycurve->ContainsRefVertex(v))
                    {
                        s[i] = polycurve->RefVertex_GetS(v);
                        if (n[i]->get_lien_topologie_reference() != v)
                        {
                            double xyzV[3]; v->get_point()->evaluer(xyzV);
                            n[i] = new MCNode (__mcEdge,v,xyzV[0],xyzV[1],xyzV[2]);
                            //n[i]->incrementer();
                        }
                        break;
                    }
                }
            }
            else
            {
                s[i] = polycurve->RefVertex_GetS(refVertex);
            }
        }
        else if (n[i]->GetRefEdgeMapping().size() == 1) // MC node is in a ref edge
        {
            polycurve->inverser(s[i], n[i]->get_coord(), 1E-6);
        }
    }
                        
    if (period != 0.0)
    {                    
        for (i=0;i<2;i++){
        if (s[i] > smax + fabs(smax-smin)*1E-6)
            s[i] -= period;
        if (s[i] < smin - fabs(smax-smin)*1E-6)
            s[i] += period; }
    }

    if (period != 0.0 && fabs(s[1] - s[0]) > .5*period)
        if (s[1]<s[0])
            s[1] += period;
        else s[1] += period;
    
    for (i=0; i<2; i++)
    {
        nodes[s[i]] = n[i];
    }

    {
        MG_ARETE * edge; double t,dt;
        double si = std::min(s[0],s[1]);
        polycurve->Parameter_SToRefEdgeT ( si+1E-6*fabs(s[1]-s[0]), &edge,&t,&dt,false);
        edges [ si ] = edge;
    }
           

    for (int i=0; i<2; i++)
    {
        OT_VECTEUR_3D v, v_node(n[i]->get_coord());
        polycurve->evaluer(s[i], v);
        double dist = (v_node-v).get_longueur();
        if (dist > .001*__mcEdge->GetPolyCurve()->get_longueur())
        {
            printf("Distance between projection and node = %f\n", dist);
            polycurve->inverser(s[i], n[i]->get_coord(), 1E-6);     
            polycurve->evaluer(s[i], v);
        }
    }

    // Create nodes in vertices between n1 and n2
    double si=0;
    for (i=0;si<std::max(s[1],s[0]);i++)
    {
        int iEdge = i%polycurve->GetRefEdgeCount();
        int iVertex = (period != 0.0)? i%(polycurve->GetRefVertexCount()-1): i%polycurve->GetRefVertexCount();
        if ((si-s[0])*(s[1]-si) > 0)
        {
            MG_SOMMET * vertex = polycurve->GetRefVertex(iVertex);
            double xyz[3];vertex->get_point()->evaluer(xyz);
            MCNode * ni = new MCNode ( __mcEdge, vertex, xyz[0], xyz[1], xyz[2] );
            nodes [ si ] = ni;
            edges [ si ] = polycurve->GetRefEdge(iEdge);
            //ni->incrementer();
        }
        si += polycurve->GetLength(iEdge);
    }

    std::map < double , MCNode * > refinedNodes = nodes;
    std::map < double , MG_ARETE * > refinedEdges = edges;

    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;
        double sii = itNii->first;
        double smid = .5*(si+sii);

        if (ni->get_lien_topologie_reference()->get_dimension() == 0)
        {
            MG_SOMMET * vi = (MG_SOMMET*)ni->get_lien_topologie_reference();
            if (nii->IsInVertex(vi))
                continue;
        }

        if (nii->get_lien_topologie_reference()->get_dimension() == 0)
        {
            MG_SOMMET * vii = (MG_SOMMET*)nii->get_lien_topologie_reference();
            if (ni->IsInVertex(vii))
                continue;
        }

        MG_ARETE * refEdge=0; double tMid; double dtMid;
        polycurve->Parameter_SToRefEdgeT(smid,&refEdge,&tMid,&dtMid,false);
        // and if the edge is periodic and tii-ti > period/2 then insert a node in the middle
        {
            if (refEdge->get_courbe()->est_periodique() != 0 )
            {
                MG_ARETE * refEdgei=0, * refEdgeii=0;
                double tii, dtii, ti, dti;
                polycurve->Parameter_SToRefEdgeT(si,&refEdgei,&ti,&dti,false);
                polycurve->Parameter_SToRefEdgeT(sii,&refEdgeii,&tii,&dtii,false);
                
                if (fabs(tii-ti) > .5*refEdge->get_courbe()->get_periode()) 
                {
                    double xyz[3];
                    refEdge->evaluer(tMid,xyz);
                    MCNode * newMCNode = new MCNode(__mcEdge, refEdge, tMid, xyz);
                    //newMCNode->incrementer();
                    refinedNodes[smid] = newMCNode;
                    refinedEdges[smid] = refEdge;
                }
            }
        }
    }
    nodes = refinedNodes;
    edges = refinedEdges;

    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;

        MG_ARETE * refEdge=edges[si];
        
        // Refine the mc segment discretization with discretization error in faces < segment length / 1000
        double deltaMax = get_longueur()*1E-2;
        if (deltaMax < 1E-7)
            deltaMax = 1E-7;
        RefineInRefEdge(__mcEdge, refEdge, ni, nii, deltaMax, 1E300, refinedNodes, refinedEdges );
    }
    nodes=refinedNodes;
    edges=refinedEdges;


    // detect edges which contain nodes, and make segments of this edge on these nodes
    for (std::map<double, MCNode * >::iterator itN = nodes.begin();
        itN != nodes.end();
        itN++)
    {
        std::map<double, MCNode * >::iterator itNii = itN;
        itNii++;
        if (itNii == nodes.end())
            break;
                           
        MCNode * ni = itN->second;
        MCNode * nii = itNii->second;
        double si = itN->first;
        std::map<double, MG_ARETE*>::iterator itRefEdge = edges.find(si);
        if (itRefEdge != edges.end())
            edge = itRefEdge->second;
        else
            edge = 0;
        if (ni->IsInEdge(edge) == false ||
            nii->IsInEdge(edge) == false )
        {
            for (MCNode::EMapIterator itNiRefEdge = ni->GetRefEdgeMapping().begin();
                itNiRefEdge != ni->GetRefEdgeMapping().end();
                itNiRefEdge++)
            {
                edge = itNiRefEdge->first;
                if ( nii->IsInEdge(edge) )
                {
                    break;
                }
            }
        }

        double params[2];
        params[0]=itN->first;
        params[1]=itNii->first;

        for (int i=0;i<2;i++)
        if ((s[1]-params[i])*(s[0]-params[i])>0)
            printf("pb\n");

        std::vector < MCNode * > segEdge;
        segEdge.push_back (ni);
        segEdge.push_back (nii);

        E.insert(std::make_pair(edge, segEdge));
        // copy the curve in adjacent faces
        for (int itF = 0; itF < edge->get_nb_mg_coarete(); itF++)
        {
            MG_FACE * face = edge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
            F.insert(std::make_pair(face, segEdge));
        }      

        if (CurveNodes.size () == 0 || CurveNodes[CurveNodes.size()-1] != ni)
            CurveNodes.push_back(ni);
        CurveNodes.push_back(nii);
        CurveTopo.push_back(edge);
    }

    if (_polylineEvaluator == NULL) _polylineEvaluator = new MCNodePolyline(liaison_topologique,CurveNodes,CurveTopo);
    else
            for (unsigned i=0;i<CurveNodes.size();i++)
                _polylineEvaluator->Add(CurveNodes[i],CurveTopo[i]);
                
    return 1;
}

MCFace  * MCSegment::Find_MCFace(MCEdge * __mcEdge, MG_SOMMET * __refVertex1, MG_SOMMET * __refVertex2)
{
    for (int itF = 0; itF != __mcEdge->get_nb_mg_coarete(); itF++)
    {
        MCFace * F = (MCFace*)__mcEdge->get_mg_coarete(itF)->get_boucle()->get_mg_face();
        bool v1Found = false, v2Found = false;
        for (int itFLoop = 0; itFLoop != F->get_nb_mg_boucle(); itFLoop ++)
        {
            MG_BOUCLE * loop = F->get_mg_boucle(itFLoop);
            for (int itFLoopEdge = 0; itFLoopEdge != loop->get_nb_mg_coarete(); itFLoopEdge++)
            {
                MCEdge * mcEdge = (MCEdge *)loop->get_mg_coarete( itFLoopEdge )->get_arete();
                if ( mcEdge->GetPolyCurve()->ContainsRefVertex(__refVertex1) )
                    v1Found = true;                                
                if ( mcEdge->GetPolyCurve()->ContainsRefVertex(__refVertex2) )
                    v2Found = true;   
                if (v1Found && v2Found)
                    return F;
            }
        }
    }
    return NULL;
}

int MCSegment::Construct_MergedVertices(MCNode * __mcNode)
{
    if (get_lien_topologie()->get_dimension() != 1)
    {
        printf("Warning: MCSegment topology is not a MCEdge ! then it's not possible to construct a mapping on a merged vertex\n");
        return 0;
    }
    MCEdge * mcEdge = (MCEdge*)get_lien_topologie();
    MG_SOMMET * polycurveRefVertex = NULL;
    MG_SOMMET * mcVertexRefVertex = (MG_SOMMET*)__mcNode->get_lien_topologie();
    MCFace * mcFace;
    for (std::set<MG_SOMMET*>::iterator itV = __mcNode->GetRefVertexMapping().begin();
        itV != __mcNode->GetRefVertexMapping().end();
        itV++ )
    {
        MG_SOMMET * v = *itV;
        if (v != mcVertexRefVertex)
        {
            polycurveRefVertex = v;
            mcFace = Find_MCFace(mcEdge, polycurveRefVertex, mcVertexRefVertex);
            if (mcFace)
            {
                MCNode * polycurve_mcNode = NULL;
                double v_xyz[3];
                v->get_point()->evaluer(v_xyz);

                polycurve_mcNode = new MCNode (get_lien_topologie(), v, v_xyz[0],v_xyz[1],v_xyz[2]);

                polycurve_mcNode->incrementer();
                if (polycurve_mcNode)
                    Construct_MCFaceByShortestPath(mcFace, __mcNode, polycurve_mcNode);
                    
                break;
            }
        }
    }

    return 1;
}

BOITE_3D MCSegment::get_boite_3D(void)
{
    return BOITE_3D(_bbox[0],_bbox[1],_bbox[2],_bbox[3],_bbox[4],_bbox[5]);
}

void MCSegment::SetSaveFormat(char __format)
{
    _saveFormat=__format;
    if (__format == 2)
    {
        MG_ELEMENT_TOPOLOGIQUE * refTopo = 0;
        switch (liaison_topologique->get_dimension())
        {
        case 1:
        {
            refTopo = ((MCEdge*)liaison_topologique)->GetPolyCurve()->GetRefEdge(0);
            break;
        }
        case 2:
        {
            refTopo = ((MCFace*)liaison_topologique)->GetPolySurface()->GetRefFace(0);
            break;
        }
        }
        change_lien_topologie(refTopo);
    }
}

void MCSegment::enregistrer(std::ostream& o,double version)
{
    if (_saveFormat == 0) // Save as CAD4FE_MC_SEGMENT in magic file
    {
        if(liaison_topologique->get_dimension()==1) o << "%" << get_id() << "=CAD4FE_MCSEGMENT($"<< get_lien_topologie()->get_id() << ",$" << noeud1->get_id() << ",$" << noeud2->get_id() <<  ");" << std::endl;
    }
    else if (_saveFormat == 1) // Save as MG_SEGMENT in magic file
    {
        MG_SEGMENT::enregistrer(o,version);
    }
    else if (_saveFormat == 2)
    {                      
        MG_SEGMENT::enregistrer(o,version);
    }
}                       
int MCSegment::get_type_entite()
{
    return IDMCSEGMENT;
}


Revision:	763
Committed:	Wed Dec 2 19:55:53 2015 UTC (9 years, 5 months ago) by francois
File size:	38557 byte(s)
Log Message:	Le fichier MAGiC est maintenant versionné. LA version actuelle est 2.0. L'ancienne version est 1.0. Tout est transparent pour l'utilisateur. Les vieilles versions sont lisibles mais les nouveaux enregistrements sont dans la version la plus récente. Changement des conditions aux limites : ajout d'un parametre pour dire si la condition numerique est une valeur ou une formule ou un lien vers une autre entité magic. Les parametres pour saisir sont maintenant -ccf -ccfi -ccff -ccft -ccfit -ccfft