DGtal  1.5.beta
dec/exampleDiscreteExteriorCalculusSolve.cpp

Example of primal and dual Helmoltz decomposition in 2D and 3D using Discrete Exterior Calculus.

See also
Helmoltz decomposition
Primal Helmoltz decomposition harmonic component.
#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::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
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...
calculus
PolyCalculus * calculus
Definition: dgtalCalculus-geodesic.cpp:70
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
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.
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
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(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
domain
Domain domain
Definition: testProjection.cpp:88
dim
unsigned int dim(const Vector &z)
Definition: viewDualSurface.cpp:174