DGtal  1.3.beta
Functions
testPolygonalCalculus.cpp File Reference
#include <iostream>
#include "DGtal/base/Common.h"
#include "ConfigTest.h"
#include "DGtalCatch.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/dec/PolygonalCalculus.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/shapes/MeshHelpers.h"
Include dependency graph for testPolygonalCalculus.cpp:

Go to the source code of this file.

Functions

RealPoint vecToRealPoint (PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
 
 TEST_CASE ("Testing PolygonalCalculus")
 

Detailed Description

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 http://www.gnu.org/licenses/.

Author
David Coeurjolly (david.nosp@m..coe.nosp@m.urjol.nosp@m.ly@l.nosp@m.iris..nosp@m.cnrs.nosp@m..fr ) Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
Date
2021/09/02

Functions for testing class PolygonalCalculus.

This file is part of the DGtal library.

Definition in file testPolygonalCalculus.cpp.

Function Documentation

◆ TEST_CASE()

TEST_CASE ( "Testing PolygonalCalculus )

Definition at line 55 of file testPolygonalCalculus.cpp.

56 {
58  std::vector<RealPoint> positions = { RealPoint( 0, 0, 0 ) ,
59  RealPoint( 1, 0, 0 ) ,
60  RealPoint( 0, 1, 0 ) ,
61  RealPoint( 1, 1, 0 ) ,
62  RealPoint( 0, 0, 1 ) ,
63  RealPoint( 1, 0, 1 ) ,
64  RealPoint( 0, 1, 1 ) ,
65  RealPoint( 1, 1, 1 ) ,
66  RealPoint( 1, 0, 2 ) ,
67  RealPoint( 0, 0, 2 ) };
68  std::vector<Mesh::Vertices> faces = { { 1, 0, 2, 3 },
69  { 0, 1, 5, 4 } ,
70  { 1, 3, 7, 5 } ,
71  { 3, 2, 6, 7 } ,
72  { 2, 0, 4, 6 } ,
73  { 4, 5, 8, 9 } };
74 
75  Mesh box(positions.cbegin(), positions.cend(),
76  faces.cbegin(), faces.cend());
77 
79 
80  SECTION("Construction and basic operators")
81  {
82  REQUIRE( boxCalculus.isValid() );
83  REQUIRE( boxCalculus.nbVertices() == positions.size() );
84 
86  auto x = boxCalculus.X(f);
87  auto d = boxCalculus.D(f);
88  auto a = boxCalculus.A(f);
89 
90  //Checking X
92  REQUIRE( vecToRealPoint(vec ) == positions[1]);
93  vec = x.row(1);
94  REQUIRE( vecToRealPoint(vec ) == positions[0]);
95  vec = x.row(2);
96  REQUIRE( vecToRealPoint(vec ) == positions[2]);
97  vec = x.row(3);
98  REQUIRE( vecToRealPoint(vec ) == positions[3]);
99 
100  trace.info()<< boxCalculus <<std::endl;
101 
102  //Some D and A
103  REQUIRE( d(1,1) == -1 );
104  REQUIRE( d(0,1) == 1 );
105  REQUIRE( d(0,2) == 0 );
106 
107  REQUIRE( a(1,1) == 0.5 );
108  REQUIRE( a(0,1) == 0.5 );
109  REQUIRE( a(0,2) == 0 );
110 
111  auto vectorArea = boxCalculus.vectorArea(f);
112 
113  //Without correction, this should match
114  for(auto f=0; f<6; ++f)
115  REQUIRE( boxCalculus.faceArea(f) == box.faceArea(f) );
116 
117  box.computeFaceNormalsFromPositions();
118  for(auto f=0; f<6; ++f)
119  {
120  auto cn = boxCalculus.faceNormalAsDGtalVector(f);
121  auto n = box.faceNormal(f);
122  REQUIRE( cn == n );
123  }
124 
125  RealPoint c = boxCalculus.centroidAsDGtalPoint(f);
126  RealPoint expected(0.5,0.5,0.0);
127  REQUIRE(c == expected);
128  }
129 
130  SECTION("Derivatives")
131  {
133  auto d = boxCalculus.D(f);
134 
135  auto nf = boxCalculus.faceDegree(f);
137  phi << 1.0, 3.0, 2.0, 6.0;
138  expected << 2,-1,4,-5;
139  auto dphi = d*phi; // n_f x 1 matrix
140  REQUIRE(dphi == expected);
141 
142  }
143 
144  SECTION("Structural propertes")
145  {
147  auto nf = boxCalculus.faceDegree(f);
149  phi << 1.0, 3.0, 2.0, 6.0;
150 
151  auto G = boxCalculus.gradient(f);
152  auto gphi = G*phi;
153  auto coG = boxCalculus.coGradient(f);
154  auto cogphi = coG*phi;
155 
156  // grad . cograd == 0
157  REQUIRE( gphi.dot(cogphi) == 0.0);
158 
159  // Gf = Uf Df
160  REQUIRE( G == boxCalculus.sharp(f)*boxCalculus.D(f));
161 
162  // UV = I - nn^t (lemma4)
163  PolygonalCalculus< RealPoint,RealVector >::Vector n = boxCalculus.faceNormal(f);
164  REQUIRE( boxCalculus.sharp(f)*boxCalculus.flat(f) == PolygonalCalculus< RealPoint,RealVector >::DenseMatrix::Identity(3,3) - n*n.transpose() );
165 
166  // P^2 = P (lemma6)
167  auto P = boxCalculus.P(f);
168  REQUIRE( P*P == P);
169 
170  // PV=0 (lemma5)
171  REQUIRE( (P*boxCalculus.flat(f)).norm() == 0.0);
172  }
173 
174  SECTION("Div / Curl")
175  {
177  auto curl = boxCalculus.curl(f);
178  //Not a great test BTW
179  REQUIRE(curl.norm() == 2.0);
180  }
181 
182  SECTION("Local Laplace-Beltrami")
183  {
185  auto nf = box.incidentVertices(f).size();
186 
187  auto L = boxCalculus.LaplaceBeltrami(f);
189  phi << 1.0, 1.0, 1.0, 1.0;
190  expected << 0,0,0,0;
191  auto lphi = L*phi;
192  REQUIRE( lphi == expected);
193  }
194  SECTION("Check lumped mass matrix")
195  {
196  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix M = boxCalculus.globalLumpedMassMatrix();
197  double a=0.0;
198  for(auto v=0; v < box.nbVertices(); ++v )
199  a += M.coeffRef(v,v);
200 
201  double fa=0.0;
202  for(auto f=0; f < box.nbFaces(); ++f )
203  fa += box.faceArea(f);
204  REQUIRE( a == fa );
205  }
206 
207  SECTION("Checking cache")
208  {
209  auto cacheU = boxCalculus.getOperatorCacheMatrix( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.sharp(f);} );
210  REQUIRE( cacheU.size() == 6 );
211 
212  auto cacheC = boxCalculus.getOperatorCacheVector( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.centroid(f);} );
213  REQUIRE( cacheC.size() == 6 );
214  }
215 
216  SECTION("Internal cache")
217  {
218  PolygonalCalculus< RealPoint,RealVector > boxCalculusCached(box,true);
219  trace.info()<< boxCalculusCached <<std::endl;
220 
221  trace.beginBlock("Without cache");
222  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix L(box.nbVertices(),box.nbVertices());
223  for(auto i=0; i < 1000 ; ++i)
224  L += i*boxCalculus.globalLaplaceBeltrami();
225  auto tps = trace.endBlock();
226 
227  trace.beginBlock("With cache");
228  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix LC(box.nbVertices(),box.nbVertices());
229  for(auto i=0; i < 1000 ; ++i)
230  LC += i*boxCalculusCached.globalLaplaceBeltrami();
231  auto tpsC = trace.endBlock();
232  REQUIRE(tpsC < tps);
233  REQUIRE(L.norm() == Approx(LC.norm()));
234 
235  }
236 
237 }

