DGtal  0.9.2
examplePropagation.cpp
1 #include <iomanip>
3 
4 #include "DECExamplesCommon.h"
5 
6 #include "DGtal/math/linalg/EigenSupport.h"
7 #include "DGtal/dec/DiscreteExteriorCalculus.h"
8 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
9 #include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
10 #include "DGtal/io/boards/Board2D.h"
11 
12 template <typename Calculus>
13 typename Calculus::Scalar
14 sample_dual_0_form(const typename Calculus::Properties& properties, const typename Calculus::DualForm0& form, const typename Calculus::Point& point)
15 {
16  const typename Calculus::Cell cell = form.myCalculus->myKSpace.uSpel(point);
17  const typename Calculus::Properties::const_iterator iter = properties.find(cell);
18  if (iter == properties.end()) return 0;
19 
20  const typename Calculus::Index index = iter->second.index;
21  return form.myContainer(index);
22 }
23 
24 template <typename Calculus>
25 void
26 set_initial_phase_dir_yy(Calculus& calculus, Eigen::VectorXcd& initial_wave)
27 {
28  typedef typename Calculus::Cell Cell;
29  typedef typename Calculus::Index Index;
30  typedef typename Calculus::Point Point;
31 
33  for (int xx=9; xx<13; xx++)
34  for (int yy=13; yy<17; yy++)
35  {
36  const Point point(xx,yy);
37  const Cell cell = calculus.myKSpace.uSpel(point);
38  const Index index = calculus.getCellIndex(cell);
39  initial_wave(index) = std::exp(std::complex<double>(0,2*M_PI*(yy-14.5)/8));
40  }
42 }
43 
44 template <typename Calculus>
45 void
46 set_initial_phase_null(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) = 1;
60  }
62 }
63 
64 using namespace DGtal;
65 using std::endl;
66 
67 void propa_2d()
68 {
69  trace.beginBlock("2d propagation");
70 
73 
74  {
75  trace.beginBlock("solving time dependent equation");
76 
77  const Calculus::Scalar cc = 8; // px/s
78  trace.info() << "cc = " << cc << endl;
79 
80  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(29,29));
81  const Calculus calculus = CalculusFactory::createFromDigitalSet(generateDiskSet(domain), false);
82 
84  const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>() + 1e-8 * calculus.identity<0, DUAL>();
86  trace.info() << "laplace = " << laplace << endl;
87 
88  trace.beginBlock("finding eigen pairs");
90  typedef Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> EigenSolverMatrix;
91  const EigenSolverMatrix eigen_solver(laplace.myContainer);
92 
93  const Eigen::VectorXd eigen_values = eigen_solver.eigenvalues();
94  const Eigen::MatrixXd eigen_vectors = eigen_solver.eigenvectors();
96  trace.info() << "eigen_values_range = " << eigen_values.minCoeff() << "/" << eigen_values.maxCoeff() << endl;
97  trace.endBlock();
98 
100  const Eigen::VectorXd angular_frequencies = cc * eigen_values.array().sqrt();
102 
103  Eigen::VectorXcd initial_wave = Eigen::VectorXcd::Zero(calculus.kFormLength(0, DUAL));
104  //set_initial_phase_null(calculus, initial_wave);
105  set_initial_phase_dir_yy(calculus, initial_wave);
106 
107  {
108  Board2D board;
109  board << domain;
110  board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
111  board << Calculus::DualForm0(calculus, initial_wave.real());
112  board.saveSVG("propagation_time_wave_initial_coarse.svg");
113  }
114 
116  Eigen::VectorXcd initial_projections = eigen_vectors.transpose() * initial_wave;
118 
119  // low pass
121  const Calculus::Scalar lambda_cutoff = 4.5;
122  const Calculus::Scalar angular_frequency_cutoff = 2*M_PI * cc / lambda_cutoff;
123  int cutted = 0;
124  for (int kk=0; kk<initial_projections.rows(); kk++)
125  {
126  const Calculus::Scalar angular_frequency = angular_frequencies(kk);
127  if (angular_frequency < angular_frequency_cutoff) continue;
128  initial_projections(kk) = 0;
129  cutted ++;
130  }
132  trace.info() << "cutted = " << cutted << "/" << initial_projections.rows() << endl;
133 
134  {
135  const Eigen::VectorXcd wave = eigen_vectors * initial_projections;
136  Board2D board;
137  board << domain;
138  board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
139  board << Calculus::DualForm0(calculus, wave.real());
140  board.saveSVG("propagation_time_wave_initial_smoothed.svg");
141  }
142 
143  const int kk_max = 60;
144  trace.progressBar(0,kk_max);
145  for (int kk=0; kk<kk_max; kk++)
146  {
148  const Calculus::Scalar time = kk/20.;
149  const Eigen::VectorXcd current_projections = (angular_frequencies * std::complex<double>(0,time)).array().exp() * initial_projections.array();
150  const Eigen::VectorXcd current_wave = eigen_vectors * current_projections;
152 
153  std::stringstream ss;
154  ss << "propagation_time_wave_solution_" << std::setfill('0') << std::setw(3) << kk << ".svg";
155 
156  Board2D board;
157  board << domain;
158  board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
159  board << Calculus::DualForm0(calculus, current_wave.real());
160  board.saveSVG(ss.str().c_str());
161 
162  trace.progressBar(kk+1,kk_max);
163  }
164  trace.info() << endl;
165 
166  trace.endBlock();
167  }
168 
169  {
170  trace.beginBlock("forced oscillations");
171 
172  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(50,50));
173  const Calculus calculus = CalculusFactory::createFromDigitalSet(generateDiskSet(domain), false);
174 
175  const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
176  trace.info() << "laplace = " << laplace << endl;
177 
178  trace.beginBlock("finding eigen pairs");
179  typedef Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> EigenSolverMatrix;
180  const EigenSolverMatrix eigen_solver(laplace.myContainer);
181 
182  const Eigen::VectorXd laplace_eigen_values = eigen_solver.eigenvalues();
183  const Eigen::MatrixXd laplace_eigen_vectors = eigen_solver.eigenvectors();
184  trace.info() << "eigen_values_range = " << laplace_eigen_values.minCoeff() << "/" << laplace_eigen_values.maxCoeff() << endl;
185  trace.endBlock();
186 
187  Calculus::DualForm0 concentration(calculus);
188  {
189  const Z2i::Point point(25,25);
190  const Calculus::Cell cell = calculus.myKSpace.uSpel(point);
191  const Calculus::Index index = calculus.getCellIndex(cell);
192  concentration.myContainer(index) = 1;
193  }
194 
195  {
196  Board2D board;
197  board << domain;
198  board << concentration;
199  board.saveSVG("propagation_forced_concentration.svg");
200  }
201 
202  for (int ll=0; ll<6; ll++)
203  {
205  const Calculus::Scalar lambda = 4*20.75/(1+2*ll);
207  trace.info() << "lambda = " << lambda << endl;
208 
210  const Eigen::VectorXd dalembert_eigen_values = (laplace_eigen_values.array() - (2*M_PI/lambda)*(2*M_PI/lambda)).array().inverse();
211  const Eigen::MatrixXd concentration_to_wave = laplace_eigen_vectors * dalembert_eigen_values.asDiagonal() * laplace_eigen_vectors.transpose();
213 
215  Calculus::DualForm0 wave(calculus, concentration_to_wave * concentration.myContainer);
217  wave.myContainer /= wave.myContainer(calculus.getCellIndex(calculus.myKSpace.uSpel(Z2i::Point(25,25))));
218 
219  {
220  trace.info() << "saving samples" << endl;
221  const Calculus::Properties properties = calculus.getProperties();
222  {
223  std::stringstream ss;
224  ss << "propagation_forced_samples_hv_" << ll << ".dat";
225  std::ofstream handle(ss.str().c_str());
226  for (int kk=0; kk<=50; kk++)
227  {
228  const Z2i::Point point_horizontal(kk,25);
229  const Z2i::Point point_vertical(25,kk);
230  const Calculus::Scalar xx = 2 * (kk/50. - .5);
231  handle << xx << " ";
232  handle << sample_dual_0_form<Calculus>(properties, wave, point_horizontal) << " ";
233  handle << sample_dual_0_form<Calculus>(properties, wave, point_vertical) << endl;
234  }
235  }
236 
237  {
238  std::stringstream ss;
239  ss << "propagation_forced_samples_diag_" << ll << ".dat";
240  std::ofstream handle(ss.str().c_str());
241  for (int kk=0; kk<=50; kk++)
242  {
243  const Z2i::Point point_diag_pos(kk,kk);
244  const Z2i::Point point_diag_neg(kk,50-kk);
245  const Calculus::Scalar xx = 2 * sqrt(2) * (kk/50. - .5);
246  handle << xx << " ";
247  handle << sample_dual_0_form<Calculus>(properties, wave, point_diag_pos) << " ";
248  handle << sample_dual_0_form<Calculus>(properties, wave, point_diag_neg) << endl;
249  }
250  }
251  }
252 
253  {
254  std::stringstream ss;
255  ss << "propagation_forced_wave_" << ll << ".svg";
256  Board2D board;
257  board << domain;
258  board << CustomStyle("KForm", new KFormStyle2D(-.5, 1));
259  board << wave;
260  board.saveSVG(ss.str().c_str());
261  }
262  }
263 
264  trace.endBlock();
265  }
266  trace.endBlock();
267 }
268 
269 int main(int argc, char* argv[])
270 {
271  propa_2d();
272  return 0;
273 }
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
void progressBar(const double currentValue, const double maximalValue)
Trace trace
Definition: Common.h:130
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