doc/stable/dec_2exampleDiscreteExteriorCalculusSolve_8cpp-example.html

#include <string>


#include <QApplication>


#include "DECExamplesCommon.h"


// always include EigenSupport.h before any other Eigen headers

#include "DGtal/math/linalg/EigenSupport.h"

#include "DGtal/dec/DiscreteExteriorCalculus.h"

#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"

#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"


#include "DGtal/io/viewers/Viewer3D.h"

#include "DGtal/io/boards/Board2D.h"

#include "DGtal/io/readers/GenericReader.h"


using namespace DGtal;

using namespace std;


void solve2d_laplace()

{

    trace.beginBlock("2d discrete exterior calculus solve laplace");


    const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(9,9));


    // create discrete exterior calculus from set

    typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;

    typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;

    Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain));

    trace.info() << calculus << endl;


    Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>() + 0.01 * calculus.identity<0, DUAL>();

    trace.info() << "laplace = " << laplace << endl;


    Calculus::DualForm0 dirac(calculus);

    dirac.myContainer(calculus.getCellIndex( calculus.myKSpace.uSpel(Z2i::Point(2,5))) ) = 1;


    {

        Board2D board;

        board << domain;

        board << dirac;

        board.saveSVG("solve_laplace_calculus.svg");

    }


    { // simplicial llt

        trace.beginBlock("simplicial llt");


        typedef EigenLinearAlgebraBackend::SolverSimplicialLLT LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(laplace);

        Calculus::DualForm0 solution = solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_simplicial_llt.svg");

    }


    { // simplicial ldlt

        trace.beginBlock("simplicial ldlt");


        typedef EigenLinearAlgebraBackend::SolverSimplicialLDLT LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(laplace);

        Calculus::DualForm0 solution = solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_simplicial_ldlt.svg");

    }


    { // conjugate gradient

        trace.beginBlock("conjugate gradient");


        typedef EigenLinearAlgebraBackend::SolverConjugateGradient LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(laplace);

        Calculus::DualForm0 solution = solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_conjugate_gradient.svg");

    }


    { // biconjugate gradient stabilized

        trace.beginBlock("biconjugate gradient stabilized (bicgstab)");


        typedef EigenLinearAlgebraBackend::SolverBiCGSTAB LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(laplace);

        Calculus::DualForm0 solution = solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_bicgstab.svg");

    }


    { // sparselu

        trace.beginBlock("sparse lu");


        typedef EigenLinearAlgebraBackend::SolverSparseLU LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(laplace);

        Calculus::DualForm0 solution = solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_sparse_lu.svg");

    }


    { // sparseqr

        trace.beginBlock("sparse qr");


        typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;

        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;


        Solver solver;

        solver.compute(-laplace);

        Calculus::DualForm0 solution = -solver.solve(dirac);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << solution << endl;

        trace.endBlock();


        Board2D board;

        board << domain;

        board << solution;

        board.saveSVG("solve_laplace_sparse_qr.svg");

    }


    trace.endBlock();

}


void solve2d_dual_decomposition()

