26#include "DECExamplesCommon.h"
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"
34#include "DGtal/io/viewers/PolyscopeViewer.h"
35#include "DGtal/io/boards/Board2D.h"
36#include "DGtal/io/readers/GenericReader.h"
43 trace.beginBlock(
"2d discrete exterior calculus solve laplace");
51 Calculus
calculus = CalculusFactory::createFromDigitalSet(generateRingSet(
domain));
58 trace.info() <<
"laplace = " << laplace << endl;
69 board.
saveSVG(
"solve_laplace_calculus.svg");
73 trace.beginBlock(
"simplicial llt");
81 Calculus::DualForm0 solution = solver.
solve(dirac);
83 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
85 trace.info() << solution << endl;
91 board.
saveSVG(
"solve_laplace_simplicial_llt.svg");
95 trace.beginBlock(
"simplicial ldlt");
103 Calculus::DualForm0 solution = solver.
solve(dirac);
105 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
107 trace.info() << solution << endl;
113 board.
saveSVG(
"solve_laplace_simplicial_ldlt.svg");
117 trace.beginBlock(
"conjugate gradient");
125 Calculus::DualForm0 solution = solver.
solve(dirac);
128 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
129 trace.info() << solution << endl;
135 board.
saveSVG(
"solve_laplace_conjugate_gradient.svg");
139 trace.beginBlock(
"biconjugate gradient stabilized (bicgstab)");
147 Calculus::DualForm0 solution = solver.
solve(dirac);
150 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
151 trace.info() << solution << endl;
157 board.
saveSVG(
"solve_laplace_bicgstab.svg");
161 trace.beginBlock(
"sparse lu");
169 Calculus::DualForm0 solution = solver.
solve(dirac);
172 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
173 trace.info() << solution << endl;
179 board.
saveSVG(
"solve_laplace_sparse_lu.svg");
183 trace.beginBlock(
"sparse qr");
191 Calculus::DualForm0 solution = -solver.
solve(dirac);
194 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
195 trace.info() << solution << endl;
201 board.
saveSVG(
"solve_laplace_sparse_qr.svg");
207void solve2d_dual_decomposition()
209 trace.beginBlock(
"2d discrete exterior calculus solve dual helmoltz decomposition");
216 Calculus
calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(
domain));
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>();
236 Calculus::DualVectorField input_vector_field(
calculus);
237 for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
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]);
243 trace.info() << input_vector_field << endl;
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;
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");
262 Calculus::DualForm0 solution_curl_free(
calculus);
264 trace.beginBlock(
"solving curl free component");
270 solution_curl_free = solver.
solve(input_one_form_anti_derivated);
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;
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");
288 Calculus::DualForm2 solution_div_free(
calculus);
290 trace.beginBlock(
"solving divergence free component");
296 solution_div_free = solver.
solve(input_one_form_derivated);
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;
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");
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;
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");
332void solve2d_primal_decomposition()
334 trace.beginBlock(
"2d discrete exterior calculus solve primal helmoltz decomposition");
341 Calculus
calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(
domain));
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>();
354 const Calculus::DualHodge1 h1p =
calculus.hodge<1,
DUAL>();
355 const Calculus::DualHodge2 h2p =
calculus.hodge<2,
DUAL>();
361 Calculus::PrimalVectorField input_vector_field(
calculus);
362 for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
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]);
368 trace.info() << input_vector_field << endl;
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;
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");
386 Calculus::PrimalForm0 solution_curl_free(
calculus);
388 trace.beginBlock(
"solving curl free component");
394 solution_curl_free = solver.
solve(input_one_form_anti_derivated);
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;
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");
412 Calculus::PrimalForm2 solution_div_free(
calculus);
414 trace.beginBlock(
"solving divergence free component");
420 solution_div_free = solver.
solve(input_one_form_derivated);
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;
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");
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;
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");
456void solve3d_decomposition()
458 trace.beginBlock(
"3d discrete exterior calculus solve helmoltz decomposition");
472 for (
int kk=2; kk<=18; kk++)
473 for (
int ll=4; ll<=36; ll++)
478 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
482 sign = Calculus::KSpace::POS;
485 sign = Calculus::KSpace::NEG;
488 sign = Calculus::KSpace::NEG;
499 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
503 sign = Calculus::KSpace::POS;
506 sign = Calculus::KSpace::POS;
509 sign = Calculus::KSpace::POS;
520 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
524 sign = Calculus::KSpace::POS;
527 sign = ( *
calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );
530 sign = Calculus::KSpace::NEG;
541 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
545 sign = Calculus::KSpace::POS;
548 sign = ( *
calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );
551 sign = Calculus::KSpace::POS;
562 for (
int kk=2; kk<=18; kk++)
563 for (
int ll=16; ll<=24; ll++)
568 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
572 sign = Calculus::KSpace::POS;
575 sign = ( *
calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );
578 sign = Calculus::KSpace::POS;
589 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
593 sign = Calculus::KSpace::POS;
596 sign = ( *
calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );
599 sign = Calculus::KSpace::NEG;
610 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
614 sign = Calculus::KSpace::POS;
617 sign = Calculus::KSpace::POS;
620 sign = Calculus::KSpace::POS;
631 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
635 sign = Calculus::KSpace::POS;
638 sign = Calculus::KSpace::NEG;
641 sign = Calculus::KSpace::NEG;
651 for (
int kk=4; kk<=36; kk++)
652 for (
int ll=0; ll<=12; ll++)
657 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
661 sign = Calculus::KSpace::POS;
664 sign = Calculus::KSpace::POS;
667 sign = Calculus::KSpace::NEG;
678 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
682 sign = Calculus::KSpace::POS;
685 sign = Calculus::KSpace::NEG;
688 sign = Calculus::KSpace::POS;
699 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
703 sign = Calculus::KSpace::POS;
706 sign = Calculus::KSpace::POS;
709 sign = Calculus::KSpace::NEG;
720 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
724 sign = Calculus::KSpace::POS;
727 sign = Calculus::KSpace::NEG;
730 sign = Calculus::KSpace::POS;
740 for (
int kk=0; kk<=12; kk++)
741 for (
int ll=16; ll<=24; ll++)
746 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
750 sign = Calculus::KSpace::POS;
753 sign = Calculus::KSpace::POS;
756 sign = Calculus::KSpace::NEG;
767 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
771 sign = Calculus::KSpace::POS;
774 sign = Calculus::KSpace::NEG;
777 sign = Calculus::KSpace::POS;
788 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
792 sign = Calculus::KSpace::POS;
795 sign = Calculus::KSpace::POS;
798 sign = Calculus::KSpace::NEG;
809 Calculus::KSpace::Sign sign = Calculus::KSpace::POS;
813 sign = Calculus::KSpace::POS;
816 sign = Calculus::KSpace::NEG;
819 sign = Calculus::KSpace::POS;
835 Viewer* viewer =
new Viewer(
calculus.myKSpace);
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>();
850 const Calculus::DualHodge2 h2p =
calculus.hodge<2,
DUAL>();
851 const Calculus::DualHodge3 h3p =
calculus.hodge<3,
DUAL>();
857 const Calculus::PrimalIdentity0 laplace =
calculus.laplace<
PRIMAL>();
858 const Eigen::VectorXd laplace_diag = laplace.myContainer.diagonal();
860 boost::array<int, 7> degrees;
861 std::fill(degrees.begin(), degrees.end(), 0);
862 for (
int kk=0; kk<laplace_diag.rows(); kk++)
864 const int degree = laplace_diag[kk];
865 ASSERT( degree >= 0 );
866 ASSERT(
static_cast<unsigned int>(degree) < degrees.size() );
870 trace.info() <<
"node degrees" << endl;
871 for (
int kk=0; kk<7; kk++)
872 trace.info() << kk <<
" " << degrees[kk] << endl;
876 Calculus::PrimalVectorField input_vector_field(
calculus);
877 for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)
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);
884 trace.info() << input_vector_field << endl;
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;
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;
902 Calculus::PrimalForm0 solution_curl_free(
calculus);
904 trace.beginBlock(
"solving curl free component");
910 solution_curl_free = solver.
solve(input_one_form_anti_derivated);
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;
920 Viewer* viewer =
new Viewer(
calculus.myKSpace);
921 (*viewer) << solution_curl_free;
922 (*viewer) <<
calculus.sharp(d0*solution_curl_free);
926 Calculus::PrimalForm2 solution_div_free(
calculus);
928 trace.beginBlock(
"solving divergence free component");
934 solution_div_free = solver.
solve(input_one_form_derivated);
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;
944 Viewer* viewer =
new Viewer(
calculus.myKSpace);
945 (*viewer) << solution_div_free;
946 (*viewer) <<
calculus.sharp(ad2*solution_div_free);
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;
957 Viewer* viewer =
new Viewer(
calculus.myKSpace);
958 (*viewer) << solution_harmonic;
959 (*viewer) <<
calculus.sharp(solution_harmonic);
966int main(
int argc,
char* argv[])
969 solve2d_dual_decomposition();
970 solve2d_primal_decomposition();
971 solve3d_decomposition();
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Structure representing an RGB triple with alpha component.
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
SolutionKForm solve(const InputKForm &input_kform) const
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
void show() override
Starts the event loop and display of elements.
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Space::RealPoint RealPoint
HyperRectDomain< Space > Domain
HyperRectDomain< Space > Domain
Space::RealPoint RealPoint
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Eigen::SimplicialLLT< SparseMatrix > SolverSimplicialLLT
Solvers on sparse matrices.
Eigen::SimplicialLDLT< SparseMatrix > SolverSimplicialLDLT
Eigen::BiCGSTAB< SparseMatrix > SolverBiCGSTAB
Eigen::SparseQR< SparseMatrix, Eigen::COLAMDOrdering< SparseMatrix::Index > > SolverSparseQR
Eigen::SparseLU< SparseMatrix > SolverSparseLU
Eigen::ConjugateGradient< SparseMatrix > SolverConjugateGradient