DGtal  1.4.beta
testPolygonalCalculus.cpp
Go to the documentation of this file.
1
33 #include <iostream>
34 #include "DGtal/base/Common.h"
35 #include "ConfigTest.h"
36 #include "DGtalCatch.h"
37 #include "DGtal/helpers/StdDefs.h"
38
39 #include "DGtal/dec/PolygonalCalculus.h"
40 #include "DGtal/shapes/SurfaceMesh.h"
41 #include "DGtal/shapes/MeshHelpers.h"
42 #include "DGtal/helpers/Shortcuts.h"
43 #include "DGtal/math/linalg/DirichletConditions.h"
45
46 using namespace std;
47 using namespace DGtal;
48 using namespace Z3i;
49
51 // Functions for testing class PolygonalCalculus.
53
55 {
56  return RealPoint(v(0),v(1),v(2));
57 }
58
59 TEST_CASE( "Testing PolygonalCalculus" )
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 ff=0; ff<6; ++ff)
119  REQUIRE( boxCalculus.faceArea(ff) == box.faceArea(ff) );
120
121  box.computeFaceNormalsFromPositions();
122  for(auto ff=0; ff<6; ++ff)
123  {
124  auto cn = boxCalculus.faceNormalAsDGtalVector(ff);
125  auto n = box.faceNormal(ff);
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
156  auto gphi = G*phi;
158  auto cogphi = coG*phi;
159
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)
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;
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() == (Eigen::Index)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  {
238  double a=0.0;
239  for( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index v=0; v < box.nbVertices(); ++v )
240  a += M.coeffRef(v,v);
241
242  double fa=0.0;
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 }
279
280 TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
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( PolygonalCalculus< RealPoint,RealVector >::MySurfaceMesh::Index 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  {
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;
361  {
362  min_phi = std::min( min_phi, phi( v ) );
363  max_phi = std::max( max_phi, phi( v ) );
364  min_u = std::min( min_u , u ( v ) );
365  max_u = std::max( max_u , u ( v ) );
366  if ( b( v ) == 0.0 )
367  {
368  min_i_u = std::min( min_i_u, u ( v ) );
369  max_i_u = std::max( max_i_u, u ( v ) );
370  }
371  }
372  REQUIRE( min_phi <= min_u );
373  REQUIRE( max_phi >= max_u );
374  REQUIRE( min_phi < min_i_u );
375  REQUIRE( max_phi > max_i_u );
376  } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
377  }
378 };
379
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
Size nbFaces() const
VertexStorage::size_type Index
Definition: Mesh.h:120
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
Implements differential operators on polygonal surfaces from .
DenseMatrix D(const Face f) const
DenseMatrix laplaceBeltrami(const Face f, const double lambda=1.0) const
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
DenseMatrix X(const Face f) const
Vector faceNormal(const Face f) const
Vector centroid(const Face f) const
DenseMatrix curl(const Face f) const
Real3dVector faceNormalAsDGtalVector(const Face f) const
Real3dPoint centroidAsDGtalPoint(const Face f) const
DenseMatrix covariantProjection(const Face f, const Vector &uf)
Vector vectorArea(const Face f) const
size_t faceDegree(Face f) const
std::vector< DenseMatrix > getOperatorCacheMatrix(const std::function< DenseMatrix(Face)> &perFaceOperator) const
DenseMatrix P(const Face f) const
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
DenseMatrix sharp(const Face f) const
double faceArea(const Face f) const
DenseMatrix connectionLaplacian(const Face &f, double lambda=1.0) const
SparseMatrix globalLumpedMassMatrix() const
MySurfaceMesh::Face Face
Face type.
std::vector< Vector > getOperatorCacheVector(const std::function< Vector(Face)> &perFaceVectorOperator) const
LinAlg::DenseVector Vector
Type of Vector.
DenseMatrix A(const Face f) const
DenseMatrix covariantGradient(const Face f, const Vector &uf)
DenseMatrix flat(const Face f) const
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
DigitalPlane::Point Vector
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:153
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:647
std::size_t Index
The type used for numbering vertices and faces.
Definition: SurfaceMesh.h:105
Size nbFaces() const
Definition: SurfaceMesh.h:296
Edges computeManifoldBoundaryEdges() const
VertexPair edgeVertices(Edge e) const
Definition: SurfaceMesh.h:338
Size nbVertices() const
Definition: SurfaceMesh.h:288
Shortcuts< KSpace > SH3
int max(int a, int b)
KSpace K
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
TEST_CASE("Testing PolygonalCalculus")
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
SECTION("Testing constant forward iterators")
REQUIRE(domain.isInside(aPoint))
PointVector< 3, double > RealPoint