{

    trace.beginBlock("2d discrete exterior calculus solve dual helmoltz decomposition");


    const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(44,29));


    // create discrete exterior calculus from set

    typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;

    typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;

    Calculus calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(domain));

    trace.info() << calculus << endl;


    // choose linear solver

    typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;


    const Calculus::DualDerivative0 d0 = calculus.derivative<0, DUAL>();

    const Calculus::DualDerivative1 d1 = calculus.derivative<1, DUAL>();

    const Calculus::PrimalDerivative0 d0p = calculus.derivative<0, PRIMAL>();

    const Calculus::PrimalDerivative1 d1p = calculus.derivative<1, PRIMAL>();

    const Calculus::DualHodge1 h1 = calculus.hodge<1, DUAL>();

    const Calculus::DualHodge2 h2 = calculus.hodge<2, DUAL>();

    const Calculus::PrimalHodge1 h1p = calculus.hodge<1, PRIMAL>();

    const Calculus::PrimalHodge2 h2p = calculus.hodge<2, PRIMAL>();

    const LinearOperator<Calculus, 1, DUAL, 0, DUAL> ad1 = h2p * d1p * h1;

    const LinearOperator<Calculus, 2, DUAL, 1, DUAL> ad2 = h1p * d0p * h2;


    Calculus::DualVectorField input_vector_field(calculus);

    for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)

    {

        const Z2i::RealPoint cell_center = Z2i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;

        input_vector_field.myCoordinates(ii, 0) = cos(-.5*cell_center[0]+ .3*cell_center[1]);

        input_vector_field.myCoordinates(ii, 1) = cos(.4*cell_center[0]+ .8*cell_center[1]);

    }

    trace.info() << input_vector_field << endl;


    const Calculus::DualForm1 input_one_form = calculus.flat(input_vector_field);

    const Calculus::DualForm0 input_one_form_anti_derivated = ad1 * input_one_form;

    const Calculus::DualForm2 input_one_form_derivated = d1 * input_one_form;


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << CustomStyle("KForm", new KFormStyle2D(-1, 1));

        board << input_one_form;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));

        board << input_vector_field;

        board.saveSVG("solve_2d_dual_decomposition_calculus.svg");

    }


    Calculus::DualForm0 solution_curl_free(calculus);

    { // solve curl free problem

        trace.beginBlock("solving curl free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, DUAL, 0, DUAL> Solver;

        Solver solver;

        solver.compute(ad1 * d0);

        solution_curl_free = solver.solve(input_one_form_anti_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_curl_free;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));

        board << calculus.sharp(d0*solution_curl_free);

        board.saveSVG("solve_2d_dual_decomposition_curl_free.svg");

    }


    Calculus::DualForm2 solution_div_free(calculus);

    { // solve divergence free problem

        trace.beginBlock("solving divergence free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 2, DUAL, 2, DUAL> Solver;

        Solver solver;

        solver.compute(d1 * ad2);

        solution_div_free = solver.solve(input_one_form_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_div_free;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(1.5));

        board << calculus.sharp(ad2*solution_div_free);

        board.saveSVG("solve_2d_dual_decomposition_div_free.svg");

    }


    const Calculus::DualForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;

    trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_harmonic;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(20));

        board << calculus.sharp(solution_harmonic);

        board.saveSVG("solve_2d_dual_decomposition_harmonic.svg");

    }


    trace.endBlock();

}


void solve2d_primal_decomposition()

