CAD4FE/src/CAD4FE_MCTriangle.cpp

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


#pragma hdrstop
                                 
#include "gestionversion.h"
                                
#include "CAD4FE_MCFace.h" 
#include "CAD4FE_PolySurface.h"
#include "CAD4FE_MCNode.h"
#include "CAD4FE_MCSegment.h"
                                
#include "CAD4FE_MCTriangle.h"
#include "CAD4FE_Geometric_Tools.h"
#include "ot_quadrature_gauss.h"

#include <math.h>

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

#pragma package(smart_init)

using namespace CAD4FE;

MCTriangle::MCTriangle(unsigned long num,MG_ELEMENT_TOPOLOGIQUE* topo,MCNode *mgnoeud1,MCNode *mgnoeud2,MCNode *mgnoeud3,MCSegment* mgsegment1,MCSegment* mgsegment2,MCSegment* mgsegment3 )
: MG_TRIANGLE(num,topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAILLEUR_AUTO), _saveFormat(0)
{
    _indexNodeA = ComputeIndexNodeA();
}

MCTriangle::MCTriangle(MG_ELEMENT_TOPOLOGIQUE* topo,MCNode *mgnoeud1,MCNode *mgnoeud2,MCNode *mgnoeud3,MCSegment* mgsegment1,MCSegment* mgsegment2,MCSegment* mgsegment3)
: MG_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAILLEUR_AUTO), _saveFormat(0)
{
    _indexNodeA = ComputeIndexNodeA();
}

MCTriangle::MCTriangle(MCTriangle& mdd)
: MG_TRIANGLE(mdd), _saveFormat(0)
{
    _indexNodeA = mdd._indexNodeA;
}

MCTriangle::~MCTriangle()
{

}

int MCTriangle::ComputeIndexNodeA()
{
    MCSegment * seg[3];
    seg[0]=(MCSegment *)get_segment1();
    seg[1]=(MCSegment *)get_segment2();
    seg[2]=(MCSegment *)get_segment3();
    double quality[3];        
    int indexSegBestQuality=0;
    for (int i=0; i<3; i++)
    {
        quality[i] = seg[i]->get_longueur()/seg[i]->get_longueur_geo();
        if (quality[i] > quality[indexSegBestQuality])
            indexSegBestQuality=i;
    }
    int indexNodeA = (indexSegBestQuality+2)%3;
    return indexNodeA;
}
                       
MCNode * MCTriangle::GetRefNode(int __index)
{
    int realIndex = (_indexNodeA + __index) % 3;
    switch (realIndex)
    {
        case 0:
            return (MCNode*)get_noeud1();
        case 1:
            return (MCNode*)get_noeud2();
        case 2:
            return (MCNode*)get_noeud3();
    }
    return NULL;
}

MCSegment * MCTriangle::GetRefSegment(int __index)
{
    int realIndex = (_indexNodeA + __index) % 3;
    switch (realIndex)
    {
        case 0:
            return (MCSegment*)get_segment1();
        case 1:
            return (MCSegment*)get_segment2();
        case 2:
            return (MCSegment*)get_segment3();
    }
    return NULL;
}
                  
void MCTriangle::SetSaveFormat(char __format)
{
    _saveFormat=__format;
    if (__format == 2)
    {
        MG_ELEMENT_TOPOLOGIQUE * refTopo = 0;
        refTopo = ((MCFace*)liaison_topologique)->GetPolySurface()->GetRefFace(0);
        change_lien_topologie(refTopo);
    }
}

void MCTriangle::enregistrer(std::ostream& o,double version)
{
    if (_saveFormat == 0)
    {
          o << "%" << get_id() << "=CAD4FE_MCTRIANGLE($"<< get_lien_topologie()->get_id() << ",$" << noeud1->get_id() << ",$" << noeud2->get_id() << ",$" << noeud3->get_id()  << ");" << std::endl;
    }
    else if (_saveFormat == 1)
    {
        MG_TRIANGLE::enregistrer(o,version);
    }
    else if (_saveFormat == 2)
    {
        MG_TRIANGLE::enregistrer(o,version);
    }
}


