Skip to content

C
CGoGN
Commit 79edfd10 authored by pitiot's avatar pitiot
Browse files
Options

adding debug

parent 22edb204
Hide whitespace changes
Inline Side-by-side
Showing with 396 additions and 396 deletions
include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp
View file @ 79edfd10
......@@ -2576,463 +2576,463 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, VertexAttribute<typename PF
// }
template <typename PFP>
Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& 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<typename PFP::MAP> mf(map); // Lock a face marker to save one dart per face
CellMarker<typename PFP::MAP, VERTEX> mv(map);
typename PFP::VEC3 volCenter;
unsigned count = 0 ;
//Store faces that are traversed and start with the face of d
std::vector<Dart> visitedFaces;
visitedFaces.reserve(512);
visitedFaces.push_back(old);
//Store the edges before the cutEdge
std::vector<Dart> oldEdges;
oldEdges.reserve(512);
mf.markOrbit<FACE>(old) ;
for(unsigned int i = 0; i < visitedFaces.size(); ++i)
{
Dart e = visitedFaces[i] ;
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);
}
//template <typename PFP>
//Dart subdivideVolumeGen(typename PFP::MAP& map, Dart d, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& 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") ;
// 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 != visitedFaces[i]) ;
}
// unsigned int vLevel = map.volumeLevel(d);
// Dart old = map.volumeOldestDart(d);
volCenter /= typename PFP::REAL(count) ;
// unsigned int cur = map.getCurrentLevel();
// map.setCurrentLevel(vLevel);
/*
* Subdivision
*/
//Store the darts from quadrangulated faces
std::vector<std::pair<Dart,Dart> > subdividedfacesQ;
subdividedfacesQ.reserve(25);
// /*
// * au niveau du volume courant i
// * stockage d'un brin de chaque face de celui-ci
// * avec calcul du centroid
// */
std::vector<std::pair<Dart,Dart> > subdividedfacesT;
subdividedfacesT.reserve(25);
// DartMarkerStore<typename PFP::MAP> mf(map); // Lock a face marker to save one dart per face
// CellMarker<typename PFP::MAP, VERTEX> mv(map);
// typename PFP::VEC3 volCenter;
// unsigned count = 0 ;
//First step : subdivide edges and faces
//creates a i+1 edge level and i+1 face level
for (std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
{
Dart d = *face;
// //Store faces that are traversed and start with the face of d
// std::vector<Dart> visitedFaces;
// visitedFaces.reserve(512);
// visitedFaces.push_back(old);
//if needed subdivide face
if(!map.faceIsSubdivided(d))
IHM::subdivideFace<PFP>(map, d, position);
// //Store the edges before the cutEdge
// std::vector<Dart> oldEdges;
// oldEdges.reserve(512);
// mf.markOrbit<FACE>(old) ;
//save a dart from the subdivided face
unsigned int cur = map.getCurrentLevel() ;
// for(unsigned int i = 0; i < visitedFaces.size(); ++i)
// {
// Dart e = visitedFaces[i] ;
// 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);
// }
unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee
map.setCurrentLevel(fLevel) ;
// // 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) ;
// }
//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
{
subdividedfacesT.push_back(std::pair<Dart,Dart>(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<Dart,Dart>(e,map.phi2(e)));
e = map.phi2(map.phi1(e));
}while (e != cf);
// e = map.phi1(e) ;
// } while(e != visitedFaces[i]) ;
// }
// volCenter /= typename PFP::REAL(count) ;
}
// /*
// * Subdivision
// */
// //Store the darts from quadrangulated faces
// std::vector<std::pair<Dart,Dart> > subdividedfacesQ;
// subdividedfacesQ.reserve(25);
map.setCurrentLevel(cur);
}
// std::vector<std::pair<Dart,Dart> > subdividedfacesT;
// subdividedfacesT.reserve(25);
map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision
std::vector<Dart> newEdges; //save darts from inner edges
newEdges.reserve(50);
// //First step : subdivide edges and faces
// //creates a i+1 edge level and i+1 face level
// for (std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
// {
// Dart d = *face;
//Second step : deconnect each corner, close each hole, subdivide each new face into 3
for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
{
Dart e = *edge;
// //if needed subdivide face
// if(!map.faceIsSubdivided(d))
// IHM::subdivideFace<PFP>(map, d, position);
std::vector<Dart> v ;
do
{
if(map.phi1(map.phi1(map.phi1(e))) != e)
v.push_back(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))
// //save a dart from the subdivided face
// unsigned int cur = map.getCurrentLevel() ;
v.push_back(map.phi1(e));
// unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee
// map.setCurrentLevel(fLevel) ;
e = map.phi2(map.phi_1(e));
}
while(e != *edge);
// //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
// {
// subdividedfacesT.push_back(std::pair<Dart,Dart>(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<Dart,Dart>(e,map.phi2(e)));
// e = map.phi2(map.phi1(e));
// }while (e != cf);
map.splitVolume(v) ;
//degree du sommet exterieur
unsigned int cornerDegree = map.Map2::vertexDegree(*edge);
// }
//tourner autour du sommet pour connaitre le brin d'un sommet de valence < cornerDegree
bool found = false;
Dart stop = e;
do
{
// map.setCurrentLevel(cur);
// }
if(map.Map2::vertexDegree(map.phi2(map.phi1(e))) < cornerDegree)
{
stop = map.phi2(map.phi1(e));
found = true;
}
// map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision
e = map.phi2(map.phi_1(e));
}
while(!found && e != *edge);
// std::vector<Dart> newEdges; //save darts from inner edges
// newEdges.reserve(50);
//si il existe un sommet de degre inferieur au degree du coin
if(found)
{
//chercher le brin de faible degree suivant
bool found2 = false;
Dart dd = map.phi1(stop);
// //Second step : deconnect each corner, close each hole, subdivide each new face into 3
// for (std::vector<Dart>::iterator edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
// {
// Dart e = *edge;
do
{
if(map.Map2::vertexDegree(dd) < cornerDegree)
found2 = true;
else
dd = map.phi1(dd);
}
while(!found2);
// std::vector<Dart> v ;
//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);
// do
// {
// if(map.phi1(map.phi1(map.phi1(e))) != e)
// v.push_back(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))
//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);
// v.push_back(map.phi1(e));
}
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);
}
// e = map.phi2(map.phi_1(e));
// }
// while(e != *edge);
}
// map.splitVolume(v) ;
}
//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);
}
// //degree du sommet exterieur
// unsigned int cornerDegree = map.Map2::vertexDegree(*edge);
}
// //tourner autour du sommet pour connaitre le brin d'un sommet de valence < cornerDegree
// bool found = false;
// Dart stop = e;
// do
// {
// 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);
// //si il existe un sommet de degre inferieur au degree du coin
// if(found)
// {
// //chercher le brin de faible degree suivant
// bool found2 = false;
// Dart dd = map.phi1(stop);
//std::cout << "1ere etape finished" << std::endl;
// do
// {
// if(map.Map2::vertexDegree(dd) < cornerDegree)
// found2 = true;
// else
// dd = map.phi1(dd);
// }
// while(!found2);
CellMarker<typename PFP::MAP, FACE> mtf(map);
// //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);
//Etape 2
for (std::vector<std::pair<Dart,Dart> >::iterator edges = subdividedfacesT.begin(); edges != subdividedfacesT.end(); ++edges)
{
// Dart f1 = (*edges).first;
Dart f2 = (*edges).second;
// //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;
//si ce n'est pas un tetrahedre
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))
{
// Dart ne = map.phi_1(stop) ;
// map.cutEdge(ne);
// position[map.phi1(ne)] = volCenter;
// stop = map.phi2(map.phi1(ne));
//map.deleteVolume(map.phi3(map.phi2(map.phi1(oldEdges.front()))));
// bool finished = false;
// Dart it = map.phi2(ne);
// 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<PFP>(map, map.phi2(f1));
// position[x] = volCenter;
// }
// else
// {
// //std::cout << "ajout d'un prisme" << std::endl;
// //Dart x = Algo::Modelisation::extrudeFace<PFP>(map,position,map.phi2(f1),5.0);
// Dart c = Algo::Modelisation::trianguleFace<PFP>(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))));
//
// }
// }
// //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);
}
}
//std::cout << "2e etape finished" << std::endl;
// {
// //Third step : 3-sew internal faces
// for (std::vector<std::pair<Dart,Dart> >::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));
//
//
// //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);
// }
//
// //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<Dart>::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<std::pair<Dart,Dart> >::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<Dart>::iterator it = newEdges.begin() ; it != newEdges.end() ; ++it)
// //sinon cas du tetraedre
// else
// {
// if(!mne.isMarked(*it))
// {
// unsigned int idedge = map.getNewEdgeId();
// map.setEdgeId(*it, idedge, EDGE);
// mne.markOrbit<EDGE>(*it);
// }
// //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);
// }
// }
//
// //cas tordu pour le prisme
// for (std::vector<std::pair<Dart,Dart> >::iterator it = subdividedfacesT.begin(); it != subdividedfacesT.end(); ++it)
// //std::cout << "1ere etape finished" << std::endl;
// CellMarker<typename PFP::MAP, FACE> mtf(map);
// //Etape 2
// for (std::vector<std::pair<Dart,Dart> >::iterator edges = subdividedfacesT.begin(); edges != subdividedfacesT.end(); ++edges)
// {
// 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))
//// Dart f1 = (*edges).first;
// Dart f2 = (*edges).second;
// //si ce n'est pas un tetrahedre
// 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))
// {
//
//
// 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.deleteVolume(map.phi3(map.phi2(map.phi1(oldEdges.front()))));
//// 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<PFP>(map, map.phi2(f1));
//// position[x] = volCenter;
//// }
//// else
//// {
//// //std::cout << "ajout d'un prisme" << std::endl;
//// //Dart x = Algo::Modelisation::extrudeFace<PFP>(map,position,map.phi2(f1),5.0);
//// Dart c = Algo::Modelisation::trianguleFace<PFP>(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);