DGtal  1.5.beta
Discrete differential calculus on polygonal surfaces
Author(s) of this documentation:
David Coeurjolly, Baptiste GENEST

Part of package DEC package.

In this documentation page, we detail the operators and tools for differential calculus computations on generic polygonal surfaces (PolygonalCalculus class). Such polygonal surfaces could be generic meshes (with non-planar, non-convex faces), or an embeddeing of a digital surface.

The overall idea, compared to the Discrete exterior calculus approach, is to define per face operators which can be later combined to define global ones.

This module documents the operators defined in [45]. An additional documentation page showcasing classical use cases of these operators for geometry processing on digital surfaces will be added.

The sign convention for the divergence and the Laplacian operator is opposite to the one of [45]. This is to match the usual mathematical convention that the Laplacian (and the Laplacian-Beltrami) has negative eigenvalues (and is the sum of second derivatives in the cartesian grid). It also follows the formal adjointness of exterior derivative and opposite of divergence as relation \( \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v \rangle \). See also https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
All illustrations below have been obtained using the DGtal+polyscope examples dgtalCalculus.cpp, dgtalCalculus-single.cpp and dgtalCalculus-poisson.cpp. To build these examples, enable the BUILD_POLYSCOPE_EXAMPLES variable (e.g. cmake .. -DBUILD_POLYSCOPE_EXAMPLES=ON).
The implementation heavily relies on implicit operators with many Eigen based small matrice constructions, which has a huge overhead in Debug mode. Please consider to build the examples in Release (e.g. CMAKE_BUILD_TYPE variable) for high performance on large geometrical objects.


The main idea of [45] is to define discrete operators per face, which can be combined to obtain a global operator acting on meshes. For vertex valued scalar or vector fields, this class implements elementary differential operators (gradient, derivatives, sharp/flat, Laplace-Beltrami...) which can be combined to define more complex ones, or to be used when solving variational geometry processing problems on meshes or digital surfaces.

The operators detailed [45] extend Mimetic Finite Difference and Virtual Element Method operators to arbitrary 3d polygons and discrete surfaces.

Per face operators

Let us consider a single (non-planar, non-convex) face example. Using

using SurfMesh = SurfaceMesh<Z3i::RealPoint,Z3i::RealVector>;
using namespace Z3i;
std::vector<RealPoint> positions={ {0,0,0},{20,0,0},{20,10,0},{10,8,5}, {0,15,1} };
//Single face
std::vector<std::vector<size_t>> faces={{ 0,1,2,3,4 }};
mesh = SurfMesh(positions.begin(),positions.end(),faces.begin(),faces.end());
PolygonalCalculus<SurfMesh> calculus(mesh);
SurfaceMesh< RealPoint, RealVector > SurfMesh
PolyCalculus * calculus

we obtain a SurfaceMesh instance with a unique face and its associated PolygonalCalculus object.

Example Example

Let us add a scalar function on vertices. For instance, using Eigen syntax, we can use;

phi << 1.0, 2.0, 0.0, 5.0 ,1.5;
LinAlg::DenseVector Vector
Type of Vector.
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
As the face is nonconvex / nonplanar, the vertex scalar quantity interpolation and the rendering in polyscope may be a bit odd.
Phi Interpolated phi

When constructing a PolygonalCalculus<SurfMesh> instance, you can provide an external function to embed the vertex positions. By default, the positions of the surface mesh vertices (positions in the previous example) is used. If you want to update the embedding (could be useful on digital surfaces), the user can specify the mapping (Face,Vertex)->RealPoint. Please refer to PolygonalCalculus<SurfMesh>::setEmbedder for an example.

As the Face id is a parameter of the embedder, a given vertex can have different embeddings for all its incident faces.

Basic operators

We first describe some standard per face operators. Note that for all extrinsic operators that require the vertex position in \( \mathbb{R}^3 \), the embedder is used. For these quantities, vector/point are encoded as Eigen 3x1 vectors and collections of n points as an \( n\times 3\) matrix.

