DGtal  1.3.beta
testPolygonalCalculus.cpp
Go to the documentation of this file.
1 
32 #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 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)
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  SECTION("Check lumped mass matrix")
199  {
201  double a=0.0;
202  for(auto v=0; v < box.nbVertices(); ++v )
203  a += M.coeffRef(v,v);
204 
205  double fa=0.0;
206  for(auto f=0; f < box.nbFaces(); ++f )
207  fa += box.faceArea(f);
208  REQUIRE( a == fa );
209  }
210 
211  SECTION("Checking cache")
212  {
213  auto cacheU = boxCalculus.getOperatorCacheMatrix( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.sharp(f);} );
214  REQUIRE( cacheU.size() == 6 );
215 
216  auto cacheC = boxCalculus.getOperatorCacheVector( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.centroid(f);} );
217  REQUIRE( cacheC.size() == 6 );
218  }
219 
220  SECTION("Internal cache")
221  {
222  PolygonalCalculus< RealPoint,RealVector > boxCalculusCached(box,true);
223  trace.info()<< boxCalculusCached <<std::endl;
224 
225  trace.beginBlock("Without cache");
226  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix L(box.nbVertices(),box.nbVertices());
227  for(auto i=0; i < 1000 ; ++i)
228  L += i*boxCalculus.globalLaplaceBeltrami();
229  auto tps = trace.endBlock();
230 
231  trace.beginBlock("With cache");
232  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix LC(box.nbVertices(),box.nbVertices());
233  for(auto i=0; i < 1000 ; ++i)
234  LC += i*boxCalculusCached.globalLaplaceBeltrami();
235  auto tpsC = trace.endBlock();
236  REQUIRE(tpsC < tps);
237  REQUIRE(L.norm() == Approx(LC.norm()));
238 
239  }
240 
241 }
242 
243 TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
244 {
245  typedef Shortcuts< KSpace > SH3;
247  typedef Mesh::Index Index;
250 
251  // Build a more complex surface.
252  auto params = SH3::defaultParameters();
253 
254  params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
255  auto implicit_shape = SH3::makeImplicitShape3D ( params );
256  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
257  auto K = SH3::getKSpace( params );
258  auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
259  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
260  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
261 
262  std::vector<std::vector< Index > > faces;
263  std::vector<RealPoint> positions = primalSurface->positions();
264  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
265  faces.push_back(primalSurface->incidentVertices( face ));
266 
267  Mesh surfmesh = Mesh( positions.begin(), positions.end(),
268  faces.begin(), faces.end() );
269  auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
270 
271  SECTION("Check surface")
272  {
273  REQUIRE( surfmesh.nbVertices() == 1364 );
274  REQUIRE( surfmesh.nbFaces() == 1279 );
275  REQUIRE( boundaryEdges.size() == 168 );
276  }
277 
278  // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
279  PolyDEC calculus( surfmesh );
280  // Laplace opeartor
282  // value on boundary
283  PolyDEC::Vector g = calculus.form0();
284  // characteristic set of boundary
285  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
286 
287  SECTION("Solve Poisson problem with boundary Dirichlet conditions")
288  {
289  for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
290  {
291  std::cout << "scale=" << scale << std::endl;
292  auto phi = [&]( Index v)
293  {
294  return cos(scale*(surfmesh.position(v)[0]))
295  * (scale*surfmesh.position(v)[1]);
296  };
297 
298  for(auto &e: boundaryEdges)
299  {
300  auto adjVertices = surfmesh.edgeVertices(e);
301  auto v1 = adjVertices.first;
302  auto v2 = adjVertices.second;
303  g(v1) = phi(v1);
304  g(v2) = phi(v2);
305  b(v1) = 1;
306  b(v2) = 1;
307  }
308  // Solve Δu=0 with g as boundary conditions
309  PolyDEC::Solver solver;
310  PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
311  solver.compute( L_dirichlet );
312  REQUIRE( solver.info() == Eigen::Success );
313  PolyDEC::Vector g_dirichlet = DC::dirichletVector( L, g, b, g );
314  PolyDEC::Vector x_dirichlet = solver.solve( g_dirichlet );
315  REQUIRE( solver.info() == Eigen::Success );
316  PolyDEC::Vector u = DC::dirichletSolution( x_dirichlet, b, g );
317  double min_phi = 0.0;
318  double max_phi = 0.0;
319  double min_u = 0.0;
320  double max_u = 0.0;
321  double min_i_u = 0.0;
322  double max_i_u = 0.0;
323  double exp_u = 0.0;
324  double exp_u2 = 0.0;
325  for ( auto v = 0; v < surfmesh.nbVertices(); ++v )
326  {
327  min_phi = std::min( min_phi, phi( v ) );
328  max_phi = std::max( max_phi, phi( v ) );
329  min_u = std::min( min_u , u ( v ) );
330  max_u = std::max( max_u , u ( v ) );
331  if ( b( v ) == 0.0 )
332  {
333  min_i_u = std::min( min_i_u, u ( v ) );
334  max_i_u = std::max( max_i_u, u ( v ) );
335  }
336  }
337  REQUIRE( min_phi <= min_u );
338  REQUIRE( max_phi >= max_u );
339  REQUIRE( min_phi < min_i_u );
340  REQUIRE( max_phi > max_i_u );
341  } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
342  }
343 };
344 
DGtal::PolygonalCalculus::Face
MySurfaceMesh::Face Face
Face type.
Definition: PolygonalCalculus.h:87
vecToRealPoint
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
Definition: testPolygonalCalculus.cpp:54
DGtal::PolygonalCalculus::getOperatorCacheMatrix
std::vector< DenseMatrix > getOperatorCacheMatrix(const std::function< DenseMatrix(Face)> &perFaceOperator) const
Definition: PolygonalCalculus.h:629
DGtal::PolygonalCalculus::vectorArea
Vector vectorArea(const Face f) const
Definition: PolygonalCalculus.h:275
DGtal::SurfaceMesh::position
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:536
DGtal::Trace::endBlock
double endBlock()
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
DGtal::PolygonalCalculus::faceNormalAsDGtalVector
Real3dVector faceNormalAsDGtalVector(const Face f) const
Definition: PolygonalCalculus.h:317
max
int max(int a, int b)
Definition: testArithmeticalDSS.cpp:1108
DGtal::PolygonalCalculus::globalLaplaceBeltrami
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Definition: PolygonalCalculus.h:551
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:103
DGtal::PolygonalCalculus::centroidAsDGtalPoint
Real3dPoint centroidAsDGtalPoint(const Face f) const
Definition: PolygonalCalculus.h:395
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::trace
Trace trace
Definition: Common.h:154
K
KSpace K
Definition: testCubicalComplex.cpp:62
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::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::PolygonalCalculus::faceArea
double faceArea(const Face f) const
Definition: PolygonalCalculus.h:299
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Shortcuts::makePrimalSurfaceMesh
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
REQUIRE
REQUIRE(domain.isInside(aPoint))
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:70
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::PolygonalCalculus::coGradient
DenseMatrix coGradient(const Face f) const
Definition: PolygonalCalculus.h:326
DGtal::PolygonalCalculus::form0
Form form0() const
Definition: PolygonalCalculus.h:524
DGtal::SurfaceMesh::nbVertices
Size nbVertices() const
Definition: SurfaceMesh.h:280
SparseMatrix
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: testHeatLaplace.cpp:50
DGtal::PolygonalCalculus::P
DenseMatrix P(const Face f) const
Definition: PolygonalCalculus.h:420
DGtal::Trace::info
std::ostream & info()
DGtal::PolygonalCalculus::D
DenseMatrix D(const Face f) const
Definition: PolygonalCalculus.h:220
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::PolygonalCalculus::sharp
DenseMatrix sharp(const Face f) const
Definition: PolygonalCalculus.h:404
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::LaplaceBeltrami
DenseMatrix LaplaceBeltrami(const Face f, const double lambda=1.0) const
Definition: PolygonalCalculus.h:504
DGtal::PolygonalCalculus::getOperatorCacheVector
std::vector< Vector > getOperatorCacheVector(const std::function< Vector(Face)> &perFaceVectorOperator) const
Definition: PolygonalCalculus.h:653
DGtal::PolygonalCalculus::gradient
DenseMatrix gradient(const Face f) const
Definition: PolygonalCalculus.h:349
DGtal::Shortcuts::makeImplicitShape3D
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::PolygonalCalculus::faceDegree
size_t faceDegree(Face f) const
Definition: PolygonalCalculus.h:693
DGtal::PolygonalCalculus::centroid
Vector centroid(const Face f) const
Definition: PolygonalCalculus.h:387
DGtal::SurfaceMesh::edgeVertices
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
DGtal::PolygonalCalculus::isValid
bool isValid() const
Definition: PolygonalCalculus.h:741
DGtal::Mesh::nbFaces
Size nbFaces() const
DGtal::PolygonalCalculus::curl
DenseMatrix curl(const Face f) const
Definition: PolygonalCalculus.h:478
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
Vector
FreemanChain< int >::Vector Vector
Definition: testCombinDSS.cpp:60
DGtal::PolygonalCalculus::nbVertices
size_t nbVertices() const
Definition: PolygonalCalculus.h:699
TEST_CASE
TEST_CASE("Testing PolygonalCalculus")
Definition: testPolygonalCalculus.cpp:59
DGtal::PolygonalCalculus::X
DenseMatrix X(const Face f) const
Definition: PolygonalCalculus.h:195
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:97
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
DGtal::PolygonalCalculus::faceNormal
Vector faceNormal(const Face f) const
Definition: PolygonalCalculus.h:307
DGtal::PolygonalCalculus::flat
DenseMatrix flat(const Face f) const
Definition: PolygonalCalculus.h:363
DGtal::PolygonalCalculus::globalLumpedMassMatrix
SparseMatrix globalLumpedMassMatrix() const
Definition: PolygonalCalculus.h:578
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288
DGtal::PolygonalCalculus::A
DenseMatrix A(const Face f) const
Definition: PolygonalCalculus.h:254