ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_MCTriangle.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (11 months ago) by francois
File size: 13699 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

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