Operator Output Description
calculus.X(f) \( n_f\times 3\) positions of the face vertices (n_f points in \( \mathbb{R}^3 \))
calculus.E(f) \( n_f\times 3\) the edge vectors of the face (n_f vectors in \( \mathbb{R}^3 \))
calculus.A(f) \( n_f\times n_f\) averaging operator that averages some values on vertices to values on edges (vectors \((p_{i+1}-p_{i})\in\mathbb{R}^3\))
calculus.B(f) \( n_f\times 3\) the edge mid-points of the face (n_f points in \( \mathbb{R}^3 \))
calculus.vectorArea(f) \( n_f\) vector the corrected vector area of the face (uses the embedder)
calculus.faceArea(f) scalar the corrected face area (uses the embedder)
calculus.faceNormalVector(f) vector in \(\mathbb{R}^3\) the corrected (unitary) normal vector of the face (uses the embedder)
calculus.faceNormalAsDGtalVector(f) RealVector the corrected (unitary) normal vector of the face (uses the embedder)
calculus.centroid(f) vector in \(\mathbb{R}^3\) the corrected centroid of the face
calculus.centroidAsDGtalVector(f) RealVector the corrected centroid of the face

Derivative operators

Derivative operators act on a scalar field defined on the vertex of the face.

Operator Output Description
calculus.D(f) \( n_f\times n_f\) derivative operator (d0 in DEC) mapping scalar values on vertices to values on edges (vertex scalars to a 1-form)
calculus.gradient(f) \( 3 \times n_f\) the gradient operator for vertex valued scalar functions of the face f
calculus.coGradient(f) \( 3 \times n_f\) the co-gradient operator for vertex valued scalar functions of the face f
A*phi D*phi (discrete 1-form)
gradient co-gradient corrected normal

Flat/Sharp and Projection Operators

Operator Output Description
calculus.flat(f) \( n_f\times 3\) Flat operator that maps an ambient vector to a local discrete 1-form
calculus.sharp(f) \( 3 \times n_f\) Sharp operator that maps a discrete 1-form to an ambient vector
calculus.P(f) \( n_f \times n_f\) Projection operator acting on a discrete 1-form that removes the associated non-tangential part of the vector reconstructed using the sharp operator

Discrete flat and sharp operators match with the DEC operators (see Discrete exterior calculus). For short, the flat operator embeds an ambient vector as a discrete 1-form (roughly values on the face edges) while the sharp operator pull-back a 1-form to reconstruct a vector in \(\mathbb{R}^3\). These operators are rank 2 and \( sharp_f flat_f s = s\) for any vector \(s\) in the face tangent plane, but \( flat_f sharp_f v = v\) for any vector v does not hold (due to rank deficiencies of the matrices). This is quantified by the projection operator \(P\) (see [45] for structural lemmas, e.g. \( normal_f\cdot (sharp_f flat_f s)=0\quad \forall s\in \mathbb{R}^3\)...).

For interested readers, [45] has a discussion about the relationship between these operators and the classical DEC or FEM ones.

a vector in \(\mathbb{R}^3\) a vector in \(\mathbb{R}^3\) Flat (1-form) Sharp of the flat (green vector) Sharp and normal vector (orange vector)

Inner Product and Laplace-Beltrami Operators

For 0-forms, the inner product is the classical one induced by the \( l_2\) norm. For 1-forms, the inner product is given by PolygonalCalculus<SurfMesh>::M() (useful to define the Laplace-Beltrami operator).

Operator Output Description
calculus.M(f,lambda) \( n_f\times n_f\) Inner product operator on discrete 1-forms operator (lambda is a regularization parameter, see [45])
calculus.laplaceBeltrami(f,lambda) \( n_f \times n_f\) Laplace-Beltrami operator \( D^tMD \) (lambda is a regularization parameter, see [45]), the matrix is PSD

Operators on directional fields and Connection Laplacian

In this section, we describe operators acting on directional fields (Levi-Civita connection, covariant gradient of a vector field –first-order derivative on VF–, and connection laplacian). Pleas refer to [45], section 5, for details.

As an example, these operators can be used to interpolate vector fields as illustrated in Vector Heat Method using discrete polygonal calculus.

