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