Commit 1d99adf6 authored by untereiner's avatar untereiner

somes changes in volumetric sqrt(3) subdivision

parent a4379092
...@@ -33,6 +33,7 @@ find_package(Boost COMPONENTS regex thread system REQUIRED) ...@@ -33,6 +33,7 @@ find_package(Boost COMPONENTS regex thread system REQUIRED)
find_package(ZLIB REQUIRED) find_package(ZLIB REQUIRED)
find_package(LibXml2 REQUIRED) find_package(LibXml2 REQUIRED)
find_package(GLEW REQUIRED) find_package(GLEW REQUIRED)
find_package(SuiteSparse REQUIRED)
IF (DEFINED ASSERTON) IF (DEFINED ASSERTON)
add_definitions(-DCGOGN_ASSERT_BOOL=${ASSERTON}) add_definitions(-DCGOGN_ASSERT_BOOL=${ASSERTON})
...@@ -77,6 +78,7 @@ SET (CGoGN_EXT_INCLUDES ...@@ -77,6 +78,7 @@ SET (CGoGN_EXT_INCLUDES
# define libs for external libs # define libs for external libs
SET (CGoGN_EXT_LIBS SET (CGoGN_EXT_LIBS
nl
${OPENGL_LIBRARY} ${OPENGL_LIBRARY}
${GLEW_LIBRARIES} ${GLEW_LIBRARIES}
${ZLIB_LIBRARIES} ${ZLIB_LIBRARIES}
...@@ -84,6 +86,7 @@ SET (CGoGN_EXT_LIBS ...@@ -84,6 +86,7 @@ SET (CGoGN_EXT_LIBS
${Boost_SYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY}
${Boost_REGEX_LIBRARY} ${Boost_REGEX_LIBRARY}
${Boost_THREAD_LIBRARY} ${Boost_THREAD_LIBRARY}
${SUITESPARSE_LIBRARIES}
) )
#optionnal libs #optionnal libs
......
...@@ -32,6 +32,7 @@ find_package(Boost COMPONENTS regex thread system REQUIRED) ...@@ -32,6 +32,7 @@ find_package(Boost COMPONENTS regex thread system REQUIRED)
find_package(ZLIB REQUIRED) find_package(ZLIB REQUIRED)
find_package(LibXml2 REQUIRED) find_package(LibXml2 REQUIRED)
find_package(GLEW REQUIRED) find_package(GLEW REQUIRED)
find_package(SuiteSparse REQUIRED)
IF (DEFINED ASSERTON) IF (DEFINED ASSERTON)
add_definitions(-DCGOGN_ASSERT_BOOL=${ASSERTON}) add_definitions(-DCGOGN_ASSERT_BOOL=${ASSERTON})
...@@ -53,6 +54,7 @@ SET(COMMON_INCLUDES ...@@ -53,6 +54,7 @@ SET(COMMON_INCLUDES
# define libs for external libs # define libs for external libs
SET (COMMON_LIBS SET (COMMON_LIBS
nl
${OPENGL_LIBRARY} ${OPENGL_LIBRARY}
${GLEW_LIBRARIES} ${GLEW_LIBRARIES}
${ZLIB_LIBRARIES} ${ZLIB_LIBRARIES}
...@@ -60,6 +62,7 @@ SET (COMMON_LIBS ...@@ -60,6 +62,7 @@ SET (COMMON_LIBS
${Boost_SYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY}
${Boost_REGEX_LIBRARY} ${Boost_REGEX_LIBRARY}
${Boost_THREAD_LIBRARY} ${Boost_THREAD_LIBRARY}
${SUITESPARSE_LIBRARIES}
) )
#optionnal libs #optionnal libs
......
...@@ -672,7 +672,7 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit ...@@ -672,7 +672,7 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
Traversor3WF<typename PFP::MAP> tWF(map, dit); Traversor3WF<typename PFP::MAP> tWF(map, dit);
for(Dart ditWF = tWF.begin() ; ditWF != tWF.end() ; ditWF = tWF.next()) for(Dart ditWF = tWF.begin() ; ditWF != tWF.end() ; ditWF = tWF.next())
{ {
if(!map.isBoundaryFace(ditWF)) if(!map.isBoundaryFace(ditWF) && !m.isMarked(ditWF))
m.markOrbit<FACE>(ditWF); m.markOrbit<FACE>(ditWF);
} }
...@@ -711,7 +711,7 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit ...@@ -711,7 +711,7 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
Traversor3WE<typename PFP::MAP> tWE(map, dit); Traversor3WE<typename PFP::MAP> tWE(map, dit);
for(Dart ditWE = tWE.begin() ; ditWE != tWE.end() ; ditWE = tWE.next()) for(Dart ditWE = tWE.begin() ; ditWE != tWE.end() ; ditWE = tWE.next())
{ {
if(map.isBoundaryEdge(ditWE)) if(map.isBoundaryEdge(ditWE) && !m.isMarked(ditWE))
m.markOrbit<EDGE>(ditWE); m.markOrbit<EDGE>(ditWE);
} }
...@@ -728,45 +728,73 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit ...@@ -728,45 +728,73 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
} }
} }
//
// edge-removal on all old boundary edges
//
TraversorE<typename PFP::MAP> tE(map,selected);
for(Dart dit = tE.begin() ; dit != tE.end() ; dit = tE.next())
{
if(m.isMarked(dit))
{
m.unmarkOrbit<EDGE>(dit);
Dart d = map.phi2(map.phi3(map.findBoundaryFaceOfEdge(dit)));
Volume::Modelisation::Tetrahedralization::swapGen3To2<PFP>(map, d);
}
}
TraversorV<typename PFP::MAP> tVg(map,selected); TraversorV<typename PFP::MAP> tVg(map,selected);
for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next()) for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
{ {
if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit)) if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
{ {
typename PFP::VEC3 P = position[dit] ; Dart db = map.findBoundaryFaceOfVertex(dit);
typename PFP::VEC3 P = position[db] ;
typename PFP::VEC3 newP(0) ; typename PFP::VEC3 newP(0) ;
unsigned int val = 0 ; unsigned int val = 0 ;
Dart vit = dit ; Dart vit = db ;
do do
{ {
newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ; newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
++val ; ++val ;
vit = map.phi2(map.phi_1(vit)) ; vit = map.phi2(map.phi_1(vit)) ;
} while(vit != dit) ; } while(vit != db) ;
typename PFP::REAL K = sqrt3_K(val) ; typename PFP::REAL K = sqrt3_K(val) ;
newP *= typename PFP::REAL(3) ; newP *= typename PFP::REAL(3) ;
newP -= typename PFP::REAL(val) * P ; newP -= typename PFP::REAL(val) * P ;
newP *= K / typename PFP::REAL(2 * val) ; newP *= K / typename PFP::REAL(2 * val) ;
newP += (typename PFP::REAL(1) - K) * P ; newP += (typename PFP::REAL(1) - K) * P ;
position[dit] = newP ; position[db] = newP ;
}
}
//
// edge-removal on all old boundary edges
//
TraversorE<typename PFP::MAP> tE(map,selected);
for(Dart dit = tE.begin() ; dit != tE.end() ; dit = tE.next())
{
if(m.isMarked(dit))
{
m.unmarkOrbit<EDGE>(dit);
Dart d = map.phi2(map.phi3(map.findBoundaryFaceOfEdge(dit)));
Volume::Modelisation::Tetrahedralization::swapGen3To2<PFP>(map, d);
} }
} }
// TraversorV<typename PFP::MAP> tVg(map,selected);
// for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
// {
// if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
// {
// Dart db = map.findBoundaryFaceOfVertex(dit);
//
// typename PFP::VEC3 P = position[db] ;
// typename PFP::VEC3 newP(0) ;
// unsigned int val = 0 ;
// Dart vit = db ;
// do
// {
// newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
// ++val ;
// vit = map.phi2(map.phi_1(vit)) ;
// } while(vit != db) ;
// typename PFP::REAL K = sqrt3_K(val) ;
// newP *= typename PFP::REAL(3) ;
// newP -= typename PFP::REAL(val) * P ;
// newP *= K / typename PFP::REAL(2 * val) ;
// newP += (typename PFP::REAL(1) - K) * P ;
// position[db] = newP ;
// }
// }
//AutoVertexAttribute laplacian qui est une copie de position //AutoVertexAttribute laplacian qui est une copie de position
// TraversorV<typename PFP::MAP> tVg2(map,selected); // TraversorV<typename PFP::MAP> tVg2(map,selected);
...@@ -781,108 +809,48 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit ...@@ -781,108 +809,48 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& posit
//echange lapaclian et position //echange lapaclian et position
VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord"); // VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord");
Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ; // Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ;
VertexAutoAttribute<unsigned int> vIndex(map, "vIndex");
unsigned int nb_vertices = map.template computeIndexCells<VERTEX>(vIndex);
CellMarker<VERTEX> lockingMarker(map);
TraversorV<typename PFP::MAP> tv(map);
for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
{
if(!lockingMarker.isMarked(dit) && map.isBoundaryVertex(dit))
lockingMarker.mark(dit);
}
NLContext nlContext = nlNewContext();
nlSolverParameteri(NL_NB_VARIABLES, nb_vertices);
nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
nlMakeCurrent(nlContext);
if(nlGetCurrentState() == NL_STATE_INITIAL)
nlBegin(NL_SYSTEM) ;
for(int coord = 0; coord < 3; ++coord)
{
LinearSolving::setupVariables<PFP>(map, vIndex, lockingMarker, position, coord) ;
nlBegin(NL_MATRIX) ;
LinearSolving::addRowsRHS_Laplacian_Topo<PFP>(map, vIndex, diffCoord, coord) ;
// LinearSolving::addRowsRHS_Laplacian_Cotan<PFP>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
nlEnd(NL_MATRIX) ;
nlEnd(NL_SYSTEM) ;
nlSolve() ;
LinearSolving::getResult<PFP>(map, vIndex, position, coord) ;
nlReset(NL_TRUE) ;
}
nlDeleteContext(nlContext);
// //
// float weight = 1.0; // VertexAutoAttribute<unsigned int> vIndex(map, "vIndex");
//
// LinearSolving::LinearSolver<PFP::REAL> ls(nbV); // unsigned int nb_vertices = map.template computeIndexCells<VERTEX>(vIndex);
// ls.set_least_squares(true); //
//
// for(unsigned int coord = 0 ; coord < 3 ; ++coord) // CellMarker<VERTEX> lockingMarker(map);
//
// TraversorV<typename PFP::MAP> tv(map);
// for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
// { // {
// std::cout << "coord " << coord << std::flush; // if(map.isBoundaryVertex(dit))
// lockingMarker.mark(dit);
// TraversorV<PFP::MAP> tv(map); // }
// for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next()) //
// { //
// ls.variable(indexV[dit]).set_value(position[dit][coord]); // NLContext nlContext = nlNewContext();
// nlSolverParameteri(NL_NB_VARIABLES, nb_vertices);
// if(map.isBoundaryVertex(dit)) // nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
// ls.variable(indexV[dit]).lock(); // nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);
// } //
// std::cout << "... variables set... " << std::flush; //
// nlMakeCurrent(nlContext);
// ls.begin_system(); // if(nlGetCurrentState() == NL_STATE_INITIAL)
// nlBegin(NL_SYSTEM) ;
// TraversorV<PFP::MAP> tv2(map); //
// for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next()) // for(int coord = 0; coord < 3; ++coord)
// { // {
// if(!map.isBoundaryVertex(dit)) // LinearSolving::setupVariables<PFP>(map, vIndex, lockingMarker, position, coord) ;
// { // nlBegin(NL_MATRIX) ;
// ls.begin_row(); // LinearSolving::addRowsRHS_Laplacian_Topo<PFP>(map, vIndex, diffCoord, coord) ;
//// LinearSolving::addRowsRHS_Laplacian_Cotan<PFP>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
// float sum = 0; // nlEnd(NL_MATRIX) ;
// Traversor3VVaE<PFP::MAP> tvvae(map, dit); // nlEnd(NL_SYSTEM) ;
// for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next()) // nlSolve() ;
// { // LinearSolving::getResult<PFP>(map, vIndex, position, coord) ;
// ls.add_coefficient(indexV[ditvvae],weight); // nlReset(NL_TRUE) ;
// sum += weight;
// }
// ls.add_coefficient(indexV[dit],-sum);
// ls.normalize_row();
// ls.end_row();
// }
// }
// ls.end_system();
// std::cout << "... system built... " << std::flush;
// ls.solve();
// std::cout << "... system solved... " << std::flush;
// TraversorV<PFP::MAP> tv3(map);
// for(Dart dit = tv3.begin() ; dit != tv3.end() ; dit = tv3.next())
// {
// position[dit][coord] = ls.variable(indexV[dit]).value();
// }
// ls.reset(false);
// std::cout << "... done" << std::endl;
// } // }
//
// nlDeleteContext(nlContext);
} }
template <typename PFP> template <typename PFP>
......
...@@ -114,6 +114,11 @@ public: ...@@ -114,6 +114,11 @@ public:
{} {}
void operator() () void operator() ()
{
}
void operator() (bool filtering)
{ {
TraversorF<typename PFP::MAP> trav(m_map) ; TraversorF<typename PFP::MAP> trav(m_map) ;
for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
...@@ -145,10 +150,6 @@ public: ...@@ -145,10 +150,6 @@ public:
} }
} }
void operator() (bool filtering)
{
}
} ; } ;
template <typename PFP> template <typename PFP>
...@@ -163,6 +164,11 @@ public: ...@@ -163,6 +164,11 @@ public:
{} {}
void operator() () void operator() ()
{
}
void operator() (bool filtering)
{ {
TraversorV<typename PFP::MAP> trav(m_map) ; TraversorV<typename PFP::MAP> trav(m_map) ;
for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
...@@ -234,10 +240,6 @@ public: ...@@ -234,10 +240,6 @@ public:
} }
} }
void operator() (bool filtering)
{
}
} ; } ;
template <typename PFP> template <typename PFP>
...@@ -252,6 +254,11 @@ public: ...@@ -252,6 +254,11 @@ public:
{} {}
void operator() () void operator() ()
{
}
void operator() (bool filtering)
{ {
TraversorV<typename PFP::MAP> trav(m_map) ; TraversorV<typename PFP::MAP> trav(m_map) ;
for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
...@@ -289,10 +296,6 @@ public: ...@@ -289,10 +296,6 @@ public:
} }
} }
void operator() (bool filtering)
{
}
} ; } ;
template <typename PFP> template <typename PFP>
...@@ -307,6 +310,11 @@ public: ...@@ -307,6 +310,11 @@ public:
{} {}
void operator() () void operator() ()
{
}
void operator() (bool filtering)
{ {
TraversorF<typename PFP::MAP> trav(m_map) ; TraversorF<typename PFP::MAP> trav(m_map) ;
for (Dart d = trav.begin(); d != trav.end(); d = trav.next()) for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
...@@ -331,10 +339,6 @@ public: ...@@ -331,10 +339,6 @@ public:
} }
} }
void operator() (bool filtering)
{
}
} ; } ;
/********************************************************************************* /*********************************************************************************
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment