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