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