Skip to content
Snippets Groups Projects
Commit 8b77e1ec authored by Emilie Sauvage's avatar Emilie Sauvage
Browse files

Update on curvature algorithm

parent 885d2954
No related branches found
No related tags found
No related merge requests found
......@@ -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
......@@ -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)
......
......@@ -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,18 +31,18 @@ Curvature::~Curvature()
void Curvature::retrievePhysicalSurfaces(const std::string & face_tag)
{
// GFaceList ptFaceList; //Pointer to face
GEntityList ptTempEntityList;
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
//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 > 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)
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);
// std::cout << "The physical tag we have stored is: "<< physicalTag << std::endl;
if (physicalTag == face_tag) //face_tag is one of the argument of "compute_curvature"
{
......@@ -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,7 +214,6 @@ 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]);
......@@ -227,7 +221,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); //Dot-product of the 2 vectors
const double k_nominator1 = -dot(_VertexNormal[V1], 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 KijAB_0 = coeff*k_nominator0;
......@@ -249,11 +243,11 @@ 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
......@@ -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
......@@ -300,7 +294,6 @@ void Curvature::computeCurvature_Simple()
_VertexCurve.resize(_VertexToInt.size());
for (unsigned int n = 0; n < _VertexToInt.size(); ++ n) //Loop over the vertex
{
vector_E = SVector3(1,0,0);
......@@ -327,7 +320,7 @@ void Curvature::computeCurvature_Simple()
Qvi(1,1) += 1.0;
Qvi(2,2) += 1.0;
//To create the transpose matrix:
//Transpose the matrix:
for (int i = 0; i<3; ++i)
{
for (int j = 0; j<3; ++j)
......@@ -351,7 +344,6 @@ void Curvature::computeCurvature_Simple()
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);
......@@ -363,14 +355,77 @@ void Curvature::computeCurvature_Simple()
//========================================================================================================
//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
void Curvature::computeRusinkiewiczNormals()
{
SVector3 vector_AB;
SVector3 vector_AC;
SVector3 vector_BC;
_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();
}
}
//========================================================================================================
// Compute per-vertex point areas
void Curvature::computePointareas()
{
// if (_pointareas.size() == _VertexToInt.size())
// return;
// need_faces();
std::cout << "Computing point areas... " << std::endl;
// std::cout << "The mesh has " << _VertexToInt.size() << " nodes" << std::endl;
......@@ -381,7 +436,6 @@ void Curvature::computePointareas()
_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"
......@@ -416,7 +470,6 @@ void Curvature::computePointareas()
_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)
{
......@@ -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,13 +660,11 @@ 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());
}
}
......@@ -633,8 +678,11 @@ void Curvature::computeCurvature_Rusinkiewicz()
// Compute curvature per face:
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 iElem = 0; iElem < entity->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
//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++)
{
MElement *E = entity->getMeshElement(iElem); //Pointer to one element
// The NEW tag of the corresponding element
......@@ -650,7 +698,6 @@ void Curvature::computeCurvature_Rusinkiewicz()
//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();
......@@ -700,7 +747,7 @@ void Curvature::computeCurvature_Rusinkiewicz()
w(1,1) = w(0,0) + w(2,2);
w(1,2) = w(0,1);
//Least Square Solution
//Least Squares Solution
double diag[3];
if (!ldltdc(w, diag))
{
......@@ -721,7 +768,6 @@ void Curvature::computeCurvature_Rusinkiewicz()
_curv1[vj] += wt*c1;
_curv12[vj] += wt*c12;
_curv2[vj] += wt*c2;
}
} //End of loop over the element (iElem)
......@@ -742,9 +788,39 @@ void Curvature::computeCurvature_Rusinkiewicz()
for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
{
//**Mean Curvature:**
_VertexCurve[ivertex] = (_curv1[ivertex]+_curv2[ivertex])*0.5;
//_VertexCurve[ivertex] = std::abs(_curv1[ivertex]) + std::abs(_curv2[ivertex]);
//**Gaussian Curvature:**
//_VertexCurve[ivertex] = (_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
......
......@@ -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.
......
......@@ -26,7 +26,7 @@ print "Model is loaded"
cv = Curvature(g)
cv.retrievePhysicalSurfaces("Wall")
cv.computeCurvature()
cv.computeCurvature_Rusinkiewicz()
cv.writeToFile("result.pos")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment