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 1f09bc895c469dd886698794d80e4fbb2f71472c..f12fdfcd4dc0580bcd62d27f435ce2667084449c 100644 --- a/include/Algo/ImplicitHierarchicalMesh/subdivision3.h +++ b/include/Algo/ImplicitHierarchicalMesh/subdivision3.h @@ -53,6 +53,9 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi template 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); diff --git a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp index 60175bf9577ef830a5309689000b945fedee8857..725b10c99528a4d58f65f02f65c93d48c013d9b4 100644 --- a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp +++ b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp @@ -168,9 +168,6 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position map.setCurrentLevel(cur) ; } - - - template Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) { @@ -364,10 +361,6 @@ Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& pos { std::cout << "pyramide" << std::endl; map.splitFace(dd, map.phi1(map.phi1(dd))); - //tester si besoin de fermer f2 (par exemple pas besoin pour hexa... mais pour tet, octa, prisme oui) - //map.closeHole(f2); - - //retenir un truc pour le sew/unsew final de l'octa au milieu } else { @@ -452,7 +445,7 @@ Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& pos // //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; @@ -534,108 +527,119 @@ Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& pos 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)) - { - //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); - } - } - } - - - 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.phi3(map.phi2(map.phi3(map.phi2(map.phi3(map.phi1(map.phi1(map.phi2(f1))))))))), - map.phi2(map.phi1(map.phi1(map.phi2(f1))))); - } - - } - - } +// { +// //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) ; @@ -651,6 +655,7 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi { 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); @@ -800,6 +805,7 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi map.closeHole(f1); + //fonction qui calcule le degree max des faces atour d'un sommet unsigned int fdeg = map.faceDegree(map.phi2(f1)); @@ -868,7 +874,7 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi else { map.closeHole(f2); - //map.sewVolumes(map.phi2(f1),map.phi2(f2)); + map.sewVolumes(map.phi2(f1),map.phi2(f2)); unsigned int idface = map.getNewFaceId(); map.setFaceId(map.phi2(f1),idface, FACE); @@ -978,203 +984,201 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi } -//template -//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") ; -// -// 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 -// { -// //add one old edge per vertex to the old edge list -// //compute volume centroid -// if(!mv.isMarked(e)) -// { -// 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))); -// -// //TODO utile ? -// //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); -// -// 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 +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") ; + + 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 + { + //add one old edge per vertex to the old edge list + //compute volume centroid + if(!mv.isMarked(e)) + { + 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); + + 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) @@ -1634,169 +1638,169 @@ void splitVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position) * */ -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); @@ -1858,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