DGtal  0.9.2
exampleDiscreteExteriorCalculusSolve.cpp
1 #include <string>
3 
4 #include <QApplication>
5 
6 #include "DECExamplesCommon.h"
7 
8 // always include EigenSupport.h before any other Eigen headers
9 #include "DGtal/math/linalg/EigenSupport.h"
10 #include "DGtal/dec/DiscreteExteriorCalculus.h"
11 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
12 #include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
13 
14 #include "DGtal/io/viewers/Viewer3D.h"
15 #include "DGtal/io/boards/Board2D.h"
16 #include "DGtal/io/readers/GenericReader.h"
17 
18 using namespace DGtal;
19 using namespace std;
20 
21 void solve2d_laplace()
22 {
23  trace.beginBlock("2d discrete exterior calculus solve laplace");
24 
25  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(9,9));
26 
27  // create discrete exterior calculus from set
31  Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain));
33  trace.info() << calculus << endl;
34 
36  Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>() + 0.01 * calculus.identity<0, DUAL>();
38  trace.info() << "laplace = " << laplace << endl;
39 
41  Calculus::DualForm0 dirac(calculus);
42  dirac.myContainer(calculus.getCellIndex( calculus.myKSpace.uSpel(Z2i::Point(2,5))) ) = 1;
44 
45  {
46  Board2D board;
47  board << domain;
48  board << dirac;
49  board.saveSVG("solve_laplace_calculus.svg");
50  }
51 
52  { // simplicial llt
53  trace.beginBlock("simplicial llt");
54 
56  typedef EigenLinearAlgebraBackend::SolverSimplicialLLT LinearAlgebraSolver;
58 
59  Solver solver;
60  solver.compute(laplace);
61  Calculus::DualForm0 solution = solver.solve(dirac);
62 
63  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
65  trace.info() << solution << endl;
66  trace.endBlock();
67 
68  Board2D board;
69  board << domain;
70  board << solution;
71  board.saveSVG("solve_laplace_simplicial_llt.svg");
72  }
73 
74  { // simplicial ldlt
75  trace.beginBlock("simplicial ldlt");
76 
78  typedef EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver;
80 
81  Solver solver;
82  solver.compute(laplace);
83  Calculus::DualForm0 solution = solver.solve(dirac);
84 
85  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
87  trace.info() << solution << endl;
88  trace.endBlock();
89 
90  Board2D board;
91  board << domain;
92  board << solution;
93  board.saveSVG("solve_laplace_simplicial_ldlt.svg");
94  }
95 
96  { // conjugate gradient
97  trace.beginBlock("conjugate gradient");
98 
100  typedef EigenLinearAlgebraBackend::SolverConjugateGradient LinearAlgebraSolver;
102 
103  Solver solver;
104  solver.compute(laplace);
105  Calculus::DualForm0 solution = solver.solve(dirac);
107 
108  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
109  trace.info() << solution << endl;
110  trace.endBlock();
111 
112  Board2D board;
113  board << domain;
114  board << solution;
115  board.saveSVG("solve_laplace_conjugate_gradient.svg");
116  }
117 
118  { // biconjugate gradient stabilized
119  trace.beginBlock("biconjugate gradient stabilized (bicgstab)");
120 
122  typedef EigenLinearAlgebraBackend::SolverBiCGSTAB LinearAlgebraSolver;
124 
125  Solver solver;
126  solver.compute(laplace);
127  Calculus::DualForm0 solution = solver.solve(dirac);
129 
130  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
131  trace.info() << solution << endl;
132  trace.endBlock();
133 
134  Board2D board;
135  board << domain;
136  board << solution;
137  board.saveSVG("solve_laplace_bicgstab.svg");
138  }
139 
140  { // sparselu
141  trace.beginBlock("sparse lu");
142 
144  typedef EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;
146 
147  Solver solver;
148  solver.compute(laplace);
149  Calculus::DualForm0 solution = solver.solve(dirac);
151 
152  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
153  trace.info() << solution << endl;
154  trace.endBlock();
155 
156  Board2D board;
157  board << domain;
158  board << solution;
159  board.saveSVG("solve_laplace_sparse_lu.svg");
160  }
161 
162  { // sparseqr
163  trace.beginBlock("sparse qr");
164 
166  typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
168 
169  Solver solver;
170  solver.compute(-laplace);
171  Calculus::DualForm0 solution = -solver.solve(dirac);
173 
174  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
175  trace.info() << solution << endl;
176  trace.endBlock();
177 
178  Board2D board;
179  board << domain;
180  board << solution;
181  board.saveSVG("solve_laplace_sparse_qr.svg");
182  }
183 
184  trace.endBlock();
185 }
186 
187 void solve2d_dual_decomposition()
188 {
189  trace.beginBlock("2d discrete exterior calculus solve dual helmoltz decomposition");
190 
191  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(44,29));
192 
193  // create discrete exterior calculus from set
196  Calculus calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(domain));
197  trace.info() << calculus << endl;
198 
199  // choose linear solver
200  typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
201 
203  const Calculus::DualDerivative0 d0 = calculus.derivative<0, DUAL>();
204  const Calculus::DualDerivative1 d1 = calculus.derivative<1, DUAL>();
205  const Calculus::PrimalDerivative0 d0p = calculus.derivative<0, PRIMAL>();
206  const Calculus::PrimalDerivative1 d1p = calculus.derivative<1, PRIMAL>();
207  const Calculus::DualHodge1 h1 = calculus.hodge<1, DUAL>();
208  const Calculus::DualHodge2 h2 = calculus.hodge<2, DUAL>();
209  const Calculus::PrimalHodge1 h1p = calculus.hodge<1, PRIMAL>();
210  const Calculus::PrimalHodge2 h2p = calculus.hodge<2, PRIMAL>();
211  const LinearOperator<Calculus, 1, DUAL, 0, DUAL> ad1 = h2p * d1p * h1;
212  const LinearOperator<Calculus, 2, DUAL, 1, DUAL> ad2 = h1p * d0p * h2;
214 
216  Calculus::DualVectorField input_vector_field(calculus);
217  for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
218  {
219  const Z2i::RealPoint cell_center = Z2i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;
220  input_vector_field.myCoordinates(ii, 0) = cos(-.5*cell_center[0]+ .3*cell_center[1]);
221  input_vector_field.myCoordinates(ii, 1) = cos(.4*cell_center[0]+ .8*cell_center[1]);
222  }
223  trace.info() << input_vector_field << endl;
224 
225  const Calculus::DualForm1 input_one_form = calculus.flat(input_vector_field);
226  const Calculus::DualForm0 input_one_form_anti_derivated = ad1 * input_one_form;
227  const Calculus::DualForm2 input_one_form_derivated = d1 * input_one_form;
229 
230  {
231  Board2D board;
232  board << domain;
233  board << calculus;
234  board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
235  board << input_one_form;
236  board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));
237  board << input_vector_field;
238  board.saveSVG("solve_2d_dual_decomposition_calculus.svg");
239  }
240 
241 
242  Calculus::DualForm0 solution_curl_free(calculus);
243  { // solve curl free problem
244  trace.beginBlock("solving curl free component");
245 
248  Solver solver;
249  solver.compute(ad1 * d0);
250  solution_curl_free = solver.solve(input_one_form_anti_derivated);
252 
253  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
254  trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;
255  trace.endBlock();
256  }
257 
258  {
259  Board2D board;
260  board << domain;
261  board << calculus;
262  board << solution_curl_free;
263  board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));
264  board << calculus.sharp(d0*solution_curl_free);
265  board.saveSVG("solve_2d_dual_decomposition_curl_free.svg");
266  }
267 
268  Calculus::DualForm2 solution_div_free(calculus);
269  { // solve divergence free problem
270  trace.beginBlock("solving divergence free component");
271 
274  Solver solver;
275  solver.compute(d1 * ad2);
276  solution_div_free = solver.solve(input_one_form_derivated);
278 
279  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
280  trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;
281  trace.endBlock();
282  }
283 
284  {
285  Board2D board;
286  board << domain;
287  board << calculus;
288  board << solution_div_free;
289  board << CustomStyle("VectorField", new VectorFieldStyle2D(1.5));
290  board << calculus.sharp(ad2*solution_div_free);
291  board.saveSVG("solve_2d_dual_decomposition_div_free.svg");
292  }
293 
295  const Calculus::DualForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;
297  trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;
298 
299  {
300  Board2D board;
301  board << domain;
302  board << calculus;
303  board << solution_harmonic;
304  board << CustomStyle("VectorField", new VectorFieldStyle2D(20));
305  board << calculus.sharp(solution_harmonic);
306  board.saveSVG("solve_2d_dual_decomposition_harmonic.svg");
307  }
308 
309  trace.endBlock();
310 }
311 
312 void solve2d_primal_decomposition()
313 {
314  trace.beginBlock("2d discrete exterior calculus solve primal helmoltz decomposition");
315 
316  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(44,29));
317 
318  // create discrete exterior calculus from set
321  Calculus calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(domain));
322  trace.info() << calculus << endl;
323 
324  // choose linear solver
325  typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
326 
328  const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
329  const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();
330  const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
331  const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();
332  const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
333  const Calculus::PrimalHodge2 h2 = calculus.hodge<2, PRIMAL>();
334  const Calculus::DualHodge1 h1p = calculus.hodge<1, DUAL>();
335  const Calculus::DualHodge2 h2p = calculus.hodge<2, DUAL>();
336  const LinearOperator<Calculus, 1, PRIMAL, 0, PRIMAL> ad1 = h2p * d1p * h1;
337  const LinearOperator<Calculus, 2, PRIMAL, 1, PRIMAL> ad2 = h1p * d0p * h2;
339 
341  Calculus::PrimalVectorField input_vector_field(calculus);
342  for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
343  {
344  const Z2i::RealPoint cell_center = Z2i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;
345  input_vector_field.myCoordinates(ii, 0) = cos(-.5*cell_center[0]+ .3*cell_center[1]);
346  input_vector_field.myCoordinates(ii, 1) = cos(.4*cell_center[0]+ .8*cell_center[1]);
347  }
348  trace.info() << input_vector_field << endl;
349 
350  const Calculus::PrimalForm1 input_one_form = calculus.flat(input_vector_field);
351  const Calculus::PrimalForm0 input_one_form_anti_derivated = ad1 * input_one_form;
352  const Calculus::PrimalForm2 input_one_form_derivated = d1 * input_one_form;
354 
355  {
356  Board2D board;
357  board << domain;
358  board << calculus;
359  board << input_one_form;
360  board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));
361  board << input_vector_field;
362  board.saveSVG("solve_2d_primal_decomposition_calculus.svg");
363  }
364 
365 
366  Calculus::PrimalForm0 solution_curl_free(calculus);
367  { // solve curl free problem
368  trace.beginBlock("solving curl free component");
369 
372  Solver solver;
373  solver.compute(ad1 * d0);
374  solution_curl_free = solver.solve(input_one_form_anti_derivated);
376 
377  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
378  trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;
379  trace.endBlock();
380  }
381 
382  {
383  Board2D board;
384  board << domain;
385  board << calculus;
386  board << solution_curl_free;
387  board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));
388  board << calculus.sharp(d0*solution_curl_free);
389  board.saveSVG("solve_2d_primal_decomposition_curl_free.svg");
390  }
391 
392  Calculus::PrimalForm2 solution_div_free(calculus);
393  { // solve divergence free problem
394  trace.beginBlock("solving divergence free component");
395 
398  Solver solver;
399  solver.compute(d1 * ad2);
400  solution_div_free = solver.solve(input_one_form_derivated);
402 
403  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
404  trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;
405  trace.endBlock();
406  }
407 
408  {
409  Board2D board;
410  board << domain;
411  board << calculus;
412  board << solution_div_free;
413  board << CustomStyle("VectorField", new VectorFieldStyle2D(1.5));
414  board << calculus.sharp(ad2*solution_div_free);
415  board.saveSVG("solve_2d_primal_decomposition_div_free.svg");
416  }
417 
419  const Calculus::PrimalForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;
421  trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;
422 
423  {
424  Board2D board;
425  board << domain;
426  board << calculus;
427  board << solution_harmonic;
428  board << CustomStyle("VectorField", new VectorFieldStyle2D(30));
429  board << calculus.sharp(solution_harmonic);
430  board.saveSVG("solve_2d_primal_decomposition_harmonic.svg");
431  }
432 
433  trace.endBlock();
434 }
435 
436 void solve3d_decomposition()
437 {
438  trace.beginBlock("3d discrete exterior calculus solve helmoltz decomposition");
439 
440  const Z3i::Domain domain(Z3i::Point(0,0,0), Z3i::Point(19,19,9));
441 
442  // choose linear solver
443  typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
444 
446  // create discrete exterior calculus from set
448  Calculus calculus;
449  calculus.initKSpace<Z3i::Domain>(domain);
450 
451  // outer ring
452  for (int kk=2; kk<=18; kk++)
453  for (int ll=4; ll<=36; ll++)
454  {
455  { // bottom
456  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4,kk));
457  const Dimension dim = calculus.myKSpace.uDim(cell);
458  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
459  switch (dim)
460  {
461  case 0:
462  sign = Calculus::KSpace::POS;
463  break;
464  case 1:
465  sign = Calculus::KSpace::NEG;
466  break;
467  case 2:
468  sign = Calculus::KSpace::NEG;
469  break;
470  default:
471  break;
472  }
473  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
474  }
475 
476  { // top
477  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,36,kk));
478  const Dimension dim = calculus.myKSpace.uDim(cell);
479  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
480  switch (dim)
481  {
482  case 0:
483  sign = Calculus::KSpace::POS;
484  break;
485  case 1:
486  sign = Calculus::KSpace::POS;
487  break;
488  case 2:
489  sign = Calculus::KSpace::POS;
490  break;
491  default:
492  break;
493  }
494  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
495  }
496 
497  { // left
498  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4,ll,kk));
499  const Dimension dim = calculus.myKSpace.uDim(cell);
500  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
501  switch (dim)
502  {
503  case 0:
504  sign = Calculus::KSpace::POS;
505  break;
506  case 1:
507  sign = ( *calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );
508  break;
509  case 2:
510  sign = Calculus::KSpace::NEG;
511  break;
512  default:
513  break;
514  }
515  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
516  }
517 
518  { // right
519  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(36,ll,kk));
520  const Dimension dim = calculus.myKSpace.uDim(cell);
521  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
522  switch (dim)
523  {
524  case 0:
525  sign = Calculus::KSpace::POS;
526  break;
527  case 1:
528  sign = ( *calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );
529  break;
530  case 2:
531  sign = Calculus::KSpace::POS;
532  break;
533  default:
534  break;
535  }
536  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
537  }
538 
539  }
540 
541  // inner ring
542  for (int kk=2; kk<=18; kk++)
543  for (int ll=16; ll<=24; ll++)
544  {
545  { // bottom
546  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,16,kk));
547  const Dimension dim = calculus.myKSpace.uDim(cell);
548  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
549  switch (dim)
550  {
551  case 0:
552  sign = Calculus::KSpace::POS;
553  break;
554  case 1:
555  sign = ( *calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );
556  break;
557  case 2:
558  sign = Calculus::KSpace::POS;
559  break;
560  default:
561  break;
562  }
563  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
564  }
565 
566  { // top
567  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24,kk));
568  const Dimension dim = calculus.myKSpace.uDim(cell);
569  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
570  switch (dim)
571  {
572  case 0:
573  sign = Calculus::KSpace::POS;
574  break;
575  case 1:
576  sign = ( *calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );
577  break;
578  case 2:
579  sign = Calculus::KSpace::NEG;
580  break;
581  default:
582  break;
583  }
584  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
585  }
586 
587  { // left
588  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(16,ll,kk));
589  const Dimension dim = calculus.myKSpace.uDim(cell);
590  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
591  switch (dim)
592  {
593  case 0:
594  sign = Calculus::KSpace::POS;
595  break;
596  case 1:
597  sign = Calculus::KSpace::POS;
598  break;
599  case 2:
600  sign = Calculus::KSpace::POS;
601  break;
602  default:
603  break;
604  }
605  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
606  }
607 
608  { // right
609  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24,ll,kk));
610  const Dimension dim = calculus.myKSpace.uDim(cell);
611  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
612  switch (dim)
613  {
614  case 0:
615  sign = Calculus::KSpace::POS;
616  break;
617  case 1:
618  sign = Calculus::KSpace::NEG;
619  break;
620  case 2:
621  sign = Calculus::KSpace::NEG;
622  break;
623  default:
624  break;
625  }
626  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
627  }
628  }
629 
630  // back and front
631  for (int kk=4; kk<=36; kk++)
632  for (int ll=0; ll<=12; ll++)
633  {
634  { // back
635  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4+ll,kk,2));
636  const Dimension dim = calculus.myKSpace.uDim(cell);
637  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
638  switch (dim)
639  {
640  case 0:
641  sign = Calculus::KSpace::POS;
642  break;
643  case 1:
644  sign = Calculus::KSpace::POS;
645  break;
646  case 2:
647  sign = Calculus::KSpace::NEG;
648  break;
649  default:
650  break;
651  }
652  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
653  }
654 
655  { // front
656  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4+ll,kk,18));
657  const Dimension dim = calculus.myKSpace.uDim(cell);
658  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
659  switch (dim)
660  {
661  case 0:
662  sign = Calculus::KSpace::POS;
663  break;
664  case 1:
665  sign = Calculus::KSpace::NEG;
666  break;
667  case 2:
668  sign = Calculus::KSpace::POS;
669  break;
670  default:
671  break;
672  }
673  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
674  }
675 
676  { // back
677  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24+ll,kk,2));
678  const Dimension dim = calculus.myKSpace.uDim(cell);
679  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
680  switch (dim)
681  {
682  case 0:
683  sign = Calculus::KSpace::POS;
684  break;
685  case 1:
686  sign = Calculus::KSpace::POS;
687  break;
688  case 2:
689  sign = Calculus::KSpace::NEG;
690  break;
691  default:
692  break;
693  }
694  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
695  }
696 
697  { // front
698  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24+ll,kk,18));
699  const Dimension dim = calculus.myKSpace.uDim(cell);
700  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
701  switch (dim)
702  {
703  case 0:
704  sign = Calculus::KSpace::POS;
705  break;
706  case 1:
707  sign = Calculus::KSpace::NEG;
708  break;
709  case 2:
710  sign = Calculus::KSpace::POS;
711  break;
712  default:
713  break;
714  }
715  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
716  }
717  }
718 
719  // back and front
720  for (int kk=0; kk<=12; kk++)
721  for (int ll=16; ll<=24; ll++)
722  {
723  { // back
724  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4+kk,2));
725  const Dimension dim = calculus.myKSpace.uDim(cell);
726  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
727  switch (dim)
728  {
729  case 0:
730  sign = Calculus::KSpace::POS;
731  break;
732  case 1:
733  sign = Calculus::KSpace::POS;
734  break;
735  case 2:
736  sign = Calculus::KSpace::NEG;
737  break;
738  default:
739  break;
740  }
741  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
742  }
743 
744  { // front
745  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4+kk,18));
746  const Dimension dim = calculus.myKSpace.uDim(cell);
747  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
748  switch (dim)
749  {
750  case 0:
751  sign = Calculus::KSpace::POS;
752  break;
753  case 1:
754  sign = Calculus::KSpace::NEG;
755  break;
756  case 2:
757  sign = Calculus::KSpace::POS;
758  break;
759  default:
760  break;
761  }
762  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
763  }
764 
765  { // back
766  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24+kk,2));
767  const Dimension dim = calculus.myKSpace.uDim(cell);
768  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
769  switch (dim)
770  {
771  case 0:
772  sign = Calculus::KSpace::POS;
773  break;
774  case 1:
775  sign = Calculus::KSpace::POS;
776  break;
777  case 2:
778  sign = Calculus::KSpace::NEG;
779  break;
780  default:
781  break;
782  }
783  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
784  }
785 
786  { // front
787  const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24+kk,18));
788  const Dimension dim = calculus.myKSpace.uDim(cell);
789  Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
790  switch (dim)
791  {
792  case 0:
793  sign = Calculus::KSpace::POS;
794  break;
795  case 1:
796  sign = Calculus::KSpace::NEG;
797  break;
798  case 2:
799  sign = Calculus::KSpace::POS;
800  break;
801  default:
802  break;
803  }
804  calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );
805  }
806  }
807 
808  calculus.updateIndexes();
810 
811  trace.info() << calculus << endl;
812 
813  {
815  Viewer* viewer = new Viewer(calculus.myKSpace);
816  viewer->show();
817  viewer->setWindowTitle("structure");
818  (*viewer) << CustomColors3D(DGtal::Color(255,0,0), DGtal::Color(0,0,0));
819  (*viewer) << domain;
821  (*viewer) << Viewer::updateDisplay;
822  }
823 
825  const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
826  const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();
827  const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();
828  const Calculus::DualDerivative2 d2p = calculus.derivative<2, DUAL>();
829  const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
830  const Calculus::PrimalHodge2 h2 = calculus.hodge<2, PRIMAL>();
831  const Calculus::DualHodge2 h2p = calculus.hodge<2, DUAL>();
832  const Calculus::DualHodge3 h3p = calculus.hodge<3, DUAL>();
833  const LinearOperator<Calculus, 1, PRIMAL, 0, PRIMAL> ad1 = h3p * d2p * h1;
834  const LinearOperator<Calculus, 2, PRIMAL, 1, PRIMAL> ad2 = h2p * d1p * h2;
836 
837  {
838  const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
839  const Eigen::VectorXd laplace_diag = laplace.myContainer.diagonal();
840 
841  boost::array<int, 7> degrees;
842  std::fill(degrees.begin(), degrees.end(), 0);
843  for (int kk=0; kk<laplace_diag.rows(); kk++)
844  {
845  const int degree = laplace_diag[kk];
846  ASSERT( degree >= 0 );
847  ASSERT( degree < degrees.size() );
848  degrees[degree] ++;
849  }
850 
851  trace.info() << "node degrees" << endl;
852  for (int kk=0; kk<7; kk++)
853  trace.info() << kk << " " << degrees[kk] << endl;
854  }
855 
857  Calculus::PrimalVectorField input_vector_field(calculus);
858  for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
859  {
860  const Z3i::RealPoint cell_center = Z3i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;
861  input_vector_field.myCoordinates(ii, 0) = -cos(-.3*cell_center[0] + .6*cell_center[1] + .8*cell_center[2]);
862  input_vector_field.myCoordinates(ii, 1) = sin(.8*cell_center[0] + .3*cell_center[1] - .4*cell_center[2]);
863  input_vector_field.myCoordinates(ii, 2) = -cos(cell_center[2]*.5);
864  }
865  trace.info() << input_vector_field << endl;
866 
867  const Calculus::PrimalForm1 input_one_form = calculus.flat(input_vector_field);
869  const Calculus::PrimalForm0 input_one_form_anti_derivated = ad1 * input_one_form;
870  const Calculus::PrimalForm2 input_one_form_derivated = d1 * input_one_form;
871 
872  {
873  typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
874  Viewer* viewer = new Viewer(calculus.myKSpace);
875  viewer->show();
876  viewer->setWindowTitle("input vector field");
877  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form);
878  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form_derivated);
879  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form_anti_derivated);
880  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_vector_field);
881  (*viewer) << Viewer::updateDisplay;
882  }
883 
884  Calculus::PrimalForm0 solution_curl_free(calculus);
885  { // solve curl free problem
886  trace.beginBlock("solving curl free component");
887 
890  Solver solver;
891  solver.compute(ad1 * d0);
892  solution_curl_free = solver.solve(input_one_form_anti_derivated);
894 
895  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
896  trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;
897  trace.endBlock();
898  }
899 
900  {
901  typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
902  Viewer* viewer = new Viewer(calculus.myKSpace);
903  viewer->show();
904  viewer->setWindowTitle("curl free solution");
905  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_curl_free);
906  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(d0*solution_curl_free));
907  (*viewer) << Viewer::updateDisplay;
908  }
909 
910  Calculus::PrimalForm2 solution_div_free(calculus);
911  { // solve divergence free problem
912  trace.beginBlock("solving divergence free component");
913 
916  Solver solver;
917  solver.compute(d1 * ad2);
918  solution_div_free = solver.solve(input_one_form_derivated);
920 
921  trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
922  trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;
923  trace.endBlock();
924  }
925 
926  {
927  typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
928  Viewer* viewer = new Viewer(calculus.myKSpace);
929  viewer->show();
930  viewer->setWindowTitle("div free solution");
931  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_div_free);
932  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(ad2*solution_div_free));
933  (*viewer) << Viewer::updateDisplay;
934  }
935 
937  const Calculus::PrimalForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;
939  trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;
940 
941  {
942  typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;
943  Viewer* viewer = new Viewer(calculus.myKSpace);
944  viewer->show();
945  viewer->setWindowTitle("harmonic");
946  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_harmonic);
947  Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(solution_harmonic), 10.);
948  (*viewer) << Viewer::updateDisplay;
949  }
950 
951  trace.endBlock();
952 }
953 
954 int main(int argc, char* argv[])
955 {
956  QApplication app(argc,argv);
957 
958  solve2d_laplace();
959  solve2d_dual_decomposition();
960  solve2d_primal_decomposition();
961  solve3d_decomposition();
962 
963  return app.exec();
964 }
965 
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Trace trace
Definition: Common.h:130
SolutionKForm solve(const InputKForm &input_kform) const
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
DGtal::uint32_t Dimension
Definition: Common.h:113
Eigen::SparseQR< SparseMatrix, Eigen::COLAMDOrdering< SparseMatrix::Index > > SolverSparseQR
Definition: EigenSupport.h:106
STL namespace.
double endBlock()
Eigen::ConjugateGradient< SparseMatrix > SolverConjugateGradient
Definition: EigenSupport.h:103
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
Eigen::SimplicialLLT< SparseMatrix > SolverSimplicialLLT
Definition: EigenSupport.h:101
Aim: This class provides static members to create DEC structures from various other DGtal structures...
Eigen::BiCGSTAB< SparseMatrix > SolverBiCGSTAB
Definition: EigenSupport.h:104
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package...
DGtal is the top-level namespace which contains all DGtal functions and types.
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
std::ostream & info()
Eigen::SimplicialLDLT< SparseMatrix > SolverSimplicialLDLT
Definition: EigenSupport.h:102
void initKSpace(ConstAlias< TDomain > domain)
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
Space::RealPoint RealPoint
Definition: StdDefs.h:97
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
Space::RealPoint RealPoint
Definition: StdDefs.h:170
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70
Eigen::SparseLU< SparseMatrix > SolverSparseLU
Definition: EigenSupport.h:105