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