Testing a curvature algorithm

Geo/Curvature.cpp

+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include "Curvature.h"
+#include "MElement.h"
+
+#include<iostream>
+#include<fstream>
+#include<cmath>
+
+//========================================================================================================
+
+//CONSTRUCTOR
+Curvature::Curvature(GModel* model)
+{
+    _model = model;
+}
+//========================================================================================================
+
+//DESTRUCTOR
+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
+    _model->getPhysicalGroups(physicals);
+    std::map<int, GEntityList > surface_physicals = physicals[2];// Here we need only the 2nd maps which represents the faces (triangles)
+
+    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"
+      {
+        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);
+
+        }
+      }
+   }
+//   //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)
+//   {
+//    std::cout << _ptFinalEntityList[i]->tag() << std::endl; //some method that prints the name of the entity
+//   }
+
+}
+
+//========================================================================================================
+
+//INITIALIZATION OF THE MAP  AND  RENUMBERING OF THE SELECTED ENTITIES:
+
+void Curvature::initializeMap()
+{
+
+  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);
+      // Tag of the corresponding element
+      const int E = e->getNum();
+      // std::cout << "We are now looking at element Nr: " << E << std::endl;
+      _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;
+
+  // 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;
+  }
+
+}
+
+//========================================================================================================
+
+void Curvature::computeVertexNormals()
+{
+  SVector3 vector_AB;
+  SVector3 vector_AC;
+
+  _TriangleArea.resize(_ElementToInt.size() );
+  _VertexArea.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() );
+
+      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;
+
+      _VertexArea[V0] += _TriangleArea[E];
+      _VertexArea[V1] += _TriangleArea[E];
+      _VertexArea[V2] += _TriangleArea[E];
+
+      _VertexNormal[V0] += cross;  //here we are actually computing the unit normal vector per vertex
+      _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();
+    }
+
+}
+
+//========================================================================================================
+
+void Curvature::curvatureTensor()
+{
+
+  STensor3 TempTensor;
+  STensor3 TijABTensorProduct;
+  SVector3 TijAB;
+  SVector3 vector_AB;
+
+  _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"
+
+    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
+      const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
+
+      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]);
+        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 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;
+
+        for (int m = 0; m<3; ++m)
+        {
+          TijAB(m) = 0.0;
+          for (int n = 0; n<3; ++n)
+          {
+            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;
+
+         for (int m = 0; m<3; ++m)
+        {
+          TijAB(m) = 0.0;
+          for (int n = 0; n<3; ++n)
+          {
+            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
+
+  }//End of loop over ptFinalEntityList
+
+}//End of method
+
+//========================================================================================================
+
+void Curvature::computeCurvature()
+{
+    SVector3 vector_E;
+    SVector3 vector_A;
+    SVector3 vector_B;
+    SVector3 vector_Wvi;
+    STensor3 Qvi;
+    STensor3 QviT;
+    STensor3 Holder;
+
+    initializeMap();
+    computeVertexNormals();
+    curvatureTensor();
+
+    _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();
+
+       if (MagB > MagA)
+       {
+          vector_Wvi = vector_B;
+       }
+       else
+       {
+          vector_Wvi = vector_A;
+       }
+
+       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;
+
+       //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;
+       }
+
+
+       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
+
+    }
+}
+
+//========================================================================================================
+
+void Curvature::writeToFile( const std::string & filename)
+{
+  std::ofstream outfile;
+  outfile.precision(18);
+  outfile.open(filename.c_str());
+  outfile << "View \"Curvature \"{" << std::endl;
+
+  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
+    {
+      MElement *e = entity->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 << "ST("; //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 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 << _VertexCurve[V1] << "," << _VertexCurve[V2] << "," << _VertexCurve[V3];
+
+      outfile << "};" << std::endl;
+
+  } //Loop over elements
+
+} // Loop over ptFinalEntityList
+
+outfile << "};" << std::endl;
+
+outfile.close();
+}
+
+//========================================================================================================
+
+

