DGtal 1.3.0
Loading...
Searching...
No Matches
testEmbedding.cpp
1#include "DGtal/math/linalg/EigenSupport.h"
2#include "DGtal/dec/DiscreteExteriorCalculus.h"
3#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
4#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
5
6#define NOVIEWER // comment to enable 3d viewers
7
8#if !defined(NOVIEWER)
9#include "DGtal/io/viewers/Viewer3D.h"
10#endif
11#include "DGtal/io/boards/Board2D.h"
12#include "DGtal/base/Common.h"
13#include "DGtal/helpers/StdDefs.h"
14
15using namespace DGtal;
16using namespace std;
17using namespace Eigen;
18
19template <typename OperatorAA, typename OperatorBB>
20bool
21equal(const OperatorAA& aa, const OperatorBB& bb)
22{
23 return MatrixXd(aa.myContainer) == MatrixXd(bb.myContainer);
24}
25
26int main(int , char** )
27{
28#if !defined(NOVIEWER)
30#endif
34
35#if !defined(NOVIEWER)
36 QApplication app(argc, argv);
37 Z3i::KSpace kspace_3d;
38#endif
39 Z2i::KSpace kspace_2d;
40
41#if !defined(NOVIEWER)
42 Viewer viewer1(kspace_3d);
43 viewer1.show();
44 viewer1.setWindowTitle("embedding_1d_calculus_3d");
45 viewer1.camera()->setPosition( Vec(2,2,2) );
46 viewer1.camera()->setUpVector( Vec(0,0,1), false );
47 viewer1.camera()->lookAt( Vec(0,0,0) );
48
49 Viewer viewer2(kspace_3d);
50 viewer2.show();
51 viewer2.setWindowTitle("embedding_2d_calculus_3d");
52 viewer2.camera()->setPosition( Vec(2,2,2) );
53 viewer2.camera()->setUpVector( Vec(0,0,2), false );
54 viewer2.camera()->lookAt( Vec(0,0,0) );
55
56 Viewer viewer3(kspace_3d);
57 viewer3.show();
58 viewer3.setWindowTitle("embedding_2d_calculus_3d_no_border");
59 viewer3.camera()->setPosition( Vec(2,2,2) );
60 viewer3.camera()->setUpVector( Vec(0,0,2), false );
61 viewer3.camera()->lookAt( Vec(0,0,0) );
62#endif
63
64 {
65 trace.beginBlock("1d manifold embedding");
66
72
74 typedef std::set<Calculus1D::SCell> SCells1D;
75 SCells1D ncells_1d_factory;
77 SCells1D cells_1d_manual;
78 Calculus1D calculus_1d_manual;
79 for (int kk=0; kk<31; kk++)
80 {
81 Calculus1D::KSpace::Point point;
82 point[0] = kk;
83 const Calculus1D::SCell cell = calculus_1d_manual.myKSpace.sCell(point);
84 calculus_1d_manual.insertSCell(cell, 1, kk == 0 || kk == 30 ? 1/2. : 1);
85 cells_1d_manual.insert(cell);
86 if (kk%2 != 0) ncells_1d_factory.insert(cell);
87 }
88 calculus_1d_manual.updateIndexes();
89 trace.info() << "calculus_1d_manual=" << calculus_1d_manual << endl;
90
91 {
92 Board2D board;
93 board << Z2i::Domain(Z2i::Point(-1,-1), Z2i::Point(15,0));
94 board.setFillColor( DGtal::Color(128,128,128) );
95 for (Calculus1D::ConstIterator ii=calculus_1d_manual.begin(), ie=calculus_1d_manual.end(); ii!=ie; ii++)
96 {
97 const Z2i::Point point(calculus_1d_manual.myKSpace.uKCoord(ii->first, 0), 0);
98 const Z2i::KSpace::SCell cell = kspace_2d.sCell(point);
100 }
101 board.saveSVG("embedding_1d_calculus_1d.svg");
102 }
103
105 const Calculus1D calculus_1d_factory = CalculusFactory::createFromNSCells<1>(ncells_1d_factory.begin(), ncells_1d_factory.end(), true);
107 trace.info() << "calculus_1d_factory=" << calculus_1d_factory << endl;
108
109 Calculus2D calculus_2d_manual;
110 {
111 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,0)), 1, 1/2. );
112 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,1), Calculus2D::KSpace::POS) );
113 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,2)) );
114 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,2), Calculus2D::KSpace::POS) );
115 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,2)) );
116 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,1), Calculus2D::KSpace::NEG) );
117 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,0)) );
118 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,-1), Calculus2D::KSpace::NEG) );
119 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,-2)) );
120 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,-2), Calculus2D::KSpace::NEG) );
121 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,-2)) );
122 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(5,-2), Calculus2D::KSpace::NEG) );
123 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(4,-2)) );
124 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(3,-2), Calculus2D::KSpace::NEG) );
125 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,-2)) );
126 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,-2), Calculus2D::KSpace::NEG) );
127 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,-2)) );
128 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,-2), Calculus2D::KSpace::NEG) );
129 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,-2)) );
130 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,-1), Calculus2D::KSpace::POS) );
131 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,0)) );
132 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,1), Calculus2D::KSpace::POS) );
133 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,2)) );
134 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,2), Calculus2D::KSpace::POS) );
135 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,2)) );
136 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,2), Calculus2D::KSpace::POS) );
137 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,2)) );
138 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,1), Calculus2D::KSpace::NEG) );
139 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,0)) );
140 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,0), Calculus2D::KSpace::NEG) );
141 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,0)), 1, 1/2.);
142 calculus_2d_manual.updateIndexes();
143 }
144 trace.info() << "calculus_2d_manual=" << calculus_2d_manual << endl;
145
146 {
147 Board2D board;
148 board << Z2i::Domain(Z2i::Point(-2,-2), Z2i::Point(4,1));
149 board << calculus_2d_manual;
150 board.saveSVG("embedding_1d_calculus_2d.svg");
151 }
152
154 typedef std::list<Calculus2D::SCell> SCells2D;
155 SCells2D ncells_2d_factory;
157
158 {
159 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,1), Calculus2D::KSpace::POS) );
160 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,2), Calculus2D::KSpace::POS) );
161 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,1), Calculus2D::KSpace::NEG) );
162 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,-1), Calculus2D::KSpace::NEG) );
163 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,-2), Calculus2D::KSpace::NEG) );
164 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(5,-2), Calculus2D::KSpace::NEG) );
165 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(3,-2), Calculus2D::KSpace::NEG) );
166 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,-2), Calculus2D::KSpace::NEG) );
167 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,-2), Calculus2D::KSpace::NEG) );
168 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,-1), Calculus2D::KSpace::POS) );
169 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,1), Calculus2D::KSpace::POS) );
170 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,2), Calculus2D::KSpace::POS) );
171 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,2), Calculus2D::KSpace::POS) );
172 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,1), Calculus2D::KSpace::NEG) );
173 ncells_2d_factory.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,0), Calculus2D::KSpace::NEG) );
174 }
176 const Calculus2D calculus_2d_factory = CalculusFactory::createFromNSCells<1>(ncells_2d_factory.begin(), ncells_2d_factory.end(), true);
178 trace.info() << "calculus_2d_factory=" << calculus_2d_factory << endl;
179
180 SCells2D cells_2d_manual;
181 {
182 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,0)) );
183 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,1), Calculus2D::KSpace::POS) );
184 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,2)) );
185 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,2), Calculus2D::KSpace::POS) );
186 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,2)) );
187 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,1), Calculus2D::KSpace::NEG) );
188 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,0)) );
189 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,-1), Calculus2D::KSpace::NEG) );
190 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(8,-2)) );
191 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(7,-2), Calculus2D::KSpace::NEG) );
192 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(6,-2)) );
193 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(5,-2), Calculus2D::KSpace::NEG) );
194 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(4,-2)) );
195 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(3,-2), Calculus2D::KSpace::NEG) );
196 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,-2)) );
197 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,-2), Calculus2D::KSpace::NEG) );
198 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,-2)) );
199 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,-2), Calculus2D::KSpace::NEG) );
200 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,-2)) );
201 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,-1), Calculus2D::KSpace::POS) );
202 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,0)) );
203 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,1), Calculus2D::KSpace::POS) );
204 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-2,2)) );
205 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(-1,2), Calculus2D::KSpace::POS) );
206 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,2)) );
207 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,2), Calculus2D::KSpace::POS) );
208 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,2)) );
209 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,1), Calculus2D::KSpace::NEG) );
210 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(2,0)) );
211 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(1,0), Calculus2D::KSpace::NEG) );
212 cells_2d_manual.push_back( calculus_2d_manual.myKSpace.sCell(Z2i::Point(0,0)) );
213 }
214
215 Calculus3D calculus_3d_manual;
216 {
217 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,0,0)), 1, 1/2. );
218 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,0,0), Calculus3D::KSpace::POS) );
219 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,0,0)) );
220 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,0,0), Calculus3D::KSpace::POS) );
221 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,0,0)) );
222 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,1,0), Calculus3D::KSpace::POS) );
223 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,2,0)) );
224 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,3,0), Calculus3D::KSpace::POS) );
225 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,4,0)) );
226 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,4,0), Calculus3D::KSpace::NEG) );
227 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,0)) );
228 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,4,0), Calculus3D::KSpace::NEG) );
229 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,4,0)) );
230 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,3,0), Calculus3D::KSpace::NEG) );
231 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,2,0)) );
232 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,2,0), Calculus3D::KSpace::POS) );
233 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,0)) );
234 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,1), Calculus3D::KSpace::POS) );
235 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,2)) );
236 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,2), Calculus3D::KSpace::POS) );
237 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,2)) );
238 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,2), Calculus3D::KSpace::POS) );
239 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,2)) );
240 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,1), Calculus3D::KSpace::NEG) );
241 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,0)) );
242 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,-1), Calculus3D::KSpace::NEG) );
243 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,-2)) );
244 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,-2), Calculus3D::KSpace::NEG) );
245 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,-2)) );
246 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,-2), Calculus3D::KSpace::NEG) );
247 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,-2)), 1, 1/2. );
248 calculus_3d_manual.updateIndexes();
249 }
250 trace.info() << "calculus_3d_manual=" << calculus_3d_manual << endl;
251
252#if !defined(NOVIEWER)
254 viewer1 << Viewer::updateDisplay;
255#endif
256
258 typedef std::vector<Calculus3D::SCell> SCells3D;
259 SCells3D ncells_3d_factory;
261 {
262 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,0,0), Calculus3D::KSpace::POS) );
263 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,0,0), Calculus3D::KSpace::POS) );
264 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,1,0), Calculus3D::KSpace::POS) );
265 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,3,0), Calculus3D::KSpace::POS) );
266 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,4,0), Calculus3D::KSpace::NEG) );
267 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,4,0), Calculus3D::KSpace::NEG) );
268 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,3,0), Calculus3D::KSpace::NEG) );
269 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,2,0), Calculus3D::KSpace::POS) );
270 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,1), Calculus3D::KSpace::POS) );
271 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,2), Calculus3D::KSpace::POS) );
272 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,2), Calculus3D::KSpace::POS) );
273 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,1), Calculus3D::KSpace::NEG) );
274 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,-1), Calculus3D::KSpace::NEG) );
275 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,-2), Calculus3D::KSpace::NEG) );
276 ncells_3d_factory.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,-2), Calculus3D::KSpace::NEG) );
277 }
279 const Calculus3D calculus_3d_factory = CalculusFactory::createFromNSCells<1>(ncells_3d_factory.begin(), ncells_3d_factory.end(), true);
281 trace.info() << "calculus_3d_factory=" << calculus_3d_factory << endl;
282
283 SCells3D cells_3d_manual;
284 {
285 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,0,0)) );
286 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,0,0), Calculus3D::KSpace::POS) );
287 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,0,0)) );
288 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,0,0), Calculus3D::KSpace::POS) );
289 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,0,0)) );
290 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,1,0), Calculus3D::KSpace::POS) );
291 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,2,0)) );
292 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,3,0), Calculus3D::KSpace::POS) );
293 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,4,0)) );
294 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(3,4,0), Calculus3D::KSpace::NEG) );
295 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,0)) );
296 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,4,0), Calculus3D::KSpace::NEG) );
297 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,4,0)) );
298 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,3,0), Calculus3D::KSpace::NEG) );
299 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(0,2,0)) );
300 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(1,2,0), Calculus3D::KSpace::POS) );
301 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,0)) );
302 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,1), Calculus3D::KSpace::POS) );
303 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,2)) );
304 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,2), Calculus3D::KSpace::POS) );
305 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,2)) );
306 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,2), Calculus3D::KSpace::POS) );
307 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,2)) );
308 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,1), Calculus3D::KSpace::NEG) );
309 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,0)) );
310 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,-1), Calculus3D::KSpace::NEG) );
311 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,6,-2)) );
312 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,5,-2), Calculus3D::KSpace::NEG) );
313 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,4,-2)) );
314 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,3,-2), Calculus3D::KSpace::NEG) );
315 cells_3d_manual.push_back( calculus_3d_manual.myKSpace.sCell(Z3i::Point(2,2,-2)) );
316 }
317 trace.beginBlock("checking operators");
318
319 const Calculus1D::PrimalIdentity0 reorder_0cell_1d = calculus_1d_manual.reorder<0, PRIMAL>(cells_1d_manual.begin(), cells_1d_manual.end());
320 const Calculus2D::PrimalIdentity0 reorder_0cell_2d = calculus_2d_manual.reorder<0, PRIMAL>(cells_2d_manual.begin(), cells_2d_manual.end());
321 const Calculus3D::PrimalIdentity0 reorder_0cell_3d = calculus_3d_manual.reorder<0, PRIMAL>(cells_3d_manual.begin(), cells_3d_manual.end());
322 const Calculus1D::PrimalIdentity1 reorder_1cell_1d = calculus_1d_manual.reorder<1, PRIMAL>(cells_1d_manual.begin(), cells_1d_manual.end());
323 const Calculus2D::PrimalIdentity1 reorder_1cell_2d = calculus_2d_manual.reorder<1, PRIMAL>(cells_2d_manual.begin(), cells_2d_manual.end());
324 const Calculus3D::PrimalIdentity1 reorder_1cell_3d = calculus_3d_manual.reorder<1, PRIMAL>(cells_3d_manual.begin(), cells_3d_manual.end());
325 const Calculus1D::DualIdentity0 reorder_0cellp_1d = calculus_1d_manual.reorder<0, DUAL>(cells_1d_manual.begin(), cells_1d_manual.end());
326 const Calculus2D::DualIdentity0 reorder_0cellp_2d = calculus_2d_manual.reorder<0, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
327 const Calculus3D::DualIdentity0 reorder_0cellp_3d = calculus_3d_manual.reorder<0, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
328 const Calculus1D::DualIdentity1 reorder_1cellp_1d = calculus_1d_manual.reorder<1, DUAL>(cells_1d_manual.begin(), cells_1d_manual.end());
329 const Calculus2D::DualIdentity1 reorder_1cellp_2d = calculus_2d_manual.reorder<1, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
330 const Calculus3D::DualIdentity1 reorder_1cellp_3d = calculus_3d_manual.reorder<1, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
331
332 { // testing primal operators
333 const Calculus1D::PrimalIdentity0 primal_laplace_1d = reorder_0cell_1d * calculus_1d_manual.laplace<PRIMAL>() * reorder_0cell_1d.transpose();
334 const Calculus2D::PrimalIdentity0 primal_laplace_2d = reorder_0cell_2d * calculus_2d_manual.laplace<PRIMAL>() * reorder_0cell_2d.transpose();
335 const Calculus3D::PrimalIdentity0 primal_laplace_3d = reorder_0cell_3d * calculus_3d_manual.laplace<PRIMAL>() * reorder_0cell_3d.transpose();
336 trace.info() << "primal_laplace_1d=" << primal_laplace_1d << endl;
337 trace.info() << "primal_laplace_2d=" << primal_laplace_2d << endl;
338 trace.info() << "primal_laplace_3d=" << primal_laplace_3d << endl;
339 trace.info() << "primal_laplace_container=" << endl << MatrixXd(primal_laplace_1d.myContainer) << endl;
340 FATAL_ERROR( equal(
341 reorder_1cellp_1d * calculus_1d_manual.hodge<0, PRIMAL>() * reorder_0cell_1d.transpose(),
342 reorder_1cellp_2d * calculus_2d_manual.hodge<0, PRIMAL>() * reorder_0cell_2d.transpose()
343 ) );
344 FATAL_ERROR( equal(
345 reorder_1cellp_1d * calculus_1d_manual.hodge<0, PRIMAL>() * reorder_0cell_1d.transpose(),
346 reorder_1cellp_3d * calculus_3d_manual.hodge<0, PRIMAL>() * reorder_0cell_3d.transpose()
347 ) );
348 FATAL_ERROR( equal(
349 reorder_0cellp_1d * calculus_1d_manual.hodge<1, PRIMAL>() * reorder_1cell_1d.transpose(),
350 reorder_0cellp_2d * calculus_2d_manual.hodge<1, PRIMAL>() * reorder_1cell_2d.transpose()
351 ) );
352 FATAL_ERROR( equal(
353 reorder_0cellp_1d * calculus_1d_manual.hodge<1, PRIMAL>() * reorder_1cell_1d.transpose(),
354 reorder_0cellp_3d * calculus_3d_manual.hodge<1, PRIMAL>() * reorder_1cell_3d.transpose()
355 ) );
356 FATAL_ERROR( equal(
357 reorder_1cell_1d * calculus_1d_manual.derivative<0, PRIMAL>() * reorder_0cell_1d.transpose(),
358 reorder_1cell_2d * calculus_2d_manual.derivative<0, PRIMAL>() * reorder_0cell_2d.transpose()
359 ) );
360 FATAL_ERROR( equal(
361 reorder_1cell_1d * calculus_1d_manual.derivative<0, PRIMAL>() * reorder_0cell_1d.transpose(),
362 reorder_1cell_3d * calculus_3d_manual.derivative<0, PRIMAL>() * reorder_0cell_3d.transpose()
363 ) );
364 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_2d) );
365 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_3d) );
366 }
367
368 { // testing dual operators
369 const Calculus1D::DualIdentity0 dual_laplace_1d = reorder_0cellp_1d * calculus_1d_manual.laplace<DUAL>() * reorder_0cellp_1d.transpose();
370 const Calculus2D::DualIdentity0 dual_laplace_2d = reorder_0cellp_2d * calculus_2d_manual.laplace<DUAL>() * reorder_0cellp_2d.transpose();
371 const Calculus3D::DualIdentity0 dual_laplace_3d = reorder_0cellp_3d * calculus_3d_manual.laplace<DUAL>() * reorder_0cellp_3d.transpose();
372 trace.info() << "dual_laplace_1d=" << dual_laplace_1d << endl;
373 trace.info() << "dual_laplace_2d=" << dual_laplace_2d << endl;
374 trace.info() << "dual_laplace_3d=" << dual_laplace_3d << endl;
375 trace.info() << "dual_laplace_container=" << endl << MatrixXd(dual_laplace_1d.myContainer) << endl;
376 FATAL_ERROR( equal(
377 reorder_1cell_1d * calculus_1d_manual.hodge<0, DUAL>() * reorder_0cellp_1d.transpose(),
378 reorder_1cell_2d * calculus_2d_manual.hodge<0, DUAL>() * reorder_0cellp_2d.transpose()
379 ) );
380 FATAL_ERROR( equal(
381 reorder_1cell_1d * calculus_1d_manual.hodge<0, DUAL>() * reorder_0cellp_1d.transpose(),
382 reorder_1cell_3d * calculus_3d_manual.hodge<0, DUAL>() * reorder_0cellp_3d.transpose()
383 ) );
384 FATAL_ERROR( equal(
385 reorder_0cell_1d * calculus_1d_manual.hodge<1, DUAL>() * reorder_1cellp_1d.transpose(),
386 reorder_0cell_2d * calculus_2d_manual.hodge<1, DUAL>() * reorder_1cellp_2d.transpose()
387 ) );
388 FATAL_ERROR( equal(
389 reorder_0cell_1d * calculus_1d_manual.hodge<1, DUAL>() * reorder_1cellp_1d.transpose(),
390 reorder_0cell_3d * calculus_3d_manual.hodge<1, DUAL>() * reorder_1cellp_3d.transpose()
391 ) );
392 FATAL_ERROR( equal(
393 reorder_1cellp_1d * calculus_1d_manual.derivative<0, DUAL>() * reorder_0cellp_1d.transpose(),
394 reorder_1cellp_2d * calculus_2d_manual.derivative<0, DUAL>() * reorder_0cellp_2d.transpose()
395 ) );
396 FATAL_ERROR( equal(
397 reorder_1cellp_1d * calculus_1d_manual.derivative<0, DUAL>() * reorder_0cellp_1d.transpose(),
398 reorder_1cellp_3d * calculus_3d_manual.derivative<0, DUAL>() * reorder_0cellp_3d.transpose()
399 ) );
400 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_2d) );
401 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_3d) );
402 }
403
404 { // checking dual laplace factory calculus vs manual calculus
405 const Calculus1D::DualIdentity0 reorder_0cellp_1d_factory = calculus_1d_factory.reorder<0, DUAL>(cells_1d_manual.begin(), cells_1d_manual.end());
406 const Calculus2D::DualIdentity0 reorder_0cellp_2d_factory = calculus_2d_factory.reorder<0, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
407 const Calculus3D::DualIdentity0 reorder_0cellp_3d_factory = calculus_3d_factory.reorder<0, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
408 FATAL_ERROR( equal(
409 reorder_0cellp_1d * calculus_1d_manual.laplace<DUAL>() * reorder_0cellp_1d.transpose(),
410 reorder_0cellp_1d_factory * calculus_1d_factory.laplace<DUAL>() * reorder_0cellp_1d_factory.transpose()
411 ) );
412 FATAL_ERROR( equal(
413 reorder_0cellp_2d * calculus_2d_manual.laplace<DUAL>() * reorder_0cellp_2d.transpose(),
414 reorder_0cellp_2d_factory * calculus_2d_factory.laplace<DUAL>() * reorder_0cellp_2d_factory.transpose()
415 ) );
416 FATAL_ERROR( equal(
417 reorder_0cellp_3d * calculus_3d_manual.laplace<DUAL>() * reorder_0cellp_3d.transpose(),
418 reorder_0cellp_3d_factory * calculus_3d_factory.laplace<DUAL>() * reorder_0cellp_3d_factory.transpose()
419 ) );
420 }
421
422 trace.endBlock();
423
424 trace.beginBlock("checking border");
425
426 { // 2d ambient border
427 const Calculus2D calculus_2d_factory_no_border = CalculusFactory::createFromNSCells<1>(ncells_2d_factory.begin(), ncells_2d_factory.end(), false);
428 trace.info() << "calculus_2d_factory_no_border=" << calculus_2d_factory_no_border << endl;
429 trace.info() << "calculus_2d_factory=" << calculus_2d_factory << endl;
430 FATAL_ERROR( calculus_2d_factory.containsCell(calculus_2d_factory.myKSpace.uCell(Z2i::Point(6,0))) );
431 FATAL_ERROR( !calculus_2d_factory_no_border.containsCell(calculus_2d_factory_no_border.myKSpace.uCell(Z2i::Point(6,0))) );
432 FATAL_ERROR( calculus_2d_factory.containsCell(calculus_2d_factory.myKSpace.uCell(Z2i::Point(0,0))) );
433 FATAL_ERROR( !calculus_2d_factory_no_border.containsCell(calculus_2d_factory_no_border.myKSpace.uCell(Z2i::Point(0,0))) );
434 }
435
436 { // 3d ambient border
437 const Calculus3D calculus_3d_factory_no_border = CalculusFactory::createFromNSCells<1>(ncells_3d_factory.begin(), ncells_3d_factory.end(), false);
438 trace.info() << "calculus_3d_factory_no_border=" << calculus_3d_factory_no_border << endl;
439 trace.info() << "calculus_3d_factory=" << calculus_3d_factory << endl;
440 FATAL_ERROR( calculus_3d_factory.containsCell(calculus_3d_factory.myKSpace.uCell(Z3i::Point(2,2,-2))) );
441 FATAL_ERROR( !calculus_3d_factory_no_border.containsCell(calculus_3d_factory_no_border.myKSpace.uCell(Z3i::Point(2,2,-2))) );
442 FATAL_ERROR( calculus_3d_factory.containsCell(calculus_3d_factory.myKSpace.uCell(Z3i::Point(0,0,0))) );
443 FATAL_ERROR( !calculus_3d_factory_no_border.containsCell(calculus_3d_factory_no_border.myKSpace.uCell(Z3i::Point(0,0,0))) );
444 }
445
446 trace.endBlock();
447
448 trace.beginBlock("checking area");
449
450 {
451 const Calculus2D::Scalar area_th = calculus_2d_factory.kFormLength(0, DUAL);
452 const Calculus2D::Scalar area_primal = (
453 calculus_2d_factory.hodge<0, PRIMAL>() *
454 Calculus2D::PrimalForm0::ones(calculus_2d_factory)
455 ).myContainer.array().sum();
456 const Calculus2D::Scalar area_dual = (
457 calculus_2d_factory.hodge<0, DUAL>() *
458 Calculus2D::DualForm0::ones(calculus_2d_factory)
459 ).myContainer.array().sum();
460 trace.info() << "area_2d_th=" << area_th << endl;
461 trace.info() << "area_2d_primal=" << area_primal << endl;
462 trace.info() << "area_2d_dual=" << area_dual << endl;
463 FATAL_ERROR( area_th == area_primal );
464 FATAL_ERROR( area_th == area_dual );
465 }
466
467 {
468 const Calculus3D::Scalar area_th = calculus_3d_factory.kFormLength(0, DUAL);
469 const Calculus3D::Scalar area_primal = (
470 calculus_3d_factory.hodge<0, PRIMAL>() *
471 Calculus3D::PrimalForm0::ones(calculus_3d_factory)
472 ).myContainer.array().sum();
473 const Calculus3D::Scalar area_dual = (
474 calculus_3d_factory.hodge<0, DUAL>() *
475 Calculus3D::DualForm0::ones(calculus_3d_factory)
476 ).myContainer.array().sum();
477 trace.info() << "area_3d_th=" << area_th << endl;
478 trace.info() << "area_3d_primal=" << area_primal << endl;
479 trace.info() << "area_3d_dual=" << area_dual << endl;
480 FATAL_ERROR( area_th == area_primal );
481 FATAL_ERROR( area_th == area_dual );
482 }
483
484 trace.endBlock();
485
486 trace.endBlock();
487 }
488
489 /*{
490 trace.beginBlock("2d manifold embedding");
491
493 typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus2D;
494 typedef DiscreteExteriorCalculus<2, 3, EigenLinearAlgebraBackend> Calculus3D;
496
497 Calculus2D calculus_2d_manual;
498 {
499 for (int xx=0; xx<=4; xx++)
500 for (int yy=0; yy<=4; yy++)
501 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy)) );
502
503 for (int xx=0; xx<=4; xx++)
504 for (int yy=0; yy<=4; yy++)
505 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+4)) );
506
507 for (int xx=0; xx<=4; xx++)
508 for (int yy=0; yy<=4; yy++)
509 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+8)) );
510
511 for (int xx=0; xx<=4; xx++)
512 for (int yy=0; yy<=4; yy++)
513 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+12)) );
514
515 for (int xx=0; xx<=4; xx++)
516 for (int yy=0; yy<=4; yy++)
517 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+16)) );
518
519 for (int xx=0; xx<=4; xx++)
520 for (int yy=0; yy<=4; yy++)
521 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+4,yy)) );
522
523 for (int xx=0; xx<=4; xx++)
524 for (int yy=0; yy<=4; yy++)
525 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+8,yy)) );
526
527 for (int xx=0; xx<=4; xx++)
528 for (int yy=0; yy<=4; yy++)
529 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+12,yy)) );
530
531 for (int xx=0; xx<=4; xx++)
532 for (int yy=0; yy<=4; yy++)
533 calculus_2d_manual.insertSCell( calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+16,yy)) );
534
535 calculus_2d_manual.updateIndexes();
536 }
537 trace.info() << "calculus_2d_manual=" << calculus_2d_manual << endl;
538
539 {
540 Board2D board;
541 board << Z2i::Domain(Z2i::Point(-1,-1), Z2i::Point(10,10));
542 board << calculus_2d_manual;
543 board.saveSVG("embedding_2d_calculus_2d.svg");
544 }
545
547 typedef std::list<Calculus2D::SCell> SCells2D;
548 SCells2D ncells_2d_factory;
550 for (int kk=0; kk<calculus_2d_manual.kFormLength(2, PRIMAL); kk++) ncells_2d_factory.push_back( calculus_2d_manual.getSCell(2, PRIMAL, kk) );
552 const Calculus2D calculus_2d_factory_weighed = CalculusFactory::createFromNSCells<2>(ncells_2d_factory.begin(), ncells_2d_factory.end(), true);
554 Calculus2D calculus_2d_factory = calculus_2d_factory_weighed;
555 calculus_2d_factory.resetSizes();
556 trace.info() << "calculus_2d_factory=" << calculus_2d_factory << endl;
557
558 SCells2D cells_2d_manual;
559 {
560 std::set<Calculus2D::SCell> already_inserted;
561
562 for (int xx=0; xx<=4; xx++)
563 for (int yy=0; yy<=4; yy++)
564 {
565 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy));
566 if (already_inserted.find(cell) != already_inserted.end()) continue;
567 already_inserted.insert(cell);
568 cells_2d_manual.push_back(cell);
569 }
570
571 for (int xx=0; xx<=4; xx++)
572 for (int yy=0; yy<=4; yy++)
573 {
574 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+4));
575 if (already_inserted.find(cell) != already_inserted.end()) continue;
576 already_inserted.insert(cell);
577 cells_2d_manual.push_back(cell);
578 }
579
580 for (int xx=0; xx<=4; xx++)
581 for (int yy=0; yy<=4; yy++)
582 {
583 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+8));
584 if (already_inserted.find(cell) != already_inserted.end()) continue;
585 already_inserted.insert(cell);
586 cells_2d_manual.push_back(cell);
587 }
588
589 for (int xx=0; xx<=4; xx++)
590 for (int yy=0; yy<=4; yy++)
591 {
592 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+12));
593 if (already_inserted.find(cell) != already_inserted.end()) continue;
594 already_inserted.insert(cell);
595 cells_2d_manual.push_back(cell);
596 }
597
598 for (int xx=0; xx<=4; xx++)
599 for (int yy=0; yy<=4; yy++)
600 {
601 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx,yy+16));
602 if (already_inserted.find(cell) != already_inserted.end()) continue;
603 already_inserted.insert(cell);
604 cells_2d_manual.push_back(cell);
605 }
606
607 for (int xx=0; xx<=4; xx++)
608 for (int yy=0; yy<=4; yy++)
609 {
610 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+4,yy));
611 if (already_inserted.find(cell) != already_inserted.end()) continue;
612 already_inserted.insert(cell);
613 cells_2d_manual.push_back(cell);
614 }
615
616 for (int xx=0; xx<=4; xx++)
617 for (int yy=0; yy<=4; yy++)
618 {
619 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+8,yy));
620 if (already_inserted.find(cell) != already_inserted.end()) continue;
621 already_inserted.insert(cell);
622 cells_2d_manual.push_back(cell);
623 }
624
625 for (int xx=0; xx<=4; xx++)
626 for (int yy=0; yy<=4; yy++)
627 {
628 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+12,yy));
629 if (already_inserted.find(cell) != already_inserted.end()) continue;
630 already_inserted.insert(cell);
631 cells_2d_manual.push_back(cell);
632 }
633
634 for (int xx=0; xx<=4; xx++)
635 for (int yy=0; yy<=4; yy++)
636 {
637 const Calculus2D::SCell cell = calculus_2d_manual.myKSpace.sCell(Z2i::Point(xx+16,yy));
638 if (already_inserted.find(cell) != already_inserted.end()) continue;
639 already_inserted.insert(cell);
640 cells_2d_manual.push_back(cell);
641 }
642 }
643
644 Calculus3D calculus_3d_manual;
645 {
646 for (int xx=0; xx<=4; xx++)
647 for (int yy=0; yy<=4; yy++)
648 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,yy,0)) );
649
650 for (int xx=0; xx<=4; xx++)
651 for (int yy=0; yy<=4; yy++)
652 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,4,yy)) );
653
654 for (int xx=0; xx<=4; xx++)
655 for (int yy=0; yy<=4; yy++)
656 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,4-yy,4),
657 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
658 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
659 Calculus3D::KSpace::POS) );
660
661 for (int xx=0; xx<=4; xx++)
662 for (int yy=0; yy<=4; yy++)
663 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,-yy,4),
664 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
665 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
666 Calculus3D::KSpace::POS) );
667
668 for (int xx=0; xx<=4; xx++)
669 for (int yy=0; yy<=4; yy++)
670 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,-4,4-yy),
671 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
672 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
673 Calculus3D::KSpace::POS) );
674
675 for (int xx=0; xx<=4; xx++)
676 for (int yy=0; yy<=4; yy++)
677 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,yy,-xx),
678 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::POS : // surfels
679 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
680 Calculus3D::KSpace::POS) );
681
682 for (int xx=0; xx<=4; xx++)
683 for (int yy=0; yy<=4; yy++)
684 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(4-xx,yy,-4),
685 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
686 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
687 Calculus3D::KSpace::POS) );
688
689 for (int xx=0; xx<=4; xx++)
690 for (int yy=0; yy<=4; yy++)
691 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(-xx,yy,-4),
692 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
693 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
694 Calculus3D::KSpace::POS) );
695
696 for (int xx=0; xx<=4; xx++)
697 for (int yy=0; yy<=4; yy++)
698 calculus_3d_manual.insertSCell( calculus_3d_manual.myKSpace.sCell(Z3i::Point(-4,yy,-4+xx),
699 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
700 xx%2 != 0 ? Calculus3D::KSpace::POS : // x-edge
701 Calculus3D::KSpace::POS) );
702
703 calculus_3d_manual.updateIndexes();
704 }
705 trace.info() << "calculus_3d_manual=" << calculus_3d_manual << endl;
706
707#if !defined(NOVIEWER)
708 Display3DFactory<Calculus3D::KSpace::Space, Calculus3D::KSpace>::draw(viewer2, calculus_3d_manual);
709 viewer2 << Viewer::updateDisplay;
710#endif
711
713 typedef std::list<Calculus3D::SCell> SCells3D;
714 SCells3D ncells_3d_factory;
716 for (int kk=0; kk<calculus_3d_manual.kFormLength(2, PRIMAL); kk++) ncells_3d_factory.push_back( calculus_3d_manual.getSCell(2, PRIMAL, kk) );
718 const Calculus3D calculus_3d_factory_weighed = CalculusFactory::createFromNSCells<2>(ncells_3d_factory.begin(), ncells_3d_factory.end(), true);
720 Calculus3D calculus_3d_factory = calculus_3d_factory_weighed;
721 calculus_3d_factory.resetSizes();
722 trace.info() << "calculus_3d_factory=" << calculus_3d_factory << endl;
723
724 SCells3D cells_3d_manual;
725 {
726 std::set<Calculus3D::SCell> already_inserted;
727
728 for (int xx=0; xx<=4; xx++)
729 for (int yy=0; yy<=4; yy++)
730 {
731 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,yy,0));
732 if (already_inserted.find(cell) != already_inserted.end()) continue;
733 already_inserted.insert(cell);
734 cells_3d_manual.push_back(cell);
735 }
736
737 for (int xx=0; xx<=4; xx++)
738 for (int yy=0; yy<=4; yy++)
739 {
740 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,4,yy));
741 if (already_inserted.find(cell) != already_inserted.end()) continue;
742 already_inserted.insert(cell);
743 cells_3d_manual.push_back(cell);
744 }
745
746 for (int xx=0; xx<=4; xx++)
747 for (int yy=0; yy<=4; yy++)
748 {
749 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,4-yy,4),
750 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
751 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
752 Calculus3D::KSpace::POS);
753 if (already_inserted.find(cell) != already_inserted.end()) continue;
754 already_inserted.insert(cell);
755 cells_3d_manual.push_back(cell);
756 }
757
758 for (int xx=0; xx<=4; xx++)
759 for (int yy=0; yy<=4; yy++)
760 {
761 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,-yy,4),
762 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
763 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
764 Calculus3D::KSpace::POS);
765 if (already_inserted.find(cell) != already_inserted.end()) continue;
766 already_inserted.insert(cell);
767 cells_3d_manual.push_back(cell);
768 }
769
770 for (int xx=0; xx<=4; xx++)
771 for (int yy=0; yy<=4; yy++)
772 {
773 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(xx,-4,4-yy),
774 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
775 yy%2 != 0 ? Calculus3D::KSpace::NEG : // y-edge
776 Calculus3D::KSpace::POS);
777 if (already_inserted.find(cell) != already_inserted.end()) continue;
778 already_inserted.insert(cell);
779 cells_3d_manual.push_back(cell);
780 }
781
782 for (int xx=0; xx<=4; xx++)
783 for (int yy=0; yy<=4; yy++)
784 {
785 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(4,yy,-xx),
786 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::POS : // surfels
787 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
788 Calculus3D::KSpace::POS);
789 if (already_inserted.find(cell) != already_inserted.end()) continue;
790 already_inserted.insert(cell);
791 cells_3d_manual.push_back(cell);
792 }
793
794 for (int xx=0; xx<=4; xx++)
795 for (int yy=0; yy<=4; yy++)
796 {
797 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(4-xx,yy,-4),
798 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
799 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
800 Calculus3D::KSpace::POS);
801 if (already_inserted.find(cell) != already_inserted.end()) continue;
802 already_inserted.insert(cell);
803 cells_3d_manual.push_back(cell);
804 }
805
806 for (int xx=0; xx<=4; xx++)
807 for (int yy=0; yy<=4; yy++)
808 {
809 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(-xx,yy,-4),
810 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
811 xx%2 != 0 ? Calculus3D::KSpace::NEG : // x-edge
812 Calculus3D::KSpace::POS);
813 if (already_inserted.find(cell) != already_inserted.end()) continue;
814 already_inserted.insert(cell);
815 cells_3d_manual.push_back(cell);
816 }
817
818 for (int xx=0; xx<=4; xx++)
819 for (int yy=0; yy<=4; yy++)
820 {
821 const Calculus3D::SCell cell = calculus_3d_manual.myKSpace.sCell(Z3i::Point(-4,yy,-4+xx),
822 xx%2 != 0 && yy%2 != 0 ? Calculus3D::KSpace::NEG : // surfels
823 xx%2 != 0 ? Calculus3D::KSpace::POS : // x-edge
824 Calculus3D::KSpace::POS);
825 if (already_inserted.find(cell) != already_inserted.end()) continue;
826 already_inserted.insert(cell);
827 cells_3d_manual.push_back(cell);
828 }
829 }
830
831 trace.beginBlock("checking operators");
832 trace.info() << calculus_3d_manual << endl;
833 trace.info() << cells_3d_manual.size() << endl;
834
835 const Calculus2D::PrimalIdentity0 reorder_0cell_2d = calculus_2d_manual.reorder<0, PRIMAL>(cells_2d_manual.begin(), cells_2d_manual.end());
836 const Calculus3D::PrimalIdentity0 reorder_0cell_3d = calculus_3d_manual.reorder<0, PRIMAL>(cells_3d_manual.begin(), cells_3d_manual.end());
837 const Calculus2D::PrimalIdentity1 reorder_1cell_2d = calculus_2d_manual.reorder<1, PRIMAL>(cells_2d_manual.begin(), cells_2d_manual.end());
838 const Calculus3D::PrimalIdentity1 reorder_1cell_3d = calculus_3d_manual.reorder<1, PRIMAL>(cells_3d_manual.begin(), cells_3d_manual.end());
839 const Calculus2D::PrimalIdentity2 reorder_2cell_2d = calculus_2d_manual.reorder<2, PRIMAL>(cells_2d_manual.begin(), cells_2d_manual.end());
840 const Calculus3D::PrimalIdentity2 reorder_2cell_3d = calculus_3d_manual.reorder<2, PRIMAL>(cells_3d_manual.begin(), cells_3d_manual.end());
841 const Calculus2D::DualIdentity0 reorder_0cellp_2d = calculus_2d_manual.reorder<0, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
842 const Calculus3D::DualIdentity0 reorder_0cellp_3d = calculus_3d_manual.reorder<0, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
843 const Calculus2D::DualIdentity1 reorder_1cellp_2d = calculus_2d_manual.reorder<1, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
844 const Calculus3D::DualIdentity1 reorder_1cellp_3d = calculus_3d_manual.reorder<1, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
845 const Calculus2D::DualIdentity2 reorder_2cellp_2d = calculus_2d_manual.reorder<2, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
846 const Calculus3D::DualIdentity2 reorder_2cellp_3d = calculus_3d_manual.reorder<2, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
847
848 { // check primal operators
849 const Calculus2D::PrimalIdentity0 primal_laplace_2d = reorder_0cell_2d * calculus_2d_manual.laplace<PRIMAL>() * reorder_0cell_2d.transpose();
850 const Calculus3D::PrimalIdentity0 primal_laplace_3d = reorder_0cell_3d * calculus_3d_manual.laplace<PRIMAL>() * reorder_0cell_3d.transpose();
851 trace.info() << "primal_laplace_2d=" << primal_laplace_2d << endl;
852 trace.info() << "primal_laplace_3d=" << primal_laplace_3d << endl;
853 trace.info() << "primal_laplace_container=" << endl << MatrixXd(primal_laplace_2d.myContainer) << endl;
854 FATAL_ERROR( equal(
855 reorder_2cellp_2d * calculus_2d_manual.hodge<0, PRIMAL>() * reorder_0cell_2d.transpose(),
856 reorder_2cellp_3d * calculus_3d_manual.hodge<0, PRIMAL>() * reorder_0cell_3d.transpose()
857 ) );
858 FATAL_ERROR( equal(
859 reorder_1cellp_2d * calculus_2d_manual.hodge<1, PRIMAL>() * reorder_1cell_2d.transpose(),
860 reorder_1cellp_3d * calculus_3d_manual.hodge<1, PRIMAL>() * reorder_1cell_3d.transpose()
861 ) );
862 FATAL_ERROR( equal(
863 reorder_0cellp_2d * calculus_2d_manual.hodge<2, PRIMAL>() * reorder_2cell_2d.transpose(),
864 reorder_0cellp_3d * calculus_3d_manual.hodge<2, PRIMAL>() * reorder_2cell_3d.transpose()
865 ) );
866 FATAL_ERROR( equal(
867 reorder_1cell_2d * calculus_2d_manual.derivative<0, PRIMAL>() * reorder_0cell_2d.transpose(),
868 reorder_1cell_3d * calculus_3d_manual.derivative<0, PRIMAL>() * reorder_0cell_3d.transpose()
869 ) );
870 FATAL_ERROR( equal(
871 reorder_2cell_2d * calculus_2d_manual.derivative<1, PRIMAL>() * reorder_1cell_2d.transpose(),
872 reorder_2cell_3d * calculus_3d_manual.derivative<1, PRIMAL>() * reorder_1cell_3d.transpose()
873 ) );
874 FATAL_ERROR( equal(primal_laplace_2d, primal_laplace_3d) );
875 }
876
877 { // check dual operators
878 const Calculus2D::DualIdentity0 dual_laplace_2d = reorder_0cellp_2d * calculus_2d_manual.laplace<DUAL>() * reorder_0cellp_2d.transpose();
879 const Calculus3D::DualIdentity0 dual_laplace_3d = reorder_0cellp_3d * calculus_3d_manual.laplace<DUAL>() * reorder_0cellp_3d.transpose();
880 trace.info() << "dual_laplace_2d=" << dual_laplace_2d << endl;
881 trace.info() << "dual_laplace_3d=" << dual_laplace_3d << endl;
882 trace.info() << "dual_laplace_container=" << endl << MatrixXd(dual_laplace_2d.myContainer) << endl;
883 FATAL_ERROR( equal(
884 reorder_2cell_2d * calculus_2d_manual.hodge<0, DUAL>() * reorder_0cellp_2d.transpose(),
885 reorder_2cell_3d * calculus_3d_manual.hodge<0, DUAL>() * reorder_0cellp_3d.transpose()
886 ) );
887 FATAL_ERROR( equal(
888 reorder_1cell_2d * calculus_2d_manual.hodge<1, DUAL>() * reorder_1cellp_2d.transpose(),
889 reorder_1cell_3d * calculus_3d_manual.hodge<1, DUAL>() * reorder_1cellp_3d.transpose()
890 ) );
891 FATAL_ERROR( equal(
892 reorder_0cell_2d * calculus_2d_manual.hodge<2, DUAL>() * reorder_2cellp_2d.transpose(),
893 reorder_0cell_3d * calculus_3d_manual.hodge<2, DUAL>() * reorder_2cellp_3d.transpose()
894 ) );
895 FATAL_ERROR( equal(
896 reorder_1cellp_2d * calculus_2d_manual.derivative<0, DUAL>() * reorder_0cellp_2d.transpose(),
897 reorder_1cellp_3d * calculus_3d_manual.derivative<0, DUAL>() * reorder_0cellp_3d.transpose()
898 ) );
899 FATAL_ERROR( equal(
900 reorder_2cellp_2d * calculus_2d_manual.derivative<1, DUAL>() * reorder_1cellp_2d.transpose(),
901 reorder_2cellp_3d * calculus_3d_manual.derivative<1, DUAL>() * reorder_1cellp_3d.transpose()
902 ) );
903 FATAL_ERROR( equal(dual_laplace_2d, dual_laplace_3d) );
904 }
905
906 { // checking dual laplace factory calculus vs manual calculus
907 const Calculus2D::DualIdentity0 reorder_0cellp_2d_factory = calculus_2d_factory.reorder<0, DUAL>(cells_2d_manual.begin(), cells_2d_manual.end());
908 const Calculus3D::DualIdentity0 reorder_0cellp_3d_factory = calculus_3d_factory.reorder<0, DUAL>(cells_3d_manual.begin(), cells_3d_manual.end());
909 FATAL_ERROR( equal(
910 reorder_0cellp_2d * calculus_2d_manual.laplace<DUAL>() * reorder_0cellp_2d.transpose(),
911 reorder_0cellp_2d_factory * calculus_2d_factory.laplace<DUAL>() * reorder_0cellp_2d_factory.transpose()
912 ) );
913 FATAL_ERROR( equal(
914 reorder_0cellp_3d * calculus_3d_manual.laplace<DUAL>() * reorder_0cellp_3d.transpose(),
915 reorder_0cellp_3d_factory * calculus_3d_factory.laplace<DUAL>() * reorder_0cellp_3d_factory.transpose()
916 ) );
917 }
918
919 trace.endBlock();
920
921 trace.beginBlock("checking border");
922
923 { // 2d ambient border
925 const Calculus2D calculus_2d_factory_no_border = CalculusFactory::createFromNSCells<2>(ncells_2d_factory.begin(), ncells_2d_factory.end(), false);
927 trace.info() << "calculus_2d_factory_no_border=" << calculus_2d_factory_no_border << endl;
928 trace.info() << "calculus_2d_factory=" << calculus_2d_factory << endl;
929
930 Board2D board;
931 board << Z2i::Domain(Z2i::Point(-1,-1), Z2i::Point(10,10));
932 board << calculus_2d_factory_no_border;
933 board.saveSVG("embedding_2d_calculus_2d_no_border.svg");
934 }
935
936 { // 3d ambient border
938 const Calculus3D calculus_3d_factory_no_border = CalculusFactory::createFromNSCells<2>(ncells_3d_factory.begin(), ncells_3d_factory.end(), false);
940 trace.info() << "calculus_3d_factory_no_border=" << calculus_3d_factory_no_border << endl;
941 trace.info() << "calculus_3d_factory=" << calculus_3d_factory << endl;
942
943#if !defined(NOVIEWER)
944 Display3DFactory<Calculus3D::KSpace::Space, Calculus3D::KSpace>::draw(viewer3, calculus_3d_factory_no_border);
945 viewer3 << Viewer::updateDisplay;
946#endif
947 }
948
949 trace.endBlock();
950
951 trace.beginBlock("checking area");
952
953 {
954 const Calculus2D::Scalar area_th = calculus_2d_factory_weighed.kFormLength(0, DUAL);
955 const Calculus2D::Scalar area_primal = (
956 calculus_2d_factory_weighed.hodge<0, PRIMAL>() *
957 Calculus2D::PrimalForm0::ones(calculus_2d_factory_weighed)
958 ).myContainer.array().sum();
959 const Calculus2D::Scalar area_dual = (
960 calculus_2d_factory_weighed.hodge<0, DUAL>() *
961 Calculus2D::DualForm0::ones(calculus_2d_factory_weighed)
962 ).myContainer.array().sum();
963 trace.info() << "area_2d_th=" << area_th << endl;
964 trace.info() << "area_2d_primal=" << area_primal << endl;
965 trace.info() << "area_2d_dual=" << area_dual << endl;
966 FATAL_ERROR( area_th == area_primal );
967 FATAL_ERROR( area_th == area_dual );
968 }
969
970 {
971 const Calculus3D::Scalar area_th = calculus_3d_factory_weighed.kFormLength(0, DUAL);
972 const Calculus3D::Scalar area_primal = (
973 calculus_3d_factory_weighed.hodge<0, PRIMAL>() *
974 Calculus3D::PrimalForm0::ones(calculus_3d_factory_weighed)
975 ).myContainer.array().sum();
976 const Calculus3D::Scalar area_dual = (
977 calculus_3d_factory_weighed.hodge<0, DUAL>() *
978 Calculus3D::DualForm0::ones(calculus_3d_factory_weighed)
979 ).myContainer.array().sum();
980 trace.info() << "area_3d_th=" << area_th << endl;
981 trace.info() << "area_3d_primal=" << area_primal << endl;
982 trace.info() << "area_3d_dual=" << area_dual << endl;
983 FATAL_ERROR( area_th == area_primal );
984 FATAL_ERROR( area_th == area_dual );
985 }
986
987 trace.endBlock();
988
989 trace.endBlock();
990 }*/
991
992#if !defined(NOVIEWER)
993 return app.exec();
994#else
995 return 0;
996#endif
997}
998
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
Aim: This class provides static members to create DEC structures from various other DGtal structures.
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
SCell sCell(const SPreCell &c) const
From a signed cell, returns a signed cell lying into this Khalismky space.
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
Board & setFillColor(const DGtal::Color &color)
Definition: Board.cpp:322
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
HyperRectDomain< Space > Domain
Definition: StdDefs.h:99
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
@ PRIMAL
Definition: Duality.h:61
@ DUAL
Definition: Duality.h:62
STL namespace.
static void drawDECSignedKhalimskyCell(DGtal::Board2D &board, const DGtal::SignedKhalimskyCell< dim, TInteger > &cell)
static void draw(Display3D< Space, KSpace > &display, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
int main()
Definition: testBits.cpp:56