MCSegment * MCTriangle::evaluer_geo_isoparam_u(double __u)
{
    MCSegment * B2C2 = NULL;
    MCNode * A = GetRefNode(0);
    MCNode * B = GetRefNode(1);
    MCSegment * BA = GetRefSegment(0);
    MCSegment * CA = GetRefSegment(2);
    MCNode * B2 = NULL;
    MCNode * C2 = NULL;

    if (__u == 0)
    {
        B2C2 = new MCSegment(*GetRefSegment(1));
        if (B2C2->get_noeud1() != B)
            B2C2->change_orientation(-1);
        else   
            B2C2->change_orientation(+1);
    }
    else
    {
        double tBA,tCA;
        if (BA->get_noeud2() != A)
            tBA = 1 - __u;
        else
            tBA = __u;

        if (CA->get_noeud2() != A)
            tCA = 1 - __u;
        else
            tCA = __u;

        B2 = new MCNode;
        C2 = new MCNode;

        BA->evaluer_geo(tBA, B2);
        CA->evaluer_geo(tCA, C2);

        MCFace * mcFace = (MCFace*) liaison_topologique;
        int nbB2RefFaceInMCFace=0;
        for (MCNode::FMapIterator itF = B2->GetRefFaceMapping().begin();
            itF != B2->GetRefFaceMapping().end();
            itF++)
        {
            MG_FACE * face = itF->first;
            if (mcFace->GetPolySurface()->Contains(face))
                nbB2RefFaceInMCFace++;
        }                       
        int nbC2RefFaceInMCFace=0;
        for (MCNode::FMapIterator itF = C2->GetRefFaceMapping().begin();
            itF != C2->GetRefFaceMapping().end();
            itF++)
        {
            MG_FACE * face = itF->first;
            if (mcFace->GetPolySurface()->Contains(face))
                nbC2RefFaceInMCFace++;
        }
        if (nbB2RefFaceInMCFace == 0 || nbC2RefFaceInMCFace == 0)
        {// nodes of the segment are outside the topology
         // E.g. along a segment connecting a ref vertex merged in a MC vertex
            B2C2 = new MCSegment( liaison_topologique,B2,C2,0,
                MCSegment::CstrNoRefTopologyMapping
            );
        }
        else
        {
            B2C2 = new MCSegment(liaison_topologique,B2,C2,0,
                MCSegment::CstrMCFaceByPlaneIntr
            //MCSegment::CstrMCFaceByShortestPath
            );
        }
    }
    
    return B2C2;
}

void MCTriangle::evaluer_geo(double __uv[2], MCNode * __result, double * tangent_u, double * tangent_v)
{
    MCSegment * B2C2 = 0, *B3C3=0;
    double du=1E-3;
    evaluer_geo(&B2C2,&B3C3,du,__uv,__result,tangent_u,tangent_v);
    delete B2C2;
    delete B3C3;
}
                
void MCTriangle::evaluer_geo(MCSegment ** B2C2, MCSegment ** B3C3, double __du, double __uv[2], MCNode * __P, double * tangent_u, double * tangent_v)
{
    double tAI;
    double du = (1-__uv[0]-__uv[1] > __du) ? __du : -__du;
    if (fabs(1-__uv[0]) < __du)
    {
        *__P = *GetRefNode(0);
        du = -fabs(__du);
    }
    // v <= u and u in 0,1 ... v(segment in [0,1]) v_seg = v/(1-u)
    if ((*B2C2) == NULL) (*B2C2) = evaluer_geo_isoparam_u(__uv[0]);
    (*B2C2)->evaluer_geo(__uv[1]/(1-__uv[0]), __P, tangent_v);
    if (tangent_u)
    {
        if ((*B3C3) == NULL)
            (*B3C3) = evaluer_geo_isoparam_u(__uv[0]+du);
        MCNode * rdu = new MCNode;
        (*B3C3)->evaluer_geo(__uv[1]/(1-(__uv[0]+du)), rdu);
        for (int i=0;i<3;i++)
            tangent_u[i] = (rdu->get_coord()[i]-__P->get_coord()[i])/du;
        delete rdu;
    }
    if (tangent_v)
    {
        for (int i=0;i<3;i++) tangent_v[i] *= 1/(1-__uv[0]);
    }
}

