From 111fef6fe0e3412e8a848bb9635864f273209710 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 14 Jul 2011 12:51:08 +0000
Subject: [PATCH] work on curvature for emi from Emi
---
Geo/Curvature.cpp | 265 +++++++--------------------
Geo/Curvature.h | 24 +--
Geo/GFaceCompound.cpp | 29 +--
Geo/OCCFace.cpp | 6 +-
Geo/discreteEdge.cpp | 15 +-
benchmarks/curvature/TorusRemesh.geo | 4 +-
6 files changed, 88 insertions(+), 255 deletions(-)
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 0622aa7742..1780772f69 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -25,16 +25,23 @@ bool Curvature::_alreadyComputedCurvature = false;
//========================================================================================================
//CONSTRUCTOR
-Curvature::Curvature()
-{
-
+Curvature::Curvature(){
}
+
+// Curvature::Curvature(std::list<GFace*> &myFaces){
+
+// std::list<GFace*>::const_iterator it = myFaces.begin();
+// for( ; it != myFaces.end() ; ++it){
+// _ptFinalEntityList.push_back(*it);
+// }
+
+// }
+
//========================================================================================================
-//DESTRUCTOR
Curvature::~Curvature()
{
- _instance = 0;
+ _instance = 0;
_destroyed = true;
}
@@ -43,31 +50,21 @@ void Curvature::onDeadReference()
{
std::cout << "Dead reference of Curvature detected" << std::endl;
}
-
-
//========================================================================================================
Curvature& Curvature::getInstance()
{
- if (!_instance)
- {
- // Check for dead reference:
- if(_destroyed)
- {
+ if (!_instance) {
+ if(_destroyed){
onDeadReference();
}
- else
- {
+ else{
create();
}
}
return *_instance;
-
}
-
//========================================================================================================
-
-
bool Curvature::valueAlreadyComputed()
{
return _alreadyComputedCurvature;
@@ -89,43 +86,26 @@ Curvature& Curvature::getInstance()
}
//========================================================================================================
+ void Curvature::retrieveCompounds() {
- void Curvature::retrieveCompounds()
- {
std::vector<GEntity*> entities;
_model->getEntities(entities);
-// std::cout << "The size of entities is " << entities.size() << std::endl;
- for(int ie = 0; ie < entities.size(); ++ie)
- {
-// std::cout << "Entity " << ie << ":" << std::endl;
-// std::cout << "\ttag = " << entities[ie]->tag() << std::endl;
-// std::cout << "\t" << entities[ie]->getInfoString() << std::endl;
-
+ for(int ie = 0; ie < entities.size(); ++ie) {
- if( entities[ie]->geomType() == GEntity::CompoundSurface )
- {
- //std::cout << "This is compound surface" << std::endl;
+ if( entities[ie]->geomType() == GEntity::CompoundSurface ) {
GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]);
-
std::list<GFace*> tempcompounds = compound->getCompounds();
- //std::cout << "The compound surface consists of " << tempcompounds.size() << " sub-surfaces " << std::endl;
-
std::list<GFace*>::iterator surfIterator;
- for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator)
- {
-// std::cout << "TAG = " << (*surfIterator)->tag() << std::endl;
-// std::cout << "THIS SURFACE HAS " << (*surfIterator)->getNumMeshElements() << std::endl;
+ for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator) {
_ptFinalEntityList.push_back(*surfIterator);
}
}
}
-
}
-
//========================================================================================================
//INITIALIZATION OF THE MAP AND RENUMBERING OF THE SELECTED ENTITIES:
@@ -170,17 +150,13 @@ void Curvature::initializeMap()
std::map<int,int>::iterator element_iterator;
for (vertex_iterator = _VertexToInt.begin() ; vertex_iterator !=_VertexToInt.end() ; ++ vertex_iterator, ++idx)
- {
(*vertex_iterator).second = idx;
- }
-
+
idx = 0;
for (element_iterator = _ElementToInt.begin() ; element_iterator !=_ElementToInt.end() ; ++ element_iterator, ++idx)
- {
(*element_iterator).second = idx;
- }
-
+
}
//========================================================================================================
@@ -646,8 +622,7 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);
// if(unlikely(kuv !=0.0f))
- if(kuv !=0.0)
- {
+ if(kuv !=0.0) {
//Jacobi rotation to diagonalize
double h= 0.5*(kv -ku)/ kuv;
tt = (h<0.0) ?
@@ -660,12 +635,10 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
k1 = ku - tt *kuv;
k2 = kv + tt *kuv;
- if(std::abs(k1) >= std::abs(k2))
- {
+ if(std::abs(k1) >= std::abs(k2)) {
pdir1 = c*r_old_u - s*r_old_v;
}
- else
- {
+ else {
std::swap(k1,k2);
pdir1 = s*r_old_u + c*r_old_v;
}
@@ -674,8 +647,7 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
//========================================================================================================
-void Curvature::computeCurvature_Rusinkiewicz()
-{
+void Curvature::computeCurvature_Rusinkiewicz(int isMax){
retrieveCompounds();
initializeMap();
computeRusinkiewiczNormals();
@@ -853,76 +825,25 @@ void Curvature::computeCurvature_Rusinkiewicz()
for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
{
- //**Mean Curvature:**
- _VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5;
-// _VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]);
- //**Gaussian Curvature:**
-// _VertexCurve[ivertex] = std::abs( _curv1[ivertex]*_curv2[ivertex] );
+
+ if (isMax){
+ _VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
+ }
+ else{
+ _VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5; //Mean curvature
+ //_VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]);
+ //_VertexCurve[ivertex] = std::abs( _curv1[ivertex]*_curv2[ivertex] ); //Gaussian
//_VertexCurve[ivertex] = std::abs(_VertexCurve[ivertex]);
+ }
}
_alreadyComputedCurvature = true;
} //End of the "computeCurvature_Rusinkiewicz" method
-//========================================================================================================
-
- void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2)
- {
- MVertex* A = triangle->getVertex(0);
- MVertex* B = triangle->getVertex(1);
- MVertex* C = triangle->getVertex(2);
-
- int V0 = 0;
- int V1 = 0;
- int V2 = 0;
-
- std::map<int,int>::iterator vertexIterator;
-
- vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V0 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
- }
-
- vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V1 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
- }
-
- vertexIterator = _VertexToInt.find( C->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V2 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
- }
-
-
-// const int V0 = _VertexToInt[A->getNum()];
-// const int V1 = _VertexToInt[B->getNum()];
-// const int V2 = _VertexToInt[C->getNum()];
-
- c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
- c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
- c2 = _VertexCurve[V2]; //Mean curvature in vertex 2
-
- }
-
//========================================================================================================
- void Curvature::triangleNodalAbsValues(MTriangle* triangle, double& c0, double& c1, double& c2)
+void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs)
{
MVertex* A = triangle->getVertex(0);
MVertex* B = triangle->getVertex(1);
@@ -933,51 +854,37 @@ void Curvature::computeCurvature_Rusinkiewicz()
int V2 = 0;
std::map<int,int>::iterator vertexIterator;
-
vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V0 = (*vertexIterator).second;
- }
+ if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
else
- {
std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
- }
-
+
vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V1 = (*vertexIterator).second;
- }
+ if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
else
- {
std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
- }
-
+
vertexIterator = _VertexToInt.find( C->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V2 = (*vertexIterator).second;
- }
+ if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
else
- {
std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
+
+ if (isAbs){
+ c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+ c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
+ c2 = std::abs(_VertexCurve[V2]); //Mean curvature in vertex 2
+ }
+ else{
+ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+ c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
+ c2 = _VertexCurve[V2]; //Mean curvature in vertex 2
}
-
-
- // const int V0 = _VertexToInt[A->getNum()];
- // const int V1 = _VertexToInt[B->getNum()];
- // const int V2 = _VertexToInt[C->getNum()];
-
- c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
- c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
- c2 = std::abs(_VertexCurve[V2]); //Mean curvature in vertex 2
}
//========================================================================================================
- void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1)
+void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
{
MVertex* A = edge->getVertex(0);
MVertex* B = edge->getVertex(1);
@@ -988,70 +895,24 @@ void Curvature::computeCurvature_Rusinkiewicz()
std::map<int,int>::iterator vertexIterator;
vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V0 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
- }
-
+ if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
+ else std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
+
vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V1 = (*vertexIterator).second;
+ if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
+ else std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
+
+ if (isAbs){
+ c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+ c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
}
- else
- {
- std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
+ else{
+ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+ c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
}
- c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
- c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
-
}
-//========================================================================================================
-
- void Curvature::edgeNodalAbsValues(MLine* edge, double& c0, double& c1)
- {
- MVertex* A = edge->getVertex(0);
- MVertex* B = edge->getVertex(1);
-
- int V0 = 0;
- int V1 = 0;
-
- std::map<int,int>::iterator vertexIterator;
-
- vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V0 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
- }
-
- vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() )
- {
- V1 = (*vertexIterator).second;
- }
- else
- {
- std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
- }
-
- c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
- c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
-
- }
-
-
-
-
//========================================================================================================
void Curvature::writeToPosFile( const std::string & filename)
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 6c7d30b4f5..59876b6b70 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -48,7 +48,6 @@ private:
//Model and list of selected entities with give physical tag:
GModel* _model;
-
GFaceList _ptFinalEntityList;
//Averaged vertex normals
@@ -67,7 +66,6 @@ private:
std::vector<double> _pointareas;
std::vector<SVector3> _cornerareas;
-
//Curvature Tensor per mesh vertex
std::vector<STensor3> _CurveTensor;
@@ -76,25 +74,18 @@ private:
//Area around a mesh vertex:
std::vector<double> _VertexArea;
-
std::vector<double> _VertexCurve;
-
-
//-----------------------------------------
// PRIVATE METHODS
- //Constructor
+ //Constructor and destructor
Curvature();
-
- //Destructor
~Curvature();
static void create();
-
static void onDeadReference();
-
void initializeMap();
void computeVertexNormals();
void curvatureTensor();
@@ -170,8 +161,8 @@ public:
static Curvature& getInstance();
static bool valueAlreadyComputed();
+
void setGModel(GModel* model);
-
void retrieveCompounds();
//void retrievePhysicalSurfaces(const std::string & face_tag);
@@ -185,15 +176,10 @@ public:
/// Estimating Curvatures and Their Derivatives on Triangle Meshes
/// Szymon Rusinkiewicz, Princeton University
/// Code taken from Rusinkiewicz' 'trimesh2' library
- void computeCurvature_Rusinkiewicz();
-
- void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2);
-
- void triangleNodalAbsValues(MTriangle* triangle, double& c0, double& c1, double& c2);
-
- void edgeNodalValues(MLine* edge, double& c0, double& c1);
+ void computeCurvature_Rusinkiewicz(int isMax=0);
- void edgeNodalAbsValues(MLine* edge, double& c0, double& c1);
+ void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs=0);
+ void edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs=0);
void writeToPosFile( const std::string & filename);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 346950716b..f63a2251a1 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1375,8 +1375,7 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
if(!oct) parametrize();
- if(trivial())
- {
+ if(trivial()) {
return (*(_compound.begin()))->curvatureMax(param);
}
@@ -1384,50 +1383,38 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
- if(!lt)
- {
+ if(!lt) {
return 0.0;
}
- if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface)
- {
- //std::cout << "I'm not in DiscreteSurface" << std::endl;
+ if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatureMax(pv);
}
- else if (lt->gf->geomType() == GEntity::DiscreteSurface)
- {
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
Curvature& curvature = Curvature::getInstance();
- if( !Curvature::valueAlreadyComputed() )
- {
+ if( !Curvature::valueAlreadyComputed() ) {
std::cout << "Need to compute discrete curvature" << std::endl;
std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() );
- curvature.computeCurvature_Rusinkiewicz();
+ int computeMax = 1;
+ curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk");
-
std::cout << " ... computing curvature finished" << std::endl;
-
}
-
-// std::cout << "I'm in DiscreteSurface" << std::endl;
double c0;
double c1;
double c2;
- curvature.triangleNodalAbsValues(lt->tri,c0, c1, c2);
+ curvature.triangleNodalValues(lt->tri,c0, c1, c2, 1);
double cv = (1-U-V)*c0 + U*c1 + V*c2;
return cv;
-//Original code:
-// double curv= 0.;
-// curv = locCurvature(lt->tri,U,V);
-// return curv;
}
return 0.;
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 238ffac0cf..c0b3b3369b 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -222,6 +222,7 @@ GEntity::GeomType OCCFace::geomType() const
double OCCFace::curvatureMax(const SPoint2 ¶m) const
{
+
const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps);
@@ -230,16 +231,17 @@ double OCCFace::curvatureMax(const SPoint2 ¶m) const
if (!prop.IsCurvatureDefined()){
return eps;
}
+
return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
-}
-
+}
double OCCFace::curvatures(const SPoint2 ¶m,
SVector3 *dirMax,
SVector3 *dirMin,
double *curvMax,
double *curvMin) const
{
+
const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps);
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 4ed94e0d98..17f6682b44 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -492,32 +492,29 @@ SVector3 discreteEdge::firstDer(double par) const
return der;
}
-double discreteEdge::curvature(double par) const
-{
-// std::cout << "Computing Curvature in Discrete Edges" << std::endl;
+double discreteEdge::curvature(double par) const{
double tLoc;
int iEdge;
if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC;
double c0, c1;
- Curvature& curvature = Curvature::getInstance();
+ Curvature& curvature = Curvature::getInstance();
- if( !Curvature::valueAlreadyComputed() )
- {
+ if( !Curvature::valueAlreadyComputed() ) {
std::cout << "Need to compute discrete curvature" << std::endl;
std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() );
- curvature.computeCurvature_Rusinkiewicz();
+ int computeMax = 1;
+ curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk");
std::cout << " ... computing curvature finished" << std::endl;
-
}
- curvature.edgeNodalAbsValues(lines[iEdge],c0, c1);
+ curvature.edgeNodalValues(lines[iEdge],c0, c1, 1);
double cv = (1-tLoc)*c0 + tLoc*c1;
diff --git a/benchmarks/curvature/TorusRemesh.geo b/benchmarks/curvature/TorusRemesh.geo
index 68624b03e8..eb97cef048 100644
--- a/benchmarks/curvature/TorusRemesh.geo
+++ b/benchmarks/curvature/TorusRemesh.geo
@@ -17,8 +17,8 @@ Mesh.CharacteristicLengthExtendFromBoundary=0;
//Merge "TorusRemeshedBAMG.stl";
//Merge "SimpleTorus.msh";
-//Merge "TorusInput.stl";
-Merge "occtorus.stl";
+Merge "TorusInput.stl";
+//Merge "occtorus.stl";
CreateTopology;
--
GitLab