DGtal 1.3.0
Loading...
Searching...
No Matches
examplePropagation.cpp
1
17
24#include <iomanip>
25#include "DECExamplesCommon.h"
26#include "DGtal/math/linalg/EigenSupport.h"
27#include "DGtal/dec/DiscreteExteriorCalculus.h"
28#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
29#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
30#include "DGtal/io/boards/Board2D.h"
31
32template <typename Calculus>
33typename Calculus::Scalar
34sample_dual_0_form(const typename Calculus::Properties& properties, const typename Calculus::DualForm0& form, const typename Calculus::Point& point)
35{
36 const typename Calculus::Cell cell = form.myCalculus->myKSpace.uSpel(point);
37 const typename Calculus::Properties::const_iterator iter = properties.find(cell);
38 if (iter == properties.end()) return 0;
39
40 const typename Calculus::Index index = iter->second.index;
41 return form.myContainer(index);
42}
43
44template <typename Calculus>
45void
46set_initial_phase_dir_yy(Calculus& calculus, Eigen::VectorXcd& initial_wave)
47{
48 typedef typename Calculus::Cell Cell;
49 typedef typename Calculus::Index Index;
50 typedef typename Calculus::Point Point;
51
53 for (int xx=9; xx<13; xx++)
54 for (int yy=13; yy<17; yy++)
55 {
56 const Point point(xx,yy);
57 const Cell cell = calculus.myKSpace.uSpel(point);
58 const Index index = calculus.getCellIndex(cell);
59 initial_wave(index) = std::exp(std::complex<double>(0,2*M_PI*(yy-14.5)/8));
60 }
62}
63
64template <typename Calculus>
65void
66set_initial_phase_null(Calculus& calculus, Eigen::VectorXcd& initial_wave)
67{
68 typedef typename Calculus::Cell Cell;
69 typedef typename Calculus::Index Index;
70 typedef typename Calculus::Point Point;
71
73 for (int xx=9; xx<13; xx++)
74 for (int yy=13; yy<17; yy++)
75 {
76 const Point point(xx,yy);
77 const Cell cell = calculus.myKSpace.uSpel(point);
78 const Index index = calculus.getCellIndex(cell);
79 initial_wave(index) = 1;
80 }
82}
83
84using namespace DGtal;
85using std::endl;
86
87void propa_2d()
88{
89 trace.beginBlock("2d propagation");
90
93
94 {
95 trace.beginBlock("solving time dependent equation");
96
97 const Calculus::Scalar cc = 8; // px/s
98 trace.info() << "cc = " << cc << endl;
99
100 const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(29,29));
101 const Calculus calculus = CalculusFactory::createFromDigitalSet(generateDiskSet(domain), false);
102
104 const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>() + 1e-8 * calculus.identity<0, DUAL>();
106 trace.info() << "laplace = " << laplace << endl;
107
108 trace.beginBlock("finding eigen pairs");
110 typedef Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> EigenSolverMatrix;
111 const EigenSolverMatrix eigen_solver(laplace.myContainer);
112
113 const Eigen::VectorXd eigen_values = eigen_solver.eigenvalues();
114 const Eigen::MatrixXd eigen_vectors = eigen_solver.eigenvectors();
116 trace.info() << "eigen_values_range = " << eigen_values.minCoeff() << "/" << eigen_values.maxCoeff() << endl;
117 trace.endBlock();
118
120 const Eigen::VectorXd angular_frequencies = cc * eigen_values.array().sqrt();
122
123 Eigen::VectorXcd initial_wave = Eigen::VectorXcd::Zero(calculus.kFormLength(0, DUAL));
124 //set_initial_phase_null(calculus, initial_wave);
125 set_initial_phase_dir_yy(calculus, initial_wave);
126
127 {
128 Board2D board;
129 board << domain;
130 board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
131 board << Calculus::DualForm0(calculus, initial_wave.real());
132 board.saveSVG("propagation_time_wave_initial_coarse.svg");
133 }
134
136 Eigen::VectorXcd initial_projections = eigen_vectors.transpose() * initial_wave;
138
139 // low pass
141 const Calculus::Scalar lambda_cutoff = 4.5;
142 const Calculus::Scalar angular_frequency_cutoff = 2*M_PI * cc / lambda_cutoff;
143 int cutted = 0;
144 for (int kk=0; kk<initial_projections.rows(); kk++)
145 {
146 const Calculus::Scalar angular_frequency = angular_frequencies(kk);
147 if (angular_frequency < angular_frequency_cutoff) continue;
148 initial_projections(kk) = 0;
149 cutted ++;
150 }
152 trace.info() << "cutted = " << cutted << "/" << initial_projections.rows() << endl;
153
154 {
155 const Eigen::VectorXcd wave = eigen_vectors * initial_projections;
156 Board2D board;
157 board << domain;
158 board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
159 board << Calculus::DualForm0(calculus, wave.real());
160 board.saveSVG("propagation_time_wave_initial_smoothed.svg");
161 }
162
163 const int kk_max = 60;
164 trace.progressBar(0,kk_max);
165 for (int kk=0; kk<kk_max; kk++)
166 {
168 const Calculus::Scalar time = kk/20.;
169 const Eigen::VectorXcd current_projections = (angular_frequencies * std::complex<double>(0,time)).array().exp() * initial_projections.array();
170 const Eigen::VectorXcd current_wave = eigen_vectors * current_projections;
172
173 std::stringstream ss;
174 ss << "propagation_time_wave_solution_" << std::setfill('0') << std::setw(3) << kk << ".svg";
175
176 Board2D board;
177 board << domain;
178 board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
179 board << Calculus::DualForm0(calculus, current_wave.real());
180 board.saveSVG(ss.str().c_str());
181
182 trace.progressBar(kk+1,kk_max);
183 }
184 trace.info() << endl;
185
186 trace.endBlock();
187 }
188
189 {
190 trace.beginBlock("forced oscillations");
191
192 const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(50,50));
193 const Calculus calculus = CalculusFactory::createFromDigitalSet(generateDiskSet(domain), false);
194
195 const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
196 trace.info() << "laplace = " << laplace << endl;
197
198 trace.beginBlock("finding eigen pairs");
199 typedef Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> EigenSolverMatrix;
200 const EigenSolverMatrix eigen_solver(laplace.myContainer);
201
202 const Eigen::VectorXd laplace_eigen_values = eigen_solver.eigenvalues();
203 const Eigen::MatrixXd laplace_eigen_vectors = eigen_solver.eigenvectors();
204 trace.info() << "eigen_values_range = " << laplace_eigen_values.minCoeff() << "/" << laplace_eigen_values.maxCoeff() << endl;
205 trace.endBlock();
206
207 Calculus::DualForm0 concentration(calculus);
208 {
209 const Z2i::Point point(25,25);
210 const Calculus::Cell cell = calculus.myKSpace.uSpel(point);
211 const Calculus::Index index = calculus.getCellIndex(cell);
212 concentration.myContainer(index) = 1;
213 }
214
215 {
216 Board2D board;
217 board << domain;
218 board << concentration;
219 board.saveSVG("propagation_forced_concentration.svg");
220 }
221
222 for (int ll=0; ll<6; ll++)
223 {
225 const Calculus::Scalar lambda = 4*20.75/(1+2*ll);
227 trace.info() << "lambda = " << lambda << endl;
228
230 const Eigen::VectorXd dalembert_eigen_values = (laplace_eigen_values.array() - (2*M_PI/lambda)*(2*M_PI/lambda)).array().inverse();
231 const Eigen::MatrixXd concentration_to_wave = laplace_eigen_vectors * dalembert_eigen_values.asDiagonal() * laplace_eigen_vectors.transpose();
233
235 Calculus::DualForm0 wave(calculus, concentration_to_wave * concentration.myContainer);
237 wave.myContainer /= wave.myContainer(calculus.getCellIndex(calculus.myKSpace.uSpel(Z2i::Point(25,25))));
238
239 {
240 trace.info() << "saving samples" << endl;
241 const Calculus::Properties properties = calculus.getProperties();
242 {
243 std::stringstream ss;
244 ss << "propagation_forced_samples_hv_" << ll << ".dat";
245 std::ofstream handle(ss.str().c_str());
246 for (int kk=0; kk<=50; kk++)
247 {
248 const Z2i::Point point_horizontal(kk,25);
249 const Z2i::Point point_vertical(25,kk);
250 const Calculus::Scalar xx = 2 * (kk/50. - .5);
251 handle << xx << " ";
252 handle << sample_dual_0_form<Calculus>(properties, wave, point_horizontal) << " ";
253 handle << sample_dual_0_form<Calculus>(properties, wave, point_vertical) << endl;
254 }
255 }
256
257 {
258 std::stringstream ss;
259 ss << "propagation_forced_samples_diag_" << ll << ".dat";
260 std::ofstream handle(ss.str().c_str());
261 for (int kk=0; kk<=50; kk++)
262 {
263 const Z2i::Point point_diag_pos(kk,kk);
264 const Z2i::Point point_diag_neg(kk,50-kk);
265 const Calculus::Scalar xx = 2 * sqrt(2) * (kk/50. - .5);
266 handle << xx << " ";
267 handle << sample_dual_0_form<Calculus>(properties, wave, point_diag_pos) << " ";
268 handle << sample_dual_0_form<Calculus>(properties, wave, point_diag_neg) << endl;
269 }
270 }
271 }
272
273 {
274 std::stringstream ss;
275 ss << "propagation_forced_wave_" << ll << ".svg";
276 Board2D board;
277 board << domain;
278 board << CustomStyle("KForm", new KFormStyle2D(-.5, 1));
279 board << wave;
280 board.saveSVG(ss.str().c_str());
281 }
282 }
283
284 trace.endBlock();
285 }
286 trace.endBlock();
287}
288
289int main()
290{
291 propa_2d();
292 return 0;
293}
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...
void beginBlock(const std::string &keyword="")
std::ostream & info()
void progressBar(const double currentValue, const double maximalValue)
double endBlock()
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
SMesh::Index Index
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
@ DUAL
Definition: Duality.h:62
int main(int argc, char **argv)
MyPointD Point
Definition: testClone2.cpp:383
KSpace::Cell Cell
Domain domain