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