Skip to content
Snippets Groups Projects
Commit 111fef6f authored by Emilie Marchandise's avatar Emilie Marchandise
Browse files

work on curvature for emi from Emi

parent 016dadb0
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/Curvature.cpp
+ 63
− 202
View file @ 111fef6f
...@@ -25,16 +25,23 @@ bool Curvature::_alreadyComputedCurvature = false; ...@@ -25,16 +25,23 @@ bool Curvature::_alreadyComputedCurvature = false;
//======================================================================================================== //========================================================================================================
//CONSTRUCTOR //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() Curvature::~Curvature()
{ {
_instance = 0; _instance = 0;
_destroyed = true; _destroyed = true;
} }
...@@ -43,31 +50,21 @@ void Curvature::onDeadReference() ...@@ -43,31 +50,21 @@ void Curvature::onDeadReference()
{ {
std::cout << "Dead reference of Curvature detected" << std::endl; std::cout << "Dead reference of Curvature detected" << std::endl;
} }
//======================================================================================================== //========================================================================================================
Curvature& Curvature::getInstance() Curvature& Curvature::getInstance()
{ {
if (!_instance) if (!_instance) {
{ if(_destroyed){
// Check for dead reference:
if(_destroyed)
{
onDeadReference(); onDeadReference();
} }
else else{
{
create(); create();
} }
} }
return *_instance; return *_instance;
} }
//======================================================================================================== //========================================================================================================
bool Curvature::valueAlreadyComputed() bool Curvature::valueAlreadyComputed()
{ {
return _alreadyComputedCurvature; return _alreadyComputedCurvature;
...@@ -89,43 +86,26 @@ Curvature& Curvature::getInstance() ...@@ -89,43 +86,26 @@ Curvature& Curvature::getInstance()
} }
//======================================================================================================== //========================================================================================================
void Curvature::retrieveCompounds() {
void Curvature::retrieveCompounds()
{
std::vector<GEntity*> entities; std::vector<GEntity*> entities;
_model->getEntities(entities); _model->getEntities(entities);
// std::cout << "The size of entities is " << entities.size() << std::endl; for(int ie = 0; ie < entities.size(); ++ie) {
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;
if( entities[ie]->geomType() == GEntity::CompoundSurface ) if( entities[ie]->geomType() == GEntity::CompoundSurface ) {
{
//std::cout << "This is compound surface" << std::endl;
GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]); GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]);
std::list<GFace*> tempcompounds = compound->getCompounds(); std::list<GFace*> tempcompounds = compound->getCompounds();
//std::cout << "The compound surface consists of " << tempcompounds.size() << " sub-surfaces " << std::endl;
std::list<GFace*>::iterator surfIterator; std::list<GFace*>::iterator surfIterator;
for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++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;
_ptFinalEntityList.push_back(*surfIterator); _ptFinalEntityList.push_back(*surfIterator);
} }
} }
} }
} }
//======================================================================================================== //========================================================================================================
//INITIALIZATION OF THE MAP AND RENUMBERING OF THE SELECTED ENTITIES: //INITIALIZATION OF THE MAP AND RENUMBERING OF THE SELECTED ENTITIES:
...@@ -170,17 +150,13 @@ void Curvature::initializeMap() ...@@ -170,17 +150,13 @@ void Curvature::initializeMap()
std::map<int,int>::iterator element_iterator; std::map<int,int>::iterator element_iterator;
for (vertex_iterator = _VertexToInt.begin() ; vertex_iterator !=_VertexToInt.end() ; ++ vertex_iterator, ++idx) for (vertex_iterator = _VertexToInt.begin() ; vertex_iterator !=_VertexToInt.end() ; ++ vertex_iterator, ++idx)
{
(*vertex_iterator).second = idx; (*vertex_iterator).second = idx;
}
idx = 0; idx = 0;
for (element_iterator = _ElementToInt.begin() ; element_iterator !=_ElementToInt.end() ; ++ element_iterator, ++idx) for (element_iterator = _ElementToInt.begin() ; element_iterator !=_ElementToInt.end() ; ++ element_iterator, ++idx)
{
(*element_iterator).second = idx; (*element_iterator).second = idx;
}
} }
//======================================================================================================== //========================================================================================================
...@@ -646,8 +622,7 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v, ...@@ -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); rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);
// if(unlikely(kuv !=0.0f)) // if(unlikely(kuv !=0.0f))
if(kuv !=0.0) if(kuv !=0.0) {
{
//Jacobi rotation to diagonalize //Jacobi rotation to diagonalize
double h= 0.5*(kv -ku)/ kuv; double h= 0.5*(kv -ku)/ kuv;
tt = (h<0.0) ? tt = (h<0.0) ?
...@@ -660,12 +635,10 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v, ...@@ -660,12 +635,10 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
k1 = ku - tt *kuv; k1 = ku - tt *kuv;
k2 = kv + 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; pdir1 = c*r_old_u - s*r_old_v;
} }
else else {
{
std::swap(k1,k2); std::swap(k1,k2);
pdir1 = s*r_old_u + c*r_old_v; 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, ...@@ -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(); retrieveCompounds();
initializeMap(); initializeMap();
computeRusinkiewiczNormals(); computeRusinkiewiczNormals();
...@@ -853,76 +825,25 @@ void Curvature::computeCurvature_Rusinkiewicz() ...@@ -853,76 +825,25 @@ void Curvature::computeCurvature_Rusinkiewicz()
for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
{ {
//**Mean Curvature:**
_VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5; if (isMax){
// _VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]); _VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
//**Gaussian Curvature:** }
// _VertexCurve[ivertex] = std::abs( _curv1[ivertex]*_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]); //_VertexCurve[ivertex] = std::abs(_VertexCurve[ivertex]);
}
} }
_alreadyComputedCurvature = true; _alreadyComputedCurvature = true;
} //End of the "computeCurvature_Rusinkiewicz" method } //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* A = triangle->getVertex(0);
MVertex* B = triangle->getVertex(1); MVertex* B = triangle->getVertex(1);
...@@ -933,51 +854,37 @@ void Curvature::computeCurvature_Rusinkiewicz() ...@@ -933,51 +854,37 @@ void Curvature::computeCurvature_Rusinkiewicz()
int V2 = 0; int V2 = 0;
std::map<int,int>::iterator vertexIterator; std::map<int,int>::iterator vertexIterator;
vertexIterator = _VertexToInt.find( A->getNum() ); vertexIterator = _VertexToInt.find( A->getNum() );
if ( vertexIterator != _VertexToInt.end() ) if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
{
V0 = (*vertexIterator).second;
}
else else
{
std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl; std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
}
vertexIterator = _VertexToInt.find( B->getNum() ); vertexIterator = _VertexToInt.find( B->getNum() );
if ( vertexIterator != _VertexToInt.end() ) if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
{
V1 = (*vertexIterator).second;
}
else else
{
std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl; std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
}
vertexIterator = _VertexToInt.find( C->getNum() ); vertexIterator = _VertexToInt.find( C->getNum() );
if ( vertexIterator != _VertexToInt.end() ) if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
{
V2 = (*vertexIterator).second;
}
else else
{
std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl; 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* A = edge->getVertex(0);
MVertex* B = edge->getVertex(1); MVertex* B = edge->getVertex(1);
...@@ -988,70 +895,24 @@ void Curvature::computeCurvature_Rusinkiewicz() ...@@ -988,70 +895,24 @@ void Curvature::computeCurvature_Rusinkiewicz()
std::map<int,int>::iterator vertexIterator; std::map<int,int>::iterator vertexIterator;
vertexIterator = _VertexToInt.find( A->getNum() ); vertexIterator = _VertexToInt.find( A->getNum() );
if ( vertexIterator != _VertexToInt.end() ) if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
{ else std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
V0 = (*vertexIterator).second;
}
else
{
std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
}
vertexIterator = _VertexToInt.find( B->getNum() ); vertexIterator = _VertexToInt.find( B->getNum() );
if ( vertexIterator != _VertexToInt.end() ) if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
{ else std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
V1 = (*vertexIterator).second;
if (isAbs){
c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
} }
else else{
{ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl; 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) void Curvature::writeToPosFile( const std::string & filename)
......
Geo/Curvature.h
+ 5
− 19
View file @ 111fef6f
...@@ -48,7 +48,6 @@ private: ...@@ -48,7 +48,6 @@ private:
//Model and list of selected entities with give physical tag: //Model and list of selected entities with give physical tag:
GModel* _model; GModel* _model;
GFaceList _ptFinalEntityList; GFaceList _ptFinalEntityList;
//Averaged vertex normals //Averaged vertex normals
...@@ -67,7 +66,6 @@ private: ...@@ -67,7 +66,6 @@ private:
std::vector<double> _pointareas; std::vector<double> _pointareas;
std::vector<SVector3> _cornerareas; std::vector<SVector3> _cornerareas;
//Curvature Tensor per mesh vertex //Curvature Tensor per mesh vertex
std::vector<STensor3> _CurveTensor; std::vector<STensor3> _CurveTensor;
...@@ -76,25 +74,18 @@ private: ...@@ -76,25 +74,18 @@ private:
//Area around a mesh vertex: //Area around a mesh vertex:
std::vector<double> _VertexArea; std::vector<double> _VertexArea;
std::vector<double> _VertexCurve; std::vector<double> _VertexCurve;
//----------------------------------------- //-----------------------------------------
// PRIVATE METHODS // PRIVATE METHODS
//Constructor //Constructor and destructor
Curvature(); Curvature();
//Destructor
~Curvature(); ~Curvature();
static void create(); static void create();
static void onDeadReference(); static void onDeadReference();
void initializeMap(); void initializeMap();
void computeVertexNormals(); void computeVertexNormals();
void curvatureTensor(); void curvatureTensor();
...@@ -170,8 +161,8 @@ public: ...@@ -170,8 +161,8 @@ public:
static Curvature& getInstance(); static Curvature& getInstance();
static bool valueAlreadyComputed(); static bool valueAlreadyComputed();
void setGModel(GModel* model); void setGModel(GModel* model);
void retrieveCompounds(); void retrieveCompounds();
//void retrievePhysicalSurfaces(const std::string & face_tag); //void retrievePhysicalSurfaces(const std::string & face_tag);
...@@ -185,15 +176,10 @@ public: ...@@ -185,15 +176,10 @@ public:
/// Estimating Curvatures and Their Derivatives on Triangle Meshes /// Estimating Curvatures and Their Derivatives on Triangle Meshes
/// Szymon Rusinkiewicz, Princeton University /// Szymon Rusinkiewicz, Princeton University
/// Code taken from Rusinkiewicz' 'trimesh2' library /// Code taken from Rusinkiewicz' 'trimesh2' library
void computeCurvature_Rusinkiewicz(); void computeCurvature_Rusinkiewicz(int isMax=0);
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 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); void writeToPosFile( const std::string & filename);
......
Geo/GFaceCompound.cpp
+ 8
− 21
View file @ 111fef6f
...@@ -1375,8 +1375,7 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const ...@@ -1375,8 +1375,7 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const
{ {
if(!oct) parametrize(); if(!oct) parametrize();
if(trivial()) if(trivial()) {
{
return (*(_compound.begin()))->curvatureMax(param); return (*(_compound.begin()))->curvatureMax(param);
} }
...@@ -1384,50 +1383,38 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const ...@@ -1384,50 +1383,38 @@ double GFaceCompound::curvatureMax(const SPoint2 &param) const
GFaceCompoundTriangle *lt; GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), &lt, U,V); getTriangle(param.x(), param.y(), &lt, U,V);
if(!lt) if(!lt) {
{
return 0.0; return 0.0;
} }
if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
{
//std::cout << "I'm not in DiscreteSurface" << std::endl;
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V; SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatureMax(pv); return lt->gf->curvatureMax(pv);
} }
else if (lt->gf->geomType() == GEntity::DiscreteSurface) else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
{
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 << "Need to compute discrete curvature" << std::endl;
std::cout << "Getting instance of curvature" << std::endl; std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() ); curvature.setGModel( model() );
curvature.computeCurvature_Rusinkiewicz(); int computeMax = 1;
curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos"); curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk"); curvature.writeToVtkFile("curvature.vtk");
std::cout << " ... computing curvature finished" << std::endl; std::cout << " ... computing curvature finished" << std::endl;
} }
// std::cout << "I'm in DiscreteSurface" << std::endl;
double c0; double c0;
double c1; double c1;
double c2; 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; double cv = (1-U-V)*c0 + U*c1 + V*c2;
return cv; return cv;
//Original code:
// double curv= 0.;
// curv = locCurvature(lt->tri,U,V);
// return curv;
} }
return 0.; return 0.;
......
Geo/OCCFace.cpp
+ 4
− 2
View file @ 111fef6f
...@@ -222,6 +222,7 @@ GEntity::GeomType OCCFace::geomType() const ...@@ -222,6 +222,7 @@ GEntity::GeomType OCCFace::geomType() const
double OCCFace::curvatureMax(const SPoint2 &param) const double OCCFace::curvatureMax(const SPoint2 &param) const
{ {
const double eps = 1.e-12; const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True); BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps); BRepLProp_SLProps prop(sf, 2, eps);
...@@ -230,16 +231,17 @@ double OCCFace::curvatureMax(const SPoint2 &param) const ...@@ -230,16 +231,17 @@ double OCCFace::curvatureMax(const SPoint2 &param) const
if (!prop.IsCurvatureDefined()){ if (!prop.IsCurvatureDefined()){
return eps; return eps;
} }
return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature())); return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
}
}
double OCCFace::curvatures(const SPoint2 &param, double OCCFace::curvatures(const SPoint2 &param,
SVector3 *dirMax, SVector3 *dirMax,
SVector3 *dirMin, SVector3 *dirMin,
double *curvMax, double *curvMax,
double *curvMin) const double *curvMin) const
{ {
const double eps = 1.e-12; const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True); BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps); BRepLProp_SLProps prop(sf, 2, eps);
......
Geo/discreteEdge.cpp
+ 6
− 9
View file @ 111fef6f
...@@ -492,32 +492,29 @@ SVector3 discreteEdge::firstDer(double par) const ...@@ -492,32 +492,29 @@ SVector3 discreteEdge::firstDer(double par) const
return der; return der;
} }
double discreteEdge::curvature(double par) const double discreteEdge::curvature(double par) const{
{
// std::cout << "Computing Curvature in Discrete Edges" << std::endl;
double tLoc; double tLoc;
int iEdge; int iEdge;
if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC; if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC;
double c0, c1; 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 << "Need to compute discrete curvature" << std::endl;
std::cout << "Getting instance of curvature" << std::endl; std::cout << "Getting instance of curvature" << std::endl;
curvature.setGModel( model() ); curvature.setGModel( model() );
curvature.computeCurvature_Rusinkiewicz(); int computeMax = 1;
curvature.computeCurvature_Rusinkiewicz(computeMax);
curvature.writeToPosFile("curvature.pos"); curvature.writeToPosFile("curvature.pos");
curvature.writeToVtkFile("curvature.vtk"); curvature.writeToVtkFile("curvature.vtk");
std::cout << " ... computing curvature finished" << std::endl; 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; double cv = (1-tLoc)*c0 + tLoc*c1;
......
benchmarks/curvature/TorusRemesh.geo
+ 2
− 2
View file @ 111fef6f
...@@ -17,8 +17,8 @@ Mesh.CharacteristicLengthExtendFromBoundary=0; ...@@ -17,8 +17,8 @@ Mesh.CharacteristicLengthExtendFromBoundary=0;
//Merge "TorusRemeshedBAMG.stl"; //Merge "TorusRemeshedBAMG.stl";
//Merge "SimpleTorus.msh"; //Merge "SimpleTorus.msh";
//Merge "TorusInput.stl"; Merge "TorusInput.stl";
Merge "occtorus.stl"; //Merge "occtorus.stl";
CreateTopology; CreateTopology;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment