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