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