void MCTriangle::evaluer_euc(double __uv[2], double * __result, double * tangent_u, double * tangent_v)
{                             
    MCNode * A = GetRefNode(0);
    MCNode * B = GetRefNode(1);
    MCNode * C = GetRefNode(2);

    double tBA=__uv[0], tCA=__uv[0];
    double xyzB2[3], xyzC2[3];
    for (int i=0; i<3; i++)
    {
        xyzB2[i] = tBA*(A->get_coord()[i]-B->get_coord()[i])+B->get_coord()[i];
        xyzC2[i] = tCA*(A->get_coord()[i]-C->get_coord()[i])+C->get_coord()[i];
    }

    double t=__uv[1]/(1-__uv[0]);
                         
    for (int i=0; i<3; i++)
        __result[i] = xyzB2[i]+(xyzC2[i]-xyzB2[i])*t;
                          
    if (tangent_u)
    {
        for (int i=0;i<3;i++) tangent_u[i] = (A->get_coord()[i]-B->get_coord()[i]);
    }

    if (tangent_v)
    {
        for (int i=0;i<3;i++) tangent_v[i] = (C->get_coord()[i]-B->get_coord()[i]);
    }
}

/*
// Tessellation with Quadrature points
void MCTriangle::Tessellate(unsigned __nbPtsPerSegment, std::vector <MCNode *> & __nodes, std::vector <unsigned> & __indexedTriangleSet)
{                           
    double I=0;
    int nx=__nbPtsPerSegment;
    int ny=__nbPtsPerSegment;

    double * x = new double[nx+1];
    double * wx = new double[nx+1];
    double * y = new double[ny+1];
    double * wy = new double[ny+1];

    OT_QUADRATURE_GAUSS::gauss_legendre_points(-1, +1, x, wx, nx);
    if (nx == ny)
        for (int i=0;i<=nx;i++) { y[i]=x[i]; wy[i]=wx[i]; }
    else  OT_QUADRATURE_GAUSS::gauss_legendre_points(-1, +1, y, wy, ny);

    std::map<double, MCSegment*> mapIso;
    std::map<double, MCSegment*>::iterator itMapIso;

    MCNode * trianglePts = new MCNode [nx*ny];
    OT_VECTEUR_3D * arrayDerU = new OT_VECTEUR_3D [nx*ny]; 
    OT_VECTEUR_3D * arrayDerV = new OT_VECTEUR_3D [nx*ny];  
    OT_VECTEUR_3D * pts = new OT_VECTEUR_3D [nx*ny];
#define idx(i,j) ny*(i-1)+(j-1)
#define node(i,j) (trianglePts+idx(i,j))
#define pts(i,j) pts[idx(i,j)]
#define derU(i,j) arrayDerU[idx(i,j)]     
#define derU2(i,j) arrayDerU2[idx(i,j)]
#define derV(i,j) arrayDerV[idx(i,j)]
    for (int i=1;i<=nx;i++)
    {
        double U=(1+x[i])*.5;
        double UU=(i<nx)?(1+x[i+1])*.5:1;
        double DU=(UU-U);
        MCSegment * isoU = NULL, * isoUU = NULL;
        itMapIso = mapIso.find(U);
        if (itMapIso!=mapIso.end())
            isoU = itMapIso->second;
        itMapIso = mapIso.find(UU);
        if (itMapIso!=mapIso.end())
            isoUU = itMapIso->second;
        for (int j=1;j<=ny;j++)
        {
            double V=.25*((1-x[i])*(1+y[j]));
            double uv[2]={U,V};
            evaluer_geo(&isoU, &isoUU, DU, uv, node(i,j), derU(i,j), derV(i,j)); 
            __nodes.push_back(node(i,j));
            pts(i,j) = OT_VECTEUR_3D (node(i,j)->get_coord());
         }
         mapIso[U] = isoU;
         mapIso[UU] = isoUU;
    }   


        for (int i=0;i+1<nx;i++)
        for (int j=0;j+1<ny;j++)
        {
            unsigned p = ny*i+j;
            unsigned north = p+ny;
            unsigned east = p+1;
            unsigned northeast = north+1;
            // triangle 1
            __indexedTriangleSet.push_back(p);
            __indexedTriangleSet.push_back(east);
            __indexedTriangleSet.push_back(northeast);
            // triangle 2
            __indexedTriangleSet.push_back(p);
            __indexedTriangleSet.push_back(northeast);
            __indexedTriangleSet.push_back(north);
        }

    for (itMapIso = mapIso.begin();itMapIso != mapIso.end(); itMapIso++)
        {
            MCSegment * iso=itMapIso->second;
            if (iso->get_nb_reference() == 0) delete iso;
        }
                                          
//    delete [] trianglePts;
    delete [] arrayDerU;  
    delete [] arrayDerV;
    delete [] pts;
#undef idx(i,j)
#undef node(i,j)
#undef pts(i,j)
#undef derU(i,j)
#undef derV(i,j)
}   */

void MCTriangle::Tessellate(unsigned __nbPtsPerSegment, std::vector <MCNode *> & __nodes, std::vector <unsigned> & __indexedTriangleSet)
{
    unsigned N = __nbPtsPerSegment;
    double D = 1.0f/(N-1);
    std::map <unsigned, unsigned> mapIndices;

        for (unsigned i=0;i<N; i++)
        {
            if (i != N-1)
            {
                MCSegment * isoparametric = evaluer_geo_isoparam_u(i*D);
                for (unsigned j=0;j<N-i; j++)
                {
                    mapIndices[i*N+j] = __nodes.size();
                    MCNode * mcNode = new MCNode;
                    isoparametric->evaluer_geo(j*D/(1-i*D),mcNode);
                    __nodes.push_back(mcNode);
                    mcNode->incrementer();
                }
                if (isoparametric->get_nb_reference() == 0)
                    delete isoparametric;
            }
            else
            {                 
                mapIndices[i*N] = __nodes.size();
                MCNode * mcNode = GetRefNode(0);
                mcNode->incrementer();      
                __nodes.push_back(mcNode);
            }
        }
        for (unsigned i=0;i+1<N;i++)
        for (unsigned j=0;j+1<N-i;j++)
        {
            unsigned p = i*N+j;
            unsigned north = p+N;
            unsigned east = p+1; 
            unsigned northeast = p+N+1;

            unsigned t[2][3]={{p,east,north},{east,northeast,north}};

            for (int k=0; k<2; k++)
            {
                if (k == 1 && j+2==N-i)
                    continue;
                for (int l=0; l<3; l++)
                    __indexedTriangleSet.push_back(mapIndices[t[k][l]]);
            }
        }
}
Revision:	763
Committed:	Wed Dec 2 19:55:53 2015 UTC (9 years, 5 months ago) by francois
File size:	12813 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
#	User	Rev	Content
1	foucault	27	//---------------------------------------------------------------------------
2
3
4			#pragma hdrstop
5
6			#include "gestionversion.h"
7
8			#include "CAD4FE_MCFace.h"
9			#include "CAD4FE_PolySurface.h"
10			#include "CAD4FE_MCNode.h"
11			#include "CAD4FE_MCSegment.h"
12
13			#include "CAD4FE_MCTriangle.h"
14	foucault	569	#include "CAD4FE_Geometric_Tools.h"
15	foucault	27	#include "ot_quadrature_gauss.h"
16
17			#include <math.h>
18
19			//---------------------------------------------------------------------------
20
21			#pragma package(smart_init)
22
23			using namespace CAD4FE;
24
25			MCTriangle::MCTriangle(unsigned long num,MG_ELEMENT_TOPOLOGIQUE* topo,MCNode mgnoeud1,MCNode mgnoeud2,MCNode mgnoeud3,MCSegment mgsegment1,MCSegment* mgsegment2,MCSegment* mgsegment3 )
26	francois	35	: MG_TRIANGLE(num,topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAILLEUR_AUTO), _saveFormat(0)
27	foucault	27	{
28			_indexNodeA = ComputeIndexNodeA();
29			}
30
31			MCTriangle::MCTriangle(MG_ELEMENT_TOPOLOGIQUE* topo,MCNode mgnoeud1,MCNode mgnoeud2,MCNode mgnoeud3,MCSegment mgsegment1,MCSegment* mgsegment2,MCSegment* mgsegment3)
32	francois	35	: MG_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAILLEUR_AUTO), _saveFormat(0)
33	foucault	27	{
34			_indexNodeA = ComputeIndexNodeA();
35			}
36
37			MCTriangle::MCTriangle(MCTriangle& mdd)
38			: MG_TRIANGLE(mdd), _saveFormat(0)
39			{
40			_indexNodeA = mdd._indexNodeA;
41			}
42
43			MCTriangle::~MCTriangle()
44			{
45
46			}
47
48			int MCTriangle::ComputeIndexNodeA()
49			{
50			MCSegment * seg[3];
51			seg[0]=(MCSegment *)get_segment1();
52			seg[1]=(MCSegment *)get_segment2();
53			seg[2]=(MCSegment *)get_segment3();
54			double quality[3];
55			int indexSegBestQuality=0;
56			for (int i=0; i<3; i++)
57			{
58			quality[i] = seg[i]->get_longueur()/seg[i]->get_longueur_geo();
59			if (quality[i] > quality[indexSegBestQuality])
60			indexSegBestQuality=i;
61			}
62			int indexNodeA = (indexSegBestQuality+2)%3;
63			return indexNodeA;
64			}
65
66			MCNode * MCTriangle::GetRefNode(int __index)
67			{
68			int realIndex = (_indexNodeA + __index) % 3;
69			switch (realIndex)
70			{
71			case 0:
72			return (MCNode*)get_noeud1();
73			case 1:
74			return (MCNode*)get_noeud2();
75			case 2:
76			return (MCNode*)get_noeud3();
77			}
78			return NULL;
79			}
80
81			MCSegment * MCTriangle::GetRefSegment(int __index)
82			{
83			int realIndex = (_indexNodeA + __index) % 3;
84			switch (realIndex)
85			{
86			case 0:
87			return (MCSegment*)get_segment1();
88			case 1:
89			return (MCSegment*)get_segment2();
90			case 2:
91			return (MCSegment*)get_segment3();
92			}
93			return NULL;
94			}
95
96			void MCTriangle::SetSaveFormat(char __format)
97			{
98			_saveFormat=__format;
99			if (__format == 2)
100			{
101			MG_ELEMENT_TOPOLOGIQUE * refTopo = 0;
102			refTopo = ((MCFace*)liaison_topologique)->GetPolySurface()->GetRefFace(0);
103			change_lien_topologie(refTopo);
104			}
105			}
106
107	francois	763	void MCTriangle::enregistrer(std::ostream& o,double version)
108	foucault	27	{
109			if (_saveFormat == 0)
110			{
111			o << "%" << get_id() << "=CAD4FE_MCTRIANGLE($"<< get_lien_topologie()->get_id() << ",$" << noeud1->get_id() << ",$" << noeud2->get_id() << ",$" << noeud3->get_id() << ");" << std::endl;
112			}
113			else if (_saveFormat == 1)
114			{
115	francois	763	MG_TRIANGLE::enregistrer(o,version);
116	foucault	27	}
117			else if (_saveFormat == 2)
118			{
119	francois	763	MG_TRIANGLE::enregistrer(o,version);
120	foucault	27	}
121			}
122
123
124			MCSegment * MCTriangle::evaluer_geo_isoparam_u(double __u)
125			{
126			MCSegment * B2C2 = NULL;
127			MCNode * A = GetRefNode(0);
128			MCNode * B = GetRefNode(1);
129			MCSegment * BA = GetRefSegment(0);
130			MCSegment * CA = GetRefSegment(2);
131			MCNode * B2 = NULL;
132			MCNode * C2 = NULL;
133
134			if (__u == 0)
135			{
136			B2C2 = new MCSegment(*GetRefSegment(1));
137			if (B2C2->get_noeud1() != B)
138			B2C2->change_orientation(-1);
139			else
140			B2C2->change_orientation(+1);
141			}
142			else
143			{
144			double tBA,tCA;
145			if (BA->get_noeud2() != A)
146			tBA = 1 - __u;
147			else
148			tBA = __u;
149
150			if (CA->get_noeud2() != A)
151			tCA = 1 - __u;
152			else
153			tCA = __u;
154
155			B2 = new MCNode;
156			C2 = new MCNode;
157
158			BA->evaluer_geo(tBA, B2);
159			CA->evaluer_geo(tCA, C2);
160
161			MCFace * mcFace = (MCFace*) liaison_topologique;
162			int nbB2RefFaceInMCFace=0;
163			for (MCNode::FMapIterator itF = B2->GetRefFaceMapping().begin();
164			itF != B2->GetRefFaceMapping().end();
165			itF++)
166			{
167			MG_FACE * face = itF->first;
168			if (mcFace->GetPolySurface()->Contains(face))
169			nbB2RefFaceInMCFace++;
170			}
171			int nbC2RefFaceInMCFace=0;
172			for (MCNode::FMapIterator itF = C2->GetRefFaceMapping().begin();
173			itF != C2->GetRefFaceMapping().end();
174			itF++)
175			{
176			MG_FACE * face = itF->first;
177			if (mcFace->GetPolySurface()->Contains(face))
178			nbC2RefFaceInMCFace++;
179			}
180			if (nbB2RefFaceInMCFace == 0 \|\| nbC2RefFaceInMCFace == 0)
181			{// nodes of the segment are outside the topology
182			// E.g. along a segment connecting a ref vertex merged in a MC vertex
183			B2C2 = new MCSegment( liaison_topologique,B2,C2,0,
184			MCSegment::CstrNoRefTopologyMapping
185			);
186			}
187			else
188			{
189			B2C2 = new MCSegment(liaison_topologique,B2,C2,0,
190			MCSegment::CstrMCFaceByPlaneIntr
191			//MCSegment::CstrMCFaceByShortestPath
192			);
193			}
194			}
195
196			return B2C2;
197			}
198
199			void MCTriangle::evaluer_geo(double __uv[2], MCNode * __result, double * tangent_u, double * tangent_v)
200			{
201			MCSegment * B2C2 = 0, *B3C3=0;
202			double du=1E-3;
203			evaluer_geo(&B2C2,&B3C3,du,__uv,__result,tangent_u,tangent_v);
204			delete B2C2;
205			delete B3C3;
206			}
207
208			void MCTriangle::evaluer_geo(MCSegment B2C2, MCSegment B3C3, double __du, double __uv[2], MCNode * __P, double * tangent_u, double * tangent_v)
209			{
210			double tAI;
211			double du = (1-__uv[0]-__uv[1] > __du) ? __du : -__du;
212			if (fabs(1-__uv[0]) < __du)
213			{
214			__P = GetRefNode(0);
215			du = -fabs(__du);
216			}
217			// v <= u and u in 0,1 ... v(segment in [0,1]) v_seg = v/(1-u)
218			if ((B2C2) == NULL) (B2C2) = evaluer_geo_isoparam_u(__uv[0]);
219			(*B2C2)->evaluer_geo(__uv[1]/(1-__uv[0]), __P, tangent_v);
220			if (tangent_u)
221			{
222			if ((*B3C3) == NULL)
223			(*B3C3) = evaluer_geo_isoparam_u(__uv[0]+du);
224			MCNode * rdu = new MCNode;
225			(*B3C3)->evaluer_geo(__uv[1]/(1-(__uv[0]+du)), rdu);
226			for (int i=0;i<3;i++)
227			tangent_u[i] = (rdu->get_coord()[i]-__P->get_coord()[i])/du;
228			delete rdu;
229			}
230			if (tangent_v)
231			{
232			for (int i=0;i<3;i++) tangent_v[i] *= 1/(1-__uv[0]);
233			}
234			}
235
236			void MCTriangle::evaluer_euc(double __uv[2], double * __result, double * tangent_u, double * tangent_v)
237			{
238			MCNode * A = GetRefNode(0);
239			MCNode * B = GetRefNode(1);
240			MCNode * C = GetRefNode(2);
241
242			double tBA=__uv[0], tCA=__uv[0];
243			double xyzB2[3], xyzC2[3];
244			for (int i=0; i<3; i++)
245			{
246			xyzB2[i] = tBA*(A->get_coord()[i]-B->get_coord()[i])+B->get_coord()[i];
247			xyzC2[i] = tCA*(A->get_coord()[i]-C->get_coord()[i])+C->get_coord()[i];
248			}
249
250			double t=__uv[1]/(1-__uv[0]);
251
252			for (int i=0; i<3; i++)
253			__result[i] = xyzB2[i]+(xyzC2[i]-xyzB2[i])*t;
254
255			if (tangent_u)
256			{
257			for (int i=0;i<3;i++) tangent_u[i] = (A->get_coord()[i]-B->get_coord()[i]);
258			}
259
260			if (tangent_v)
261			{
262			for (int i=0;i<3;i++) tangent_v[i] = (C->get_coord()[i]-B->get_coord()[i]);
263			}
264			}
265
266			/*
267			// Tessellation with Quadrature points
268			void MCTriangle::Tessellate(unsigned __nbPtsPerSegment, std::vector <MCNode *> & __nodes, std::vector <unsigned> & __indexedTriangleSet)
269			{
270			double I=0;
271			int nx=__nbPtsPerSegment;
272			int ny=__nbPtsPerSegment;
273
274			double * x = new double[nx+1];
275			double * wx = new double[nx+1];
276			double * y = new double[ny+1];
277			double * wy = new double[ny+1];
278
279			OT_QUADRATURE_GAUSS::gauss_legendre_points(-1, +1, x, wx, nx);
280			if (nx == ny)
281			for (int i=0;i<=nx;i++) { y[i]=x[i]; wy[i]=wx[i]; }
282			else OT_QUADRATURE_GAUSS::gauss_legendre_points(-1, +1, y, wy, ny);
283
284			std::map<double, MCSegment*> mapIso;
285			std::map<double, MCSegment*>::iterator itMapIso;
286
287			MCNode * trianglePts = new MCNode [nx*ny];
288			OT_VECTEUR_3D * arrayDerU = new OT_VECTEUR_3D [nx*ny];
289			OT_VECTEUR_3D * arrayDerV = new OT_VECTEUR_3D [nx*ny];
290			OT_VECTEUR_3D * pts = new OT_VECTEUR_3D [nx*ny];
291			#define idx(i,j) ny*(i-1)+(j-1)
292			#define node(i,j) (trianglePts+idx(i,j))
293			#define pts(i,j) pts[idx(i,j)]
294			#define derU(i,j) arrayDerU[idx(i,j)]
295			#define derU2(i,j) arrayDerU2[idx(i,j)]
296			#define derV(i,j) arrayDerV[idx(i,j)]
297			for (int i=1;i<=nx;i++)
298			{
299			double U=(1+x[i])*.5;
300			double UU=(i<nx)?(1+x[i+1])*.5:1;
301			double DU=(UU-U);
302			MCSegment * isoU = NULL, * isoUU = NULL;
303			itMapIso = mapIso.find(U);
304			if (itMapIso!=mapIso.end())
305			isoU = itMapIso->second;
306			itMapIso = mapIso.find(UU);
307			if (itMapIso!=mapIso.end())
308			isoUU = itMapIso->second;
309			for (int j=1;j<=ny;j++)
310			{
311			double V=.25((1-x[i])(1+y[j]));
312			double uv[2]={U,V};
313			evaluer_geo(&isoU, &isoUU, DU, uv, node(i,j), derU(i,j), derV(i,j));
314			__nodes.push_back(node(i,j));
315			pts(i,j) = OT_VECTEUR_3D (node(i,j)->get_coord());
316			}
317			mapIso[U] = isoU;
318			mapIso[UU] = isoUU;
319			}
320
321
322			for (int i=0;i+1<nx;i++)
323			for (int j=0;j+1<ny;j++)
324			{
325			unsigned p = ny*i+j;
326			unsigned north = p+ny;
327			unsigned east = p+1;
328			unsigned northeast = north+1;
329			// triangle 1
330			__indexedTriangleSet.push_back(p);
331			__indexedTriangleSet.push_back(east);
332			__indexedTriangleSet.push_back(northeast);
333			// triangle 2
334			__indexedTriangleSet.push_back(p);
335			__indexedTriangleSet.push_back(northeast);
336			__indexedTriangleSet.push_back(north);
337			}
338
339			for (itMapIso = mapIso.begin();itMapIso != mapIso.end(); itMapIso++)
340			{
341			MCSegment * iso=itMapIso->second;
342			if (iso->get_nb_reference() == 0) delete iso;
343			}
344
345			// delete [] trianglePts;
346			delete [] arrayDerU;
347			delete [] arrayDerV;
348			delete [] pts;
349			#undef idx(i,j)
350			#undef node(i,j)
351			#undef pts(i,j)
352			#undef derU(i,j)
353			#undef derV(i,j)
354			} */
355
356			void MCTriangle::Tessellate(unsigned __nbPtsPerSegment, std::vector <MCNode *> & __nodes, std::vector <unsigned> & __indexedTriangleSet)
357			{
358			unsigned N = __nbPtsPerSegment;
359			double D = 1.0f/(N-1);
360			std::map <unsigned, unsigned> mapIndices;
361
362	francois	763	for (unsigned i=0;i<N; i++)
363			{
364			if (i != N-1)
365			{
366			MCSegment * isoparametric = evaluer_geo_isoparam_u(i*D);
367			for (unsigned j=0;j<N-i; j++)
368			{
369			mapIndices[i*N+j] = __nodes.size();
370			MCNode * mcNode = new MCNode;
371			isoparametric->evaluer_geo(jD/(1-iD),mcNode);
372			__nodes.push_back(mcNode);
373			mcNode->incrementer();
374			}
375			if (isoparametric->get_nb_reference() == 0)
376			delete isoparametric;
377			}
378			else
379			{
380			mapIndices[i*N] = __nodes.size();
381			MCNode * mcNode = GetRefNode(0);
382			mcNode->incrementer();
383			__nodes.push_back(mcNode);
384			}
385	foucault	27	}
386			for (unsigned i=0;i+1<N;i++)
387			for (unsigned j=0;j+1<N-i;j++)
388			{
389			unsigned p = i*N+j;
390			unsigned north = p+N;
391			unsigned east = p+1;
392			unsigned northeast = p+N+1;
393
394			unsigned t[2][3]={{p,east,north},{east,northeast,north}};
395
396			for (int k=0; k<2; k++)
397			{
398			if (k == 1 && j+2==N-i)
399			continue;
400			for (int l=0; l<3; l++)
401			__indexedTriangleSet.push_back(mapIndices[t[k][l]]);
402			}
403			}
404			}