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