ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_MCTriangle.cpp
Revision: 791
Committed: Fri Mar 18 21:23:00 2016 UTC (9 years, 1 month ago) by francois
File size: 12845 byte(s)
Log Message:
Transferts de la definition des origines dans le namespace MAGiC

File Contents

# 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 791 : MG_TRIANGLE(num,topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::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 791 : MG_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::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(j*D/(1-i*D),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     }