diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 358daf5f2dacd1771c6f7b6016e49fcfe18af371..9a41dee7f44986e64c4403fef6f2e99742cd458d 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -3,9 +3,6 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#include <iostream>
-#include <fstream>
-#include <cmath>
#include "Curvature.h"
#include "MElement.h"
#include "MTriangle.h"
@@ -15,26 +12,53 @@
#include "GRbf.h"
#include "OS.h"
#include "SBoundingBox3d.h"
+#include "discreteEdge.h"
+
+#include<iostream>
+#include<fstream>
+#include<cmath>
#define NEXT(i) ((i)<2 ? (i)+1 : (i)-2)
#define PREV(i) ((i)>0 ? (i)-1 : (i)+2)
+//========================================================================================================
+//Initialization of the static variables:
Curvature* Curvature::_instance = 0;
bool Curvature::_destroyed = false;
bool Curvature::_alreadyComputedCurvature = false;
+//========================================================================================================
+
+//CONSTRUCTOR
+Curvature::Curvature() : _isMapInitialized(false)
+{
+}
+
+// Curvature::Curvature(std::list<GFace*> &myFaces){
+
+// std::list<GFace*>::const_iterator it = myFaces.begin();
+// for( ; it != myFaces.end() ; ++it){
+// _ptFinalEntityList.push_back(*it);
+// }
+
+// }
+
+//========================================================================================================
+
Curvature::~Curvature()
{
- _instance = 0;
+ _instance = 0;
_destroyed = true;
-}
+}
+//========================================================================================================
void Curvature::onDeadReference()
{
std::cout << "Dead reference of Curvature detected" << std::endl;
}
+//========================================================================================================
-Curvature &Curvature::getInstance()
+Curvature& Curvature::getInstance()
{
if (!_instance) {
if(_destroyed){
@@ -46,83 +70,313 @@ Curvature &Curvature::getInstance()
}
return *_instance;
}
+//========================================================================================================
+ bool Curvature::valueAlreadyComputed()
+ {
+ return _alreadyComputedCurvature;
+ }
-bool Curvature::valueAlreadyComputed()
-{
- return _alreadyComputedCurvature;
-}
+//========================================================================================================
-void Curvature::create()
-{
- static Curvature instance;
- _instance = & instance;
-}
+ void Curvature::create()
+ {
+ static Curvature instance;
+ _instance = & instance;
+ }
-void Curvature::retrieveCompounds()
-{
+ //========================================================================================================
+ void Curvature::retrieveCompounds() {
#if defined(HAVE_SOLVER)
- std::vector<GEntity*> entities;
- _model->getEntities(entities);
-
- for(unsigned ie = 0; ie < entities.size(); ++ie) {
-
- if(entities[ie]->geomType() == GEntity::CompoundSurface) {
- GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]);
- std::list<GFace*> tempcompounds = compound->getCompounds();
- std::list<GFace*>::iterator surfIterator;
- for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end();
- ++surfIterator) {
- if ((*surfIterator)->geomType() == GEntity::DiscreteSurface) {
- _ptFinalEntityList.push_back(*surfIterator);
- }
- }
- }
- else if (entities[ie]->geomType() == GEntity::DiscreteSurface) {
- _ptFinalEntityList.push_back(dynamic_cast<GFace*>(entities[ie]));
- }
- }
+
+ /// -------------------------------------------------------------------------------------------
+ // Loop over all faces. Check if the face is a compound. If it is, store all of its discrete
+ // faces in _EntityArray
+
+ std::list<GFace*> global_face_list;
+
+ for(GModel::fiter face_iter = _model->firstFace(); face_iter != _model->lastFace(); ++face_iter)
+ {
+ if ( (*face_iter)->geomType() == GEntity::CompoundSurface )
+ {
+
+ GFaceCompound* compound = dynamic_cast<GFaceCompound*>(*face_iter);
+ std::list<GFace*> tempcompounds = compound->getCompounds();
+ std::list<GFace*>::iterator surfIterator;
+
+ for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator)
+ {
+ if ((*surfIterator)->geomType() == GEntity::DiscreteSurface)
+ {
+ global_face_list.push_back(*surfIterator);
+ }
+ }
+ }
+
+ } // Loop over faces of the model
+
+ global_face_list.sort();
+ global_face_list.unique();
+ _EntityArray.resize(global_face_list.size());
+ std::copy(global_face_list.begin(),global_face_list.end(),_EntityArray.begin());
+
+// std::cout << "Retrieve compounds: some preliminary info:" << std::endl;
+// std::cout << "The size of _ptFinalEntityList = " << _EntityArray.size() << std::endl;
+// std::cout << "Number of elements in all faces: " << std::endl;
+// for(int i = 0; i < _EntityArray.size(); ++i)
+// {
+// std::cout << "\t" << _EntityArray[i]->getNumMeshElements() << " elements" << std::endl;
+// }
+
#endif
-}
+ }
+
+ //========================================================================================================
+
+ // Method to detect boundary edges of the mesh. We want to find which edges are on the boundary in order
+ // to correct the curvature field close to boundaries
+
+ void Curvature::correctOnEdges()
+ {
+#if defined(HAVE_SOLVER)
+
+ if (! _alreadyComputedCurvature )
+ {
+ Msg::Error("Cannot correct the curvature values at the edges because the curvatures weren't computed yet");
+ return;
+ }
+ // Remark: this can only be used after the call to initializeMap() !
+
+ buildEdgeList();
+
+ std::list<MeshEdgeInfo>::iterator VertToEdgeListIter;
+
+ _isOnBoundary.resize(_VertexToInt.size());
+ for (int i = 0; i < _VertexToInt.size(); ++i)
+ {
+ _isOnBoundary[i] = 0;
+ }
+
+ // To detect the nodes on the egdes of a geometry, we create a list of all edges on the mesh. The
+ // edges which are shared by only one mesh element are boundary edges. Their nodes are tagged by 1
+ for (int i = 0; i < _VertexToInt.size(); ++i)
+ {
+ for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
+ ++VertToEdgeListIter)
+ {
+ if ((*VertToEdgeListIter).NbElemNeighbour == 1)
+ {
+ _isOnBoundary[(*VertToEdgeListIter).StartV] = 1;
+ _isOnBoundary[(*VertToEdgeListIter).EndV] = 1;
+ }
+ }
+
+ }
+
+ // We want to find the nodes that are the immediate neighbours of the boundary nodes. Those nodes are
+ // considered nodes with 'level 2'. The neighbours of neighbours have 'level 3' etc.
+ // We want to construct levels 1,2,3. Nodes with level 0 are internal nodes of the mesh
+ //Loop over the vector of edgelist
+ for (int level = 1; level < 3; ++level)
+ {
+ for (int i = 0; i < _VertexToEdgeList.size(); ++i)
+ {
+ for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
+ ++VertToEdgeListIter)
+ {
+ if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level && _isOnBoundary[(*VertToEdgeListIter).EndV] == 0)
+ {
+ _isOnBoundary[(*VertToEdgeListIter).EndV] = level+1;
+ }
+ if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level && _isOnBoundary[(*VertToEdgeListIter).StartV] == 0)
+ {
+ _isOnBoundary[(*VertToEdgeListIter).StartV] = level+1;
+ }
+ }
+ }
+ }//Loop over the level of the ring
+
+
+ // Check test to see the tag of the nodes on the boundary:
+// for (int i = 0; i< _EntityArray.size(); ++i)
+// {
+// GFace* face = _EntityArray[i]; //List of all the discrete face on which we compute the curvature
+
+// for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+// {
+// MElement *e = face->getMeshElement(iElem);
+
+// // The NEW tag of the corresponding element
+// const int E = _ElementToInt[e->getNum()];
+
+// // Pointers to vertices of triangle
+// MVertex* A = e->getVertex(0);
+// MVertex* B = e->getVertex(1);
+// MVertex* C = e->getVertex(2);
+
+// V[0] = _VertexToInt[A->getNum()]; //The new number of the vertex
+// V[1] = _VertexToInt[B->getNum()];
+// V[2] = _VertexToInt[C->getNum()];
+
+// if (_isOnBoundary[V[0]] == 3)
+// {
+// std::cout << "Vertex: " << A->getNum() << " -- " << A->x() << "; " << A->y() << "; " << A->z() << std::endl;
+// }
+// if (_isOnBoundary[V[1]] == 3)
+// {
+// std::cout << "Vertex: " << B->getNum() << " -- " << B->x() << "; " << B->y() << "; " << B->z() << std::endl;
+// }
+// if (_isOnBoundary[V[2]] == 3)
+// {
+// std::cout << "Vertex: " << C->getNum() << " -- " << C->x() << "; " << C->y() << "; " << C->z() << std::endl;
+// }
+// }
+// }
+
+
+ // Now we'll propagate the cuvature values from inside nodes with level = 3 close to the boundary - first
+ // to nodes with level 2, then from nodes with level 2 to nodes with level 1 (on the boundary)
+ _NbNeighbour.resize(_VertexToInt.size());
+ for(int i = 0; i < _NbNeighbour.size(); ++i)
+ {
+ _NbNeighbour[i] = 0;
+ }
+
+ for (int level = 2; level > 0 ; --level)
+ {
+
+ for (int i = 0; i< _NbNeighbour.size(); ++i)
+ {
+ _NbNeighbour[i] = 0;
+ if (_isOnBoundary[i] == level)
+ {
+ _VertexCurve[i] = 0.0;
+ }
+ }
+
+
+ for (int i = 0; i < _VertexToEdgeList.size(); ++i)
+ {
+ for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
+ ++VertToEdgeListIter)
+ {
+
+ if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level && _isOnBoundary[(*VertToEdgeListIter).EndV] == level+1)
+ {
+ _VertexCurve[(*VertToEdgeListIter).StartV] += _VertexCurve[(*VertToEdgeListIter).EndV];
+ _NbNeighbour[(*VertToEdgeListIter).StartV] ++;
+
+ }
+ if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level && _isOnBoundary[(*VertToEdgeListIter).StartV] == level+1)
+ {
+ _VertexCurve[(*VertToEdgeListIter).EndV] += _VertexCurve[(*VertToEdgeListIter).StartV];
+ _NbNeighbour[(*VertToEdgeListIter).EndV] ++;
+ }
+
+ }
+ }
+
+ // Correction for a degenerate case when a node has neighbours with the same or lower level, but zero
+ // neighbours with level+1
+
+ for (int i = 0; i < _VertexToEdgeList.size(); ++i)
+ {
+ for (VertToEdgeListIter = _VertexToEdgeList[i].begin(); VertToEdgeListIter != _VertexToEdgeList[i].end();
+ ++VertToEdgeListIter)
+ {
+
+ if (_isOnBoundary[(*VertToEdgeListIter).StartV] == level
+ && _isOnBoundary[(*VertToEdgeListIter).EndV] == level
+ && _NbNeighbour[(*VertToEdgeListIter).StartV] == 0)
+ {
+ _VertexCurve[(*VertToEdgeListIter).StartV] += _VertexCurve[(*VertToEdgeListIter).EndV];
+ _NbNeighbour[(*VertToEdgeListIter).StartV] = _NbNeighbour[(*VertToEdgeListIter).EndV];
+ }
+ if (_isOnBoundary[(*VertToEdgeListIter).EndV] == level
+ && _isOnBoundary[(*VertToEdgeListIter).StartV] == level
+ && _NbNeighbour[(*VertToEdgeListIter).EndV] == 0)
+ {
+ _VertexCurve[(*VertToEdgeListIter).EndV] += _VertexCurve[(*VertToEdgeListIter).StartV];
+ _NbNeighbour[(*VertToEdgeListIter).EndV] = _NbNeighbour[(*VertToEdgeListIter).StartV];
+ }
+ }
+ }
+
+ for (int i = 0; i < _isOnBoundary.size(); ++i)
+ {
+ if (_isOnBoundary[i] == level)
+ {
+ if (_NbNeighbour[i] == 0)
+ {
+ std::cout << "ERROR: VERTEX " << i << " WITH LEVEL " << level << " HAS 0 NEIGHBOURS" << std::endl;
+ }
+ _VertexCurve[i] = _VertexCurve[i]/_NbNeighbour[i];
+ }
+ }
+
+ }//Loop over the levels of the ring
+
+// for (int i = 0; i < _isOnBoundary.size(); ++i)
+// {
+// _VertexCurve[i] = _isOnBoundary[i];
+// }
+
+#endif
+ }
+
+//========================================================================================================
+
+//INITIALIZATION OF THE MAP AND RENUMBERING OF THE SELECTED ENTITIES:
-// initialization of the map and renumbering of the selected entities
void Curvature::initializeMap()
{
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
- GFace* face = _ptFinalEntityList[i];
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
- MElement *e = face->getMeshElement(iElem);
- const int E = e->getNum();
- _ElementToInt[E] = 1;
- const int A = e->getVertex(0)->getNum(); // pointer to 1st vertex of the triangle
- const int B = e->getVertex(1)->getNum();
- const int C = e->getVertex(2)->getNum();
- _VertexToInt[A] = 1;
- _VertexToInt[B] = 1;
- _VertexToInt[C] = 1;
+ if (! _isMapInitialized)
+ {
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i];
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
+ MElement *e = face->getMeshElement(iElem);
+ const int E = e->getNum();
+ _ElementToInt[E] = 1;
+
+ const int A = e->getVertex(0)->getNum(); //Pointer to 1st vertex of the triangle
+ const int B = e->getVertex(1)->getNum();
+ const int C = e->getVertex(2)->getNum();
+
+ _VertexToInt[A] = 1;
+ _VertexToInt[B] = 1;
+ _VertexToInt[C] = 1;
+ }
}
- }
- // set up a new numbering of chosen vertices and triangles
- int idx = 0;
+ /// Set up a new numbering of chosen vertices and triangles
+ int idx = 0;
- // map between the pointer to vertex and the new numbering of the vertex
- std::map<int,int>::iterator vertex_iterator;
- // map between the pointer to element and the new numbering of the element
- std::map<int,int>::iterator element_iterator;
+ // map between the pointer to vertex and the new numbering of the vertex
+ std::map<int,int>::iterator vertex_iterator;
+ // map between the pointer to element and the new numbering of the element
+ 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;
+
+ _isMapInitialized = true;
+
+ }
- 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;
-
}
-// compute the normal at the vertex and the area around
+//========================================================================================================
+
+//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
+
void Curvature::computeVertexNormals()
{
SVector3 vector_AB;
@@ -132,12 +386,14 @@ void Curvature::computeVertexNormals()
_VertexArea.resize(_ElementToInt.size() );
_VertexNormal.resize(_VertexToInt.size());
- for (unsigned i = 0; i < _ptFinalEntityList.size(); ++i){
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
// face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ GFace* face = _EntityArray[i];
//Loop over the element all the element of the "myTag"-surface
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
// Pointer to one element
MElement *e = face->getMeshElement(iElem);
// The NEW tag of the corresponding element
@@ -160,15 +416,13 @@ void Curvature::computeVertexNormals()
// Area of the triangles:
_TriangleArea[E] = 0.5*cross.norm();
- // std::cout << "The area of the triangle nr: " << e->getNum()
- // << " is: "<< TriangleArea[E] << std::endl;
+ // std::cout << "The area of the triangle nr: " << e->getNum() << " is: "<< TriangleArea[E] << std::endl;
_VertexArea[V0] += _TriangleArea[E];
_VertexArea[V1] += _TriangleArea[E];
_VertexArea[V2] += _TriangleArea[E];
- // here we are actually computing the unit normal vector per vertex
- _VertexNormal[V0] += cross;
+ _VertexNormal[V0] += cross; //here we are actually computing the unit normal vector per vertex
_VertexNormal[V1] += cross;
_VertexNormal[V2] += cross;
@@ -176,14 +430,18 @@ void Curvature::computeVertexNormals()
} //Loop over _ptFinalEntityList
- // Normalize the vertex-normals.
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
- _VertexNormal[n].normalize();
- }
+ ///////Normalize the vertex-normals.
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
+ {
+ _VertexNormal[n].normalize();
+ }
}
+//========================================================================================================
+
void Curvature::curvatureTensor()
{
+
STensor3 TempTensor;
STensor3 TijABTensorProduct;
SVector3 TijAB;
@@ -191,21 +449,23 @@ void Curvature::curvatureTensor()
_CurveTensor.resize(_VertexToInt.size());
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
- // face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"
- //Loop over the element all the element of the "myTag"-surface
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
MElement *e = face->getMeshElement(iElem); //Pointer to one element
const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
- // Loop over the 3 edges of each element
- for (unsigned int i = 0; i<3; ++i){
- MVertex* A = e->getVertex(i); //Pointer to 1st vertex of the edge A-to-B
- MVertex* B = e->getVertex((i+1)%3); //Pointer to 2nd vertex of the edge A-to-B
- const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the edge A-to-B
- const int V1 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the edge A-to-B
+ for (unsigned int i = 0; i<3; ++i) // Loop over the 3 edges of each element
+ {
+
+ MVertex* A = e->getVertex(i); //Pointer to 1st vertex of the edge A-to-B
+ MVertex* B = e->getVertex((i+1)%3); //Pointer to 2nd vertex of the edge A-to-B
+
+ const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the edge A-to-B
+ const int V1 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the edge A-to-B
//Weight for triangle-i-th's contribution to the shape tensor:
const double Wij0 = _TriangleArea[E] / (2 * _VertexArea[V0]);
@@ -213,7 +473,7 @@ void Curvature::curvatureTensor()
//Approximate Curvature "kappa" along some tangential vector T:
vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
- const double k_nominator0 = dot(_VertexNormal[V0], vector_AB);
+ const double k_nominator0 = dot(_VertexNormal[V0], vector_AB); //Dot-product of the 2 vectors
const double k_nominator1 = -dot(_VertexNormal[V1], vector_AB);
const double coeff = 2.0/vector_AB.normSq(); //normSq is the norm to the power 2
@@ -228,9 +488,11 @@ void Curvature::curvatureTensor()
TempTensor(1,1) += 1.0;
TempTensor(2,2) += 1.0;
- for (int m = 0; m<3; ++m){
+ for (int m = 0; m<3; ++m)
+ {
TijAB(m) = 0.0;
- for (int n = 0; n<3; ++n){
+ for (int n = 0; n<3; ++n)
+ {
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
@@ -246,9 +508,11 @@ void Curvature::curvatureTensor()
TempTensor(1,1) += 1.0;
TempTensor(2,2) += 1.0;
- for (int m = 0; m<3; ++m){
+ for (int m = 0; m<3; ++m)
+ {
TijAB(m) = 0.0;
- for (int n = 0; n<3; ++n){
+ for (int n = 0; n<3; ++n)
+ {
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
@@ -265,6 +529,8 @@ void Curvature::curvatureTensor()
}//End of method
+//========================================================================================================
+
void Curvature::computeCurvature_Simple()
{
SVector3 vector_E;
@@ -281,8 +547,8 @@ void Curvature::computeCurvature_Simple()
_VertexCurve.resize(_VertexToInt.size());
- //Loop over the vertex
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n) //Loop over the vertex
+ {
vector_E = SVector3(1,0,0);
vector_A = vector_E + _VertexNormal[n];
vector_B = vector_E - _VertexNormal[n];
@@ -290,10 +556,12 @@ void Curvature::computeCurvature_Simple()
const double MagA = vector_A.norm();
const double MagB = vector_B.norm();
- if (MagB > MagA){
+ if (MagB > MagA)
+ {
vector_Wvi = vector_B;
}
- else{
+ else
+ {
vector_Wvi = vector_A;
}
@@ -306,12 +574,14 @@ void Curvature::computeCurvature_Simple()
Qvi(2,2) += 1.0;
//Transpose the matrix:
- for (int i = 0; i<3; ++i){
- for (int j = 0; j<3; ++j){
+ for (int i = 0; i<3; ++i)
+ {
+ for (int j = 0; j<3; ++j)
+ {
QviT(i,j) = Qvi(j,i);
- }
- }
-
+ }
+ }
+
QviT *= _CurveTensor[n];
QviT *=Qvi;
Holder = QviT;
@@ -322,7 +592,8 @@ void Curvature::computeCurvature_Simple()
const double C = Holder(1,1)*Holder(2,2) - Holder(1,2)*Holder(2,1);
const double Delta = std::sqrt(B*B-4*A*C);
- if((B*B-4.*A*C) < 0.0){
+ if((B*B-4.*A*C) < 0.0)
+ {
std::cout << "WARNING: negative discriminant: " << B*B-4.*A*C << std::endl;
}
@@ -335,6 +606,10 @@ void Curvature::computeCurvature_Simple()
}
}
+//========================================================================================================
+
+//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
+
void Curvature::computeRusinkiewiczNormals()
{
SVector3 vector_AB;
@@ -344,12 +619,14 @@ void Curvature::computeRusinkiewiczNormals()
_TriangleArea.resize(_ElementToInt.size() );
_VertexNormal.resize(_VertexToInt.size());
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
// face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ GFace* face = _EntityArray[i];
//Loop over the element all the element of the "myTag"-surface
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
// Pointer to one element
MElement *e = face->getMeshElement(iElem);
const int E = _ElementToInt[e->getNum()];
@@ -359,7 +636,7 @@ void Curvature::computeRusinkiewiczNormals()
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V0 = _VertexToInt[A->getNum()]; //The new number of the vertex
+ const int V0 = _VertexToInt[A->getNum()]; //The new number of the vertex
const int V1 = _VertexToInt[B->getNum()];
const int V2 = _VertexToInt[C->getNum()];
@@ -383,30 +660,32 @@ void Curvature::computeRusinkiewiczNormals()
} // end of loop over elements of one face
} //Loop over _ptFinalEntityList
-
- // Normalize the vertex-normals.
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
- _VertexNormal[n].normalize();
- }
+
+ ///////Normalize the vertex-normals.
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
+ {
+ _VertexNormal[n].normalize();
+ }
}
+//========================================================================================================
// Compute per-vertex point areas
-void Curvature::computePointareas()
-{
+void Curvature::computePointareas(){
+
SVector3 e[3];
SVector3 l2;
SVector3 ew;
+ double factor[3];
_pointareas.resize(_VertexToInt.size());
_cornerareas.resize(_ElementToInt.size());
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i)
+ for (int i = 0; i< _EntityArray.size(); ++i)
{
- //face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"
- //Loop over the element all the element of the "myTag"-surface
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
MElement *E = face->getMeshElement(iElem); //Pointer to one element
// The NEW tag of the corresponding element
const int EIdx = _ElementToInt[E->getNum()];
@@ -415,8 +694,30 @@ void Curvature::computePointareas()
MVertex* B = E->getVertex(1);
MVertex* C = E->getVertex(2);
+ const int V0 = _VertexToInt[A->getNum()];
+ const int V1 = _VertexToInt[B->getNum()];
+ const int V2 = _VertexToInt[C->getNum()];
+
+ if (_isOnBoundary[V0])
+ {
+ factor[0] = 1.0;
+ }
+ else {factor[0] = 1.0;}
+
+ if (_isOnBoundary[V1])
+ {
+ factor[1] = 1.0;
+ }
+ else {factor[1] = 1.0;}
+
+ if (_isOnBoundary[V2])
+ {
+ factor[2] = 1.0;
+ }
+ else {factor[2] = 1.0;}
+
//Edges
- e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
+ e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangular element
e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
@@ -430,57 +731,62 @@ void Curvature::computePointareas()
l2[1] * (l2[2] + l2[0] - l2[1]),
l2[2] * (l2[0] + l2[1] - l2[2]) );
- if (ew[0] <= 0.0){
+ if (ew[0] <= 0.0)
+ {
_cornerareas[EIdx][1] = -0.25 * l2[2] * area / dot(e[0], e[2]);
_cornerareas[EIdx][2] = -0.25 * l2[1] * area / dot(e[0], e[1]);
_cornerareas[EIdx][0] = area - _cornerareas[EIdx][1] - _cornerareas[EIdx][2];
}
- else if (ew[1] <= 0.0){
+ else if (ew[1] <= 0.0)
+ {
_cornerareas[EIdx][2] = -0.25 * l2[0] * area / dot(e[1], e[0]);
_cornerareas[EIdx][0] = -0.25 * l2[2] * area / dot(e[1], e[2]);
_cornerareas[EIdx][1] = area - _cornerareas[EIdx][2] - _cornerareas[EIdx][0];
}
- else if (ew[2] <= 0.0){
+ else if (ew[2] <= 0.0)
+ {
_cornerareas[EIdx][0] = -0.25 * l2[1] * area / dot(e[2], e[1]);
_cornerareas[EIdx][1] = -0.25 * l2[0] * area / dot(e[2], e[0]);
_cornerareas[EIdx][2] = area - _cornerareas[EIdx][0] - _cornerareas[EIdx][1];
}
- else{
+ else
+ {
float ewscale = 0.5 * area / (ew[0] + ew[1] + ew[2]);
for (int j = 0; j < 3; j++)
_cornerareas[EIdx][j] = ewscale * (ew[(j+1)%3] + ew[(j+2)%3]);
}
- const int V0 = _VertexToInt[A->getNum()];
- const int V1 = _VertexToInt[B->getNum()];
- const int V2 = _VertexToInt[C->getNum()];
-
_pointareas[V0] += _cornerareas[EIdx][0];
- //_pointareas[faces[i][0]] += _cornerareas[i][0];
_pointareas[V1] += _cornerareas[EIdx][1];
_pointareas[V2] += _cornerareas[EIdx][2];
+ for (int j = 0; j < 3; j++)
+ {
+ _cornerareas[EIdx][j] = factor[j]*_cornerareas[EIdx][j];
+ }
+
+
} //End of loop over iElem
- // std::cout << "_pointareas.size = " << _pointareas.size() << std::endl;
- // std::cout << "_cornerareas.size = " << _cornerareas.size() << std::endl;
+// std::cout << "_pointareas.size = " << _pointareas.size() << std::endl;
+// std::cout << "_cornerareas.size = " << _cornerareas.size() << std::endl;
} //End of loop over _ptFinalEntityList
-}
+} //End of the method "computePointareas"
+
+//========================================================================================================
//Rotate a coordinate system to be perpendicular to the given normal
-void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
- const SVector3 &new_norm, SVector3 &new_u,
- SVector3 &new_v)
-{
+void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v){
+
new_u = old_u;
new_v = old_v;
SVector3 old_norm = crossprod(old_u, old_v);
double ndot = dot(old_norm, new_norm);
-// if (unlikely(ndot <= -1.0f))
+// if (unlikely(ndot <= -1.0f))
if (ndot <= -1.0f)
{
new_u = -1.0*new_u;
@@ -494,14 +800,16 @@ void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
new_v -= dperp*dot(new_v, perp_old);
}
-// Project a curvature tensor from the basis spanned by old_u and
-// old_v (which are assumed to be unit-length and perpendicular) to
-// the new_u and new_v basis
-void Curvature::proj_curv(const SVector3 &old_u, const SVector3 &old_v,
+//========================================================================================================
+
+//Project a curvature tensor from the basis spanned by old_u and old_v
+//(which are assumed to be unit-length and perpendicular) to the new_u
+//and new_v basis
+
+void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
double old_ku, double old_kuv, double old_kv,
const SVector3 &new_u, const SVector3 &new_v,
- double &new_ku, double &new_kuv, double &new_kv)
-{
+ double &new_ku, double &new_kuv, double &new_kv){
SVector3 r_new_u;
SVector3 r_new_v;
rot_coord_sys(new_u, new_v, crossprod(old_u,old_v), r_new_u, r_new_v);
@@ -516,14 +824,16 @@ void Curvature::proj_curv(const SVector3 &old_u, const SVector3 &old_v,
new_kv = old_ku*u2*u2 + old_kuv*(2.0f * u2*v2) + old_kv*v2*v2;
}
-// Given a curvature tensor, find principal directions and curvatures
-// Makes sure that pdir1 and pdir2 are perpendicular to normal
+
+//========================================================================================================
+
+//Given a curvature tensor, find principal directions and curvatures
+//Makes sure that pdir1 and pdir2 are perpendicular to normal
+
void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
- double ku, double kuv, double kv,
- const SVector3 &new_norm,
- SVector3 &pdir1, SVector3 &pdir2,
- double &k1, double &k2)
-{
+ double ku, double kuv, double kv,
+ const SVector3 &new_norm,
+ SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2){
SVector3 r_old_u;
SVector3 r_old_v;
@@ -533,7 +843,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(unlikely(kuv !=0.0f))
if(kuv !=0.0) {
//Jacobi rotation to diagonalize
double h= 0.5*(kv -ku)/ kuv;
@@ -547,18 +857,20 @@ 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;
}
pdir2 = crossprod(new_norm, pdir1);
}
+//========================================================================================================
void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
{
+
_model = model;
double t0 = Cpu();
@@ -575,12 +887,153 @@ void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
writeToPosFile("curvature.pos");
writeToVtkFile("curvature.vtk");
+
}
+//========================================================================================================
+
+void Curvature::buildEdgeList()
+{
+
+ int V[3];
+ bool edgefound;
+
+ _VertexToEdgeList.clear();
+
+ std::list<MeshEdgeInfo>::iterator VertToEdgeListIter;
+ MeshEdgeInfo TempEdge;
+
+ _VertexToEdgeList.resize(_VertexToInt.size());
+
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i]; //List of all the discrete face on which we compute the curvature
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
+ MElement *e = face->getMeshElement(iElem);
+
+ // The NEW tag of the corresponding element
+ const int E = _ElementToInt[e->getNum()];
+
+ // Pointers to vertices of triangle
+ MVertex* A = e->getVertex(0);
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ V[0] = _VertexToInt[A->getNum()]; //The new number of the vertex
+ V[1] = _VertexToInt[B->getNum()];
+ V[2] = _VertexToInt[C->getNum()];
+
+ // Try to add edge [V0,V1] to the global list
+ for (int j = 0; j < 3; ++j)
+ {
+ const int StartV = std::min(V[j], V[(j+1)%3]);
+ const int EndV = std::max(V[j], V[(j+1)%3]);
+
+ edgefound = false;
+ for (VertToEdgeListIter = _VertexToEdgeList[StartV].begin(); VertToEdgeListIter != _VertexToEdgeList[StartV].end();
+ ++VertToEdgeListIter)
+ {
+ if(StartV == (*VertToEdgeListIter).StartV && EndV == (*VertToEdgeListIter).EndV)
+ {
+ (*VertToEdgeListIter).NbElemNeighbour ++;
+ edgefound = true;
+ }
+
+ }
+ if (false == edgefound)
+ {
+ TempEdge.StartV = StartV;
+ TempEdge.EndV = EndV;
+ TempEdge.NbElemNeighbour = 1;
+ _VertexToEdgeList[StartV].push_back(TempEdge);
+ }
+
+
+ } // Loop over j
+
+
+ } // Loop over the elements (triangles) of the face
+ }
+
+ int NbEdges = 0;
+ for(int i = 0; i < _VertexToEdgeList.size(); ++i)
+ {
+ NbEdges += _VertexToEdgeList[i].size();
+ }
+
+ std::cout << "Euler formula:" << std::endl;
+ std::cout << "Edges + 2 = " << NbEdges + 2 << std::endl;
+ std::cout << "Elements + Nodes = " << _VertexToInt.size() + _ElementToInt.size() << std::endl;
+
+}
+
+//========================================================================================================
+
+void Curvature::smoothCurvatureField(const int NbIter)
+{
+ if ( _VertexToEdgeList.size() == 0 ) { buildEdgeList(); }
+
+ std::vector<double> smoothedCurvature;
+ smoothedCurvature.resize( _VertexToInt.size() );
+
+ _NbNeighbour.resize(_VertexToInt.size());
+ for(int i = 0; i < _VertexToInt.size(); ++i)
+ {
+ _NbNeighbour[i] = 0;
+ }
+
+
+ // Smoothing iterations
+ for(int iter = 0; iter < NbIter; ++iter)
+ {
+
+ for(int i = 0; i < smoothedCurvature.size(); ++i)
+ {
+ smoothedCurvature[i] = 0.0;
+ }
+
+ std::list<MeshEdgeInfo>::const_iterator edgeIterator;
+
+
+ for(int i = 0; i < _VertexToEdgeList.size(); ++i)
+ {
+ for(edgeIterator = _VertexToEdgeList[i].begin(); edgeIterator != _VertexToEdgeList[i].end(); ++edgeIterator)
+ {
+ const int N0 = (*edgeIterator).StartV;
+ const int N1 = (*edgeIterator).EndV;
+
+ const int V0 = _VertexToInt[N0];
+ const int V1 = _VertexToInt[N1];
+
+ smoothedCurvature[V0] += _VertexCurve[V1];
+ smoothedCurvature[V1] += _VertexCurve[V0];
+ _NbNeighbour[V0]++;
+ _NbNeighbour[V1]++;
+ }
+ }
+
+
+
+ const double Lambda = 0.3;
+ for(int i = 0; i < _VertexCurve.size(); ++i)
+ {
+ _VertexCurve[i] = Lambda*_VertexCurve[i] + (1-Lambda)*smoothedCurvature[i] / _NbNeighbour[i];
+ }
+ }
+
+}
+
+//========================================================================================================
+
+
+
void Curvature::computeCurvature_Rusinkiewicz(int isMax)
{
retrieveCompounds();
initializeMap();
+
computeRusinkiewiczNormals();
computePointareas();
@@ -592,8 +1045,6 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
SVector3 m;
SVector3 dn;
- int vj;
-
double u;
double v;
double dnu;
@@ -610,43 +1061,46 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
_curv2.resize(_VertexToInt.size());
_curv12.resize(_VertexToInt.size());
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ for (int i = 0; i< _EntityArray.size(); ++i)
{
- //face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"
- //Loop over the element all the element of the "myTag"-surface
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
- MElement *E = face->getMeshElement(iElem); // Pointer to one element
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
+ MElement *E = face->getMeshElement(iElem); //Pointer to one element
- MVertex* A = E->getVertex(0); // Pointers to vertices of triangle
+ MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
MVertex* B = E->getVertex(1);
MVertex* C = E->getVertex(2);
- const int V0 = _VertexToInt[A->getNum()]; // Tag of the 1st vertex of the triangle
+ const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
const int V1 = _VertexToInt[B->getNum()];
const int V2 = _VertexToInt[C->getNum()];
- //Set up an initial coordinate system per vertex:
+ ///Set up an initial coordinate system per vertex:
+
_pdir1[V0] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
_pdir1[V1] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
_pdir1[V2] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
}
}
- for (unsigned ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
+ {
_pdir1[ivertex] = crossprod(_pdir1[ivertex], _VertexNormal[ivertex]);
_pdir1[ivertex].normalize();
_pdir2[ivertex] = crossprod(_VertexNormal[ivertex], _pdir1[ivertex]);
}
// Compute curvature per face:
- for (int i = 0; i< _ptFinalEntityList.size(); ++i){
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
//face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ GFace* face = _EntityArray[i];
//Loop over all elements of this face:
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
MElement *E = face->getMeshElement(iElem); //Pointer to one element
// The NEW tag of the corresponding element
const int EIdx = _ElementToInt[E->getNum()];
@@ -659,31 +1113,31 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
- // SVector3 e[3] = {SVector3(C->x() - B->x(), C->y() - B->y(),
- // C->z() - B->z()), SVector3(A->x() - C->x(), A->y() - C->y(),
- // A->z() - C->z()), SVector3(B->x() - A->x(), B->y() - A->y(),
- // B->z() - A->z()) };
+ //SVector3 e[3] = {SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()), SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z()), SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z()) };
- // N-T-B coordinate system per face
+ //N-T-B coordinate system per face
t = e[0];
t.normalize();
n = crossprod( e[0], e[1]);
b = crossprod(n, t);
b.normalize();
- // Estimate curvature based on variations of normals along
- // edges: intialization:
+ //Estimate curvature based on variations of normals along edges:
+ //intialization:
m = SVector3(0.0, 0.0, 0.0);
//maybe double m[3] = { 0.0, 0.0, 0.0 };
// w *= 0.0; //Reset w to zero
- for (int i = 0; i< 3; ++i){
- for (int j = 0; j<3; ++j){
+ for (int i = 0; i< 3; ++i)
+ {
+ for (int j = 0; j<3; ++j)
+ {
w(i,j) = 0.0;
}
}
//filling:
- for (int j = 0; j< 3; ++j){
+ for (int j = 0; j< 3; ++j)
+ {
u = dot(e[j], t);
v = dot(e[j], b);
@@ -710,21 +1164,27 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
w(1,1) = w(0,0) + w(2,2);
w(1,2) = w(0,1);
- // Least Squares Solution
+ //Least Squares Solution
double diag[3];
- if (!ldltdc(w, diag)){
+ if (!ldltdc(w, diag))
+ {
std::cout << "ldltdc failed" << std::endl;
continue;
}
ldltsl(w, diag, m, m);
- // Push it back out to the vertices
- for (int j = 0; j< 3; ++j){
+ //Push it back out to the vertices
+ for (int j = 0; j< 3; ++j)
+ {
const int old_vj = E->getVertex(j)->getNum();
const int vj = _VertexToInt[old_vj];
proj_curv(t, b, m[0], m[1], m[2], _pdir1[vj], _pdir2[vj], c1, c12, c2);
wt = _cornerareas[EIdx][j]/_pointareas[vj];
- // wt = 1.0;
+// wt = 1.0;
+// if (_isOnBoundary[vj])
+// {
+// wt = 2.0*wt;
+// }
_curv1[vj] += wt*c1;
_curv12[vj] += wt*c12;
@@ -733,19 +1193,19 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
} //End of loop over the element (iElem)
- } //End of loop over "_ptFinalEntityList"
+ } //End of loop over "_EntityArray"
+
//Compute principal directions and curvatures at each vertex
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) {
- diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex],
- _curv12[ivertex], _curv2[ivertex],
- _VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex],
- _curv1[ivertex], _curv2[ivertex]);
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) {
+ diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv12[ivertex], _curv2[ivertex],
+ _VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv2[ivertex]);
}
-
- _VertexCurve.resize(_VertexToInt.size());
+
+ _VertexCurve.resize( _VertexToInt.size() );
for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){
+
if (isMax){
_VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
}
@@ -757,19 +1217,26 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
}
}
+
+// smoothCurvatureField(1);
_alreadyComputedCurvature = true;
-}
+// Propagate the value of curvature from nodes close the edges of the geometry onto the edges
+ correctOnEdges();
+
+} //End of the "computeCurvature_Rusinkiewicz" method
+
+ //========================================================================================================
void Curvature::computeCurvature_RBF()
{
retrieveCompounds();
initializeMap();
- // fill set of MVertex
+ //fill set of MVertex
std::set<MVertex*> allNodes;
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i) {
- GFaceCompound* face = (GFaceCompound*)_ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i) {
+ GFaceCompound* face = (GFaceCompound*)_EntityArray[i];
for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++) {
MElement *e = face->getMeshElement(iElem);
for(int j = 0; j < e->getNumVertices(); j++){
@@ -778,7 +1245,7 @@ void Curvature::computeCurvature_RBF()
}
}
- // bounding box
+ //bounding box
SBoundingBox3d bb;
std::vector<SPoint3> vertices;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -788,18 +1255,18 @@ void Curvature::computeCurvature_RBF()
}
double sizeBox = norm(SVector3(bb.max(), bb.min()));
- // compure curvature RBF
+ //compure curvature RBF
std::map<MVertex*, SVector3> _normals;
std::vector<MVertex*> _ordered;
std::map<MVertex*, double> curvRBF;
- // GLOBAL
+ //GLOBAL
GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered);
_rbf->computeCurvature(_rbf->getXYZ(),curvRBF);
- // LOCAL FD
- // GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true);
- // _rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
+ //LOCAL FD
+ //GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true);
+ //_rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
- // fill vertex curve
+ //fill vertex curve
_VertexCurve.resize( _VertexToInt.size() );
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
@@ -808,58 +1275,57 @@ void Curvature::computeCurvature_RBF()
vertexIterator = _VertexToInt.find( v->getNum() );
if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
_VertexCurve[V0] = curvRBF[v];
- }
-
+ }
+
_alreadyComputedCurvature = true;
-}
+} //End of the "computeCurvature_RBF" method
+
+
+ //========================================================================================================
+
+void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs)
+ {
+ 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;
+
+ 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
+ }
-void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1,
- double& c2, int isAbs)
-{
- 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;
-
- 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
}
-
-}
-void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax,
- SVector3* dMin, double* cMax,
- double* cMin, int isAbs)
+//========================================================================================================
+
+void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs)
{
MVertex* A = triangle->getVertex(0);
MVertex* B = triangle->getVertex(1);
@@ -871,22 +1337,19 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
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;
+ 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;
+ 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;
+ std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
if (isAbs){
dMax[0] = _pdir1[V0];
@@ -904,8 +1367,10 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
cMin[0] = std::abs(_curv2[V0]);
cMin[1] = std::abs(_curv2[V1]);
cMin[2] = std::abs(_curv2[V2]);
+
}
else{
+
dMax[0] = _pdir1[V0];
dMax[1] = _pdir1[V1];
dMax[2] = _pdir1[V2];
@@ -924,42 +1389,45 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
}
}
+
+
+ //========================================================================================================
+
void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
-{
- 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;
+ {
+ 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;
-
- if (isAbs){
- c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
- c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
- }
- else{
- c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
- c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
- }
-}
+ 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;
+
+ if (isAbs){
+ c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+ c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
+ }
+ else{
+ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+ c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
+ }
-double Curvature::getAtVertex(const MVertex *v) const
-{
+ }
+
+//========================================================================================================
+
+double Curvature::getAtVertex(const MVertex *v) const {
std::map<int,int>::const_iterator it = _VertexToInt.find(v->getNum());
if (it == _VertexToInt.end()) {
- Msg::Error("curvature has not been computed for vertex %i (%i)",
- v->getNum(), _VertexToInt.size());
+ Msg::Error("curvature has not been computed for vertex %i (%i)", v->getNum(), _VertexToInt.size());
return 1;
}
return _VertexCurve[it->second];
@@ -971,11 +1439,11 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile.precision(18);
outfile.open(filename.c_str());
outfile << "View \"Curvature \"{" << std::endl;
-
+
int idxelem = 0;
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i) {
- GFace* face = _ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i) {
+ GFace* face = _EntityArray[i];
for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem);
@@ -986,11 +1454,12 @@ void Curvature::writeToPosFile( const std::string & filename)
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
- const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
- const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
-
- // Here is printing the triplet X-Y-Z of each vertex:
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+
+ //Here is printing the triplet X-Y-Z of each vertex:
+ //*************************************************
outfile << "ST("; //VT = vector triangles //ST = scalar triangle
outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
@@ -1000,31 +1469,33 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile << ")";
outfile <<"{";
- // Here is printing the 3 components of the normal vector for each vertex:
+ //Here is printing the 3 components of the normal vector for each vertex:
+ //**********************************************************************
+
+// outfile << _VertexNormal[V1].x() << ","<< _VertexNormal[V1].y() << ","<< _VertexNormal[V1].z() << ",";
+// outfile << _VertexNormal[V2].x() << ","<< _VertexNormal[V2].y() << ","<< _VertexNormal[V2].z() << ",";
+// outfile << _VertexNormal[V3].x() << ","<< _VertexNormal[V3].y() << ","<< _VertexNormal[V3].z();
- //outfile << VertexNormal[V1].x() << ","<< VertexNormal[V1].y() << ","
- // << VertexNormal[V1].z() << ",";
- //outfile << VertexNormal[V2].x() << ","<< VertexNormal[V2].y() << ","
- // << VertexNormal[V2].z() << ",";
- //outfile << VertexNormal[V3].x() << ","<< VertexNormal[V3].y() << ","
- // << VertexNormal[V3].z();
-
outfile << _VertexCurve[V1] << "," << _VertexCurve[V2] << "," << _VertexCurve[V3];
outfile << "};" << std::endl;
idxelem++;
- }
- }
+ } //Loop over elements
+
+} // Loop over ptFinalEntityList
- outfile << "};" << std::endl;
-
- outfile.close();
+outfile << "};" << std::endl;
+
+outfile.close();
}
+//========================================================================================================
+
void Curvature::writeToVtkFile( const std::string & filename)
{
+
std::ofstream outfile;
outfile.precision(18);
outfile.open(filename.c_str());
@@ -1037,20 +1508,22 @@ void Curvature::writeToVtkFile( const std::string & filename)
outfile << "POINTS " << npoints << " double" << std::endl;
- // Build a table of coordinates
- // - Loop over all elements and look at the 'old' (not necessarily
- // continuous) numbers of vertices
- // - Get the 'new' index of each vertex through _VertexToInt and the
- // [x,y,z] coordinates of this vertex
- // - Store them in coordx,coordy and coordz
+ /// Build a table of coordinates
+ /// Loop over all elements and look at the 'old' (not necessarily continuous) numbers of vertices
+ /// Get the 'new' index of each vertex through _VertexToInt and the [x,y,z] coordinates of this vertex
+ /// Store them in coordx,coordy and coordz
+
std::vector<VtkPoint> coord;
+
coord.resize(npoints);
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
- GFace* face = _ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
MElement *e = face->getMeshElement(iElem);
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
@@ -1079,23 +1552,25 @@ void Curvature::writeToVtkFile( const std::string & filename)
}
}
- for (int v = 0; v < npoints; ++v){
+ for (int v = 0; v < npoints; ++v)
+ {
outfile << coord[v].x << " " << coord[v].y << " " << coord[v].z << std::endl;
}
- // Empty the array 'coord' to free the memory
- // Point coordinates will not be needed anymore
+ /// Empty the array 'coord' to free the memory
+ /// Point coordinates will not be needed anymore
coord.clear();
- // Write the cell connectivity
+ /// Write the cell connectivity
- outfile << std::endl << "CELLS " << _ElementToInt.size() << " "
- << 4*_ElementToInt.size() << std::endl;
+ outfile << std::endl << "CELLS " << _ElementToInt.size() << " " << 4*_ElementToInt.size() << std::endl;
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
- GFace* face = _ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
+ {
MElement *e = face->getMeshElement(iElem);
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
@@ -1113,25 +1588,32 @@ void Curvature::writeToVtkFile( const std::string & filename)
outfile << "3 " << newIdxA << " " << newIdxB << " " << newIdxC << std::endl;
}
}
-
+
outfile << std::endl << "CELL_TYPES " << _ElementToInt.size() << std::endl;
- for(int ie = 0; ie < _ElementToInt.size(); ++ie){
+ for(int ie = 0; ie < _ElementToInt.size(); ++ie)
+ {
outfile << "5" << std::endl; //Triangle is element type 5 in vtk
+
}
- // Write the curvature values as vtk 'point data'
+ /// Write the curvature values as vtk 'point data'
outfile << std::endl << "POINT_DATA " << npoints << std::endl;
outfile << "SCALARS curvature float 1" << std::endl;
outfile << "LOOKUP_TABLE default" << std::endl;
- for (int iv = 0; iv < npoints; ++iv){
+ for (int iv = 0; iv < npoints; ++iv)
+ {
outfile << _VertexCurve[iv] << std::endl;
}
-
+
outfile.close();
+
+
}
+//========================================================================================================
+
void Curvature::writeDirectionsToPosFile( const std::string & filename)
{
std::ofstream outfile;
@@ -1139,26 +1621,28 @@ void Curvature::writeDirectionsToPosFile( const std::string & filename)
outfile.open(filename.c_str());
outfile << "View \"Curvature_DirMax \"{" << std::endl;
- for (unsigned i = 0; i< _ptFinalEntityList.size(); ++i){
- //face is a pointer to one surface of the group "FinalEntityList"
- GFace* face = _ptFinalEntityList[i];
+ for (int i = 0; i< _EntityArray.size(); ++i)
+ {
+ GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"
- //Loop over the element all the element of the "myTag"-surface
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
MElement *e = face->getMeshElement(iElem); //Pointer to one element
const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
-
+
//std::cout << "We are now looking at element Nr: " << E << std::endl;
-
+
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
- const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
- const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
//Here is printing the triplet X-Y-Z of each vertex:
+ //*************************************************
+
outfile << "VT("; //VT = vector triangles //ST = scalar triangle
outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
@@ -1168,54 +1652,75 @@ void Curvature::writeDirectionsToPosFile( const std::string & filename)
outfile <<"{";
//Here is printing the 3 components of the curvature max direction for each vertex:
- outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
- outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
- outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();
+ //*********************************************************************************
+
+ outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
+ outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
+ outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();
+
outfile << "};" << std::endl;
- }
- }
-
- outfile << "};" << std::endl;
-
- outfile << "View \"Curvature_DirMin \"{" << std::endl;
-
- for (int i = 0; i< _ptFinalEntityList.size(); ++i){
- GFace* face = _ptFinalEntityList[i];
-
- for (unsigned iElem = 0; iElem < face->getNumMeshElements(); iElem++){
- MElement *e = face->getMeshElement(iElem); //Pointer to one element
- const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
-
- //std::cout << "We are now looking at element Nr: " << E << std::endl;
-
- MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
- MVertex* B = e->getVertex(1);
- MVertex* C = e->getVertex(2);
-
- const int V1 = _VertexToInt[A->getNum()];
- const int V2 = _VertexToInt[B->getNum()];
- const int V3 = _VertexToInt[C->getNum()];
-
- //Here is printing the triplet X-Y-Z of each vertex:
- outfile << "VT("; //VT = vector triangles //ST = scalar triangle
- outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
- outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
- outfile << C->x() << ","<< C->y() << "," << C->z();
-
- outfile << ")";
- outfile <<"{";
-
- //Here is printing the 3 components of the curvature min direction for each vertex:
- outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
- outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
- outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();
-
- outfile << "};" << std::endl;
- }
- }
-
- outfile << "};" << std::endl;
-
- outfile.close();
+
+ } //Loop over elements
+
+} // Loop over ptFinalEntityList
+
+outfile << "};" << std::endl;
+
+
+//----------------------------------------------------------------------------------------------
+
+outfile << "View \"Curvature_DirMin \"{" << std::endl;
+
+for (int i = 0; i< _EntityArray.size(); ++i)
+{
+ GFace* face = _EntityArray[i]; //face is a pointer to one surface of the group "FinalEntityList"
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
+ MElement *e = face->getMeshElement(iElem); //Pointer to one element
+ const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
+
+ //std::cout << "We are now looking at element Nr: " << E << std::endl;
+
+ MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+
+ //Here is printing the triplet X-Y-Z of each vertex:
+ //*************************************************
+
+ outfile << "VT("; //VT = vector triangles //ST = scalar triangle
+ outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
+ outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
+ outfile << C->x() << ","<< C->y() << "," << C->z();
+
+ outfile << ")";
+ outfile <<"{";
+
+ //Here is printing the 3 components of the curvature min direction for each vertex:
+ //*********************************************************************************
+
+ outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
+ outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
+ outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();
+
+
+ outfile << "};" << std::endl;
+
+} //Loop over elements
+
+} // Loop over ptFinalEntityList
+
+outfile << "};" << std::endl;
+
+outfile.close();
}
+//========================================================================================================
+
+
+
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index c34a909ce87d476f31e5121f02cc0c63a07ec547..095d23f6bd694123d10d4826bc3f3f79ca724f35 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -3,160 +3,225 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _CURVATURE_H_
+#ifndef _CURVATUREL_H_
#define _CURVATURE_H_
-#include <map>
-#include <vector>
#include "GModel.h"
#include "STensor3.h"
+#include<map>
+#include<vector>
+#include<list>
+
class Curvature {
- private:
- typedef std::vector<GFace*> GFaceList;
- typedef std::map<int, GFaceList >::iterator GFaceIter;
- // helper type for writing VTK files
- struct VtkPoint
- {
- double x;
- double y;
- double z;
- };
-
- // static member variable to implement the class curvature as a Singleton
- static Curvature *_instance;
- static bool _destroyed;
+private:
+ //----------------------------------------
+ //TYPEDEFS:
- // boolean to check if the curvature has already been computed
- static bool _alreadyComputedCurvature;
+ typedef std::vector<GFace*> GFaceList;
+ //typedef std::map<int, GEntityList >::iterator GEntityIter;
+ typedef std::map<int, GFaceList >::iterator GFaceIter;
- // map definition
- std::map<int, int> _VertexToInt;
- std::map<int, int> _ElementToInt;
+ //-----------------------------------------
+ // HELPER TYPE FOR WRITING TO VTK FILES
+ struct VtkPoint
+ {
+ double x;
+ double y;
+ double z;
+ };
- // model and list of selected entities with give physical tag:
- GModel* _model;
- GFaceList _ptFinalEntityList;
+ //-----------------------------------------
+ // HELPER TYPE FOR DETECTING BOUNDARY NODES
+ struct MeshEdgeInfo
+ {
+ int StartV;
+ int EndV;
+ int NbElemNeighbour;
+ };
- // averaged vertex normals
- std::vector<SVector3> _VertexNormal;
+ //-----------------------------------------
+ // MEMBER VARIABLES
- // vector of principal dircections
- std::vector<SVector3> _pdir1;
- std::vector<SVector3> _pdir2;
+ //Static member variable to implement the class curvature as a Singleton
+ static Curvature *_instance;
+ static bool _destroyed;
- // vector of principal curvature
- std::vector<double> _curv1;
- std::vector<double> _curv2;
- std::vector<double> _curv12;
-
- // area around point
- std::vector<double> _pointareas;
- std::vector<SVector3> _cornerareas;
-
- // curvature Tensor per mesh vertex
- std::vector<STensor3> _CurveTensor;
-
- // area of a triangle element:
- std::vector<double> _TriangleArea;
-
- // area around a mesh vertex:
- std::vector<double> _VertexArea;
- std::vector<double> _VertexCurve;
-
- // constructor and destructor
- Curvature(){}
- ~Curvature();
- static void create();
- static void onDeadReference();
- void initializeMap();
- void computeVertexNormals();
- void curvatureTensor();
- static void rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
- const SVector3 &new_norm, SVector3 &new_u,
- SVector3 &new_v);
- void proj_curv(const SVector3 &old_u, const SVector3 &old_v, double old_ku,
- double old_kuv, double old_kv, const SVector3 &new_u,
- const SVector3 &new_v, double &new_ku, double &new_kuv,
- double &new_kv);
- void diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
- double ku, double kuv, double kv,
- const SVector3 &new_norm,
- SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2);
- void computePointareas();
- void computeRusinkiewiczNormals();
-
- // Perform LDL^T decomposition of a symmetric positive definite matrix.
- // Like Cholesky, but no square roots. Overwrites lower triangle of matrix.
- static inline bool ldltdc(STensor3& A, double rdiag[3])
- {
- double v[2];
- for (int i = 0; i < 3; i++){
- for (int k = 0; k < i; k++)
- v[k] = A(i,k) * rdiag[k];
- for (int j = i; j < 3; j++){
- double sum = A(i,j);
+ //Boolean to check if the curvature has already been computed
+ static bool _alreadyComputedCurvature;
+
+ //Map definition
+ std::map<int, int> _VertexToInt;
+ std::map<int, int> _ElementToInt;
+
+ bool _isMapInitialized;
+
+ //Model and list of selected entities with give physical tag:
+ GModel* _model;
+ GFaceList _EntityArray;
+
+ //Averaged vertex normals
+ std::vector<SVector3> _VertexNormal;
+
+ // Vector of principal dircections
+ std::vector<SVector3> _pdir1;
+ std::vector<SVector3> _pdir2;
+
+ // Vector of principal curvature
+ std::vector<double> _curv1;
+ std::vector<double> _curv2;
+ std::vector<double> _curv12;
+
+ // Area around point
+ std::vector<double> _pointareas;
+ std::vector<SVector3> _cornerareas;
+
+ //Curvature Tensor per mesh vertex
+ std::vector<STensor3> _CurveTensor;
+
+ //Area of a triangle element:
+ std::vector<double> _TriangleArea;
+
+ //Area around a mesh vertex:
+ std::vector<double> _VertexArea;
+ std::vector<double> _VertexCurve;
+
+
+ //Number of neighbour for each vertices in the mesh:
+ std::vector<int> _NbNeighbour;
+
+ //Vector of 0/1 to check which nodes are on the boundary:
+ std::vector<int> _isOnBoundary;
+
+ std::vector<std::list<MeshEdgeInfo> > _VertexToEdgeList;
+ //-----------------------------------------
+ // PRIVATE METHODS
+
+ //Constructor and destructor
+ Curvature();
+ ~Curvature();
+
+ static void create();
+ static void onDeadReference();
+
+ void initializeMap();
+ void computeVertexNormals();
+ void curvatureTensor();
+ static void rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
+ const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v);
+ void proj_curv( const SVector3 &old_u, const SVector3 &old_v, double old_ku, double old_kuv,
+ double old_kv, const SVector3 &new_u, const SVector3 &new_v,
+ double &new_ku, double &new_kuv, double &new_kv);
+ void diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
+ double ku, double kuv, double kv,
+ const SVector3 &new_norm,
+ SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2);
+ void computePointareas();
+ void computeRusinkiewiczNormals();
+
+ // Perform LDL^T decomposition of a symmetric positive definite matrix.
+ // Like Cholesky, but no square roots. Overwrites lower triangle of matrix.
+ static inline bool ldltdc(STensor3& A, double rdiag[3])
+ {
+ double v[2];
+ for (int i = 0; i < 3; i++)
+ {
for (int k = 0; k < i; k++)
- sum -= v[k] * A(j,k);
- if (i == j){
- //if (unlikely(sum <= T(0)))
- if (sum <= 0.0)
- return false;
- rdiag[i] = 1.0 / sum;
- }
- else{
- A(j,i) = sum;
+ v[k] = A(i,k) * rdiag[k];
+ for (int j = i; j < 3; j++)
+ {
+ double sum = A(i,j);
+ for (int k = 0; k < i; k++)
+ sum -= v[k] * A(j,k);
+ if (i == j)
+ {
+ //if (unlikely(sum <= T(0)))
+ if (sum <= 0.0)
+ return false;
+ rdiag[i] = 1.0 / sum;
+ }
+ else
+ {
+ A(j,i) = sum;
+ }
}
}
+
+ return true;
}
- return true;
- }
- // Solve Ax=B after ldltdc
- static inline void ldltsl(STensor3& A, double rdiag[3], double B[3], double x[3])
- {
- int i;
- for (i = 0; i < 3; i++){
- double sum = B[i];
- for (int k = 0; k < i; k++)
- sum -= A(i,k) * x[k];
- x[i] = sum * rdiag[i];
- }
- for (i = 3 - 1; i >= 0; i--){
- double sum = 0;
- for (int k = i + 1; k < 3; k++)
- sum += A(k,i) * x[k];
- x[i] -= sum * rdiag[i];
+ // Solve Ax=B after ldltdc
+ static inline void ldltsl(STensor3& A, double rdiag[3], double B[3], double x[3])
+ {
+ int i;
+ for (i = 0; i < 3; i++)
+ {
+ double sum = B[i];
+ for (int k = 0; k < i; k++)
+ sum -= A(i,k) * x[k];
+ x[i] = sum * rdiag[i];
+ }
+ for (i = 3 - 1; i >= 0; i--)
+ {
+ double sum = 0;
+ for (int k = i + 1; k < 3; k++)
+ sum += A(k,i) * x[k];
+ x[i] -= sum * rdiag[i];
+ }
}
- }
-
- public:
- typedef enum {RUSIN=1, RBF=2, SIMPLE=3} typeOfCurvature;
+
+ //Build the list of all edges:
+ void buildEdgeList();
+
+ //Smooth the computed curvature field via Laplacian smoothing method:
+ void smoothCurvatureField(const int NbIter);
+
+public:
+
+ typedef enum {RUSIN=1,RBF=2, SIMPLE=3} typeOfCurvature;
static Curvature& getInstance();
static bool valueAlreadyComputed();
+
+ inline void setGModel(GModel* model)
+ {
+ _model = model;
+ }
+
void retrieveCompounds();
+ void correctOnEdges();
double getAtVertex(const MVertex *v) const;
+ //void retrievePhysicalSurfaces(const std::string & face_tag);
+
void computeCurvature(GModel* model, typeOfCurvature typ);
+
/// The following function implements algorithm from:
/// Implementation of an Algorithm for Approximating the Curvature Tensor
/// on a Triangular Surface Mesh in the Vish Environment
/// Edwin Matthews, Werner Benger, Marcel Ritter
void computeCurvature_Simple();
+
/// The following function implements algorithm from:
/// Estimating Curvatures and Their Derivatives on Triangle Meshes
/// Szymon Rusinkiewicz, Princeton University
/// Code taken from Rusinkiewicz' 'trimesh2' library
void computeCurvature_Rusinkiewicz(int isMax=0);
+
void computeCurvature_RBF();
- void triangleNodalValues(MTriangle* triangle, double& c0, double& c1,
- double& c2, int isAbs=0);
- void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax,
- SVector3* dMin, double* cMax, double* cMin,
- int isAbs=0);
+
+ void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs=0);
+ void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs=0);
+
void edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs=0);
+
void writeToPosFile( const std::string & filename);
+
void writeToVtkFile( const std::string & filename);
+
void writeDirectionsToPosFile( const std::string & filename);
+
+
+
};
+
#endif