PolygonalCalculus.h

 
#pragma once
 
// Inclusions
#include <iostream>
#include <functional>
#include <vector>
#include <string>
#include <map>
#include <unordered_map>
#include "DGtal/base/ConstAlias.h"
#include "DGtal/base/Common.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/math/linalg/EigenSupport.h"
 
namespace DGtal
{
 
// template class PolygonalCalculus
template <typename TRealPoint, typename TRealVector>
class PolygonalCalculus
{
  // ----------------------- Standard services ------------------------------
public:
  
  static const Dimension dimension = TRealPoint::dimension;
  BOOST_STATIC_ASSERT( ( dimension == 3 ) );
  
  typedef PolygonalCalculus<TRealPoint, TRealVector> Self;
  
  typedef SurfaceMesh<TRealPoint, TRealVector> MySurfaceMesh;
  typedef typename MySurfaceMesh::Vertex Vertex;
  typedef typename MySurfaceMesh::Face Face;
  typedef typename MySurfaceMesh::RealPoint Real3dPoint;
  typedef typename MySurfaceMesh::RealVector Real3dVector;
  
  typedef EigenLinearAlgebraBackend LinAlg;
  typedef LinAlg::DenseVector Vector;
  typedef Vector Form;
  typedef LinAlg::DenseMatrix DenseMatrix;
  typedef LinAlg::SparseMatrix SparseMatrix;
  typedef LinAlg::Triplet Triplet;
 
  typedef LinAlg::SolverSimplicialLDLT Solver;
 
  
  PolygonalCalculus(const ConstAlias<MySurfaceMesh> surf,
                    bool globalInternalCacheEnabled = false):
          mySurfaceMesh(&surf),  myGlobalCacheEnabled(globalInternalCacheEnabled)
  {
    myEmbedder =[&](Face f,Vertex v){ (void)f; return mySurfaceMesh->position(v);};
    myVertexNormalEmbedder = [&](Vertex v){ return toReal3dVector(computeVertexNormal(v));};
    init();
  };
  
  PolygonalCalculus(const ConstAlias<MySurfaceMesh> surf,
                    const std::function<Real3dPoint(Face,Vertex)> &embedder,
                    bool globalInternalCacheEnabled = false):
  mySurfaceMesh(&surf), myEmbedder(embedder), myGlobalCacheEnabled(globalInternalCacheEnabled)
  {
    myVertexNormalEmbedder = [&](Vertex v){ return toReal3dVector(computeVertexNormal(v)); };
    init();
  };
  
  PolygonalCalculus(const ConstAlias<MySurfaceMesh> surf,
                    const std::function<Real3dVector(Vertex)> &embedder,
                    bool globalInternalCacheEnabled = false):
      mySurfaceMesh(&surf), myVertexNormalEmbedder(embedder),
      myGlobalCacheEnabled(globalInternalCacheEnabled)
  {
    myEmbedder = [&](Face f,Vertex v){(void)f; return mySurfaceMesh->position(v); };
    init();
  };
 
  PolygonalCalculus(const ConstAlias<MySurfaceMesh> surf,
                    const std::function<Real3dPoint(Face,Vertex)> &pos_embedder,
                    const std::function<Vector(Vertex)> &normal_embedder,
                    bool globalInternalCacheEnabled = false) :
  mySurfaceMesh(&surf), myEmbedder(pos_embedder),
  myVertexNormalEmbedder(normal_embedder),
  myGlobalCacheEnabled(globalInternalCacheEnabled)
  {
    init();
  };
 
  PolygonalCalculus() = delete;
  
  ~PolygonalCalculus() = default;
  
  PolygonalCalculus ( const PolygonalCalculus & other ) = delete;
  
  PolygonalCalculus ( PolygonalCalculus && other ) = delete;
  
  PolygonalCalculus & operator= ( const PolygonalCalculus & other ) = delete;
  
