DGtal 1.3.0
Loading...
Searching...
No Matches
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
46using namespace std;
47using namespace DGtal;
48using namespace Z3i;
49
51// Functions for testing class PolygonalCalculus.
53
55{
56 return RealPoint(v(0),v(1),v(2));
57}
58
59TEST_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
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 {
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}
279
280TEST_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(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
318 PolyDEC::SparseMatrix L = calculus.globalLaplaceBeltrami();
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};
381
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
std::vector< DenseMatrix > getOperatorCacheMatrix(const std::function< DenseMatrix(Face)> &perFaceOperator) 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
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
std::vector< Vector > getOperatorCacheVector(const std::function< Vector(Face)> &perFaceVectorOperator) 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.
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
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
static Parameters defaultParameters()
Definition: Shortcuts.h:203
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
SMesh::Index Index
Space::RealPoint RealPoint
Definition: StdDefs.h:170
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
Size nbFaces() const
Definition: SurfaceMesh.h:288
Edges computeManifoldBoundaryEdges() const
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
Size nbVertices() const
Definition: SurfaceMesh.h:280
Shortcuts< KSpace > SH3
KSpace K
TEST_CASE("Testing PolygonalCalculus")
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
SECTION("Testing constant forward iterators")
REQUIRE(domain.isInside(aPoint))