DGtal  0.9.2
exampleDiscreteExteriorCalculusChladni.cpp
1 #include <sstream>
3 #include <string>
4 
5 // always include EigenSupport.h before any other Eigen headers
6 #include "DGtal/math/linalg/EigenSupport.h"
7 #include <Eigen/Eigenvalues>
8 
9 #include "DGtal/dec/DiscreteExteriorCalculus.h"
10 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
11 
12 #include "DGtal/base/Common.h"
13 #include "DGtal/helpers/StdDefs.h"
14 #include "DGtal/io/boards/Board2D.h"
15 #include "DGtal/math/linalg/EigenSupport.h"
16 #include "DGtal/dec/DiscreteExteriorCalculus.h"
17 #include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
18 
19 using namespace std;
20 using namespace DGtal;
21 using namespace Eigen;
22 
23 double standard_deviation(const VectorXd& xx)
24 {
25  const double mean = xx.mean();
26  double accum = 0;
27  for (int kk=0, kk_max=xx.rows(); kk<kk_max; kk++)
28  accum += (xx(kk)-mean)*(xx(kk)-mean);
29  return sqrt(accum)/(xx.size()-1);
30 }
31 
32 int main()
33 {
34  trace.beginBlock("building calculus");
35 
36  const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(10,10));
37 
39  Calculus calculus;
40  calculus.initKSpace<Z2i::Domain>(domain);
41 
42  // bottom linear structure
43  // left and right Dirichlet boundary condition
44  for (int kk=2; kk<=20; kk++)
45  calculus.insertSCell( calculus.myKSpace.sCell(Z2i::Point(kk, 1)) );
46 
47  // top linear structure
48  // left Neumann boundary condition, right Dirichlet boundary condition
49  for (int kk=3; kk<=20; kk++)
50  calculus.insertSCell( calculus.myKSpace.sCell(Z2i::Point(kk, 21)) );
51 
52  for (int kk=3; kk<20; kk++)
53  for (int ll=3; ll<20; ll++)
54  {
55  if (kk==11 && ll==11) continue;
56  calculus.insertSCell( calculus.myKSpace.sCell(Z2i::Point(kk, ll)) );
57  }
58 
59  calculus.updateIndexes();
60  trace.info() << calculus << endl;
61 
62  trace.endBlock();
63 
64  trace.beginBlock("building laplace");
65 
66  const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>();
67  trace.info() << "laplace=" << laplace << endl;
68 
69  {
70  Board2D board;
71  board << domain;
72  board << calculus;
73  board.saveSVG("chladni_calculus.svg");
74  }
75 
76  trace.endBlock();
77 
78  trace.beginBlock("finding laplace eigen vectors");
79 
80  typedef SelfAdjointEigenSolver<MatrixXd> EigenSolverMatrix;
81  const EigenSolverMatrix eigen_solver(laplace.myContainer);
82 
83  const VectorXd eigen_values = eigen_solver.eigenvalues();
84  const MatrixXd eigen_vectors = eigen_solver.eigenvectors();
85  for (int kk=0; kk<laplace.myContainer.rows(); kk++)
86  {
87  Calculus::Scalar eigen_value = eigen_values(kk, 0);
88 
89  const VectorXd eigen_vector = eigen_vectors.col(kk);
90  const Calculus::DualForm0 eigen_form = Calculus::DualForm0(calculus, eigen_vector);
91  std::stringstream ss;
92  ss << "chladni_eigen_" << kk << ".svg";
93  const std::string filename = ss.str();
94  ss << "chladni_eigen_vector_" << kk << ".svg";
95  trace.info() << kk << " " << eigen_value << " " << sqrt(eigen_value) << " " << eigen_vector.minCoeff() << " " << eigen_vector.maxCoeff() << " " << standard_deviation(eigen_vector) << endl;
96 
97  Board2D board;
98  board << domain;
99  board << calculus;
100  board << CustomStyle("KForm", new KFormStyle2D(eigen_vectors.minCoeff(),eigen_vectors.maxCoeff()));
101  board << eigen_form;
102  board.saveSVG(filename.c_str());
103  }
104 
105  trace.endBlock();
106 
107  return 0;
108 }
109 
void beginBlock(const std::string &keyword="")
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Trace trace
Definition: Common.h:130
STL namespace.
double endBlock()
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()
void initKSpace(ConstAlias< TDomain > domain)
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70