  PolygonalCalculus & operator= ( PolygonalCalculus && other ) = delete;
 
  
  // ----------------------- embedding  services  --------------------------
  //MARK: Embedding services
  
  void setEmbedder(const std::function<Real3dPoint(Face,Vertex)> &externalFunctor)
  {
    myEmbedder = externalFunctor;
  }
  
 
  // ----------------------- Per face operators --------------------------------------
  //MARK: Per face operator on scalars
  
  DenseMatrix X(const Face f) const
  {
    if (checkCache(X_,f))
      return myGlobalCache[X_][f];
    
    const auto vertices = mySurfaceMesh->incidentVertices(f);
    const auto nf = myFaceDegree[f];
    DenseMatrix Xt(nf,3);
    size_t cpt=0;
    for(auto v: vertices)
    {
      Xt(cpt,0) = myEmbedder(f,v)[0];
      Xt(cpt,1) = myEmbedder(f,v)[1];
      Xt(cpt,2) = myEmbedder(f,v)[2];
      ++cpt;
    }
    
    setInCache(X_,f,Xt);
    return  Xt;
  }
  
 
  DenseMatrix D(const Face f) const
  {
    if (checkCache(D_,f))
      return myGlobalCache[D_][f];
    
    const auto nf = myFaceDegree[f];
    DenseMatrix d = DenseMatrix::Zero(nf ,nf);
    for(auto i=0u; i < nf; ++i)
    {
      d(i,i) = -1.;
      d(i, (i+1)%nf) = 1.;
    }
    
    setInCache(D_,f,d);
    return d;
  }
  
  DenseMatrix E(const Face f) const
  {
    if (checkCache(E_,f))
      return myGlobalCache[E_][f];
 
    DenseMatrix op = D(f)*X(f);
 
    setInCache(E_,f,op);
    return op;
  }
  
  DenseMatrix A(const Face f) const
  {
    if (checkCache(A_,f))
      return myGlobalCache[A_][f];
    
    const auto nf = myFaceDegree[f];
    DenseMatrix a = DenseMatrix::Zero(nf ,nf);
    for(auto i=0u; i < nf; ++i)
    {
      a(i, (i+1)%nf) = 0.5;
      a(i,i) = 0.5;
    }
 
    setInCache(A_,f,a);
    return a;
  }
  
  
  Vector vectorArea(const Face f) const
  {
    Real3dPoint af(0.0,0.0,0.0);
    const auto vertices = mySurfaceMesh->incidentVertices(f);
    auto it     = vertices.cbegin();
    auto itnext = vertices.cbegin();
    ++itnext;
    while (it != vertices.cend())
    {
      auto xi  = myEmbedder(f,*it);
      auto xip = myEmbedder(f,*itnext);
      af += xi.crossProduct(xip);
      ++it;
      ++itnext;
      if (itnext == vertices.cend())
        itnext = vertices.cbegin();
    }
    Eigen::Vector3d output = {af[0],af[1],af[2]};
    return 0.5*output;
  }
  
  double faceArea(const Face f) const
  {
    return vectorArea(f).norm();
  }
  
  Vector faceNormal(const Face f) const
  {
    Vector v = vectorArea(f);
    v.normalize();
    return v;
  }
  
  Real3dVector faceNormalAsDGtalVector(const Face f) const
  {
    Vector v = faceNormal(f);
    return {v(0),v(1),v(2)};
  }
  
  DenseMatrix coGradient(const Face f) const
  {
    if (checkCache(COGRAD_,f))
      return myGlobalCache[COGRAD_][f];
    DenseMatrix op = E(f).transpose() * A(f);
    setInCache(COGRAD_, f, op);
    return op;
  }
  
  DenseMatrix bracket(const Vector &n) const
  {
    DenseMatrix brack(3,3);
    brack << 0.0 , -n(2), n(1),
             n(2), 0.0 , -n(0),
            -n(1) , n(0),0.0 ;
    return brack;
  }
  