Geo/Curvature.h

+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _CURVATUREL_H_
+#define _CURVATURE_H_
+
+#include "GModel.h"
+#include "STensor3.h"
+
+#include<map>
+#include<vector>
+
+class Curvature
+{
+private:
+
+    //-----------------------------------------
+    // TYPEDEFS
+
+    typedef std::vector<GEntity*> GEntityList;
+    typedef std::vector<GFace*> GFaceList;
+    typedef std::map<int, GEntityList >::iterator GEntityIter;
+    typedef std::map<int, GFaceList >::iterator GFaceIter;
+
+    //-----------------------------------------
+    // MEMBER VARIABLES
+
+    //Map definition
+    std::map<int, int> _VertexToInt;
+    std::map<int, int> _ElementToInt;
+
+    //Model and list of selected entities with give physical tag:
+    GModel* _model;    
+
+    GEntityList _ptFinalEntityList;
+
+    //Averaged vertex normals
+    std::vector<SVector3> _VertexNormal;
+
+    //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;
+
+
+    //-----------------------------------------
+    // PRIVATE METHODS
+
+    void initializeMap();
+    void computeVertexNormals();
+    void curvatureTensor();
+
+
+
+public:
+
+  //-----------------------------------------
+  // PUBLIC METHODS
+  //Constructor
+  Curvature(GModel* model);
+
+  //Destructor
+  ~Curvature();
+
+  void retrievePhysicalSurfaces(const std::string & face_tag);
+  void computeCurvature();
+  void writeToFile( const std::string & filename);
+
+
+
+};
+
+
+#endif

benchmarks/curvature/Aorta.geo

+Merge "Aorta_Source.msh";
+
+Physical Surface("Wall") = {1};

benchmarks/curvature/Bifurcation.geo

+Merge "Bifurcation_Source.msh";
+
+Physical Surface("Wall") = {1};

benchmarks/curvature/ComputeCurvature.py

+from dgpy import *
+from gmshpy import *
+import os
+import sys
+
+#if (os.path.isfile("Torus.msh")) :
+  #print('---- The mesh file already exists ----');
+
+#else :
+  #cmd = "gmsh -3 Torus.geo"
+  #os.system(cmd)
+  #print('---- Mesh has been executed ----');
+  
+#meshName = "Torus"
+#meshName = "Sphere"
+meshName = "Bifurcation"
+
+g = GModel()
+
+#g.load(meshName+".geo")
+g.load(meshName+".msh")
+print "Model is loaded"
+
+#g.compute_curvature("Wall","curvature.pos")
+
+
+cv = Curvature(g)
+cv.retrievePhysicalSurfaces("Wall")
+cv.computeCurvature()
+cv.writeToFile("result.pos")
+
+
+  
+print "Youpiii Exit :-)"
+sys.exit(0)

benchmarks/curvature/Sphere.geo

