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 
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)
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() == (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  {
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;
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 coGradient(const Face f) const
DenseMatrix connectionLaplacian(const Face &f, double lambda=1.0) const
SparseMatrix globalLumpedMassMatrix() const
DenseMatrix gradient(const Face f) 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:637
std::size_t Index
The type used for numbering vertices and faces.
Definition: SurfaceMesh.h:105
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)
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