diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index f54d09e2d3afbb3b4e89bc647005ed380935340c..5b35cfa4316d849b523784f4416aa010de6803eb 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -51,6 +51,7 @@
#include "GPoint.h"
#include "GmshDefines.h"
#include "JacobianBasis.h"
+ #include "Curvature.h"
#if defined(HAVE_FLTK)
#include "FlGui.h"
#endif
@@ -146,6 +147,7 @@ namespace std {
%include "Generator.h"
%include "GmshDefines.h"
%include "JacobianBasis.h"
+%include "Curvature.h"
#if defined(HAVE_FLTK)
%include "FlGui.h"
#endif
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index a4c44720bbaecbf0116185a245247c7ee0e8bb8d..25b77e8f59afab40b132765ad21938df69a0dc60 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -38,6 +38,7 @@ set(SRC
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp
Cell.cpp CellComplex.cpp ChainComplex.cpp Homology.cpp
+ Curvature.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 35f66c4fd7e23b6f76eea9ae30fa82313a3fb423..b78a6528abd8d5b802bf47fcc308c4d3ae44b492 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -13,7 +13,6 @@
#define NEXT(i) ((i)<2 ? (i)+1 : (i)-2)
#define PREV(i) ((i)>0 ? (i)-1 : (i)+2)
-//std::cout << "NEXT(1) = " << NEXT(1) << std::endl;
//========================================================================================================
//CONSTRUCTOR
@@ -32,32 +31,32 @@ Curvature::~Curvature()
void Curvature::retrievePhysicalSurfaces(const std::string & face_tag)
{
-// GFaceList ptFaceList; //Pointer to face
- GEntityList ptTempEntityList;
+ GEntityList ptTempEntityList;
+ //This is a vector of 4 maps: 1 for 0D, one for 1D, one for 2D, one for 3D
+ std::map<int, GEntityList > physicals[4];
+ _model->getPhysicalGroups(physicals);
- std::map<int, GEntityList > physicals[4]; //The right vector is a vector of 4 maps: 1 for 0D, one for 1D, one for 2D, one for 3D
- _model->getPhysicalGroups(physicals);
- std::map<int, GEntityList > surface_physicals = physicals[2];// Here we need only the 2nd maps which represents the faces (triangles)
+ // Here we need only the 2nd map which represents the faces (triangles)
- for (GEntityIter it = surface_physicals.begin(); it != surface_physicals.end(); it++) //Loop over physical names
+ std::map<int, GEntityList > surface_physicals = physicals[2];
+ for (GEntityIter it = surface_physicals.begin(); it != surface_physicals.end(); it++) //Loop over physical names
+ {
+ std::string physicalTag = _model->getPhysicalName(2, it->first);
+
+ if (physicalTag == face_tag) //face_tag is one of the argument of "compute_curvature"
{
- std::string physicalTag = _model->getPhysicalName(2, it->first);
-// std::cout << "The physical tag we have stored is: "<< physicalTag << std::endl;
+ ptTempEntityList.clear(); // the vector is cleared before the for-loop where is it filled
+ ptTempEntityList = it->second;
- if (physicalTag == face_tag) //face_tag is one of the argument of "compute_curvature"
+ for (unsigned int iEnt = 0; iEnt < ptTempEntityList.size(); iEnt++)
{
- ptTempEntityList.clear(); // the vector is cleared before the for-loop where is it filled
- ptTempEntityList = it->second;
-
- for (unsigned int iEnt = 0; iEnt < ptTempEntityList.size(); iEnt++)
- {
- GEntity* entity = ptTempEntityList[iEnt];
- _ptFinalEntityList.push_back(entity);
+ GEntity* entity = ptTempEntityList[iEnt];
+ _ptFinalEntityList.push_back(entity);
- }
}
- }
+ }
+ }
// //To check that all the surfaces with the physical-tag "face_tag", are stored in ptFinalEntityList:
// std::cout << "Stored physical tags:" << std::endl;
// for(unsigned int i =0; i < _ptFinalEntityList.size(); ++i)
@@ -174,18 +173,15 @@ void Curvature::computeVertexNormals()
_VertexNormal[V1] += cross;
_VertexNormal[V2] += cross;
-
} // end of loop over elements of one entity
} //Loop over _ptFinalEntityList
-
///////Normalize the vertex-normals.
for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
{
_VertexNormal[n].normalize();
}
-
}
//========================================================================================================
@@ -200,7 +196,6 @@ void Curvature::curvatureTensor()
_CurveTensor.resize(_VertexToInt.size());
-
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
{
GEntity* entity = _ptFinalEntityList[i]; //entity is a pointer to one surface of the group "FinalEntityList"
@@ -219,27 +214,26 @@ void Curvature::curvatureTensor()
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]);
const double Wij1 = _TriangleArea[E] / (2 * _VertexArea[V1]);
//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); //Dot-product of the 2 vectors
- const double k_nominator1 = -dot(_VertexNormal[V1], vector_AB); //Dot-product of the 2 vectors
+ 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
+ const double coeff = 2.0/vector_AB.normSq(); //normSq is the norm to the power 2
const double KijAB_0 = coeff*k_nominator0;
const double KijAB_1 = coeff*k_nominator1;
//Projection of the edge vector AB on the tangential plane:
//** For Vertex 0
tensprod(_VertexNormal[V0], _VertexNormal[V0], TempTensor);
- TempTensor *= -1.0; //-N*Nt
- TempTensor(0,0) += 1.0; //I-N*Nt
- TempTensor(1,1) += 1.0;
- TempTensor(2,2) += 1.0;
+ TempTensor *= -1.0; //-N*Nt
+ TempTensor(0,0) += 1.0; //I-N*Nt
+ TempTensor(1,1) += 1.0;
+ TempTensor(2,2) += 1.0;
for (int m = 0; m<3; ++m)
{
@@ -249,17 +243,17 @@ void Curvature::curvatureTensor()
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
+
TijAB.normalize();
tensprod(TijAB, TijAB, TijABTensorProduct);
_CurveTensor[V0] += Wij0*KijAB_0*TijABTensorProduct;
-
//** For Vertex 1
tensprod(_VertexNormal[V1], _VertexNormal[V1], TempTensor);
- TempTensor *= -1.0; //-N*Nt
- TempTensor(0,0) += 1.0; //I-N*Nt
- TempTensor(1,1) += 1.0;
- TempTensor(2,2) += 1.0;
+ TempTensor *= -1.0; //-N*Nt
+ TempTensor(0,0) += 1.0; //I-N*Nt
+ TempTensor(1,1) += 1.0;
+ TempTensor(2,2) += 1.0;
for (int m = 0; m<3; ++m)
{
@@ -269,11 +263,11 @@ void Curvature::curvatureTensor()
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
+
TijAB.normalize();
tensprod(TijAB, TijAB, TijABTensorProduct);
_CurveTensor[V1] += Wij1*KijAB_1*TijABTensorProduct;
-
}//end of loop over the edges
}//end of loop over all elements of one GEntity
@@ -286,78 +280,143 @@ void Curvature::curvatureTensor()
void Curvature::computeCurvature_Simple()
{
- SVector3 vector_E;
- SVector3 vector_A;
- SVector3 vector_B;
- SVector3 vector_Wvi;
- STensor3 Qvi;
- STensor3 QviT;
- STensor3 Holder;
+ SVector3 vector_E;
+ SVector3 vector_A;
+ SVector3 vector_B;
+ SVector3 vector_Wvi;
+ STensor3 Qvi;
+ STensor3 QviT;
+ STensor3 Holder;
- initializeMap();
- computeVertexNormals();
- curvatureTensor();
+ initializeMap();
+ computeVertexNormals();
+ curvatureTensor();
+
+ _VertexCurve.resize(_VertexToInt.size());
- _VertexCurve.resize(_VertexToInt.size());
+ 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];
+ const double MagA = vector_A.norm();
+ const double MagB = vector_B.norm();
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n) //Loop over the vertex
+ if (MagB > MagA)
+ {
+ vector_Wvi = vector_B;
+ }
+ else
{
- vector_E = SVector3(1,0,0);
- vector_A = vector_E + _VertexNormal[n];
- vector_B = vector_E - _VertexNormal[n];
+ vector_Wvi = vector_A;
+ }
- const double MagA = vector_A.norm();
- const double MagB = vector_B.norm();
+ vector_Wvi.normalize();
+ //to obtain the Qvi = Id -2*Wvi*Wvi^t
+ tensprod(vector_Wvi, vector_Wvi, Qvi); //Qvi = Wvi*Wvi^t
+ Qvi *= -2.0; //-2*Wvi*Wvi^t
+ Qvi(0,0) += 1.0; //I - 2*Wvi*Wvi^t ==> Householder transformation
+ Qvi(1,1) += 1.0;
+ Qvi(2,2) += 1.0;
- if (MagB > MagA)
- {
- vector_Wvi = vector_B;
- }
- else
- {
- vector_Wvi = vector_A;
+ //Transpose the matrix:
+ for (int i = 0; i<3; ++i)
+ {
+ for (int j = 0; j<3; ++j)
+ {
+ QviT(i,j) = Qvi(j,i);
}
+ }
- vector_Wvi.normalize();
- //to obtain the Qvi = Id -2*Wvi*Wvi^t
- tensprod(vector_Wvi, vector_Wvi, Qvi); //Qvi = Wvi*Wvi^t
- Qvi *= -2.0; //-2*Wvi*Wvi^t
- Qvi(0,0) += 1.0; //I - 2*Wvi*Wvi^t ==> Householder transformation
- Qvi(1,1) += 1.0;
- Qvi(2,2) += 1.0;
-
- //To create the transpose matrix:
- 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;
- QviT *= _CurveTensor[n];
- QviT *=Qvi;
- Holder = QviT;
+ //Eigenvalues of the 2*2 minor from the Householder matrix
+ const double A = 1.0;
+ const double B = -(Holder(1,1) + Holder(2,2));
+ 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);
- //Eigenvalues of the 2*2 minor from the Householder matrix
- const double A = 1.0;
- const double B = -(Holder(1,1) + Holder(2,2));
- 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)
+ {
+ std::cout << "WARNING: negative discriminant: " << B*B-4.*A*C << std::endl;
+ }
- if((B*B-4.*A*C) < 0.0)
- {
- std::cout << "WARNING: negative discriminant: " << B*B-4.*A*C << std::endl;
- }
+ const double m11 = (-B + Delta)/(2*A); //Eigenvalue of Householder submatrix
+ const double m22 = (-B - Delta)/(2*A);
+
+ //_VertexCurve[n] = (3*m11-m22)*(3*m22-m11); //Gaussian Curvature
+ _VertexCurve[n] = ((3*m11-m22) + (3*m22-m11))*0.5; //Mean Curvature
+
+ }
+}
+
+//========================================================================================================
+//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
- const double m11 = (-B + Delta)/(2*A); //Eigenvalue of Householder submatrix
- const double m22 = (-B - Delta)/(2*A);
+void Curvature::computeRusinkiewiczNormals()
+{
+ SVector3 vector_AB;
+ SVector3 vector_AC;
+ SVector3 vector_BC;
- //_VertexCurve[n] = (3*m11-m22)*(3*m22-m11); //Gaussian Curvature
- _VertexCurve[n] = ((3*m11-m22) + (3*m22-m11))*0.5; //Mean Curvature
+ _TriangleArea.resize(_ElementToInt.size() );
+ _VertexNormal.resize(_VertexToInt.size());
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ {
+ // entity is a pointer to one surface of the group "FinalEntityList"
+ GEntity* entity = _ptFinalEntityList[i];
+
+ //Loop over the element all the element of the "myTag"-surface
+ for (int iElem = 0; iElem < entity->getNumMeshElements(); iElem++)
+ {
+ // Pointer to one element
+ MElement *e = entity->getMeshElement(iElem);
+ // The NEW tag of the corresponding element
+ const int E = _ElementToInt[e->getNum()];
+ // std::cout << "We are now looking at element Nr: " << E << std::endl;
+
+ // Pointers to vertices of triangle
+ MVertex* A = e->getVertex(0);
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ const int V0 = _VertexToInt[A->getNum()]; //The new number of the vertex
+ const int V1 = _VertexToInt[B->getNum()];
+ const int V2 = _VertexToInt[C->getNum()];
+
+ vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
+ vector_AC = SVector3(C->x() - A->x(), C->y() - A->y(), C->z() - A->z() );
+ vector_BC = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z() );
+
+
+ const SVector3 cross = crossprod(vector_AB, vector_AC);
+
+ // 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;
+
+
+ const double l_AB = vector_AB.normSq();
+ const double l_AC = vector_AC.normSq();
+ const double l_BC = vector_BC.normSq();
+
+ _VertexNormal[V0] += cross * (1.0/ (l_AB*l_AC));
+ _VertexNormal[V1] += cross * (1.0/ (l_AB*l_BC));
+ _VertexNormal[V2] += cross * (1.0/ (l_AC*l_BC));
+
+ } // end of loop over elements of one entity
+
+ } //Loop over _ptFinalEntityList
+
+ ///////Normalize the vertex-normals.
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
+ {
+ _VertexNormal[n].normalize();
}
}
@@ -367,96 +426,90 @@ void Curvature::computeCurvature_Simple()
void Curvature::computePointareas()
{
-// if (_pointareas.size() == _VertexToInt.size())
-// return;
-// need_faces();
-
- std::cout << "Computing point areas... " << std::endl;
+ std::cout << "Computing point areas... " << std::endl;
// std::cout << "The mesh has " << _VertexToInt.size() << " nodes" << std::endl;
- SVector3 e[3];
- SVector3 l2;
- SVector3 ew;
+ SVector3 e[3];
+ SVector3 l2;
+ SVector3 ew;
- _pointareas.resize(_VertexToInt.size());
- _cornerareas.resize(_ElementToInt.size());
+ _pointareas.resize(_VertexToInt.size());
+ _cornerareas.resize(_ElementToInt.size());
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ {
+ GEntity* entity = _ptFinalEntityList[i]; //entity is a pointer to one surface of the group "FinalEntityList"
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ for (int iElem = 0; iElem < entity->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
{
- GEntity* entity = _ptFinalEntityList[i]; //entity is a pointer to one surface of the group "FinalEntityList"
+ MElement *E = entity->getMeshElement(iElem); //Pointer to one element
+ // The NEW tag of the corresponding element
+ const int EIdx = _ElementToInt[E->getNum()];
- for (int iElem = 0; iElem < entity->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
- MElement *E = entity->getMeshElement(iElem); //Pointer to one element
- // The NEW tag of the corresponding element
- const int EIdx = _ElementToInt[E->getNum()];
-
- MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
- MVertex* B = E->getVertex(1);
- MVertex* C = E->getVertex(2);
-
- //Edges
- e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangilar 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());
-
- // Compute corner weights
- //SVector3 len = crossprod(e[0], e[1]);
- //len = norm
- //len2 = normSq
- double area = 0.5 * norm(crossprod(e[0], e[1])); //area of a triangle
- l2 = SVector3( normSq(e[0]), normSq(e[1]), normSq(e[2]) );
- ew = SVector3( l2[0] * (l2[1] + l2[2] - l2[0]),
- l2[1] * (l2[2] + l2[0] - l2[1]),
- l2[2] * (l2[0] + l2[1] - l2[2]) );
-
- 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];
+ MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = E->getVertex(1);
+ MVertex* C = E->getVertex(2);
- }
- 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)
- {
- _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
- {
- 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]);
- }
+ //Edges
+ e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangilar 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());
+
+ // Compute corner weights
+ //SVector3 len = crossprod(e[0], e[1]);
+ //len = norm
+ //len2 = normSq
+ double area = 0.5 * norm(crossprod(e[0], e[1])); //area of a triangle
+ l2 = SVector3( normSq(e[0]), normSq(e[1]), normSq(e[2]) );
+ ew = SVector3( l2[0] * (l2[1] + l2[2] - l2[0]),
+ l2[1] * (l2[2] + l2[0] - l2[1]),
+ l2[2] * (l2[0] + l2[1] - l2[2]) );
+
+ 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)
+ {
+ _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)
+ {
+ _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
+ {
+ 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()];
+ 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[V0] += _cornerareas[EIdx][0];
+ //_pointareas[faces[i][0]] += _cornerareas[i][0];
- _pointareas[V1] += _cornerareas[EIdx][1];
+ _pointareas[V1] += _cornerareas[EIdx][1];
- _pointareas[V2] += _cornerareas[EIdx][2];
+ _pointareas[V2] += _cornerareas[EIdx][2];
- }//End of loop over iElem
+ } //End of loop over iElem
- std::cout << "Done computing pointareas" << std::endl;
+ std::cout << "Done computing pointareas" << 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 loop over _ptFinalEntityList
-}//End of the method "computePointareas"
+} //End of the method "computePointareas"
//========================================================================================================
@@ -475,16 +528,14 @@ void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, cons
new_u = -1.0*new_u;
new_v = -1.0*new_v;
return;
- }//end of if-condtion
+ }
SVector3 perp_old = new_norm - ndot*old_norm;
SVector3 dperp = 1.0f/(1 + ndot) * (old_norm + new_norm);
new_u -= dperp*dot(new_u, perp_old);
new_v -= dperp*dot(new_v, perp_old);
-
}
-
//========================================================================================================
//Project a curvature tensor from the basis spanned by old_u and old_v
@@ -500,16 +551,14 @@ void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
SVector3 r_new_v;
rot_coord_sys(new_u, new_v, crossprod(old_u,old_v), r_new_u, r_new_v);
- double u1 = dot(r_new_u, old_u);
- double v1 = dot(r_new_u, old_v);
- double u2 = dot(r_new_v, old_u);
- double v2 = dot(r_new_v, old_v);
+ const double u1 = dot(r_new_u, old_u);
+ const double v1 = dot(r_new_u, old_v);
+ const double u2 = dot(r_new_v, old_u);
+ const double v2 = dot(r_new_v, old_v);
new_ku = old_ku*u1*u1 + old_kuv*(2.0f * u1*v1) + old_kv*v1*v1;
new_kuv = old_ku*u1*u2 + old_kuv*(u1*v2 * u2*v1) + old_kv*v1*v2;
new_kv = old_ku*u2*u2 + old_kuv*(2.0f * u2*v2) + old_kv*v2*v2;
-
-
}
@@ -564,7 +613,7 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
void Curvature::computeCurvature_Rusinkiewicz()
{
initializeMap();
- computeVertexNormals();
+ computeRusinkiewiczNormals();
computePointareas();
SVector3 e[3];
@@ -593,8 +642,6 @@ void Curvature::computeCurvature_Rusinkiewicz()
_curv2.resize(_VertexToInt.size());
_curv12.resize(_VertexToInt.size());
-
-
std::cout << "Computing Curvature.." << std::endl;
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
@@ -613,138 +660,167 @@ void Curvature::computeCurvature_Rusinkiewicz()
const int V1 = _VertexToInt[B->getNum()];
const int V2 = _VertexToInt[C->getNum()];
-
///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 (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
- {
- _pdir1[ivertex] = crossprod(_pdir1[ivertex], _VertexNormal[ivertex]);
- _pdir1[ivertex].normalize();
- _pdir2[ivertex] = crossprod(_VertexNormal[ivertex], _pdir1[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)
+ // Compute curvature per face:
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ {
+ //entity is a pointer to one surface of the group "FinalEntityList"
+ GEntity* entity = _ptFinalEntityList[i];
+
+ //Loop over all elements of this entity:
+ for (int iElem = 0; iElem < entity->getNumMeshElements(); iElem++)
{
- GEntity* entity = _ptFinalEntityList[i]; //entity is a pointer to one surface of the group "FinalEntityList"
- for (int iElem = 0; iElem < entity->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ MElement *E = entity->getMeshElement(iElem); //Pointer to one element
+ // The NEW tag of the corresponding element
+ const int EIdx = _ElementToInt[E->getNum()];
+
+ MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = E->getVertex(1);
+ MVertex* C = E->getVertex(2);
+
+ e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
+ 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()) };
+
+ //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:
+ 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)
{
- MElement *E = entity->getMeshElement(iElem); //Pointer to one element
- // The NEW tag of the corresponding element
- const int EIdx = _ElementToInt[E->getNum()];
-
- MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
- MVertex* B = E->getVertex(1);
- MVertex* C = E->getVertex(2);
-
- e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
- 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()) };
-
-
- //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:
- 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 j = 0; j<3; ++j)
- {
- w(i,j) = 0.0;
- }
+ w(i,j) = 0.0;
}
+ }
- //filling:
- for (int j = 0; j< 3; ++j)
- {
- u = dot(e[j], t);
- v = dot(e[j], b);
+ //filling:
+ for (int j = 0; j< 3; ++j)
+ {
+ u = dot(e[j], t);
+ v = dot(e[j], b);
- w(0,0) += u*u;
- w(0,1) += u*v;
- w(2,2) += v*v;
+ w(0,0) += u*u;
+ w(0,1) += u*v;
+ w(2,2) += v*v;
- MVertex* U = E->getVertex(PREV(j));
- MVertex* V = E->getVertex(NEXT(j));
+ MVertex* U = E->getVertex(PREV(j));
+ MVertex* V = E->getVertex(NEXT(j));
- const int UIdx = _VertexToInt[U->getNum()];
- const int VIdx = _VertexToInt[V->getNum()];
+ const int UIdx = _VertexToInt[U->getNum()];
+ const int VIdx = _VertexToInt[V->getNum()];
- dn = _VertexNormal[UIdx] - _VertexNormal[VIdx];
+ dn = _VertexNormal[UIdx] - _VertexNormal[VIdx];
- dnu = dot(dn, t);
- dnv = dot(dn, b);
+ dnu = dot(dn, t);
+ dnv = dot(dn, b);
- m[0] += dnu*u;
- m[1] += dnu*v + dnv*u;
- m[2] += dnv*v;
- }
+ m[0] += dnu*u;
+ m[1] += dnu*v + dnv*u;
+ m[2] += dnv*v;
+ }
- w(1,1) = w(0,0) + w(2,2);
- w(1,2) = w(0,1);
+ w(1,1) = w(0,0) + w(2,2);
+ w(1,2) = w(0,1);
- //Least Square Solution
- double diag[3];
- if (!ldltdc(w, diag))
- {
- std::cout << "ldltdc failed" << std::endl;
- continue;
- }
- ldltsl(w, diag, m, m);
+ //Least Squares Solution
+ double diag[3];
+ 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)
- {
- 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];
+ //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;
- _curv1[vj] += wt*c1;
- _curv12[vj] += wt*c12;
- _curv2[vj] += wt*c2;
-
- }
+ _curv1[vj] += wt*c1;
+ _curv12[vj] += wt*c12;
+ _curv2[vj] += wt*c2;
+ }
} //End of loop over the element (iElem)
} //End of loop over "_ptFinalEntityList"
- //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]);
- }
+ //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]);
+ }
- std::cout << "Done" << std::endl;
+ std::cout << "Done" << std::endl;
- _VertexCurve.resize( _VertexToInt.size() );
+ _VertexCurve.resize( _VertexToInt.size() );
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
- {
- _VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5;
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
+ {
+ //**Mean Curvature:**
+ _VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5;
+ //_VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]);
+ //**Gaussian Curvature:**
+ //_VertexCurve[ivertex] = (_curv1[ivertex]*_curv2[ivertex]);
+ //_VertexCurve[ivertex] = std::abs(_VertexCurve[ivertex]);
+
+ }
+
+/*
+ std::vector<double>::iterator MaxVal = std::max_element(_VertexCurve.begin(), _VertexCurve.end());
+
+ const double max_double_value = *MaxVal;
+ std::cout << "The max value is: " << max_double_value << std::endl;
+
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
+ {
+ _VertexCurve[ivertex]= 1.1*max_double_value - _VertexCurve[ivertex];
+// _VertexCurve[ivertex]= MaxVal - _VertexCurve[ivertex];
+
+ }
+*/
+
+ //Emilie Marchandise:
+
+ const int N = 5;
+
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
+ {
+ //_VertexCurve[ivertex] = 2.*3.14159/( std::abs(_VertexCurve[ivertex]) * N);
+ _VertexCurve[ivertex] = 1./ std::pow( std::exp(std::abs(_VertexCurve[ivertex])), 0.25 );
+ }
- }
} //End of the "computeCurvature_Rusinkiewicz" method
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index a11dbfca6501e4c18e7cb946e0ec9fd1d866d710..c856b49ae0de5ee63bf70422364beb964c578809 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -82,6 +82,7 @@ private:
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.
diff --git a/benchmarks/curvature/ComputeCurvature.py b/benchmarks/curvature/ComputeCurvature.py
index 9c3bc2193836e003ad17535c9d9291c745bdd1d0..53704e844ea168bfc2ff1370e2105606d3c5f9bc 100644
--- a/benchmarks/curvature/ComputeCurvature.py
+++ b/benchmarks/curvature/ComputeCurvature.py
@@ -26,7 +26,7 @@ print "Model is loaded"
cv = Curvature(g)
cv.retrievePhysicalSurfaces("Wall")
-cv.computeCurvature()
+cv.computeCurvature_Rusinkiewicz()
cv.writeToFile("result.pos")