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"
9#include "DGtal/io/viewers/Viewer3D.h"
11#include "DGtal/io/boards/Board2D.h"
12#include "DGtal/base/Common.h"
13#include "DGtal/helpers/StdDefs.h"
19template <
typename OperatorAA,
typename OperatorBB>
21equal(
const OperatorAA& aa,
const OperatorBB& bb)
23 return MatrixXd(aa.myContainer) == MatrixXd(bb.myContainer);
26int main(
int ,
char** )
36 QApplication app(argc, argv);
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) );
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) );
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) );
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++)
81 Calculus1D::KSpace::Point point;
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);
88 calculus_1d_manual.updateIndexes();
89 trace.
info() <<
"calculus_1d_manual=" << calculus_1d_manual << endl;
95 for (Calculus1D::ConstIterator ii=calculus_1d_manual.begin(), ie=calculus_1d_manual.end(); ii!=ie; ii++)
97 const Z2i::Point point(calculus_1d_manual.myKSpace.uKCoord(ii->first, 0), 0);
101 board.
saveSVG(
"embedding_1d_calculus_1d.svg");
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;
109 Calculus2D calculus_2d_manual;
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();
144 trace.
info() <<
"calculus_2d_manual=" << calculus_2d_manual << endl;
149 board << calculus_2d_manual;
150 board.
saveSVG(
"embedding_1d_calculus_2d.svg");
154 typedef std::list<Calculus2D::SCell> SCells2D;
155 SCells2D ncells_2d_factory;
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) );
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;
180 SCells2D cells_2d_manual;
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)) );
215 Calculus3D calculus_3d_manual;
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();
250 trace.
info() <<
"calculus_3d_manual=" << calculus_3d_manual << endl;
252#if !defined(NOVIEWER)
258 typedef std::vector<Calculus3D::SCell> SCells3D;
259 SCells3D ncells_3d_factory;
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) );
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;
283 SCells3D cells_3d_manual;
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)) );
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());
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;
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()
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()
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()
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()
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()
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()
364 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_2d) );
365 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_3d) );
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;
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()
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()
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()
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()
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()
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()
400 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_2d) );
401 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_3d) );
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());
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()
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()
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()
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))) );
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))) );
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 );
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 );
493 typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus2D;
494 typedef DiscreteExteriorCalculus<2, 3, EigenLinearAlgebraBackend> Calculus3D;
497 Calculus2D calculus_2d_manual;
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
535 calculus_2d_manual.updateIndexes();
537 trace.info() << "calculus_2d_manual=" << calculus_2d_manual << endl;
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");
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;
558 SCells2D cells_2d_manual;
560 std::set<Calculus2D::SCell> already_inserted;
562 for (int xx=0; xx<=4; xx++)
563 for (int yy=0; yy<=4; yy++)
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);
571 for (int xx=0; xx<=4; xx++)
572 for (int yy=0; yy<=4; yy++)
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);
580 for (int xx=0; xx<=4; xx++)
581 for (int yy=0; yy<=4; yy++)
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);
589 for (int xx=0; xx<=4; xx++)
590 for (int yy=0; yy<=4; yy++)
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);
598 for (int xx=0; xx<=4; xx++)
599 for (int yy=0; yy<=4; yy++)
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);
607 for (int xx=0; xx<=4; xx++)
608 for (int yy=0; yy<=4; yy++)
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);
616 for (int xx=0; xx<=4; xx++)
617 for (int yy=0; yy<=4; yy++)
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);
625 for (int xx=0; xx<=4; xx++)
626 for (int yy=0; yy<=4; yy++)
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);
634 for (int xx=0; xx<=4; xx++)
635 for (int yy=0; yy<=4; yy++)
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);
644 Calculus3D calculus_3d_manual;
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)) );
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)) );
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) );
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) );
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) );
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) );
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) );
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) );
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) );
703 calculus_3d_manual.updateIndexes();
705 trace.info() << "calculus_3d_manual=" << calculus_3d_manual << endl;
707#if !defined(NOVIEWER)
708 Display3DFactory<Calculus3D::KSpace::Space, Calculus3D::KSpace>::draw(viewer2, calculus_3d_manual);
709 viewer2 << Viewer::updateDisplay;
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;
724 SCells3D cells_3d_manual;
726 std::set<Calculus3D::SCell> already_inserted;
728 for (int xx=0; xx<=4; xx++)
729 for (int yy=0; yy<=4; yy++)
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);
737 for (int xx=0; xx<=4; xx++)
738 for (int yy=0; yy<=4; yy++)
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);
746 for (int xx=0; xx<=4; xx++)
747 for (int yy=0; yy<=4; yy++)
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);
758 for (int xx=0; xx<=4; xx++)
759 for (int yy=0; yy<=4; yy++)
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);
770 for (int xx=0; xx<=4; xx++)
771 for (int yy=0; yy<=4; yy++)
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);
782 for (int xx=0; xx<=4; xx++)
783 for (int yy=0; yy<=4; yy++)
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);
794 for (int xx=0; xx<=4; xx++)
795 for (int yy=0; yy<=4; yy++)
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);
806 for (int xx=0; xx<=4; xx++)
807 for (int yy=0; yy<=4; yy++)
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);
818 for (int xx=0; xx<=4; xx++)
819 for (int yy=0; yy<=4; yy++)
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);
831 trace.beginBlock("checking operators");
832 trace.info() << calculus_3d_manual << endl;
833 trace.info() << cells_3d_manual.size() << endl;
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());
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;
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()
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()
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()
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()
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()
874 FATAL_ERROR( equal(primal_laplace_2d, primal_laplace_3d) );
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;
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()
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()
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()
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()
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()
903 FATAL_ERROR( equal(dual_laplace_2d, dual_laplace_3d) );
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());
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()
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()
921 trace.beginBlock("checking border");
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;
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");
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;
943#if !defined(NOVIEWER)
944 Display3DFactory<Calculus3D::KSpace::Space, Calculus3D::KSpace>::draw(viewer3, calculus_3d_factory_no_border);
945 viewer3 << Viewer::updateDisplay;
951 trace.beginBlock("checking area");
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 );
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 );
992#if !defined(NOVIEWER)
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Structure representing an RGB triple with alpha component.
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="")
Board & setFillColor(const DGtal::Color &color)
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
HyperRectDomain< Space > Domain
DGtal is the top-level namespace which contains all DGtal functions and types.
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.