Transport of a vertex based vector field to face tangent plane

Covariant Gradient and Projection

Since these operators are themselves matrices they cannot be constructed in the same way as the others ( Operator builder Matrix * function restricted to face Vector ), you need the specify the local vector field as well, in the same format as in [45], i.e.

\begin{eqnarray*} u_{f} = [u_{v_1}^t \ldots u_{v_{nf}}^t]^t \end{eqnarray*}

where \(u_{v_i}\) are the intrinsic vectors at each adjacent vertex expressed in their own tangent basis.

Operator Output Description
calculus.covariantGradient(f,uf) \( 2\times 2\) Covariant Gradient defined for intrinsic vector valued 0-forms at face f
calculus.covariantProjection(f,uf) \( n_f \times 2\) Covariant Projection operator defined for intrinsic vector valued 0-forms at face f

Vector Laplacian

Using the same format for vector fields as above, we can define a Vector Laplacian operator per face.

Operator Output Description
calculus.connectionLaplacian(f,lambda) \( 2n_f\times 2n_f\) Vector Laplacian defined as the associated Matrix with the Dirichlet energy for vector valued 0-forms (lambda is a regularization parameter, see [45]) at face f, PSD matrix

Global calculus

Given a scalar function defined on a generic surface mesh vertices, all previously mentioned operators can be applied to obtain consistent quantities on the overall mesh. For instance, from the dgtalCalculus.cpp example using Shortcuts and ShortcutsGeometry to set up the surface:

Surface Phi Gradient Gradient+co-gradient

To solve some global PDE (e.g. Laplace/Poisson problems, see below), one can combine the local operators into a global one, gathering the contributions of each face.

For example, the PolygonalCalculus<SurfMesh>::globalLaplaceBeltrami() method outputs a global (sparse) Laplace-Beltrami operator which can later be used for diffusion.

You can also retrieve the global mass matrix via PolygonalCalculus::globalLumpedMassMatrix, which is used to define an inner product for 0-forms (i.e. scalar functions), as well as its inverse PolygonalCalculus::globalInverseLumpedMassMatrix(), generally used when one solves a weak problem and wishes to get a pointwise per-vertex solution.

Corrected Calculus using Estimated Normal Vectors

On digital surfaces, solving PDE on original embedding with axis aligned quad surfaces may fail to correctly capture the surface metric. As discussed in [26], given an estimation of the tangent bundle of the discrete surface (for instance using estimated normal vectors from Integral invariant curvature estimator 2D/3D, or Digital Voronoi Covariance Measure and geometry estimation, cf Shortcuts (for the impatient developper)), one can implicitly project each face to a prescribed tangent plane and perform the computations on this new embedding of the geometry.

Geodesic distances without correction Geodesic distances with correction

The functor functors::EmbedderFromNormalVectors can be used to implicitly project Face vertices onto the prescribed tangent plane. A classical usage is the following one:

//A surface mesh. Eg. primal surface of a digital surface
SurfaceMesh< Z3i::RealPoint, Z3i::RealVector > surfmesh(...);
//Per face normal vector estimation
std::vector<Z3i::RealVector> ii_normals = ....
//New embedder using tangent plane projection of face vertices
functors::EmbedderFromNormalVectors<Z3i::RealPoint, Z3i::RealVector> embedderFromNormals(ii_normals,surfmesh);
//The calculus instance with the new embedder
PolygonalCalculus<Z3i::RealPoint,23i::RealVector> calculus(surfmesh);
calculus->setEmbedder( embedderFromNormals );
SurfMesh surfmesh

A complete code is given in the dgtalCalculus-geodesic.cpp example.

Example: Solving a Laplace problem

Let suppose we want to solve the following Laplace problem for data interpolation:

\begin{eqnarray*} \Delta_\Omega u& = 0 \\ & s.t. u = g \text{ on } \partial\Omega \end{eqnarray*}

We want to solve that problem on a polygonal mesh \(\Omega\) (digital surface here) with a boundary and some scalar values attached to boundary vertices, or sampled on the object surface.

Furthermore, the discrete version of the Laplace problem boils down to a simple linear problem using on the discrete Laplace-Beltrami sparse matrix.