  DenseMatrix gradient(const Face f) const
  {
    if (checkCache(GRAD_,f))
      return myGlobalCache[GRAD_][f];
    
    DenseMatrix op = -1.0/faceArea(f) * bracket( faceNormal(f) ) * coGradient(f);
    
    setInCache(GRAD_,f,op);
    return op;
  }
  
  DenseMatrix  flat(const Face f) const
  {
    if (checkCache(FLAT_,f))
      return myGlobalCache[FLAT_][f];
    DenseMatrix n = faceNormal(f);
    DenseMatrix op = E(f)*( DenseMatrix::Identity(3,3) - n*n.transpose());
    setInCache(FLAT_,f,op);
    return op;
  }
  
  DenseMatrix B(const Face f) const
  {
    if (checkCache(B_,f))
      return myGlobalCache[B_][f];
    DenseMatrix res = A(f) * X(f);
    setInCache(B_,f,res);
    return res;
  }
  
  Vector centroid(const Face f) const
  {
    const auto nf = myFaceDegree[f];
    return 1.0/(double)nf * X(f).transpose() * Vector::Ones(nf);
  }
  
  Real3dPoint centroidAsDGtalPoint(const Face f) const
  {
    const Vector c = centroid(f);
    return {c(0),c(1),c(2)};
  }
  
  DenseMatrix sharp(const Face f) const
  {
    if (checkCache(SHARP_,f))
      return myGlobalCache[SHARP_][f];
 
    const auto  nf = myFaceDegree[f];
    DenseMatrix op = 1.0/faceArea(f) * bracket(faceNormal(f)) *
          ( B(f).transpose() - centroid(f)* Vector::Ones(nf).transpose() );
 
    setInCache(SHARP_,f,op);
    return op;
  }
  
  DenseMatrix P(const Face f) const
  {
    if (checkCache(P_,f))
      return myGlobalCache[P_][f];
    
    const auto  nf = myFaceDegree[f];
    DenseMatrix op = DenseMatrix::Identity(nf,nf) - flat(f)*sharp(f);
 
    setInCache(P_, f, op);
    return op;
  }
  
  DenseMatrix M(const Face f, const double lambda=1.0) const
  {
    if (checkCache(M_,f))
      return myGlobalCache[M_][f];
    
    DenseMatrix Uf = sharp(f);
    DenseMatrix Pf = P(f);
    DenseMatrix op = faceArea(f) * Uf.transpose()*Uf + lambda * Pf.transpose()*Pf;
    
    setInCache(M_,f,op);
    return op;
  }
  
  DenseMatrix divergence(const Face f) const
  {
    if (checkCache(DIVERGENCE_,f))
      return myGlobalCache[DIVERGENCE_][f];
 
    DenseMatrix op = -1.0 * D(f).transpose() * M(f);
    setInCache(DIVERGENCE_,f,op);
    
    return op;
  }
  
  DenseMatrix curl(const Face f) const
  {
    if (checkCache(CURL_,f))
      return myGlobalCache[CURL_][f];
    
    DenseMatrix op = DenseMatrix::Identity(myFaceDegree[f],myFaceDegree[f]);
 
    setInCache(CURL_,f,op);
    return op;
  }
  
  
  DenseMatrix laplaceBeltrami(const Face f, const double lambda=1.0) const
  {
    if (checkCache(L_,f))
      return myGlobalCache[L_][f];
 
    DenseMatrix Df = D(f);
    // Laplacian is a negative operator.
    DenseMatrix op = -1.0 * Df.transpose() * M(f,lambda) * Df;
    
    setInCache(L_, f, op);
    return op;
  }
  
