diff --git a/include/Algo/ImplicitHierarchicalMesh/ihm3.h b/include/Algo/ImplicitHierarchicalMesh/ihm3.h index 933889674d636d7e86b45ea9d0bdd0da0c608133..643687675b1a728d4617658186ed2b6fcc40f805 100644 --- a/include/Algo/ImplicitHierarchicalMesh/ihm3.h +++ b/include/Algo/ImplicitHierarchicalMesh/ihm3.h @@ -264,6 +264,10 @@ public: */ bool volumeIsSubdividedOnce(Dart d); + /** + * + */ + bool neighborhoodLevelDiffersByOne(Dart d); } ; diff --git a/include/Algo/ImplicitHierarchicalMesh/subdivision3.h b/include/Algo/ImplicitHierarchicalMesh/subdivision3.h index 40de7c90a8e7b31cbb2c200aa39c81176d3e1161..f12fdfcd4dc0580bcd62d27f435ce2667084449c 100644 --- a/include/Algo/ImplicitHierarchicalMesh/subdivision3.h +++ b/include/Algo/ImplicitHierarchicalMesh/subdivision3.h @@ -45,13 +45,19 @@ template void subdivideEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) ; template -void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType = S_QUAD); +void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); template Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); template -void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); +Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); + +template +Dart subdivideVolumeOld(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); + +//template +//void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); template void coarsenEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position); diff --git a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp index 6b1feed3c029691dacc2006764eb51ff53398b58..725b10c99528a4d58f65f02f65c93d48c013d9b4 100644 --- a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp +++ b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp @@ -24,6 +24,7 @@ #include "Algo/Geometry/centroid.h" #include "Algo/Modelisation/subdivision.h" +#include "Algo/Modelisation/extrusion.h" namespace CGoGN { @@ -66,7 +67,7 @@ void subdivideEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position } template -void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType) +void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) { assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; assert(!map.faceIsSubdivided(d) || !"Trying to subdivide an already subdivided face") ; @@ -94,7 +95,7 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position Dart res; - if(degree == 3 && sType == S_TRI) //subdiviser une face triangulaire + if(degree == 3) //subdiviser une face triangulaire { Dart dd = map.phi1(old) ; Dart e = map.phi1(map.phi1(dd)) ; @@ -168,7 +169,7 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position } template -Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) +Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) { assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; @@ -179,7 +180,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi unsigned int cur = map.getCurrentLevel(); map.setCurrentLevel(vLevel); - /* * au niveau du volume courant i * stockage d'un brin de chaque face de celui-ci @@ -235,10 +235,13 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi /* * Subdivision */ - //Store the darts from quadrangulated faces - std::vector > subdividedfaces; - subdividedfaces.reserve(25); + std::vector > subdividedfacesQ; + subdividedfacesQ.reserve(25); + + std::vector > subdividedfacesT; + subdividedfacesT.reserve(25); + //First step : subdivide edges and faces //creates a i+1 edge level and i+1 face level @@ -250,6 +253,7 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi if(!map.faceIsSubdivided(d)) Algo::IHM::subdivideFace(map, d, position); + //save a dart from the subdivided face unsigned int cur = map.getCurrentLevel() ; @@ -257,14 +261,32 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi map.setCurrentLevel(fLevel) ; - //le brin est forcement du niveau cur - Dart cf = map.phi1(d); - Dart e = cf; - do + + //test si la face est triangulaire ou non + if(map.phi1(map.phi1(map.phi1(d))) == d) { - subdividedfaces.push_back(std::pair(e,map.phi2(e))); - e = map.phi2(map.phi1(e)); - }while (e != cf); + std::cout << "trian" << std::endl; + Dart cf = map.phi2(map.phi1(d)); + Dart e = cf; + do + { + subdividedfacesT.push_back(std::pair(e,map.phi2(e))); + e = map.phi1(e); + }while (e != cf); + } + else + { + std::cout << "quad" << std::endl; + Dart cf = map.phi1(d); + Dart e = cf; + do + { + subdividedfacesQ.push_back(std::pair(e,map.phi2(e))); + e = map.phi2(map.phi1(e)); + }while (e != cf); + + + } map.setCurrentLevel(cur); } @@ -277,17 +299,21 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi //Second step : deconnect each corner, close each hole, subdivide each new face into 3 for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) { + //std::vector::iterator edge = oldEdges.begin(); Dart e = *edge; Dart f1 = map.phi1(*edge); + Dart f2 = map.phi2(f1); do { - map.unsewFaces(map.phi1(map.phi1(e))); + if(map.phi1(map.phi1(map.phi1(e))) != e) + { + map.unsewFaces(map.phi1(map.phi1(e))); //remplacer par une boucle qui découd toute la face et non juste une face carre (jusqu'a phi_1(e)) + } + + map.unsewFaces(map.phi1(e)); - //TODO utile ? - //if(map.phi2(map.phi1(e)) != map.phi1(e)) - map.unsewFaces(map.phi1(e)); e = map.phi2(map.phi_1(e)); } @@ -295,89 +321,347 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi map.closeHole(f1); - Dart old = map.phi2(map.phi1(e)); - Dart dd = map.phi1(map.phi1(old)) ; - map.splitFace(old,dd) ; + //degree du sommet exterieur + unsigned int cornerDegree = map.Map2::vertexDegree(*edge); - Dart ne = map.phi1(map.phi1(old)) ; + //tourner autour du sommet pour connaitre le brin d'un sommet de valence < cornerDegree + bool found = false; + Dart stop = e; + do + { - map.cutEdge(ne); - position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? - newEdges.push_back(ne); - newEdges.push_back(map.phi1(ne)); + if(map.Map2::vertexDegree(map.phi2(map.phi1(e))) < cornerDegree) + { + stop = map.phi2(map.phi1(e)); + found = true; + } + e = map.phi2(map.phi_1(e)); + } + while(!found && e != *edge); - Dart stop = map.phi2(map.phi1(ne)); - ne = map.phi2(ne); - do + //si il existe un sommet de degre inferieur au degree du coin + if(found) { - dd = map.phi1(map.phi1(map.phi1(ne))); + //chercher le brin de faible degree suivant + bool found2 = false; + Dart dd = map.phi1(stop); - //A Verifier !! - map.splitFace(ne, dd) ; + do + { + if(map.Map2::vertexDegree(dd) < cornerDegree) + found2 = true; + else + dd = map.phi1(dd); + } + while(!found2); - newEdges.push_back(map.phi1(dd)); + //cas de la pyramide + if(dd == stop) + { + std::cout << "pyramide" << std::endl; + map.splitFace(dd, map.phi1(map.phi1(dd))); + } + else + { + map.splitFace(dd, stop); - ne = map.phi2(map.phi_1(ne)); - dd = map.phi1(map.phi1(dd)); - } - while(dd != stop); - } + //calcul de la taille des faces de chaque cote de stop + if(!( (map.Map2::faceDegree(map.phi_1(stop)) == 3 && map.Map2::faceDegree(map.phi2(map.phi_1(stop))) == 4) || + (map.Map2::faceDegree(map.phi_1(stop)) == 4 && map.Map2::faceDegree(map.phi2(map.phi_1(stop))) == 3) )) + { + std::cout << "octaedre ou hexaedre" << std::endl; + + Dart ne = map.phi_1(stop) ; + map.cutEdge(ne); + position[map.phi1(ne)] = volCenter; + stop = map.phi2(map.phi1(ne)); + + bool finished = false; + Dart it = map.phi2(ne); + + do + { + //chercher le brin de faible degree suivant + bool found2 = false; + Dart dd = map.phi1(it); + + do + { + if(dd == stop) + finished = true; + else if(map.Map2::vertexDegree(dd) < cornerDegree) + found2 = true; + else + dd = map.phi1(dd); + } + while(!found2 & !finished); + + if(found2) + { + map.splitFace(it,dd); + } + + it = map.phi_1(dd); + + if(it == stop) + finished = true; + + } + while(!finished); - //Third step : 3-sew internal faces - for (std::vector >::iterator it = subdividedfaces.begin(); it != subdividedfaces.end(); ++it) - { - Dart f1 = (*it).first; - Dart f2 = (*it).second; + } + else + { + std::cout << "prisme" << std::endl; + //tester si besoin de fermer f2 (par exemple pas besoin pour hexa... mais pour tet, octa, prisme oui) + //map.closeHole(f2); + } - if(map.phi3(map.phi2(f1)) == map.phi2(f1) && map.phi3(map.phi2(f2)) == map.phi2(f2)) - { - //id pour toutes les faces interieures - map.sewVolumes(map.phi2(f1), map.phi2(f2)); + } - //Fais a la couture !!!!! - unsigned int idface = map.getNewFaceId(); - map.setFaceId(map.phi2(f1),idface, FACE); + } + //sinon cas du tetraedre + else + { + std::cout << "tetraedre" << std::endl; + //tester si besoin de fermer f2 (par exemple pas besoin pour hexa... mais pour tet, octa, prisme oui) + //map.closeHole(f2); } - //FAIS a la couture !!!!!!! - //id pour toutes les aretes exterieurs des faces quadrangulees - unsigned int idedge = map.getEdgeId(f1); - map.setEdgeId(map.phi2(f1), idedge, DART); - map.setEdgeId( map.phi2(f2), idedge, DART); } - //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! - //id pour les aretes interieurs : (i.e. 6 pour un hexa) - DartMarker mne(map); - for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) + std::cout << "1ere etape finished" << std::endl; + + CellMarker mtf(map, FACE); + + //Etape 2 + for (std::vector >::iterator edges = subdividedfacesT.begin(); edges != subdividedfacesT.end(); ++edges) { - if(!mne.isMarked(*it)) + Dart f1 = (*edges).first; + Dart f2 = (*edges).second; + +//Fonction isTetrahedron ?? +// //if(Algo::Modelisation::Tetrahedron::isTetrahedron(map,f2)) + if( (map.Map2::faceDegree(f2) == 3 && map.Map2::faceDegree(map.phi2(f2)) == 3 && + map.Map2::faceDegree(map.phi2(map.phi_1(f2))) == 3) && map.Map2::vertexDegree(f2) == 3) + { //cas du tetrahedre + + std::cout << "ajout d'une face" << std::endl; + + if(map.phi3(map.phi2(f2)) == map.phi2(f2)) + { + Dart nf = map.newFace(3); + map.sewVolumes(map.phi2(f2),nf); + } + + if(map.phi2(map.phi3(map.phi2(f2))) == map.phi3(map.phi2(f2))) + { + map.sewFaces(map.phi3(map.phi2(f2)), f1); + } + } + else { - unsigned int idedge = map.getNewEdgeId(); - map.setEdgeId(*it, idedge, EDGE); - mne.markOrbit(EDGE,*it); + if(!mtf.isMarked(f1)) + { + mtf.mark(f1); + + map.closeHole(f1); + + if(map.Map2::faceDegree(map.phi2(f2)) == 3) + { + std::cout << "ajout d'un tetraedre" << std::endl; + Dart x = Algo::Modelisation::trianguleFace(map, map.phi2(f1)); + position[x] = volCenter; + } + else + { + std::cout << "ajout d'un prisme" << std::endl; + //Dart x = Algo::Modelisation::extrudeFace(map,position,map.phi2(f1),5.0); + Dart c = Algo::Modelisation::trianguleFace(map, map.phi2(f1)); + + Dart cc = c; + // cut edges + do + { + + typename PFP::VEC3 p1 = position[cc] ; + typename PFP::VEC3 p2 = position[map.phi1(cc)] ; + + map.cutEdge(cc); + + position[map.phi1(cc)] = (p1 + p2) * typename PFP::REAL(0.5) ; + + cc = map.phi2(map.phi_1(cc)); + }while (cc != c); + + // cut faces + do + { + Dart d1 = map.phi1(cc); + Dart d2 = map.phi_1(cc); + map.splitFace(d1,d2); + cc = map.phi2(map.phi_1(cc));//map.Map2::alpha1(cc); + }while (cc != c); + + //merge central faces by removing edges + bool notFinished=true; + do + { + Dart d1 = map.Map2::alpha1(cc); + if (d1 == cc) // last edge is pending edge inside of face + notFinished = false; + map.deleteFace(cc); + cc = d1; + } while (notFinished); + + + map.closeHole(map.phi1(map.phi1(map.phi2(f1)))); + + } + } } + } + std::cout << "2e etape finished" << std::endl; + + +// { +// //Third step : 3-sew internal faces +// for (std::vector >::iterator it = subdividedfacesT.begin(); it != subdividedfacesT.end(); ++it) +// { +// Dart f1 = (*it).first; +// Dart f2 = (*it).second; +// +// +// +// if(map.phi3(map.phi2(f1)) == map.phi2(f1) && map.phi3(map.phi2(f2)) == map.phi2(f2)) +// { +// if(map.getEmbedding(VERTEX, map.phi_1(map.phi2(f2))) == map.getEmbedding(VERTEX, map.phi_1(map.phi2(f1)))) +// { +// map.Map3::sewVolumes(map.phi2(f2), map.phi2(f1)); +// } +// else +// { +// +// //id pour toutes les faces interieures +// map.sewVolumes(map.phi2(f2), map.phi2(f1)); +// +// +// } +// +// //Fais a la couture !!!!! +// unsigned int idface = map.getNewFaceId(); +// map.setFaceId(map.phi2(f1),idface, FACE); +// } +// +// +// //FAIS a la couture !!!!!!! +// //id pour toutes les aretes exterieurs des faces quadrangulees +// unsigned int idedge = map.getEdgeId(f1); +// map.setEdgeId(map.phi2(f1), idedge, DART); +// map.setEdgeId( map.phi2(f2), idedge, DART); +// +// } +// +// //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! +// //id pour les aretes interieurs : (i.e. 16 pour un octa) +// DartMarker mne(map); +// for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) +// { +// if(!mne.isMarked(*it)) +// { +// unsigned int idedge = map.getNewEdgeId(); +// map.setEdgeId(*it, idedge, EDGE); +// mne.markOrbit(EDGE,*it); +// } +// } +// } +// +// { +// //Third step : 3-sew internal faces +// for (std::vector >::iterator it = subdividedfacesQ.begin(); it != subdividedfacesQ.end(); ++it) +// { +// Dart f1 = (*it).first; +// Dart f2 = (*it).second; +// +// if(map.phi3(map.phi2(f1)) == map.phi2(f1) && map.phi3(map.phi2(f2)) == map.phi2(f2)) +// { +// //id pour toutes les faces interieures +// map.sewVolumes(map.phi2(f2), map.phi2(f1)); +// +// //Fais a la couture !!!!! +// unsigned int idface = map.getNewFaceId(); +// map.setFaceId(map.phi2(f1),idface, FACE); +// } +// +// //FAIS a la couture !!!!!!! +// //id pour toutes les aretes exterieurs des faces quadrangulees +// unsigned int idedge = map.getEdgeId(f1); +// map.setEdgeId(map.phi2(f1), idedge, DART); +// map.setEdgeId( map.phi2(f2), idedge, DART); +// } +// //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! +// //id pour les aretes interieurs : (i.e. 16 pour un octa) +// DartMarker mne(map); +// for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) +// { +// if(!mne.isMarked(*it)) +// { +// unsigned int idedge = map.getNewEdgeId(); +// map.setEdgeId(*it, idedge, EDGE); +// mne.markOrbit(EDGE,*it); +// } +// } +// } +// +// +// //cas tordu pour le prisme +// for (std::vector >::iterator it = subdividedfacesT.begin(); it != subdividedfacesT.end(); ++it) +// { +// Dart f1 = (*it).first; +// Dart f2 = (*it).second; +// +// if( !(map.Map2::faceDegree(f2) == 3 && map.Map2::faceDegree(map.phi2(f2)) == 3 && +// map.Map2::faceDegree(map.phi2(map.phi1(f2))) == 3 && map.Map2::faceDegree(map.phi2(map.phi_1(f2))) == 3)) +// { +// +// +// if(map.phi2(map.phi1(map.phi1(map.phi2(f1)))) == map.phi3(map.phi2(map.phi1(map.phi1(map.phi2(f1))))) && +// map.phi2(map.phi3(map.phi2(map.phi3(map.phi2(map.phi3(map.phi2(map.phi2(map.phi1(map.phi1(map.phi2(f1))))))))))) == +// map.phi3(map.phi2(map.phi3(map.phi2(map.phi3(map.phi2(map.phi3(map.phi2(map.phi2(map.phi1(map.phi1(map.phi2(f1)))))))))))) +// ) +// { +// map.sewVolumes(map.phi2(map.phi1(map.phi1(map.phi2(f1)))), +// map.phi2(map.phi3(map.phi2(map.phi3(map.phi2(map.phi3(map.phi1(map.phi1(map.phi2(f1)))))))))); +// } +// +// } +// +// } + + map.setCurrentLevel(cur) ; - return subdividedfaces.begin()->first; + return subdividedfacesQ.begin()->first; } + + template -void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) +Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) { assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; + assert(!map.neighborhoodLevelDiffersByOne(d) || !"Trying to subdivide a volume with neighborhood Level difference greater than 1"); - unsigned int vLevel = map.volumeLevel(d) ; - Dart old = map.volumeOldestDart(d) ; - - unsigned int cur = map.getCurrentLevel() ; - map.setCurrentLevel(vLevel) ; // go to the level of the face to subdivide its edges + unsigned int vLevel = map.volumeLevel(d); + Dart old = map.volumeOldestDart(d); + unsigned int cur = map.getCurrentLevel(); + map.setCurrentLevel(vLevel); /* * au niveau du volume courant i @@ -386,7 +670,7 @@ void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position */ DartMarkerStore mf(map); // Lock a face marker to save one dart per face - DartMarkerStore mv(map); // Lock a vertex marker to compute volume center + CellMarker mv(map, VERTEX); typename PFP::VEC3 volCenter; unsigned count = 0 ; @@ -411,9 +695,9 @@ void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position //compute volume centroid if(!mv.isMarked(e)) { + mv.mark(e); volCenter += position[e]; ++count; - mv.markOrbit(VERTEX, e); oldEdges.push_back(e); } @@ -429,168 +713,698 @@ void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position } while(e != *face) ; } - volCenter /= typename PFP::REAL(count) ; + volCenter /= typename PFP::REAL(count) ; /* * Subdivision */ //Store the darts from quadrangulated faces - std::vector > subdividedfaces; - subdividedfaces.reserve(25); + std::vector > subdividedfacesQ; + subdividedfacesQ.reserve(25); + + std::vector > subdividedfacesT; + subdividedfacesT.reserve(25); + + + //First step : subdivide edges and faces + //creates a i+1 edge level and i+1 face level + for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) + { + Dart d = *face; + + //if needed subdivide face + if(!map.faceIsSubdivided(d)) + Algo::IHM::subdivideFace(map, d, position); - //First step : subdivide edges and faces - //creates a i+1 edge level and i+1 face level - for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) - { - Dart d = *face; - //if needed subdivide face - if(!map.faceIsSubdivided(d)) - Algo::IHM::subdivideFace(map, d, position, Algo::IHM::S_TRI); + //save a dart from the subdivided face + unsigned int cur = map.getCurrentLevel() ; + + unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee + map.setCurrentLevel(fLevel) ; - //save a dart from the subdivided face - unsigned int cur = map.getCurrentLevel() ; - unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee - map.setCurrentLevel(fLevel) ; - //le brin est forcement du niveau cur + //test si la face est triangulaire ou non + if(map.phi1(map.phi1(map.phi1(d))) == d) + { + std::cout << "trian" << std::endl; Dart cf = map.phi2(map.phi1(d)); Dart e = cf; do { - subdividedfaces.push_back(std::pair(e,map.phi2(e))); + subdividedfacesT.push_back(std::pair(e,map.phi2(e))); e = map.phi1(e); }while (e != cf); - - map.setCurrentLevel(cur); } - - map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision - - std::vector newEdges; //save darts from inner edges - newEdges.reserve(50); - - bool istet = true; - - //Second step : deconnect each corner, close each hole, subdivide each new face into 3 - for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) + else { - Dart e = *edge; - - Dart f1 = map.phi1(*edge); - Dart f2 = map.phi2(f1); - + std::cout << "quad" << std::endl; + Dart cf = map.phi1(d); + Dart e = cf; do { + subdividedfacesQ.push_back(std::pair(e,map.phi2(e))); + e = map.phi2(map.phi1(e)); + }while (e != cf); - map.unsewFaces(map.phi1(e)); - e = map.phi2(map.phi_1(e)); - } - while(e != *edge); + } - map.closeHole(f1); - map.closeHole(f2); - map.sewVolumes(map.phi2(f1),map.phi2(f2)); + map.setCurrentLevel(cur); - unsigned int fdeg = map.faceDegree(map.phi2(f1)); + } - if(fdeg > 3) - { - istet = false; - Dart old = map.phi2(map.phi1(*edge)); - Dart dd = map.phi1(old) ; - map.splitFace(old,dd) ; + map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision - Dart ne = map.phi1(old); + std::vector newEdges; //save darts from inner edges + newEdges.reserve(50); - map.cutEdge(ne); - position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? - newEdges.push_back(ne); - newEdges.push_back(map.phi1(ne)); + bool istet = true; - Dart stop = map.phi2(map.phi1(ne)); - ne = map.phi2(ne); - do - { - dd = map.phi1(map.phi1(ne)); + //Second step : deconnect each corner, close each hole, subdivide each new face into 3 + for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) + { + Dart e = *edge; - map.splitFace(ne, dd) ; + Dart f1 = map.phi1(*edge); + Dart f2 = map.phi2(f1); - newEdges.push_back(map.phi1(dd)); + do + { + if(map.phi1(map.phi1(map.phi1(e))) != e) + map.unsewFaces(map.phi1(map.phi1(e))); - ne = map.phi2(map.phi_1(ne)); - dd = map.phi1(dd); - } - while(dd != stop); + map.unsewFaces(map.phi1(e)); + + + e = map.phi2(map.phi_1(e)); + } + while(e != *edge); + + map.closeHole(f1); + + + //fonction qui calcule le degree max des faces atour d'un sommet + unsigned int fdeg = map.faceDegree(map.phi2(f1)); + + + if(fdeg > 4) + { + Dart old = map.phi2(map.phi1(e)); + Dart dd = map.phi1(map.phi1(old)) ; + map.splitFace(old,dd) ; + + Dart ne = map.phi1(map.phi1(old)) ; + + map.cutEdge(ne); + position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? + newEdges.push_back(ne); + newEdges.push_back(map.phi1(ne)); + + + Dart stop = map.phi2(map.phi1(ne)); + ne = map.phi2(ne); + do + { + dd = map.phi1(map.phi1(map.phi1(ne))); + + //A Verifier !! + map.splitFace(ne, dd) ; + + newEdges.push_back(map.phi1(dd)); + + ne = map.phi2(map.phi_1(ne)); + dd = map.phi1(map.phi1(dd)); } - else + while(dd != stop); + } + else if(fdeg > 3) + { + map.closeHole(f2); + map.sewVolumes(map.phi2(f1),map.phi2(f2)); + + istet = false; + Dart old = map.phi2(map.phi1(*edge)); + Dart dd = map.phi1(old) ; + map.splitFace(old,dd) ; + + Dart ne = map.phi1(old); + + map.cutEdge(ne); + position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? + newEdges.push_back(ne); + newEdges.push_back(map.phi1(ne)); + + Dart stop = map.phi2(map.phi1(ne)); + ne = map.phi2(ne); + do { - unsigned int idface = map.getNewFaceId(); - map.setFaceId(map.phi2(f1),idface, FACE); + dd = map.phi1(map.phi1(ne)); + + map.splitFace(ne, dd) ; + + newEdges.push_back(map.phi1(dd)); + + ne = map.phi2(map.phi_1(ne)); + dd = map.phi1(dd); } + while(dd != stop); + } + else + { + map.closeHole(f2); + map.sewVolumes(map.phi2(f1),map.phi2(f2)); + unsigned int idface = map.getNewFaceId(); + map.setFaceId(map.phi2(f1),idface, FACE); } - if(!istet) + } + + if(!istet) + { + for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) { - for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) + Dart e = *edge; + + Dart x = map.phi_1(map.phi2(map.phi1(*edge))); + Dart f = x; + + do { - Dart e = *edge; + Dart f3 = map.phi3(f); + Dart tmp = map.phi_1(map.phi2(map.phi_1(map.phi2(map.phi_1(f3))))); - Dart x = map.phi_1(map.phi2(map.phi1(*edge))); - Dart f = x; + map.unsewFaces(f3); + map.unsewFaces(tmp); + map.sewFaces(f3, tmp); - do - { - Dart f3 = map.phi3(f); - Dart tmp = map.phi_1(map.phi2(map.phi_1(map.phi2(map.phi_1(f3))))); + unsigned int idface = map.getNewFaceId(); + map.setFaceId(map.phi2(f3),idface, FACE); - map.unsewFaces(f3); - map.unsewFaces(tmp); - map.sewFaces(f3, tmp); - unsigned int idface = map.getNewFaceId(); - map.setFaceId(map.phi2(f3),idface, FACE); + f = map.phi2(map.phi_1(f)); + }while(f != x); + } - f = map.phi2(map.phi_1(f)); - }while(f != x); + } - } + { + //Third step : 3-sew internal faces + for (std::vector >::iterator it = subdividedfacesT.begin(); it != subdividedfacesT.end(); ++it) + { + Dart f1 = (*it).first; + Dart f2 = (*it).second; + + //FAIS a la couture !!!!!!! + //id pour toutes les aretes exterieurs des faces quadrangulees + unsigned int idedge = map.getEdgeId(f1); + map.setEdgeId(map.phi2(f1), idedge, DART); + map.setEdgeId( map.phi2(f2), idedge, DART); + + } + + //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! + //id pour les aretes interieurs : (i.e. 16 pour un octa) + DartMarker mne(map); + for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) + { + if(!mne.isMarked(*it)) + { + unsigned int idedge = map.getNewEdgeId(); + map.setEdgeId(*it, idedge, EDGE); + mne.markOrbit(EDGE,*it); } + } + } + + { + //Third step : 3-sew internal faces + for (std::vector >::iterator it = subdividedfacesQ.begin(); it != subdividedfacesQ.end(); ++it) + { + Dart f1 = (*it).first; + Dart f2 = (*it).second; - //Third step : 3-sew internal faces - for (std::vector >::iterator it = subdividedfaces.begin(); it != subdividedfaces.end(); ++it) + if(map.phi3(map.phi2(f1)) == map.phi2(f1) && map.phi3(map.phi2(f2)) == map.phi2(f2)) { - Dart f1 = (*it).first; - Dart f2 = (*it).second; + //id pour toutes les faces interieures + map.sewVolumes(map.phi2(f1), map.phi2(f2)); + + //Fais a la couture !!!!! + unsigned int idface = map.getNewFaceId(); + map.setFaceId(map.phi2(f1),idface, FACE); + } - //FAIS a la couture !!!!!!! - //id pour toutes les aretes exterieurs des faces quadrangulees - unsigned int idedge = map.getEdgeId(f1); - map.setEdgeId(map.phi2(f1), idedge, DART); - map.setEdgeId( map.phi2(f2), idedge, DART); + //FAIS a la couture !!!!!!! + //id pour toutes les aretes exterieurs des faces quadrangulees + unsigned int idedge = map.getEdgeId(f1); + map.setEdgeId(map.phi2(f1), idedge, DART); + map.setEdgeId( map.phi2(f2), idedge, DART); + } + //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! + //id pour les aretes interieurs : (i.e. 6 pour un hexa) + DartMarker mne(map); + for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) + { + if(!mne.isMarked(*it)) + { + unsigned int idedge = map.getNewEdgeId(); + map.setEdgeId(*it, idedge, EDGE); + mne.markOrbit(EDGE,*it); } + } + } + + map.setCurrentLevel(cur) ; + + return subdividedfacesQ.begin()->first; +} + + +template +Dart subdivideVolumeOld(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) +{ + assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; + assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; - //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! - //id pour les aretes interieurs : (i.e. 16 pour un octa) - DartMarker mne(map); - for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) + unsigned int vLevel = map.volumeLevel(d); + Dart old = map.volumeOldestDart(d); + + unsigned int cur = map.getCurrentLevel(); + map.setCurrentLevel(vLevel); + + + /* + * au niveau du volume courant i + * stockage d'un brin de chaque face de celui-ci + * avec calcul du centroid + */ + + DartMarkerStore mf(map); // Lock a face marker to save one dart per face + CellMarker mv(map, VERTEX); + + typename PFP::VEC3 volCenter; + unsigned count = 0 ; + + //Store faces that are traversed and start with the face of d + std::vector visitedFaces; + visitedFaces.reserve(20); + visitedFaces.push_back(old); + + //Store the edges before the cutEdge + std::vector oldEdges; + oldEdges.reserve(20); + + mf.markOrbit(FACE, old) ; + + for(std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) + { + Dart e = *face ; + do { - if(!mne.isMarked(*it)) + //add one old edge per vertex to the old edge list + //compute volume centroid + if(!mv.isMarked(e)) { - unsigned int idedge = map.getNewEdgeId(); - map.setEdgeId(*it, idedge, EDGE); - mne.markOrbit(EDGE,*it); + mv.mark(e); + volCenter += position[e]; + ++count; + oldEdges.push_back(e); } + + // add all face neighbours to the table + Dart ee = map.phi2(e) ; + if(!mf.isMarked(ee)) // not already marked + { + visitedFaces.push_back(ee) ; + mf.markOrbit(FACE, ee) ; + } + + e = map.phi1(e) ; + } while(e != *face) ; + } + + volCenter /= typename PFP::REAL(count) ; + + /* + * Subdivision + */ + + //Store the darts from quadrangulated faces + std::vector > subdividedfaces; + subdividedfaces.reserve(25); + + //First step : subdivide edges and faces + //creates a i+1 edge level and i+1 face level + for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) + { + Dart d = *face; + + //if needed subdivide face + if(!map.faceIsSubdivided(d)) + Algo::IHM::subdivideFace(map, d, position); + + //save a dart from the subdivided face + unsigned int cur = map.getCurrentLevel() ; + + unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee + map.setCurrentLevel(fLevel) ; + + + //le brin est forcement du niveau cur + Dart cf = map.phi1(d); + Dart e = cf; + do + { + subdividedfaces.push_back(std::pair(e,map.phi2(e))); + e = map.phi2(map.phi1(e)); + }while (e != cf); + + map.setCurrentLevel(cur); + } + + map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision + + std::vector newEdges; //save darts from inner edges + newEdges.reserve(50); + + //Second step : deconnect each corner, close each hole, subdivide each new face into 3 + for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) + { + Dart e = *edge; + + Dart f1 = map.phi1(*edge); + + do + { + map.unsewFaces(map.phi1(map.phi1(e))); + + map.unsewFaces(map.phi1(e)); + + e = map.phi2(map.phi_1(e)); } + while(e != *edge); + map.closeHole(f1); - map.setCurrentLevel(cur) ; + Dart old = map.phi2(map.phi1(e)); + Dart dd = map.phi1(map.phi1(old)) ; + map.splitFace(old,dd) ; + + Dart ne = map.phi1(map.phi1(old)) ; + + map.cutEdge(ne); + position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? + newEdges.push_back(ne); + newEdges.push_back(map.phi1(ne)); + + + Dart stop = map.phi2(map.phi1(ne)); + ne = map.phi2(ne); + do + { + dd = map.phi1(map.phi1(map.phi1(ne))); + + //A Verifier !! + map.splitFace(ne, dd) ; + + newEdges.push_back(map.phi1(dd)); + + ne = map.phi2(map.phi_1(ne)); + dd = map.phi1(map.phi1(dd)); + } + while(dd != stop); + } + + //Third step : 3-sew internal faces + for (std::vector >::iterator it = subdividedfaces.begin(); it != subdividedfaces.end(); ++it) + { + Dart f1 = (*it).first; + Dart f2 = (*it).second; + + if(map.phi3(map.phi2(f1)) == map.phi2(f1) && map.phi3(map.phi2(f2)) == map.phi2(f2)) + { + //id pour toutes les faces interieures + map.sewVolumes(map.phi2(f1), map.phi2(f2)); + + //Fais a la couture !!!!! + unsigned int idface = map.getNewFaceId(); + map.setFaceId(map.phi2(f1),idface, FACE); + } + + //FAIS a la couture !!!!!!! + //id pour toutes les aretes exterieurs des faces quadrangulees + unsigned int idedge = map.getEdgeId(f1); + map.setEdgeId(map.phi2(f1), idedge, DART); + map.setEdgeId( map.phi2(f2), idedge, DART); + } + + //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! + //id pour les aretes interieurs : (i.e. 6 pour un hexa) + DartMarker mne(map); + for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) + { + if(!mne.isMarked(*it)) + { + unsigned int idedge = map.getNewEdgeId(); + map.setEdgeId(*it, idedge, EDGE); + mne.markOrbit(EDGE,*it); + } + } + + map.setCurrentLevel(cur) ; + + return subdividedfaces.begin()->first; } +// +//template +//void subdivideLoop(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) +//{ +// assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; +// assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; +// +// unsigned int vLevel = map.volumeLevel(d) ; +// Dart old = map.volumeOldestDart(d) ; +// +// unsigned int cur = map.getCurrentLevel() ; +// map.setCurrentLevel(vLevel) ; // go to the level of the face to subdivide its edges +// +// +// /* +// * au niveau du volume courant i +// * stockage d'un brin de chaque face de celui-ci +// * avec calcul du centroid +// */ +// +// DartMarkerStore mf(map); // Lock a face marker to save one dart per face +// DartMarkerStore mv(map); // Lock a vertex marker to compute volume center +// +// typename PFP::VEC3 volCenter; +// unsigned count = 0 ; +// +// //Store faces that are traversed and start with the face of d +// std::vector visitedFaces; +// visitedFaces.reserve(20); +// visitedFaces.push_back(old); +// +// //Store the edges before the cutEdge +// std::vector oldEdges; +// oldEdges.reserve(20); +// +// mf.markOrbit(FACE, old) ; +// +// for(std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) +// { +// Dart e = *face ; +// do +// { +// //add one old edge per vertex to the old edge list +// //compute volume centroid +// if(!mv.isMarked(e)) +// { +// volCenter += position[e]; +// ++count; +// mv.markOrbit(VERTEX, e); +// oldEdges.push_back(e); +// } +// +// // add all face neighbours to the table +// Dart ee = map.phi2(e) ; +// if(!mf.isMarked(ee)) // not already marked +// { +// visitedFaces.push_back(ee) ; +// mf.markOrbit(FACE, ee) ; +// } +// +// e = map.phi1(e) ; +// } while(e != *face) ; +// } +// +// volCenter /= typename PFP::REAL(count) ; +// +// /* +// * Subdivision +// */ +// //Store the darts from quadrangulated faces +// std::vector > subdividedfaces; +// subdividedfaces.reserve(25); +// +// //First step : subdivide edges and faces +// //creates a i+1 edge level and i+1 face level +// for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) +// { +// Dart d = *face; +// +// //if needed subdivide face +// if(!map.faceIsSubdivided(d)) +// Algo::IHM::subdivideFace(map, d, position, Algo::IHM::S_TRI); +// +// //save a dart from the subdivided face +// unsigned int cur = map.getCurrentLevel() ; +// +// unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee +// map.setCurrentLevel(fLevel) ; +// +// //le brin est forcement du niveau cur +// Dart cf = map.phi2(map.phi1(d)); +// Dart e = cf; +// do +// { +// subdividedfaces.push_back(std::pair(e,map.phi2(e))); +// e = map.phi1(e); +// }while (e != cf); +// +// map.setCurrentLevel(cur); +// } +// +// map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision +// +// std::vector newEdges; //save darts from inner edges +// newEdges.reserve(50); +// +// bool istet = true; +// +// //Second step : deconnect each corner, close each hole, subdivide each new face into 3 +// for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) +// { +// Dart e = *edge; +// +// Dart f1 = map.phi1(*edge); +// Dart f2 = map.phi2(f1); +// +// do +// { +// +// map.unsewFaces(map.phi1(e)); +// +// e = map.phi2(map.phi_1(e)); +// } +// while(e != *edge); +// +// map.closeHole(f1); +// map.closeHole(f2); +// map.sewVolumes(map.phi2(f1),map.phi2(f2)); +// +// unsigned int fdeg = map.faceDegree(map.phi2(f1)); +// +// if(fdeg > 3) +// { +// istet = false; +// Dart old = map.phi2(map.phi1(*edge)); +// Dart dd = map.phi1(old) ; +// map.splitFace(old,dd) ; +// +// Dart ne = map.phi1(old); +// +// map.cutEdge(ne); +// position[map.phi1(ne)] = volCenter; //plonger a la fin de la boucle ???? +// newEdges.push_back(ne); +// newEdges.push_back(map.phi1(ne)); +// +// Dart stop = map.phi2(map.phi1(ne)); +// ne = map.phi2(ne); +// do +// { +// dd = map.phi1(map.phi1(ne)); +// +// map.splitFace(ne, dd) ; +// +// newEdges.push_back(map.phi1(dd)); +// +// ne = map.phi2(map.phi_1(ne)); +// dd = map.phi1(dd); +// } +// while(dd != stop); +// } +// else +// { +// unsigned int idface = map.getNewFaceId(); +// map.setFaceId(map.phi2(f1),idface, FACE); +// } +// +// } +// +// if(!istet) +// { +// for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) +// { +// Dart e = *edge; +// +// Dart x = map.phi_1(map.phi2(map.phi1(*edge))); +// Dart f = x; +// +// do +// { +// Dart f3 = map.phi3(f); +// Dart tmp = map.phi_1(map.phi2(map.phi_1(map.phi2(map.phi_1(f3))))); +// +// map.unsewFaces(f3); +// map.unsewFaces(tmp); +// map.sewFaces(f3, tmp); +// +// unsigned int idface = map.getNewFaceId(); +// map.setFaceId(map.phi2(f3),idface, FACE); +// +// +// f = map.phi2(map.phi_1(f)); +// }while(f != x); +// +// } +// } +// +// //Third step : 3-sew internal faces +// for (std::vector >::iterator it = subdividedfaces.begin(); it != subdividedfaces.end(); ++it) +// { +// Dart f1 = (*it).first; +// Dart f2 = (*it).second; +// +// //FAIS a la couture !!!!!!! +// //id pour toutes les aretes exterieurs des faces quadrangulees +// unsigned int idedge = map.getEdgeId(f1); +// map.setEdgeId(map.phi2(f1), idedge, DART); +// map.setEdgeId( map.phi2(f2), idedge, DART); +// +// } +// +// //LA copie de L'id est a gerer avec le sewVolumes normalement !!!!!! +// //id pour les aretes interieurs : (i.e. 16 pour un octa) +// DartMarker mne(map); +// for(std::vector::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it) +// { +// if(!mne.isMarked(*it)) +// { +// unsigned int idedge = map.getNewEdgeId(); +// map.setEdgeId(*it, idedge, EDGE); +// mne.markOrbit(EDGE,*it); +// } +// } +// +// +// map.setCurrentLevel(cur) ; +//} template @@ -811,173 +1625,182 @@ void splitVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) map.setCurrentLevel(cur) ; } -/************************************************************************************* - * - */ - -template -void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType) -{ - assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; - assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; - - unsigned int vLevel = map.volumeLevel(d) ; - Dart old = map.volumeOldestDart(d) ; - - unsigned int cur = map.getCurrentLevel() ; - map.setCurrentLevel(vLevel) ; // go to the level of the face to subdivide its edges - - - /* - * au niveau du volume courant i - * stockage d'un brin de chaque face de celui-ci - * avec calcul du centroid - */ - - DartMarkerStore mf(map); // Lock a face marker to save one dart per face - DartMarkerStore mv(map); // Lock a vertex marker to compute volume center - - typename PFP::VEC3 volCenter; - unsigned count = 0 ; - - //Store faces that are traversed and start with the face of d - std::vector visitedFaces; - visitedFaces.reserve(20); - visitedFaces.push_back(old); - - //Store the edges before the cutEdge - std::vector oldEdges; - oldEdges.reserve(20); - - mf.markOrbit(FACE, old) ; - - for(std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) - { - Dart e = *face ; - do - { - //add one old edge per vertex to the old edge list - //compute volume centroid - if(!mv.isMarked(e)) - { - volCenter += position[e]; - ++count; - mv.markOrbit(VERTEX, e); - oldEdges.push_back(e); - } - - // add all face neighbours to the table - Dart ee = map.phi2(e) ; - if(!mf.isMarked(ee)) // not already marked - { - visitedFaces.push_back(ee) ; - mf.markOrbit(FACE, ee) ; - } - - e = map.phi1(e) ; - } while(e != *face) ; - } - - volCenter /= typename PFP::REAL(count) ; - - /* - * Subdivision - */ - - //Store the darts from quadrangulated faces - std::vector > subdividedTriFaces; - subdividedTriFaces.reserve(15); - - std::vector > subdividedQuadFaces; - subdividedQuadFaces.reserve(15); - - //First step : subdivide edges and faces - //creates a i+1 edge level and i+1 face level - for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) - { - Dart d = *face; - - //if needed subdivide face - if(!map.faceIsSubdivided(d)) - Algo::IHM::subdivideFace(map, d, position, Algo::IHM::S_TRI); - - //save a dart from the subdivided face - unsigned int cur = map.getCurrentLevel() ; - - unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee - map.setCurrentLevel(fLevel) ; - - //Face tri - if(map.phi1(map.phi1(map.phi1(d))) == d) - { - - //le brin est forcement du niveau cur - Dart cf = map.phi2(map.phi1(d)); - Dart e = cf; - do - { - subdividedTriFaces.push_back(std::pair(e,map.phi2(e))); - e = map.phi1(e); - }while (e != cf); - } - else - { - //le brin est forcement du niveau cur - Dart cf = map.phi1(d); - Dart e = cf; - do - { - subdividedQuadFaces.push_back(std::pair(e,map.phi2(e))); - e = map.phi2(map.phi1(e)); - }while (e != cf); - } - - map.setCurrentLevel(cur); - } - - map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision - //Second step : create corners - for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) - { - Dart e = *edge; - Dart f1 = map.phi1(*edge); - Dart f2 = map.phi2(map.phi1(*edge)); - do - { - if(map.phi1(map.phi1(map.phi1(e))) != e) - map.unsewFaces(map.phi1(map.phi1(e))); - if(map.phi2(map.phi1(e)) != map.phi1(e)) - map.unsewFaces(map.phi1(e)); - e = map.phi2(map.phi_1(e)); - }while(e != *edge); - map.closeHole(f1); - map.closeHole(f2); - map.sewVolumes(map.phi2(f1),map.phi2(f2)); - //un passage sur chaque nouvelle face pour recopier l'id d'arete - Dart d = map.phi2(f1); - e = d; - do - { - unsigned int idedge = map.getEdgeId(map.phi2(e)); - map.setEdgeId(e, idedge, EDGE); - e = map.phi1(e); - }while(e != d); - //Identifiant de face pour les 2 faces orientees - unsigned int idface = map.getNewFaceId(); - map.setFaceId(map.phi2(f1),idface, FACE); - } +/************************************************************************************* + * + */ - std::vector inFaces; - inFaces.reserve(10); +//template +//void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position, SubdivideType sType) +//{ +// assert(map.getDartLevel(d) <= map.getCurrentLevel() || !"Access to a dart introduced after current level") ; +// assert(!map.volumeIsSubdivided(d) || !"Trying to subdivide an already subdivided volume") ; +// +// unsigned int vLevel = map.volumeLevel(d) ; +// Dart old = map.volumeOldestDart(d) ; +// +// unsigned int cur = map.getCurrentLevel() ; +// map.setCurrentLevel(vLevel) ; // go to the level of the face to subdivide its edges +// +// +// /* +// * au niveau du volume courant i +// * stockage d'un brin de chaque face de celui-ci +// * avec calcul du centroid +// */ +// +// DartMarkerStore mf(map); // Lock a face marker to save one dart per face +// DartMarkerStore mv(map); // Lock a vertex marker to compute volume center +// +// typename PFP::VEC3 volCenter; +// unsigned count = 0 ; +// +// //Store faces that are traversed and start with the face of d +// std::vector visitedFaces; +// visitedFaces.reserve(20); +// visitedFaces.push_back(old); +// +// //Store the edges before the cutEdge +// std::vector oldEdges; +// oldEdges.reserve(20); +// +// mf.markOrbit(FACE, old) ; +// +// for(std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) +// { +// Dart e = *face ; +// do +// { +// //add one old edge per vertex to the old edge list +// //compute volume centroid +// if(!mv.isMarked(e)) +// { +// volCenter += position[e]; +// ++count; +// mv.markOrbit(VERTEX, e); +// oldEdges.push_back(e); +// } +// +// // add all face neighbours to the table +// Dart ee = map.phi2(e) ; +// if(!mf.isMarked(ee)) // not already marked +// { +// visitedFaces.push_back(ee) ; +// mf.markOrbit(FACE, ee) ; +// } +// +// e = map.phi1(e) ; +// } while(e != *face) ; +// } +// +// volCenter /= typename PFP::REAL(count) ; +// +// /* +// * Subdivision +// */ +// +// //Store the darts from quadrangulated faces +// std::vector > subdividedTriFaces; +// subdividedTriFaces.reserve(15); +// +// std::vector > subdividedQuadFaces; +// subdividedQuadFaces.reserve(15); +// +// //First step : subdivide edges and faces +// //creates a i+1 edge level and i+1 face level +// for (std::vector::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face) +// { +// Dart d = *face; +// +// //if needed subdivide face +// if(!map.faceIsSubdivided(d)) +// Algo::IHM::subdivideFaceQuad(map, d, position); +// +// //save a dart from the subdivided face +// unsigned int cur = map.getCurrentLevel() ; +// +// unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee +// map.setCurrentLevel(fLevel) ; +// +// //Face tri +// if(map.phi1(map.phi1(map.phi1(d))) == d) +// { +// +// //le brin est forcement du niveau cur +// Dart cf = map.phi2(map.phi1(d)); +// Dart e = cf; +// do +// { +// subdividedTriFaces.push_back(std::pair(e,map.phi2(e))); +// e = map.phi1(e); +// }while (e != cf); +// } +// else +// { +// //le brin est forcement du niveau cur +// Dart cf = map.phi1(d); +// Dart e = cf; +// do +// { +// subdividedQuadFaces.push_back(std::pair(e,map.phi2(e))); +// e = map.phi2(map.phi1(e)); +// }while (e != cf); +// } +// +// map.setCurrentLevel(cur); +// } +// +// map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision +// +// //Second step : create corners +// for (std::vector::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge) +// { +// Dart e = *edge; +// +// Dart f1 = map.phi1(*edge); +// Dart f2 = map.phi2(map.phi1(*edge)); +// +// do +// { +// if(map.phi1(map.phi1(map.phi1(e))) != e) +// map.unsewFaces(map.phi1(map.phi1(e))); +// if(map.phi2(map.phi1(e)) != map.phi1(e)) +// map.unsewFaces(map.phi1(e)); +// +// e = map.phi2(map.phi_1(e)); +// +// }while(e != *edge); +// +// map.closeHole(f1); +// map.closeHole(f2); +// map.sewVolumes(map.phi2(f1),map.phi2(f2)); +// +// //un passage sur chaque nouvelle face pour recopier l'id d'arete +// Dart d = map.phi2(f1); +// e = d; +// do +// { +// unsigned int idedge = map.getEdgeId(map.phi2(e)); +// map.setEdgeId(e, idedge, EDGE); +// e = map.phi1(e); +// }while(e != d); +// +// +// //Identifiant de face pour les 2 faces orientees +// unsigned int idface = map.getNewFaceId(); +// map.setFaceId(map.phi2(f1),idface, FACE); +// } +// +// std::vector inFaces; +// inFaces.reserve(10); //DartMarkerStore mv(map); @@ -1039,8 +1862,8 @@ void subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi // // } - map.setCurrentLevel(cur); -} +// map.setCurrentLevel(cur); +//} // Dart old = map.phi1(map.phi2(map.phi1(e))); diff --git a/src/Algo/ImplicitHierarchicalMesh/ihm3.cpp b/src/Algo/ImplicitHierarchicalMesh/ihm3.cpp index 74fe96da1a7f29cd60c48399df2b5dcbf7db9f9e..10944fbf8402437aedfbffe6a2391790741924aa 100644 --- a/src/Algo/ImplicitHierarchicalMesh/ihm3.cpp +++ b/src/Algo/ImplicitHierarchicalMesh/ihm3.cpp @@ -515,6 +515,51 @@ bool ImplicitHierarchicalMap3:: volumeIsSubdividedOnce(Dart d) return true; } +bool ImplicitHierarchicalMap3::neighborhoodLevelDiffersByOne(Dart d) +{ + assert(m_dartLevel[d] <= m_curLevel || !"Access to a dart introduced after current level") ; + + bool found = false; + + unsigned int vLevel = volumeLevel(d) + 1; + + Dart old = volumeOldestDart(d); + + DartMarkerStore mf(*this); // Lock a face marker to save one dart per face + + //Store faces that are traversed and start with the face of d + std::vector visitedFaces; + visitedFaces.reserve(20); + visitedFaces.push_back(old); + + mf.markOrbit(FACE, old) ; + + for(std::vector::iterator face = visitedFaces.begin(); !found && face != visitedFaces.end(); ++face) + { + Dart e = *face ; + do + { + // add all face neighbours to the table + + if(phi3(e) != e && (abs(volumeLevel(phi3(e)) - vLevel) > 1)) + { + found = true; + } + + Dart ee = phi2(e) ; + if(!mf.isMarked(ee)) // not already marked + { + visitedFaces.push_back(ee) ; + mf.markOrbit(FACE, ee) ; + } + + e = phi1(e) ; + } while(e != *face) ; + } + + return found; +} + } //namespace IHM } //namespace Algo