DGtal  0.9.2
exampleDiscreteExteriorCalculusUsage.cpp
1 #include <string>
3 
4 #include "DECExamplesCommon.h"
5 
7 // always include EigenSupport.h before any other Eigen headers
8 #include "DGtal/math/linalg/EigenSupport.h"
9 #include "DGtal/dec/DiscreteExteriorCalculus.h"
10 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
11 #include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
13 
14 #include "DGtal/io/boards/Board2D.h"
15 #include "DGtal/io/readers/GenericReader.h"
16 
17 using namespace std;
18 using namespace DGtal;
19 
20 void usage2d()
21 {
22  trace.beginBlock("2d discrete exterior calculus usage");
23 
24  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(9,9));
25 
30 
31  // create discrete exterior calculus from set without border
32  {
34  Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain), false);
35 
36  calculus.eraseCell(calculus.myKSpace.uSpel(Z2i::Point(8, 5)));
37 
38  calculus.updateIndexes();
40 
41  trace.info() << calculus << endl;
42 
43  Board2D board;
44  board << domain;
45  board << calculus;
46  board.saveSVG("usage_calculus_without_border.svg");
47  }
48 
49  // create discrete exterior calculus from set with border
51  Calculus calculus = CalculusFactory::createFromDigitalSet(generateRingSet(domain), true);
52 
53  calculus.eraseCell(calculus.myKSpace.uSpel(Z2i::Point(8, 5)));
54  calculus.eraseCell(calculus.myKSpace.uCell(Z2i::Point(18, 11)));
55 
56  calculus.updateIndexes();
58 
59  trace.info() << calculus << endl;
60 
61  {
62  Board2D board;
63  board << domain;
64  board << calculus;
65  board.saveSVG("usage_calculus_with_border.svg");
66  }
67 
68  const Z2i::Point center(13,7);
69 
70  // primal path
71  {
72  trace.info() << "primal path" << endl;
73 
74  // create primal 0-form and fill it with euclidian metric
76  Calculus::PrimalForm0 primal_zero_form(calculus);
77  for (Calculus::Index index=0; index<primal_zero_form.length(); index++)
78  {
79  const Calculus::SCell& cell = primal_zero_form.getSCell(index);
80  const Calculus::Scalar& value = Z2i::l2Metric(calculus.myKSpace.sKCoords(cell), center)/2;
81  primal_zero_form.myContainer(index) = value;
82  }
84  // one can do linear algebra operation between equaly typed kforms
86  const Calculus::PrimalForm0 foo = 2 * primal_zero_form + primal_zero_form;
88 
89  {
90  Board2D board;
91  board << domain;
92  board << calculus;
93  board << primal_zero_form;
94  board.saveSVG("usage_primal_zero_form.svg");
95  }
96 
97  // create primal gradient vector field and primal derivative one form
99  const Calculus::PrimalDerivative0 primal_zero_derivative = calculus.derivative<0, PRIMAL>();
100  const Calculus::PrimalForm1 primal_one_form = primal_zero_derivative * primal_zero_form;
101  const Calculus::PrimalVectorField primal_vector_field = calculus.sharp(primal_one_form);
103 
104  {
105  Board2D board;
106  board << domain;
107  board << calculus;
108  board << primal_one_form;
109  board << primal_vector_field;
110  board.saveSVG("usage_primal_one_form.svg");
111  }
112 
113  // test primal flat and sharp
115  const Calculus::PrimalForm1 flat_sharp_primal_one_form = calculus.flat(primal_vector_field);
116  const Calculus::PrimalVectorField sharp_flat_primal_vector_field = calculus.sharp(flat_sharp_primal_one_form);
118 
119  {
120  Board2D board;
121  board << domain;
122  board << calculus;
123  board << flat_sharp_primal_one_form;
124  board << sharp_flat_primal_vector_field;
125  board.saveSVG("usage_primal_one_form_sharp_flat.svg");
126  }
127 
128  // create dual gradient vector field and hodge*d dual one form
130  const Calculus::PrimalHodge1 primal_one_hodge = calculus.hodge<1, PRIMAL>();
131  const Calculus::DualForm1 dual_one_form = primal_one_hodge * primal_zero_derivative * primal_zero_form;
132  const Calculus::DualVectorField dual_vector_field = calculus.sharp(dual_one_form);
134 
135  {
136  Board2D board;
137  board << domain;
138  board << calculus;
139  board << dual_one_form;
140  board << dual_vector_field;
141  board << primal_vector_field;
142  board.saveSVG("usage_primal_one_form_hodge.svg");
143  }
144  }
145 
146  // dual path
147  {
148  trace.info() << "dual path" << endl;
149 
150  // create dual 0-form and fill it with euclidian metric
151  Calculus::DualForm0 dual_zero_form(calculus);
152  for (Calculus::Index index=0; index<dual_zero_form.length(); index++)
153  {
154  const Calculus::SCell& cell = dual_zero_form.getSCell(index);
155  const Calculus::Scalar& value = Z2i::l2Metric(calculus.myKSpace.sKCoords(cell), center)/2;
156  dual_zero_form.myContainer(index) = value;
157  }
158 
159  {
160  Board2D board;
161  board << domain;
162  board << calculus;
163  board << dual_zero_form;
164  board.saveSVG("usage_dual_zero_form.svg");
165  }
166 
167  // create dual gradient vector field and dual derivative one form
168  const Calculus::DualDerivative0 dual_zero_derivative = calculus.derivative<0, DUAL>();
169  const Calculus::DualForm1 dual_one_form = dual_zero_derivative * dual_zero_form;
170  const Calculus::DualVectorField dual_vector_field = calculus.sharp(dual_one_form);
171 
172  {
173  Board2D board;
174  board << domain;
175  board << calculus;
176  board << dual_one_form;
177  board << dual_vector_field;
178  board.saveSVG("usage_dual_one_form.svg");
179  }
180 
181  // test primal flat and sharp
182  const Calculus::DualForm1 flat_sharp_dual_one_form = calculus.flat(dual_vector_field);
183  const Calculus::DualVectorField sharp_flat_dual_vector_field = -calculus.sharp(flat_sharp_dual_one_form);
184 
185  {
186  Board2D board;
187  board << domain;
188  board << calculus;
189  board << flat_sharp_dual_one_form;
190  board << -sharp_flat_dual_vector_field;
191  board.saveSVG("usage_dual_one_form_sharp_flat.svg");
192  }
193 
194  // create primal gradient vector field and hodge*d primal one form
195  const Calculus::DualHodge1 dual_one_hodge = calculus.hodge<1, DUAL>();
196  const Calculus::PrimalForm1 primal_one_form = dual_one_hodge * dual_zero_derivative * dual_zero_form;
197  const Calculus::PrimalVectorField primal_vector_field = calculus.sharp(primal_one_form);
198 
199  {
200  Board2D board;
201  board << domain;
202  board << calculus;
203  board << primal_one_form;
204  board << primal_vector_field;
205  board << dual_vector_field;
206  board.saveSVG("usage_dual_one_form_hodge.svg");
207  }
208  }
209 
210  trace.endBlock();
211 }
212 
213 int main()
214 {
215  usage2d();
216 
217  return 0;
218 }
219 
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Trace trace
Definition: Common.h:130
STL namespace.
double endBlock()
Aim: This class provides static members to create DEC structures from various other DGtal structures...
DGtal is the top-level namespace which contains all DGtal functions and types.
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
std::ostream & info()
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70