positions = primalSurface->positions();
+ for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
+ faces.push_back(primalSurface->incidentVertices( face ));
+
+ surfmesh = SurfMesh(positions.begin(),
+ positions.end(),
+ faces.begin(),
+ faces.end());
+ surfmesh.computeFaceNormalsFromPositions();
+ psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
+
+ // Initialize polyscope
+ polyscope::init();
+
+ // Set the callback function
+ polyscope::state::userCallback = myCallback;
+ polyscope::show();
+ return EXIT_SUCCESS;
+}
diff --git a/src/DGtal/dec/InterpolatedCorrectedCalculus.h b/src/DGtal/dec/InterpolatedCorrectedCalculus.h
new file mode 100644
index 0000000000..75dd7d6666
--- /dev/null
+++ b/src/DGtal/dec/InterpolatedCorrectedCalculus.h
@@ -0,0 +1,811 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+#pragma once
+
+/**
+ * @file
+ * @author
+ *
+ * @date 2024/06/21
+ *
+ * Header file for module SurfaceDEC.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+
+#if !defined InterpolatedCorrectedCalculus_h
+/** Prevents repeated inclusion of headers. */
+#define InterpolatedCorrectedCalculus_h
+
+#include "DGtal/base/Common.h"
+#include "DGtal/dec/SurfaceDEC.h"
+#include "DGtal/shapes/SurfaceMesh.h"
+
+namespace DGtal
+{
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class InterpolatedCorrectedCalculus
+ /**
+ Description of template class 'InterpolatedCorrectedCalculus'
+ \brief Implement differential operators on digital surfaces with
+ a vertex normal field.
+
+ Works digital surfaces.
+
+ Requires vertex normals on attached surface mesh.
+
+ Implements SurfaceDEC methods. Use them to access the corresponding
+ operators.
+
+ * @tparam TLinearAlgebraBackend linear algebra backend used (i.e.
+ EigenSparseLinearAlgebraBackend).
+ @tparam TRealPoint an arbitrary model of 3D RealPoint.
+ @tparam TRealVector an arbitrary model of 3D RealVector.
+ */
+ template
+ class InterpolatedCorrectedCalculus
+ : public SurfaceDEC,
+ TLinearAlgebraBackend>
+ {
+ friend class SurfaceDEC;
+
+ public:
+ typedef TLinearAlgebraBackend LinearAlgebraBackend;
+ typedef TRealPoint RealPoint;
+ typedef TRealVector RealVector;
+ typedef typename LinearAlgebraBackend::DenseVector::Index Index;
+ typedef typename LinearAlgebraBackend::DenseVector::Scalar Scalar;
+ typedef typename LinearAlgebraBackend::SparseMatrix LinearOperator;
+ typedef
+ typename LinearAlgebraBackend::SparseMatrix::StorageIndex StorageIndex;
+ typedef typename LinearAlgebraBackend::DenseMatrix DenseMatrix;
+ typedef typename LinearAlgebraBackend::DenseVector DenseVector;
+ typedef SurfaceMesh Mesh;
+ typedef typename Mesh::Vertex Vertex;
+ typedef typename Mesh::Edge Edge;
+ typedef typename Mesh::Face Face;
+ typedef typename Mesh::Face Corner;
+ typedef typename Mesh::Size Size;
+
+ /// A triplet (row, col, value), useful for initializaing sparse matrices.
+ typedef typename LinearAlgebraBackend::Triplet Triplet;
+ /// A range of triplets (row,col,value).
+ typedef std::vector Triplets;
+
+ private:
+ const Mesh * myMesh;
+ bool use2ndOrder = false;
+ double lambda = 1.;
+
+ DenseMatrix quadrangleSharpLocalMatrix( Face f ) const
+ {
+ DenseMatrix Basis = DenseMatrix::Zero( 3, 3 );
+ DenseMatrix Gavg = DenseMatrix::Zero( 3, 4 );
+ DenseMatrix Gedge = DenseMatrix::Zero( 3, 4 );
+ DenseMatrix S = DenseMatrix::Zero( 3, 4 );
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto i = vtcs[ 0 ];
+ const auto j = vtcs[ 1 ];
+ const auto k = vtcs[ 2 ];
+ const auto l = vtcs[ 3 ];
+ const auto dxij = myMesh->position( j ) - myMesh->position( i );
+ const auto dxjk = myMesh->position( k ) - myMesh->position( j );
+ const auto xij = dxij.getNormalized();
+ const auto xjk = dxjk.getNormalized();
+
+ const RealVector n = xij.crossProduct( xjk );
+ const Dimension x = fabs( xij[ 0 ] ) > fabs( xij[ 1 ] )
+ ? ( fabs( xij[ 0 ] ) > fabs( xij[ 2 ] ) ? 0 : 2 )
+ : ( fabs( xij[ 1 ] ) > fabs( xij[ 2 ] ) ? 1 : 2 );
+ const Dimension y = fabs( xjk[ 0 ] ) > fabs( xjk[ 1 ] )
+ ? ( fabs( xjk[ 0 ] ) > fabs( xjk[ 2 ] ) ? 0 : 2 )
+ : ( fabs( xjk[ 1 ] ) > fabs( xjk[ 2 ] ) ? 1 : 2 );
+ const Dimension z = fabs( n[ 0 ] ) > fabs( n[ 1 ] )
+ ? ( fabs( n[ 0 ] ) > fabs( n[ 2 ] ) ? 0 : 2 )
+ : ( fabs( n[ 1 ] ) > fabs( n[ 2 ] ) ? 1 : 2 );
+ if ( ( x == y ) || ( x == z ) || ( y == z ) || ( x + y + z > 3 ) )
+ trace.warning() << "Bad basis: " << x << y << z << std::endl;
+ const Scalar sign_x = xij[ x ] > 0.0 ? 1.0 : -1.0;
+ const Scalar sign_y = xjk[ y ] > 0.0 ? 1.0 : -1.0;
+ const Scalar sign_z = n[ z ] > 0.0 ? 1.0 : -1.0;
+
+ if ( use2ndOrder == false )
+ {
+ const auto ui = myMesh->vertexNormal( i );
+ const auto uj = myMesh->vertexNormal( j );
+ const auto uk = myMesh->vertexNormal( k );
+ const auto ul = myMesh->vertexNormal( l );
+ const auto u0 = 0.5 * ( ui + uj );
+ const auto u1 = 0.5 * ( uj + uk );
+ const auto u2 = 0.5 * ( uk + ul );
+ const auto u3 = 0.5 * ( ul + ui );
+ const auto u = 0.5 * ( u0 + u2 );
+ // We put u into the "natural" frame of the surfel
+ const auto ux = xij.dot( u );
+ const auto uy = xjk.dot( u );
+ const auto uz = n.dot( u );
+ Gavg.coeffRef( 0, 0 ) = uz;
+ Gavg.coeffRef( 0, 2 ) = -uz;
+ Gavg.coeffRef( 1, 1 ) = uz;
+ Gavg.coeffRef( 1, 3 ) = -uz;
+ Gavg.coeffRef( 2, 0 ) = -ux;
+ Gavg.coeffRef( 2, 1 ) = -uy;
+ Gavg.coeffRef( 2, 2 ) = ux;
+ Gavg.coeffRef( 2, 3 ) = uy;
+ Gavg /= 3.0;
+ Gedge.coeffRef( 0, 0 ) = n.dot( u0 ); // u0[z];
+ Gedge.coeffRef( 0, 2 ) = -n.dot( u2 ); // -u2[z];
+ Gedge.coeffRef( 1, 1 ) = n.dot( u1 ); // u1[z];
+ Gedge.coeffRef( 1, 3 ) = -n.dot( u3 ); // -u3[z];
+ Gedge.coeffRef( 2, 0 ) = -xij.dot( u0 ); // -u0[x];
+ Gedge.coeffRef( 2, 1 ) = -xjk.dot( u1 ); // -u1[y];
+ Gedge.coeffRef( 2, 2 ) = xij.dot( u2 ); // u2[x];
+ Gedge.coeffRef( 2, 3 ) = xjk.dot( u3 ); // u3[y];
+ Gedge /= 6.0;
+ // We compute the change of basis to come back to the canonic frame of
+ // R3.
+ S = Gavg + Gedge;
+ }
+ else
+ { // use2ndOrder
+ const auto u00 = myMesh->vertexNormal( i );
+ const auto u10 = myMesh->vertexNormal( j );
+ const auto u11 = myMesh->vertexNormal( k );
+ const auto u01 = myMesh->vertexNormal( l );
+ const auto ux00 = xij.dot( u00 );
+ const auto ux10 = xij.dot( u10 );
+ const auto ux11 = xij.dot( u11 );
+ const auto ux01 = xij.dot( u01 );
+ const auto uy00 = xjk.dot( u00 );
+ const auto uy10 = xjk.dot( u10 );
+ const auto uy11 = xjk.dot( u11 );
+ const auto uy01 = xjk.dot( u01 );
+ const auto uz00 = n.dot( u00 );
+ const auto uz10 = n.dot( u10 );
+ const auto uz11 = n.dot( u11 );
+ const auto uz01 = n.dot( u01 );
+ const auto a = 1.0 / ( u00 + u10 ).norm();
+ const auto b = 1.0 / ( u10 + u11 ).norm();
+ const auto c = 1.0 / ( u11 + u01 ).norm();
+ const auto d = 1.0 / ( u01 + u00 ).norm();
+ const auto e = 1.0 / ( u00 + u10 + u11 + u01 ).norm();
+ S.coeffRef( 0, 0 ) =
+ ( 4 * a + 2 * d + 8 * e + 1 ) * uz00 + 2 * ( d + 4 * e ) * uz01 +
+ ( 4 * a + 2 * b + 8 * e + 1 ) * uz10 + 2 * ( b + 4 * e ) * uz11;
+ S.coeffRef( 0, 2 ) =
+ -2 * ( d + 4 * e ) * uz00 - ( 4 * c + 2 * d + 8 * e + 1 ) * uz01 -
+ 2 * ( b + 4 * e ) * uz10 - ( 2 * b + 4 * c + 8 * e + 1 ) * uz11;
+ S.coeffRef( 1, 1 ) = 2 * ( a + 4 * e ) * uz00 +
+ 2 * ( c + 4 * e ) * uz01 +
+ ( 2 * a + 4 * b + 8 * e + 1 ) * uz10 +
+ ( 4 * b + 2 * c + 8 * e + 1 ) * uz11;
+ S.coeffRef( 1, 3 ) = -( 2 * a + 4 * d + 8 * e + 1 ) * uz00 -
+ ( 2 * c + 4 * d + 8 * e + 1 ) * uz01 -
+ 2 * ( a + 4 * e ) * uz10 -
+ 2 * ( c + 4 * e ) * uz11;
+ S.coeffRef( 2, 0 ) =
+ -( 4 * a + 2 * d + 8 * e + 1 ) * ux00 - 2 * ( d + 4 * e ) * ux01 -
+ ( 4 * a + 2 * b + 8 * e + 1 ) * ux10 - 2 * ( b + 4 * e ) * ux11;
+ S.coeffRef( 2, 1 ) = -2 * ( a + 4 * e ) * uy00 -
+ 2 * ( c + 4 * e ) * uy01 -
+ ( 2 * a + 4 * b + 8 * e + 1 ) * uy10 -
+ ( 4 * b + 2 * c + 8 * e + 1 ) * uy11;
+ S.coeffRef( 2, 2 ) =
+ 2 * ( d + 4 * e ) * ux00 + ( 4 * c + 2 * d + 8 * e + 1 ) * ux01 +
+ 2 * ( b + 4 * e ) * ux10 + ( 2 * b + 4 * c + 8 * e + 1 ) * ux11;
+ S.coeffRef( 2, 3 ) = ( 2 * a + 4 * d + 8 * e + 1 ) * uy00 +
+ ( 2 * c + 4 * d + 8 * e + 1 ) * uy01 +
+ 2 * ( a + 4 * e ) * uy10 +
+ 2 * ( c + 4 * e ) * uy11;
+ S /= 36.0;
+ }
+ const Scalar h = fabs( dxij[ x ] );
+ Basis.coeffRef( 0, x ) = sign_x * h;
+ Basis.coeffRef( 1, y ) = sign_y * h;
+ Basis.coeffRef( 2, z ) = sign_z * h;
+ const auto a_f = quadrangleArea( f );
+ return ( Basis.transpose() / a_f ) * S;
+ }
+
+ DenseMatrix quadrangleFlatLocalMatrix( Face f ) const
+ {
+ DenseMatrix Basis = DenseMatrix::Zero( 3, 3 );
+ DenseMatrix V = DenseMatrix::Zero( 4, 3 );
+ const auto vtcs = myMesh->incidentVertices( f );
+ if ( vtcs.size() != 4 )
+ {
+ trace.error()
+ << "[SurfaceMeshDEC::quadrangleCorrectedDigitalFlatLocalMatrix]"
+ << " restricted to quadrangle faces." << std::endl;
+ return V;
+ }
+ const auto i = vtcs[ 0 ];
+ const auto j = vtcs[ 1 ];
+ const auto k = vtcs[ 2 ];
+ const auto l = vtcs[ 3 ];
+ const auto dxij = myMesh->position( j ) - myMesh->position( i );
+ const auto dxjk = myMesh->position( k ) - myMesh->position( j );
+ const auto xij = dxij.getNormalized();
+ const auto xjk = dxjk.getNormalized();
+
+ const RealVector n = xij.crossProduct( xjk );
+ const Dimension x = fabs( xij[ 0 ] ) > fabs( xij[ 1 ] )
+ ? ( fabs( xij[ 0 ] ) > fabs( xij[ 2 ] ) ? 0 : 2 )
+ : ( fabs( xij[ 1 ] ) > fabs( xij[ 2 ] ) ? 1 : 2 );
+ const Dimension y = fabs( xjk[ 0 ] ) > fabs( xjk[ 1 ] )
+ ? ( fabs( xjk[ 0 ] ) > fabs( xjk[ 2 ] ) ? 0 : 2 )
+ : ( fabs( xjk[ 1 ] ) > fabs( xjk[ 2 ] ) ? 1 : 2 );
+ const Dimension z = fabs( n[ 0 ] ) > fabs( n[ 1 ] )
+ ? ( fabs( n[ 0 ] ) > fabs( n[ 2 ] ) ? 0 : 2 )
+ : ( fabs( n[ 1 ] ) > fabs( n[ 2 ] ) ? 1 : 2 );
+ if ( ( x == y ) || ( x == z ) || ( y == z ) || ( x + y + z > 3 ) )
+ trace.warning() << "Bad basis: " << x << y << z << std::endl;
+ const Scalar sign_x = xij[ x ] > 0.0 ? 1.0 : -1.0;
+ const Scalar sign_y = xjk[ y ] > 0.0 ? 1.0 : -1.0;
+ const Scalar sign_z = n[ z ] > 0.0 ? 1.0 : -1.0;
+ const auto ui = myMesh->vertexNormal( i );
+ const auto uj = myMesh->vertexNormal( j );
+ const auto uk = myMesh->vertexNormal( k );
+ const auto ul = myMesh->vertexNormal( l );
+ // We put all the u into the "natural" frame of the surfel
+ const auto uix = xij.dot( ui );
+ const auto uiy = xjk.dot( ui );
+ const auto uiz = n.dot( ui );
+ const auto ujx = xij.dot( uj );
+ const auto ujy = xjk.dot( uj );
+ const auto ujz = n.dot( uj );
+ const auto ukx = xij.dot( uk );
+ const auto uky = xjk.dot( uk );
+ const auto ukz = n.dot( uk );
+ const auto ulx = xij.dot( ul );
+ const auto uly = xjk.dot( ul );
+ const auto ulz = n.dot( ul );
+ if ( !use2ndOrder )
+ {
+ V.coeffRef( 0, 0 ) =
+ 6.0 - 2.0 * ( 2 * uix + ujx ) * uix - 2.0 * ( uix + 2 * ujx ) * ujx;
+ V.coeffRef( 0, 1 ) =
+ -2.0 * ( 2 * uix + ujx ) * uiy - 2.0 * ( uix + 2 * ujx ) * ujy;
+ V.coeffRef( 0, 2 ) =
+ -2.0 * ( 2 * uix + ujx ) * uiz - 2.0 * ( uix + 2 * ujx ) * ujz;
+ V.coeffRef( 1, 0 ) =
+ -2.0 * ( 2 * ujy + uky ) * ujx - 2.0 * ( ujy + 2 * uky ) * ukx;
+ V.coeffRef( 1, 1 ) =
+ 6.0 - 2.0 * ( 2 * ujy + uky ) * ujy - 2.0 * ( ujy + 2 * uky ) * uky;
+ V.coeffRef( 1, 2 ) =
+ -2.0 * ( 2 * ujy + uky ) * ujz - 2.0 * ( ujy + 2 * uky ) * ukz;
+ V.coeffRef( 2, 0 ) =
+ -6.0 + 2.0 * ( 2 * ukx + ulx ) * ukx + 2.0 * ( ukx + 2 * ulx ) * ulx;
+ V.coeffRef( 2, 1 ) =
+ 2.0 * ( 2 * ukx + ulx ) * uky + 2.0 * ( ukx + 2 * ulx ) * uly;
+ V.coeffRef( 2, 2 ) =
+ 2.0 * ( 2 * ukx + ulx ) * ukz + 2.0 * ( ukx + 2 * ulx ) * ulz;
+ V.coeffRef( 3, 0 ) =
+ 2.0 * ( 2 * uly + uiy ) * ulx + 2.0 * ( uly + 2 * uiy ) * uix;
+ V.coeffRef( 3, 1 ) =
+ -6.0 + 2.0 * ( 2 * uly + uiy ) * uly + 2.0 * ( uly + 2 * uiy ) * uiy;
+ V.coeffRef( 3, 2 ) =
+ 2.0 * ( 2 * uly + uiy ) * ulz + 2.0 * ( uly + 2 * uiy ) * uiz;
+ V /= 6.0;
+ }
+ else
+ {
+ const auto a = 1.0 / ( ui + uj ).norm();
+ const auto b = 1.0 / ( uj + uk ).norm();
+ const auto c = 1.0 / ( uk + ul ).norm();
+ const auto d = 1.0 / ( ul + ui ).norm();
+ const auto e = 1.0 / ( ui + uj + uk + ul ).norm();
+ const auto ap = 4 * a * a + a + 1;
+ const auto app = 16 * a * a + 4 * a - 1;
+ const auto bp = 4 * b * b + b + 1;
+ const auto bpp = 16 * b * b + 4 * b - 1;
+ const auto cp = 4 * c * c + c + 1;
+ const auto cpp = 16 * c * c + 4 * c - 1;
+ const auto dp = 4 * d * d + d + 1;
+ const auto dpp = 16 * d * d + 4 * d - 1;
+ V.coeffRef( 0, 0 ) =
+ -4 * ap * uix * uix - 2 * app * uix * ujx - 4 * ap * ujx * ujx + 30.0;
+ V.coeffRef( 0, 1 ) = -( 4 * ap * uix + app * ujx ) * uiy -
+ ( app * uix + 4 * ap * ujx ) * ujy;
+ V.coeffRef( 0, 2 ) = -( 4 * ap * uix + app * ujx ) * uiz -
+ ( app * uix + 4 * ap * ujx ) * ujz;
+ V.coeffRef( 1, 0 ) = -( 4 * bp * ujx + bpp * ukx ) * ujy -
+ ( bpp * ujx + 4 * bp * ukx ) * uky;
+ V.coeffRef( 1, 1 ) =
+ -4 * bp * ujy * ujy - 2 * bpp * ujy * uky - 4 * bp * uky * uky + 30.0;
+ V.coeffRef( 1, 2 ) = -( 4 * bp * ujy + bpp * uky ) * ujz -
+ ( bpp * ujy + 4 * bp * uky ) * ukz;
+ V.coeffRef( 2, 0 ) =
+ 4 * cp * ulx * ulx + 2 * cpp * ulx * ukx + 4 * cp * ukx * ukx - 30.0;
+ V.coeffRef( 2, 1 ) =
+ ( 4 * cp * ulx + cpp * ukx ) * uly + ( cpp * ulx + 4 * cp * ukx ) * uky;
+ V.coeffRef( 2, 2 ) =
+ ( 4 * cp * ulx + cpp * ukx ) * ulz + ( cpp * ulx + 4 * cp * ukx ) * ukz;
+ V.coeffRef( 3, 0 ) =
+ ( 4 * dp * uix + dpp * ulx ) * uiy + ( dpp * uix + 4 * dp * ulx ) * uly;
+ V.coeffRef( 3, 1 ) =
+ 4 * dp * uiy * uiy + 2 * dpp * uiy * uly + 4 * dp * uly * uly - 30.0;
+ V.coeffRef( 3, 2 ) =
+ ( 4 * dp * uiy + dpp * uly ) * uiz + ( dpp * uiy + 4 * dp * uly ) * ulz;
+ V /= 30.0;
+ }
+
+ // We compute the change of basis to come back to the canonic frame of R3.
+ const Scalar h = fabs( dxij[ x ] );
+ Basis.coeffRef( 0, x ) = sign_x * h;
+ Basis.coeffRef( 1, y ) = sign_y * h;
+ Basis.coeffRef( 2, z ) = sign_z * h;
+ return V * Basis;
+ }
+
+ //-----------------------------------------------------------------------------
+ DenseMatrix quadrangleInnerProduct0( Face f ) const
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto i = vtcs[ 0 ];
+ const auto j = vtcs[ 1 ];
+ const auto k = vtcs[ 2 ];
+ const auto l = vtcs[ 3 ];
+ const auto xij = myMesh->position( j ) - myMesh->position( i );
+ const auto xil = myMesh->position( l ) - myMesh->position( i );
+ const auto xlk = myMesh->position( k ) - myMesh->position( l );
+ const auto xjk = myMesh->position( k ) - myMesh->position( j );
+ const RealVector nn[ 4 ] = { /* ni */ xij.crossProduct( xil ),
+ /* nj */ xij.crossProduct( xjk ),
+ /* nk */ xlk.crossProduct( xjk ),
+ /* nl */ xlk.crossProduct( xil ) };
+ const Scalar factor = 1.0 / 3600.0;
+ const Vertex id[ 4 ] = { i, j, k, l };
+ DenseMatrix M = DenseMatrix::Zero( 4, 4 );
+ if ( !use2ndOrder )
+ {
+ // All coefs are computed from polynomial of the form:
+ // s^a (1-s)^a t^b (1-t)^b
+ // Indices are 0,1,2,3 == i,j,k,l
+ const int pa[ 4 ] = { 0, 1, 1, 0 };
+ const int pb[ 4 ] = { 0, 0, 1, 1 };
+ const double co[ 5 ][ 5 ] = { { 144.0, 36.0, 24.0, 36.0, 144.0 },
+ { 36.0, 9.0, 6.0, 9.0, 36.0 },
+ { 24.0, 6.0, 4.0, 6.0, 24.0 },
+ { 36.0, 9.0, 6.0, 9.0, 36.0 },
+ { 144.0, 36.0, 24.0, 36.0, 144.0 } };
+ for ( int g = 0; g < 4; ++g ) // left 0-form at vertex
+ for ( int h = 0; h < 4; ++h ) // right 0-form at vertex
+ { // Computing coefficients for ( id[ g ], id[ h ] )
+ for ( int u = 0; u < 4; ++u ) // corrected normal at vertex
+ for ( int n = 0; n < 4; ++n ) // naive normal at vertex
+ {
+ const int a = pa[ g ] + pa[ h ] + pa[ u ] + pa[ n ];
+ const int b = pb[ g ] + pb[ h ] + pb[ u ] + pb[ n ];
+ const double coeff =
+ factor * co[ a ][ b ] *
+ myMesh->vertexNormal( id[ u ] ).dot( nn[ n ] );
+ M( g, h ) += coeff;
+ }
+ }
+ }
+ else
+ { // use2ndOrder
+ // # 1ere ligne de M0 (avec kle facteur 3600*)
+ const auto ui = myMesh->vertexNormal( i );
+ const auto uj = myMesh->vertexNormal( j );
+ const auto uk = myMesh->vertexNormal( k );
+ const auto ul = myMesh->vertexNormal( l );
+ // formula is theoretically valid for constant geometric normals nn
+ const auto uiz = nn[ 0 ].dot( myMesh->vertexNormal( i ) );
+ const auto ujz = nn[ 1 ].dot( myMesh->vertexNormal( j ) );
+ const auto ukz = nn[ 2 ].dot( myMesh->vertexNormal( k ) );
+ const auto ulz = nn[ 3 ].dot( myMesh->vertexNormal( l ) );
+ const auto a = 1.0 / ( ui + uj ).norm();
+ const auto b = 1.0 / ( uj + uk ).norm();
+ const auto c = 1.0 / ( uk + ul ).norm();
+ const auto d = 1.0 / ( ul + ui ).norm();
+ const auto e = 1.0 / ( ui + uj + uk + ul ).norm();
+ M.coeffRef( 0, 0 ) = ( 108 * a + 108 * d + 144 * e + 81 ) * uiz +
+ ( 108 * a - 12 * b + 144 * e - 9 ) * ujz +
+ ( -12 * b - 12 * c + 144 * e + 1 ) * ukz +
+ ( -12 * c + 108 * d + 144 * e - 9 ) * ulz;
+ M.coeffRef( 0, 1 ) = ( 72 * a + 12 * d + 96 * e + 9 ) * uiz +
+ ( 72 * a + 12 * b + 96 * e + 9 ) * ujz +
+ ( 12 * b - 8 * c + 96 * e - 1 ) * ukz +
+ ( -8 * c + 12 * d + 96 * e - 1 ) * ulz;
+ M.coeffRef( 0, 2 ) = ( 8 * a + 8 * d + 64 * e + 1 ) * uiz +
+ ( 8 * a + 8 * b + 64 * e + 1 ) * ujz +
+ ( 8 * b + 8 * c + 64 * e + 1 ) * ukz +
+ ( 8 * c + 8 * d + 64 * e + 1 ) * ulz;
+ M.coeffRef( 0, 3 ) = ( 12 * a + 72 * d + 96 * e + 9 ) * uiz +
+ ( 12 * a - 8 * b + 96 * e - 1 ) * ujz +
+ ( -8 * b + 12 * c + 96 * e - 1 ) * ukz +
+ ( 12 * c + 72 * d + 96 * e + 9 ) * ulz;
+ M.coeffRef( 1, 0 ) = M( 0, 1 );
+ M.coeffRef( 1, 1 ) = ( 108 * a - 12 * d + 144 * e - 9 ) * uiz +
+ ( 108 * a + 108 * b + 144 * e + 81 ) * ujz +
+ ( 108 * b - 12 * c + 144 * e - 9 ) * ukz +
+ ( -12 * c - 12 * d + 144 * e + 1 ) * ulz;
+ M.coeffRef( 1, 2 ) = ( 12 * a - 8 * d + 96 * e - 1 ) * uiz +
+ ( 12 * a + 72 * b + 96 * e + 9 ) * ujz +
+ ( 72 * b + 12 * c + 96 * e + 9 ) * ukz +
+ ( 12 * c - 8 * d + 96 * e - 1 ) * ulz;
+ M.coeffRef( 1, 3 ) = ( 8 * a + 8 * d + 64 * e + 1 ) * uiz +
+ ( 8 * a + 8 * b + 64 * e + 1 ) * ujz +
+ ( 8 * b + 8 * c + 64 * e + 1 ) * ukz +
+ ( 8 * c + 8 * d + 64 * e + 1 ) * ulz;
+ M.coeffRef( 2, 0 ) = M( 0, 2 );
+ M.coeffRef( 2, 1 ) = M( 1, 2 );
+ M.coeffRef( 2, 2 ) = ( -12 * a - 12 * d + 144 * e + 1 ) * uiz +
+ ( -12 * a + 108 * b + 144 * e - 9 ) * ujz +
+ ( 108 * b + 108 * c + 144 * e + 81 ) * ukz +
+ ( 108 * c - 12 * d + 144 * e - 9 ) * ulz;
+ M.coeffRef( 2, 3 ) = ( -8 * a + 12 * d + 96 * e - 1 ) * uiz +
+ ( -8 * a + 12 * b + 96 * e - 1 ) * ujz +
+ ( 12 * b + 72 * c + 96 * e + 9 ) * ukz +
+ ( 72 * c + 12 * d + 96 * e + 9 ) * ulz;
+ M.coeffRef( 3, 0 ) = M( 0, 3 );
+ M.coeffRef( 3, 1 ) = M( 1, 3 );
+ M.coeffRef( 3, 2 ) = M( 2, 3 );
+ M.coeffRef( 3, 3 ) = ( -12 * a + 108 * d + 144 * e - 9 ) * uiz +
+ ( -12 * a - 12 * b + 144 * e + 1 ) * ujz +
+ ( -12 * b + 108 * c + 144 * e - 9 ) * ukz +
+ ( 108 * c + 108 * d + 144 * e + 81 ) * ulz;
+ M *= factor;
+ }
+ return M;
+ }
+
+ Scalar quadrangleArea( Face f ) const
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ ASSERT( vtcs.size() == 4 );
+ const auto i = vtcs[ 0 ];
+ const auto j = vtcs[ 1 ];
+ const auto k = vtcs[ 2 ];
+ const auto l = vtcs[ 3 ];
+ const auto xij = myMesh->position( j ) - myMesh->position( i );
+ const auto xil = myMesh->position( l ) - myMesh->position( i );
+ const auto xlk = myMesh->position( k ) - myMesh->position( l );
+ const auto xjk = myMesh->position( k ) - myMesh->position( j );
+ const RealVector ni = xij.crossProduct( xil );
+ const RealVector nj = xij.crossProduct( xjk );
+ const RealVector nk = xlk.crossProduct( xjk );
+ const RealVector nl = xlk.crossProduct( xil );
+ const auto & ui = myMesh->vertexNormal( i );
+ const auto & uj = myMesh->vertexNormal( j );
+ const auto & uk = myMesh->vertexNormal( k );
+ const auto & ul = myMesh->vertexNormal( l );
+ const Scalar factor = 1.0 / 36.0;
+ if ( !use2ndOrder )
+ return factor * ( ni.dot( 4 * ui + 2 * uj + uk + 2 * ul ) +
+ nj.dot( 2 * ui + 4 * uj + 2 * uk + ul ) +
+ nk.dot( ui + 2 * uj + 4 * uk + 2 * ul ) +
+ nl.dot( 2 * ui + uj + 2 * uk + 4 * ul ) );
+ else
+ {
+ // Works only for surfels ! i.e. xij==xlk and xjk==xil
+ // If we assume f00=f10=f01=f11=1, the corrected area of a surfel is,
+ // after *36 factor A = (4*a+4*d+16*e+1) < u00 | n >
+ // + (4*a+4*b+16*e+1) < u10 | n >
+ // + (4*b+4*c+16*e+1) < u11 | n >
+ // + (4*c+4*d+16*e+1) < u01 | n >
+ const auto a = 1.0 / ( ui + uj ).norm();
+ const auto b = 1.0 / ( uj + uk ).norm();
+ const auto c = 1.0 / ( uk + ul ).norm();
+ const auto d = 1.0 / ( ul + ui ).norm();
+ const auto e = 1.0 / ( ui + uj + uk + ul ).norm();
+ return factor * ( ( 4 * a + 4 * d + 16 * e + 1 ) * ni.dot( ui ) +
+ ( 4 * a + 4 * b + 16 * e + 1 ) * nj.dot( uj ) +
+ ( 4 * b + 4 * c + 16 * e + 1 ) * nk.dot( uk ) +
+ ( 4 * c + 4 * d + 16 * e + 1 ) * nl.dot( ul ) );
+ }
+ }
+
+ public:
+ /// @name Initialization services
+ /// @{
+
+ /// Default constructor. The object is invalid.
+ InterpolatedCorrectedCalculus() : myMesh( nullptr )
+ {
+ }
+
+ /// Constructor from surface mesh \a smesh.
+ /// @param smesh any surface mesh
+ /// @param use2ndOrder wether to use 2nd order method for interpolating
+ /// normals (slightly more precise)
+ /// @param lambda value used for M1 positive definitness following @cite
+ /// degoes2020discrete
+ InterpolatedCorrectedCalculus( const ConstAlias smesh,
+ bool use2ndOrder = false,
+ double lambda = 0.1 )
+ : myMesh( &smesh ), use2ndOrder( use2ndOrder ), lambda( lambda )
+ {
+ if ( myMesh->vertexNormals().empty() )
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::init]"
+ << " vertex normals should be provided." << std::endl;
+ return;
+ }
+ }
+ ///@}
+
+ private:
+ DenseMatrix buildLocalM0( Index f ) const
+ {
+ if ( myMesh->incidentVertices( f ).size() == 4 )
+ return quadrangleInnerProduct0( f );
+ else
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildLocalM0]"
+ << " face should be quadrangles" << std::endl;
+ return DenseMatrix();
+ }
+ }
+
+ LinearOperator buildM0() const
+ {
+ Triplets triplets;
+ LinearOperator myM0 =
+ LinearOperator( myMesh->nbVertices(), myMesh->nbVertices() );
+ DenseMatrix values;
+ for ( Face f = 0; f < myMesh->nbFaces(); ++f )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto nv = vtcs.size();
+ if ( nv == 4 )
+ values = quadrangleInnerProduct0( f );
+ else
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::init]"
+ << " faces should be quadrangles" << std::endl;
+ return myM0;
+ }
+ for ( int i = 0; i < values.rows(); i++ )
+ for ( int j = 0; j < values.cols(); j++ )
+ triplets.push_back( { (StorageIndex)vtcs[ i ],
+ (StorageIndex)vtcs[ j ], values( i, j ) } );
+ }
+ myM0.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return myM0;
+ }
+
+ LinearOperator buildLumpedM0() const
+ {
+ return buildM0();
+ }
+
+ DenseMatrix buildLocalM1( Index f ) const
+ {
+ const auto a_f = quadrangleArea( f );
+ DenseMatrix CA = DenseMatrix::Identity( 3, 3 ) * a_f;
+ DenseMatrix sharp = buildLocalSharp( f );
+ DenseMatrix Id = DenseMatrix::Identity( 4, 4 );
+ DenseMatrix P = Id - buildLocalFlat( f ) * sharp;
+ // Compute inner product for 1-forms
+ DenseMatrix myM1 =
+ sharp.transpose() * CA * sharp + lambda * P.transpose() * P;
+ return myM1;
+ }
+
+ LinearOperator buildM1() const
+ {
+
+ LinearOperator CA( 3 * myMesh->nbFaces(), 3 * myMesh->nbFaces() );
+ Triplets triplets;
+ for ( Face f = 0; f < myMesh->nbFaces(); ++f )
+ {
+ const auto a_f = quadrangleArea( f );
+ if ( a_f <= 1e-8 )
+ trace.warning() << "Bad corrected area " << f << " " << a_f
+ << std::endl;
+ triplets.push_back( { (StorageIndex)3 * f, (StorageIndex)3 * f, a_f } );
+ triplets.push_back(
+ { (StorageIndex)3 * f + 1, (StorageIndex)3 * f + 1, a_f } );
+ triplets.push_back(
+ { (StorageIndex)3 * f + 2, (StorageIndex)3 * f + 2, a_f } );
+ }
+ CA.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ // Compute projection operator
+ LinearOperator mySharp = buildSharp();
+ LinearOperator Id1( myMesh->nbEdges(), myMesh->nbEdges() );
+ Id1.setIdentity();
+ LinearOperator P = Id1 - buildFlat() * mySharp;
+ // Compute inner product for 1-forms
+ LinearOperator myM1 =
+ mySharp.transpose() * CA * mySharp + lambda * P.transpose() * P;
+ return myM1;
+ }
+
+ DenseMatrix buildLocalM2( Index f ) const
+ {
+ const auto nv = myMesh->incidentVertices( f ).size();
+ DenseMatrix res;
+ if ( nv == 4 )
+ res = DenseMatrix::Identity( 1, 1 ) / quadrangleArea( f );
+ else
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildLocalM2]"
+ << " face should be quadrangles" << std::endl;
+ }
+ return res;
+ }
+
+ LinearOperator buildM2() const
+ {
+ Triplets triplets;
+ LinearOperator myM2 =
+ LinearOperator( myMesh->nbFaces(), myMesh->nbFaces() );
+ for ( Face f = 0; f < myMesh->nbFaces(); ++f )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ double area = 0.;
+ if ( vtcs.size() == 4 )
+ area = quadrangleArea( f );
+ else
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildM2]"
+ << " faces should be quadrangular." << std::endl;
+ }
+ triplets.push_back( { (StorageIndex)f, (StorageIndex)f, 1.0 / area } );
+ }
+ myM2.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return myM2;
+ }
+
+ DenseMatrix buildLocalSharp( Index f ) const
+ {
+ if ( myMesh->incidentVertices( f ).size() != 4 )
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildLocalSharp]"
+ << " restricted to quadrangle faces." << std::endl;
+ return DenseMatrix();
+ }
+ return quadrangleSharpLocalMatrix( f );
+ }
+
+ LinearOperator buildSharp() const
+ {
+ Triplets triplets;
+ LinearOperator mySharp =
+ LinearOperator( 3 * myMesh->nbFaces(), myMesh->nbEdges() );
+ for ( Face f = 0; f < myMesh->nbFaces(); ++f )
+ {
+ // std::cout << "Face " << f << std::endl;
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto nc = vtcs.size();
+ if ( nc != 4 )
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildSharp]"
+ << " restricted to quadrangle faces." << std::endl;
+ return mySharp;
+ }
+ Edge edge_idcs[ 4 ];
+ Scalar edge_sign[ 4 ];
+ for ( Size i = 0; i < nc; i++ )
+ {
+ edge_idcs[ i ] =
+ myMesh->makeEdge( vtcs[ i ], vtcs[ ( i + 1 ) % nc ] );
+ edge_sign[ i ] = vtcs[ i ] < vtcs[ ( i + 1 ) % nc ] ? 1.0 : -1.0;
+ }
+ const auto G_f = buildLocalSharp( f );
+ for ( Dimension x = 0; x < 3; ++x )
+ for ( Dimension n = 0; n < nc; ++n )
+ triplets.push_back( { (StorageIndex)3 * f + x,
+ (StorageIndex)edge_idcs[ n ],
+ edge_sign[ n ] * G_f( x, n ) } );
+ }
+ mySharp.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return mySharp;
+ }
+
+ DenseMatrix buildLocalFlat( Index f ) const
+ {
+ if ( myMesh->incidentVertices( f ).size() != 4 )
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildLocalFlat]"
+ << " restricted to quadrangle faces." << std::endl;
+ return DenseMatrix();
+ }
+ return quadrangleFlatLocalMatrix( f );
+ }
+
+ LinearOperator buildFlat() const
+ {
+ Triplets triplets;
+ LinearOperator myFlat =
+ LinearOperator( myMesh->nbEdges(), 3 * myMesh->nbFaces() );
+ for ( Face f = 0; f < myMesh->nbFaces(); ++f )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto nc = vtcs.size();
+ if ( nc != 4 )
+ {
+ trace.error() << "[InterpolatedCorrectedCalculus::buildFlat]"
+ << " restricted to quadrangle faces." << std::endl;
+ return myFlat;
+ }
+ Edge edge_idcs[ 4 ];
+ Scalar edge_sign[ 4 ];
+ for ( Size i = 0; i < nc; i++ )
+ {
+ edge_idcs[ i ] =
+ myMesh->makeEdge( vtcs[ i ], vtcs[ ( i + 1 ) % nc ] );
+ edge_sign[ i ] = vtcs[ i ] < vtcs[ ( i + 1 ) % nc ] ? 1.0 : -1.0;
+ }
+ const auto V_f = buildLocalFlat( f );
+ for ( Dimension x = 0; x < 3; ++x )
+ for ( Dimension n = 0; n < nc; ++n )
+ {
+ const double nbFaces = myMesh->edgeFaces( edge_idcs[ n ] ).size();
+ triplets.push_back( { (StorageIndex)edge_idcs[ n ],
+ (StorageIndex)3 * f + x,
+ edge_sign[ n ] * V_f( n, x ) / nbFaces } );
+ }
+ }
+ myFlat.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return myFlat;
+ }
+
+ DenseMatrix buildLocalD0( Index f ) const
+ {
+ const auto nc = myMesh->incidentVertices( f ).size();
+ DenseMatrix res = DenseMatrix::Zero( nc, nc );
+ for ( int i = 0; i < nc; i++ )
+ {
+ res( ( i - 1 + nc ) % nc, i ) = 1;
+ res( i, i ) = -1;
+ }
+ return res;
+ }
+
+ LinearOperator buildD0() const
+ {
+ Triplets triplets;
+ // trace.beginBlock( "Init derivative operator D0" );
+ Edge e = 0;
+ for ( auto && vtcs : myMesh->allEdgeVertices() )
+ {
+ triplets.push_back( { (StorageIndex)e, (StorageIndex)vtcs.first, -1 } );
+ triplets.push_back( { (StorageIndex)e, (StorageIndex)vtcs.second, 1 } );
+ e++;
+ }
+ LinearOperator myD0 =
+ LinearOperator( myMesh->nbEdges(), myMesh->nbVertices() );
+ myD0.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return myD0;
+ }
+
+ DenseMatrix buildLocalL0( Index f ) const
+ {
+ DenseMatrix D0 = buildLocalD0( f );
+ DenseMatrix M1 = buildLocalM1( f );
+ return -1.0 * D0.transpose() * M1 * D0;
+ }
+
+ LinearOperator buildL0() const
+ {
+ LinearOperator myD0 = buildD0();
+ LinearOperator myM1 = buildM1();
+ return -1.0 * myD0.transpose() * myM1 * myD0;
+ }
+ };
+
+} // namespace DGtal
+
+#endif
diff --git a/src/DGtal/dec/NormalCorrectedFEM.h b/src/DGtal/dec/NormalCorrectedFEM.h
new file mode 100644
index 0000000000..37e3a3c892
--- /dev/null
+++ b/src/DGtal/dec/NormalCorrectedFEM.h
@@ -0,0 +1,467 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+ #pragma once
+
+
+/**
+ * @file
+ * @author
+ *
+ * @date 2024/06/21
+ *
+ * Header file for module SurfaceDEC.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if !defined NormalCorrectedFEM_h
+/** Prevents repeated inclusion of headers. */
+#define NormalCorrectedFEM_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include "DGtal/base/Common.h"
+#include "DGtal/helpers/StdDefs.h"
+#include "DGtal/shapes/SurfaceMesh.h"
+#include
+#include
+#include
+
+namespace DGtal
+{
+ /////////////////////////////////////////////////////////////////////////////
+ // template class CorrectedFEM
+ /**
+ Description of template class 'NormalCorrectedFEM' \brief Aim:
+ Provides methods for building a stiffness and mass matrix to solve
+ Poisson problems using the Finite Element Method using a
+ corrected normal field.
+
+ Works with triangular and quadrangular surfaces.
+
+ Requires face normals on attached surface mesh.
+
+ * @tparam TLinearAlgebraBackend linear algebra backend used (i.e.
+ EigenSparseLinearAlgebraBackend).
+ @tparam TRealPoint an arbitrary model of 3D RealPoint.
+ @tparam TRealVector an arbitrary model of 3D RealVector.
+ */
+ template
+ struct NormalCorrectedFEM
+ {
+ typedef TLinearAlgebraBackend LinearAlgebraBackend;
+ typedef TRealPoint RealPoint;
+ typedef TRealVector RealVector;
+ typedef typename LinearAlgebraBackend::SparseMatrix LinearOperator;
+ typedef typename LinearAlgebraBackend::SparseMatrix::StorageIndex StorageIndex;
+ typedef typename LinearAlgebraBackend::DenseMatrix DenseMatrix;
+ typedef typename LinearAlgebraBackend::DenseVector DenseVector;
+ typedef SurfaceMesh Mesh;
+ typedef typename Mesh::Vertex Vertex;
+ typedef typename Mesh::Edge Edge;
+ typedef typename Mesh::Face Face;
+ typedef typename Mesh::Face Corner;
+ typedef typename Mesh::Size Size;
+
+ /// A triplet (row, col, value), useful for initializaing sparse matrices.
+ typedef typename LinearAlgebraBackend::Triplet Triplet;
+ /// A range of triplets (row,col,value).
+ typedef std::vector Triplets;
+
+ private:
+ // ------------------------- Initialization services
+ // ----------------------------
+
+ void buildTriangleFaceMass( Triplets & triplets, Triplets & lumpedTriplets,
+ Face f ) const
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto v0 = vtcs[ 0 ];
+ const auto v1 = vtcs[ 1 ];
+ const auto v2 = vtcs[ 2 ];
+ const auto xij = myMesh->position( v1 ) - myMesh->position( v0 );
+ const auto xik = myMesh->position( v2 ) - myMesh->position( v0 );
+ const auto n = xij.crossProduct( xik );
+ const auto nn = n.getNormalized();
+ const auto face_n = myMesh->faceNormal( f );
+ const double z = face_n.dot( nn );
+ const double doubleArea = n.norm();
+ triplets.push_back( { (StorageIndex)v0, (StorageIndex)v0, 1. / 12. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v0, (StorageIndex)v1, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v0, (StorageIndex)v2, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v1, (StorageIndex)v0, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v1, (StorageIndex)v1, 1. / 12. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v1, (StorageIndex)v2, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v2, (StorageIndex)v0, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v2, (StorageIndex)v1, 1. / 24. * z * doubleArea } );
+ triplets.push_back( { (StorageIndex)v2, (StorageIndex)v2, 1. / 12. * z * doubleArea } );
+ lumpedTriplets.push_back( { (StorageIndex)v0, (StorageIndex)v0, 1. / 6. * z * doubleArea } );
+ lumpedTriplets.push_back( { (StorageIndex)v1, (StorageIndex)v1, 1. / 6. * z * doubleArea } );
+ lumpedTriplets.push_back( { (StorageIndex)v2, (StorageIndex)v2, 1. / 6. * z * doubleArea } );
+ }
+
+ void buildQuadrangleFaceMass( Triplets & triplets,
+ Triplets & lumpedTriplets, Face f ) const
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto v0 = vtcs[ 0 ];
+ const auto v1 = vtcs[ 1 ];
+ const auto v2 = vtcs[ 2 ];
+ const auto v3 = vtcs[ 3 ];
+ const auto xij = myMesh->position( v1 ) - myMesh->position( v0 );
+ const auto xil = myMesh->position( v3 ) - myMesh->position( v0 );
+ const auto nn = xij.crossProduct( xil ).getNormalized();
+ const auto face_n = myMesh->faceNormal( f );
+ const double z = nn.dot( face_n );
+ if ( z < 0. )
+ {
+ trace.warning() << "Bad corrected normals at face " << f << std::endl;
+ }
+
+ double area = xij.crossProduct( xil ).norm();
+ int v[ 2 ][ 2 ] = { { (int)v0, (int)v1 }, { (int)v3, (int)v2 } };
+ for ( int i = 0; i < 2; i++ )
+ {
+ for ( int j = 0; j < 2; j++ )
+ {
+ for ( int i2 = 0; i2 < 2; i2++ )
+ {
+ for ( int j2 = 0; j2 < 2; j2++ )
+ {
+ int v_i = v[ i ][ j ];
+ int v_j = v[ i2 ][ j2 ];
+ if ( i == i2 && j == j2 )
+ {
+ triplets.push_back( { (StorageIndex)v_i, (StorageIndex)v_j, area * z * 1. / 9. } );
+ lumpedTriplets.push_back( { (StorageIndex)v_i, (StorageIndex)v_j, 1. / 4. * z * area } );
+ }
+ else if ( i == i2 && j == 1 - j2 )
+ {
+ triplets.push_back( { (StorageIndex)v_i, (StorageIndex)v_j, area * z * 1. / 18. } );
+ }
+ else if ( j == j2 && i == 1 - i2 )
+ {
+ triplets.push_back( { (StorageIndex)v_i, (StorageIndex)v_j, area * z * 1. / 18. } );
+ }
+ else
+ {
+ triplets.push_back( { (StorageIndex)v_i, (StorageIndex)v_j, area * z * 1. / 36. } );
+ }
+ }
+ }
+ }
+ }
+ }
+
+ DenseMatrix buildTriangleFaceStiffness( Face f ) const
+ {
+ DenseMatrix res = DenseMatrix( 3, 3 );
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto v0 = vtcs[ 0 ];
+ const auto v1 = vtcs[ 1 ];
+ const auto v2 = vtcs[ 2 ];
+ const auto xij = myMesh->position( v1 ) - myMesh->position( v0 );
+ const auto xik = myMesh->position( v2 ) - myMesh->position( v0 );
+ const auto n = xij.crossProduct( xik );
+ const auto nn = n.getNormalized();
+ const auto xil = nn.crossProduct( xij ).getNormalized();
+ const auto face_n = myMesh->faceNormal( f );
+ const double x = face_n.dot( xij.getNormalized() );
+ const double y = face_n.dot( xil );
+ const double z = face_n.dot( nn );
+ const double l1 = xij.norm();
+ const double l2 = xik.norm();
+ const double cos_t = xij.dot( xik ) / ( l1 * l2 );
+ const double sin_t = n.norm() / ( l1 * l2 );
+ const double x2 = x * x;
+ const double y2 = y * y;
+ const double z2 = z * z;
+ const double l12 = l1 * l1;
+ const double l22 = l2 * l2;
+ const double cos_t2 = cos_t * cos_t;
+ const double sin_t2 = sin_t * sin_t;
+ res( 0, 0 ) =
+ 1. / 2 *
+ ( 2 * ( l1 * l2 - cos_t * l22 ) * sin_t * x * y - 2 * cos_t * l1 * l2 +
+ l22 * sin_t2 * z2 - l22 * sin_t2 +
+ ( 2 * cos_t * l1 * l2 + 2 * l22 * sin_t2 - l12 - l22 ) * x2 + l12 +
+ l22 ) /
+ ( l1 * l2 * sin_t * z );
+ res( 1, 0 ) =
+ 1. / 2 *
+ ( ( 2 * cos_t * l2 - l1 ) * sin_t * x * y + l2 * sin_t2 * y2 +
+ cos_t * l1 - ( cos_t * l1 + l2 * sin_t2 - l2 ) * x2 - l2 ) /
+ ( l1 * sin_t * z );
+ res( 2, 0 ) =
+ -1. / 2 *
+ ( l2 * sin_t * x * y - cos_t * l2 + ( cos_t * l2 - l1 ) * x2 + l1 ) /
+ ( l2 * sin_t * z );
+ res( 0, 1 ) =
+ 1. / 2 *
+ ( ( 2 * cos_t * l2 - l1 ) * sin_t * x * y + l2 * sin_t2 * y2 +
+ cos_t * l1 - ( cos_t * l1 + l2 * sin_t2 - l2 ) * x2 - l2 ) /
+ ( l1 * sin_t * z );
+ res( 1, 1 ) = -1. / 2 *
+ ( 2 * cos_t * l2 * sin_t * x * y + cos_t2 * l2 * x2 -
+ ( cos_t2 - 1 ) * l2 * y2 - l2 ) /
+ ( l1 * sin_t * z );
+ res( 2, 1 ) =
+ 1. / 2 * ( sin_t * x * y + cos_t * x2 - cos_t ) / ( sin_t * z );
+ res( 0, 2 ) =
+ -1. / 2 *
+ ( l2 * sin_t * x * y - cos_t * l2 + ( cos_t * l2 - l1 ) * x2 + l1 ) /
+ ( l2 * sin_t * z );
+ res( 1, 2 ) =
+ 1. / 2 * ( sin_t * x * y + cos_t * x2 - cos_t ) / ( sin_t * z );
+ res( 2, 2 ) = -1. / 2 * ( l1 * x2 - l1 ) / ( l2 * sin_t * z );
+ return res;
+ }
+
+ DenseMatrix buildQuadrangleFaceStiffness( Face f ) const
+ {
+ DenseMatrix res = DenseMatrix( 4, 4 );
+ const auto vtcs = myMesh->incidentVertices( f );
+ const auto v0 = vtcs[ 0 ];
+ const auto v1 = vtcs[ 1 ];
+ const auto v2 = vtcs[ 2 ];
+ const auto v3 = vtcs[ 3 ];
+ const auto xij = myMesh->position( v1 ) - myMesh->position( v0 );
+ const auto xil = myMesh->position( v3 ) - myMesh->position( v0 );
+ const auto nn = xij.crossProduct( xil ).getNormalized();
+ const auto face_n = myMesh->faceNormal( f );
+ const double x = face_n.dot( xij.getNormalized() );
+ const double y = face_n.dot( xil.getNormalized() );
+ const double z = face_n.dot( nn );
+ if ( z < 0. )
+ {
+ trace.warning() << "Bad corrected normals at face " << f << std::endl;
+ }
+ const double lap = 0.5 * ( x * y / z );
+ const double sqx = 1. - x * x;
+ const double sqy = 1. - y * y;
+ const double self = lap + 1. / 3. * ( sqx + sqy ) / z;
+ const double op1 = -1. / 3. * sqy / z + 1. / 6. * sqx / z;
+ const double op2 = -1. / 3. * sqx / z + 1. / 6. * sqy / z;
+ const double op = -lap - 1. / 6. * ( sqx + sqy ) / z;
+ res( 0, 0 ) = self;
+ res( 0, 1 ) = op1;
+ res( 0, 2 ) = op;
+ res( 0, 3 ) = op2;
+
+ res( 1, 0 ) = op1;
+ res( 1, 1 ) = self - 2. * lap;
+ res( 1, 2 ) = op2;
+ res( 1, 3 ) = op + 2. * lap;
+
+ res( 2, 0 ) = op;
+ res( 2, 1 ) = op2;
+ res( 2, 2 ) = self;
+ res( 2, 3 ) = op1;
+
+ res( 3, 0 ) = op2;
+ res( 3, 1 ) = op + 2. * lap;
+ res( 3, 2 ) = op1;
+ res( 3, 3 ) = self - 2. * lap;
+ return res;
+ }
+
+ LinearOperator buildMass()
+ {
+ return true;
+ }
+
+ bool init( const Mesh * pMesh )
+ {
+ myMesh = pMesh;
+ if ( order != 1 )
+ {
+ trace.error() << "[NormalCorrectedFEM::init]"
+ << " order should be 1" << std::endl;
+ return false;
+ }
+ if ( myMesh->faceNormals().empty() )
+ {
+ trace.error() << "[NormalCorrectedFEM::init]"
+ << " face normals should be provided." << std::endl;
+ return false;
+ }
+ return true;
+ }
+
+ LinearOperator buildL0() const
+ {
+ LinearOperator L( myMesh->nbVertices(), myMesh->nbVertices() );
+ Triplets triplets;
+ DenseMatrix values;
+ for ( Face f = 0; f < myMesh->nbFaces(); f++ )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ if ( vtcs.size() == 3 )
+ values = buildTriangleFaceStiffness( f );
+ else if ( vtcs.size() == 4 )
+ values = buildQuadrangleFaceStiffness( f );
+ else
+ {
+ trace.error() << "[NormalCorrectedFEM::init]"
+ << " faces should be triangle or quadrangles"
+ << std::endl;
+
+ return L;
+ }
+ for ( int i = 0; i < values.rows(); i++ )
+ for ( int j = 0; j < values.cols(); j++ )
+ triplets.push_back( { (StorageIndex)vtcs[ i ], (StorageIndex)vtcs[ j ], values( i, j ) } );
+ }
+ L.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return -L;
+ }
+
+ LinearOperator buildM0() const
+ {
+ LinearOperator M( myMesh->nbVertices(), myMesh->nbVertices() );
+
+ Triplets triplets;
+ Triplets lumpedTriplets;
+ for ( Face f = 0; f < myMesh->nbFaces(); f++ )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ if ( vtcs.size() == 3 )
+ buildTriangleFaceMass( triplets, lumpedTriplets, f );
+ else if ( vtcs.size() == 4 )
+ buildQuadrangleFaceMass( triplets, lumpedTriplets, f );
+ else
+ {
+ trace.error() << "[NormalCorrectedFEM::buildM0]"
+ << " faces should be triangles or quadrangles"
+ << std::endl;
+ return M;
+ }
+ }
+ M.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return M;
+ }
+
+ LinearOperator buildLumpedM0() const
+ {
+ LinearOperator lumpedM( myMesh->nbVertices(), myMesh->nbVertices() );
+
+ Triplets triplets;
+ Triplets lumpedTriplets;
+ for ( Face f = 0; f < myMesh->nbFaces(); f++ )
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ if ( vtcs.size() == 3 )
+ buildTriangleFaceMass( triplets, lumpedTriplets, f );
+ else if ( vtcs.size() == 4 )
+ buildQuadrangleFaceMass( triplets, lumpedTriplets, f );
+ else
+ {
+ trace.error() << "[NormalCorrectedFEM::buildLumpedM0]"
+ << " faces should be triangles or quadrangles"
+ << std::endl;
+ return lumpedM;
+ }
+ }
+ lumpedM.setFromTriplets( lumpedTriplets.cbegin(), lumpedTriplets.cend() );
+ return lumpedM;
+ }
+
+ private:
+ const Mesh * myMesh;
+ int order;
+
+ public:
+ /// @name Initialization services
+ /// @{
+
+ /// Default constructor. The object is invalid.
+ NormalCorrectedFEM() : myMesh( nullptr )
+ {
+ }
+
+ /// Constructor from surface mesh \a smesh.
+ /// @param smesh any surface mesh
+ /// @param order the order (default 1)
+ NormalCorrectedFEM( ConstAlias smesh, int order = 1 )
+ : myMesh( nullptr ), order( order )
+ {
+ init( &smesh );
+ }
+
+ /// @}
+ /// @name Standard services
+ /// @{
+
+ /// Return the stiffness matrix degree x degree of the face.
+ /// @param f a face
+ /// @return the degree x degree stiffness matrix
+ DenseMatrix localL0( Face f ) const
+ {
+ const auto vtcs = myMesh->incidentVertices( f );
+ if ( vtcs.size() == 3 )
+ return -buildTriangleFaceStiffness( f );
+ else if ( vtcs.size() == 4 )
+ return -buildQuadrangleFaceStiffness( f );
+ else
+ {
+ trace.error() << "[NormalCorrectedFEM::buildLocalL0]"
+ << " faces should be triangles or quadrangles"
+ << std::endl;
+ return DenseMatrix();
+ }
+ }
+
+ /// Return the global stiffness matrix n_v x n_v.
+ /// @return the n_v x n_v stiffness matrix
+ ///
+ /// @note The sign convention for the divergence and the Laplacian
+ /// operator is opposite to the one of @cite degoes2020discrete .
+ /// 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ LinearOperator L0() const
+ {
+ return buildL0();
+ }
+
+ /// Return the global mass matrix n_v x n_v.
+ /// This matrix is not diagonal. To obtain a (less precise)
+ /// diagonal version use lumpedM0()
+ /// @return the n_v x n_v mass matrix
+ LinearOperator M0() const
+ {
+ return buildM0();
+ }
+
+ /// Return the global lumped mass matrix n_v x n_v.
+ /// This matrix is diagonal. To obtain a (more precise)
+ /// non diagonal version use M0()
+ /// @return the n_v x n_v mass matrix
+ LinearOperator lumpedM0() const
+ {
+ return buildLumpedM0();
+ }
+ /// @}
+ };
+} // namespace DGtal
+
+#endif
diff --git a/src/DGtal/dec/PolygonalCalculus.h b/src/DGtal/dec/PolygonalCalculus.h
index 116a29d2ae..da1fb8da86 100644
--- a/src/DGtal/dec/PolygonalCalculus.h
+++ b/src/DGtal/dec/PolygonalCalculus.h
@@ -19,7 +19,8 @@
/**
* @file
* @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
- * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS,
+ * UMR 5205), CNRS, France
*
* @date 2021/09/02
*
@@ -38,1273 +39,1522 @@
#include "DGtal/base/ConstAlias.h"
#include "DGtal/base/Common.h"
#include "DGtal/shapes/SurfaceMesh.h"
+#include "DGtal/dec/SurfaceDEC.h"
#include "DGtal/math/linalg/EigenSupport.h"
//////////////////////////////////////////////////////////////////////////////
namespace DGtal
{
- namespace functors {
+ namespace functors
+ {
/**
*
- * \brief Functor that projects a face vertex of a surface mesh onto the tangent plane
- * given by a per-face normal vector.
- * This functor can be used in PolygonalCalculus to correct the embedding of
- * digital surfaces using an estimated normal vector field (see @cite coeurjolly2022simple).
+ * \brief Functor that projects a face vertex of a surface mesh onto the
+ * tangent plane given by a per-face normal vector. This functor can be used
+ * in PolygonalCalculus to correct the embedding of digital surfaces using
+ * an estimated normal vector field (see @cite coeurjolly2022simple).
*
- * @note when used in PolygonalCalculus, all operators being invariant by translation, all
- * tangent planes pass through the origin (0,0,0) (no offest).
+ * @note when used in PolygonalCalculus, all operators being invariant by
+ * translation, all tangent planes pass through the origin (0,0,0) (no
+ * offest).
*
- * @tparam TRealPoint a model of points @f$\mathbb{R}^3@f$ (e.g. PointVector).
- * @tparam TRealVector a model of vectors in @f$\mathbb{R}^3@f$ (e.g. PointVector).
- */ template
+ * @tparam TRealPoint a model of points @f$\mathbb{R}^3@f$ (e.g.
+ * PointVector).
+ * @tparam TRealVector a model of vectors in @f$\mathbb{R}^3@f$ (e.g.
+ * PointVector).
+ */
+ template
struct EmbedderFromNormalVectors
{
- ///Type of SurfaceMesh
+ /// Type of SurfaceMesh
typedef SurfaceMesh MySurfaceMesh;
- ///Vertex type
+ /// Vertex type
typedef typename MySurfaceMesh::Vertex Vertex;
- ///Face type
+ /// Face type
typedef typename MySurfaceMesh::Face Face;
- ///Position type
+ /// Position type
typedef typename MySurfaceMesh::RealPoint Real3dPoint;
-
+
EmbedderFromNormalVectors() = delete;
-
- /// Constructor from an array of normal vectors and a surface mesh instance.
- /// @param normals a vector of per face normal vectors (same ordering as the SurfaceMesh face indicies).
+
+ /// Constructor from an array of normal vectors and a surface mesh
+ /// instance.
+ /// @param normals a vector of per face normal vectors (same ordering as
+ /// the SurfaceMesh face indicies).
/// @param surfmesh an instance of SurfaceMesh
- EmbedderFromNormalVectors(ConstAlias> normals,
- ConstAlias surfmesh)
+ EmbedderFromNormalVectors( ConstAlias> normals,
+ ConstAlias surfmesh )
{
myNormals = &normals;
mySurfaceMesh = &surfmesh;
}
-
- /// Project a face vertex onto its tangent plane (given by the per-face estimated
- /// normal vector).
+
+ /// Project a face vertex onto its tangent plane (given by the per-face
+ /// estimated normal vector).
///
/// @param f the face that contains the vertex
/// @param v the vertex to project
- Real3dPoint operator()(const Face &f,const Vertex &v)
+ Real3dPoint operator()( const Face & f, const Vertex & v )
{
- const auto nn = (*myNormals)[f];
- Real3dPoint p = mySurfaceMesh->position(v);
- return p - nn.dot(p)*nn;
+ const auto nn = ( *myNormals )[ f ];
+ Real3dPoint p = mySurfaceMesh->position( v );
+ return p - nn.dot( p ) * nn;
}
-
- ///Alias to the normal vectors
- const std::vector *myNormals;
- ///Alias to the surface mesh
- const MySurfaceMesh *mySurfaceMesh;
+
+ /// Alias to the normal vectors
+ const std::vector * myNormals;
+ /// Alias to the surface mesh
+ const MySurfaceMesh * mySurfaceMesh;
};
- }
-
-
-
-/////////////////////////////////////////////////////////////////////////////
-// template class PolygonalCalculus
-/**
- * Description of template class 'PolygonalCalculus'
- * \brief Implements differential operators on polygonal surfaces from
- * @cite degoes2020discrete
- *
- * See @ref modulePolygonalCalculus for details.
- *
- * @note The sign convention for the divergence and the Laplacian
- * operator is opposite to the one of @cite degoes2020discrete. 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 \f$ \langle \mathrm{d} u, v \rangle = -
- * \langle u, \mathrm{div} v \rangle \f$. See also
- * https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- *
- * @tparam TRealPoint a model of points @f$\mathbb{R}^3@f$ (e.g. PointVector).
- * @tparam TRealVector a model of vectors in @f$\mathbb{R}^3@f$ (e.g. PointVector).
- */
-template
-class PolygonalCalculus
-{
- // ----------------------- Standard services ------------------------------
-public:
-
- ///Concept checking
- static const Dimension dimension = TRealPoint::dimension;
- BOOST_STATIC_ASSERT( ( dimension == 3 ) );
-
- ///Self type
- typedef PolygonalCalculus Self;
-
- ///Type of SurfaceMesh
- typedef SurfaceMesh MySurfaceMesh;
- ///Vertex type
- typedef typename MySurfaceMesh::Vertex Vertex;
- ///Face type
- typedef typename MySurfaceMesh::Face Face;
- ///Position type
- typedef typename MySurfaceMesh::RealPoint Real3dPoint;
- ///Real vector type
- typedef typename MySurfaceMesh::RealVector Real3dVector;
-
- ///Linear Algebra Backend from Eigen
- typedef EigenLinearAlgebraBackend LinAlg;
- ///Type of Vector
- typedef LinAlg::DenseVector Vector;
- ///Global 0-form, 1-form, 2-form are Vector
- typedef Vector Form;
- ///Type of dense matrix
- typedef LinAlg::DenseMatrix DenseMatrix;
- ///Type of sparse matrix
- typedef LinAlg::SparseMatrix SparseMatrix;
- ///Type of sparse matrix triplet
- typedef LinAlg::Triplet Triplet;
-
- ///Type of a sparse matrix solver
- typedef LinAlg::SolverSimplicialLDLT Solver;
-
- /// @name Standard services
- /// @{
-
- /// Create a Polygonal DEC structure from a surface mesh (@a surf)
- /// using an default identity embedder.
- /// @param surf an instance of SurfaceMesh
- /// @param globalInternalCacheEnabled enable the internal cache for all operators (default: false)
- PolygonalCalculus(const ConstAlias 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();
- };
-
- /// Create a Polygonal DEC structure from a surface mesh (@a surf)
- /// and an embedder for the vertex position: function with two parameters, a face and a vertex
- /// which outputs the embedding in R^3 of the vertex w.r.t. to the face.
- /// @param surf an instance of SurfaceMesh
- /// @param embedder an embedder for the vertex position
- /// @param globalInternalCacheEnabled if true, the global operator cache is enabled
- PolygonalCalculus(const ConstAlias surf,
- const std::function &embedder,
- bool globalInternalCacheEnabled = false):
- mySurfaceMesh(&surf), myEmbedder(embedder), myGlobalCacheEnabled(globalInternalCacheEnabled)
- {
- myVertexNormalEmbedder = [&](Vertex v){ return toReal3dVector(computeVertexNormal(v)); };
- init();
- };
-
- /// Create a Polygonal DEC structure from a surface mesh (@a surf)
- /// and an embedder for the vertex normal: function with a vertex as
- /// parameter which outputs the embedding in R^3 of the vertex normal.
- /// @param surf an instance of SurfaceMesh
- /// @param embedder an embedder for the vertex position
- /// @param globalInternalCacheEnabled if true, the global operator cache is enabled
- PolygonalCalculus(const ConstAlias surf,
- const std::function &embedder,
- bool globalInternalCacheEnabled = false):
- mySurfaceMesh(&surf), myVertexNormalEmbedder(embedder),
- myGlobalCacheEnabled(globalInternalCacheEnabled)
- {
- myEmbedder = [&](Face f,Vertex v){(void)f; return mySurfaceMesh->position(v); };
- init();
- };
-
- /// Create a Polygonal DEC structure from a surface mesh (@a surf)
- /// and an embedder for the vertex position: function with two parameters, a
- /// face and a vertex which outputs the embedding in R^3 of the vertex
- /// w.r.t. to the face. and an embedder for the vertex normal: function with
- /// a vertex as parameter which outputs the embedding in R^3 of the vertex
- /// normal.
- /// @param surf an instance of SurfaceMesh
- /// @param pos_embedder an embedder for the position
- /// @param normal_embedder an embedder for the position
- /// @param globalInternalCacheEnabled
- PolygonalCalculus(const ConstAlias surf,
- const std::function &pos_embedder,
- const std::function &normal_embedder,
- bool globalInternalCacheEnabled = false) :
- mySurfaceMesh(&surf), myEmbedder(pos_embedder),
- myVertexNormalEmbedder(normal_embedder),
- myGlobalCacheEnabled(globalInternalCacheEnabled)
- {
- init();
- };
+ } // namespace functors
+ /////////////////////////////////////////////////////////////////////////////
+ // template class PolygonalCalculus
/**
- * Deleted default constructor.
+ * Description of template class 'PolygonalCalculus'
+ * \brief Implements differential operators on polygonal surfaces from
+ * @cite degoes2020discrete
+ *
+ * See @ref modulePolygonalCalculus for details.
+ *
+ * @note The sign convention for the divergence and the Laplacian
+ * operator is opposite to the one of @cite degoes2020discrete. 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 \f$ \langle \mathrm{d} u, v \rangle = -
+ * \langle u, \mathrm{div} v \rangle \f$. See also
+ * https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ *
+ * @tparam TRealPoint a model of points @f$\mathbb{R}^3@f$ (e.g. PointVector).
+ * @tparam TRealVector a model of vectors in @f$\mathbb{R}^3@f$ (e.g.
+ * PointVector).
*/
- PolygonalCalculus() = delete;
-
- /**
- * Destructor (default).
- */
- ~PolygonalCalculus() = default;
-
- /**
- * Deleted copy constructor.
- * @param other the object to clone.
- */
- PolygonalCalculus ( const PolygonalCalculus & other ) = delete;
-
- /**
- * Deleted move constructor.
- * @param other the object to move.
- */
- PolygonalCalculus ( PolygonalCalculus && other ) = delete;
-
- /**
- * Deleted copy assignment operator.
- * @param other the object to copy.
- * @return a reference on 'this'.
- */
- PolygonalCalculus & operator= ( const PolygonalCalculus & other ) = delete;
-
- /**
- * Deleted move assignment operator.
- * @param other the object to move.
- * @return a reference on 'this'.
- */
- PolygonalCalculus & operator= ( PolygonalCalculus && other ) = delete;
-
- /// @}
-
- // ----------------------- embedding services --------------------------
- //MARK: Embedding services
- /// @name Embedding services
- /// @{
-
- /// Update the embedding function.
- /// @param externalFunctor a new embedding functor (Face,Vertex)->RealPoint.
- void setEmbedder(const std::function &externalFunctor)
- {
- myEmbedder = externalFunctor;
- }
- /// @}
-
- // ----------------------- Per face operators --------------------------------------
- //MARK: Per face operator on scalars
- /// @name Per face operators on scalars
- /// @{
-
- /// Return the vertex position matrix degree x 3 of the face.
- /// @param f a face
- /// @return the n_f x 3 position matrix
- DenseMatrix X(const Face f) const
+ template
+ class PolygonalCalculus
+ : public SurfaceDEC,
+ EigenLinearAlgebraBackend>
{
- 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;
- }
-
+ friend class SurfaceDEC;
- /// Derivative operator (d_0) of a face.
- /// @param f the face
- /// @return a degree x degree matrix
- 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;
- }
-
- /// Edge vector operator per face.
- /// @param f the face
- /// @return degree x 3 matrix
- DenseMatrix E(const Face f) const
- {
- if (checkCache(E_,f))
- return myGlobalCache[E_][f];
+ // ----------------------- Standard services ------------------------------
+ public:
+ /// Concept checking
+ static const Dimension dimension = TRealPoint::dimension;
+ BOOST_STATIC_ASSERT( ( dimension == 3 ) );
- DenseMatrix op = D(f)*X(f);
-
- setInCache(E_,f,op);
- return op;
- }
-
- /// Average operator to average, per edge, its vertex values.
- /// @param f the face
- /// @return a degree x degree matrix
- 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)
+ /// Self type
+ typedef PolygonalCalculus Self;
+
+ /// Type of SurfaceMesh
+ typedef SurfaceMesh MySurfaceMesh;
+ /// Vertex type
+ typedef typename MySurfaceMesh::Vertex Vertex;
+ /// Face type
+ typedef typename MySurfaceMesh::Face Face;
+ /// Position type
+ typedef typename MySurfaceMesh::RealPoint Real3dPoint;
+ /// Real vector type
+ typedef typename MySurfaceMesh::RealVector Real3dVector;
+
+ /// Linear Algebra Backend from Eigen
+ typedef EigenLinearAlgebraBackend LinAlg;
+ /// Type of Vector
+ typedef LinAlg::DenseVector Vector;
+ /// Global 0-form, 1-form, 2-form are Vector
+ typedef Vector Form;
+ /// Type of dense matrix
+ typedef LinAlg::DenseMatrix DenseMatrix;
+ /// Type of sparse matrix
+ typedef LinAlg::SparseMatrix SparseMatrix;
+ /// Type of sparse matrix triplet
+ typedef LinAlg::Triplet Triplet;
+
+ /// Type of a sparse matrix solver
+ typedef LinAlg::SolverSimplicialLDLT Solver;
+
+ /// @name Standard services
+ /// @{
+
+ /// Create a Polygonal DEC structure from a surface mesh (@a surf)
+ /// using an default identity embedder.
+ /// @param surf an instance of SurfaceMesh
+ /// @param globalInternalCacheEnabled enable the internal cache for all
+ /// operators (default: false)
+ PolygonalCalculus( const ConstAlias surf,
+ bool globalInternalCacheEnabled = false )
+ : mySurfaceMesh( &surf )
+ , myGlobalCacheEnabled( globalInternalCacheEnabled )
{
- a(i, (i+1)%nf) = 0.5;
- a(i,i) = 0.5;
+ myEmbedder = [ & ]( Face f, Vertex v )
+ {
+ (void)f;
+ return mySurfaceMesh->position( v );
+ };
+ myVertexNormalEmbedder = [ & ]( Vertex v )
+ { return toReal3dVector( computeVertexNormal( v ) ); };
+ init();
+ };
+
+ /// Create a Polygonal DEC structure from a surface mesh (@a surf)
+ /// and an embedder for the vertex position: function with two parameters, a
+ /// face and a vertex which outputs the embedding in R^3 of the vertex
+ /// w.r.t. to the face.
+ /// @param surf an instance of SurfaceMesh
+ /// @param embedder an embedder for the vertex position
+ /// @param globalInternalCacheEnabled if true, the global operator cache is
+ /// enabled
+ PolygonalCalculus(
+ const ConstAlias surf,
+ const std::function & embedder,
+ bool globalInternalCacheEnabled = false )
+ : mySurfaceMesh( &surf )
+ , myEmbedder( embedder )
+ , myGlobalCacheEnabled( globalInternalCacheEnabled )
+ {
+ myVertexNormalEmbedder = [ & ]( Vertex v )
+ { return toReal3dVector( computeVertexNormal( v ) ); };
+ init();
+ };
+
+ /// Create a Polygonal DEC structure from a surface mesh (@a surf)
+ /// and an embedder for the vertex normal: function with a vertex as
+ /// parameter which outputs the embedding in R^3 of the vertex normal.
+ /// @param surf an instance of SurfaceMesh
+ /// @param embedder an embedder for the vertex position
+ /// @param globalInternalCacheEnabled if true, the global operator cache is
+ /// enabled
+ PolygonalCalculus( const ConstAlias surf,
+ const std::function & embedder,
+ bool globalInternalCacheEnabled = false )
+ : mySurfaceMesh( &surf )
+ , myVertexNormalEmbedder( embedder )
+ , myGlobalCacheEnabled( globalInternalCacheEnabled )
+ {
+ myEmbedder = [ & ]( Face f, Vertex v )
+ {
+ (void)f;
+ return mySurfaceMesh->position( v );
+ };
+ init();
+ };
+
+ /// Create a Polygonal DEC structure from a surface mesh (@a surf)
+ /// and an embedder for the vertex position: function with two parameters, a
+ /// face and a vertex which outputs the embedding in R^3 of the vertex
+ /// w.r.t. to the face. and an embedder for the vertex normal: function with
+ /// a vertex as parameter which outputs the embedding in R^3 of the vertex
+ /// normal.
+ /// @param surf an instance of SurfaceMesh
+ /// @param pos_embedder an embedder for the position
+ /// @param normal_embedder an embedder for the position
+ /// @param globalInternalCacheEnabled
+ PolygonalCalculus(
+ const ConstAlias surf,
+ const std::function & pos_embedder,
+ const std::function & normal_embedder,
+ bool globalInternalCacheEnabled = false )
+ : mySurfaceMesh( &surf )
+ , myEmbedder( pos_embedder )
+ , myVertexNormalEmbedder( normal_embedder )
+ , myGlobalCacheEnabled( globalInternalCacheEnabled )
+ {
+ init();
+ };
+
+ /**
+ * Deleted default constructor.
+ */
+ PolygonalCalculus() = delete;
+
+ /**
+ * Destructor (default).
+ */
+ ~PolygonalCalculus() = default;
+
+ /**
+ * Deleted copy constructor.
+ * @param other the object to clone.
+ */
+ PolygonalCalculus( const PolygonalCalculus & other ) = delete;
+
+ /**
+ * Deleted move constructor.
+ * @param other the object to move.
+ */
+ PolygonalCalculus( PolygonalCalculus && other ) = delete;
+
+ /**
+ * Deleted copy assignment operator.
+ * @param other the object to copy.
+ * @return a reference on 'this'.
+ */
+ PolygonalCalculus & operator=( const PolygonalCalculus & other ) = delete;
+
+ /**
+ * Deleted move assignment operator.
+ * @param other the object to move.
+ * @return a reference on 'this'.
+ */
+ PolygonalCalculus & operator=( PolygonalCalculus && other ) = delete;
+
+ /// @}
+
+ // ----------------------- embedding services --------------------------
+ // MARK: Embedding services
+ /// @name Embedding services
+ /// @{
+
+ /// Update the embedding function.
+ /// @param externalFunctor a new embedding functor (Face,Vertex)->RealPoint.
+ void setEmbedder(
+ const std::function & externalFunctor )
+ {
+ myEmbedder = externalFunctor;
}
+ /// @}
- setInCache(A_,f,a);
- return a;
- }
-
-
- /// Polygonal (corrected) vector area.
- /// @param f the face
- /// @return a vector oriented in the (corrected) normal direction and with length equal to the (corrected) area of the face \a f.
- 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;
+ // ----------------------- Per face operators
+ // --------------------------------------
+ // MARK: Per face operator on scalars
+ /// @name Per face operators on scalars
+ /// @{
+
+ /// Return the vertex position matrix degree x 3 of the face.
+ /// @param f a face
+ /// @return the n_f x 3 position matrix
+ 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;
+ }
+
+ /// Derivative operator (d_0) of a face.
+ /// @param f the face
+ /// @return a degree x degree matrix
+ 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;
+ }
+
+ /// Edge vector operator per face.
+ /// @param f the face
+ /// @return degree x 3 matrix
+ 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;
+ }
+
+ /// Average operator to average, per edge, its vertex values.
+ /// @param f the face
+ /// @return a degree x degree matrix
+ 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;
+ }
+
+ /// Polygonal (corrected) vector area.
+ /// @param f the face
+ /// @return a vector oriented in the (corrected) normal direction and with
+ /// length equal to the (corrected) area of the face \a f.
+ 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;
- if (itnext == vertices.cend())
- itnext = vertices.cbegin();
+ 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;
}
- Eigen::Vector3d output = {af[0],af[1],af[2]};
- return 0.5*output;
- }
-
- /// Area of a face from the vector area.
- /// @param f the face
- /// @return the corrected area of the face
- double faceArea(const Face f) const
- {
- return vectorArea(f).norm();
- }
-
- /// Corrected normal vector of a face.
- /// @param f the face
- /// @return a vector (Eigen vector)
- Vector faceNormal(const Face f) const
- {
- Vector v = vectorArea(f);
- v.normalize();
- return v;
- }
-
- /// Corrected normal vector of a face.
- /// @param f the face
- /// @return a vector (DGtal RealVector/RealPoint)
- Real3dVector faceNormalAsDGtalVector(const Face f) const
- {
- Vector v = faceNormal(f);
- return {v(0),v(1),v(2)};
- }
-
- /// co-Gradient operator of the face
- /// @param f the face
- /// @return a 3 x degree matrix
- 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;
- }
-
- ///Return [n] as the 3x3 operator such that [n]q = n x q
- ///@param n a vector
- 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;
- }
-
- /// Gradient operator of the face.
- /// @param f the face
- /// @return 3 x degree matrix
- 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;
- }
-
- /// Flat operator for the face.
- /// @param f the face
- /// @return a degree x 3 matrix
- 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;
- }
-
- /// Edge mid-point operator of the face.
- /// @param f the face
- /// @return a degree x 3 matrix
- 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;
- }
-
- /// @returns the centroid of the face
- /// @param f the face
- Vector centroid(const Face f) const
- {
- const auto nf = myFaceDegree[f];
- return 1.0/(double)nf * X(f).transpose() * Vector::Ones(nf);
- }
-
- /// @returns the centroid of the face as a DGtal RealPoint
- /// @param f the face
- Real3dPoint centroidAsDGtalPoint(const Face f) const
- {
- const Vector c = centroid(f);
- return {c(0),c(1),c(2)};
- }
-
- /// Sharp operator for the face.
- /// @param f the face
- /// @return a 3 x degree matrix
- 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() );
+ /// Area of a face from the vector area.
+ /// @param f the face
+ /// @return the corrected area of the face
+ double faceArea( const Face f ) const
+ {
+ return vectorArea( f ).norm();
+ }
- setInCache(SHARP_,f,op);
- return op;
- }
-
- /// Projection operator for the face.
- /// @param f the face
- /// @return a degree x degree matrix
- 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;
- }
-
- /// Inner product on 1-forms associated with the face
- /// @param f the face
- /// @param lambda the regularization parameter
- /// @return a degree x degree matrix
- 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;
- }
-
- /// Divergence operator of a one-form.
- /// @param f the face
- /// @return a degree x degree matrix
- ///
- /// @note The sign convention for the divergence and the Laplacian
- /// operator is opposite to the one of @cite degoes2020discrete .
- /// 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 \f$
- /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
- /// \rangle \f$. See also
- /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- 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;
- }
-
- /// Curl operator of a one-form (identity matrix).
- /// @param f the face
- /// @return a degree x degree matrix
- DenseMatrix curl(const Face f) const
- {
- if (checkCache(CURL_,f))
- return myGlobalCache[CURL_][f];
-
- DenseMatrix op = DenseMatrix::Identity(myFaceDegree[f],myFaceDegree[f]);
+ /// Corrected normal vector of a face.
+ /// @param f the face
+ /// @return a vector (Eigen vector)
+ Vector faceNormal( const Face f ) const
+ {
+ Vector v = vectorArea( f );
+ v.normalize();
+ return v;
+ }
- setInCache(CURL_,f,op);
- return op;
- }
-
-
- /// (weak) Laplace-Beltrami operator for the face.
- /// @param f the face
- /// @param lambda the regularization parameter
- /// @return a degree x degree matrix
- ///
- /// @note The sign convention for the divergence and the Laplacian
- /// operator is opposite to the one of @cite degoes2020discrete .
- /// 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 \f$
- /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
- /// \rangle \f$. See also https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- 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
- ///@name Per face operators on vector fields
- ///@{
-
-public:
- ///@return 3x2 matrix defining the tangent space at vertex v, with basis
- ///vectors in columns
- 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;
- }
+ /// Corrected normal vector of a face.
+ /// @param f the face
+ /// @return a vector (DGtal RealVector/RealPoint)
+ Real3dVector faceNormalAsDGtalVector( const Face f ) const
+ {
+ Vector v = faceNormal( f );
+ return { v( 0 ), v( 1 ), v( 2 ) };
+ }
- ///@return 3x2 matrix defining the tangent space at face f, with basis
- ///vectors in columns
- 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;
- }
+ /// co-Gradient operator of the face
+ /// @param f the face
+ /// @return a 3 x degree matrix
+ 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;
+ }
- /// \brief toExtrinsicVector
- /// \param v the vertex
- /// \param I the intrinsic vector at Tv
- /// \return 3D extrinsic vector from intrinsic 2D vector I expressed from
- /// tangent frame at vertex v
- Vector toExtrinsicVector(const Vertex v, const Vector & I) const
- {
- DenseMatrix T = Tv(v);
- return T.col(0) * I(0) + T.col(1) * I(1);
- }
+ /// Return [n] as the 3x3 operator such that [n]q = n x q
+ ///@param n a vector
+ 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;
+ }
- /// \param I set of intrinsic vectors, vectors indices must be the same as
- /// their associated vertex
- ///@return converts a set of intrinsic vectors to their extrinsic
- ///equivalent, expressed in correponding tangent frame
- std::vector toExtrinsicVectors(const std::vector & I) const
- {
- std::vector ext(mySurfaceMesh->nbVertices());
- for (auto v = 0; v < mySurfaceMesh->nbVertices(); v++)
- ext[v] = toExtrinsicVector(v,I[v]);
- return ext;
- }
+ /// Gradient operator of the face.
+ /// @param f the face
+ /// @return 3 x degree matrix
+ DenseMatrix gradient( const Face f ) const
+ {
+ if ( checkCache( GRAD_, f ) )
+ return myGlobalCache[ GRAD_ ][ f ];
- /// https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
- ///@return 3x3 Rotation matrix to align n_v to n_f
- 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 op =
+ -1.0 / faceArea( f ) * bracket( faceNormal( f ) ) * coGradient( f );
- ///@return Levi-Civita connection from vertex v tangent space to face f
- ///tangent space (2x2 rotation matrix)
- DenseMatrix Rvf(const Vertex & v, const Face & f) const
- {
- return Tf(f).transpose() * Qvf(v,f) * Tv(v);
- }
+ setInCache( GRAD_, f, op );
+ return op;
+ }
- /// Shape operator on the face @a f (2x2 operator).
- ///@return the shape operator at face f
- DenseMatrix shapeOperator(const Face f) const
- {
- DenseMatrix N(myFaceDegree[f],3);
- uint cpt = 0;
- for (Vertex v : mySurfaceMesh->incidentVertices(f))
+ /// Flat operator for the face.
+ /// @param f the face
+ /// @return a degree x 3 matrix
+ DenseMatrix flat( const Face f ) const
{
- N.block(cpt,0,3,1) = n_v(v).transpose();
- cpt++;
+ 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 GN = gradient(f) * N, Tf = T_f(f);
- return 0.5 * Tf.transpose() * (GN + GN.transpose()) * Tf;
- }
+ /// Edge mid-point operator of the face.
+ /// @param f the face
+ /// @return a degree x 3 matrix
+ 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;
+ }
- /// @return to fit @cite degoes2020discrete paper's notations,
- /// this function maps all the per vertex vectors (expressed in the (2*nf)
- /// vector form) into the nfx2 matrix with transported vectors (to face f) in
- /// each row.
- /// @note Unlike the rest of the per face operators, the
- /// covariant operators need to be applied directly to the restriction of
- /// the vector field to the face,
- DenseMatrix transportAndFormatVectorField(const Face f, const Vector & uf)
- {
- DenseMatrix uf_nabla(myFaceDegree[f], 2);
- size_t cpt = 0;
- for (auto v : mySurfaceMesh->incidentVertices(f))
+ /// @returns the centroid of the face
+ /// @param f the face
+ Vector centroid( const Face f ) const
{
- uf_nabla.block(cpt,0,1,2) =
- (Rvf(v,f) * uf.block(2 * cpt,0,2,1)).transpose();
- ++cpt;
+ const auto nf = myFaceDegree[ f ];
+ return 1.0 / (double)nf * X( f ).transpose() * Vector::Ones( nf );
}
- return uf_nabla;
- }
- /// Covarient gradient at a face a @a f of intrinsic vectors @a uf.
- /// @param uf list of all intrinsic vectors per vertex concatenated in a
- /// column vector
- /// @param f the face
- /// @return the covariant gradient of the given vector field uf (expressed
- /// in corresponding vertex tangent frames), wrt face f
- /// @note Unlike the rest of the per face operators, the
- /// covariant operators need to be applied directly to the restriction of
- /// the vector field to the face,
- DenseMatrix covariantGradient(const Face f, const Vector & uf)
- {
- return Tf(f).transpose() * gradient(f) *
- transportAndFormatVectorField(f,uf);
- }
+ /// @returns the centroid of the face as a DGtal RealPoint
+ /// @param f the face
+ Real3dPoint centroidAsDGtalPoint( const Face f ) const
+ {
+ const Vector c = centroid( f );
+ return { c( 0 ), c( 1 ), c( 2 ) };
+ }
- /// Compute the covariance projection at a face @a f of intrinsic vectors @a uf.
- /// @param uf list of all intrinsic vectors per vertex concatenated in a
- /// column vector
- /// @param f the face
- /// @return the covariant projection of the given vector field uf (
- /// restricted to face f and expressed in corresponding vertex tangent
- /// frames)
- /// @note Unlike the rest of the per face operators, the
- /// covariant operators need to be applied directly to the restriction of
- /// the vector field to the face
- DenseMatrix covariantProjection(const Face f, const Vector & uf)
- {
- return P(f) * D(f) * transportAndFormatVectorField(f,uf);
- }
-
- /// @return Covariant Gradient Operator, returns the operator that acts on
- /// the concatenated vectors. When applied, gives the associated 2x2 matrix
- /// in the isomorphic vector form (a b c d)^t to be used in the dirichlet
- /// energy (vector laplacian) G∇_f.
- /// Used to define the connection Laplacian.
- DenseMatrix covariantGradient_f(const Face & f) const
- {
- return kroneckerWithI2(Tf(f).transpose() * gradient(f)) * blockConnection(f);
- }
-
- /// @return Projection Gradient Operator, returns the operator that acts on
- /// the concatenated vectors. When applied, gives the associated nfx2 matrix
- /// in the isomorphic vector form (a b c d ...)^t to be used in the
- /// dirichlet energy (vector laplacian) P∇_f.
- /// Used to define the connection Laplacian.
- DenseMatrix covariantProjection_f(const Face & f) const
- {
- return kroneckerWithI2(P(f) * D(f)) * blockConnection(f);
- ;
- }
-
- /// L∇ := -(afG∇tG∇+λP∇tP∇)
- /// @return Connection/Vector laplacian at face f
- /// @note The sign convention for the divergence and the Laplacian
- /// operator is opposite to the one of @cite degoes2020discrete.
- /// 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 \f$
- /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
- /// \rangle \f$. See also
- /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- 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
- /// @name Global operators
- /// @{
-
- /// @return a 0-form initialized to zero
- Form form0() const
- {
- return Form::Zero( nbVertices() );
- }
- /// @return the identity linear operator for 0-forms
- SparseMatrix identity0() const
- {
- SparseMatrix Id0( nbVertices(), nbVertices() );
- Id0.setIdentity();
- return Id0;
- }
-
- /// Computes the global Laplace-Beltrami operator by assembling the
- /// per face operators.
- ///
- /// @param lambda the regularization parameter for the local Laplace-Beltrami operators
- /// @return a sparse nbVertices x nbVertices matrix
- ///
- /// @note The sign convention for the divergence is opposite to the
- /// one of @cite degoes2020discrete. This is also true for the
- /// Laplacian operator. 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 \f$
- /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
- /// \rangle \f$. See also https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
- {
- SparseMatrix lapGlobal(mySurfaceMesh->nbVertices(), mySurfaceMesh->nbVertices());
- std::vector 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 ) ) );
- }
+ /// Sharp operator for the face.
+ /// @param f the face
+ /// @return a 3 x degree matrix
+ 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;
}
- lapGlobal.setFromTriplets(triplets.begin(), triplets.end());
- return lapGlobal;
- }
-
- /// Compute and returns the global lumped mass matrix
- /// (diagonal matrix with Max's weights for each vertex).
- /// M(i,i) = ∑_{adjface f} faceArea(f)/degree(f) ;
- ///
- /// @return the global lumped mass matrix.
- SparseMatrix globalLumpedMassMatrix() const
- {
- SparseMatrix M(mySurfaceMesh->nbVertices(), mySurfaceMesh->nbVertices());
- std::vector triplets;
- for ( typename MySurfaceMesh::Index v = 0; v < mySurfaceMesh->nbVertices(); ++v )
+
+ /// Projection operator for the face.
+ /// @param f the face
+ /// @return a degree x degree matrix
+ 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;
+ }
+
+ /// Inner product on 1-forms associated with the face
+ /// @param f the face
+ /// @param lambda the regularization parameter
+ /// @return a degree x degree matrix
+ 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;
+ }
+
+ /// Divergence operator of a one-form.
+ /// @param f the face
+ /// @return a degree x degree matrix
+ ///
+ /// @note The sign convention for the divergence and the Laplacian
+ /// operator is opposite to the one of @cite degoes2020discrete .
+ /// 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ 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;
+ }
+
+ /// Curl operator of a one-form (identity matrix).
+ /// @param f the face
+ /// @return a degree x degree matrix
+ 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;
+ }
+
+ /// (weak) Laplace-Beltrami operator for the face.
+ /// @param f the face
+ /// @param lambda the regularization parameter
+ /// @return a degree x degree matrix
+ ///
+ /// @note The sign convention for the divergence and the Laplacian
+ /// operator is opposite to the one of @cite degoes2020discrete .
+ /// 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ 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
+ ///@name Per face operators on vector fields
+ ///@{
+
+ public:
+ ///@return 3x2 matrix defining the tangent space at vertex v, with basis
+ /// vectors in columns
+ 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;
+ }
+
+ ///@return 3x2 matrix defining the tangent space at face f, with basis
+ /// vectors in columns
+ 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;
+ }
+
+ /// \brief toExtrinsicVector
+ /// \param v the vertex
+ /// \param I the intrinsic vector at Tv
+ /// \return 3D extrinsic vector from intrinsic 2D vector I expressed from
+ /// tangent frame at vertex v
+ Vector toExtrinsicVector( const Vertex v, const Vector & I ) const
+ {
+ DenseMatrix T = Tv( v );
+ return T.col( 0 ) * I( 0 ) + T.col( 1 ) * I( 1 );
+ }
+
+ /// \param I set of intrinsic vectors, vectors indices must be the same as
+ /// their associated vertex
+ ///@return converts a set of intrinsic vectors to their extrinsic
+ /// equivalent, expressed in correponding tangent frame
+ std::vector
+ toExtrinsicVectors( const std::vector & I ) const
+ {
+ std::vector ext( mySurfaceMesh->nbVertices() );
+ for ( auto v = 0; v < mySurfaceMesh->nbVertices(); v++ )
+ ext[ v ] = toExtrinsicVector( v, I[ v ] );
+ return ext;
+ }
+
+ /// https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
+ ///@return 3x3 Rotation matrix to align n_v to n_f
+ 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;
+ }
+
+ ///@return Levi-Civita connection from vertex v tangent space to face f
+ /// tangent space (2x2 rotation matrix)
+ DenseMatrix Rvf( const Vertex & v, const Face & f ) const
+ {
+ return Tf( f ).transpose() * Qvf( v, f ) * Tv( v );
+ }
+
+ /// Shape operator on the face @a f (2x2 operator).
+ ///@return the shape operator at face f
+ DenseMatrix shapeOperator( const Face f ) const
+ {
+ DenseMatrix N( myFaceDegree[ f ], 3 );
+ uint cpt = 0;
+ for ( Vertex v : mySurfaceMesh->incidentVertices( f ) )
{
- auto faces = mySurfaceMesh->incidentFaces(v);
+ 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;
+ }
+
+ /// @return to fit @cite degoes2020discrete paper's notations,
+ /// this function maps all the per vertex vectors (expressed in the (2*nf)
+ /// vector form) into the nfx2 matrix with transported vectors (to face f)
+ /// in each row.
+ /// @note Unlike the rest of the per face operators, the
+ /// covariant operators need to be applied directly to the restriction of
+ /// the vector field to the face,
+ 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;
+ }
+
+ /// Covarient gradient at a face a @a f of intrinsic vectors @a uf.
+ /// @param uf list of all intrinsic vectors per vertex concatenated in a
+ /// column vector
+ /// @param f the face
+ /// @return the covariant gradient of the given vector field uf (expressed
+ /// in corresponding vertex tangent frames), wrt face f
+ /// @note Unlike the rest of the per face operators, the
+ /// covariant operators need to be applied directly to the restriction of
+ /// the vector field to the face,
+ DenseMatrix covariantGradient( const Face f, const Vector & uf )
+ {
+ return Tf( f ).transpose() * gradient( f ) *
+ transportAndFormatVectorField( f, uf );
+ }
+
+ /// Compute the covariance projection at a face @a f of intrinsic vectors @a
+ /// uf.
+ /// @param uf list of all intrinsic vectors per vertex concatenated in a
+ /// column vector
+ /// @param f the face
+ /// @return the covariant projection of the given vector field uf (
+ /// restricted to face f and expressed in corresponding vertex tangent
+ /// frames)
+ /// @note Unlike the rest of the per face operators, the
+ /// covariant operators need to be applied directly to the restriction of
+ /// the vector field to the face
+ DenseMatrix covariantProjection( const Face f, const Vector & uf )
+ {
+ return P( f ) * D( f ) * transportAndFormatVectorField( f, uf );
+ }
+
+ /// @return Covariant Gradient Operator, returns the operator that acts on
+ /// the concatenated vectors. When applied, gives the associated 2x2 matrix
+ /// in the isomorphic vector form (a b c d)^t to be used in the dirichlet
+ /// energy (vector laplacian) G∇_f.
+ /// Used to define the connection Laplacian.
+ DenseMatrix covariantGradient_f( const Face & f ) const
+ {
+ return kroneckerWithI2( Tf( f ).transpose() * gradient( f ) ) *
+ blockConnection( f );
+ }
+
+ /// @return Projection Gradient Operator, returns the operator that acts on
+ /// the concatenated vectors. When applied, gives the associated nfx2 matrix
+ /// in the isomorphic vector form (a b c d ...)^t to be used in the
+ /// dirichlet energy (vector laplacian) P∇_f.
+ /// Used to define the connection Laplacian.
+ DenseMatrix covariantProjection_f( const Face & f ) const
+ {
+ return kroneckerWithI2( P( f ) * D( f ) ) * blockConnection( f );
+ ;
+ }
+
+ /// L∇ := -(afG∇tG∇+λP∇tP∇)
+ /// @return Connection/Vector laplacian at face f
+ /// @note The sign convention for the divergence and the Laplacian
+ /// operator is opposite to the one of @cite degoes2020discrete.
+ /// 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ 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
+ /// @name Global operators
+ /// @{
+
+ /// @return a 0-form initialized to zero
+ Form form0() const
+ {
+ return Form::Zero( nbVertices() );
+ }
+ /// @return the identity linear operator for 0-forms
+ SparseMatrix identity0() const
+ {
+ SparseMatrix Id0( nbVertices(), nbVertices() );
+ Id0.setIdentity();
+ return Id0;
+ }
+
+ /// Computes the global Laplace-Beltrami operator by assembling the
+ /// per face operators.
+ ///
+ /// @param lambda the regularization parameter for the local
+ /// Laplace-Beltrami operators
+ /// @return a sparse nbVertices x nbVertices matrix
+ ///
+ /// @note The sign convention for the divergence is opposite to the
+ /// one of @cite degoes2020discrete. This is also true for the
+ /// Laplacian operator. 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ SparseMatrix globalLaplaceBeltrami( const double lambda = 1.0 ) const
+ {
+ SparseMatrix lapGlobal( mySurfaceMesh->nbVertices(),
+ mySurfaceMesh->nbVertices() );
+ std::vector 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;
+ }
+
+ /// Compute and returns the global lumped mass matrix
+ /// (diagonal matrix with Max's weights for each vertex).
+ /// M(i,i) = ∑_{adjface f} faceArea(f)/degree(f) ;
+ ///
+ /// @return the global lumped mass matrix.
+ SparseMatrix globalLumpedMassMatrix() const
+ {
+ SparseMatrix M( mySurfaceMesh->nbVertices(),
+ mySurfaceMesh->nbVertices() );
+ std::vector 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));
+ 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;
- }
+ M.setFromTriplets( triplets.begin(), triplets.end() );
+ return M;
+ }
- /// Compute and returns the inverse of the global lumped mass matrix
- /// (diagonal matrix with Max's weights for each vertex).
- ///
- /// @return the inverse of the global lumped mass matrix.
- 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;
- }
+ /// Compute and returns the inverse of the global lumped mass matrix
+ /// (diagonal matrix with Max's weights for each vertex).
+ ///
+ /// @return the inverse of the global lumped mass matrix.
+ 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;
+ }
- /// Computes the global Connection-Laplace-Beltrami operator by accumulating
- /// the per face operators.
- ///
- /// @param lambda the regualrization parameter for the local
- /// Connection-Laplace-Beltrami operators
- /// @return a sparse 2*nbVertices x 2*nbVertices matrix
- ///
- /// @note The sign convention for the divergence is opposite to the
- /// one of @cite degoes2020discrete. This is also true for the
- /// Laplacian operator. 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 \f$
- /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
- /// \rangle \f$. See also
- /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
- 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 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;
- }
+ /// Computes the global Connection-Laplace-Beltrami operator by accumulating
+ /// the per face operators.
+ ///
+ /// @param lambda the regualrization parameter for the local
+ /// Connection-Laplace-Beltrami operators
+ /// @return a sparse 2*nbVertices x 2*nbVertices matrix
+ ///
+ /// @note The sign convention for the divergence is opposite to the
+ /// one of @cite degoes2020discrete. This is also true for the
+ /// Laplacian operator. 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 \f$
+ /// \langle \mathrm{d} u, v \rangle = - \langle u, \mathrm{div} v
+ /// \rangle \f$. See also
+ /// https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+ 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 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;
+ }
- /// Compute and returns the global lumped mass matrix tensorized with Id_2
- /// (used for connection laplacian) (diagonal matrix with Max's weights for
- /// each vertex).
- /// M(2*i,2*i) = ∑_{adjface f} faceArea(f)/degree(f) ;
- /// M(2*i+1,2*i+1) = M(2*i,2*i)
- /// @return the global lumped mass matrix.
- SparseMatrix doubledGlobalLumpedMassMatrix() const
- {
- auto nv = mySurfaceMesh->nbVertices();
- SparseMatrix M(2 * nv, 2 * nv);
- std::vector 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 --------------------------------------
- /// @name Cache mechanism
- /// @{
-
- /// Generic method to compute all the per face DenseMatrices and store them in an
- /// indexed container.
- ///
- /// Usage example:
- /// @code
- ///auto opM = [&](const PolygonalCalculus::Face f){ return calculus.M(f);};
- ///auto cacheM = boxCalculus.getOperatorCacheMatrix(opM);
- ///...
- /////Now you have access to the cached values and mixed them with un-cached ones
- /// Face f = ...;
- /// auto res = cacheM[f] * calculus.D(f) * phi;
- /// ...
- ///@endcode
- ///
- /// @param perFaceOperator the per face operator
- /// @return an indexed container of all DenseMatrix operators (indexed per Face).
- std::vector getOperatorCacheMatrix(const std::function &perFaceOperator) const
- {
- std::vector cache;
- for( typename MySurfaceMesh::Index f=0; f < mySurfaceMesh->nbFaces(); ++f)
- cache.push_back(perFaceOperator(f));
- return cache;
- }
-
- /// Generic method to compute all the per face Vector and store them in an
- /// indexed container.
- ///
- /// Usage example:
- /// @code
- ///auto opCentroid = [&](const PolygonalCalculus::Face f){ return calculus.centroid(f);};
- ///auto cacheCentroid = boxCalculus.getOperatorCacheVector(opCentroid);
- ///...
- /////Now you have access to the cached values and mixed them with un-cached ones
- /// Face f = ...;
- /// auto res = calculus.P(f) * cacheCentroid[f] ;
- /// ...
- ///@endcode
- ///
- /// @param perFaceVectorOperator the per face operator
- /// @return an indexed container of all Vector quantities (indexed per Face).
- std::vector getOperatorCacheVector(const std::function &perFaceVectorOperator) const
- {
- std::vector cache;
- for( typename MySurfaceMesh::Index f=0; f < mySurfaceMesh->nbFaces(); ++f)
- cache.push_back(perFaceVectorOperator(f));
- return cache;
- }
+ /// Compute and returns the global lumped mass matrix tensorized with Id_2
+ /// (used for connection laplacian) (diagonal matrix with Max's weights for
+ /// each vertex).
+ /// M(2*i,2*i) = ∑_{adjface f} faceArea(f)/degree(f) ;
+ /// M(2*i+1,2*i+1) = M(2*i,2*i)
+ /// @return the global lumped mass matrix.
+ SparseMatrix doubledGlobalLumpedMassMatrix() const
+ {
+ auto nv = mySurfaceMesh->nbVertices();
+ SparseMatrix M( 2 * nv, 2 * nv );
+ std::vector 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;
+ }
+ /// @}
- /// Enable the internal global cache for operators.
- ///
- void enableInternalGlobalCache()
- {
- myGlobalCacheEnabled = true;
- }
-
- /// Disable the internal global cache for operators.
- /// This method will also clean up the
- void disableInternalGlobalCache()
- {
- myGlobalCacheEnabled = false;
- myGlobalCache.clear();
- }
+ // ----------------------- Cache mechanism
+ // --------------------------------------
+ /// @name Cache mechanism
+ /// @{
- /// @}
-
- // ----------------------- Common --------------------------------------
-public:
- /// @name Common services
- /// @{
-
- /// Update the internal cache structures
- /// (e.g. degree of each face).
- void init()
- {
- updateFaceDegree();
- }
-
- /// Helper to retrieve the degree of the face from the cache.
- /// @param f the face
- /// @return the number of vertices of the face.
- size_t faceDegree(Face f) const
- {
- return myFaceDegree[f];
- }
-
- /// @return the number of vertices of the underlying surface mesh.
- size_t nbVertices() const
- {
- return mySurfaceMesh->nbVertices();
- }
-
- /// @return the number of faces of the underlying surface mesh.
- size_t nbFaces() const
- {
- return mySurfaceMesh->nbFaces();
- }
-
- /// @returns the degree of the face f (number of vertices)
- /// @param f the face
- size_t degree(const Face f) const
- {
- return myFaceDegree[f];
- }
-
- /// @returns an pointer to the underlying SurfaceMash object.
- const MySurfaceMesh * getSurfaceMeshPtr() const
- {
- return mySurfaceMesh;
- }
-
- /**
- * Writes/Displays the object on an output stream.
- * @param out the output stream where the object is written.
- */
- void selfDisplay ( std::ostream & out ) const
- {
- out << "[PolygonalCalculus]: ";
- if (myGlobalCacheEnabled)
- out<< "internal cache enabled, ";
- else
- out<<"internal cache disabled, ";
- out <<"SurfaceMesh="<<*mySurfaceMesh;
- }
-
- /**
- * Checks the validity/consistency of the object.
- * @return 'true' if the object is valid, 'false' otherwise.
- */
- bool isValid() const
- {
- return true;
- }
+ /// Generic method to compute all the per face DenseMatrices and store them
+ /// in an indexed container.
+ ///
+ /// Usage example:
+ /// @code
+ /// auto opM = [&](const PolygonalCalculus::Face f){ return
+ /// calculus.M(f);}; auto cacheM = boxCalculus.getOperatorCacheMatrix(opM);
+ ///...
+ /////Now you have access to the cached values and mixed them with un-cached
+ /// ones
+ /// Face f = ...;
+ /// auto res = cacheM[f] * calculus.D(f) * phi;
+ /// ...
+ ///@endcode
+ ///
+ /// @param perFaceOperator the per face operator
+ /// @return an indexed container of all DenseMatrix operators (indexed per
+ /// Face).
+ std::vector getOperatorCacheMatrix(
+ const std::function & perFaceOperator ) const
+ {
+ std::vector cache;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ cache.push_back( perFaceOperator( f ) );
+ return cache;
+ }
- /// @}
-
- // ------------------------- Protected Datas ------------------------------
- //MARK: Protected
-
-protected:
- /// @name Protected services and types
- /// @{
-
- ///Enum for operators in the internal cache strategy
- enum OPERATOR { X_, D_, E_, A_, COGRAD_, GRAD_, FLAT_, B_, SHARP_, P_, M_, DIVERGENCE_, CURL_, L_,CON_L_};
-
- /// Update the face degree cache
- void updateFaceDegree()
- {
- myFaceDegree.resize(mySurfaceMesh->nbFaces());
- for(typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces(); ++f)
+ /// Generic method to compute all the per face Vector and store them in an
+ /// indexed container.
+ ///
+ /// Usage example:
+ /// @code
+ /// auto opCentroid = [&](const PolygonalCalculus::Face f){ return
+ /// calculus.centroid(f);}; auto cacheCentroid =
+ /// boxCalculus.getOperatorCacheVector(opCentroid);
+ ///...
+ /////Now you have access to the cached values and mixed them with un-cached
+ /// ones
+ /// Face f = ...;
+ /// auto res = calculus.P(f) * cacheCentroid[f] ;
+ /// ...
+ ///@endcode
+ ///
+ /// @param perFaceVectorOperator the per face operator
+ /// @return an indexed container of all Vector quantities (indexed per
+ /// Face).
+ std::vector getOperatorCacheVector(
+ const std::function & perFaceVectorOperator ) const
{
- auto vertices = mySurfaceMesh->incidentVertices(f);
- auto nf = vertices.size();
- myFaceDegree[f] = nf;
+ std::vector cache;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ cache.push_back( perFaceVectorOperator( f ) );
+ return cache;
}
- }
-
- /// Check internal cache if enabled.
- /// @param key the operator name
- /// @param f the face
- /// @returns true if the operator "key" for the face f has been computed.
- bool checkCache(OPERATOR key, const Face f) const
- {
- if (myGlobalCacheEnabled)
- if (myGlobalCache[key].find(f) != myGlobalCache[key].end())
- return true;
- return false;
- }
- /// Set an operator in the internal cache.
- /// @param key the operator name
- /// @param f the face
- /// @param ope the operator to store
- void setInCache(OPERATOR key, const Face f,
- const DenseMatrix &ope) const
- {
- if (myGlobalCacheEnabled)
- myGlobalCache[key][f] = ope;
- }
-
- /// Project u on the orthgonal of n
- /// \param u vector to project
- /// \param n vector to build orthogonal space from
- /// \return projected vector
- static Vector project(const Vector & u, const Vector & n)
- {
- return u - (u.dot(n) / n.squaredNorm()) * n;
- }
-
- /// Conversion routines.
- /// \brief toVector convert Real3dPoint to Eigen::VectorXd
- /// \param x the vector
- /// \return the same vector in eigen type
- static Vector toVector(const Eigen::Vector3d & x)
- {
- Vector X(3);
- for (int i = 0; i < 3; i++)
- X(i) = x(i);
- return X;
- }
-
- /// \brief toVec3 convert Real3dPoint to Eigen::Vector3d
- /// \param x the vector
- /// \return the same vector in eigen type
- static Eigen::Vector3d toVec3(const Real3dPoint & x)
- {
- return Eigen::Vector3d(x(0),x(1),x(2));
- }
-
- /// Conversion routines.
+ /// Enable the internal global cache for operators.
+ ///
+ void enableInternalGlobalCache()
+ {
+ myGlobalCacheEnabled = true;
+ }
+
+ /// Disable the internal global cache for operators.
+ /// This method will also clean up the
+ void disableInternalGlobalCache()
+ {
+ myGlobalCacheEnabled = false;
+ myGlobalCache.clear();
+ }
+
+ /// @}
+
+ // ----------------------- Common --------------------------------------
+ public:
+ /// @name Common services
+ /// @{
+
+ /// Update the internal cache structures
+ /// (e.g. degree of each face).
+ void init()
+ {
+ updateFaceDegree();
+ }
+
+ /// Helper to retrieve the degree of the face from the cache.
+ /// @param f the face
+ /// @return the number of vertices of the face.
+ size_t faceDegree( Face f ) const
+ {
+ return myFaceDegree[ f ];
+ }
+
+ /// @return the number of vertices of the underlying surface mesh.
+ size_t nbVertices() const
+ {
+ return mySurfaceMesh->nbVertices();
+ }
+
+ /// @return the number of faces of the underlying surface mesh.
+ size_t nbFaces() const
+ {
+ return mySurfaceMesh->nbFaces();
+ }
+
+ /// @returns the degree of the face f (number of vertices)
+ /// @param f the face
+ size_t degree( const Face f ) const
+ {
+ return myFaceDegree[ f ];
+ }
+
+ /// @returns an pointer to the underlying SurfaceMash object.
+ const MySurfaceMesh * getSurfaceMeshPtr() const
+ {
+ return mySurfaceMesh;
+ }
+
+ /**
+ * Writes/Displays the object on an output stream.
+ * @param out the output stream where the object is written.
+ */
+ void selfDisplay( std::ostream & out ) const
+ {
+ out << "[PolygonalCalculus]: ";
+ if ( myGlobalCacheEnabled )
+ out << "internal cache enabled, ";
+ else
+ out << "internal cache disabled, ";
+ out << "SurfaceMesh=" << *mySurfaceMesh;
+ }
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const
+ {
+ return true;
+ }
+
+ /// @}
+
+ // ------------------------- Protected Datas ------------------------------
+ // MARK: Protected
+
+ protected:
+ /// @name Protected services and types
+ /// @{
+
+ /// Enum for operators in the internal cache strategy
+ enum OPERATOR
+ {
+ X_,
+ D_,
+ E_,
+ A_,
+ COGRAD_,
+ GRAD_,
+ FLAT_,
+ B_,
+ SHARP_,
+ P_,
+ M_,
+ DIVERGENCE_,
+ CURL_,
+ L_,
+ CON_L_
+ };
+
+ /// Update the face degree cache
+ 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;
+ }
+ }
+
+ /// Check internal cache if enabled.
+ /// @param key the operator name
+ /// @param f the face
+ /// @returns true if the operator "key" for the face f has been computed.
+ bool checkCache( OPERATOR key, const Face f ) const
+ {
+ if ( myGlobalCacheEnabled )
+ if ( myGlobalCache[ key ].find( f ) != myGlobalCache[ key ].end() )
+ return true;
+ return false;
+ }
+
+ /// Set an operator in the internal cache.
+ /// @param key the operator name
+ /// @param f the face
+ /// @param ope the operator to store
+ void setInCache( OPERATOR key, const Face f, const DenseMatrix & ope ) const
+ {
+ if ( myGlobalCacheEnabled )
+ myGlobalCache[ key ][ f ] = ope;
+ }
+
+ /// Project u on the orthgonal of n
+ /// \param u vector to project
+ /// \param n vector to build orthogonal space from
+ /// \return projected vector
+ static Vector project( const Vector & u, const Vector & n )
+ {
+ return u - ( u.dot( n ) / n.squaredNorm() ) * n;
+ }
+
+ /// Conversion routines.
+ /// \brief toVector convert Real3dPoint to Eigen::VectorXd
+ /// \param x the vector
+ /// \return the same vector in eigen type
+ static Vector toVector( const Eigen::Vector3d & x )
+ {
+ Vector X( 3 );
+ for ( int i = 0; i < 3; i++ )
+ X( i ) = x( i );
+ return X;
+ }
+
+ /// \brief toVec3 convert Real3dPoint to Eigen::Vector3d
+ /// \param x the vector
+ /// \return the same vector in eigen type
+ static Eigen::Vector3d toVec3( const Real3dPoint & x )
+ {
+ return Eigen::Vector3d( x( 0 ), x( 1 ), x( 2 ) );
+ }
+
+ /// Conversion routines.
/// \brief toReal3dVector converts Eigen::Vector3d to Real3dVector.
/// \param x the vector
/// \return the same vector in DGtal type
- static Real3dVector toReal3dVector(const Eigen::Vector3d & x)
- {
- return { x(0), x(1), x(2)};
- }
-
-
- /// Compute the (normalized) normal vector at a Vertex by averaging
- /// the adjacent face normal vectors.
- /// \param v the vertex to compute the normal from
- /// \return 3D normal vector at vertex v
- ///
- 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."<incidentVertices(f))
+ /// Compute the (normalized) normal vector at a Vertex by averaging
+ /// the adjacent face normal vectors.
+ /// \param v the vertex to compute the normal from
+ /// \return 3D normal vector at vertex v
+ ///
+ Vector computeVertexNormal( const Vertex & v ) const
{
- auto Rv = Rvf(v,f);
- RU_fO.block<2,2>(2 * cpt,2 * cpt) = Rv;
- ++cpt;
+ 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;
}
- return RU_fO;
- }
-
- /// @return the tensor-kronecker product of M with 2x2 identity matrix
- 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++)
+
+ ///@return the normal vector at vertex v, if no normal vertex embedder is
+ /// set, the normal will be computed
+ Eigen::Vector3d n_v( const Vertex & v ) const
+ {
+ return toVec3( myVertexNormalEmbedder( v ) );
+ }
+
+ // Covariant operators routines
+ /// @return Block Diagonal matrix with Rvf for each vertex v in face f
+ 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 ) )
{
- MK(2 * j, 2 * i) = M(j, i);
- MK(2 * j + 1, 2 * i + 1) = M(j, i);
+ auto Rv = Rvf( v, f );
+ RU_fO.block<2, 2>( 2 * cpt, 2 * cpt ) = Rv;
+ ++cpt;
}
- return MK;
- }
+ return RU_fO;
+ }
-
-
- /// @}
-
- // ------------------------- Internals ------------------------------------
- //MARK: Internals
-private:
-
- ///Underlying SurfaceMesh
- const MySurfaceMesh *mySurfaceMesh;
-
- ///Embedding function (face,vertex)->R^3 for the vertex position wrt. the face.
- std::function myEmbedder;
-
- ///Embedding function (vertex)->R^3 for the vertex normal.
- std::function myVertexNormalEmbedder;
-
- ///Cache containing the face degree
- std::vector myFaceDegree;
-
- ///Global cache
- bool myGlobalCacheEnabled;
- mutable std::array, 15> myGlobalCache;
-
-}; // end of class PolygonalCalculus
+ /// @return the tensor-kronecker product of M with 2x2 identity matrix
+ 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;
+ }
-/**
- * Overloads 'operator<<' for displaying objects of class 'PolygonalCalculus'.
- * @param out the output stream where the object is written.
- * @param object the object of class 'PolygonalCalculus' to write.
- * @return the output stream after the writing.
- */
-template
-std::ostream&
-operator<< ( std::ostream & out, const PolygonalCalculus & object )
-{
- object.selfDisplay( out );
- return out;
-}
+ /// @}
+
+ // ------------------------- Internals ------------------------------------
+ // MARK: Internals
+ private:
+ /// Underlying SurfaceMesh
+ const MySurfaceMesh * mySurfaceMesh;
+
+ /// Embedding function (face,vertex)->R^3 for the vertex position wrt. the
+ /// face.
+ std::function myEmbedder;
+
+ /// Embedding function (vertex)->R^3 for the vertex normal.
+ std::function myVertexNormalEmbedder;
+
+ /// Cache containing the face degree
+ std::vector myFaceDegree;
+
+ /// Global cache
+ bool myGlobalCacheEnabled;
+ mutable std::array, 15> myGlobalCache;
+
+ DenseMatrix buildLocalD0( Face f ) const
+ {
+ return D( f );
+ }
+
+ DenseMatrix buildLocalM0( Face f ) const
+ {
+ auto nf = degree( f );
+ return DenseMatrix::Identity( nf, nf ) * faceArea( f ) / nf;
+ }
+
+ DenseMatrix buildLocalM1( Face f ) const
+ {
+ return M( f );
+ }
+
+ DenseMatrix buildLocalM2( Face f ) const
+ {
+ return DenseMatrix::Identity( 1, 1 ) / faceArea( f );
+ }
+
+ DenseMatrix buildLocalSharp( Face f ) const
+ {
+ return sharp( f );
+ }
+
+ DenseMatrix buildLocalFlat( Face f ) const
+ {
+ return flat( f );
+ }
+
+ DenseMatrix buildLocalL0( Face f ) const
+ {
+ return laplaceBeltrami( f );
+ }
+
+ SparseMatrix buildD0() const
+ {
+ std::vector triplets;
+ // trace.beginBlock( "Init derivative operator D0" );
+ typename MySurfaceMesh::Index e = 0;
+ for ( auto && vtcs : mySurfaceMesh->allEdgeVertices() )
+ {
+ triplets.push_back( {(SparseMatrix::StorageIndex)e,(SparseMatrix::StorageIndex)vtcs.first, -1 } );
+ triplets.push_back( {(SparseMatrix::StorageIndex)e,(SparseMatrix::StorageIndex)vtcs.second, 1 } );
+ e++;
+ }
+ SparseMatrix myD0 =
+ SparseMatrix( mySurfaceMesh->nbEdges(), mySurfaceMesh->nbVertices() );
+ myD0.setFromTriplets( triplets.cbegin(), triplets.cend() );
+ return myD0;
+ }
+
+ SparseMatrix buildM0() const
+ {
+ return globalLumpedMassMatrix();
+ }
+
+ SparseMatrix buildM1() const
+ {
+ SparseMatrix mGlobal( mySurfaceMesh->nbEdges(),
+ mySurfaceMesh->nbEdges() );
+ std::vector triplets;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ {
+ auto nf = myFaceDegree[ f ];
+ DenseMatrix M = this->M( f );
+ const auto vertices = mySurfaceMesh->incidentVertices( f );
+ for ( auto i = 0u; i < nf; ++i )
+ for ( auto j = 0u; j < nf; ++j )
+ {
+ auto v = M( i, j );
+ auto edge_i = mySurfaceMesh->makeEdge( vertices[ i ],
+ vertices[ ( i + 1 ) % nf ] );
+ auto edge_j = mySurfaceMesh->makeEdge( vertices[ j ],
+ vertices[ ( j + 1 ) % nf ] );
+ if ( v != 0.0 )
+ triplets.emplace_back(
+ Triplet( (SparseMatrix::StorageIndex)edge_i,
+ (SparseMatrix::StorageIndex)edge_j, M( i, j ) ) );
+ }
+ }
+ mGlobal.setFromTriplets( triplets.begin(), triplets.end() );
+ return mGlobal;
+ }
+
+ SparseMatrix buildM2() const
+ {
+ SparseMatrix mGlobal( mySurfaceMesh->nbFaces(),
+ mySurfaceMesh->nbFaces() );
+ std::vector triplets;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ {
+ triplets.emplace_back( Triplet( (SparseMatrix::StorageIndex)f,
+ (SparseMatrix::StorageIndex)f,
+ this->faceArea( f ) ) );
+ }
+ mGlobal.setFromTriplets( triplets.begin(), triplets.end() );
+ return mGlobal;
+ }
+
+ SparseMatrix buildSharp() const
+ {
+ SparseMatrix mySharp( 3 * mySurfaceMesh->nbFaces(),
+ mySurfaceMesh->nbEdges() );
+
+ std::vector triplets;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ {
+ auto nf = myFaceDegree[ f ];
+ DenseMatrix sharp = this->sharp( f );
+ const auto vertices = mySurfaceMesh->incidentVertices( f );
+ for ( auto i = 0u; i < nf; ++i )
+ for ( auto j = 0; j < 3; ++j )
+ {
+ auto v = sharp( j, i );
+ auto edge_i = mySurfaceMesh->makeEdge( vertices[ i ],
+ vertices[ ( i + 1 ) % nf ] );
+ if(vertices[i] < vertices[(i + 1) % nf])
+ triplets.emplace_back( Triplet(
+ (SparseMatrix::StorageIndex)3 * f + j,
+ (SparseMatrix::StorageIndex)edge_i, v ) );
+ else
+ triplets.emplace_back( Triplet(
+ (SparseMatrix::StorageIndex)3 * f + j,
+ (SparseMatrix::StorageIndex)edge_i, -v ) );
+ }
+ }
+ mySharp.setFromTriplets( triplets.begin(), triplets.end() );
+ return mySharp;
+ }
+
+ SparseMatrix buildFlat() const
+ {
+ SparseMatrix myFlat( mySurfaceMesh->nbEdges(),
+ 3 * mySurfaceMesh->nbFaces() );
+ std::vector triplets;
+ for ( typename MySurfaceMesh::Index f = 0; f < mySurfaceMesh->nbFaces();
+ ++f )
+ {
+ auto nf = myFaceDegree[ f ];
+ DenseMatrix flat = this->flat( f );
+ const auto vertices = mySurfaceMesh->incidentVertices( f );
+ for ( auto i = 0u; i < nf; ++i )
+ for ( auto j = 0; j < 3; ++j )
+ {
+ auto edge_i = mySurfaceMesh->makeEdge( vertices[ i ],
+ vertices[ ( i + 1 ) % nf ] );
+ auto v = flat( i, j ) / (double) mySurfaceMesh->edgeFaces(edge_i).size();
+ if(vertices[i] < vertices[(i + 1) % nf])
+ triplets.emplace_back( Triplet(
+ (SparseMatrix::StorageIndex)edge_i,
+ (SparseMatrix::StorageIndex)3 * f + j,
+ v ) );
+ else
+ triplets.emplace_back( Triplet(
+ (SparseMatrix::StorageIndex)edge_i,
+ (SparseMatrix::StorageIndex)3 * f + j,
+ -v ) );
+ }
+ }
+ myFlat.setFromTriplets( triplets.begin(), triplets.end() );
+ return myFlat;
+ }
+
+ SparseMatrix buildL0() const
+ {
+ return globalLaplaceBeltrami();
+ }
+
+ }; // end of class PolygonalCalculus
+
+ /**
+ * Overloads 'operator<<' for displaying objects of class 'PolygonalCalculus'.
+ * @param out the output stream where the object is written.
+ * @param object the object of class 'PolygonalCalculus' to write.
+ * @return the output stream after the writing.
+ */
+ template
+ std::ostream & operator<<( std::ostream & out,
+ const PolygonalCalculus & object )
+ {
+ object.selfDisplay( out );
+ return out;
+ }
} // namespace DGtal
///////////////////////////////////////////////////////////////////////////////
diff --git a/src/DGtal/dec/SurfaceDEC.h b/src/DGtal/dec/SurfaceDEC.h
new file mode 100644
index 0000000000..fccd7e21cf
--- /dev/null
+++ b/src/DGtal/dec/SurfaceDEC.h
@@ -0,0 +1,182 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+#pragma once
+
+/**
+ * @file
+ * @author
+ *
+ * @date 2024/06/21
+ *
+ * Header file for module SurfaceDEC.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if !defined SurfaceDEC_h
+/** Prevents repeated inclusion of headers. */
+#define SurfaceDEC_h
+
+#include "DGtal/math/linalg/EigenSupport.h"
+namespace DGtal
+{
+
+ /**
+ Description of template class 'SurfaceDEC' \brief
+ Provides methods and operators definitions for DEC.
+
+ * @tparam T An implementation that also inherit `SurfaceDEC`
+ @tparam TLinearAlgebraBackend linear algebra backend used (i.e.
+ EigenSparseLinearAlgebraBackend).
+ */
+ template
+ class SurfaceDEC
+ {
+ typedef typename TLinearAlgebraBackend::DenseVector::Index Index;
+ typedef typename TLinearAlgebraBackend::DenseMatrix DenseMatrix;
+ typedef typename TLinearAlgebraBackend::SparseMatrix LinearOperator;
+
+ public:
+ /// @name Global operators
+ /// @{
+
+ /// \brief Global differential
+ /// @return the n_e x n_v diffential matrix
+ LinearOperator D0() const
+ {
+ return static_cast( this )->buildD0();
+ }
+
+ /// \brief Global Laplace-Beltrami
+ /// @return the n_v x n_v striffness matrix
+ LinearOperator L0() const
+ {
+ return static_cast( this )->buildL0();
+ }
+
+ /// \brief Global inner product between 0 forms
+ ///
+ /// For a diagonal version see lumpedM0()
+ /// @return the n_v x n_v mass matrix
+ LinearOperator M0() const
+ {
+ return static_cast( this )->buildM0();
+ }
+
+ /// \brief Global diagonal inner product between 0 forms
+ ///
+ /// @param f a face
+ /// @return the n_v x n_v diagonal mass matrix
+ LinearOperator lumpedM0() const
+ {
+ return static_cast( this )->buildLumpedM0();
+ }
+
+ /// \brief Global inner product between 1 forms
+ /// @return the n_e x n_e mass matrix
+ LinearOperator M1() const
+ {
+ return static_cast( this )->buildM1();
+ }
+
+ /// \brief Global inner product between 2 forms
+ /// @return the n_f x n_f mass matrix
+ LinearOperator M2() const
+ {
+ return static_cast( this )->buildM2();
+ }
+
+ /// \brief Global Sharp operator
+ /// @return a 3*n_f x n_e matrix
+ LinearOperator Sharp() const
+ {
+ return static_cast( this )->buildSharp();
+ }
+
+ /// \brief Global Flat operator
+ /// @return a n_e x 3*n_f matrix
+ LinearOperator Flat() const
+ {
+ return static_cast( this )->buildFlat();
+ }
+
+ /// @}
+
+ /// @name Per face operators
+ /// @{
+
+ /// \brief Differential inside a face
+ /// @param f a face
+ /// @return the degree x degree differential matrix
+ DenseMatrix localD0( Index f ) const
+ {
+ return static_cast( this )->buildLocalD0( f );
+ }
+
+ /// \brief Laplace-Beltrami inside a face
+ /// @param f a face
+ /// @return the degree x degree stiffness matrix
+ DenseMatrix localL0( Index f ) const
+ {
+ return static_cast( this )->buildLocalL0( f );
+ }
+
+ /// \brief Inner product between 0 forms inside a face
+ /// @param f a face
+ /// @return the degree x degree mass matrix
+ DenseMatrix localM0( Index f ) const
+ {
+ return static_cast( this )->buildLocalM0( f );
+ }
+
+ /// \brief Inner product between 1 forms inside a face
+ /// @param f a face
+ /// @return the degree x degree mass matrix
+ DenseMatrix localM1( Index f ) const
+ {
+ return static_cast( this )->buildLocalM1( f );
+ }
+
+ /// \brief Inner product between 2 forms inside a face
+ /// @param f a face
+ /// @return the 1 x 1 mass matrix
+ DenseMatrix localM2( Index f ) const
+ {
+ return static_cast( this )->buildLocalM2( f );
+ }
+
+ /// \brief Sharp operator for inside a face
+ /// @param f a face
+ /// @return a 3 x degree matrix
+ DenseMatrix localSharp( Index f ) const
+ {
+ return static_cast( this )->buildLocalSharp( f );
+ }
+
+ /// \brief Flat operator inside a face
+ /// @param f a face
+ /// @return a degree x 3 matrix
+ DenseMatrix localFlat( Index f ) const
+ {
+ return static_cast( this )->buildLocalFlat( f );
+ }
+
+ /// @}
+ };
+
+} // namespace DGtal
+
+#endif
diff --git a/src/DGtal/dec/doc/images/cc/poisson-cc-g.png b/src/DGtal/dec/doc/images/cc/poisson-cc-g.png
new file mode 100644
index 0000000000..5508269702
Binary files /dev/null and b/src/DGtal/dec/doc/images/cc/poisson-cc-g.png differ
diff --git a/src/DGtal/dec/doc/images/cc/poisson-cc-surf.png b/src/DGtal/dec/doc/images/cc/poisson-cc-surf.png
new file mode 100644
index 0000000000..4ae40abf0a
Binary files /dev/null and b/src/DGtal/dec/doc/images/cc/poisson-cc-surf.png differ
diff --git a/src/DGtal/dec/doc/images/cc/poisson-cc-u.png b/src/DGtal/dec/doc/images/cc/poisson-cc-u.png
new file mode 100644
index 0000000000..7c0bf50b02
Binary files /dev/null and b/src/DGtal/dec/doc/images/cc/poisson-cc-u.png differ
diff --git a/src/DGtal/dec/doc/images/fem/poisson-fem-g.png b/src/DGtal/dec/doc/images/fem/poisson-fem-g.png
new file mode 100644
index 0000000000..e00722508b
Binary files /dev/null and b/src/DGtal/dec/doc/images/fem/poisson-fem-g.png differ
diff --git a/src/DGtal/dec/doc/images/fem/poisson-fem-surf.png b/src/DGtal/dec/doc/images/fem/poisson-fem-surf.png
new file mode 100644
index 0000000000..3a7a613177
Binary files /dev/null and b/src/DGtal/dec/doc/images/fem/poisson-fem-surf.png differ
diff --git a/src/DGtal/dec/doc/images/fem/poisson-fem-u.png b/src/DGtal/dec/doc/images/fem/poisson-fem-u.png
new file mode 100644
index 0000000000..5466a4a860
Binary files /dev/null and b/src/DGtal/dec/doc/images/fem/poisson-fem-u.png differ
diff --git a/src/DGtal/dec/doc/moduleInterpolatedCorrectedCalculus.dox b/src/DGtal/dec/doc/moduleInterpolatedCorrectedCalculus.dox
new file mode 100644
index 0000000000..f0608d9104
--- /dev/null
+++ b/src/DGtal/dec/doc/moduleInterpolatedCorrectedCalculus.dox
@@ -0,0 +1,192 @@
+/**
+ * @file
+ * @author Colin Weill--Duflos (\c colin.weill-duflos@univ-smb.fr )
+ * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
+ *
+ * @date 2024/06/05
+ *
+ * Documentation file for feature InterpolatedCorrectedCalculus
+ *
+ * This file is part of the DGtal library.
+ */
+
+/*
+ * Useful to avoid writing DGtal:: in front of every class.
+ * Do not forget to add an entry in src/DGtal/base/Config.h.in !
+ */
+namespace DGtal {
+//----------------------------------------
+/*!
+@page moduleInterpolatedCorrectedCalculus Interpolated Normals corrected DEC on digital surfaces
+@writers Colin Weill--Duflos
+
+[TOC]
+
+@since 1.5
+
+ Part of package \ref packageDEC.
+
+In this documentation page, we detail the operators and tools for differential
+ calculus computations on digital surfaces equipped with a normal field at vertices.
+
+The method use for building operators is similar to \ref modulePolygonalCalculus ,
+except that it assumes that the surface it is built on is a digital surface and
+that it has corrected normal on vertices that are then interpolated when computing operators.
+
+@note The sign convention for the divergence and the Laplacian
+operator is opposite to the one of @cite degoes2020discrete. 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 \f$ \langle \mathrm{d} u, v \rangle = - \langle u,
+\mathrm{div} v \rangle \f$. See also
+https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+
+@note All illustrations below have been obtained using the
+DGtal+[polyscope](https://polyscope.run) example \ref dgtalCC-poisson.cpp. To build these examples, enable the `BUILD_POLYSCOPE_EXAMPLES` variable (e.g. `cmake .. -DBUILD_POLYSCOPE_EXAMPLES=ON`).
+
+@warning 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.
+
+
+
+\section sectInterpolatedCorrectedCalcIntro Introduction
+
+We provide the same operators as those described in \cite degoes2020discrete . All of them are accessible as
+ global and local operators through the implementation of SurfaceDEC. See \ref moduleDECIntroduction for a more thorough introduction to DEC concepts.
+
+The main difference compared to other provided DEC implementation is the fact that normals are not considered constant through the faces,
+but as defined at vertices and the bilinearly interpolated during calculations. One of the key difference of this approach
+ is that the corrected geometry is continuous at edges.
+
+\section sectInterpolatedCorrectedCalcFace Operators
+
+\subsection sectInterpolatedCorrectedbuilding Building the object
+
+
+Let us consider a single face digital surface (i.e. : a square). Using :
+
+@code
+using CC = InterpolatedCorrectedCalculus;
+using SurfMesh = SurfaceMesh;
+using namespace Z3i;
+
+//Vertices
+std::vector positions={ {0,0,0},{1,0,0},{1,1,0},{0,0,1} };
+//Single face
+std::vector> faces={{ 0,1,2,3 }};
+
+mesh = SurfMesh(positions.begin(),positions.end(),faces.begin(),faces.end());
+// We require normal at vertices. For digital surfaces, prefer using a convergent estimator
+mesh.computeFaceNormalsFromPositions();
+mesh.computeVertexNormalsFromFaceNormals();
+
+CC calculus(mesh);
+@endcode
+
+we obtain a SurfaceMesh instance with a unique face, we add the naturally defined normal
+to it and we define its associated InterpolatedCorrectedCalculus object.
+
+\subsection sectInterpolatedCorrectedLambda Second order normal interpolation, regularization parameter
+
+Note that the constructor takes two optional parameters : wether or not to use second order normal interpolation (defaults to false),
+and a regularization parameter lambda (defaults to 0.1).
+
+Second normal interpolation uses average values of normals, normalized, at edge midpoints to use a second order
+normal interpolation instead of the bilinear interpolation used in computations. This does not require more
+values of normals and thus the normal field used is not, in itself, more precise, but this avoid interpolated
+values of normal with a norm too far from 1.
+
+Lambda is the same regularization parameter that appears in
+\cite degoes2020discrete when building the product between 1 forms and the Lapalce-Beltrami operator.
+
+We can use second order interpolation as well as a different value for this lambda, here 0.5 :
+@code
+CC calculus(mesh, true, 0.5);
+@endcode
+
+\subsection sectInterpolatedCorrectedbuildingop Building the operators
+
+The following operators can be built, both using local (per face) or global variants. Here @f$ n_v@f$, @f$ n_e @f$ and @f$ n_f @f$
+denote the number of vertices, edges and faces of the surface and @f$ deg(f) @f$ the degree of face f.
+
+Operator | Local Variant | Local shape | Global variant | Global shape
+--|--|--|--|--
+D0 : differential for 0 forms | `calculus.localD0(f);` | @f$ deg(f) \times deg(f) @f$ | `calculus.D0()` | @f$ n_e \times n_v @f$
+M0 : inner product between 0 forms | `calculus.localM0(f);` | @f$ deg(f) \times deg(f) @f$ | `calculus.M0()` | @f$ n_v \times n_v @f$
+lumpedM0 : diagonal version of M0 (less accurate, easier to inverse) | | | `calculus.lumpedM0()`| @f$ n_v \times n_v @f$
+M1 : inner product between 1 forms | `calculus.localM1(f);` | @f$ deg(f) \times deg(f) @f$ | `calculus.M1()` | @f$ n_e \times n_e @f$
+M2 : inner product between 2 forms | `calculus.localM2(f);` | @f$ 1 \times 1 @f$ | `calculus.M2()` | @f$ n_f \times n_f @f$
+Sharp | `calculus.localSharp(f);` | @f$ 3 \times deg(f) @f$ | `calculus.Sharp()` | @f$ 3n_f \times n_e @f$
+Flat | `calculus.localFlat(f);` | @f$ deg(f) \times 3 @f$ | `calculus.Flat()` | @f$ n_e \times 3n_f @f$
+L0 : integrated Laplace-Beltrami | `calculus.localL0(f);` | @f$ deg(f) \times deg(f) @f$ | `calculus.L0()` | @f$ n_v \times n_v @f$
+
+Vectors fields are represented by one vector of size @f$ 3 n_f @f$ (for a field defined at faces), whith the first
+@f$n_f@f$ values containing the x coordinates, then the next @f$n_f@f$ the y coordinates and the last @f$n_f@f$ the z
+coordinates.
+
+Let us add a scalar function on vertices. For instance, using Eigen syntax, we can use;
+
+@code
+CC::DenseVector phi(4); //4 vertices
+phi << 1.0, 2.0, 0.0, 5.0;
+@endcode
+
+We can compute the gradient of this function using DEC formulas : @f$\nabla f = (d f)^\sharp @f$.
+
+Since we have a unique face, we can compute it locally or globally and obtain the same result.
+
+The local approach would be :
+@code
+CC::DenseVector localGrad = calculus.localSharp(0) * calculus.localD0(0) * phi;
+@endcode
+
+The global approach would be :
+@code
+CC::DenseVector grad = calculus.Sharp() * calculus.D0() * phi;
+@endcode
+
+In a similar manner, we can compute the (integrated) divergence of this gradient either locally or globally, using @f$\nabla . \mathbf(u) = *d*\mathbf(u)^\flat @f$.
+Since we use the integrated version (to havoid having to inverse M0), we only have to compute @f$ d*\mathbf(u)^\flat @f$
+
+Locally :
+@code
+CC::DenseVector localDiv = calculus.localD0(0).transpose() * calculus.localM1(0) * calculus.localFlat(0) * localGrad;
+@endcode
+
+Globally :
+@code
+CC::DenseVector div = calculus.D0().transpose() * calculus.M1() * calculus.Flat() * localGrad;
+@endcode
+
+\section sectInterpolatedCorrectedCalcPoisson Example: Solving a Laplace problem
+
+Let suppose we want to solve the following Laplace problem for data interpolation:
+\f{eqnarray*}{
+ \Delta_\Omega u& = 0 \\
+ & s.t. u = g \text{ on } \partial\Omega
+\f}
+
+We want to solve that problem on a digital surface @f$\Omega@f$
+ 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:
+\snippet dgtalCC-poisson.cpp CC-init
+
+Leading to the following results (see \ref dgtalCC-poisson.cpp):
+
+Surface | Boundary condition @f$ g@f$ | Solution @f$ u @f$
+--|--|--
+@image html images/cc/poisson-cc-surf.png "" | @image html images/cc/poisson-cc-g.png "" | @image html images/cc/poisson-cc-u.png ""
+
+*/
+
+}
diff --git a/src/DGtal/dec/doc/moduleNormalCorrectedFEM.dox b/src/DGtal/dec/doc/moduleNormalCorrectedFEM.dox
new file mode 100644
index 0000000000..06c3ea6084
--- /dev/null
+++ b/src/DGtal/dec/doc/moduleNormalCorrectedFEM.dox
@@ -0,0 +1,153 @@
+
+/**
+ * @file
+ * @author Colin Weill--Duflos (\c colin.weill-duflos@univ-smb.fr )
+ * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
+ *
+ * @date 2024/06/05
+ *
+ * Documentation file for feature NormalCorrectedFEM
+ *
+ * This file is part of the DGtal library.
+ */
+
+/*
+ * Useful to avoid writing DGtal:: in front of every class.
+ * Do not forget to add an entry in src/DGtal/base/Config.h.in !
+ */
+namespace DGtal {
+//----------------------------------------
+/*!
+@page moduleNormalCorrectedFEM Normal corrected Finite Element Method for Poisson problems
+@writers Colin Weill--Duflos
+
+[TOC]
+
+@since 1.5
+
+ Part of package \ref packageDEC.
+
+In this documentation page, we detail a laplacian operator
+ computations on surface mesh equipped with faces normals provided by the
+ NormalCorrectedFEM class.
+
+@note The sign convention for the divergence and the Laplacian
+operator is opposite to the one of @cite degoes2020discrete. 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 \f$ \langle \mathrm{d} u, v \rangle = - \langle u,
+\mathrm{div} v \rangle \f$. See also
+https://en.wikipedia.org/wiki/Laplace–Beltrami_operator
+
+@note All illustrations below have been obtained using the
+DGtal+[polyscope](https://polyscope.run) example dgtalFEM-poisson.cpp. To build these examples, enable the `BUILD_POLYSCOPE_EXAMPLES` variable (e.g. `cmake .. -DBUILD_POLYSCOPE_EXAMPLES=ON`).
+
+@warning 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.
+
+
+
+\section sectNormalCorrectedFEMIntro Introduction
+
+We aim at solving a problem of the form @f$\Delta u = f@f$.
+We follows the Finite Element Method to derive two matrices @f$L, M@f$
+such that @f$L \mathbf{u} = M \mathbf{f}@f$. We use a normal based
+metric to make our operator work on digital surfaces.
+
+The metric used is @f$ G = \begin{bmatrix} 1 - (\mathbf{u}_x)^2 & -\mathbf{u}_x \mathbf{u}_y \\
+-\mathbf{u}_x \mathbf{u}_y & 1-(\mathbf{u}_y)^2 \end{bmatrix} @f$
+
+Where @f$\mathbf{u} = (\mathbf{u}_x, \mathbf{u}_y, \mathbf{u}_z)@f$ is the provided normal
+for a surfel in the natural coordinate surface of the surfel.
+
+\section sectNormalCorrectedFEMface Per face operators
+
+Let us consider a single face digital surface (i.e. : a square). Using :
+
+@code
+using ncFEM = NormalCorrectedFEM;
+using SurfMesh = SurfaceMesh;
+using namespace Z3i;
+
+//Vertices
+std::vector positions={ {0,0,0},{1,0,0},{1,1,0},{0,0,1} };
+//Single face
+std::vector> faces={{ 0,1,2,3 }};
+
+mesh = SurfMesh(positions.begin(),positions.end(),faces.begin(),faces.end());
+// We require normal at faces. For digital surfaces, prefer using a convergent estimator
+mesh.computeFaceNormalsFromPositions();
+
+ncFEM calculus(mesh);
+@endcode
+
+we obtain a SurfaceMesh instance with a unique face, we add the naturally defined normal
+to it and we define its associated NormalCorrectedFEM object.
+
+Let us add a scalar function on vertices. For instance, using Eigen syntax, we can use;
+
+@code
+ncFEM::DenseVector phi(4); //4 vertices
+phi << 1.0, 2.0, 0.0, 5.0;
+@endcode
+
+\subsection sectNormalCorrectedMMstiff Mass matrix and stiffness matrix
+
+The expression for the integral of the product of two functions @f$f@f$ and @f$g@f$ inside a surfel is given by the following expression :
+@f[ \int_{\square} \sqrt{det(G)}fg @f]
+
+The expression for the integral of the dot product of the gradient of two functions @f$f@f$ and @f$g@f$ inside a surfel is given by the following expression :
+@f[ \int_{\square} \sqrt{det(G)} \nabla f^T G^{-1} \nabla g \rangle @f]
+
+By evaluating this expression with a linear basis, we can build matrices letting us evaluate these integrals for any two linear functions inside a face.
+The results can be obtained in the following way :
+
+Operator | Output | Description
+------------- | ------------- | ----
+`calculus.localM0(f)` | \f$ n_v\times n_v\f$ | mass matrix (corresponds to the inner product between 0 forms in @cite degoes2020discrete)
+`calculus.localL0(f)`| \f$ n_v \times n_v\f$ | stiffness matrix (corresponds to Laplace-Beltrami operator in @cite degoes2020discrete), the matrix is PSD
+
+\section sectNormalCorrectedFEMGlo Global operators
+
+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.
+
+
+Operator | Output | Description
+------------- | ------------- | ----
+`calculus.M0()` | \f$ n_v \times n_v\f$ | global mass matrix (PD)
+`calculus.lumpedM0()` | \f$ n_v \times n_v\f$ | global lumped mass matrix (Diagonal, positive)
+`calculus.L0()`| \f$ n_v \times n_v\f$ | global stiffness matrix (SPD)
+
+
+\section secNormalCorrectedFEMLap Example: Solving a Laplace problem
+
+Let suppose we want to solve the following Laplace problem for data interpolation:
+\f{eqnarray*}{
+ \Delta_\Omega u& = 0 \\
+ & s.t. u = g \text{ on } \partial\Omega
+\f}
+
+We want to solve that problem on a digital surface @f$\Omega@f$
+ 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:
+\snippet dgtalFEM-poisson.cpp FEM-init
+
+Leading to the following results (see \ref dgtalFEM-poisson.cpp):
+
+Surface | Boundary condition @f$ g@f$ | Solution @f$ u @f$
+--|--|--
+@image html images/fem/poisson-fem-surf.png "" | @image html images/fem/poisson-fem-g.png "" | @image html images/fem/poisson-fem-u.png ""
+
+*/
+
+}
diff --git a/src/DGtal/dec/doc/modulePolygonalCalculus.dox b/src/DGtal/dec/doc/modulePolygonalCalculus.dox
index 643ee65900..cf34205ed7 100644
--- a/src/DGtal/dec/doc/modulePolygonalCalculus.dox
+++ b/src/DGtal/dec/doc/modulePolygonalCalculus.dox
@@ -1,5 +1,5 @@
/**
- * @file
+ * @file
* @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
* Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
*
@@ -52,8 +52,8 @@ DGtal+[polyscope](https://polyscope.run) examples \ref
dgtalCalculus.cpp, \ref dgtalCalculus-single.cpp and \ref dgtalCalculus-poisson.cpp. To build these examples, enable the `BUILD_POLYSCOPE_EXAMPLES` variable (e.g. `cmake .. -DBUILD_POLYSCOPE_EXAMPLES=ON`).
@warning 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.
-
-
+
+
\section sectPolygonalCalculusIntro Introduction
@@ -84,13 +84,13 @@ std::vector positions={ {0,0,0},{20,0,0},{20,10,0},{10,8,5}, {0,15,1}
std::vector> faces={{ 0,1,2,3,4 }};
mesh = SurfMesh(positions.begin(),positions.end(),faces.begin(),faces.end());
-
+
PolygonalCalculus calculus(mesh);
@endcode
we obtain a SurfaceMesh instance with a unique face and its associated PolygonalCalculus object.
- Example| Example
+ Example| Example
--|--
@image html images/poly/face.png "" | @image html images/poly/face_rot.png ""
@@ -114,17 +114,17 @@ default, the positions of the surface mesh vertices (@a positions @a 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::setEmbedder for an example.
+PolygonalCalculus::setEmbedder for an example.
@note As the Face id is a parameter of the embedder, a given vertex
can have different embeddings for all its incident faces.
-
+
\subsection susub1 Basic operators
-
-
-
+
+
+
We first describe some standard per face operators. Note that for all
extrinsic operators that require the vertex position in @f$
@@ -148,7 +148,7 @@ Operator | Output | Description
\subsubsection subderiv Derivative operators
-Derivative operators act on a scalar field defined on the vertex of the face.
+Derivative operators act on a scalar field defined on the vertex of the face.
Operator | Output | Description
@@ -160,7 +160,7 @@ Operator | Output | Description
A*phi | D*phi (discrete 1-form)
--| --
-@image html images/poly/A_phi.png "" | @image html images/poly/D_phi.png ""
+@image html images/poly/A_phi.png "" | @image html images/poly/D_phi.png ""
gradient | co-gradient | corrected normal
@@ -194,10 +194,10 @@ a vector in @f$\mathbb{R}^3@f$ | a vector in @f$\mathbb{R}^3@f$ | Flat (1-form)
-\subsection sublap Inner Product and Laplace-Beltrami Operators
+\subsection subdlap Inner Product and Laplace-Beltrami Operators
-For 0-forms, the inner product is the classical one induced by the @f$ l_2@f$ norm. For 1-forms, the inner product is given by PolygonalCalculus::M() (useful to define the Laplace-Beltrami operator).
+For 0-forms, the inner product is the classical one induced by the @f$ l_2@f$ norm. For 1-forms, the inner product is given by PolygonalCalculus::M() (useful to define the Laplace-Beltrami operator).
Operator | Output | Description
------------- | ------------- | ----
@@ -208,15 +208,15 @@ Operator | Output | Description
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 @cite degoes2020discrete, section 5, for details.
-
+
As an example, these operators can be used to interpolate vector fields as illustrated in @ref moduleVectorsInHeat.
-
+
@image html images/poly/transport_to_face.png "Transport of a vertex based vector field to face tangent plane"
\subsubsection sectcovOP 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 @cite degoes2020discrete,
-i.e.
+i.e.
\f{eqnarray*}{
u_{f} = [u_{v_1}^t \ldots u_{v_{nf}}^t]^t
\f}
@@ -232,7 +232,7 @@ Operator | Output | Description
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)` | \f$ 2n_f\times 2n_f\f$ | Vector Laplacian defined as the associated Matrix with the Dirichlet energy for vector valued 0-forms (lambda is a regularization parameter, see @cite degoes2020discrete) at face f, PSD matrix
+`calculus.connectionLaplacian(f,lambda)` | \f$ 2n_f\times 2n_f\f$ | Vector Laplacian defined as the associated Matrix with the Dirichlet energy for vector valued 0-forms (lambda is a regularization parameter, see @cite degoes2020discrete) at face f, PSD matrix
\section sectPolygonalCalculusGlo Global calculus
@@ -240,12 +240,12 @@ Operator | Output | Description
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 \ref dgtalCalculus.cpp example using Shortcuts and ShortcutsGeometry to set up the surface:
-Surface | Phi | Gradient | Gradient+co-gradient
+Surface | Phi | Gradient | Gradient+co-gradient
--|--|--|--
@image html images/poly/init.png "" | @image html images/poly/goursat_phi.png "" | @image html images/poly/goursat_grad.png "" | @image html images/poly/goursat_gradcograd.png ""
-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.
+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::globalLaplaceBeltrami() method outputs a global (sparse) Laplace-Beltrami operator which can later be used for diffusion.
@@ -257,9 +257,9 @@ used when one solves a weak problem and wishes to get a pointwise
per-vertex solution.
-
+
\section sectCorrected 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 @cite coeurjolly2022simple, given an estimation of the tangent bundle of the discrete surface (for instance using
estimated normal vectors from @ref moduleIntegralInvariant, or @ref moduleVCM, cf @ref moduleShortcuts),
@@ -269,16 +269,16 @@ Geodesic distances without correction | Geodesic distances with correction
--|--
@image html images/poly/corrected-without.png "" | @image html images/poly/corrected-with.png ""
-
+
The functor functors::EmbedderFromNormalVectors can be used to implicitly project Face vertices onto the prescribed tangent plane. A classical usage is the following one:
-
+
@code
//A surface mesh. Eg. primal surface of a digital surface
SurfaceMesh< Z3i::RealPoint, Z3i::RealVector > surfmesh(...);
//Per face normal vector estimation
std::vector ii_normals = ....
-
+
//New embedder using tangent plane projection of face vertices
functors::EmbedderFromNormalVectors embedderFromNormals(ii_normals,surfmesh);
@@ -287,58 +287,58 @@ The functor functors::EmbedderFromNormalVectors can be used to implicitly projec
calculus->setEmbedder( embedderFromNormals );
@endcode
A complete code is given in the @ref dgtalCalculus-geodesic.cpp example.
-
+
\section secLap Example: Solving a Laplace problem
-
+
Let suppose we want to solve the following Laplace problem for data interpolation:
\f{eqnarray*}{
\Delta_\Omega u& = 0 \\
& s.t. u = g \text{ on } \partial\Omega
\f}
-
+
We want to solve that problem on a polygonal mesh @f$\Omega@f$
(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:
\snippet dgtalCalculus-poisson.cpp PolyDEC-init
-
+
Leading to the following results (see \ref dgtalCalculus-poisson.cpp):
-
+
Surface | Boundary condition @f$ g@f$ | Solution @f$ u @f$
--|--|--
@image html images/poly/poisson-surf.png "" | @image html images/poly/poisson-g.png "" | @image html images/poly/poisson-u.png ""
@image html images/poly/bunny-init.png "" | @image html images/poly/bunny-g.png "" | @image html images/poly/bunny-u.png ""
@image html images/poly/cat-init.png "" | @image html images/poly/cat-g.png "" | @image html images/poly/cat-u.png ""
-
+
\subsection Global Vector Calculus
-
+
Global Vector Laplace/Poisson problems can also be solved by the same way, using instead PolygonalCalculus::globalConnectionLaplace() and PolygonalCalculus::doubledGlobalLumpedMassMatrix(). One can find examples of such use in the \ref VectorsInHeat class.
-
-
-
-
-
+
+
+
+
+
\section sectMisc Miscellaneous
\subsection sectPolygonalCalculusHP 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 @b 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.
@code
PolygonalCalculus calculus(surfmesh,true); //global internal cache enabled.
@endcode
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 \ref 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 @e M @e operator per face is given by
@code
DenseMatrix M(const Face f, const double lambda=1.0) const
@@ -347,12 +347,12 @@ DenseMatrix M(const Face f, const double lambda=1.0) const
auto Pf=P(f);
return faceArea(f) * Uf.transpose()*Uf + lambda * Pf.transpose()*Pf;
}
-@endcode
-which could be time consuming as the internal operators may be computed several times.
+@endcode
+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
+A typical usage is
@code
auto cacheSharp = getOperatorCacheMatrix( [&](Face f){ return(calculus.sharp(f);} );
auto cachefaceArea = getOperatorCacheMatrix( [&](Face f){ return(calculus.faceArea(f);} );
@@ -362,9 +362,9 @@ auto cacheP = getOperatorCacheMatrix( [&](Face f){ return(calculus.P(f);}
Then, cached operators can be accessed and combined:
@code
auto Mf = cachefaceArea[f] * cacheU[f].transpose()*cacheU[f] + lambda * cacheP[f].transpose() * cacheP[f];
-@endcode
+@endcode
+
-
*/
diff --git a/src/DGtal/dec/doc/packageDEC.dox b/src/DGtal/dec/doc/packageDEC.dox
index 63c88e8602..54236e07ef 100644
--- a/src/DGtal/dec/doc/packageDEC.dox
+++ b/src/DGtal/dec/doc/packageDEC.dox
@@ -59,10 +59,13 @@ Basic operators, such as Hodge duality operator or exterior derivative, can be c
- \subpage moduleHeatLaplaceOperator (Thomas Caissard)
- \subpage moduleGenericAT (Jacques-Olivier Lachaud, Marion Foare, David Coeurjolly, Pierre Gueth)
+
- Discrete Corrected Polygonal Calculus
- \subpage modulePolygonalCalculus (David Coeurjolly, Jacques-Olivier Lachaud, Baptiste Genest)
- \subpage moduleGeodesicsInHeat (David Coeurjolly, Jacques-Olivier Lachaud)
- \subpage moduleVectorsInHeat (Baptiste Genest, David Coeurjolly)
+ - \subpage moduleNormalCorrectedFEM (Colin Weill--Duflos)
+ - \subpage moduleInterpolatedCorrectedCalculus (Colin Weill--Duflos)
@b Package @b Concepts @b Overview
- \subpage packageDECConcepts
@@ -84,6 +87,8 @@ Basic operators, such as Hodge duality operator or exterior derivative, can be c
- exampleVectorHeatMethod.cpp
- exampleHarmonicParametrization.cpp
- exampleBunnyHead.cpp
+- dgtalFEM-poisson.cpp
+- dgtalCC-poisson.cpp
*/
diff --git a/tests/dec/CMakeLists.txt b/tests/dec/CMakeLists.txt
index ddeb51b911..cef4a4b324 100644
--- a/tests/dec/CMakeLists.txt
+++ b/tests/dec/CMakeLists.txt
@@ -1,10 +1,12 @@
set(DGTAL_TESTS_SRC
- testDiscreteExteriorCalculus
+ testDiscreteExteriorCalculus
testEmbedding
testHeatLaplace
testPolygonalCalculus
testGeodesicsInHeat
testVectorsInHeat
+ testInterpolatedCorrectedCalculus
+ testCorrectedFEM
)
# add_test is disabled for the following sources
@@ -21,4 +23,3 @@ if(WITH_EIGEN)
DGtal_add_test(${FILE} ONLY_ADD_EXECUTABLE)
endforeach()
endif()
-
diff --git a/tests/dec/testCorrectedFEM.cpp b/tests/dec/testCorrectedFEM.cpp
new file mode 100644
index 0000000000..c68ae8e4cb
--- /dev/null
+++ b/tests/dec/testCorrectedFEM.cpp
@@ -0,0 +1,185 @@
+
+///////////////////////////////////////////////////////////////////////////////
+#include
+#include "DGtal/base/Common.h"
+#include "ConfigTest.h"
+#include "DGtal/geometry/helpers/PlaneProbingEstimatorHelper.h"
+#include "DGtalCatch.h"
+#include "DGtal/helpers/StdDefs.h"
+
+#include "DGtal/dec/NormalCorrectedFEM.h"
+#include "DGtal/shapes/SurfaceMesh.h"
+#include "DGtal/shapes/MeshHelpers.h"
+#include "DGtal/helpers/Shortcuts.h"
+#include "DGtal/math/linalg/DirichletConditions.h"
+#include "DGtal/math/linalg/EigenSupport.h"
+///////////////////////////////////////////////////////////////////////////////
+
+using namespace std;
+using namespace DGtal;
+using namespace Z3i;
+
+///////////////////////////////////////////////////////////////////////////////
+// Functions for testing class InterpolatedCorrectedCalculus.
+///////////////////////////////////////////////////////////////////////////////
+
+TEST_CASE( "Testing PolygonalCalculus" )
+{
+ typedef SurfaceMesh< RealPoint,RealPoint > Mesh;
+ typedef NormalCorrectedFEM CFEM;
+ std::vector positions = { RealPoint( 0, 0, 0 ) ,
+ RealPoint( 1, 0, 0 ) ,
+ RealPoint( 0, 1, 0 ) ,
+ RealPoint( 1, 1, 0 ) ,
+ RealPoint( 0, 0, 1 ) ,
+ RealPoint( 1, 0, 1 ) ,
+ RealPoint( 0, 1, 1 ) ,
+ RealPoint( 1, 1, 1 ) ,
+ RealPoint( 1, 0, 2 ) ,
+ RealPoint( 0, 0, 2 ) };
+ std::vector faces = { { 1, 0, 2, 3 },
+ { 0, 1, 5, 4 } ,
+ { 1, 3, 7, 5 } ,
+ { 3, 2, 6, 7 } ,
+ { 2, 0, 4, 6 } ,
+ { 4, 5, 8, 9 } };
+
+ Mesh box(positions.cbegin(), positions.cend(),
+ faces.cbegin(), faces.cend());
+
+ box.computeFaceNormalsFromPositions();
+ CFEM boxCalculus(box);
+ SECTION("Local Laplace-Beltrami")
+ {
+ CFEM::Face f = 0;
+ auto nf = box.incidentVertices(f).size();
+
+ auto L = boxCalculus.localL0(f);
+ CFEM::LinearAlgebraBackend::DenseVector phi(nf),expected(nf);
+ phi << 1.0, 1.0, 1.0, 1.0;
+ expected << 0,0,0,0;
+ auto lphi = L*phi;
+ for(int i = 0; i < nf; i++) {
+ REQUIRE(abs(lphi(i)) < 1e-15);
+ }
+ }
+
+ /*
+ SECTION("Check lumped mass matrix")
+ {
+ PolygonalCalculus< RealPoint,RealVector >::SparseMatrix M = boxCalculus.globalLumpedMassMatrix();
+ double a=0.0;
+ for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index v=0; v < box.nbVertices(); ++v )
+ a += M.coeffRef(v,v);
+
+ double fa=0.0;
+ for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index f=0; f < box.nbFaces(); ++f )
+ fa += box.faceArea(f);
+ REQUIRE( a == fa );
+ }
+ */
+}
+
+TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
+{
+ typedef Shortcuts< KSpace > SH3;
+ typedef SurfaceMesh< RealPoint,RealPoint > Mesh;
+ typedef Mesh::Index Index;
+ typedef NormalCorrectedFEM CFEM;
+ typedef CFEM::LinearAlgebraBackend::DenseVector DenseVector;
+ typedef DirichletConditions< EigenLinearAlgebraBackend > DC;
+
+ // Build a more complex surface.
+ auto params = SH3::defaultParameters();
+
+ params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
+ auto implicit_shape = SH3::makeImplicitShape3D ( params );
+ auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
+ auto K = SH3::getKSpace( params );
+ auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
+ auto surface = SH3::makeDigitalSurface( binary_image, K, params );
+ auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
+
+ std::vector > faces;
+ std::vector positions = primalSurface->positions();
+ for( Index face= 0 ; face < primalSurface->nbFaces(); ++face)
+ faces.push_back(primalSurface->incidentVertices( face ));
+
+ Mesh surfmesh = Mesh( positions.begin(), positions.end(),
+ faces.begin(), faces.end() );
+ auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
+ surfmesh.computeFaceNormalsFromPositions();
+
+ SECTION("Check surface")
+ {
+ REQUIRE( surfmesh.nbVertices() == 1364 );
+ REQUIRE( surfmesh.nbFaces() == 1279 );
+ REQUIRE( boundaryEdges.size() == 168 );
+ }
+
+ // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
+ CFEM calculus( surfmesh );
+ // Laplace opeartor
+ CFEM::LinearOperator L = calculus.L0();
+ // value on boundary
+ CFEM::LinearAlgebraBackend::DenseVector g = CFEM::LinearAlgebraBackend::DenseVector::Zero(surfmesh.nbVertices());
+ // characteristic set of boundary
+ DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
+
+ SECTION("Solve Poisson problem with boundary Dirichlet conditions")
+ {
+ for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
+ {
+ std::cout << "scale=" << scale << std::endl;
+ auto phi = [&]( Index v)
+ {
+ return cos(scale*(surfmesh.position(v)[0]))
+ * (scale*surfmesh.position(v)[1]);
+ };
+
+ for(auto &e: boundaryEdges)
+ {
+ auto adjVertices = surfmesh.edgeVertices(e);
+ auto v1 = adjVertices.first;
+ auto v2 = adjVertices.second;
+ g(v1) = phi(v1);
+ g(v2) = phi(v2);
+ b(v1) = 1;
+ b(v2) = 1;
+ }
+ // Solve Δu=0 with g as boundary conditions
+ EigenLinearAlgebraBackend::SolverSimplicialLDLT solver;
+ CFEM::LinearOperator L_dirichlet = DC::dirichletOperator( L, b );
+ solver.compute( L_dirichlet );
+ REQUIRE( solver.info() == Eigen::Success );
+ DenseVector g_dirichlet = DC::dirichletVector( L, g, b, g );
+ DenseVector x_dirichlet = solver.solve( g_dirichlet );
+ REQUIRE( solver.info() == Eigen::Success );
+ DenseVector u = DC::dirichletSolution( x_dirichlet, b, g );
+ double min_phi = 0.0;
+ double max_phi = 0.0;
+ double min_u = 0.0;
+ double max_u = 0.0;
+ double min_i_u = 0.0;
+ double max_i_u = 0.0;
+ for ( Index v = 0; v < surfmesh.nbVertices(); ++v )
+ {
+ min_phi = std::min( min_phi, phi( v ) );
+ max_phi = std::max( max_phi, phi( v ) );
+ min_u = std::min( min_u , u ( v ) );
+ max_u = std::max( max_u , u ( v ) );
+ if ( b( v ) == 0.0 )
+ {
+ min_i_u = std::min( min_i_u, u ( v ) );
+ max_i_u = std::max( max_i_u, u ( v ) );
+ }
+ }
+ REQUIRE( min_phi <= min_u );
+ REQUIRE( max_phi >= max_u );
+ REQUIRE( min_phi < min_i_u );
+ REQUIRE( max_phi > max_i_u );
+ } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
+ }
+};
+
+/** @ingroup Tests **/
diff --git a/tests/dec/testInterpolatedCorrectedCalculus.cpp b/tests/dec/testInterpolatedCorrectedCalculus.cpp
new file mode 100644
index 0000000000..01ae077f8d
--- /dev/null
+++ b/tests/dec/testInterpolatedCorrectedCalculus.cpp
@@ -0,0 +1,229 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+
+/**
+ * @file testPolygonalCalculus.cpp
+ * @ingroup Tests
+ * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
+ * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
+ *
+ * @date 2021/09/02
+ *
+ * Functions for testing class PolygonalCalculus.
+ *
+ * This file is part of the DGtal library.
+ */
+
+///////////////////////////////////////////////////////////////////////////////
+#include
+#include "DGtal/base/Common.h"
+#include "ConfigTest.h"
+#include "DGtalCatch.h"
+#include "DGtal/helpers/StdDefs.h"
+
+#include "DGtal/dec/InterpolatedCorrectedCalculus.h"
+#include "DGtal/shapes/SurfaceMesh.h"
+#include "DGtal/shapes/MeshHelpers.h"
+#include "DGtal/helpers/Shortcuts.h"
+#include "DGtal/math/linalg/DirichletConditions.h"
+#include "DGtal/math/linalg/EigenSupport.h"
+///////////////////////////////////////////////////////////////////////////////
+
+using namespace std;
+using namespace DGtal;
+using namespace Z3i;
+
+///////////////////////////////////////////////////////////////////////////////
+// Functions for testing class InterpolatedCorrectedCalculus.
+///////////////////////////////////////////////////////////////////////////////
+
+TEST_CASE( "Testing PolygonalCalculus" )
+{
+ typedef SurfaceMesh< RealPoint,RealPoint > Mesh;
+ typedef InterpolatedCorrectedCalculus CC;
+ std::vector positions = { RealPoint( 0, 0, 0 ) ,
+ RealPoint( 1, 0, 0 ) ,
+ RealPoint( 0, 1, 0 ) ,
+ RealPoint( 1, 1, 0 ) ,
+ RealPoint( 0, 0, 1 ) ,
+ RealPoint( 1, 0, 1 ) ,
+ RealPoint( 0, 1, 1 ) ,
+ RealPoint( 1, 1, 1 ) ,
+ RealPoint( 1, 0, 2 ) ,
+ RealPoint( 0, 0, 2 ) };
+ std::vector faces = { { 1, 0, 2, 3 },
+ { 0, 1, 5, 4 } ,
+ { 1, 3, 7, 5 } ,
+ { 3, 2, 6, 7 } ,
+ { 2, 0, 4, 6 } ,
+ { 4, 5, 8, 9 } };
+
+ Mesh box(positions.cbegin(), positions.cend(),
+ faces.cbegin(), faces.cend());
+
+ box.computeFaceNormalsFromPositions();
+ box.computeVertexNormalsFromFaceNormals();
+ CC boxCalculus(box);
+
+ SECTION("Derivatives")
+ {
+ CC::Face f = 0;
+ auto d = boxCalculus.localD0(f);
+
+ auto nf = box.incidentVertices(f).size();
+ CC::DenseVector phi(nf),expected(nf);
+ phi << 1.0, 3.0, 2.0, 6.0;
+ expected << 2,-1,4,-5;
+ auto dphi = d*phi; // n_f x 1 matrix
+ REQUIRE(dphi == expected);
+
+ }
+
+ SECTION("Local Laplace-Beltrami")
+ {
+ CC::Face f = 0;
+ auto nf = box.incidentVertices(f).size();
+
+ auto L = boxCalculus.localL0(f);
+ CC::DenseVector phi(nf),expected(nf);
+ phi << 1.0, 1.0, 1.0, 1.0;
+ expected << 0,0,0,0;
+ auto lphi = L*phi;
+ for(int i = 0; i < nf; i++) {
+ REQUIRE(abs(lphi(i)) < 1e-15);
+ }
+ }
+ /*
+ SECTION("Check lumped mass matrix")
+ {
+ PolygonalCalculus< RealPoint,RealVector >::SparseMatrix M = boxCalculus.globalLumpedMassMatrix();
+ double a=0.0;
+ for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index v=0; v < box.nbVertices(); ++v )
+ a += M.coeffRef(v,v);
+
+ double fa=0.0;
+ for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index f=0; f < box.nbFaces(); ++f )
+ fa += box.faceArea(f);
+ REQUIRE( a == fa );
+ }
+ */
+}
+
+TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
+{
+ typedef Shortcuts< KSpace > SH3;
+ typedef SurfaceMesh< RealPoint,RealPoint > Mesh;
+ typedef Mesh::Index Index;
+ typedef InterpolatedCorrectedCalculus CC;
+ typedef DirichletConditions< EigenLinearAlgebraBackend > DC;
+
+ // Build a more complex surface.
+ auto params = SH3::defaultParameters();
+
+ params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
+ auto implicit_shape = SH3::makeImplicitShape3D ( params );
+ auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
+ auto K = SH3::getKSpace( params );
+ auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
+ auto surface = SH3::makeDigitalSurface( binary_image, K, params );
+ auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
+
+ std::vector > faces;
+ std::vector positions = primalSurface->positions();
+ for( Index face= 0 ; face < primalSurface->nbFaces(); ++face)
+ faces.push_back(primalSurface->incidentVertices( face ));
+
+ Mesh surfmesh = Mesh( positions.begin(), positions.end(),
+ faces.begin(), faces.end() );
+ auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
+ surfmesh.computeFaceNormalsFromPositions();
+ surfmesh.computeVertexNormalsFromFaceNormals();
+
+ SECTION("Check surface")
+ {
+ REQUIRE( surfmesh.nbVertices() == 1364 );
+ REQUIRE( surfmesh.nbFaces() == 1279 );
+ REQUIRE( boundaryEdges.size() == 168 );
+ }
+
+ // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
+ CC calculus( surfmesh );
+ // Laplace opeartor
+ CC::LinearOperator L = calculus.L0();
+ // value on boundary
+ CC::DenseVector g = CC::DenseVector::Zero(surfmesh.nbVertices());
+ // characteristic set of boundary
+ DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
+
+ SECTION("Solve Poisson problem with boundary Dirichlet conditions")
+ {
+ for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
+ {
+ std::cout << "scale=" << scale << std::endl;
+ auto phi = [&]( Index v)
+ {
+ return cos(scale*(surfmesh.position(v)[0]))
+ * (scale*surfmesh.position(v)[1]);
+ };
+
+ for(auto &e: boundaryEdges)
+ {
+ auto adjVertices = surfmesh.edgeVertices(e);
+ auto v1 = adjVertices.first;
+ auto v2 = adjVertices.second;
+ g(v1) = phi(v1);
+ g(v2) = phi(v2);
+ b(v1) = 1;
+ b(v2) = 1;
+ }
+ // Solve Δu=0 with g as boundary conditions
+ EigenLinearAlgebraBackend::SolverSimplicialLDLT solver;
+ CC::LinearOperator L_dirichlet = DC::dirichletOperator( L, b );
+ solver.compute( L_dirichlet );
+ REQUIRE( solver.info() == Eigen::Success );
+ CC::DenseVector g_dirichlet = DC::dirichletVector( L, g, b, g );
+ CC::DenseVector x_dirichlet = solver.solve( g_dirichlet );
+ REQUIRE( solver.info() == Eigen::Success );
+ CC::DenseVector u = DC::dirichletSolution( x_dirichlet, b, g );
+ double min_phi = 0.0;
+ double max_phi = 0.0;
+ double min_u = 0.0;
+ double max_u = 0.0;
+ double min_i_u = 0.0;
+ double max_i_u = 0.0;
+ for ( Index v = 0; v < surfmesh.nbVertices(); ++v )
+ {
+ min_phi = std::min( min_phi, phi( v ) );
+ max_phi = std::max( max_phi, phi( v ) );
+ min_u = std::min( min_u , u ( v ) );
+ max_u = std::max( max_u , u ( v ) );
+ if ( b( v ) == 0.0 )
+ {
+ min_i_u = std::min( min_i_u, u ( v ) );
+ max_i_u = std::max( max_i_u, u ( v ) );
+ }
+ }
+ REQUIRE( min_phi <= min_u );
+ REQUIRE( max_phi >= max_u );
+ REQUIRE( min_phi < min_i_u );
+ REQUIRE( max_phi > max_i_u );
+ } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
+ }
+};
+
+/** @ingroup Tests **/
diff --git a/tests/dec/testPolygonalCalculus.cpp b/tests/dec/testPolygonalCalculus.cpp
index c0a4adf149..5ba78c6abb 100644
--- a/tests/dec/testPolygonalCalculus.cpp
+++ b/tests/dec/testPolygonalCalculus.cpp
@@ -75,12 +75,12 @@ TEST_CASE( "Testing PolygonalCalculus" )
{ 3, 2, 6, 7 } ,
{ 2, 0, 4, 6 } ,
{ 4, 5, 8, 9 } };
-
+
Mesh box(positions.cbegin(), positions.cend(),
faces.cbegin(), faces.cend());
-
+
PolygonalCalculus< RealPoint,RealVector > boxCalculus(box);
-
+
SECTION("Construction and basic operators")
{
REQUIRE( boxCalculus.isValid() );
@@ -90,7 +90,7 @@ TEST_CASE( "Testing PolygonalCalculus" )
auto x = boxCalculus.X(f);
auto d = boxCalculus.D(f);
auto a = boxCalculus.A(f);
-
+
//Checking X
PolygonalCalculus< RealPoint,RealVector >::Vector vec = x.row(0);
REQUIRE( vecToRealPoint(vec ) == positions[1]);
@@ -102,7 +102,7 @@ TEST_CASE( "Testing PolygonalCalculus" )
REQUIRE( vecToRealPoint(vec ) == positions[3]);
trace.info()<< boxCalculus <::Face f = 0;
auto d = boxCalculus.D(f);
-
+
auto nf = boxCalculus.faceDegree(f);
PolygonalCalculus< RealPoint,RealVector >::Vector phi(nf),expected(nf);
phi << 1.0, 3.0, 2.0, 6.0;
expected << 2,-1,4,-5;
auto dphi = d*phi; // n_f x 1 matrix
REQUIRE(dphi == expected);
-
+
}
-
+
SECTION("Structural propertes")
{
PolygonalCalculus< RealPoint,RealVector >::Face f = 0;
auto nf = boxCalculus.faceDegree(f);
PolygonalCalculus< RealPoint,RealVector >::Vector phi(nf);
phi << 1.0, 3.0, 2.0, 6.0;
-
+
auto G = boxCalculus.gradient(f);
auto gphi = G*phi;
auto coG = boxCalculus.coGradient(f);
auto cogphi = coG*phi;
-
+
// grad . cograd == 0
REQUIRE( gphi.dot(cogphi) == 0.0);
-
+
// Gf = Uf Df
REQUIRE( G == boxCalculus.sharp(f)*boxCalculus.D(f));
-
+
// UV = I - nn^t (lemma4)
PolygonalCalculus< RealPoint,RealVector >::Vector n = boxCalculus.faceNormal(f);
REQUIRE( boxCalculus.sharp(f)*boxCalculus.flat(f) == PolygonalCalculus< RealPoint,RealVector >::DenseMatrix::Identity(3,3) - n*n.transpose() );
-
+
// P^2 = P (lemma6)
auto P = boxCalculus.P(f);
REQUIRE( P*P == P);
-
+
// PV=0 (lemma5)
REQUIRE( (P*boxCalculus.flat(f)).norm() == 0.0);
}
-
+
SECTION("Div / Curl")
{
PolygonalCalculus< RealPoint,RealVector >::Face f = 0;
@@ -182,12 +182,12 @@ TEST_CASE( "Testing PolygonalCalculus" )
//Not a great test BTW
REQUIRE(curl.norm() == 2.0);
}
-
+
SECTION("Local Laplace-Beltrami")
{
PolygonalCalculus< RealPoint,RealVector >::Face f = 0;
auto nf = box.incidentVertices(f).size();
-
+
auto L = boxCalculus.laplaceBeltrami(f);
PolygonalCalculus< RealPoint,RealVector >::Vector phi(nf),expected(nf);
phi << 1.0, 1.0, 1.0, 1.0;
@@ -212,7 +212,7 @@ TEST_CASE( "Testing PolygonalCalculus" )
REQUIRE( det == Approx(1.0));
REQUIRE( lphi[2] == Approx(-3.683));
}
-
+
SECTION("Covariant Operators")
{
PolygonalCalculus< RealPoint,RealVector >::Face f = 0;
@@ -238,33 +238,33 @@ TEST_CASE( "Testing PolygonalCalculus" )
double a=0.0;
for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index v=0; v < box.nbVertices(); ++v )
a += M.coeffRef(v,v);
-
+
double fa=0.0;
for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index f=0; f < box.nbFaces(); ++f )
fa += box.faceArea(f);
REQUIRE( a == fa );
}
-
+
SECTION("Checking cache")
{
auto cacheU = boxCalculus.getOperatorCacheMatrix( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.sharp(f);} );
REQUIRE( cacheU.size() == 6 );
-
+
auto cacheC = boxCalculus.getOperatorCacheVector( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.centroid(f);} );
REQUIRE( cacheC.size() == 6 );
}
-
+
SECTION("Internal cache")
{
PolygonalCalculus< RealPoint,RealVector > boxCalculusCached(box,true);
trace.info()<< boxCalculusCached <::SparseMatrix L(box.nbVertices(),box.nbVertices());
for(auto i=0; i < 1000 ; ++i)
L += i*boxCalculus.globalLaplaceBeltrami();
auto tps = trace.endBlock();
-
+
trace.beginBlock("With cache");
PolygonalCalculus< RealPoint,RealVector >::SparseMatrix LC(box.nbVertices(),box.nbVertices());
for(auto i=0; i < 1000 ; ++i)
@@ -272,9 +272,9 @@ TEST_CASE( "Testing PolygonalCalculus" )
auto tpsC = trace.endBlock();
REQUIRE(tpsC < tps);
REQUIRE(L.norm() == Approx(LC.norm()));
-
+
}
-
+
}
TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
@@ -287,7 +287,7 @@ TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
// Build a more complex surface.
auto params = SH3::defaultParameters();
-
+
params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
auto implicit_shape = SH3::makeImplicitShape3D ( params );
auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
@@ -295,16 +295,16 @@ TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
auto surface = SH3::makeDigitalSurface( binary_image, K, params );
auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
-
+
std::vector > faces;
std::vector positions = primalSurface->positions();
for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index face= 0 ; face < primalSurface->nbFaces(); ++face)
faces.push_back(primalSurface->incidentVertices( face ));
-
+
Mesh surfmesh = Mesh( positions.begin(), positions.end(),
faces.begin(), faces.end() );
auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
-
+
// Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
PolyDEC calculus( surfmesh );
// Laplace opeartor
@@ -324,7 +324,7 @@ TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
return cos(scale*(surfmesh.position(v)[0]))
* (scale*surfmesh.position(v)[1]);
};
-
+
for(auto &e: boundaryEdges)
{
auto adjVertices = surfmesh.edgeVertices(e);