  // ----------------------- Vector calculus----------------------------------
  //MARK: Vector Field Calculus
  
public:
  DenseMatrix Tv(const Vertex & v) const
  {
    Eigen::Vector3d nv = n_v(v);
    ASSERT(std::abs(nv.norm() - 1.0) < 0.001);
    const auto & N            = getSurfaceMeshPtr()->neighborVertices(v);
    auto neighbor             = *N.begin();
    Real3dPoint tangentVector = getSurfaceMeshPtr()->position(v) -
                                getSurfaceMeshPtr()->position(neighbor);
    Eigen::Vector3d w  = toVec3(tangentVector);
    Eigen::Vector3d uu = project(w,nv).normalized();
    Eigen::Vector3d vv = nv.cross(uu);
 
    DenseMatrix tanB(3,2);
    tanB.col(0) = uu;
    tanB.col(1) = vv;
    return tanB;
  }
 
  DenseMatrix Tf(const Face & f) const
  {
    Eigen::Vector3d nf = faceNormal(f);
    ASSERT(std::abs(nf.norm() - 1.0) < 0.001);
    const auto & N = getSurfaceMeshPtr()->incidentVertices(f);
    auto v1        = *(N.begin());
    auto v2        = *(N.begin() + 1);
    Real3dPoint tangentVector =
    getSurfaceMeshPtr()->position(v2) - getSurfaceMeshPtr()->position(v1);
    Eigen::Vector3d w  = toVec3(tangentVector);
    Eigen::Vector3d uu = project(w,nf).normalized();
    Eigen::Vector3d vv = nf.cross(uu);
 
    DenseMatrix tanB(3,2);
    tanB.col(0) = uu;
    tanB.col(1) = vv;
    return tanB;
  }
 
  Vector toExtrinsicVector(const Vertex v, const Vector & I) const
  {
    DenseMatrix T = Tv(v);
    return T.col(0) * I(0) + T.col(1) * I(1);
  }
 
  std::vector<Vector> toExtrinsicVectors(const std::vector<Vector> & I) const
  {
    std::vector<Vector> ext(mySurfaceMesh->nbVertices());
    for (auto v = 0; v < mySurfaceMesh->nbVertices(); v++)
      ext[v] = toExtrinsicVector(v,I[v]);
    return ext;
  }
 
  DenseMatrix Qvf(const Vertex & v, const Face & f) const
  {
    Eigen::Vector3d nf = faceNormal(f);
    Eigen::Vector3d nv = n_v(v);
    double c           = nv.dot(nf);
    ASSERT(std::abs( c + 1.0) > 0.0001);
    //Special case for opposite nv and nf vectors.
    if (std::abs( c + 1.0) < 0.00001)
      return -Eigen::Matrix3d::Identity();
  
    auto vv          = nv.cross(nf);
    DenseMatrix skew = bracket(vv);
    return Eigen::Matrix3d::Identity() + skew +
           1.0 / (1.0 + c) * skew * skew;
  }
 
  DenseMatrix Rvf(const Vertex & v, const Face & f) const
  {
    return Tf(f).transpose() * Qvf(v,f) * Tv(v);
  }
 
  DenseMatrix shapeOperator(const Face f) const
  {
    DenseMatrix N(myFaceDegree[f],3);
    uint cpt = 0;
    for (Vertex v : mySurfaceMesh->incidentVertices(f))
    {
      N.block(cpt,0,3,1) = n_v(v).transpose();
      cpt++;
    }
    DenseMatrix GN = gradient(f) * N, Tf = T_f(f);
 
    return 0.5 * Tf.transpose() * (GN + GN.transpose()) * Tf;
  }
 
  DenseMatrix transportAndFormatVectorField(const Face f, const Vector & uf)
  {
    DenseMatrix uf_nabla(myFaceDegree[f], 2);
    size_t cpt = 0;
    for (auto v : mySurfaceMesh->incidentVertices(f))
    {
      uf_nabla.block(cpt,0,1,2) =
      (Rvf(v,f) * uf.block(2 * cpt,0,2,1)).transpose();
      ++cpt;
    }
    return uf_nabla;
  }
 
  DenseMatrix covariantGradient(const Face f, const Vector & uf)
  {
    return Tf(f).transpose() * gradient(f) *
           transportAndFormatVectorField(f,uf);
  }
 
