DGtal  1.1.0
examplePropagation.cpp
1 
16 
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 
32 template <typename Calculus>
33 typename Calculus::Scalar
34 sample_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 
44 template <typename Calculus>
45 void
46 set_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 
64 template <typename Calculus>
65 void
66 set_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 
84 using namespace DGtal;
85 using std::endl;
86 
87 void 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 
289 int main()
290 {
291  propa_2d();
292  return 0;
293 }
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition: DiscreteExteriorCalculus.h:98
DGtal::HyperRectDomain< Space >
DGtal::Trace::endBlock
double endBlock()
DGtal::trace
Trace trace
Definition: Common.h:150
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::CustomStyle
Definition: Board2D.h:217
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::progressBar
void progressBar(const double currentValue, const double maximalValue)
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition: ClosedIntegerHalfPlane.h:49
DGtal::DiscreteExteriorCalculusFactory
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Definition: DiscreteExteriorCalculusFactory.h:62
LibBoard::Board::saveSVG
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::Board2D
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::PointVector< dim, Integer >
DGtal::DUAL
@ DUAL
Definition: Duality.h:62
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::KhalimskyCell< dim, Integer >