DGtal  0.9.4beta
dec/exampleDiscreteExteriorCalculusUsage.cpp

Example of vector field manipulation using Discrete Exterior Calculus.

See also
Introduction
usage_primal_one_form_sharp_flat.png
Flat sharped primal gradient vector field.
#include <string>
#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/boards/Board2D.h"
#include "DGtal/io/readers/GenericReader.h"
using namespace std;
using namespace DGtal;
void usage2d()
{
trace.beginBlock("2d discrete exterior calculus usage");
const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(9,9));
// create discrete exterior calculus from set without border
{
Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain), false);
calculus.eraseCell(calculus.myKSpace.uSpel(Z2i::Point(8, 5)));
calculus.updateIndexes();
trace.info() << calculus << endl;
Board2D board;
board << domain;
board << calculus;
board.saveSVG("usage_calculus_without_border.svg");
}
// create discrete exterior calculus from set with border
Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain), true);
calculus.eraseCell(calculus.myKSpace.uSpel(Z2i::Point(8, 5)));
calculus.eraseCell(calculus.myKSpace.uCell(Z2i::Point(18, 11)));
calculus.updateIndexes();
trace.info() << calculus << endl;
{
Board2D board;
board << domain;
board << calculus;
board.saveSVG("usage_calculus_with_border.svg");
}
const Z2i::Point center(13,7);
// primal path
{
trace.info() << "primal path" << endl;
// create primal 0-form and fill it with euclidian metric
Calculus::PrimalForm0 primal_zero_form(calculus);
for (Calculus::Index index=0; index<primal_zero_form.length(); index++)
{
const Calculus::SCell& cell = primal_zero_form.getSCell(index);
const Calculus::Scalar& value = Z2i::l2Metric(calculus.myKSpace.sKCoords(cell), center)/2;
primal_zero_form.myContainer(index) = value;
}
// one can do linear algebra operation between equaly typed kforms
const Calculus::PrimalForm0 foo = 2 * primal_zero_form + primal_zero_form;
{
Board2D board;
board << domain;
board << calculus;
board << primal_zero_form;
board.saveSVG("usage_primal_zero_form.svg");
}
// create primal gradient vector field and primal derivative one form
const Calculus::PrimalDerivative0 primal_zero_derivative = calculus.derivative<0, PRIMAL>();
const Calculus::PrimalForm1 primal_one_form = primal_zero_derivative * primal_zero_form;
const Calculus::PrimalVectorField primal_vector_field = calculus.sharp(primal_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << primal_one_form;
board << primal_vector_field;
board.saveSVG("usage_primal_one_form.svg");
}
// test primal flat and sharp
const Calculus::PrimalForm1 flat_sharp_primal_one_form = calculus.flat(primal_vector_field);
const Calculus::PrimalVectorField sharp_flat_primal_vector_field = calculus.sharp(flat_sharp_primal_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << flat_sharp_primal_one_form;
board << sharp_flat_primal_vector_field;
board.saveSVG("usage_primal_one_form_sharp_flat.svg");
}
// create dual gradient vector field and hodge*d dual one form
const Calculus::PrimalHodge1 primal_one_hodge = calculus.hodge<1, PRIMAL>();
const Calculus::DualForm1 dual_one_form = primal_one_hodge * primal_zero_derivative * primal_zero_form;
const Calculus::DualVectorField dual_vector_field = calculus.sharp(dual_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << dual_one_form;
board << dual_vector_field;
board << primal_vector_field;
board.saveSVG("usage_primal_one_form_hodge.svg");
}
}
// dual path
{
trace.info() << "dual path" << endl;
// create dual 0-form and fill it with euclidian metric
Calculus::DualForm0 dual_zero_form(calculus);
for (Calculus::Index index=0; index<dual_zero_form.length(); index++)
{
const Calculus::SCell& cell = dual_zero_form.getSCell(index);
const Calculus::Scalar& value = Z2i::l2Metric(calculus.myKSpace.sKCoords(cell), center)/2;
dual_zero_form.myContainer(index) = value;
}
{
Board2D board;
board << domain;
board << calculus;
board << dual_zero_form;
board.saveSVG("usage_dual_zero_form.svg");
}
// create dual gradient vector field and dual derivative one form
const Calculus::DualDerivative0 dual_zero_derivative = calculus.derivative<0, DUAL>();
const Calculus::DualForm1 dual_one_form = dual_zero_derivative * dual_zero_form;
const Calculus::DualVectorField dual_vector_field = calculus.sharp(dual_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << dual_one_form;
board << dual_vector_field;
board.saveSVG("usage_dual_one_form.svg");
}
// test primal flat and sharp
const Calculus::DualForm1 flat_sharp_dual_one_form = calculus.flat(dual_vector_field);
const Calculus::DualVectorField sharp_flat_dual_vector_field = -calculus.sharp(flat_sharp_dual_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << flat_sharp_dual_one_form;
board << -sharp_flat_dual_vector_field;
board.saveSVG("usage_dual_one_form_sharp_flat.svg");
}
// create primal gradient vector field and hodge*d primal one form
const Calculus::DualHodge1 dual_one_hodge = calculus.hodge<1, DUAL>();
const Calculus::PrimalForm1 primal_one_form = dual_one_hodge * dual_zero_derivative * dual_zero_form;
const Calculus::PrimalVectorField primal_vector_field = calculus.sharp(primal_one_form);
{
Board2D board;
board << domain;
board << calculus;
board << primal_one_form;
board << primal_vector_field;
board << dual_vector_field;
board.saveSVG("usage_dual_one_form_hodge.svg");
}
}
}
int main()
{
usage2d();
return 0;
}