DGtal  1.5.beta
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 
38 using namespace DGtal;
39 using namespace Z2i;
40 using std::endl;
41 
42 #if BOOST_VERSION >= 105000
43 #define TEST_HARDCODED_ORDER
44 #endif
45 
46 void 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 
279 template <typename Operator>
280 void 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 
285 void 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 
336 void 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 
449 void 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 
1014 int
1015 main()
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::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::info
std::ostream & info()
DGtal::Trace::warning
std::ostream & warning()
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
SCell
Z3i::SCell SCell
Definition: digitalPolyhedronBuilder3D.cpp:79
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
ConstIterator
MyDigitalSurface::ConstIterator ConstIterator
Definition: greedy-plane-segmentation-ex2.cpp:93
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.
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
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
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
main
int main()
Definition: testLinearStructure.cpp:1015
test_manual_operators_3d
void test_manual_operators_3d()
Definition: testLinearStructure.cpp:336
domain
Domain domain
Definition: testProjection.cpp:88
dim
unsigned int dim(const Vector &z)
Definition: viewDualSurface.cpp:174