References DGtal::PolygonalCalculus< TRealPoint, TRealVector >::A(), DGtal::Trace::beginBlock(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::centroid(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::centroidAsDGtalPoint(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::coGradient(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::curl(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::D(), DGtal::Trace::endBlock(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceArea(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceDegree(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceNormal(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceNormalAsDGtalVector(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::flat(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::getOperatorCacheMatrix(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::getOperatorCacheVector(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::globalLaplaceBeltrami(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::globalLumpedMassMatrix(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::gradient(), DGtal::Trace::info(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::isValid(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::LaplaceBeltrami(), DGtal::Mesh< TPoint >::nbFaces(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::nbVertices(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::P(), REQUIRE(), SECTION(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::sharp(), DGtal::trace, DGtal::PolygonalCalculus< TRealPoint, TRealVector >::vectorArea(), vecToRealPoint(), and DGtal::PolygonalCalculus< TRealPoint, TRealVector >::X().

◆ vecToRealPoint()

RealPoint vecToRealPoint ( PolygonalCalculus< RealPoint, RealPoint >::Vector v)

Definition at line 50 of file testPolygonalCalculus.cpp.

51 {
52  return RealPoint(v(0),v(1),v(2));
53 }

Referenced by TEST_CASE().

DGtal::PolygonalCalculus::Face
MySurfaceMesh::Face Face
Face type.
Definition: PolygonalCalculus.h:77
vecToRealPoint
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
Definition: testPolygonalCalculus.cpp:50
DGtal::Trace::endBlock
double endBlock()
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:91
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::PolygonalCalculus::faceArea
double faceArea(const Face f) const
Definition: PolygonalCalculus.h:279
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
REQUIRE
REQUIRE(domain.isInside(aPoint))
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:60
DGtal::Trace::info
std::ostream & info()
DGtal::Mesh
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Definition: Mesh.h:91
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:87
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
SECTION
SECTION("Testing constant forward iterators")
Definition: testSimpleRandomAccessRangeFromPoint.cpp:66
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46