  DenseMatrix covariantProjection(const Face f, const Vector & uf)
  {
    return P(f) * D(f) * transportAndFormatVectorField(f,uf);
  }
  
  DenseMatrix covariantGradient_f(const Face & f) const
  {
    return kroneckerWithI2(Tf(f).transpose() * gradient(f)) * blockConnection(f);
  }
  
  DenseMatrix covariantProjection_f(const Face & f) const
  {
    return kroneckerWithI2(P(f) * D(f)) * blockConnection(f);
    ;
  }
  
  DenseMatrix connectionLaplacian(const Face & f, double lambda = 1.0) const
  {
    if (checkCache(CON_L_,f))
      return myGlobalCache[CON_L_][f];
    DenseMatrix G = covariantGradient_f(f);
    DenseMatrix P = covariantProjection_f(f);
    DenseMatrix L = -(faceArea(f) * G.transpose() * G + lambda * P.transpose() * P);
    setInCache(CON_L_,f,L);
    return L;
  }
  
  // ----------------------- Global operators --------------------------------------
  //MARK: Global Operators
 
  Form form0() const
  {
    return Form::Zero( nbVertices() );
  }
  SparseMatrix identity0() const
  {
    SparseMatrix Id0( nbVertices(), nbVertices() );
    Id0.setIdentity();
    return Id0; 
  }
  
  SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
  {
    SparseMatrix lapGlobal(mySurfaceMesh->nbVertices(), mySurfaceMesh->nbVertices());
    std::vector<Triplet> triplets;
    for( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces(); ++f )
    {
      auto nf = myFaceDegree[f];
      DenseMatrix Lap = this->laplaceBeltrami(f,lambda);
      const auto vertices = mySurfaceMesh->incidentVertices(f);
      for(auto i=0u; i < nf; ++i)
        for(auto j=0u; j < nf; ++j)
        {
          auto v = Lap(i,j);
          if (v!= 0.0)
            triplets.emplace_back( Triplet( (SparseMatrix::StorageIndex)vertices[ i ], (SparseMatrix::StorageIndex)vertices[ j ],
                                            Lap( i, j ) ) );
        }
    }
    lapGlobal.setFromTriplets(triplets.begin(), triplets.end());
    return lapGlobal;
  }
  
  SparseMatrix globalLumpedMassMatrix() const
  {
    SparseMatrix M(mySurfaceMesh->nbVertices(), mySurfaceMesh->nbVertices());
    std::vector<Triplet> triplets;
    for ( typename MySurfaceMesh::Index v = 0; v < mySurfaceMesh->nbVertices(); ++v )
      {
        auto faces = mySurfaceMesh->incidentFaces(v);
        auto varea = 0.0;
        for(auto f: faces)
          varea += faceArea(f) /(double)myFaceDegree[f];
        triplets.emplace_back(Triplet(v,v,varea));
      }
    M.setFromTriplets(triplets.begin(),triplets.end());
    return M;
  }
 
  SparseMatrix globalInverseLumpedMassMatrix() const
  {
    SparseMatrix iM0 = globalLumpedMassMatrix();
    for ( int k = 0; k < iM0.outerSize(); ++k )
      for ( typename SparseMatrix::InnerIterator it( iM0, k ); it; ++it )
        it.valueRef() = 1.0 / it.value();
    return iM0;
  }
 
  SparseMatrix globalConnectionLaplace(const double lambda = 1.0) const
  {
    auto nv = mySurfaceMesh->nbVertices();
    SparseMatrix lapGlobal(2 * nv, 2 * nv);
    SparseMatrix local(2 * nv, 2 * nv);
    std::vector<Triplet> triplets;
    for (typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces(); f++)
    {
      auto nf             = degree(f);
      DenseMatrix Lap     = connectionLaplacian(f,lambda);
      const auto vertices = mySurfaceMesh->incidentVertices(f);
      for (auto i = 0u; i < nf; ++i)
        for (auto j = 0u; j < nf; ++j)
          for (short k1 = 0; k1 < 2; k1++)
            for (short k2 = 0; k2 < 2; k2++)
            {
              auto v = Lap(2 * i + k1, 2 * j + k2);
              if (v != 0.0)
                triplets.emplace_back(Triplet(2 * vertices[i] + k1,
                                                2 * vertices[j] + k2, v));
            }
    }
    lapGlobal.setFromTriplets(triplets.begin(), triplets.end());
    return lapGlobal;
  }
 
