DGtal 1.4.0
Loading...
Searching...
No Matches
exampleDiscreteExteriorCalculusSolve.cpp
1
17
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
40using namespace DGtal;
41using namespace std;
42
43void 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
209void 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
334void 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
458void 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
976int 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::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:68
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::PolygonalCalculus::sharp
DenseMatrix sharp(const Face f) const
Definition PolygonalCalculus.h:503
DGtal::PolygonalCalculus::flat
DenseMatrix flat(const Face f) const
Definition PolygonalCalculus.h:462
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
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...
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::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::Dimension
DGtal::uint32_t Dimension
Definition Common.h:136
DGtal::trace
Trace trace
Definition Common.h:153
DGtal::PRIMAL
@ PRIMAL
Definition Duality.h:61
DGtal::DUAL
@ DUAL
Definition Duality.h:62
std
STL namespace.
DGtal::Display3DFactory::draw
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
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