DGtal  1.5.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  // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
309  PolyDEC calculus( surfmesh );
310  // Laplace opeartor
312  // value on boundary
314  // characteristic set of boundary
315  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
316 
317  SECTION("Solve Poisson problem with boundary Dirichlet conditions")
318  {
319  for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
320  {
321  std::cout << "scale=" << scale << std::endl;
322  auto phi = [&]( Index v)
323  {
324  return cos(scale*(surfmesh.position(v)[0]))
325  * (scale*surfmesh.position(v)[1]);
326  };
327 
328  for(auto &e: boundaryEdges)
329  {
330  auto adjVertices = surfmesh.edgeVertices(e);
331  auto v1 = adjVertices.first;
332  auto v2 = adjVertices.second;
333  g(v1) = phi(v1);
334  g(v2) = phi(v2);
335  b(v1) = 1;
336  b(v2) = 1;
337  }
338  // Solve Δu=0 with g as boundary conditions
339  PolyDEC::Solver solver;
340  PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
341  solver.compute( L_dirichlet );
342  REQUIRE( solver.info() == Eigen::Success );
343  PolyDEC::Vector g_dirichlet = DC::dirichletVector( L, g, b, g );
344  PolyDEC::Vector x_dirichlet = solver.solve( g_dirichlet );
345  REQUIRE( solver.info() == Eigen::Success );
346  PolyDEC::Vector u = DC::dirichletSolution( x_dirichlet, b, g );
347  double min_phi = 0.0;
348  double max_phi = 0.0;
349  double min_u = 0.0;
350  double max_u = 0.0;
351  double min_i_u = 0.0;
352  double max_i_u = 0.0;
354  {
355  min_phi = std::min( min_phi, phi( v ) );
356  max_phi = std::max( max_phi, phi( v ) );
357  min_u = std::min( min_u , u ( v ) );
358  max_u = std::max( max_u , u ( v ) );
359  if ( b( v ) == 0.0 )
360  {
361  min_i_u = std::min( min_i_u, u ( v ) );
362  max_i_u = std::max( max_i_u, u ( v ) );
363  }
364  }
365  REQUIRE( min_phi <= min_u );
366  REQUIRE( max_phi >= max_u );
367  REQUIRE( min_phi < min_i_u );
368  REQUIRE( max_phi > max_i_u );
369  } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
370  }
371 };
372 
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()
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi(const Face f)
SurfMesh surfmesh
PolyCalculus * calculus
CountedPtr< SH3::DigitalSurface > surface
CountedPtr< SH3::BinaryImage > binary_image
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
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