ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/CAD4FE/src/CAD4FE_MakeLoops.cpp
Revision: 569
Committed: Thu Oct 16 14:36:31 2014 UTC (10 years, 6 months ago) by foucault
File size: 14361 byte(s)
Log Message:
Mise à jour pour CAD4FE (Gilles) : operation 1 (tentative)

File Contents

# User Rev Content
1 foucault 27 //---------------------------------------------------------------------------
2    
3     #include <sstream>
4     #include <string>
5    
6     #pragma hdrstop
7    
8     #include "gestionversion.h"
9    
10     #include "ot_mathematique.h"
11     #include "ot_algorithme_geometrique.h"
12     #include "mg_coarete.h"
13     #include "mg_face.h"
14     #include "mg_arete.h"
15     #include "mg_boucle.h"
16     #include "mg_maillage.h"
17    
18 foucault 569 #include "CAD4FE_Geometric_Tools.h"
19 foucault 27
20     #include "CAD4FE_MakeLoops.h"
21     #include "CAD4FE_MCVertex.h"
22     #include "CAD4FE_MCEdge.h"
23     #include "CAD4FE_MCFace.h"
24    
25    
26     //---------------------------------------------------------------------------
27     #pragma package(smart_init)
28    
29     using namespace CAD4FE;
30    
31    
32     //---------------------------------------------------------------------------
33     MakeLoops::CoEdge::CoEdge(MG_FACE * __face, MG_ARETE * __e, int __sense):f(__face),sense(__sense),e(__e),prev(0),next(0){}
34     MakeLoops::CoEdge::~CoEdge()
35     {
36    
37     }
38     int MakeLoops::CoEdge::Sense(){return sense;}
39     MG_ARETE * MakeLoops::CoEdge::Edge(){return e;}
40     MG_FACE * MakeLoops::CoEdge::Face(){return f;}
41     bool MakeLoops::CoEdge::IsInverse(CoEdge & __mcEdgeSense)
42     {
43     return (__mcEdgeSense.e == e && __mcEdgeSense.Sense() == -Sense() ) ;
44     }
45     MG_SOMMET * MakeLoops::CoEdge::StartVertex()
46     {
47     if (Sense()==1)
48     return e->get_cosommet1()->get_sommet();
49     else
50     return e->get_cosommet2()->get_sommet();
51     }
52     MG_SOMMET * MakeLoops::CoEdge::EndVertex()
53     {
54     if (Sense()==1)
55     return e->get_cosommet2()->get_sommet();
56     else
57     return e->get_cosommet1()->get_sommet();
58     }
59     MG_COSOMMET * MakeLoops::CoEdge::StartCoVertex()
60     {
61     if (Sense()==1)
62     return e->get_cosommet1();
63     else
64     return e->get_cosommet2();
65     }
66     MG_COSOMMET * MakeLoops::CoEdge::EndCoVertex()
67     {
68     if (Sense()==1)
69     return e->get_cosommet2();
70     else
71     return e->get_cosommet1();
72     }
73     bool MakeLoops::CoEdge::IsAfter(CoEdge & __mcEdgeSense)
74     {
75     return this->StartVertex() == __mcEdgeSense.EndVertex();
76     }
77     bool MakeLoops::CoEdge::IsBefore(CoEdge & __mcEdgeSense)
78     {
79     return this->EndVertex() == __mcEdgeSense.StartVertex();
80     }
81     OT_VECTEUR_3D MakeLoops::CoEdge::StartDir()
82     {
83     OT_VECTEUR_3D dir;
84     double t=StartCoVertex()->get_t();
85     e->deriver(t,dir);
86     dir=Sense()*dir;
87     dir.norme();
88     return dir;
89     }
90    
91     OT_VECTEUR_3D MakeLoops::CoEdge::EndDir()
92     {
93     OT_VECTEUR_3D dir;
94     double t=EndCoVertex()->get_t();
95     e->deriver(t,dir);
96     dir=Sense()*dir;
97     dir.norme();
98     return dir;
99     }
100    
101     OT_VECTEUR_3D MakeLoops::CoEdge::StartMeshDir(MG_MAILLAGE * __mesh)
102     {
103     OT_VECTEUR_3D result;
104     MG_SOMMET * v1=StartVertex(), *v2=EndVertex();
105     int meshSense = 0; // if meshSense = sense Then use segment which direction is equal to sense
106     // if meshSense = -sense Then use segment which direction is opposite to sense
107     // if meshSense = 0 Then use any segment to compute mesh derivative
108     if (v1==v2)
109     meshSense = Sense();
110     else
111     meshSense = 0;
112    
113     result=MeshDir(__mesh, v1, meshSense);
114     if (result.get_longueur2()>1E100)
115     {
116     result = StartDir();
117     }
118    
119     return result;
120     }
121     OT_VECTEUR_3D MakeLoops::CoEdge::EndMeshDir(MG_MAILLAGE * __mesh)
122     {
123     OT_VECTEUR_3D result;
124     MG_SOMMET * v1=StartVertex(), *v2=EndVertex();
125     int meshSense = 0;
126     if (v1==v2) // if meshSense = sense Then use segment which direction is equal to sense
127     // if meshSense = -sense Then use segment which direction is opposite to sense
128     // if meshSense = 0 Then use any segment to compute mesh derivative
129     meshSense = -Sense();
130     else
131     meshSense = 0;
132    
133     result=MeshDir(__mesh, v2, meshSense);
134     if (result.get_longueur2()>1E100)
135     {
136     result = EndDir();
137     }
138    
139     return result;
140     }
141     OT_VECTEUR_3D MakeLoops::CoEdge::MeshDir(MG_MAILLAGE * __mesh, MG_SOMMET * __v, int __meshSense)
142     {
143     OT_VECTEUR_3D dir;
144     TPL_SET < MG_ELEMENT_MAILLAGE *> * mesh = __v->get_lien_maillage();
145     TPL_SET < MG_ELEMENT_MAILLAGE *>::ITERATEUR it;
146     MG_NOEUD * n;
147     MG_SEGMENT * vertexSeg = NULL;
148     double t1,t2,*xyz1,*xyz2, t1o, o, t2o;
149     double period=e->get_courbe()->get_periode();
150     for (n = (MG_NOEUD*)mesh->get_premier(it); n; n = (MG_NOEUD*)mesh->get_suivant(it))
151     if (!__mesh || __mesh->contient(n))
152     {
153     int nSeg = n->get_lien_segment()->get_nb();
154     TPL_LISTE_ENTITE< MG_SEGMENT * > * lstSeg = n->get_lien_segment();
155     for (int i=0; i<nSeg; i++)
156     {
157    
158     MG_SEGMENT * seg = lstSeg->get(i);
159     if (seg->get_lien_topologie() == e)
160     {
161     bool valid_sense;
162    
163     if (__meshSense == +1 || __meshSense == -1)
164     {
165     if (seg->get_noeud1()->get_lien_topologie()==__v)
166     {
167     xyz1=seg->get_noeud1()->get_coord();
168     xyz2=seg->get_noeud2()->get_coord();
169     }
170     else
171     {
172     xyz1=seg->get_noeud2()->get_coord();
173     xyz2=seg->get_noeud1()->get_coord();
174     }
175    
176     e->inverser(t1, xyz1);
177     e->inverser(t2, xyz2);
178    
179     t1o=.5*period;
180     o=t1o-t1;
181     t2o=t2+o;
182     if (t2o>period)
183     t2o -= period;
184     if (t2o<0)
185     t2o += period;
186    
187     if (__meshSense == +1)
188     valid_sense = (t1o<t2o);
189     else if (__meshSense == -1)
190     valid_sense = (t1o>t2o);
191     }
192     else
193     valid_sense = true;
194    
195     if (valid_sense)
196     {
197     if (vertexSeg == NULL || seg->get_longueur() > vertexSeg->get_longueur())
198     vertexSeg = seg;
199     }
200     }
201     }
202     }
203    
204     MG_NOEUD * vertexSegN[2];
205    
206     if (vertexSeg == NULL)
207     {
208     return OT_VECTEUR_3D(1E308,1E308,1E308);
209     }
210    
211     xyz1=vertexSeg->get_noeud1()->get_coord();
212     xyz2=vertexSeg->get_noeud2()->get_coord();
213     e->inverser(t1, xyz1);
214     e->inverser(t2, xyz2);
215    
216     t1o=.5*e->get_courbe()->get_periode();
217     o=t1o-t1;
218     t2o=t2+o;
219     if (t2o>period)
220     t2o -= period;
221     if (t2o<0)
222     t2o += period;
223    
224     dir = OT_VECTEUR_3D(xyz2);
225     dir -= OT_VECTEUR_3D(xyz1);
226     dir /= Sense()*(t2o-t1o);
227     return dir;
228     }
229    
230     double MakeLoops::CoEdge::MeshAngleInPlane(CoEdge & __coEdge, OT_VECTEUR_3D & __normal, MG_MAILLAGE * __mesh)
231     {
232     OT_MATRICE_3D repereNormal, transform3DToRepereNormal;
233     OT_VECTEUR_3D dir1_3D, dir2_3D, dir1_2D, dir2_2D, normal_2D;
234     if (IsInverse(__coEdge))
235     return -M_PI;
236     if (IsAfter(__coEdge))
237     {
238     dir2_3D=StartMeshDir(__mesh);
239     dir1_3D=__coEdge.EndMeshDir(__mesh);
240     }
241     else if (IsBefore(__coEdge))
242     {
243     dir1_3D=EndMeshDir(__mesh);
244     dir2_3D=__coEdge.StartMeshDir(__mesh);
245     }
246     else return 1E308;
247     double angle = AngleInPlane(dir1_3D, dir2_3D, __normal);
248     printf("Angle Between %d and %d = %f\n", e->get_id(), __coEdge.e->get_id(), angle);
249     return angle;
250     }
251    
252     double MakeLoops::CoEdge::AngleInPlane(CoEdge & __coEdge, OT_VECTEUR_3D & __normal)
253     {
254     OT_VECTEUR_3D dir1_3D, dir2_3D;
255     if (IsInverse(__coEdge))
256     return -M_PI;
257     if (IsAfter(__coEdge))
258     {
259     dir2_3D=StartDir();
260     dir1_3D=__coEdge.EndDir();
261     }
262     else if (IsBefore(__coEdge))
263     {
264     dir1_3D=EndDir();
265     dir2_3D=__coEdge.StartDir();
266     }
267     else return 1E308;
268     double angle = AngleInPlane(dir1_3D, dir2_3D, __normal);
269     printf("Angle Between %d and %d = %f\n", e->get_id(), __coEdge.e->get_id(), angle);
270     return angle;
271     }
272    
273     double MakeLoops::CoEdge::AngleInPlane(OT_VECTEUR_3D & dir1_3D, OT_VECTEUR_3D & dir2_3D, OT_VECTEUR_3D & __normal)
274     {
275     OT_VECTEUR_3D dir1_2D, dir2_2D;
276     OT_MATRICE_3D repereNormal, transform3DToRepereNormal;
277     repereNormal = GeometricTools::GetPlaneFrame(__normal);
278     repereNormal.transpose(transform3DToRepereNormal);
279     dir1_2D=transform3DToRepereNormal*dir1_3D;
280     dir2_2D=transform3DToRepereNormal*dir2_3D;
281     dir2_2D[2]=dir1_2D[2]=0;
282     dir1_2D = dir1_2D;
283     double x=1/(dir1_2D.get_longueur()*dir2_2D.get_longueur());
284     double cs = x*(dir1_2D*dir2_2D);
285     OT_VECTEUR_3D Y=OT_VECTEUR_3D(0,0,1)&dir1_2D;
286     double sn = dir2_2D*Y;
287     if (cs>1)cs=1;else if (cs<-1)cs=-1;
288     double angle = acos (cs);
289     if (sn<0) angle = -angle;
290     return angle;
291     }
292    
293     //---------------------------------------------------------------------------
294     MakeLoops::MakeLoops(std::vector <MG_FACE*> __faceList, std::vector <MG_ARETE*> __edgeList, std::vector <int> __senseList)
295     {
296     for (unsigned i=0; i<__edgeList.size(); i++)
297     {
298     MG_FACE * face = __faceList[i];
299     MG_ARETE * edge = __edgeList[i];
300     int sense = __senseList[i];
301     if (_mapFaceCoEdge.find(face) == _mapFaceCoEdge.end())
302     {
303     std::set < CoEdge * > faceCoEdges;
304     faceCoEdges.insert(new CoEdge(face, edge, sense));
305     _mapFaceCoEdge[face] = faceCoEdges;
306     }
307     else
308     {
309     _mapFaceCoEdge[face].insert(new CoEdge(face, edge, sense));
310     }
311     }
312     }
313     //---------------------------------------------------------------------------
314     MakeLoops::~MakeLoops()
315     {
316     for (std::map < MG_FACE * , std::set < CoEdge * > >::iterator itLst = _mapFaceCoEdge.begin();
317     itLst != _mapFaceCoEdge.end();
318     itLst++)
319     {
320     std::set < CoEdge * > & lst = itLst->second;
321     for ( std::set < CoEdge * >::iterator itCoEdge = lst.begin(); itCoEdge != lst.end(); itCoEdge++)
322     {
323     CoEdge * coedge = *itCoEdge;
324     delete coedge;
325     }
326     }
327     }
328     std::string MakeLoops::PrintFaceNormalAtVertices(MG_FACE * __face)
329     {
330     std::stringstream out;
331     unsigned N=0;
332     MCFace * mcFace = (MCFace*) __face;
333     out << "BaseColor { rgb 0.5 0 0 }\n";
334     out << "\n Coordinate3 {\n point [ \n";
335     std::set < CoEdge * > & lst = _mapFaceCoEdge[__face];
336     for (std::set < CoEdge * >::iterator itCoEdge = lst.begin();
337     itCoEdge != lst.end();
338     itCoEdge++)
339     {
340     CoEdge * current = *itCoEdge;
341     MCVertex * vertex = (MCVertex *)current->StartVertex();
342     OT_VECTEUR_3D normal(0,0,0);
343     int nbRefFaces;
344     mcFace->calcul_normale_unitaire(vertex,normal,&nbRefFaces);
345     double xyz[3],xyz2[3];
346     vertex->get_point()->evaluer(xyz);
347     for (int i=0; i<3; i++) xyz2[i] = xyz[i]+normal[i]*.005;
348     out << xyz[0] <<" "<< xyz[1] <<" "<< xyz[2]<<",\n";
349     out << xyz2[0] <<" "<< xyz2[1] <<" "<< xyz2[2]<<",\n";
350     N+=2;
351     }
352     out << "\n]\n}\n";
353     out << "\nIndexedLineSet {";
354     out << "\ncoordIndex\n [ \n";
355     for (unsigned int j=0; j+1<N; j+=2)
356     {
357     out << j << ", ";
358     out << j+1 << ", ";
359     out << "-1,\n";
360     }
361     out << "] \n}\n";
362    
363     return out.str();
364     }
365     //---------------------------------------------------------------------------
366     void MakeLoops::GetFaceLoops(MG_FACE* __face, std::vector < std::vector < CoEdge *> > & __loops )
367     {
368     std::set < CoEdge * > & lst = _mapFaceCoEdge[__face];
369     std::set < CoEdge * > unvisited = lst;
370    
371 foucault 64 if (lst.size()==0) return;
372    
373 foucault 27 bool isMCT;
374     {
375     std::string MCSTR("MC");
376     std::string idorig = __face->get_idoriginal();
377     string::size_type loc = idorig.find( MCSTR, 0 );
378     isMCT = ( loc == 0 );
379     }
380    
381     std::vector<CoEdge*> currentLoop;
382     CoEdge * current = 0;
383    
384     while (unvisited.size())
385     {
386     if (current == 0 || current->next != 0 )
387     {
388     if (currentLoop.size())
389     {
390     __loops.push_back(currentLoop);
391     currentLoop.clear();
392     }
393     printf("Loop %d of Face %d\n", __loops.size(), __face->get_id());
394     current = *(unvisited.begin());
395     }
396     currentLoop.push_back(current);
397     unvisited.erase(current);
398     printf("Edge %d\n", current->e->get_id());
399    
400     double score_max = -(+M_PI+.001);
401     OT_VECTEUR_3D normal;
402     if (isMCT)
403     {
404     int nbRefFaceNormal;
405     MCVertex * mcVertex = (MCVertex*) current->EndVertex();
406     MCFace * mcFace = (MCFace*)__face;
407     mcFace->calcul_normale_unitaire(mcVertex, normal, &nbRefFaceNormal);
408     }
409     else
410     {
411     MG_SOMMET * vertex = current->EndVertex();
412     double xyzVertex[3]; vertex->get_point()->evaluer(xyzVertex);
413     double uvVertex[2]; __face->inverser(uvVertex, xyzVertex);
414     __face->calcul_normale_unitaire(uvVertex, normal);
415     }
416    
417    
418     for (std::set<CoEdge*>::const_iterator it = unvisited.begin(); it != unvisited.end(); it++)
419     {
420     CoEdge * candidate = *it;
421     if (candidate->IsAfter(*current))
422     {
423     double score = current->MeshAngleInPlane(*candidate, normal,0);
424     if (score > score_max)
425     {
426     score_max = score;
427     current->next = candidate;
428     }
429     }
430     }
431    
432     current = current->next;
433     }
434    
435     __loops.push_back(currentLoop);
436     int N=0;
437     for (unsigned k=0; k<__loops.size(); k++)
438     N+=__loops[k].size();
439     printf("N=%d lst=%d\n", N, lst.size());
440     }