DGtal  1.4.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 "DGtal/helpers/Shortcuts.h"
#include "DGtal/math/linalg/DirichletConditions.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")
 
 TEST_CASE ("Testing PolygonalCalculus and DirichletConditions")
 

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
Jacques-Olivier Lachaud (jacqu.nosp@m.es-o.nosp@m.livie.nosp@m.r.la.nosp@m.chaud.nosp@m.@uni.nosp@m.v-sav.nosp@m.oie..nosp@m.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.

Definition in file testPolygonalCalculus.cpp.

Function Documentation

◆ TEST_CASE() [1/2]

TEST_CASE ( "Testing PolygonalCalculus and DirichletConditions )

Definition at line 280 of file testPolygonalCalculus.cpp.

281 {
282  typedef Shortcuts< KSpace > SH3;
284  typedef Mesh::Index Index;
287 
288  // Build a more complex surface.
289  auto params = SH3::defaultParameters();
290 
291  params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
292  auto implicit_shape = SH3::makeImplicitShape3D ( params );
293  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
294  auto K = SH3::getKSpace( params );
295  auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
296  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
297  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
298 
299  std::vector<std::vector< Index > > faces;
300  std::vector<RealPoint> positions = primalSurface->positions();
301  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
302  faces.push_back(primalSurface->incidentVertices( face ));
303 
304  Mesh surfmesh = Mesh( positions.begin(), positions.end(),
305  faces.begin(), faces.end() );
306  auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
307 
308  SECTION("Check surface")
309  {
310  REQUIRE( surfmesh.nbVertices() == 1364 );
311  REQUIRE( surfmesh.nbFaces() == 1279 );
312  REQUIRE( boundaryEdges.size() == 168 );
313  }
314 
315  // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
316  PolyDEC calculus( surfmesh );
317  // Laplace opeartor
319  // value on boundary
320  PolyDEC::Vector g = calculus.form0();
321  // characteristic set of boundary
322  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
323 
324  SECTION("Solve Poisson problem with boundary Dirichlet conditions")
325  {
326  for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
327  {
328  std::cout << "scale=" << scale << std::endl;
329  auto phi = [&]( Index v)
330  {
331  return cos(scale*(surfmesh.position(v)[0]))
332  * (scale*surfmesh.position(v)[1]);
333  };
334 
335  for(auto &e: boundaryEdges)
336  {
337  auto adjVertices = surfmesh.edgeVertices(e);
338  auto v1 = adjVertices.first;
339  auto v2 = adjVertices.second;
340  g(v1) = phi(v1);
341  g(v2) = phi(v2);
342  b(v1) = 1;
343  b(v2) = 1;
344  }
345  // Solve Δu=0 with g as boundary conditions
346  PolyDEC::Solver solver;
347  PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
348  solver.compute( L_dirichlet );
349  REQUIRE( solver.info() == Eigen::Success );
350  PolyDEC::Vector g_dirichlet = DC::dirichletVector( L, g, b, g );
351  PolyDEC::Vector x_dirichlet = solver.solve( g_dirichlet );
352  REQUIRE( solver.info() == Eigen::Success );
353  PolyDEC::Vector u = DC::dirichletSolution( x_dirichlet, b, g );
354  double min_phi = 0.0;
355  double max_phi = 0.0;
356  double min_u = 0.0;
357  double max_u = 0.0;
358  double min_i_u = 0.0;
359  double max_i_u = 0.0;
360  double exp_u = 0.0;
361  double exp_u2 = 0.0;
362  for ( auto v = 0; v < surfmesh.nbVertices(); ++v )
363  {
364  min_phi = std::min( min_phi, phi( v ) );
365  max_phi = std::max( max_phi, phi( v ) );
366  min_u = std::min( min_u , u ( v ) );
367  max_u = std::max( max_u , u ( v ) );
368  if ( b( v ) == 0.0 )
369  {
370  min_i_u = std::min( min_i_u, u ( v ) );
371  max_i_u = std::max( max_i_u, u ( v ) );
372  }
373  }
374  REQUIRE( min_phi <= min_u );
375  REQUIRE( max_phi >= max_u );
376  REQUIRE( min_phi < min_i_u );
377  REQUIRE( max_phi > max_i_u );
378  } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
379  }
380 };
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Definition: Mesh.h:92
VertexStorage::size_type Index
Definition: Mesh.h:120
Implements differential operators on polygonal surfaces from .
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
SMesh::Index Index
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
Size nbFaces() const
Definition: SurfaceMesh.h:288
Edges computeManifoldBoundaryEdges() const
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
Size nbVertices() const
Definition: SurfaceMesh.h:280
Shortcuts< KSpace > SH3
int max(int a, int b)
FreemanChain< int >::Vector Vector
KSpace K
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
SECTION("Testing constant forward iterators")
REQUIRE(domain.isInside(aPoint))

References DGtal::SurfaceMesh< TRealPoint, TRealVector >::computeManifoldBoundaryEdges(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::edgeVertices(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::form0(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::globalLaplaceBeltrami(), K, max(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbFaces(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbVertices(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::position(), REQUIRE(), and SECTION().

◆ TEST_CASE() [2/2]

TEST_CASE ( "Testing PolygonalCalculus )

Definition at line 59 of file testPolygonalCalculus.cpp.

60 {
62  std::vector<RealPoint> positions = { RealPoint( 0, 0, 0 ) ,
63  RealPoint( 1, 0, 0 ) ,
64  RealPoint( 0, 1, 0 ) ,
65  RealPoint( 1, 1, 0 ) ,
66  RealPoint( 0, 0, 1 ) ,
67  RealPoint( 1, 0, 1 ) ,
68  RealPoint( 0, 1, 1 ) ,
69  RealPoint( 1, 1, 1 ) ,
70  RealPoint( 1, 0, 2 ) ,
71  RealPoint( 0, 0, 2 ) };
72  std::vector<Mesh::Vertices> faces = { { 1, 0, 2, 3 },
73  { 0, 1, 5, 4 } ,
74  { 1, 3, 7, 5 } ,
75  { 3, 2, 6, 7 } ,
76  { 2, 0, 4, 6 } ,
77  { 4, 5, 8, 9 } };
78 
79  Mesh box(positions.cbegin(), positions.cend(),
80  faces.cbegin(), faces.cend());
81 
83 
84  SECTION("Construction and basic operators")
85  {
86  REQUIRE( boxCalculus.isValid() );
87  REQUIRE( boxCalculus.nbVertices() == positions.size() );
88 
90  auto x = boxCalculus.X(f);
91  auto d = boxCalculus.D(f);
92  auto a = boxCalculus.A(f);
93 
94  //Checking X
96  REQUIRE( vecToRealPoint(vec ) == positions[1]);
97  vec = x.row(1);
98  REQUIRE( vecToRealPoint(vec ) == positions[0]);
99  vec = x.row(2);
100  REQUIRE( vecToRealPoint(vec ) == positions[2]);
101  vec = x.row(3);
102  REQUIRE( vecToRealPoint(vec ) == positions[3]);
103 
104  trace.info()<< boxCalculus <<std::endl;
105 
106  //Some D and A
107  REQUIRE( d(1,1) == -1 );
108  REQUIRE( d(0,1) == 1 );
109  REQUIRE( d(0,2) == 0 );
110 
111  REQUIRE( a(1,1) == 0.5 );
112  REQUIRE( a(0,1) == 0.5 );
113  REQUIRE( a(0,2) == 0 );
114 
115  auto vectorArea = boxCalculus.vectorArea(f);
116 
117  //Without correction, this should match
118  for(auto f=0; f<6; ++f)
119  REQUIRE( boxCalculus.faceArea(f) == box.faceArea(f) );
120 
121  box.computeFaceNormalsFromPositions();
122  for(auto f=0; f<6; ++f)
123  {
124  auto cn = boxCalculus.faceNormalAsDGtalVector(f);
125  auto n = box.faceNormal(f);
126  REQUIRE( cn == n );
127  }
128 
129  RealPoint c = boxCalculus.centroidAsDGtalPoint(f);
130  RealPoint expected(0.5,0.5,0.0);
131  REQUIRE(c == expected);
132  }
133 
134  SECTION("Derivatives")
135  {
137  auto d = boxCalculus.D(f);
138 
139  auto nf = boxCalculus.faceDegree(f);
141  phi << 1.0, 3.0, 2.0, 6.0;
142  expected << 2,-1,4,-5;
143  auto dphi = d*phi; // n_f x 1 matrix
144  REQUIRE(dphi == expected);
145 
146  }
147 
148  SECTION("Structural propertes")
149  {
151  auto nf = boxCalculus.faceDegree(f);
153  phi << 1.0, 3.0, 2.0, 6.0;
154 
155  auto G = boxCalculus.gradient(f);
156  auto gphi = G*phi;
157  auto coG = boxCalculus.coGradient(f);
158  auto cogphi = coG*phi;
159 
160  // grad . cograd == 0
161  REQUIRE( gphi.dot(cogphi) == 0.0);
162 
163  // Gf = Uf Df
164  REQUIRE( G == boxCalculus.sharp(f)*boxCalculus.D(f));
165 
166  // UV = I - nn^t (lemma4)
167  PolygonalCalculus< RealPoint,RealVector >::Vector n = boxCalculus.faceNormal(f);
168  REQUIRE( boxCalculus.sharp(f)*boxCalculus.flat(f) == PolygonalCalculus< RealPoint,RealVector >::DenseMatrix::Identity(3,3) - n*n.transpose() );
169 
170  // P^2 = P (lemma6)
171  auto P = boxCalculus.P(f);
172  REQUIRE( P*P == P);
173 
174  // PV=0 (lemma5)
175  REQUIRE( (P*boxCalculus.flat(f)).norm() == 0.0);
176  }
177 
178  SECTION("Div / Curl")
179  {
181  auto curl = boxCalculus.curl(f);
182  //Not a great test BTW
183  REQUIRE(curl.norm() == 2.0);
184  }
185 
186  SECTION("Local Laplace-Beltrami")
187  {
189  auto nf = box.incidentVertices(f).size();
190 
191  auto L = boxCalculus.laplaceBeltrami(f);
193  phi << 1.0, 1.0, 1.0, 1.0;
194  expected << 0,0,0,0;
195  auto lphi = L*phi;
196  REQUIRE( lphi == expected);
197  }
198 
199  SECTION("Local Connection-Laplace-Beltrami")
200  {
202  auto nf = box.incidentVertices(f).size();
203 
204  auto L = boxCalculus.connectionLaplacian(f);
206  phi << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0;
207  //since connection laplacian transports the phi vectors to the face,
208  //it's not expected to have 0 since these vectors aren't actually the same
209  auto lphi = L*phi;
210  //but we can still check that L is semi-definite
211  double det = L.determinant()+1.;
212  REQUIRE( det == Approx(1.0));
213  REQUIRE( lphi[2] == Approx(-3.683));
214  }
215 
216  SECTION("Covariant Operators")
217  {
219  auto nf = box.incidentVertices(f).size();
220 
222  phi << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0;
223  auto CG = boxCalculus.covariantGradient(f,phi);
224  auto CP = boxCalculus.covariantProjection(f,phi);
225 
226  //check sizes
227  REQUIRE( CG.rows() == 2);
228  REQUIRE( CG.cols() == 2);
229  REQUIRE( CP.rows() == faces[f].size());
230  REQUIRE( CP.cols() == 2);
231 
232  REQUIRE( CG(0,0) == Approx(0.707106));
233  REQUIRE( CP(0,0) == Approx(1.224744));
234  }
235  SECTION("Check lumped mass matrix")
236  {
237  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix M = boxCalculus.globalLumpedMassMatrix();
238  double a=0.0;
239  for(auto v=0; v < box.nbVertices(); ++v )
240  a += M.coeffRef(v,v);
241 
242  double fa=0.0;
243  for(auto f=0; f < box.nbFaces(); ++f )
244  fa += box.faceArea(f);
245  REQUIRE( a == fa );
246  }
247 
248  SECTION("Checking cache")
249  {
250  auto cacheU = boxCalculus.getOperatorCacheMatrix( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.sharp(f);} );
251  REQUIRE( cacheU.size() == 6 );
252 
253  auto cacheC = boxCalculus.getOperatorCacheVector( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.centroid(f);} );
254  REQUIRE( cacheC.size() == 6 );
255  }
256 
257  SECTION("Internal cache")
258  {
259  PolygonalCalculus< RealPoint,RealVector > boxCalculusCached(box,true);
260  trace.info()<< boxCalculusCached <<std::endl;
261 
262  trace.beginBlock("Without cache");
263  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix L(box.nbVertices(),box.nbVertices());
264  for(auto i=0; i < 1000 ; ++i)
265  L += i*boxCalculus.globalLaplaceBeltrami();
266  auto tps = trace.endBlock();
267 
268  trace.beginBlock("With cache");
269  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix LC(box.nbVertices(),box.nbVertices());
270  for(auto i=0; i < 1000 ; ++i)
271  LC += i*boxCalculusCached.globalLaplaceBeltrami();
272  auto tpsC = trace.endBlock();
273  REQUIRE(tpsC < tps);
274  REQUIRE(L.norm() == Approx(LC.norm()));
275 
276  }
277 
278 }
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
double faceArea(const Face f) const
MySurfaceMesh::Face Face
Face type.
LinAlg::DenseVector Vector
Type of Vector.
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
Trace trace
Definition: Common.h:154
Z2i::RealPoint RealPoint
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)

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 >::connectionLaplacian(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::covariantGradient(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::covariantProjection(), 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 54 of file testPolygonalCalculus.cpp.

55 {
56  return RealPoint(v(0),v(1),v(2));
57 }

Referenced by TEST_CASE().