+// Gmsh project created on Fri Feb  1 12:17:44 2008
+
+//Delete All;
+
+RSphere = 1.;
+lcSphere = .25;
+
+RDom = 10;
+lcDom = 0.7; //1.0
+
+/*
+Point(1) = {0,0,0,lcSphere};
+Point(2) = {RSphere,0,0,lcSphere};
+Point(3) = {0,RSphere,0,lcSphere};
+Circle(1) = {2,1,3};
+Point(4) = {-RSphere,0,0.0,lcSphere};
+Point(5) = {0,-RSphere,0.0,lcSphere};
+Circle(2) = {3,1,4};
+Circle(3) = {4,1,5};
+Circle(4) = {5,1,2};
+Point(6) = {0,0,-RSphere,lcSphere};
+Point(7) = {0,0,RSphere,lcSphere};
+Circle(5) = {3,1,6};
+Circle(6) = {6,1,5};
+Circle(7) = {5,1,7};
+Circle(8) = {7,1,3};
+Circle(9) = {2,1,7};
+Circle(10) = {7,1,4};
+Circle(11) = {4,1,6};
+Circle(12) = {6,1,2};
+Line Loop(13) = {2,8,-10};
+Ruled Surface(14) = {13};
+Line Loop(15) = {10,3,7};
+Ruled Surface(16) = {15};
+Line Loop(17) = {-8,-9,1};
+Ruled Surface(18) = {17};
+Line Loop(19) = {-11,-2,5};
+Ruled Surface(20) = {19};
+Line Loop(21) = {-5,-12,-1};
+Ruled Surface(22) = {21};
+Line Loop(23) = {-3,11,6};
+Ruled Surface(24) = {23};
+Line Loop(25) = {-7,4,9};
+Ruled Surface(26) = {25};
+Line Loop(27) = {-4,12,-6};
+Ruled Surface(28) = {27};
+Surface Loop(29) = {28,26,16,14,20,24,22,18};
+*/
+
+Point(101) = {0,0,0,lcDom};
+Point(102) = {RDom,0,0,lcDom};
+Point(103) = {0,RDom,0,lcDom};
+Circle(101) = {102,101,103};
+Point(104) = {-RDom,0,0.0,lcDom};
+Point(105) = {0,-RDom,0.0,lcDom};
+Circle(102) = {103,101,104};
+Circle(103) = {104,101,105};
+Circle(104) = {105,101,102};
+Point(106) = {0,0,-RDom,lcDom};
+Point(107) = {0,0,RDom,lcDom};
+Circle(105) = {103,101,106};
+Circle(106) = {106,101,105};
+Circle(107) = {105,101,107};
+Circle(108) = {107,101,103};
+Circle(109) = {102,101,107};
+Circle(110) = {107,101,104};
+Circle(111) = {104,101,106};
+Circle(112) = {106,101,102};
+Line Loop(113) = {102,108,-110};
+Ruled Surface(114) = {113};
+Line Loop(115) = {110,103,107};
+Ruled Surface(116) = {115};
+Line Loop(117) = {-108,-109,101};
+Ruled Surface(118) = {117};
+Line Loop(119) = {-111,-102,105};
+Ruled Surface(120) = {119};
+Line Loop(121) = {-105,-112,-101};
+Ruled Surface(122) = {121};
+Line Loop(123) = {-103,111,106};
+Ruled Surface(124) = {123};
+Line Loop(125) = {-107,104,109};
+Ruled Surface(126) = {125};
+Line Loop(127) = {-104,112,-106};
+Ruled Surface(128) = {127};
+Surface Loop(129) = {128,126,116,114,120,124,122,118};
+
+//Volume(200) = {129,29};
+Volume(200) = {129};
+
+
+// Wall ID = 2100
+//Physical Surface(2100) = {28,26,16,14,20,24,22,18};
+// Farfield CBC ID = 2222
+Physical Surface("Wall") = {128,126,116,114,120,124,122,118};
+// Interior ID = 1000
+Physical Volume(1000) = {200};
+
+

benchmarks/curvature/Torus.geo

+lc = 0.1;          
+Point(1) = {2.0,0.0,0.0,lc};          
+Point(2) = {2.0,1,0.0,lc};          
+Point(3) = {1,0,0.0,lc};          
+Point(4) = {3,0,0.0,lc};          
+Point(5) = {2,-1,0.0,lc}; 
+         
+Circle(1) = {4,1,2};                 
+Circle(2) = {2,1,3};         
+Circle(3) = {3,1,5};         
+Circle(4) = {5,1,4};         
+Line Loop(5) = {4,1,2,3};   //Inlet      
+Plane Surface(6) = {5};         
+       
+//out[] = Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159};
+out[] = Extrude{ {0.0, 1.0, 0.0},{ 0.0, 0.0, 0.0 }, Pi } { Surface{6}; };         
+//Recombine Surface {6, 27, 23, 15, 19, 28};
+
+Physical Surface("Wall") = {out[2],out[3],out[4],out[5]};