{

    trace.beginBlock("2d discrete exterior calculus solve primal helmoltz decomposition");


    const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(44,29));


    // create discrete exterior calculus from set

    typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;

    typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;

    Calculus calculus = CalculusFactory::createFromDigitalSet(generateDoubleRingSet(domain));

    trace.info() << calculus << endl;


    // choose linear solver

    typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;


    const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();

    const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();

    const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();

    const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();

    const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();

    const Calculus::PrimalHodge2 h2 = calculus.hodge<2, PRIMAL>();

    const Calculus::DualHodge1 h1p = calculus.hodge<1, DUAL>();

    const Calculus::DualHodge2 h2p = calculus.hodge<2, DUAL>();

    const LinearOperator<Calculus, 1, PRIMAL, 0, PRIMAL> ad1 = h2p * d1p * h1;

    const LinearOperator<Calculus, 2, PRIMAL, 1, PRIMAL> ad2 = h1p * d0p * h2;


    Calculus::PrimalVectorField input_vector_field(calculus);

    for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)

    {

        const Z2i::RealPoint cell_center = Z2i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;

        input_vector_field.myCoordinates(ii, 0) = cos(-.5*cell_center[0]+ .3*cell_center[1]);

        input_vector_field.myCoordinates(ii, 1) = cos(.4*cell_center[0]+ .8*cell_center[1]);

    }

    trace.info() << input_vector_field << endl;


    const Calculus::PrimalForm1 input_one_form = calculus.flat(input_vector_field);

    const Calculus::PrimalForm0 input_one_form_anti_derivated = ad1 * input_one_form;

    const Calculus::PrimalForm2 input_one_form_derivated = d1 * input_one_form;


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << input_one_form;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));

        board << input_vector_field;

        board.saveSVG("solve_2d_primal_decomposition_calculus.svg");

    }


    Calculus::PrimalForm0 solution_curl_free(calculus);

    { // solve curl free problem

        trace.beginBlock("solving curl free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;

        Solver solver;

        solver.compute(ad1 * d0);

        solution_curl_free = solver.solve(input_one_form_anti_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_curl_free;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(.75));

        board << calculus.sharp(d0*solution_curl_free);

        board.saveSVG("solve_2d_primal_decomposition_curl_free.svg");

    }


    Calculus::PrimalForm2 solution_div_free(calculus);

    { // solve divergence free problem

        trace.beginBlock("solving divergence free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL> Solver;

        Solver solver;

        solver.compute(d1 * ad2);

        solution_div_free = solver.solve(input_one_form_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_div_free;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(1.5));

        board << calculus.sharp(ad2*solution_div_free);

        board.saveSVG("solve_2d_primal_decomposition_div_free.svg");

    }


    const Calculus::PrimalForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;

    trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;


    {

        Board2D board;

        board << domain;

        board << calculus;

        board << solution_harmonic;

        board << CustomStyle("VectorField", new VectorFieldStyle2D(30));

        board << calculus.sharp(solution_harmonic);

        board.saveSVG("solve_2d_primal_decomposition_harmonic.svg");

    }


    trace.endBlock();

}


void solve3d_decomposition()

{

    trace.beginBlock("3d discrete exterior calculus solve helmoltz decomposition");


    const Z3i::Domain domain(Z3i::Point(0,0,0), Z3i::Point(19,19,9));


    // choose linear solver

    typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;


    // create discrete exterior calculus from set

    typedef DiscreteExteriorCalculus<3, 3, EigenLinearAlgebraBackend> Calculus;

    Calculus calculus;

    calculus.initKSpace<Z3i::Domain>(domain);


    // outer ring

    for (int kk=2; kk<=18; kk++)

        for (int ll=4; ll<=36; ll++)

        {

            { // bottom

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // top

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,36,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // left

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4,ll,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = ( *calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // right

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(36,ll,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = ( *calculus.myKSpace.uDirs(cell) == 2 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


        }


    // inner ring

    for (int kk=2; kk<=18; kk++)

        for (int ll=16; ll<=24; ll++)

        {

            { // bottom

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,16,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = ( *calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS );

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // top

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = ( *calculus.myKSpace.uDirs(cell) == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG );

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // left

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(16,ll,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // right

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24,ll,kk));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }

        }


    // back and front

    for (int kk=4; kk<=36; kk++)

        for (int ll=0; ll<=12; ll++)

        {

            { // back

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4+ll,kk,2));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // front

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(4+ll,kk,18));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // back

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24+ll,kk,2));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // front

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(24+ll,kk,18));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }

        }


    // back and front

    for (int kk=0; kk<=12; kk++)

        for (int ll=16; ll<=24; ll++)

        {

            { // back

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4+kk,2));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // front

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,4+kk,18));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // back

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24+kk,2));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 2:

                        sign = Calculus::KSpace::NEG;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }


            { // front

                const Calculus::Cell cell = calculus.myKSpace.uCell(Z3i::Point(ll,24+kk,18));

                const Dimension dim = calculus.myKSpace.uDim(cell);

                Calculus::KSpace::Sign sign = Calculus::KSpace::POS;

                switch (dim)

                {

                    case 0:

                        sign = Calculus::KSpace::POS;

                        break;

                    case 1:

                        sign = Calculus::KSpace::NEG;

                        break;

                    case 2:

                        sign = Calculus::KSpace::POS;

                        break;

                    default:

                        break;

                }

                calculus.insertSCell( calculus.myKSpace.signs(cell, sign) );

            }

        }


    calculus.updateIndexes();


    trace.info() << calculus << endl;


    {

        typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;

        Viewer* viewer = new Viewer(calculus.myKSpace);

        viewer->show();

        viewer->setWindowTitle("structure");

        (*viewer) << CustomColors3D(DGtal::Color(255,0,0), DGtal::Color(0,0,0));

        (*viewer) << domain;

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus);

        (*viewer) << Viewer::updateDisplay;

    }


    const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();

    const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();

    const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();

    const Calculus::DualDerivative2 d2p = calculus.derivative<2, DUAL>();

    const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();

    const Calculus::PrimalHodge2 h2 = calculus.hodge<2, PRIMAL>();

    const Calculus::DualHodge2 h2p = calculus.hodge<2, DUAL>();

    const Calculus::DualHodge3 h3p = calculus.hodge<3, DUAL>();

    const LinearOperator<Calculus, 1, PRIMAL, 0, PRIMAL> ad1 = h3p * d2p * h1;

    const LinearOperator<Calculus, 2, PRIMAL, 1, PRIMAL> ad2 = h2p * d1p * h2;


    {

        const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();

        const Eigen::VectorXd laplace_diag = laplace.myContainer.diagonal();


        boost::array<int, 7> degrees;

        std::fill(degrees.begin(), degrees.end(), 0);

        for (int kk=0; kk<laplace_diag.rows(); kk++)

        {

            const int degree = laplace_diag[kk];

            ASSERT( degree >= 0 );

            ASSERT( static_cast<unsigned int>(degree) < degrees.size() );

            degrees[degree] ++;

        }


        trace.info() << "node degrees" << endl;

        for (int kk=0; kk<7; kk++)

            trace.info() << kk << " " << degrees[kk] << endl;

    }


    Calculus::PrimalVectorField input_vector_field(calculus);

    for (Calculus::Index ii=0; ii<input_vector_field.length(); ii++)

    {

        const Z3i::RealPoint cell_center = Z3i::RealPoint(input_vector_field.getSCell(ii).preCell().coordinates)/2.;

        input_vector_field.myCoordinates(ii, 0) = -cos(-.3*cell_center[0] + .6*cell_center[1] + .8*cell_center[2]);

        input_vector_field.myCoordinates(ii, 1) = sin(.8*cell_center[0] + .3*cell_center[1] - .4*cell_center[2]);

        input_vector_field.myCoordinates(ii, 2) = -cos(cell_center[2]*.5);

    }

    trace.info() << input_vector_field << endl;


    const Calculus::PrimalForm1 input_one_form = calculus.flat(input_vector_field);

    const Calculus::PrimalForm0 input_one_form_anti_derivated = ad1 * input_one_form;

    const Calculus::PrimalForm2 input_one_form_derivated = d1 * input_one_form;


    {

        typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;

        Viewer* viewer = new Viewer(calculus.myKSpace);

        viewer->show();

        viewer->setWindowTitle("input vector field");

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form_derivated);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_one_form_anti_derivated);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, input_vector_field);

        (*viewer) << Viewer::updateDisplay;

    }


    Calculus::PrimalForm0 solution_curl_free(calculus);

    { // solve curl free problem

        trace.beginBlock("solving curl free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;

        Solver solver;

        solver.compute(ad1 * d0);

        solution_curl_free = solver.solve(input_one_form_anti_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_curl_free.myContainer.minCoeff() << " max=" << solution_curl_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;

        Viewer* viewer = new Viewer(calculus.myKSpace);

        viewer->show();

        viewer->setWindowTitle("curl free solution");

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_curl_free);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(d0*solution_curl_free));

        (*viewer) << Viewer::updateDisplay;

    }


    Calculus::PrimalForm2 solution_div_free(calculus);

    { // solve divergence free problem

        trace.beginBlock("solving divergence free component");


        typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 2, PRIMAL, 2, PRIMAL> Solver;

        Solver solver;

        solver.compute(d1 * ad2);

        solution_div_free = solver.solve(input_one_form_derivated);


        trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;

        trace.info() << "min=" << solution_div_free.myContainer.minCoeff() << " max=" << solution_div_free.myContainer.maxCoeff() << endl;

        trace.endBlock();

    }


    {

        typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;

        Viewer* viewer = new Viewer(calculus.myKSpace);

        viewer->show();

        viewer->setWindowTitle("div free solution");

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_div_free);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(ad2*solution_div_free));

        (*viewer) << Viewer::updateDisplay;

    }


    const Calculus::PrimalForm1 solution_harmonic = input_one_form - d0*solution_curl_free - ad2*solution_div_free;

    trace.info() << "min=" << solution_harmonic.myContainer.minCoeff() << " max=" << solution_harmonic.myContainer.maxCoeff() << endl;


    {

        typedef Viewer3D<Z3i::Space, Z3i::KSpace> Viewer;

        Viewer* viewer = new Viewer(calculus.myKSpace);

        viewer->show();

        viewer->setWindowTitle("harmonic");

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, solution_harmonic);

        Display3DFactory<Z3i::Space, Z3i::KSpace>::draw(*viewer, calculus.sharp(solution_harmonic), 10.);

        (*viewer) << Viewer::updateDisplay;

    }


    trace.endBlock();

}


int main(int argc, char* argv[])

{

    QApplication app(argc,argv);


    solve2d_laplace();

    solve2d_dual_decomposition();

    solve2d_primal_decomposition();

    solve3d_decomposition();


    return app.exec();

}

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::DiscreteExteriorCalculus::initKSpace
void initKSpace(ConstAlias< TDomain > domain)

DGtal::Display3D::updateDisplay
@ updateDisplay
Definition Display3D.h:249

DGtal::HyperRectDomain< Space >

DGtal::LinearOperator
Aim: LinearOperator represents discrete linear operator between discrete kforms in the DEC package.
Definition LinearOperator.h:69

DGtal::PointVector< dim, Integer >

DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")

DGtal::Trace::info
std::ostream & info()

DGtal::Trace::endBlock
double endBlock()

DGtal::Viewer3D
Definition Viewer3D.h:132

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
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::CustomColors3D
Definition DrawWithDisplay3DModifier.h:130

DGtal::CustomStyle
Definition Board2D.h:217

DGtal::Display3DFactory
Factory for GPL Display3D:
Definition Display3DFactory.h:87

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

dim
#define dim
Definition testImageContainerByHashTree.cpp:40

domain
Domain domain
Definition testProjection.cpp:88