  SparseMatrix doubledGlobalLumpedMassMatrix() const
  {
    auto nv = mySurfaceMesh->nbVertices();
    SparseMatrix M(2 * nv, 2 * nv);
    std::vector<Triplet> triplets;
    for (typename MySurfaceMesh::Index v = 0; v < mySurfaceMesh->nbVertices(); ++v)
    {
      auto faces = mySurfaceMesh->incidentFaces(v);
      auto varea = 0.0;
      for (auto f : faces)
        varea += faceArea(f) / (double)myFaceDegree[f];
      triplets.emplace_back(Triplet(2 * v, 2 * v, varea));
      triplets.emplace_back(Triplet(2 * v + 1, 2 * v + 1, varea));
    }
    M.setFromTriplets(triplets.begin(), triplets.end());
    return M;
  }
  
  // ----------------------- Cache mechanism --------------------------------------
  
  std::vector<DenseMatrix> getOperatorCacheMatrix(const std::function<DenseMatrix(Face)> &perFaceOperator) const
  {
    std::vector<DenseMatrix> cache;
    for( typename MySurfaceMesh::Index f=0; f < mySurfaceMesh->nbFaces(); ++f)
      cache.push_back(perFaceOperator(f));
    return cache;
  }
  
  std::vector<Vector> getOperatorCacheVector(const std::function<Vector(Face)> &perFaceVectorOperator) const
  {
    std::vector<Vector> cache;
    for( typename MySurfaceMesh::Index f=0; f < mySurfaceMesh->nbFaces(); ++f)
      cache.push_back(perFaceVectorOperator(f));
    return cache;
  }
 
  void enableInternalGlobalCache()
  {
    myGlobalCacheEnabled = true;
  }
  
  void disableInternalGlobalCache()
  {
    myGlobalCacheEnabled = false;
    myGlobalCache.clear();
  }
 
  
  // ----------------------- Common --------------------------------------
public:
  
  void init()
  {
    updateFaceDegree();
  }
  
  size_t faceDegree(Face f) const
  {
    return myFaceDegree[f];
  }
  
  size_t nbVertices() const
  {
    return mySurfaceMesh->nbVertices();
  }
  
  size_t nbFaces() const
  {
    return mySurfaceMesh->nbFaces();
  }
  
  size_t degree(const Face f) const
  {
    return myFaceDegree[f];
  }
  
  const MySurfaceMesh * getSurfaceMeshPtr() const
  {
    return mySurfaceMesh;
  }
  
  void selfDisplay ( std::ostream & out ) const
  {
    out << "[PolygonalCalculus]: ";
    if (myGlobalCacheEnabled)
      out<< "internal cache enabled, ";
    else
      out<<"internal cache disabled, ";
    out <<"SurfaceMesh="<<*mySurfaceMesh;
  }
  
  bool isValid() const
  {
    return true;
  }
 
  
  // ------------------------- Protected Datas ------------------------------
  //MARK: Protected
  
protected:
  
  enum OPERATOR { X_, D_, E_, A_, COGRAD_, GRAD_, FLAT_, B_, SHARP_, P_, M_, DIVERGENCE_, CURL_, L_,CON_L_};
  
  void updateFaceDegree()
  {
    myFaceDegree.resize(mySurfaceMesh->nbFaces());
    for(typename MySurfaceMesh::Index f = 0; f <  mySurfaceMesh->nbFaces(); ++f)
    {
      auto vertices = mySurfaceMesh->incidentVertices(f);
      auto nf = vertices.size();
      myFaceDegree[f] = nf;
    }
  }
  
