DGtal 1.4.0
Loading...
Searching...
No Matches
testLinearStructure.cpp
Go to the documentation of this file.
1
27#include "DGtal/math/linalg/EigenSupport.h"
28#include "DGtal/dec/DiscreteExteriorCalculus.h"
29#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
30
31//#include "DGtal/io/viewers/Viewer3D.h"
32#include "DGtal/io/boards/Board2D.h"
33#include "DGtal/base/Common.h"
34#include "DGtal/helpers/StdDefs.h"
35
36#include "boost/version.hpp"
37
38using namespace DGtal;
39using namespace Z2i;
40using std::endl;
41
42#if BOOST_VERSION >= 105000
43#define TEST_HARDCODED_ORDER
44#endif
45
46void test_linear_structure()
47{
48 trace.beginBlock("creating dec problem with neumann border condition");
49
51 typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
52
53 SCells cells;
54
55 // fill cells container
56 {
57 const Calculus::KSpace kspace;
58
59 for (int kk=20; kk>0; kk--) cells.push_back( kspace.sCell(Point(0,kk), kk%2 != 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS) );
60 for (int kk=0; kk<10; kk++) cells.push_back( kspace.sCell(Point(kk,0)) );
61 for (int kk=0; kk<10; kk++) cells.push_back( kspace.sCell(Point(10,kk)) );
62 cells.push_back( kspace.sCell(Point(10,10)) );
63 cells.push_back( kspace.sCell(Point(9,10), Calculus::KSpace::NEG) );
64 for (int kk=10; kk<20; kk++) cells.push_back( kspace.sCell(Point(8,kk)) );
65 for (int kk=8; kk<12; kk++) cells.push_back( kspace.sCell(Point(kk,20)) );
66 for (int kk=20; kk>0; kk--) cells.push_back( kspace.sCell(Point(12,kk), kk%2 != 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS) );
67 cells.push_back( kspace.sCell(Point(12,0)) );
68 }
69
70 // fill calculus
71 const Domain domain(Point(-1,-1), Point(10,10));
72
73 // create DEC structure
74 Calculus calculus;
75 calculus.initKSpace<Domain>(domain);
76 for (SCells::const_iterator ci=cells.begin(), ce=cells.end(); ci!=ce; ci++) calculus.insertSCell( *ci );
77 calculus.updateIndexes();
79
80 trace.info() << calculus << endl;
81
83 Calculus::Cell dirac_cell = calculus.myKSpace.uCell(Point(10,4));
84 Calculus::PrimalForm0 dirac = Calculus::PrimalForm0::dirac(calculus, dirac_cell);
86 trace.info() << "dirac_cell_index=" << calculus.getCellIndex(dirac_cell) << endl;
87
88 {
89 Board2D board;
90 board << domain;
91 board << calculus;
92 board << dirac;
93 board.saveSVG("linear_structure_neumann_dirac.svg");
94 }
95
96 trace.endBlock();
97
98 {
99 trace.beginBlock("solving problem with neumann border condition using sparse qr solver");
100
102 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
104 trace.info() << "laplace = " << laplace << endl;
105
106 const Calculus::PrimalIdentity0 reorder = calculus.reorder<0, PRIMAL>(cells.begin(), cells.end());
107 const Calculus::PrimalIdentity0 laplace_ordered = reorder * laplace * reorder.transpose();
108 trace.info() << Eigen::MatrixXd(laplace_ordered.myContainer) << endl;
109
111 typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
112 typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;
113
114 Solver solver;
115 solver.compute(laplace);
116 Calculus::PrimalForm0 solved_solution = solver.solve(dirac);
118 Calculus::PrimalForm0 solved_solution_ordered = reorder * solved_solution;
119
120 Calculus::PrimalForm0 analytic_solution(calculus);
121 {
122 const Calculus::Index dirac_position = 17;
123 const Calculus::Index length = analytic_solution.length();
124 for (Calculus::Index kk=0; kk<length; kk++)
125 {
126 Calculus::Scalar alpha = 1. * (kk)/dirac_position * (kk+1.)/(dirac_position+1.);
127 if (kk>dirac_position)
128 {
129 alpha = 1. * (length-kk)/dirac_position * (length-kk-1.)/(dirac_position+1);
130 alpha -= 1. * (length-dirac_position)/dirac_position * (length-dirac_position-1.)/(dirac_position+1);
131 alpha += 1;
132 }
133 analytic_solution.myContainer(kk) = alpha;
134 }
135 }
136
137 const double shift = solved_solution_ordered.myContainer[0]-analytic_solution.myContainer[0];
138 for (Calculus::Index index=0; index<solved_solution_ordered.length(); index++) solved_solution_ordered.myContainer[index] -= shift;
139 solved_solution_ordered.myContainer /= solved_solution_ordered.myContainer.maxCoeff();
140
141 trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
142
143 {
144 std::ofstream handle("linear_structure_neumann.dat");
145 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
146 {
147 handle << solved_solution_ordered.myContainer(kk) << " " << analytic_solution.myContainer(kk) << endl;
148 }
149 }
150
151 const double error = (solved_solution_ordered-analytic_solution).myContainer.array().abs().maxCoeff();
152 trace.info() << "error=" << error << endl;
153 FATAL_ERROR( error < 1e-5 );
154
155 {
156 Board2D board;
157 board << domain;
158 board << calculus;
159 board << solved_solution;
160 board.saveSVG("linear_structure_neumann_solution.svg");
161 }
162
163 {
164 Calculus::PrimalForm1 solved_solution_gradient = calculus.derivative<0, PRIMAL>() * solved_solution;
165 Board2D board;
166 board << domain;
167 board << calculus;
168 board << solved_solution_gradient;
169 board << CustomStyle("VectorField", new VectorFieldStyle2D(1));
170 board << calculus.sharp(solved_solution_gradient);
171 board.saveSVG("linear_structure_neumann_solution_gradient.svg");
172 }
173
174 trace.endBlock();
175 }
176
177 trace.beginBlock("creating dec problem with dirichlet border condition");
178
180 calculus.insertSCell( calculus.myKSpace.sCell(Point(13,0)) );
181 calculus.insertSCell( calculus.myKSpace.sCell(Point(1,20), Calculus::KSpace::NEG) );
182 calculus.updateIndexes();
184
185 dirac = Calculus::PrimalForm0::dirac(calculus, dirac_cell);
186 trace.info() << calculus << endl;
187 trace.info() << "dirac_cell_index=" << calculus.getCellIndex(dirac_cell) << endl;
188
189 {
190 Board2D board;
191 board << domain;
192 board << calculus;
193 board << dirac;
194 board.saveSVG("linear_structure_dirichlet_dirac.svg");
195 }
196
197 trace.endBlock();
198
199 {
200 trace.beginBlock("solving problem with dirichlet border condition using sparse qr solver");
201
203 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
205 trace.info() << "laplace = " << laplace << endl;
206
207 const Calculus::PrimalIdentity0 reorder = calculus.reorder<0, PRIMAL>(cells.begin(), cells.end());
208 const Calculus::PrimalIdentity0 laplace_ordered = reorder * laplace * reorder.transpose();
209 trace.info() << Eigen::MatrixXd(laplace_ordered.myContainer) << endl;
210
212 typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
213 typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;
214
215 Solver solver;
216 solver.compute(laplace);
217 Calculus::PrimalForm0 solved_solution = solver.solve(dirac);
219 solved_solution.myContainer.array() /= solved_solution.myContainer.maxCoeff();
220 Calculus::PrimalForm0 solved_solution_ordered = reorder * solved_solution;
221
222 Calculus::PrimalForm0 analytic_solution(calculus);
223 {
224 const Calculus::Index dirac_position = 17;
225 const Calculus::Index length = analytic_solution.length();
226 for (Calculus::Index kk=0; kk<length; kk++)
227 {
228 Calculus::Scalar alpha = (kk+1.)/(dirac_position+1.);
229 if (kk>dirac_position)
230 {
231 alpha = 1. - (kk-dirac_position)/(1.*length-dirac_position);
232 }
233 analytic_solution.myContainer(kk) = alpha;
234 }
235 }
236
237 const double shift = solved_solution_ordered.myContainer[0]-analytic_solution.myContainer[0];
238 for (Calculus::Index index=0; index<solved_solution_ordered.length(); index++) solved_solution_ordered.myContainer[index] -= shift;
239 solved_solution_ordered.myContainer /= solved_solution_ordered.myContainer.maxCoeff();
240
241 trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
242
243 {
244 std::ofstream handle("linear_structure_dirichlet.dat");
245 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
246 {
247 handle << solved_solution_ordered.myContainer(kk) << " " << analytic_solution.myContainer(kk) << endl;
248 }
249 }
250
251 const double error = (solved_solution_ordered-analytic_solution).myContainer.array().abs().maxCoeff();
252 trace.info() << "error=" << error << endl;
253 FATAL_ERROR( error < 1e-5 );
254
255 {
256 Board2D board;
257 board << domain;
258 board << calculus;
259 board << solved_solution;
260 board.saveSVG("linear_structure_dirichlet_solution.svg");
261 }
262
263 {
264 Calculus::PrimalForm1 solved_solution_gradient = calculus.derivative<0, PRIMAL>() * solved_solution;
265
266 Board2D board;
267 board << domain;
268 board << calculus;
269 board << solved_solution_gradient;
270 board << calculus.sharp(solved_solution_gradient);
271 board.saveSVG("linear_structure_dirichlet_solution_gradient.svg");
272 }
273
274 trace.endBlock();
275 }
276
277}
278
279template <typename Operator>
280void display_operator_info(const std::string& name, const Operator& op)
281{
282 trace.info() << name << endl << op << endl << Eigen::MatrixXd(op.myContainer) << endl;
283}
284
285void test_linear_ring()
286{
287 trace.beginBlock("linear ring");
288
289 const Domain domain(Point(-5,-5), Point(5,5));
290
291 typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
292 Calculus calculus;
293 calculus.initKSpace<Domain>(domain);
294
295 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(-8,kk), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
296 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,10), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
297 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(10,kk)) );
298 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,-8)) );
299 calculus.updateIndexes();
300
301 {
302 trace.info() << calculus << endl;
303 Board2D board;
304 board << domain;
305 board << calculus;
306 board.saveSVG("ring_structure.svg");
307 }
308
309 const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
310 display_operator_info("d0", d0);
311
312 const Calculus::PrimalHodge0 h0 = calculus.hodge<0, PRIMAL>();
313 display_operator_info("h0", h0);
314
315 const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
316 display_operator_info("d0p", d0p);
317
318 const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
319 display_operator_info("h1", h1);
320
321 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
322 display_operator_info("laplace", laplace);
323
324 const auto laplace_size = calculus.kFormLength(0, PRIMAL);
325 const Eigen::MatrixXd laplace_dense(laplace.myContainer);
326
327 for (int ii=0; ii<laplace_size; ii++)
328 FATAL_ERROR( laplace_dense(ii,ii) == 2 );
329
330 FATAL_ERROR( laplace_dense.array().rowwise().sum().abs().sum() == 0 );
331 FATAL_ERROR( laplace_dense.transpose() == laplace_dense );
332
333 trace.endBlock();
334}
335
336void test_manual_operators_3d()
337{
338 trace.beginBlock("testing 3d operators");
339
340 const Z3i::Domain domain(Z3i::Point(0,0,0), Z3i::Point(1,1,1));
341
342 typedef DiscreteExteriorCalculus<3, 3, EigenLinearAlgebraBackend> Calculus;
343
344 Calculus calculus;
345 calculus.initKSpace<Z3i::Domain>(domain);
346
347 { // filling primal calculus
348 // 0-cells
349 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,0,0)) );
350 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,0,0)) );
351
352 // 1-cells
353 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,1,0), Calculus::KSpace::POS) );
354 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,0,1), Calculus::KSpace::POS) );
355 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,0,0), Calculus::KSpace::POS) );
356 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,1,0), Calculus::KSpace::NEG) );
357 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,0,1), Calculus::KSpace::NEG) );
358
359 // 2-cells
360 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,1,1), Calculus::KSpace::NEG) );
361 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,0,1), Calculus::KSpace::POS) ); //FIXME strange cell orientation
362 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,1,1), Calculus::KSpace::POS) );
363 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,1,0), Calculus::KSpace::NEG) );
364
365 // 3-cells
366 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,1,1)) );
367
368 calculus.updateIndexes();
369
370 trace.info() << calculus << endl;
371 }
372
373 trace.beginBlock("base operators");
374
375 {
376 const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
377 const Calculus::DualDerivative2 d2p = calculus.derivative<2, DUAL>();
378 display_operator_info("d0", d0);
379 display_operator_info("d2p", d2p);
380
381#if defined(TEST_HARDCODED_ORDER)
382 Eigen::MatrixXd d0_th(5, 2);
383 d0_th <<
384 -1, 0,
385 -1, 0,
386 -1, 1,
387 0, 1,
388 0, 1;
389
390 FATAL_ERROR( Eigen::MatrixXd(d0.myContainer) == d0_th );
391#endif
392 FATAL_ERROR( Eigen::MatrixXd(d2p.transpose().myContainer) == Eigen::MatrixXd(d0.myContainer) );
393 }
394
395 {
396 const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();
397 const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();
398 display_operator_info("d1", d1);
399 display_operator_info("d1p", d1p);
400
401#if defined(TEST_HARDCODED_ORDER)
402 Eigen::MatrixXd d1_th(4, 5);
403 d1_th <<
404 0, -1, 1, 0, -1,
405 1, 0, -1, 1, 0,
406 0, 0, 0, -1, 1,
407 -1, 1, 0, 0, 0;
408
409 FATAL_ERROR( Eigen::MatrixXd(d1.myContainer) == d1_th );
410#endif
411 FATAL_ERROR( Eigen::MatrixXd(d1p.transpose().myContainer) == Eigen::MatrixXd(d1.myContainer) );
412 }
413
414 {
415 const Calculus::PrimalDerivative2 d2 = calculus.derivative<2, PRIMAL>();
416 const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
417 display_operator_info("d2", d2);
418 display_operator_info("d0p", d0p);
419
420#if defined(TEST_HARDCODED_ORDER)
421 Eigen::MatrixXd d2_th(1, 4);
422 d2_th << -1, -1, -1, -1;
423
424 FATAL_ERROR( Eigen::MatrixXd(d2.myContainer) == d2_th );
425#endif
426 FATAL_ERROR( Eigen::MatrixXd(d0p.transpose().myContainer) == Eigen::MatrixXd(d2.myContainer) );
427 }
428
429 trace.endBlock();
430
431 /*
432 QApplication app(0, NULL);
433
434 typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
435 Viewer* viewer = new Viewer(calculus.myKSpace);
436 viewer->show();
437 viewer->setWindowTitle("structure");
438 (*viewer) << CustomColors3D(DGtal::Color(255,0,0), DGtal::Color(0,0,0));
439 (*viewer) << domain;
440 Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus);
441 (*viewer) << Viewer::updateDisplay;
442
443 app.exec();
444 */
445
446 trace.endBlock();
447}
448
449void test_manual_operators_2d()
450{
451 trace.beginBlock("testing 2d operators");
452
453 const Domain domain(Point(0,0), Point(5,4));
454
455 typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;
456
457 Calculus primal_calculus;
458 primal_calculus.initKSpace<Domain>(domain);
459
460 { // filling primal calculus
461 // 0-cells
462 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,2)) );
463 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,2)) );
464 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,4)) );
465 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,4)) );
466 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,6)) );
467 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,6)) );
468
469 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(1,2)) ); // insert cell
470
471 // 1-cells
472 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,2), Calculus::KSpace::NEG) ); // insert negative cell
473 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,2), Calculus::KSpace::POS) ); // then reinserting negative cell in structure
474 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,3), Calculus::KSpace::NEG) );
475 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,3), Calculus::KSpace::POS) );
476 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,4), Calculus::KSpace::NEG) );
477 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,5), Calculus::KSpace::POS) );
478 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,5), Calculus::KSpace::NEG) );
479 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,6), Calculus::KSpace::POS) );
480
481 primal_calculus.eraseCell( primal_calculus.myKSpace.uCell(Point(1,2)) ); // then remove it
482
483 // 2-cells
484 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,3)) );
485 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,5)) );
486
487 primal_calculus.updateIndexes();
488
489 trace.beginBlock("primal calculus");
490 trace.info() << primal_calculus << endl;
491 for (Calculus::ConstIterator iter_property=primal_calculus.begin(), iter_property_end=primal_calculus.end(); iter_property!=iter_property_end; iter_property++)
492 {
493 const Calculus::Cell cell = iter_property->first;
494 const Calculus::Property property = iter_property->second;
495 const Dimension dim = primal_calculus.myKSpace.uDim(cell);
496 const Calculus::SCell signed_cell = primal_calculus.myKSpace.signs(cell, property.flipped ? Calculus::KSpace::NEG : Calculus::KSpace::POS);
497
498 ASSERT( signed_cell == primal_calculus.getSCell(dim, PRIMAL, property.index) );
499
500 trace.info() << cell
501 << " " << dim
502 << " " << signed_cell
503 << " " << property.primal_size
504 << " " << property.dual_size
505 << " " << property.index
506 << " " << (property.flipped ? "negative" : "positive")
507 << endl;
508 }
509 trace.endBlock();
510 }
511
512
513 Calculus dual_calculus;
514 dual_calculus.initKSpace<Domain>(domain);
515
516 { // filling dual calculus
517 // 2-cells
518 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,3)) );
519 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,3)) );
520 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,5)) );
521 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,5)) );
522 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,7)) );
523 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,7)) );
524
525 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(6,3)) ); // insert cell
526
527 // 1-cells
528 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,3), Calculus::KSpace::NEG) ); // insert negative cell
529 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,3), Calculus::KSpace::POS) ); // then reinserting negative cell in structure
530 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,4), Calculus::KSpace::POS) );
531 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,4), Calculus::KSpace::NEG) );
532 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,5), Calculus::KSpace::NEG) );
533 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,6), Calculus::KSpace::NEG) );
534 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,6), Calculus::KSpace::POS) );
535 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,7), Calculus::KSpace::POS) );
536
537 dual_calculus.eraseCell( dual_calculus.myKSpace.uCell(Point(6,3)) ); // then remove it
538
539 // 0-cells
540 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,4)) );
541 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,6)) );
542
543 dual_calculus.updateIndexes();
544
545 trace.beginBlock("dual calculus");
546 trace.info() << dual_calculus << endl;
547 for (Calculus::ConstIterator iter_property=dual_calculus.begin(), iter_property_end=dual_calculus.end(); iter_property!=iter_property_end; iter_property++)
548 {
549 const Calculus::Cell cell = iter_property->first;
550 const Calculus::Property property = iter_property->second;
551 const Dimension dim = dual_calculus.myKSpace.uDim(cell);
552 const Calculus::SCell signed_cell = dual_calculus.myKSpace.signs(cell, property.flipped ? Calculus::KSpace::NEG : Calculus::KSpace::POS);
553
554 ASSERT( signed_cell == dual_calculus.getSCell(dim, PRIMAL, property.index) );
555
556 trace.info() << cell
557 << " " << dim
558 << " " << signed_cell
559 << " " << property.primal_size
560 << " " << property.dual_size
561 << " " << property.index
562 << " " << (property.flipped ? "negative" : "positive")
563 << endl;
564 }
565 trace.endBlock();
566 }
567
568 {
569 Board2D board;
570 board << domain;
571 board << primal_calculus;
572 board << dual_calculus;
573 board.saveSVG("operators_structure.svg");
574 }
575
576 trace.beginBlock("base operators");
577
578 const Calculus::PrimalDerivative0 primal_d0 = primal_calculus.derivative<0, PRIMAL>();
579 const Calculus::DualDerivative0 dual_d0p = dual_calculus.derivative<0, DUAL>();
580 {
581 display_operator_info("primal d0", primal_d0);
582 display_operator_info("dual d0p", dual_d0p);
583
584#if defined(TEST_HARDCODED_ORDER)
585 Eigen::MatrixXd d0_th(7, 6);
586 d0_th <<
587 -1, 1, 0, 0, 0, 0,
588 1, 0, 0, -1, 0, 0,
589 0, 0, 1, 0, 0, -1,
590 0, 0, 0, 0, -1, 1,
591 0, -1, 1, 0, 0, 0,
592 0, 0, -1, 1, 0, 0,
593 0, 0, 0, -1, 1, 0;
594 FATAL_ERROR( Eigen::MatrixXd(primal_d0.myContainer) == d0_th );
595
596 Eigen::MatrixXd d0p_th(7, 6);
597 d0p_th <<
598 1, -1, 0, 0, 0, 0,
599 0, 0, -1, 1, 0, 0,
600 0, 1, 0, -1, 0, 0,
601 0, 0, 1, 0, -1, 0,
602 0, 0, 0, 0, 1, -1,
603 -1, 0, 1, 0, 0, 0,
604 0, 0, 0, -1, 0, 1;
605 FATAL_ERROR( Eigen::MatrixXd(dual_d0p.myContainer) == d0p_th );
606#endif
607 }
608
609 const Calculus::PrimalDerivative1 primal_d1 = primal_calculus.derivative<1, PRIMAL>();
610 const Calculus::DualDerivative1 dual_d1p = dual_calculus.derivative<1, DUAL>();
611 {
612 display_operator_info("primal d1", primal_d1);
613 display_operator_info("dual d1p", dual_d1p);
614
615#if defined(TEST_HARDCODED_ORDER)
616 Eigen::MatrixXd d1_th(2, 7);
617 d1_th <<
618 1, 1, 0, 0, 1, 1, 0,
619 0, 0, -1, -1, 0, -1, -1;
620 FATAL_ERROR( Eigen::MatrixXd(primal_d1.myContainer) == d1_th );
621
622 Eigen::MatrixXd d1p_th(2, 7);
623 d1p_th <<
624 -1, -1, -1, 0, 0, -1, 0,
625 0, 1, 0, 1, 1, 0, 1;
626 FATAL_ERROR( Eigen::MatrixXd(dual_d1p.myContainer) == d1p_th );
627#endif
628 }
629
630 {
631 display_operator_info("primal d1*d0", primal_d1*primal_d0);
632 display_operator_info("dual d1p*d0p", dual_d1p*dual_d0p);
633
634 FATAL_ERROR( Eigen::MatrixXd((primal_d1*primal_d0).myContainer) == Eigen::MatrixXd::Zero(2,6) );
635 FATAL_ERROR( Eigen::MatrixXd((dual_d1p*dual_d0p).myContainer) == Eigen::MatrixXd::Zero(2,6) );
636 }
637
638 const Calculus::PrimalHodge0 primal_h0 = primal_calculus.hodge<0, PRIMAL>();
639 const Calculus::DualHodge0 dual_h0p = dual_calculus.hodge<0, DUAL>();
640 const Calculus::DualHodge2 primal_h2p = primal_calculus.hodge<2, DUAL>();
641 const Calculus::PrimalHodge2 dual_h2 = dual_calculus.hodge<2, PRIMAL>();
642 {
643 display_operator_info("primal h0", primal_h0);
644 display_operator_info("dual h0p", dual_h0p);
645 display_operator_info("primal h2p", primal_h2p);
646 display_operator_info("dual h2", dual_h2);
647
648 FATAL_ERROR( Eigen::MatrixXd(primal_h0.myContainer) == Eigen::MatrixXd::Identity(6,6) );
649 FATAL_ERROR( Eigen::MatrixXd(dual_h0p.myContainer) == Eigen::MatrixXd::Identity(6,6) );
650 FATAL_ERROR( Eigen::MatrixXd(primal_h2p.myContainer) == Eigen::MatrixXd::Identity(6,6) );
651 FATAL_ERROR( Eigen::MatrixXd(dual_h2.myContainer) == Eigen::MatrixXd::Identity(6,6) );
652 }
653
654 const Calculus::PrimalHodge2 primal_h2 = primal_calculus.hodge<2, PRIMAL>();
655 const Calculus::DualHodge2 dual_h2p = dual_calculus.hodge<2, DUAL>();
656 const Calculus::DualHodge0 primal_h0p = primal_calculus.hodge<0, DUAL>();
657 const Calculus::PrimalHodge0 dual_h0 = dual_calculus.hodge<0, PRIMAL>();
658 {
659 display_operator_info("primal h2", primal_h2);
660 display_operator_info("dual h2p", dual_h2p);
661 display_operator_info("primal h0p", primal_h0p);
662 display_operator_info("dual h0", dual_h0);
663
664 FATAL_ERROR( Eigen::MatrixXd(primal_h2.myContainer) == Eigen::MatrixXd::Identity(2,2) );
665 FATAL_ERROR( Eigen::MatrixXd(dual_h2p.myContainer) == Eigen::MatrixXd::Identity(2,2) );
666 FATAL_ERROR( Eigen::MatrixXd(primal_h0p.myContainer) == Eigen::MatrixXd::Identity(2,2) );
667 FATAL_ERROR( Eigen::MatrixXd(dual_h0.myContainer) == Eigen::MatrixXd::Identity(2,2) );
668 }
669
670 const Calculus::DualDerivative0 primal_d0p = primal_calculus.derivative<0, DUAL>();
671 const Calculus::PrimalDerivative0 dual_d0 = dual_calculus.derivative<0, PRIMAL>();
672 {
673 display_operator_info("primal d0p", primal_d0p);
674 display_operator_info("dual d0", dual_d0);
675
676#if defined(TEST_HARDCODED_ORDER)
677 Eigen::MatrixXd d0p_th_transpose(2, 7);
678 d0p_th_transpose <<
679 1, 1, 0, 0, 1, 1, 0,
680 0, 0, -1, -1, 0, -1, -1;
681 FATAL_ERROR( Eigen::MatrixXd(primal_d0p.myContainer) == d0p_th_transpose.transpose() );
682
683 Eigen::MatrixXd minus_d0_th_transpose(2, 7);
684 minus_d0_th_transpose <<
685 -1, -1, -1, 0, 0, -1, 0,
686 0, 1, 0, 1, 1, 0, 1;
687 FATAL_ERROR( Eigen::MatrixXd(dual_d0.myContainer) == -minus_d0_th_transpose.transpose() );
688#endif
689 }
690
691 const Calculus::DualDerivative1 primal_d1p = primal_calculus.derivative<1, DUAL>();
692 const Calculus::PrimalDerivative1 dual_d1 = dual_calculus.derivative<1, PRIMAL>();
693 {
694 display_operator_info("primal d1p", primal_d1p);
695 display_operator_info("dual d1", dual_d1);
696
697#if defined(TEST_HARDCODED_ORDER)
698 Eigen::MatrixXd minus_d1p_th_transpose(7, 6);
699 minus_d1p_th_transpose <<
700 -1, 1, 0, 0, 0, 0,
701 1, 0, 0, -1, 0, 0,
702 0, 0, 1, 0, 0, -1,
703 0, 0, 0, 0, -1, 1,
704 0, -1, 1, 0, 0, 0,
705 0, 0, -1, 1, 0, 0,
706 0, 0, 0, -1, 1, 0;
707 FATAL_ERROR( Eigen::MatrixXd(primal_d1p.myContainer) == -minus_d1p_th_transpose.transpose() );
708
709 Eigen::MatrixXd d1_th_transpose(7, 6);
710 d1_th_transpose <<
711 1, -1, 0, 0, 0, 0,
712 0, 0, -1, 1, 0, 0,
713 0, 1, 0, -1, 0, 0,
714 0, 0, 1, 0, -1, 0,
715 0, 0, 0, 0, 1, -1,
716 -1, 0, 1, 0, 0, 0,
717 0, 0, 0, -1, 0, 1;
718 FATAL_ERROR( Eigen::MatrixXd(dual_d1.myContainer) == d1_th_transpose.transpose() );
719#endif
720 }
721
722 const Calculus::PrimalHodge1 primal_h1 = primal_calculus.hodge<1, PRIMAL>();
723 const Calculus::DualHodge1 dual_h1p = dual_calculus.hodge<1, DUAL>();
724 const Calculus::DualHodge1 primal_h1p = primal_calculus.hodge<1, DUAL>();
725 const Calculus::PrimalHodge1 dual_h1 = dual_calculus.hodge<1, PRIMAL>();
726 {
727 display_operator_info("primal h1", primal_h1);
728 display_operator_info("dual h1p", dual_h1p);
729 display_operator_info("primal h1p", primal_h1p);
730 display_operator_info("dual h1", dual_h1);
731
732 FATAL_ERROR( Eigen::MatrixXd(primal_h1.myContainer) == Eigen::MatrixXd::Identity(7,7) );
733 FATAL_ERROR( Eigen::MatrixXd(dual_h1p.myContainer) == -Eigen::MatrixXd::Identity(7,7) );
734 FATAL_ERROR( Eigen::MatrixXd((primal_h1p*primal_h1).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
735 FATAL_ERROR( Eigen::MatrixXd((dual_h1*dual_h1p).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
736 FATAL_ERROR( Eigen::MatrixXd((primal_h1*primal_h1p).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
737 FATAL_ERROR( Eigen::MatrixXd((dual_h1p*dual_h1).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
738 }
739
740 trace.endBlock();
741
742 trace.beginBlock("laplace operators");
743 display_operator_info("primal laplace", primal_calculus.laplace<PRIMAL>());
744 display_operator_info("dual laplacep", dual_calculus.laplace<DUAL>());
745 display_operator_info("primal laplacep", primal_calculus.laplace<DUAL>());
746 display_operator_info("dual laplace", dual_calculus.laplace<PRIMAL>());
747 trace.endBlock();
748
749 trace.beginBlock("sharp operators");
750
751 { // primal sharp
752 display_operator_info("primal #x", primal_calculus.sharpDirectional<PRIMAL>(0));
753 display_operator_info("dual #xp", dual_calculus.sharpDirectional<DUAL>(0));
754 display_operator_info("primal #y", primal_calculus.sharpDirectional<PRIMAL>(1));
755 display_operator_info("dual #yp", dual_calculus.sharpDirectional<DUAL>(1));
756
757 {
758 Calculus::PrimalForm1::Container dx_container(7);
759 dx_container << -1, 0, 0, -1, 0, 1, 0;
760 const Calculus::PrimalForm1 primal_dx(primal_calculus, dx_container);
761 const Calculus::PrimalVectorField primal_dx_field = primal_calculus.sharp(primal_dx);
762
763 Calculus::PrimalForm1::Container dxp_container(7);
764 dxp_container << 1, -1, 0, 0, 1, 0, 0;
765 const Calculus::DualForm1 dual_dx(dual_calculus, dxp_container);
766 const Calculus::DualVectorField dual_dx_field = dual_calculus.sharp(dual_dx);
767
768 {
769 Board2D board;
770 board << domain;
771 board << primal_calculus;
772 board << primal_dx << primal_dx_field;
773 board << dual_calculus;
774 board << dual_dx << dual_dx_field;
775 board.saveSVG("operators_sharp_dx_primal.svg");
776 }
777
778#if defined(TEST_HARDCODED_ORDER)
779 FATAL_ERROR( primal_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(6) );
780 FATAL_ERROR( primal_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(6) );
781 FATAL_ERROR( dual_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(6) );
782 FATAL_ERROR( dual_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(6) );
783#endif
784 }
785
786 {
787 Calculus::PrimalForm1::Container dy_container(7);
788 dy_container << 0, 1, 1, 0, -1, 0, -1;
789 const Calculus::PrimalForm1 primal_dy(primal_calculus, dy_container);
790 const Calculus::PrimalVectorField primal_dy_field = primal_calculus.sharp(primal_dy);
791
792 Calculus::PrimalForm1::Container dyp_container(7);
793 dyp_container << 0, 0, 1, 1, 0, -1, -1;
794 const Calculus::DualForm1 dual_dy(dual_calculus, dyp_container);
795 const Calculus::DualVectorField dual_dy_field = dual_calculus.sharp(dual_dy);
796
797 {
798 Board2D board;
799 board << domain;
800 board << primal_calculus;
801 board << primal_dy << primal_dy_field;
802 board << dual_calculus;
803 board << dual_dy << dual_dy_field;
804 board.saveSVG("operators_sharp_dy_primal.svg");
805 }
806
807#if defined(TEST_HARDCODED_ORDER)
808 FATAL_ERROR( primal_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(6) );
809 FATAL_ERROR( primal_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(6) );
810 FATAL_ERROR( dual_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(6) );
811 FATAL_ERROR( dual_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(6) );
812#endif
813 }
814 }
815
816 { // dual sharp
817 display_operator_info("primal #xp", primal_calculus.sharpDirectional<DUAL>(0));
818 display_operator_info("dual #x", dual_calculus.sharpDirectional<PRIMAL>(0));
819 display_operator_info("primal #yp", primal_calculus.sharpDirectional<DUAL>(1));
820 display_operator_info("dual #y", dual_calculus.sharpDirectional<PRIMAL>(1));
821
822 {
823 Calculus::DualForm1::Container dx_container(7);
824 dx_container << 0, -1, -1, 0, 1, 0, 1;
825 const Calculus::DualForm1 primal_dx(primal_calculus, dx_container);
826 const Calculus::DualVectorField primal_dx_field = primal_calculus.sharp(primal_dx);
827
828 Calculus::DualForm1::Container dxp_container(7);
829 dxp_container << 0, 0, 1, 1, 0, -1, -1;
830 const Calculus::PrimalForm1 dual_dx(dual_calculus, dxp_container);
831 const Calculus::PrimalVectorField dual_dx_field = dual_calculus.sharp(dual_dx);
832
833 {
834 Board2D board;
835 board << domain;
836 board << primal_calculus;
837 board << primal_dx << primal_dx_field;
838 board << dual_calculus;
839 board << dual_dx << dual_dx_field;
840 board.saveSVG("operators_sharp_dx_dual.svg");
841 }
842
843#if defined(TEST_HARDCODED_ORDER)
844 FATAL_ERROR( primal_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(2) );
845 FATAL_ERROR( primal_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(2) );
846 FATAL_ERROR( dual_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(2) );
847 FATAL_ERROR( dual_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(2) );
848#endif
849 }
850
851 {
852 Calculus::DualForm1::Container dy_container(7);
853 dy_container << -1, 0, 0, -1, 0 , 1, 0;
854 const Calculus::DualForm1 primal_dy(primal_calculus, dy_container);
855 const Calculus::DualVectorField primal_dy_field = primal_calculus.sharp(primal_dy);
856
857 Calculus::DualForm1::Container dyp_container(7);
858 dyp_container << -1, 1, 0, 0, -1, 0, 0;
859 const Calculus::PrimalForm1 dual_dy(dual_calculus, dyp_container);
860 const Calculus::PrimalVectorField dual_dy_field = dual_calculus.sharp(dual_dy);
861
862 {
863 Board2D board;
864 board << domain;
865 board << primal_calculus;
866 board << primal_dy << primal_dy_field;
867 board << dual_calculus;
868 board << dual_dy << dual_dy_field;
869 board.saveSVG("operators_sharp_dy_dual.svg");
870 }
871
872#if defined(TEST_HARDCODED_ORDER)
873 FATAL_ERROR( primal_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(2) );
874 FATAL_ERROR( primal_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(2) );
875 FATAL_ERROR( dual_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(2) );
876 FATAL_ERROR( dual_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(2) );
877#endif
878 }
879 }
880
881 trace.endBlock();
882
883 trace.beginBlock("flat operators");
884
885 { // primal flat
886 display_operator_info("primal bx", primal_calculus.flatDirectional<PRIMAL>(0));
887 display_operator_info("dual bxp", dual_calculus.flatDirectional<DUAL>(0));
888 display_operator_info("primal by", primal_calculus.flatDirectional<PRIMAL>(1));
889 display_operator_info("dual byp", dual_calculus.flatDirectional<DUAL>(1));
890
891 Calculus::PrimalVectorField::Coordinates dx_coords(6,2);
892 dx_coords.col(0) = Eigen::VectorXd::Ones(6);
893 dx_coords.col(1) = Eigen::VectorXd::Zero(6);
894
895 Calculus::PrimalVectorField::Coordinates dy_coords(6,2);
896 dy_coords.col(0) = Eigen::VectorXd::Zero(6);
897 dy_coords.col(1) = Eigen::VectorXd::Ones(6);
898
899 const Calculus::PrimalVectorField primal_dx_field(primal_calculus, dx_coords);
900 const Calculus::PrimalForm1 primal_dx = primal_calculus.flat(primal_dx_field);
901 const Calculus::DualVectorField dual_dx_field(dual_calculus, dx_coords);
902 const Calculus::DualForm1 dual_dx = dual_calculus.flat(dual_dx_field);
903
904 {
905 Board2D board;
906 board << domain;
907 board << primal_calculus;
908 board << primal_dx << primal_dx_field;
909 board << dual_calculus;
910 board << dual_dx << dual_dx_field;
911 board.saveSVG("operators_flat_dx_primal.svg");
912 }
913
914 const Calculus::PrimalVectorField primal_dy_field(primal_calculus, dy_coords);
915 const Calculus::PrimalForm1 primal_dy = primal_calculus.flat(primal_dy_field);
916 const Calculus::DualVectorField dual_dy_field(dual_calculus, dy_coords);
917 const Calculus::DualForm1 dual_dy = dual_calculus.flat(dual_dy_field);
918
919 {
920 Board2D board;
921 board << domain;
922 board << primal_calculus;
923 board << primal_dy << primal_dy_field;
924 board << dual_calculus;
925 board << dual_dy << dual_dy_field;
926 board.saveSVG("operators_flat_dy_primal.svg");
927 }
928
929#if defined(TEST_HARDCODED_ORDER)
930 Calculus::PrimalForm1::Container dx_container(7);
931 dx_container << -1, 0, 0, -1, 0, 1, 0;
932 Calculus::PrimalForm1::Container dxp_container(7);
933 dxp_container << 1, -1, 0, 0, 1, 0, 0;
934 FATAL_ERROR( primal_dx.myContainer == dx_container );
935 FATAL_ERROR( dual_dx.myContainer == dxp_container );
936
937 Calculus::PrimalForm1::Container dy_container(7);
938 dy_container << 0, 1, 1, 0, -1, 0, -1;
939 Calculus::PrimalForm1::Container dyp_container(7);
940 dyp_container << 0, 0, 1, 1, 0, -1, -1;
941 FATAL_ERROR( primal_dy.myContainer == dy_container );
942 FATAL_ERROR( dual_dy.myContainer == dyp_container );
943#endif
944 }
945
946 { // dual flat
947 display_operator_info("primal bxp", primal_calculus.flatDirectional<DUAL>(0));
948 display_operator_info("dual bx", dual_calculus.flatDirectional<PRIMAL>(0));
949 display_operator_info("primal byp", primal_calculus.flatDirectional<DUAL>(1));
950 display_operator_info("dual by", dual_calculus.flatDirectional<PRIMAL>(1));
951
952 Calculus::PrimalVectorField::Coordinates dx_coords(2,2);
953 dx_coords.col(0) = Eigen::VectorXd::Ones(2);
954 dx_coords.col(1) = Eigen::VectorXd::Zero(2);
955
956 Calculus::PrimalVectorField::Coordinates dy_coords(2,2);
957 dy_coords.col(0) = Eigen::VectorXd::Zero(2);
958 dy_coords.col(1) = Eigen::VectorXd::Ones(2);
959
960 const Calculus::DualVectorField primal_dx_field(primal_calculus, dx_coords);
961 const Calculus::DualForm1 primal_dx = primal_calculus.flat(primal_dx_field);
962 const Calculus::PrimalVectorField dual_dx_field(dual_calculus, dx_coords);
963 const Calculus::PrimalForm1 dual_dx = dual_calculus.flat(dual_dx_field);
964
965 {
966 Board2D board;
967 board << domain;
968 board << primal_calculus;
969 board << primal_dx << primal_dx_field;
970 board << dual_calculus;
971 board << dual_dx << dual_dx_field;
972 board.saveSVG("operators_flat_dx_dual.svg");
973 }
974
975 const Calculus::DualVectorField primal_dy_field(primal_calculus, dy_coords);
976 const Calculus::DualForm1 primal_dy = primal_calculus.flat(primal_dy_field);
977 const Calculus::PrimalVectorField dual_dy_field(dual_calculus, dy_coords);
978 const Calculus::PrimalForm1 dual_dy = dual_calculus.flat(dual_dy_field);
979
980 {
981 Board2D board;
982 board << domain;
983 board << primal_calculus;
984 board << primal_dy << primal_dy_field;
985 board << dual_calculus;
986 board << dual_dy << dual_dy_field;
987 board.saveSVG("operators_flat_dy_dual.svg");
988 }
989
990#if defined(TEST_HARDCODED_ORDER)
991 Calculus::PrimalForm1::Container dx_container(7);
992 dx_container << 0, -1, -1, 0, 1, 0, 1;
993 Calculus::PrimalForm1::Container dxp_container(7);
994 dxp_container << 0, 0, 1, 1, 0, -1, -1;
995 FATAL_ERROR( primal_dx.myContainer == dx_container );
996 FATAL_ERROR( dual_dx.myContainer == dxp_container );
997
998 Calculus::PrimalForm1::Container dy_container(7);
999 dy_container << -1, 0, 0, -1, 0, 1, 0;
1000 Calculus::PrimalForm1::Container dyp_container(7);
1001 dyp_container << -1, 1, 0, 0, -1, 0, 0;
1002 FATAL_ERROR( primal_dy.myContainer == dy_container );
1003 FATAL_ERROR( dual_dy.myContainer == dyp_container );
1004#endif
1005 }
1006
1007 trace.endBlock();
1008
1009
1010 trace.endBlock();
1011}
1012
1013
1014int
1015main()
1016{
1017#if !defined(TEST_HARDCODED_ORDER)
1018 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1019#endif
1020 test_manual_operators_3d();
1021 test_manual_operators_2d();
1022 test_linear_ring();
1023 test_linear_structure();
1024#if !defined(TEST_HARDCODED_ORDER)
1025 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1026#endif
1027 return 0;
1028}
1029
DGtal::Board2D
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition Board2D.h:71
DGtal::DiscreteExteriorCalculusSolver
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
Definition DiscreteExteriorCalculusSolver.h:70
DGtal::DiscreteExteriorCalculusSolver::solve
SolutionKForm solve(const InputKForm &input_kform) const
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition DiscreteExteriorCalculus.h:98
DGtal::DiscreteExteriorCalculus::initKSpace
void initKSpace(ConstAlias< TDomain > domain)
DGtal::KForm::myContainer
Container myContainer
Definition KForm.h:131
DGtal::PolygonalCalculus::sharp
DenseMatrix sharp(const Face f) const
Definition PolygonalCalculus.h:503
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::warning
std::ostream & warning()
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
LibBoard::Board::saveSVG
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition Board.cpp:1011
calculus
PolyCalculus * calculus
Definition dgtalCalculus-geodesic.cpp:70
DGtal::Z2i::SCells
KSpace::SCells SCells
Definition StdDefs.h:82
DGtal::Z2i::Point
Space::Point Point
Definition StdDefs.h:95
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition ClosedIntegerHalfPlane.h:49
DGtal::Dimension
DGtal::uint32_t Dimension
Definition Common.h:136
DGtal::trace
Trace trace
Definition Common.h:153
DGtal::PRIMAL
@ PRIMAL
Definition Duality.h:61
DGtal::DUAL
@ DUAL
Definition Duality.h:62
DGtal::EigenLinearAlgebraBackend::SolverSparseQR
Eigen::SparseQR< SparseMatrix, Eigen::COLAMDOrdering< SparseMatrix::Index > > SolverSparseQR
Definition EigenSupport.h:110
test_linear_structure
void test_linear_structure()
Definition testLinearStructure.cpp:46
test_linear_ring
void test_linear_ring()
Definition testLinearStructure.cpp:285
display_operator_info
void display_operator_info(const std::string &name, const Operator &op)
Definition testLinearStructure.cpp:280
test_manual_operators_2d
void test_manual_operators_2d()
Definition testLinearStructure.cpp:449
test_manual_operators_3d
void test_manual_operators_3d()
Definition testLinearStructure.cpp:336
domain
Domain domain
Definition testProjection.cpp:88