We also use class DirichletConditions to enforce Dirichlet boundary conditions on the system.

The overall code is:

typedef PolygonalCalculus<SH3::RealPoint,SH3::RealVector> PolyDEC;
typedef DirichletConditions< EigenLinearAlgebraBackend > DC;
PolyDEC calculus(surfmesh);
PolyDEC::SparseMatrix L = calculus.globalLaplaceBeltrami();
PolyDEC::Form g = calculus.form0();
DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
//We set values on the boundary
auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
std::cout<< "Number of boundary edges= "<<boundaryEdges.size()<<std::endl;
auto pihVertex=[&](const SurfMesh::Vertex &v){return cos(scale*(surfmesh.position(v)[0]))*(scale*surfmesh.position(v)[1]);};
for(auto &e: boundaryEdges)
auto adjVertices = surfmesh.edgeVertices(e);
g(adjVertices.first) = pihVertex(adjVertices.first);
g(adjVertices.second) = pihVertex(adjVertices.second);
b(adjVertices.first) = 1;
b(adjVertices.second) = 1;
// Solve Δu=0 with g as boundary conditions
PolyDEC::Solver solver;
PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
solver.compute( L_dirichlet );
PolyDEC::Form g_dirichlet = DC::dirichletVector( L, g, b, g );
PolyDEC::Form x_dirichlet = solver.solve( g_dirichlet );
PolyDEC::Form u = DC::dirichletSolution( x_dirichlet, b, g );
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix

Leading to the following results (see dgtalCalculus-poisson.cpp):

Surface Boundary condition \( g\) Solution \( u \)

Vector Calculus

Global Vector Laplace/Poisson problems can also be solved by the same way, using instead PolygonalCalculus<SurfMesh>::globalConnectionLaplace() and PolygonalCalculus<SurfMesh>::doubledGlobalLumpedMassMatrix(). One can find examples of such use in the VectorsInHeat class.


Cache mechanisms and high-performance computing

The PolygonalCalculus class has two cache mechanisms:

  • An external cache strategy to store a given operator into a compact container. Typical use case is when the user wants to precompute a given operator, store it and efficiently reuse it while iterating over the faces. We detail this construction below.
  • The second one is a global internal cache strategy that will store all per face operators on the fly. In that case, each operator returning a DenseMatrix is stored in a cache the first time the calculus.operator(f) is called. Typical use case is when the user wants to use many times a large set of different operators. To enable this strategy, you can use the calculus.enableInternalGlobalCache(), or from the last parameter of the class constructor (boolean set to true). E.g.
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector> calculus(surfmesh,true); //global internal cache enabled.
    By default, this behavior is disabled as it is memory expensive (all operators are explicitly stored when used for the first time), and may not have a huge running time impact for some applications. An example is given in the dgtalCalculus-bunny.cpp. Once enabled, the class API remains the same, everything is transperent to the user.

We describe here the first external cache strategy. For the sake of readability, each operator has been implemented implicitly. For example, the M operator per face is given by

DenseMatrix M(const Face f, const double lambda=1.0) const
auto Uf=sharp(f);
auto Pf=P(f);
return faceArea(f) * Uf.transpose()*Uf + lambda * Pf.transpose()*Pf;
EigenLinearAlgebraBackend::DenseMatrix DenseMatrix
TriMesh::Face Face

which could be time consuming as the internal operators may be computed several times.

For high performance computations, we provide a generic cache mechanism to explicitly store all per face operators of a surface mesh (stored in a random access container).

A typical usage is

auto cacheSharp = getOperatorCacheMatrix( [&](Face f){ return(calculus.sharp(f);} );
auto cachefaceArea = getOperatorCacheMatrix( [&](Face f){ return(calculus.faceArea(f);} );
auto cacheP = getOperatorCacheMatrix( [&](Face f){ return(calculus.P(f);} );

Then, cached operators can be accessed and combined:

auto Mf = cachefaceArea[f] * cacheU[f].transpose()*cacheU[f] + lambda * cacheP[f].transpose() * cacheP[f];