  bool checkCache(OPERATOR key, const Face f) const
  {
    if (myGlobalCacheEnabled)
      if (myGlobalCache[key].find(f) != myGlobalCache[key].end())
        return true;
    return false;
  }
 
  void setInCache(OPERATOR key, const Face f,
                  const DenseMatrix &ope) const
  {
    if (myGlobalCacheEnabled)
      myGlobalCache[key][f]  = ope;
  }
  
  static Vector project(const Vector & u, const Vector & n)
  {
    return u - (u.dot(n) / n.squaredNorm()) * n;
  }
  
  static Vector toVector(const Eigen::Vector3d & x)
  {
    Vector X(3);
    for (int i = 0; i < 3; i++)
      X(i) = x(i);
    return X;
  }
  
  static Eigen::Vector3d toVec3(const Real3dPoint & x)
  {
    return Eigen::Vector3d(x(0),x(1),x(2));
  }
  
    static Real3dVector toReal3dVector(const Eigen::Vector3d & x)
    {
      return { x(0), x(1), x(2)};
    }
    
  
  Vector computeVertexNormal(const Vertex & v) const
  {
    Vector n(3);
    n(0) = 0.;
    n(1) = 0.;
    n(2) = 0.;
   /* for (auto f : mySurfaceMesh->incidentFaces(v))
      n += vectorArea(f);
    return n.normalized();
    */
    auto faces = mySurfaceMesh->incidentFaces(v);
    for (auto f : faces)
      n += vectorArea(f);
    
    if (fabs(n.norm() - 0.0) < 0.00001)
    {
      //On non-manifold edges touching the boundary, n may be null.
      trace.warning()<<"[PolygonalCalculus] Trying to compute the normal vector at a boundary vertex incident to pnon-manifold edge, we return a random vector."<<std::endl;
      n << Vector::Random(3);
    }
    n = n.normalized();
    return n;
  }
  
  Eigen::Vector3d n_v(const Vertex & v) const
  {
    return toVec3(myVertexNormalEmbedder(v));
  }
 
  //Covariant operators routines
  DenseMatrix blockConnection(const Face & f) const
  {
    auto nf           = degree(f);
    DenseMatrix RU_fO = DenseMatrix::Zero(nf * 2,nf * 2);
    size_t cpt        = 0;
    for (auto v : getSurfaceMeshPtr()->incidentVertices(f))
    {
      auto Rv                               = Rvf(v,f);
      RU_fO.block<2,2>(2 * cpt,2 * cpt) = Rv;
      ++cpt;
    }
    return RU_fO;
  }
  
  DenseMatrix kroneckerWithI2(const DenseMatrix & M) const
  {
    size_t h       = M.rows();
    size_t w       = M.cols();
    DenseMatrix MK = DenseMatrix::Zero(h * 2,w * 2);
    for (size_t j = 0; j < h; j++)
      for (size_t i = 0; i < w; i++)
      {
        MK(2 * j, 2 * i)         = M(j, i);
        MK(2 * j + 1, 2 * i + 1) = M(j, i);
      }
    return MK;
  }
 
  
  
  
  // ------------------------- Internals ------------------------------------
  //MARK: Internals
private:
  
  const MySurfaceMesh *mySurfaceMesh;
  
  std::function<Real3dPoint(Face, Vertex)> myEmbedder;
 
  std::function<Real3dVector(Vertex)> myVertexNormalEmbedder;
  
  std::vector<size_t> myFaceDegree;
    
  bool myGlobalCacheEnabled;
  mutable std::array<std::unordered_map<Face,DenseMatrix>, 15> myGlobalCache;
  
}; // end of class PolygonalCalculus
 
template <typename TP, typename TV>
std::ostream&
operator<< ( std::ostream & out, const PolygonalCalculus<TP,TV> & object )
{
  object.selfDisplay( out );
  return out;
}
 
} // namespace DGtal