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