MAGiC  V5.0
Mailleurs Automatiques de Géometries intégrés à la Cao
CAD4FE_MCTriangle.cpp
Aller à la documentation de ce fichier.
1 //####//------------------------------------------------------------
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 
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 #include "CAD4FE_Geometric_Tools.h"
35 #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 : MG_TRIANGLE(num,topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::MAILLEUR_AUTO), _saveFormat(0)
46 {
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 : MG_TRIANGLE(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,MAGIC::ORIGINE::MAILLEUR_AUTO), _saveFormat(0)
52 {
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 void MCTriangle::enregistrer(std::ostream& o,double version)
127 {
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  MG_TRIANGLE::enregistrer(o,version);
135  }
136  else if (_saveFormat == 2)
137  {
138  MG_TRIANGLE::enregistrer(o,version);
139  }
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  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  }
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 }
MG_TRIANGLE::get_segment1
virtual MG_SEGMENT * get_segment1(void)
Definition: mg_triangle.cpp:142
gestionversion.h
CAD4FE::MCSegment::CstrMCFaceByPlaneIntr
@ CstrMCFaceByPlaneIntr
Definition: CAD4FE_MCSegment.h:105
MG_SEGMENT::get_noeud2
virtual MG_NOEUD * get_noeud2(void)
Definition: mg_segment.cpp:113
CAD4FE::PolySurface::Contains
bool Contains(MG_FACE *__refFace)
Definition: CAD4FE_PolySurface.cpp:167
MG_IDENTIFICATEUR::get_id
unsigned long get_id()
Definition: mg_identificateur.cpp:53
CAD4FE::MCSegment
Definition: CAD4FE_MCSegment.h:50
CAD4FE::MCFace::GetPolySurface
PolySurface * GetPolySurface()
Definition: CAD4FE_MCFace.cpp:113
CAD4FE::MCTriangle::MCTriangle
MCTriangle(unsigned long num, MG_ELEMENT_TOPOLOGIQUE *topo, MCNode *mgnoeud1, MCNode *mgnoeud2, MCNode *mgnoeud3, MCSegment *mgsegment1, MCSegment *mgsegment2, MCSegment *mgsegment3)
Definition: CAD4FE_MCTriangle.cpp:44
MG_TRIANGLE::get_segment2
virtual MG_SEGMENT * get_segment2(void)
Definition: mg_triangle.cpp:147
CAD4FE::MCTriangle::_saveFormat
char _saveFormat
Definition: CAD4FE_MCTriangle.h:62
MG_ELEMENT_MAILLAGE::change_lien_topologie
void change_lien_topologie(MG_ELEMENT_TOPOLOGIQUE *topo)
Definition: mg_element_maillage.cpp:56
MG_TRIANGLE
Definition: mg_triangle.h:38
MG_TRIANGLE::noeud3
class MG_NOEUD * noeud3
Definition: mg_triangle.h:76
MG_TRIANGLE::get_segment3
virtual MG_SEGMENT * get_segment3(void)
Definition: mg_triangle.cpp:152
CAD4FE_MCTriangle.h
MG_ELEMENT_TOPOLOGIQUE
Definition: mg_element_topologique.h:51
CAD4FE::MCTriangle::Tessellate
void Tessellate(unsigned __nbPtsPerSegment, std::vector< MCNode * > &__result, std::vector< unsigned > &__indexedTriangleSet)
Definition: CAD4FE_MCTriangle.cpp:373
OT_REFERENCE::get_nb_reference
int get_nb_reference(void)
Definition: ot_reference.cpp:50
CAD4FE::MCTriangle::evaluer_geo_isoparam_u
MCSegment * evaluer_geo_isoparam_u(double __u)
Definition: CAD4FE_MCTriangle.cpp:143
MAGIC
Definition: mg_fast_marching.cpp:40
CAD4FE::MCTriangle::evaluer_euc
void evaluer_euc(double __uv[2], double *__result, double *tangent_u, double *tangent_v)
Definition: CAD4FE_MCTriangle.cpp:255
CAD4FE::MCSegment::CstrNoRefTopologyMapping
@ CstrNoRefTopologyMapping
Definition: CAD4FE_MCSegment.h:108
CAD4FE::MCTriangle::~MCTriangle
virtual ~MCTriangle()
Definition: CAD4FE_MCTriangle.cpp:62
MG_SEGMENT::get_noeud1
virtual MG_NOEUD * get_noeud1(void)
Definition: mg_segment.cpp:108
CAD4FE::MCTriangle::GetRefSegment
MCSegment * GetRefSegment(int __index)
Definition: CAD4FE_MCTriangle.cpp:100
CAD4FE::MCTriangle
Definition: CAD4FE_MCTriangle.h:39
MG_ELEMENT_MAILLAGE::liaison_topologique
MG_ELEMENT_TOPOLOGIQUE * liaison_topologique
Definition: mg_element_maillage.h:69
MG_NOEUD::get_coord
virtual double * get_coord(void)
Definition: mg_noeud.cpp:92
MG_TRIANGLE::noeud2
class MG_NOEUD * noeud2
Definition: mg_triangle.h:75
CAD4FE::MCTriangle::GetRefNode
MCNode * GetRefNode(int __index)
Definition: CAD4FE_MCTriangle.cpp:85
CAD4FE_PolySurface.h
CAD4FE::MCNode::FMapIterator
FMap::iterator FMapIterator
Definition: CAD4FE_MCNode.h:51
MG_TRIANGLE::get_noeud2
virtual MG_NOEUD * get_noeud2(void)
Definition: mg_triangle.cpp:131
MG_TRIANGLE::enregistrer
virtual void enregistrer(std::ostream &o, double version)
Definition: mg_triangle.cpp:288
CAD4FE_MCSegment.h
CAD4FE_MCNode.h
MG_TRIANGLE::get_noeud1
virtual MG_NOEUD * get_noeud1(void)
Definition: mg_triangle.cpp:126
CAD4FE_Geometric_Tools.h
CAD4FE::MCTriangle::evaluer_geo
void evaluer_geo(double __uv[2], MCNode *__result, double *tangent_u, double *tangent_v)
Definition: CAD4FE_MCTriangle.cpp:218
MG_ELEMENT_MAILLAGE::get_lien_topologie
MG_ELEMENT_TOPOLOGIQUE * get_lien_topologie(void)
Definition: mg_element_maillage.cpp:51
CAD4FE::MCTriangle::ComputeIndexNodeA
int ComputeIndexNodeA()
Definition: CAD4FE_MCTriangle.cpp:67
MG_TRIANGLE::get_noeud3
virtual MG_NOEUD * get_noeud3(void)
Definition: mg_triangle.cpp:137
CAD4FE
Definition: CAD4FE_ClosestPoint_Segment_MG_ARETE.h:34
CAD4FE::MCSegment::change_orientation
void change_orientation(int __sense)
Definition: CAD4FE_MCSegment.cpp:150
OT_REFERENCE::incrementer
void incrementer(void)
Definition: ot_reference.cpp:40
MG_TRIANGLE::noeud1
class MG_NOEUD * noeud1
Definition: mg_triangle.h:74
CAD4FE::MCNode
Definition: CAD4FE_MCNode.h:47
CAD4FE::MCNode::GetRefFaceMapping
FMap & GetRefFaceMapping()
Definition: CAD4FE_MCNode.cpp:134
MG_SEGMENT::get_longueur
virtual double get_longueur(void)
Definition: mg_segment.cpp:125
MG_FACE
Definition: mg_face.h:34
CAD4FE::MCTriangle::enregistrer
virtual void enregistrer(std::ostream &o, double version)
Definition: CAD4FE_MCTriangle.cpp:126
ot_quadrature_gauss.h
CAD4FE::MCFace
Definition: CAD4FE_MCFace.h:50
CAD4FE::MCTriangle::_indexNodeA
int _indexNodeA
Definition: CAD4FE_MCTriangle.h:61
CAD4FE_MCFace.h
CAD4FE::MCTriangle::SetSaveFormat
void SetSaveFormat(char __format)
Definition: CAD4FE_MCTriangle.cpp:115
CAD4FE::MCSegment::get_longueur_geo
double get_longueur_geo()
Definition: CAD4FE_MCSegment.cpp:114
CAD4FE::MCSegment::evaluer_geo
void evaluer_geo(double __t, MCNode *__result, double *tangent=0, double *curvature=0)
Definition: CAD4FE_MCSegment.cpp:119