DGtal 1.3.0
Searching...
No Matches
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"
37
38using namespace std;
39using namespace DGtal;
40
41void 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
234int 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:1012
static const L2Metric l2Metric
Definition: StdDefs.h:123
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:154
@ PRIMAL
Definition: Duality.h:61
@ DUAL
Definition: Duality.h:62
STL namespace.
int main(int argc, char **argv)
Domain domain