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