Moka libraries
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
corefine-3d.cc
Go to the documentation of this file.
1 /*
2  * lib-corefinement : Opérations de corafinement.
3  * Copyright (C) 2004, Moka Team, Université de Poitiers, Laboratoire SIC
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * Copyright (C) 2009, Guillaume Damiand, CNRS, LIRIS,
6  * guillaume.damiand@liris.cnrs.fr, http://liris.cnrs.fr/
7  *
8  * This file is part of lib-corefinement
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU Lesser General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public License
21  * along with this program. If not, see <http://www.gnu.org/licenses/>.
22  */
23 
24 #include <cmath>
25 #include <sstream>
26 #include "corefine-3d.hh"
27 #include "time.hh"
28 
29 #include INCLUDE_NON_INLINE("corefine-3d.icc")
30 
31 using namespace std;
32 using namespace GMap3d;
33 
34 //******************************************************************************
35 
36 #define a0 FMap->alpha0
37 #define a1 FMap->alpha1
38 #define a2 FMap->alpha2
39 #define a3 FMap->alpha3
40 
41 #define VTX(d) (AVertexDI < 0 ? FMap->findVertex(d) : \
42  (CAttributeVertex*)FMap->getDirectInfo(d, AVertexDI))
43 
44 //******************************************************************************
45 
46 // Résolutions miminale et maximale de la grille selon une direction
47 #define MIN_GRID_RES 2
48 #define MAX_GRID_RES 512
49 
50 // La grille est de taille fixe si l'option est activée
51 // #define FIXED_GRID_RES 32
52 
53 // Les sous-grilles sont de tailles fixes si l'option est activée
54 // #define FIXED_SUB_GRID_RES 8
55 
56 // Version du co-raffinement (optionnel)
57 #define COREFINEMENT_VERSION "Co-raffinement v11"
58 
59 //******************************************************************************
60 
61 /* Permet d'extraire les lignes de coupes.
62  */
63 // #define EXTRACT_LINES
64 
65 /* Permet de sauver les points d'intersection dans un fichier nommé
66  * intersections.dat
67  */
68 // #define SAVE_INTERSECTION_POINTS
69 
70 /* Permet de sauvegarder la progression des arêtes d'intersection dans des
71  * fichiers portant comme nom le numéro de l'arête.
72  */
73 // #define FOLLOW_DEBUG_MODE
74 
75 /* Permet d'afficher dans la console le détail des calculs et la progression
76  * de l'algorithme. Très utile lorsque ça plante pour connaître l'endroit exact
77  * de l'erreur.
78  */
79 // #define DEBUG_MESSAGES
80 #define WARNING_MESSAGES
81 #include "message-macros.hh"
82 
83 //******************************************************************************
84 
85 #ifdef SAVE_INTERSECTION_POINTS
86 
87 FILE * point_file = NULL;
88 
89 #define SAVE_POINT(p) (fprintf(point_file ? point_file : stdout, \
90  "%.6f # %.6f # %.6f\n", \
91  (p).getX()*10, (p).getY()*10, (p).getZ()*10))
92 
93 #else // SAVE_INTERSECTION_POINTS
94 
95 #define SAVE_POINT(p)
96 
97 #endif // SAVE_INTERSECTION_POINTS
98 
99 //******************************************************************************
100 
101 CCorefine3d::CCorefine3d(CGMapVertex * AMap, bool ACalculateOrientation,
102  TCoordinate AEpsilon, int AVertexDI)
103  : CCorefine(AMap, AEpsilon), FCalculateOrientation(ACalculateOrientation),
104  FComputeOnlyFirstIntersection(false), FGridResolution(70),
105  FTools(AMap, AEpsilon)
106 {
107  ENTER;
108 
109  if (FEps > 0)
110  CBoundingBox::setEpsilon(2*FEps);
111 
112  if (AVertexDI < 0) {
113  FVertexDI = FMap->getNewDirectInfo();
114  FLocalVertexDirectInfo = true;
115  }
116  else {
117  FVertexDI = AVertexDI;
118  FLocalVertexDirectInfo = false;
119  }
120 
121  EXIT;
122 }
123 
124 //******************************************************************************
125 
126 CCorefine3d::~CCorefine3d()
127 {
128  ENTER;
129 
130  if (FLocalVertexDirectInfo)
131  FMap->freeDirectInfo(FVertexDI);
132 
133  EXIT;
134 }
135 
136 //******************************************************************************
137 
138 void CCorefine3d::computeOnlyFirstIntersection(bool ABoolean)
139 {
140  FComputeOnlyFirstIntersection = ABoolean;
141 }
142 
143 //******************************************************************************
144 
145 void CCorefine3d::setGridResolution(int ARes)
146 {
147  FGridResolution = ARes;
148 }
149 
150 //******************************************************************************
151 
152 int CCorefine3d::corefine(CDart *& AMesh1, CDart *& AMesh2,
153  bitset<NB_MARKS> ACopyMarks)
154 {
155  ENTER;
156 
157  list<CDart*> vertex_list;
158  CEdgeIntersection inter;
159  unsigned long nb_faces;
160  CDart *vertex=NULL, *edge=NULL, *tmp_edge=NULL, *tmp_face=NULL;
161  int edge_mark = FMap->getNewMark();
162  int vertex_mark = FMap->getNewMark();
163  CVertex *pt1, *pt2, tmp_pt;
164  CPlane plane;
165  CTime start, end;
166 
167  FFaceDI = FMap->getNewDirectInfo();
168  FAlpha2DI = FMap->getNewDirectInfo();
169 
170  FOrientMark = FMap->getNewMark();
171  FFictiveMark = FMap->getNewMark();
172  FIntersectionMark = FMap->getNewMark();
173 
174  FCopyMarks = ACopyMarks;
175 
176  FNumberOfIntersectionLines = 0;
177  FNumberOfIntersectionEdges = 0;
178 
179 // if (!FTools.isMeshOk(AMesh1) || !FTools.isMeshOk(AMesh2)) {
180 // cout << "\033[1;31m"
181 // << "Co-raffinement impossible !\n"
182 // << "Les maillages contiennent des arêtes de longueur nulle."
183 // << "\033[0m\n";
184 
185 // FMap->freeMark(edge_mark);
186 // FMap->freeMark(vertex_mark);
187 // return;
188 // }
189 
190 #ifdef COREFINEMENT_VERSION
191  cout << COREFINEMENT_VERSION << endl;
192 #endif
193 
194 #ifdef SAVE_INTERSECTION_POINTS
195  point_file = fopen("intersections.dat", "w");
196 #endif
197 
198  cout << "Initialisation des maillages" << endl;
199  start.updateTime();
200 
201  nb_faces = initMesh(AMesh1);
202  cout << "Nombre de faces sur le premier maillage : " << nb_faces << endl;
203 
204  nb_faces = initMesh(AMesh2);
205  cout << "Nombre de faces sur le second maillage : " << nb_faces << endl;
206 
207  end.updateTime();
208  FInitialisationTime = end - start;
209  cout << "Durée de l'initialisation : " << FInitialisationTime << endl;
210 
211  cout << "Orientation des maillages" << endl;
212  start.updateTime();
213 
214  // Recherche d'un brin correctement orienté sur chacun des maillages
215  if (FCalculateOrientation) {
216  AMesh1 = FTools.findWellOrientedDart(AMesh1, FVertexDI);
217  AMesh2 = FTools.findWellOrientedDart(AMesh2, FVertexDI);
218  }
219 
220  // Orientation des maillages
221  FMap->halfMarkOrbit(AMesh1, ORBIT_CC, FOrientMark);
222  FMap->halfMarkOrbit(AMesh2, ORBIT_CC, FOrientMark);
223 
224  end.updateTime();
225  cout << "Durée de l'orientation : " << end - start << endl;
226 
227  cout << "Création de la grille régulière pour l'optimisation" << endl;
228  start.updateTime();
229 
230  createGrid(AMesh2, nb_faces);
231 
232 // int depth = 0;
233 // while (depth < 5 && refineGrid(MAX_GRID_RES / 16, FNumberOfFacesPerCell) > 0)
234 // depth++;
235 // cout << "Profondeur de la grille : " << depth << endl;
236 
237  end.updateTime();
238  FGridCreationTime = end - start;
239  cout << "Durée de la création : " << FGridCreationTime << endl;
240 
241  FLineCreationTime.setTime(0);
242 
243  cout << "Traitement de la première intersection" << endl;
244  start.updateTime();
245 
246  // Gestion de la première intersection :
247  int tmp_dim;
248  tmp_pt = *getVertex(AMesh1);
249  tmp_face = findFirstIntersectionInGrid(tmp_pt, &tmp_dim);
250 
251  if (tmp_face != NULL) {
252  switch (tmp_dim) {
253  case 0:
254  *getVertex(AMesh1) = *getVertex(tmp_face);
255  break;
256 
257  case 1:
258  tmp_face = splitEdge(tmp_face, tmp_pt);
259  break;
260 
261  case 2:
262  tmp_face = insertVertexInFace(tmp_face, tmp_pt);
263  break;
264  }
265 
266  followIntersection(AMesh1, tmp_face, vertex_mark);
267  }
268 
269  cout << "Parcours des maillages et création des intersections" << endl;
270 
271  // Ajoute le sommet courant dans la liste
272  vertex_list.push_back(AMesh1);
273 
274  // Tant qu'il y a des éléments dans la liste, on continue le parcours
275  while(!vertex_list.empty()) {
276 
277  // Récupération du prochain sommet à traiter
278  vertex = vertex_list.front(), vertex_list.pop_front();
279 
280  // Parcours des arêtes incidentes au sommet
281  for (CStaticCoverageVertex dcv(FMap, vertex); dcv.cont(); ++dcv) {
282  edge = *dcv;
283 
284  if (!FMap->isMarked(edge, edge_mark) &&
285  FMap->isMarked(edge, FOrientMark)) {
286  pt1 = getVertex(edge);
287  pt2 = getVertex(a0(edge));
288 
289  MSG("Arête parcourue = [" << *pt1 << "," << *pt2 << "]");
290 
291  if (!FMap->isMarked(edge, FFictiveMark))
292  inter = findNearestIntersectionInGrid(*pt1, *pt2);
293  else
294  inter.setCell(NULL);
295 
296  if (inter.getCell() != NULL && // FNumberOfIntersectionLines == 0 &&
297 // (inter.getPosition() != EP_OnSecondVertex ||
298  ((*pt2 - inter.getPoint()).norm() > 2.0 * FEps ||
299  !FMap->isMarked(a0(edge), vertex_mark))) {
300 
301  assert(inter.getPosition() != EP_OnFirstVertex);
302 
303 // cout << "Une intersection trouvée sur l'arête ["
304 // << *pt1 << "," << *pt2 << "] de longueur "
305 // << (*pt2 - *pt1).norm() << " avec la face :" << endl;
306 // FTools.displayFaceVertices(inter.getCell(), FVertexDI);
307 // cout << "Propriétés de l'intersection :" << endl << inter;
308 
309  switch (inter.getPosition()) {
310  case EP_OnFirstVertex:
311  tmp_edge = edge;
312  *pt1 = inter.getPoint();
313  break;
314 
315  case EP_OnSecondVertex:
316  tmp_edge = a3(a0(edge));
317  *pt2 = inter.getPoint();
318  break;
319 
320  case EP_OnEdge:
321  tmp_edge = splitEdge(edge, inter.getPoint());
322  break;
323 
324  default:
325  break;
326  }
327 
328  switch (inter.getCellDimension()) {
329  case 0:
330  tmp_face = inter.getCell();
331  *getVertex(inter.getCell()) = inter.getPoint();
332  break;
333 
334  case 1:
335  tmp_face = splitEdge(inter.getCell(), inter.getPoint());
336  break;
337 
338  case 2:
339  tmp_face = insertVertexInFace(inter.getCell(), inter.getPoint());
340  break;
341 
342  default:
343  tmp_face = NULL;
344  break;
345  }
346 
347  followIntersection(tmp_edge, tmp_face, vertex_mark);
348  }
349  else {
350  tmp_edge = a3(a0(edge));
351  }
352 
353  if (!FMap->isMarked(tmp_edge, vertex_mark) ||
354  !FComputeOnlyFirstIntersection)
355  vertex_list.push_back(tmp_edge);
356 
357  FMap->markOrbit(edge, ORBIT_EDGE, edge_mark);
358  }
359  }
360  }
361 
362  FMap->unmarkOrbit(AMesh1, ORBIT_CC, edge_mark);
363  FMap->unmarkOrbit(AMesh1, ORBIT_CC, vertex_mark);
364 
365  if (FNumberOfIntersectionEdges == 0)
366  FMap->unmarkOrbit(AMesh2, ORBIT_CC, vertex_mark);
367 
368  FMap->freeMark(edge_mark);
369  FMap->freeMark(vertex_mark);
370 
371  end.updateTime();
372  FResearchTime = end - start;
373  cout << "Durée de la recherche des intersections : "
374  << FResearchTime << endl;
375 
376  FResearchTime -= FLineCreationTime;
377 
378  cout << "Assemblage des maillages" << endl;
379  start.updateTime();
380 
381  applyModifications(AMesh1);
382 
383  end.updateTime();
384  FUpdateTime = end - start;
385  cout << "Durée de l'assemblage : " << FUpdateTime << endl;
386 
387 #ifdef EXTRACT_LINES
388  cout << "Extraction des lignes de coupes" << endl;
389  extractIntersectionLines(AMesh1);
390  FMap->unmarkOrbit(AMesh1, ORBIT_CC, FIntersectionMark);
391 #endif // EXTRACT_LINES
392 
393  cout << "Destruction de la grille" << endl;
394  start.updateTime();
395 
396  destroyGrid();
397 
398  end.updateTime();
399  cout << "Durée de la destruction : " << end - start << endl;
400 
401  cout << "Nettoyage des maillages" << endl;
402  start.updateTime();
403 
404  cleanMesh(AMesh1);
405 
406  if (FNumberOfIntersectionEdges == 0)
407  cleanMesh(AMesh2);
408 
409  end.updateTime();
410  cout << "Durée du nettoyage : " << end - start << endl;
411 
412  cout << "Il y a eu " << FNumberOfIntersectionLines << " ligne(s) de coupe et "
413  << FNumberOfIntersectionEdges << " arête(s) d'intersection" << endl;
414 
415 #ifdef SAVE_INTERSECTION_POINTS
416  if (point_file != NULL) {
417  fclose(point_file);
418  point_file = NULL;
419  }
420 #endif
421 
422  FMap->freeDirectInfo(FFaceDI);
423  FMap->freeDirectInfo(FAlpha2DI);
424 
425  FMap->freeMark(FOrientMark);
426  FMap->freeMark(FFictiveMark);
427  FMap->freeMark(FIntersectionMark);
428 
429  EXIT;
430 
431  return FNumberOfIntersectionLines;
432 }
433 
434 //******************************************************************************
435 
436 unsigned long CCorefine3d::initMesh(CDart * AMesh)
437 {
438  ENTER;
439 
440  int mark = FMap->getNewMark();
441  unsigned long int nb_faces = 0;
442 
443  MSG("mark = " << mark);
444 
445  // Fermeture des 3-bords
446  for (CStaticCoverageCC scc(FMap, AMesh); scc.cont(); ++scc) {
447  if (FMap->isFree1(*scc))
448  cerr << "<initMesh> Le brin " << *scc
449  << " n'est pas lié par alpha1 !!!" << endl;
450 
451  if (FMap->isFree3(*scc))
452  FMap->stopUp(*scc, 3);
453  }
454 
455  // Fermeture des 2-bords
456 // for (CStaticCoverageCC scc(FMap, AMesh); scc.cont(); ++scc) {
457 // if (FMap->isFree2(*scc) && FMap->canSew2(*scc, a3(*scc)))
458 // FMap->sew2(*scc, a3(*scc));
459 // }
460 
461  CDynamicCoverageCC dcc(FMap, AMesh);
462 
463 // if (FLocalVertexDirectInfo) {
464  // Mise à jour des pointeurs vers les plongements sommets
465  for (dcc.reinit(); dcc.cont(); ++dcc) {
466  if (!FMap->isMarked(*dcc, mark)) {
467  FMap->markOrbit(*dcc, ORBIT_VERTEX, mark);
468  updateVertexLinks(*dcc);
469  }
470  }
471 // }
472 // else {
473 // for (dcc.reinit(); dcc.cont(); ++dcc)
474 // FMap->setMark(*dcc, mark);
475 // }
476 
477  // Mise à jour des pointeurs vers un brin des faces
478  for (dcc.reinit(); dcc.cont(); ++dcc) {
479  FMap->setDirectInfo(*dcc, FAlpha2DI, NULL);
480 
481  if (FMap->isMarked(*dcc, mark)) {
482  nb_faces++;
483  FMap->unmarkOrbit(*dcc, ORBIT_FACE, mark);
484  updateFaceLinks(*dcc);
485  }
486  }
487 
488  FMap->freeMark(mark);
489 
490  EXIT;
491 
492  return nb_faces;
493 }
494 
495 //******************************************************************************
496 
497 void CCorefine3d::cleanMesh(CDart * /*AMesh*/)
498 {
499  ENTER;
500 
501  int delete_mark = FMap->getNewMark();
502 
503  MSG("delete_mark = " << delete_mark);
504 
505  for (CDynamicCoverageAll dca(FMap); dca.cont(); ++dca) {
506  if (FMap->isMarked(*dca, FOrientMark)) {
507  FMap->unmarkOrbit(*dca, ORBIT_EDGE, FOrientMark);
508 
509  if (FMap->isMarked(*dca, FFictiveMark)) {
510  if (getFace(*dca) != getFace(a2(*dca))) {
511  removeEdge(*dca, delete_mark);
512  }
513  else
514  FMap->unmarkOrbit(*dca, ORBIT_EDGE, FFictiveMark);
515  }
516  }
517  }
518  FMap->deleteMarkedDarts(delete_mark);
519  FMap->freeMark(delete_mark);
520 
521  EXIT;
522 }
523 
524 //******************************************************************************
525 
526 // void CCorefine3d::createGrid(CDart * AMesh, int AMaxDiv)
527 // {
528 // ENTER;
529 
530 // int face_mark = FMap->getNewMark();
531 // CBoundingBox bb = FTools.orbitBoundingBox(AMesh, ORBIT_CC, FVertexDI);
532 // CVertex size;
533 // TCoordinate max;
534 // int size_i, size_j, size_k;
535 
536 // // Calcul de la résolution de la grille
537 // size = bb.getEpsMaxBound() - bb.getEpsMinBound();
538 // max = (size.getX() > size.getY() ?
539 // (size.getX() > size.getZ() ? size.getX() : size.getZ()) :
540 // (size.getY() > size.getZ() ? size.getY() : size.getZ()));
541 // size /= max;
542 // size_i = (int)(size.getX() * AMaxDiv);
543 // size_j = (int)(size.getY() * AMaxDiv);
544 // size_k = (int)(size.getZ() * AMaxDiv);
545 // if (size_i == 0) size_i = 1;
546 // if (size_j == 0) size_j = 1;
547 // if (size_k == 0) size_k = 1;
548 
549 // cout << "Résolution de la grille = "
550 // << size_i << "x" << size_j << "x" << size_k << endl;
551 
552 // FGrid = new CGrid3d<list<CDart*>*>(size_i, size_j, size_k, bb);
553 // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid) ; gi.cont() ; ++gi)
554 // FGrid->setCell(gi, new list<CDart*>);
555 
556 // for (CDynamicCoverageCC dcc(FMap, AMesh); dcc.cont(); ++dcc) {
557 // if (!FMap->isMarked(*dcc, face_mark) &&
558 // FMap->isMarked(*dcc, FOrientMark)) {
559 // FMap->markOrbit(*dcc, ORBIT_FACE, face_mark);
560 // addFaceToGrid(getFace(*dcc));
561 // }
562 // }
563 
564 // FMap->unmarkOrbit(AMesh, ORBIT_CC, face_mark);
565 // FMap->freeMark(face_mark);
566 
567 // // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid) ; gi.cont() ; ++gi)
568 // // cout << gi << " : " << (*gi)->size() << endl;
569 
570 // EXIT;
571 // }
572 
573 //******************************************************************************
574 
575 // void CCorefine3d::destroyGrid()
576 // {
577 // ENTER;
578 
579 // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid) ; gi.cont() ; ++gi)
580 // delete FGrid->getCell(gi);
581 // delete FGrid;
582 
583 // EXIT;
584 // }
585 
586 //******************************************************************************
587 
588 void CCorefine3d::createGrid(CDart * AMesh, unsigned long ANbFaces)
589 {
590  ENTER;
591 
592  int face_mark = FMap->getNewMark();
593  CBoundingBox bb = FTools.orbitBoundingBox(AMesh, ORBIT_CC, FVertexDI);
594  CVertex size = bb.getEpsMaxBound() - bb.getEpsMinBound();
595  unsigned int size_i, size_j, size_k;
596 
597  MSG("face_mark = " << face_mark);
598 
599  // Calcul de la résolution de la grille
600  if (FGridResolution > 0)
601  getGridResolution(normalizeGridSize(size), 1, FGridResolution,
602  &size_i, &size_j, &size_k);
603  else
604  computeGridResolution(normalizeGridSize(size),
605  ANbFaces, 0.05,
606  &size_i, &size_j, &size_k,
607  MIN_GRID_RES, MAX_GRID_RES);
608 
609  cout << "Résolution de la grille = "
610  << size_i << "x" << size_j << "x" << size_k << endl;
611 
612  FGrid = new TFaceGrid(size_i, size_j, size_k, bb);
613  for (TFaceGridIter gi(*FGrid); gi.cont(); ++gi)
614  FGrid->setCell(gi, new list<CDart*>);
615 
616  for (CDynamicCoverageCC dcc(FMap, AMesh); dcc.cont(); ++dcc) {
617  if (!FMap->isMarked(*dcc, face_mark) &&
618  FMap->isMarked(*dcc, FOrientMark)) {
619  FMap->markOrbit(*dcc, ORBIT_FACE, face_mark);
620  addFaceToGrid(FGrid, *dcc);
621  }
622  }
623 
624  FMap->unmarkOrbit(AMesh, ORBIT_CC, face_mark);
625  FMap->freeMark(face_mark);
626 
627 // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid); gi.cont(); ++gi)
628 // cout << gi << " : " << (*gi)->size() << endl;
629 
630  EXIT;
631 }
632 
633 //******************************************************************************
634 
635 void CCorefine3d::destroyGrid()
636 {
637  ENTER;
638 
639  for (TFaceGridIter gi(*FGrid); gi.cont(); ++gi)
640  delete *gi;
641  delete FGrid;
642 
643  EXIT;
644 }
645 
646 //******************************************************************************
647 
648 // void CCorefine3d::addFaceToGrid(CDart * AFace)
649 // {
650 // ENTER;
651 
652 // CBoundingBox bb = FTools.orbitBoundingBox(AFace, ORBIT_01, FVertexDI);
653 
654 // // int x, y, z;
655 // // FGrid->getCellPosition(bb->getEpsMinBound(), &x, &y, &z);
656 // // cout << "min = (" << x << "," << y << "," << z << ")" << endl;
657 // // FGrid->getCellPosition(bb->getEpsMaxBound(), &x, &y, &z);
658 // // cout << "max = (" << x << "," << y << "," << z << ")" << endl;
659 
660 // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid, &bb); gi.cont(); ++gi) {
661 // (*gi)->push_back(AFace);
662 // }
663 
664 // EXIT;
665 // }
666 
667 //******************************************************************************
668 
669 // void CCorefine3d::removeFaceFromGrid(CDart * AFace)
670 // {
671 // ENTER;
672 
673 // CBoundingBox bb = FTools.orbitBoundingBox(AFace, ORBIT_01, FVertexDI);
674 // list<CDart*> *dart_list;
675 // list<CDart*>::iterator li, old_li;
676 
677 // for (CGrid3dIterator<list<CDart*>*> gi(*FGrid, &bb); gi.cont(); ++gi) {
678 // dart_list = *gi;
679 
680 // for (li = dart_list->begin() ; li != dart_list->end() ; ) {
681 // old_li = li;
682 // ++li;
683 // if (*old_li == AFace)
684 // dart_list->erase(old_li);
685 // }
686 // }
687 
688 // EXIT;
689 // }
690 
691 //******************************************************************************
692 
693 unsigned int CCorefine3d::getMaxVerticesDegree(list<CDart*> * AList)
694 {
695  ENTER;
696 
697  list<CDart*>::iterator it;
698  unsigned int max_deg = 0, current_deg;
699 
700  for (it = AList->begin(); it != AList->end(); it++) {
701  for (CDynamicCoverage01 dcf(FMap, *it); dcf.cont(); dcf++) {
702  if (FMap->isMarked(*dcf, FOrientMark)) {
703  current_deg = 0;
704 
705  for (CDynamicCoverageVertex dcv(FMap, *dcf); dcv.cont(); dcv++)
706  current_deg++;
707 
708  if (current_deg > max_deg)
709  max_deg = current_deg;
710  }
711  }
712  }
713 
714  EXIT;
715 
716  return max_deg / 4;
717 }
718 
719 //******************************************************************************
720 
721 CVertex CCorefine3d::normalizeGridSize(const CVertex & AGridSize)
722 {
723  TCoordinate max;
724 
725  max = (AGridSize.getX() > AGridSize.getY() ?
726  (AGridSize.getX() > AGridSize.getZ() ?
727  AGridSize.getX() : AGridSize.getZ()) :
728  (AGridSize.getY() > AGridSize.getZ() ?
729  AGridSize.getY() : AGridSize.getZ()));
730 
731  return AGridSize / max;
732 }
733 
734 //******************************************************************************
735 
736 void CCorefine3d::getGridResolution(const CVertex & AGridSize,
737  unsigned int AMinRes,
738  unsigned int AMaxRes,
739  unsigned int * AResX,
740  unsigned int * AResY,
741  unsigned int * AResZ)
742 {
743  assert(AGridSize.getX() >= 0.0 && AGridSize.getX() <= 1.0 &&
744  AGridSize.getY() >= 0.0 && AGridSize.getY() <= 1.0 &&
745  AGridSize.getZ() >= 0.0 && AGridSize.getZ() <= 1.0);
746 
747  *AResX = (int)(AGridSize.getX() * AMaxRes + 0.5);
748  *AResY = (int)(AGridSize.getY() * AMaxRes + 0.5);
749  *AResZ = (int)(AGridSize.getZ() * AMaxRes + 0.5);
750  if (*AResX < AMinRes) *AResX = AMinRes;
751  if (*AResY < AMinRes) *AResY = AMinRes;
752  if (*AResZ < AMinRes) *AResZ = AMinRes;
753 }
754 
755 //******************************************************************************
756 
757 void CCorefine3d::computeGridResolution(const CVertex & AGridSize,
758  unsigned long ANbFaces,
759  TCoordinate ANbFacesPerCell,
760  unsigned int * AResX,
761  unsigned int * AResY,
762  unsigned int * AResZ,
763  unsigned int AMinRes,
764  unsigned int /*AMaxRes*/)
765 {
766  assert(AGridSize.getX() >= 0.0 && AGridSize.getX() <= 1.0 &&
767  AGridSize.getY() >= 0.0 && AGridSize.getY() <= 1.0 &&
768  AGridSize.getZ() >= 0.0 && AGridSize.getZ() <= 1.0);
769 
770  double res = pow((double)ANbFaces /
771  (double)(ANbFacesPerCell * AGridSize.getX() *
772  AGridSize.getY() * AGridSize.getZ()),
773  1.0/3.0);
774 
775  getGridResolution(AGridSize, AMinRes, (int)(res+0.5), AResX, AResY, AResZ);
776 }
777 
778 //******************************************************************************
779 
780 unsigned int CCorefine3d::refineGrid(unsigned int AMaxSubDiv,
781  unsigned int AMaxNumberOfFaces)
782 {
783  ENTER;
784 
785  list<CDart*> *l;
786  list<CDart*>::iterator it;
787  TFaceGrid *sub_grid;
788  unsigned int max_deg, size, result = 0;
789  unsigned int div_x, div_y, div_z;
790 
791  for (TFaceGridIter gi(*FGrid); gi.cont(); ++gi) {
792  l = *gi;
793  size = l->size();
794 
795  if (size > AMaxNumberOfFaces) {
796  max_deg = getMaxVerticesDegree(l);
797 
798  if (max_deg <= AMaxNumberOfFaces || size > 2 * max_deg) {
799  result++;
800 
801 #ifdef FIXED_SUB_GRID_RES
802  div_x = div_y = div_z = FIXED_SUB_GRID_RES;
803 #else
804  computeGridResolution(CVertex(1.0, 1.0, 1.0),
805  size, 0.05,
806  &div_x, &div_y, &div_z,
807  MIN_GRID_RES, AMaxSubDiv);
808 #endif
809 
810 // cout << "Cellule " << gi << " éclatée en une grille de taille "
811 // << sub_div << "^3" << endl;
812 
813  sub_grid = FGrid->splitCellInGrid(gi, div_x, div_y, div_z);
814 
815  for (TFaceGridIter gi2(*sub_grid); gi2.cont(); ++gi2)
816  sub_grid->setCell(gi2, new list<CDart*>);
817 
818  for (it = l->begin(); it != l->end(); ++it)
819  addFaceToGrid(sub_grid, *it);
820 
821  delete l;
822  }
823  }
824  }
825 
826  EXIT;
827 
828  return result;
829 }
830 
831 //******************************************************************************
832 
833 void CCorefine3d::addFaceToGrid(TFaceGrid * AGrid, CDart * AFace)
834 {
835  ENTER;
836 
837  assert(AGrid != NULL);
838  assert(AFace != NULL);
839 
840  CBoundingBox bb = FTools.orbitBoundingBox(AFace, ORBIT_01, FVertexDI);
841 
842  for (TFaceGridIter gi(*AGrid, bb); gi.cont(); ++gi)
843  (*gi)->push_back(getFace(AFace));
844 
845  EXIT;
846 }
847 
848 //******************************************************************************
849 
850 void CCorefine3d::removeFaceFromGrid(TFaceGrid * AGrid,
851  CDart * AFace)
852 {
853  ENTER;
854 
855  AFace = getFace(AFace);
856 
857  CBoundingBox bb = FTools.orbitBoundingBox(AFace, ORBIT_01, FVertexDI);
858  list<CDart*> *dart_list;
859  list<CDart*>::iterator li, old_li;
860 
861  for (TFaceGridIter gi(*AGrid, bb); gi.cont(); ++gi) {
862  dart_list = *gi;
863 
864  for (li = dart_list->begin(); li != dart_list->end(); ) {
865  old_li = li;
866  ++li;
867  if (*old_li == AFace)
868  dart_list->erase(old_li);
869  }
870  }
871 
872  EXIT;
873 }
874 
875 //******************************************************************************
876 
877 void CCorefine3d::updateVertexLinks(CDart * ADart, CAttributeVertex * AVertex)
878 {
879  ENTER;
880 
881  for (CDynamicCoverageVertex dcv(FMap, ADart); dcv.cont(); ++dcv) {
882  FMap->setDirectInfo(*dcv, FVertexDI, AVertex);
883  }
884 
885  EXIT;
886 }
887 
888 //******************************************************************************
889 
890 void CCorefine3d::updateVertexLinks(CDart * ADart)
891 {
892  ENTER;
893 
894  CAttributeVertex * att = FMap->findVertex(ADart);
895 
896  for (CDynamicCoverageVertex dcv(FMap, ADart); dcv.cont(); ++dcv) {
897  FMap->setDirectInfo(*dcv, FVertexDI, att);
898  }
899 
900  EXIT;
901 }
902 
903 //******************************************************************************
904 
905 void CCorefine3d::updateFaceLinks(CDart * AFace)
906 {
907  ENTER;
908 
909  for (CDynamicCoverageFace dcf(FMap, AFace); dcf.cont(); ++dcf) {
910  FMap->setDirectInfo(*dcf, FFaceDI, AFace);
911  }
912 
913  EXIT;
914 }
915 
916 //******************************************************************************
917 
918 CDart * CCorefine3d::createEdge()
919 {
920  ENTER;
921 
922  CDart *d[8];
923 
924  for (int i=0; i<8; ++i) {
925  d[i] = FMap->addMapDart();
926  FMap->setDirectInfo(d[i], FAlpha2DI, NULL);
927  }
928 
929  FMap->linkAlpha0(d[0], d[1]);
930  FMap->linkAlpha0(d[2], d[3]);
931  FMap->linkAlpha0(d[4], d[5]);
932  FMap->linkAlpha0(d[6], d[7]);
933 
934  FMap->linkAlpha2(d[0], d[2]);
935  FMap->linkAlpha2(d[1], d[3]);
936  FMap->linkAlpha2(d[4], d[6]);
937  FMap->linkAlpha2(d[5], d[7]);
938 
939  FMap->linkAlpha3(d[0], d[4]);
940  FMap->linkAlpha3(d[1], d[5]);
941  FMap->linkAlpha3(d[2], d[6]);
942  FMap->linkAlpha3(d[3], d[7]);
943 
944  FMap->halfMarkOrbit(d[0], ORBIT_EDGE, FOrientMark);
945 
946 // FMap->setMark(d[0], FOrientMark);
947 // FMap->setMark(d[3], FOrientMark);
948 // FMap->setMark(d[5], FOrientMark);
949 // FMap->setMark(d[6], FOrientMark);
950 
951  EXIT;
952 
953  return d[0];
954 }
955 
956 //******************************************************************************
957 
958 CDart * CCorefine3d::insertVertexInFace(CDart * AFace, const CVertex & APoint)
959 {
960  ENTER;
961 
962  assert(FMap->isMarked(AFace, FOrientMark));
963 
964  CDart *vertex = a1(a0(insertEdgeInFace(AFace, APoint)));
965 
966  FMap->markOrbit(vertex, ORBIT_EDGE, FFictiveMark);
967 
968  EXIT;
969 
970  return vertex;
971 }
972 
973 //******************************************************************************
974 
975 CDart * CCorefine3d::insertEdgeInFace(CDart * AVertex1,
976  const CVertex & AVertex2)
977 {
978  ENTER;
979 
980  assert(FMap->isMarked(AVertex1, FOrientMark));
981 
982  CDart *edge;
983 
984 // if (FTools.isPointOnLine(AVertex2, *getVertex(AVertex1),
985 // *getVertex(a0(AVertex1)))) {
986 // edge = AVertex1;
987 // splitEdge(edge, AVertex2);
988 // if (FMap->isMarked(edge, FFictiveMark))
989 // FMap->unmarkOrbit(edge, ORBIT_EDGE, FFictiveMark);
990 // }
991 // else if (FTools.isPointOnLine(AVertex2, *getVertex(AVertex1),
992 // *getVertex(a0(a1(AVertex1))))) {
993 // edge = a3(a1(AVertex1));
994 // splitEdge(edge, AVertex2);
995 // if (FMap->isMarked(edge, FFictiveMark))
996 // FMap->unmarkOrbit(edge, ORBIT_EDGE, FFictiveMark);
997 // }
998 // else {
999  CDart *old_a1, *face = getFace(AVertex1);
1000 
1001  CPlane plane = FTools.facePlane(AVertex1, FVertexDI);
1002 
1003  assert(FTools.isPointOnPlane(AVertex2, plane));
1004 
1005  if (FTools.isVectorNull(plane.getNormal())) {
1006  cout << "\033[1;34m";
1007  FTools.displayFaceVertices(AVertex1, FVertexDI);
1008  cout << "<insertEdgeInFace> La normale à la face est nulle : "
1009  << plane.getNormal() << "\033[0m\n";
1010  }
1011 
1012  // Ici, on s'assure d'insérer l'arête au bon endroit afin de ne pas obtenir
1013  // une topologie en désaccord avec la géométrie.
1014  AVertex1 = FTools.findSectorOfVector(AVertex2 - *getVertex(AVertex1),
1015  AVertex1, plane, FVertexDI);
1016 
1017  AVertex1 = a1(AVertex1);
1018 
1019  edge = createEdge();
1020 
1021  assert(!FMap->isFree1(AVertex1));
1022 
1023  old_a1 = a1(AVertex1);
1024  FMap->unsew1(AVertex1);
1025 
1026  FMap->sew1(AVertex1, edge);
1027  FMap->sew1(old_a1, a2(edge));
1028  FMap->sew1(a0(edge), a0(a2(edge)));
1029 
1030  updateVertexLinks(AVertex1);
1031 
1032  FMap->setVertex(a0(edge), AVertex2);
1033 
1034  for (CDynamicCoverage23 dc(FMap, edge); dc.cont(); ++dc) {
1035  FMap->setDirectInfo(*dc, FFaceDI, face);
1036  FMap->setDirectInfo(a0(*dc), FFaceDI, face);
1037 
1038  for (int m = 0; m < NB_MARKS; m++)
1039  if (FCopyMarks[m] && FMap->isMarked(a1(*dc), m)) {
1040  FMap->setMark(*dc, m);
1041  FMap->setMark(a0(*dc), m);
1042  }
1043  }
1044 
1045  updateVertexLinks(a0(edge));
1046 // }
1047 
1048  EXIT;
1049 
1050  return edge;
1051 }
1052 
1053 //******************************************************************************
1054 
1055 CDart * CCorefine3d::splitFace(CDart * AVertex1, CDart * AVertex2)
1056 {
1057  ENTER;
1058 
1059  assert(AVertex1 != AVertex2);
1060  assert(FMap->isMarked(AVertex1, FOrientMark));
1061  assert(FMap->isMarked(AVertex2, FOrientMark));
1062 
1063  CDart *edge = NULL;
1064 
1065  // Gestion du cas ou l'arête est déjà existante
1066  for (CDynamicCoverage12 dc1(FMap, AVertex1); dc1.cont() && !edge; ++dc1)
1067  for (CDynamicCoverage12 dc2(FMap, AVertex2); dc2.cont() && !edge; ++dc2)
1068  if (a0(*dc1) == *dc2) {
1069  if (FMap->isMarked(*dc1, FOrientMark))
1070  edge = *dc1;
1071  else
1072  edge = a3(*dc1);
1073 
1074  if (FMap->isMarked(edge, FFictiveMark))
1075  FMap->unmarkOrbit(edge, ORBIT_EDGE, FFictiveMark);
1076  }
1077 
1078  // Cas où les deux sommets ne sont pas reliés par une arête
1079  if (edge == NULL) {
1080  CVertex v = *getVertex(AVertex2) - *getVertex(AVertex1);
1081 
1082  MSG("Création d'une arête entre " << *getVertex(AVertex1)
1083  << " et " << *getVertex(AVertex2));
1084 
1085  CPlane plane = FTools.facePlane(AVertex1, FVertexDI);
1086 
1087  if (FTools.isVectorNull(plane.getNormal())) {
1088  cout << "\033[1;34m";
1089  FTools.displayFaceVertices(AVertex1, FVertexDI);
1090  cout << "<splitFace> La normale à la face est nulle : "
1091  << plane.getNormal() << "\033[0m\n";
1092  }
1093 
1094  AVertex1 = FTools.findSectorOfVector(v, AVertex1, plane, FVertexDI);
1095  AVertex2 = FTools.findSectorOfVector(-v, AVertex2, plane, FVertexDI);
1096 
1097  CDart *old_a1[2], *face = getFace(AVertex1);
1098 
1099  edge = createEdge();
1100 
1101  // On prend le brin voisin qui est orienté dans le mauvais sens car le
1102  // premier brin de la nouvelle arête est orienté dans le bon sens
1103  AVertex1 = a1(AVertex1);
1104 
1105  old_a1[0] = a1(AVertex1);
1106  old_a1[1] = a1(AVertex2);
1107 
1108  assert(!FMap->isFree1(AVertex1));
1109  FMap->unsew1(AVertex1);
1110 
1111  assert(!FMap->isFree1(AVertex2));
1112  FMap->unsew1(AVertex2);
1113 
1114  FMap->sew1(AVertex1, edge);
1115  FMap->sew1(old_a1[0], a2(edge));
1116 
1117  FMap->sew1(AVertex2, a0(edge));
1118  FMap->sew1(old_a1[1], a2(a0(edge)));
1119 
1120  updateVertexLinks(edge);
1121  updateVertexLinks(a0(edge));
1122 
1123  if (!FMap->isSameOrbit(edge, a2(edge), ORBIT_01)) {
1124  updateFaceLinks(edge);
1125  updateFaceLinks(a2(edge));
1126  }
1127  else {
1128  for (CDynamicCoverageEdge dce(FMap, edge); dce.cont(); ++dce)
1129  FMap->setDirectInfo(*dce, FFaceDI, face);
1130  }
1131 
1132 // cleanFictiveEdges(edge);
1133 // if (!FMap->is_same_orbit(edge, a2(edge), ORBIT_01))
1134 // cleanFictiveEdges(a2(edge));
1135  }
1136 
1137  for (CDynamicCoverageEdge dce(FMap, edge); dce.cont(); ++dce) {
1138  for (int m = 0; m < NB_MARKS; m++)
1139  if (FCopyMarks[m] && FMap->isMarked(a1(*dce), m)) {
1140  FMap->setMark(*dce, m);
1141  }
1142  }
1143 
1144  EXIT;
1145 
1146  return edge;
1147 }
1148 
1149 //******************************************************************************
1150 
1151 CDart * CCorefine3d::splitEdge(CDart * AVertex, const CVertex & APoint)
1152 {
1153  ENTER;
1154 
1155  assert(FMap->isMarked(AVertex, FOrientMark));
1156 // assert(FTools.isPointOnLine(APoint, *getVertex(AVertex),
1157 // *getVertex(a0(AVertex))));
1158 
1159  MSG("Arête [" << *getVertex(AVertex) << "," << *getVertex(a0(AVertex))
1160  << "] éclatée en " << APoint);
1161 
1162  CDart *d, *face = getFace(AVertex);
1163 
1164  FMap->CGMapGeneric::insertVertex(AVertex);
1165 
1166  d = a1(a0(AVertex));
1167 
1168  FMap->halfMarkOrbit(d, ORBIT_VERTEX, FOrientMark);
1169 
1170  if (FMap->isMarked(AVertex, FFictiveMark))
1171  FMap->markOrbit(d, ORBIT_VERTEX, FFictiveMark);
1172 
1173  FMap->setVertex(d, APoint);
1174  updateVertexLinks(d);
1175 
1176  for (CDynamicCoverageVertex dcv(FMap, d); dcv.cont(); ++dcv) {
1177  FMap->setDirectInfo(*dcv, FFaceDI, face);
1178  FMap->setDirectInfo(*dcv, FAlpha2DI, NULL);
1179 
1180  for (int m = 0; m < NB_MARKS; m++)
1181  if (FCopyMarks[m] && FMap->isMarked(a0(*dcv), m)) {
1182  FMap->setMark(*dcv, m);
1183  }
1184  }
1185 
1186  EXIT;
1187 
1188  return d;
1189 }
1190 
1191 //******************************************************************************
1192 
1193 CDart * CCorefine3d::removeEdge(CDart * AEdge, int ADeleteMark)
1194 {
1195  CDart *d1 = a1(AEdge);
1196  CDart *d2 = a1(a0(AEdge));
1197 
1198  if (ADeleteMark < 0)
1199  FMap->merge(AEdge, a2(AEdge), 2, true);
1200  else {
1201  FMap->markOrbit(AEdge, ORBIT_EDGE, ADeleteMark);
1202  FMap->merge(AEdge, a2(AEdge), 2, false);
1203  }
1204 
1205  updateVertexLinks(d1);
1206  updateVertexLinks(d2);
1207 
1208  return d1;
1209 }
1210 
1211 //******************************************************************************
1212 
1213 CVertex CCorefine3d::getProjectionOnPlane(CDart * AVertex,
1214  const CPlane & APlane)
1215 {
1216  ENTER;
1217 
1218  CVertex pt, v1, v2;
1219  TCoordinate t=0;
1220 
1221  pt = *getVertex(AVertex);
1222  v1 = *getVertex(a0(AVertex)) - pt;
1223  v2 = *getVertex(a0(a1(AVertex))) - pt;
1224 
1225  APlane.getLineIntersection(pt, v1, &t);
1226 
1227  if (t > 0)
1228  return pt + t * v1;
1229 
1230  APlane.getLineIntersection(pt, v2, &t);
1231 
1232  EXIT;
1233 
1234  return pt + t * v2;
1235 }
1236 
1237 //******************************************************************************
1238 
1239 bool CCorefine3d::isSameEdge_Naive(CDart * AEdge1, CDart * AEdge2)
1240 {
1241  ENTER;
1242 
1243  assert(FMap->isMarked(AEdge1, FOrientMark));
1244  assert(FMap->isMarked(AEdge2, FOrientMark));
1245 
1246  bool found = false;
1247 
1248  for (CDynamicCoverage23 dc1(FMap, AEdge1); dc1.cont() && !found; ++dc1)
1249  for (CDynamicCoverage23 dc2(FMap, AEdge2); dc2.cont() && !found; ++dc2)
1250  found = (FMap->getDirectInfo(*dc1, FAlpha2DI) == *dc2);
1251 
1252  EXIT;
1253 
1254  return found;
1255 }
1256 
1257 //******************************************************************************
1258 
1259 bool CCorefine3d::isSameEdge_Optimized(CDart * AEdge1, CDart * AEdge2)
1260 {
1261  ENTER;
1262 
1263  assert(FMap->isMarked(AEdge1, FOrientMark));
1264  assert(FMap->isMarked(AEdge2, FOrientMark));
1265 
1266  CDart *d = AEdge1;
1267  bool same_edge = false;
1268 
1269  do {
1270  if (FMap->getDirectInfo(d, FAlpha2DI) != NULL)
1271  d = a3((CDart*)FMap->getDirectInfo(d, FAlpha2DI));
1272  else
1273  d = a3(FTools.alpha2(d));
1274 
1275  same_edge = (d == AEdge2) || (d == a3(a0(AEdge2)));
1276  }
1277  while (d != AEdge1 && !same_edge);
1278 
1279  EXIT;
1280 
1281  return same_edge;
1282 }
1283 
1284 //******************************************************************************
1285 #define isSameEdge isSameEdge_Optimized
1286 
1287 bool CCorefine3d::areFacesLinked(CDart * AFace1, CDart * AFace2)
1288 {
1289  assert(FMap->isMarked(AFace1, FOrientMark));
1290  assert(FMap->isMarked(AFace2, FOrientMark));
1291 
1292  return (isSameEdge(AFace1, AFace2) ||
1293  isSameEdge(AFace1, a3(a1(AFace2))) ||
1294  isSameEdge(a3(a1(AFace1)), AFace2) ||
1295  isSameEdge(a3(a1(AFace1)), a3(a1(AFace2))));
1296 }
1297 
1298 //******************************************************************************
1299 
1300 CEdgeIntersection
1301 CCorefine3d::findNearestIntersectionInGrid(const CVertex & AVertex1,
1302  const CVertex & AVertex2)
1303 {
1304  ENTER;
1305 
1306  CEdgeIntersection inter;
1307  TPositionOnEdge pos;
1308  TCoordinate t1 = 0.0, t2 = 0.0;
1309  CVertex *pt1, *pt2, pt;
1310  CPlane plane;
1311  CBoundingBox line_bb(AVertex1, AVertex2);
1312  CDart *face;
1313  int mark = FMap->getNewMark();
1314  int edge_mark = FMap->getNewMark();
1315  int face_mark = FMap->getNewMark();
1316 
1317  MSG("mark = " << mark);
1318  MSG("edge_mark = " << edge_mark);
1319  MSG("face_mark = " << face_mark);
1320 
1321  TFaceGridIter gi(*FGrid, line_bb);
1322  list<CDart*>::iterator li;
1323 
1324  inter.setCell(NULL);
1325  inter.setPosition(EP_OutOfEdge);
1326  inter.setParameter(1.0);
1327 
1328  MSG("Recherche d'une éventuelle intersection avec les sommets");
1329 
1330  for (gi.reinit(); gi.cont(); ++gi) {
1331  for (li = (*gi)->begin(); li != (*gi)->end(); ++li) {
1332  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont(); ++dcf) {
1333  if (!FMap->isMarked(*dcf, mark) && FMap->isMarked(*dcf, FOrientMark)) {
1334  // Marquage du sommet comme traité
1335  FMap->markOrbit(*dcf, ORBIT_VERTEX, mark);
1336 
1337  pt = *getVertex(*dcf);
1338 
1339  if (line_bb.isInBox(pt)) {
1340  pos = FTools.localizePointOnEdge(pt, AVertex1, AVertex2, &t2);
1341 
1342  if ((pos == EP_OnEdge && t2 < inter.getParameter()) ||
1343  (pos == EP_OnSecondVertex && t2 <= inter.getParameter())) {
1344  MSG("Une intersection trouvée avec le sommet " << pt);
1345 
1346  inter.setCell(*dcf);
1347  inter.setCellDimension(0);
1348  inter.setPosition(pos);
1349  inter.setPoint(pt);
1350  inter.setParameter(t2);
1351  }
1352 
1353  if (pos != EP_OutOfEdge) {
1354  for (CDynamicCoverageVertex dcv(FMap, *dcf); dcv.cont(); ++dcv) {
1355  // Marquage des arêtes incidentes pour ne pas les tester
1356  if (!FMap->isMarked(*dcv, edge_mark))
1357  FMap->markOrbit(*dcv, ORBIT_EDGE, edge_mark);
1358 
1359  // Marquage des faces incidentes pour ne pas les tester
1360  if (!FMap->isMarked(*dcv, face_mark))
1361  FMap->markOrbit(*dcv, ORBIT_FACE, face_mark);
1362  }
1363  }
1364  }
1365  }
1366  }
1367  }
1368  }
1369 
1370  MSG("Recherche d'une éventuelle intersection avec les arêtes");
1371 
1372  for (gi.reinit(); gi.cont(); ++gi) {
1373  for (li = (*gi)->begin(); li != (*gi)->end(); ++li) {
1374  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont(); ++dcf) {
1375  if (FMap->isMarked(*dcf, edge_mark)) {
1376  FMap->unmarkOrbit(*dcf, ORBIT_EDGE, edge_mark);
1377  FMap->unmarkOrbit(*dcf, ORBIT_EDGE, mark);
1378  }
1379  else if (FMap->isMarked(*dcf, mark) &&
1380 // !FMap->isMarked(*dcf, FFictiveMark) &&
1381  FMap->isMarked(*dcf, FOrientMark)) {
1382  FMap->unmarkOrbit(*dcf, ORBIT_EDGE, mark);
1383 
1384  pt1 = getVertex(*dcf);
1385  pt2 = getVertex(a0(*dcf));
1386 
1387  if (CBoundingBox(*pt1, *pt2) * line_bb) {
1388  pos = FTools.localizeEdgesIntersection(*pt1, *pt2,
1389  AVertex1, AVertex2,
1390  &t1, &t2);
1391 
1392  pt = *pt1 + t1 * (*pt2 - *pt1);
1393 
1394  if (/*FGrid->getCellBoundingBox(gi.getI(),
1395  gi.getJ(),
1396  gi.getK()).isInBox(pt) &&*/
1397  (t1 >= 0.0 && t1 <= 1.0)) {
1398  if ((pos == EP_OnEdge && t2 < inter.getParameter()) ||
1399  (pos == EP_OnSecondVertex && t2 <= inter.getParameter())) {
1400  MSG("Une intersection trouvée avec l'arête "
1401  << "[" << *pt1 << "," << *pt2 << "]");
1402 
1403  inter.setCell(*dcf);
1404  inter.setCellDimension(1);
1405  inter.setPosition(pos);
1406  inter.setPoint(pt);
1407  inter.setParameter(t2);
1408  }
1409 
1410  if (pos != EP_OutOfEdge) {
1411  for (CDynamicCoverageEdge dce(FMap, *dcf); dce.cont(); ++dce) {
1412  // Marquage des faces incidentes pour ne pas les tester
1413  if (!FMap->isMarked(*dce, face_mark))
1414  FMap->markOrbit(*dce, ORBIT_FACE, face_mark);
1415  }
1416  }
1417  }
1418  }
1419  }
1420  }
1421  }
1422  }
1423 
1424  MSG("Recherche d'une éventuelle intersection avec les faces");
1425 
1426  for (gi.reinit(); gi.cont(); ++gi) {
1427  for (li = (*gi)->begin(); li != (*gi)->end(); ++li) {
1428  if (FMap->isMarked(*li, face_mark)) {
1429  FMap->unmarkOrbit(*li, ORBIT_FACE, face_mark);
1430  FMap->markOrbit(*li, ORBIT_FACE, mark);
1431  }
1432  else if (!FMap->isMarked(*li, mark)) {
1433  face = (FMap->isMarked(*li, FOrientMark)) ? *li : a3(*li);
1434  FMap->markOrbit(face, ORBIT_FACE, mark);
1435 
1436 // if (FTools.orbitBoundingBox(face, ORBIT_01, FVertexDI) * line_bb) {
1437  plane = FTools.facePlane(face, FVertexDI);
1438 
1439  pos = FTools.localizeEdgeAndPlaneIntersection(AVertex1, AVertex2,
1440  plane, &t2);
1441 
1442  pt = AVertex1 + t2 * (AVertex2 - AVertex1);
1443 
1444  if (/*FGrid->getCellBoundingBox(gi.getI(),
1445  gi.getJ(),
1446  gi.getK()).isInBox(pt) &&*/
1447  FTools.isPointInFace(pt, face, &plane, FVertexDI) &&
1448  ((pos == EP_OnEdge && t2 < inter.getParameter()) ||
1449  (pos == EP_OnSecondVertex && t2 <= inter.getParameter()))) {
1450  MSG("Une intersection avec une face trouvée");
1451 
1452  inter.setCell(face);
1453  inter.setCellDimension(2);
1454  inter.setPosition(pos);
1455  inter.setPoint(pt);
1456  inter.setParameter(t2);
1457  }
1458 // }
1459  }
1460  }
1461  }
1462 
1463  MSG("Démarquage des faces");
1464 
1465  for (gi.reinit(); gi.cont(); ++gi)
1466  for (li = (*gi)->begin(); li != (*gi)->end(); ++li)
1467  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont(); ++dcf)
1468  if (FMap->isMarked(*dcf, mark))
1469  FMap->unmarkOrbit(*dcf, ORBIT_VERTEX, mark);
1470 
1471  FMap->freeMark(mark);
1472  FMap->freeMark(edge_mark);
1473  FMap->freeMark(face_mark);
1474 
1475  EXIT;
1476 
1477  return inter;
1478 }
1479 
1480 //******************************************************************************
1481 
1482 CDart *
1483 CCorefine3d::findFirstIntersectionInGrid(const CVertex & AVertex, int * ADim)
1484 {
1485  ENTER;
1486 
1487  CVertex pt1, pt2, pt;
1488  CDart *face, *result = NULL;
1489  CPlane plane;
1490  TCoordinate t;
1491  int mark = FMap->getNewMark();
1492 
1493  MSG("mark = " << mark);
1494 
1495  TFaceGridIter gi(*FGrid, CBoundingBox(AVertex));
1496  list<CDart*>::iterator li;
1497 
1498  // Recherche d'une éventuelle intersection avec les sommets
1499  for (gi.reinit(); gi.cont() && !result; ++gi)
1500  for (li = (*gi)->begin(); li != (*gi)->end() && !result; ++li)
1501  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont() && !result; ++dcf)
1502  if (!FMap->isMarked(*dcf, mark) && FMap->isMarked(*dcf, FOrientMark)) {
1503  // Marquage du sommet comme traité
1504  FMap->markOrbit(*dcf, ORBIT_VERTEX, mark);
1505 
1506  pt = *getVertex(*dcf);
1507 
1508  if (FTools.arePointsEqual(AVertex, pt)) {
1509  MSG("Une intersection trouvée avec le sommet " << pt);
1510  result = *dcf;
1511  *ADim = 0;
1512  }
1513  }
1514 
1515  // Recherche d'une éventuelle intersection avec les arêtes
1516  for (gi.reinit(); gi.cont() && !result; ++gi)
1517  for (li = (*gi)->begin(); li != (*gi)->end() && !result; ++li)
1518  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont() && !result; ++dcf)
1519  if (FMap->isMarked(*dcf, mark) && FMap->isMarked(*dcf, FOrientMark)) {
1520  FMap->unmarkOrbit(*dcf, ORBIT_EDGE, mark);
1521 
1522  pt1 = *getVertex(*dcf);
1523  pt2 = *getVertex(a0(*dcf));
1524 
1525  if (FTools.isPointOnLine(AVertex, pt1, pt2)) {
1526  t = FTools.pointParameterOnLine(AVertex, pt1, pt2);
1527 
1528  if (t >= 0.0 && t <= 1.0) {
1529  MSG("Une intersection trouvée avec l'arête "
1530  << "[" << pt1 << "," << pt2 << "]");
1531  result = *dcf;
1532  *ADim = 1;
1533  }
1534  }
1535  }
1536 
1537  // Recherche d'une éventuelle intersection avec les faces
1538  for (gi.reinit(); gi.cont() && !result; ++gi)
1539  for (li = (*gi)->begin(); li != (*gi)->end() && !result; ++li)
1540  if (!FMap->isMarked(*li, mark)) {
1541  face = (FMap->isMarked(*li, FOrientMark)) ? *li : a3(*li);
1542  FMap->markOrbit(face, ORBIT_FACE, mark);
1543 
1544  plane = FTools.facePlane(face, FVertexDI);
1545 
1546  if (FTools.isPointOnPlane(AVertex, plane) &&
1547  FTools.isPointInFace(AVertex, face, &plane, FVertexDI)) {
1548  MSG("Une intersection avec une face trouvée");
1549  result = face;
1550  *ADim = 2;
1551  }
1552  }
1553 
1554  for (gi.reinit(); gi.cont(); ++gi)
1555  for (li = (*gi)->begin(); li != (*gi)->end(); ++li)
1556  for (CDynamicCoverageFace dcf(FMap, *li); dcf.cont(); ++dcf)
1557  if (FMap->isMarked(*dcf, mark))
1558  FMap->unmarkOrbit(*dcf, ORBIT_VERTEX, mark);
1559 
1560  FMap->freeMark(mark);
1561 
1562  EXIT;
1563 
1564  return result;
1565 }
1566 
1567 //******************************************************************************
1568 
1569 void CCorefine3d::followIntersection(CDart * AVertex1, CDart * AVertex2,
1570  int AMark)
1571 {
1572  ENTER;
1573 
1574  /* Le marquage n'est pas encore au point car il teste plusieurs fois des
1575  * intersections déjà traitées. Cependant, l'algorithme fonctionne quand
1576  * même car il se contente de traiter l'intersection une deuxième fois comme
1577  * s'il s'agissait d'arêtes confondues...
1578  */
1579 
1580  assert(FMap->isMarked(AVertex1, FOrientMark));
1581  assert(FMap->isMarked(AVertex2, FOrientMark));
1582 
1583  list<CDart*> inter_list;
1584  CDart *d1, *d2, *current, *first, *last;
1585  bool ok;
1586  int mark = FMap->getNewMark();
1587  CTime start, end;
1588  int nb_vertices = 0;
1589  CPlane plane1, plane2;
1590 
1591  MSG("mark = " << mark);
1592 
1593  FNumberOfIntersectionLines++;
1594 
1595  cout << "--> Suivi d'une ligne de coupe" << endl;
1596  start.updateTime();
1597 
1598  SAVE_POINT(*getVertex(AVertex1));
1599 
1600  inter_list.push_back(AVertex1);
1601  inter_list.push_back(AVertex2);
1602  while (!inter_list.empty()) {
1603  d1 = inter_list.front(), inter_list.pop_front();
1604  d2 = inter_list.front(), inter_list.pop_front();
1605 
1606  nb_vertices++;
1607 
1608  for (CStaticCoverageVertex dcv1(FMap, d1); dcv1.cont(); ++dcv1)
1609  if (!FMap->isMarked(*dcv1, mark) &&
1610  FMap->isMarked(*dcv1, FOrientMark)) {
1611  FMap->markOrbit(*dcv1, ORBIT_13, mark);
1612 
1613  first = *dcv1;
1614  last = FTools.alpha2(a1(*dcv1));
1615 
1616  for (CStaticCoverageVertex dcv2(FMap, d2); dcv2.cont(); ++dcv2)
1617  if (!FMap->isMarked(*dcv2, mark) &&
1618  FMap->isMarked(*dcv2, FOrientMark)) {
1619  FMap->markOrbit(*dcv2, ORBIT_13, mark);
1620 
1621  /* Comme la première face peut être éclatée par chacune des
1622  * autres faces testées, il faut à chaque tour parcourir les
1623  * sous-faces de la première face
1624  */
1625  ok = false;
1626  current = first;
1627  do {
1628  plane1 = FTools.facePlane(current, FVertexDI);
1629  plane2 = FTools.facePlane(*dcv2, FVertexDI);
1630 
1631  if (!areFacesLinked(current, *dcv2)) {
1632  if ((plane1.getNormal() * plane2.getNormal()).norm() > FEps) {
1633  ok = manageFacesIntersection(current, plane1,
1634  *dcv2, plane2,
1635  AMark, &inter_list);
1636  }
1637  else {
1638  WARN_BB("<followIntersection> Faces coplanaires !");
1639  }
1640  }
1641 
1642  current = FTools.alpha2(a1(current));
1643  }
1644  while (current != last && !ok);
1645  }
1646 
1647  FMap->unmarkOrbit(d2, ORBIT_VERTEX, mark);
1648  }
1649 
1650  FMap->markOrbit(d1, ORBIT_VERTEX, AMark);
1651  FMap->unmarkOrbit(d1, ORBIT_VERTEX, mark);
1652  }
1653 
1654  FMap->freeMark(mark);
1655 
1656  end.updateTime();
1657 // cout << " Il y a eu " << nb_vertices << " étape(s) lors du suivi\n";
1658  cout << " Durée du suivi : " << end - start << endl;
1659  FLineCreationTime += end - start;
1660 
1661 // cout << " Mise à jour de la grille" << endl;
1662 // start.updateTime();
1663 // int nb_cells = refineGrid(MAX_GRID_RES / 16 , MAX_NUMBER_OF_FACES);
1664 // end.updateTime();
1665 // cout << " " << nb_cells << " mise(s) à jour effectuée(s) en "
1666 // << end - start << endl;
1667 
1668  EXIT;
1669 }
1670 
1671 //******************************************************************************
1672 
1673 bool CCorefine3d::manageFacesIntersection(CDart * AFace1,
1674  const CPlane & APlane1,
1675  CDart * AFace2,
1676  const CPlane & APlane2,
1677  int AMark,
1678  list<CDart*> * AList)
1679 {
1680  ENTER;
1681 
1682  assert(FMap->isMarked(AFace1, FOrientMark));
1683  assert(FMap->isMarked(AFace2, FOrientMark));
1684 
1685  CEdgeIntersection inter1, inter2;
1686  CVertex pt, dir;
1687  bool result = false;
1688 
1689 #ifdef DEBUG_MESSAGES
1690  FTools.displayFaceVertices(AFace1, FVertexDI);
1691  FTools.displayFaceVertices(AFace2, FVertexDI);
1692 #endif
1693 
1694  if (FTools.isVectorNull(APlane1.getNormal())) {
1695  cout << "\033[1;34m";
1696  FTools.displayFaceVertices(AFace1, FVertexDI);
1697  cout << "<manageFacesIntersection> La normale à la face 1 est nulle : "
1698  << APlane1.getNormal() << "\033[0m\n";
1699  }
1700  if (FTools.isVectorNull(APlane2.getNormal())) {
1701  cout << "\033[1;34m";
1702  FTools.displayFaceVertices(AFace2, FVertexDI);
1703  cout << "<manageFacesIntersection> La normale à la face 2 est nulle : "
1704  << APlane2.getNormal() << "\033[0m\n";
1705  }
1706 
1707  pt = *getVertex(AFace1);
1708 
1709  // Calcul de la ligne de coupe
1710  dir = APlane1.getNormal() * APlane2.getNormal();
1711 
1712 // dir = FTools.normalizeVector(dir);
1713 
1714  if (!FTools.isVectorInSector(dir, AFace1, APlane1, true, FVertexDI) ||
1715  !FTools.isVectorInSector(dir, AFace2, APlane2, true, FVertexDI))
1716  dir = -dir;
1717 
1718  if (FTools.isVectorInSector(dir, AFace1, APlane1, true, FVertexDI) &&
1719  FTools.isVectorInSector(dir, AFace2, APlane2, true, FVertexDI)) {
1720  MSG("Ligne de coupe = (" << pt << "," << pt + dir << ")");
1721 
1722  // On cherche la plus proche intersection dans chaque face
1723  inter1 = FTools.findNearestIntersectionInFace(pt, dir, AFace1, APlane1,
1724  true, FVertexDI);
1725  inter2 = FTools.findNearestIntersectionInFace(pt, dir, AFace2, APlane2,
1726  true, FVertexDI);
1727 
1728  if (inter1.getCell() != AFace1 && inter2.getCell() != AFace2) {
1729  if (inter1.getCell() != NULL && inter2.getCell() != NULL) {
1730  createIntersectionEdge(AFace1, AFace2, APlane1, APlane2,
1731  inter1, inter2, AMark, AList);
1732  result = true;
1733  }
1734  }
1735  else {
1736  WARN_BB("<manageFacesIntersection> Intersections trouvées invalides !");
1737  }
1738  }
1739 
1740  EXIT;
1741 
1742  return result;
1743 }
1744 
1745 //******************************************************************************
1746 #define sortFacesAroundEdges sortFacesAroundEdges_SuperNaive
1747 
1748 void CCorefine3d::createIntersectionEdge(CDart * AFace1, CDart * AFace2,
1749  const CPlane & APlane1,
1750  const CPlane & APlane2,
1751  const CEdgeIntersection & AInter1,
1752  const CEdgeIntersection & AInter2,
1753  int AMark, list<CDart*> * AList)
1754 {
1755  ENTER;
1756 
1757  CDart *d1, *d2, *edge1, *edge2;
1758  bool add_to_list = true;
1759  CVertex pt;
1760 
1761  removeFaceFromGrid(FGrid, AFace2);
1762 
1763  // L'intersection se trouve à la même distance sur chacune des faces
1764  if (FTools.arePointsEqual(AInter1.getPoint(), AInter2.getPoint())) {
1765  MSG("L'intersection sur chaque face se trouve à égale distance");
1766 
1767  pt = (AInter1.getPoint() + AInter2.getPoint()) / 2.0;
1768 
1769  if (AInter1.getCellDimension() == 0) {
1770  d1 = AInter1.getCell();
1771  *getVertex(AInter1.getCell()) = pt;
1772  }
1773  else
1774  d1 = splitEdge(AInter1.getCell(), pt);
1775 
1776  if (AInter2.getCellDimension() == 0) {
1777  d2 = AInter2.getCell();
1778  *getVertex(AInter2.getCell()) = pt;
1779  }
1780  else
1781  d2 = splitEdge(AInter2.getCell(), pt);
1782 
1783  if (FMap->isMarked(d1, AMark))
1784  add_to_list = false;
1785 
1786  edge1 = splitFace(AFace1, d1);
1787  edge2 = splitFace(AFace2, d2);
1788 
1789 // for (CDynamicCoverage23 dc(FMap, a0(edge1)); dc.cont(); ++dc) {
1790 // FMap->setMark(*dc, AMark);
1791 // FMap->setMark(a1(*dc), AMark);
1792 // FMap->setMark(a0(*dc), AMark);
1793 // FMap->setMark(a1(a0(*dc)), AMark);
1794 // }
1795 
1796  SAVE_POINT(pt);
1797  }
1798  // L'intersection la plus proche se trouve sur la même face
1799  else if (AInter1.getParameter() < AInter2.getParameter()) {
1800  MSG("L'intersection la plus proche se trouve sur la même face");
1801 
1802  if (AInter1.getCellDimension() == 0)
1803  d1 = AInter1.getCell();
1804  else
1805  d1 = splitEdge(AInter1.getCell(), AInter1.getPoint());
1806 
1807  pt = getProjectionOnPlane(d1, APlane2);
1808  *getVertex(d1) = pt;
1809 
1810  if (FMap->isMarked(d1, AMark))
1811  add_to_list = false;
1812 
1813  edge1 = splitFace(AFace1, d1);
1814  edge2 = insertEdgeInFace(AFace2, pt);
1815 
1816  d2 = a1(a0(edge2));
1817 
1818 // for (CDynamicCoverage23 dc(FMap, a0(edge1)); dc.cont(); ++dc) {
1819 // FMap->setMark(*dc, AMark);
1820 // FMap->setMark(a1(*dc), AMark);
1821 // FMap->setMark(a0(*dc), AMark);
1822 // FMap->setMark(a1(a0(*dc)), AMark);
1823 // }
1824 
1825  SAVE_POINT(pt);
1826  }
1827  // L'intersection la plus proche se trouve sur l'autre face
1828  else { // (AInter2.getParameter() < AInter1.getParameter())
1829  MSG("L'intersection la plus proche se trouve sur l'autre face");
1830 
1831  if (AInter2.getCellDimension() == 0)
1832  d2 = AInter2.getCell();
1833  else
1834  d2 = splitEdge(AInter2.getCell(), AInter2.getPoint());
1835 
1836  pt = getProjectionOnPlane(d2, APlane1);
1837  *getVertex(d2) = pt;
1838 
1839  edge1 = insertEdgeInFace(AFace1, pt);
1840  edge2 = splitFace(AFace2, d2);
1841 
1842  d1 = a1(a0(edge1));
1843 
1844 // for (CDynamicCoverage23 dc(FMap, edge1); dc.cont(); ++dc) {
1845 // FMap->setMark(*dc, AMark);
1846 // FMap->setMark(a1(*dc), AMark);
1847 // }
1848 
1849  SAVE_POINT(pt);
1850  }
1851 
1852  FMap->markOrbit(a0(edge1), ORBIT_VERTEX, AMark);
1853 
1854  assert(FMap->isMarked(edge1, FOrientMark));
1855  assert(FMap->isMarked(edge2, FOrientMark));
1856 
1857  if (!isSameEdge(edge1, edge2)) {
1858  FNumberOfIntersectionEdges++;
1859 
1860  MSG("Tri des faces autour de l'arête d'intersection n°"
1861  << FNumberOfIntersectionEdges);
1862 
1863  FTools.sortFacesAroundEdges(edge1, edge2, FAlpha2DI, FVertexDI);
1864 
1865 #ifdef EXTRACT_LINES //---------------------------------------------------------
1866  FMap->markOrbit(edge1, ORBIT_EDGE, FIntersectionMark);
1867  FMap->markOrbit(edge2, ORBIT_EDGE, FIntersectionMark);
1868 #endif // EXTRACT_LINES --------------------------------------------------------
1869  }
1870 
1871  if (getFace(edge2) == getFace(a2(edge2))) {
1872  addFaceToGrid(FGrid, edge2);
1873  }
1874  else {
1875  addFaceToGrid(FGrid, edge2);
1876  addFaceToGrid(FGrid, a2(edge2));
1877  }
1878 
1879  if (add_to_list) {
1880  AList->push_back(d1);
1881  AList->push_back(d2);
1882  }
1883 
1884 #ifdef FOLLOW_DEBUG_MODE //-----------------------------------------------------
1885  CGMapVertex * tmp_map = new CGMapVertex;
1886  char file_name[1024] = "\0";
1887 
1888  stringstream s;
1889  s<<"edge";
1890  s.width(4); s.fill(0); s.flags (ios::right);
1891  s<<FNumberOfIntersectionEdges;
1892  s.width(0);
1893  s<<".map";
1894 
1895  d1 = tmp_map->addMapDart();
1896  d2 = tmp_map->addMapDart();
1897 
1898  tmp_map->linkAlpha0(d1, d2);
1899 
1900  tmp_map->setVertex(d1, *getVertex(edge1));
1901  tmp_map->setVertex(d2, *getVertex(a0(edge1)));
1902 
1903  tmp_map->save(s.str(), AsciiFormat);
1904 
1905  delete tmp_map;
1906 #endif // FOLLOW_DEBUG_MODE ----------------------------------------------------
1907 
1908  EXIT;
1909 }
1910 
1911 #undef sortFacesAroundEdges
1912 //******************************************************************************
1913 
1914 void CCorefine3d::applyModifications(CDart * AMesh)
1915 {
1916  ENTER;
1917 
1918  CDart *d;
1919 
1920  for (CStaticCoverageCC scc(FMap, AMesh); scc.cont(); ++scc) {
1921  d = (CDart*)FMap->getDirectInfo(*scc, FAlpha2DI);
1922 
1923  if (d != NULL) {
1924  if (!FMap->isFree2(*scc))
1925  FMap->unsew2(*scc);
1926  if (!FMap->isFree2(d))
1927  FMap->unsew2(d);
1928 
1929  FMap->sew2(*scc, d);
1930 
1931  FMap->setDirectInfo(*scc, FAlpha2DI, NULL);
1932  FMap->setDirectInfo(a0(*scc), FAlpha2DI, NULL);
1933  FMap->setDirectInfo(d, FAlpha2DI, NULL);
1934  FMap->setDirectInfo(a0(d), FAlpha2DI, NULL);
1935  }
1936  }
1937 
1938  EXIT;
1939 }
1940 
1941 //******************************************************************************
1942 #define COPY(d) ((CDart*)FMap->getDirectInfo(d, copy_direct_info))
1943 
1944 void CCorefine3d::extractIntersectionLines(CDart * AMesh)
1945 {
1946  ENTER;
1947 
1948  int copy_direct_info = FMap->getNewDirectInfo();
1949  CDart *d;
1950  CDynamicCoverageCC dcc(FMap, AMesh);
1951 
1952  // Duplication des brins appartenant aux lignes de coupe
1953  for (dcc.reinit(); dcc.cont(); ++dcc) {
1954  if (FMap->isMarked(*dcc, FIntersectionMark)) {
1955  FMap->setDirectInfo(*dcc, copy_direct_info, FMap->addMapDart());
1956  FMap->setMark(COPY(*dcc), FIntersectionMark);
1957 // FMap->setVertex(COPY(*dcc), *FMap->findVertex(*dcc));
1958  }
1959  }
1960 
1961  // Créations des liaisons entre ces brins
1962  for (dcc.reinit(); dcc.cont(); ++dcc) {
1963  if (FMap->isMarked(*dcc, FIntersectionMark)) {
1964  if (FMap->canSew0(COPY(*dcc), COPY(a0(*dcc))))
1965  FMap->sew0(COPY(*dcc), COPY(a0(*dcc)));
1966 
1967  d = a1(*dcc);
1968  while (!FMap->isMarked(d, FIntersectionMark))
1969  d = a1(a2(d));
1970 
1971  if (FMap->canSew1(COPY(*dcc), COPY(d)))
1972  FMap->sew1(COPY(*dcc), COPY(d));
1973  }
1974  }
1975 
1976  // Créations des liaisons entre ces brins
1977  for (dcc.reinit(); dcc.cont(); ++dcc) {
1978  if (FMap->isMarked(*dcc, FIntersectionMark)) {
1979  if (FMap->canSew2(COPY(*dcc), COPY(a2(*dcc))))
1980  FMap->sew2(COPY(*dcc), COPY(a2(*dcc)));
1981  if (FMap->canSew3(COPY(*dcc), COPY(a3(*dcc))))
1982  FMap->sew3(COPY(*dcc), COPY(a3(*dcc)));
1983  }
1984  }
1985 
1986  // Copie des plongements
1987  for (dcc.reinit(); dcc.cont(); ++dcc) {
1988  if (FMap->isMarked(COPY(*dcc), FIntersectionMark)&&
1989  FMap->isMarked(COPY(*dcc), FIntersectionMark)) {
1990  FMap->unmarkOrbit(COPY(*dcc), ORBIT_VERTEX, FIntersectionMark);
1991  FMap->setVertex(COPY(*dcc), *FMap->findVertex(*dcc));
1992  }
1993  }
1994 
1995  FMap->freeDirectInfo(copy_direct_info);
1996 
1997  EXIT;
1998 }
1999 
2000 #undef COPY
2001 //******************************************************************************
2002 
2003 #undef a0
2004 #undef a1
2005 #undef a2
2006 #undef a3
2007 
2008 //******************************************************************************