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/PolyscopeViewer.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** )
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++)
63 Calculus1D::KSpace::Point point;
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);
70 calculus_1d_manual.updateIndexes();
71 trace.
info() <<
"calculus_1d_manual=" << calculus_1d_manual << endl;
77 for (Calculus1D::ConstIterator ii=calculus_1d_manual.begin(), ie=calculus_1d_manual.end(); ii!=ie; ii++)
79 const Z2i::Point point(calculus_1d_manual.myKSpace.uKCoord(ii->first, 0), 0);
83 board.
saveSVG(
"embedding_1d_calculus_1d.svg");
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;
91 Calculus2D calculus_2d_manual;
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();
126 trace.
info() <<
"calculus_2d_manual=" << calculus_2d_manual << endl;
131 board << calculus_2d_manual;
132 board.
saveSVG(
"embedding_1d_calculus_2d.svg");
136 typedef std::list<Calculus2D::SCell> SCells2D;
137 SCells2D ncells_2d_factory;
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) );
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;
162 SCells2D cells_2d_manual;
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)) );
197 Calculus3D calculus_3d_manual;
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();
232 trace.
info() <<
"calculus_3d_manual=" << calculus_3d_manual << endl;
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});
241 typedef std::vector<Calculus3D::SCell> SCells3D;
242 SCells3D ncells_3d_factory;
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) );
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;
266 SCells3D cells_3d_manual;
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)) );
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());
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;
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()
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()
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()
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()
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()
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()
347 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_2d) );
348 FATAL_ERROR( equal(primal_laplace_1d, primal_laplace_3d) );
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;
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()
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()
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()
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()
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()
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()
383 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_2d) );
384 FATAL_ERROR( equal(dual_laplace_1d, dual_laplace_3d) );
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());
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()
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()
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()
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))) );
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))) );
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 );
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 );
476 typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus2D;
477 typedef DiscreteExteriorCalculus<2, 3, EigenLinearAlgebraBackend> Calculus3D;
480 Calculus2D calculus_2d_manual;
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
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)) );
518 calculus_2d_manual.updateIndexes();
520 trace.info() << "calculus_2d_manual=" << calculus_2d_manual << endl;
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");
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;
541 SCells2D cells_2d_manual;
543 std::set<Calculus2D::SCell> already_inserted;
545 for (int xx=0; xx<=4; xx++)
546 for (int yy=0; yy<=4; yy++)
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);
554 for (int xx=0; xx<=4; xx++)
555 for (int yy=0; yy<=4; yy++)
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);
563 for (int xx=0; xx<=4; xx++)
564 for (int yy=0; yy<=4; yy++)
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);
572 for (int xx=0; xx<=4; xx++)
573 for (int yy=0; yy<=4; yy++)
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);
581 for (int xx=0; xx<=4; xx++)
582 for (int yy=0; yy<=4; yy++)
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);
590 for (int xx=0; xx<=4; xx++)
591 for (int yy=0; yy<=4; yy++)
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);
599 for (int xx=0; xx<=4; xx++)
600 for (int yy=0; yy<=4; yy++)
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);
608 for (int xx=0; xx<=4; xx++)
609 for (int yy=0; yy<=4; yy++)
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);
617 for (int xx=0; xx<=4; xx++)
618 for (int yy=0; yy<=4; yy++)
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);
627 Calculus3D calculus_3d_manual;
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)) );
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)) );
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) );
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) );
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) );
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) );
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) );
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) );
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) );
686 calculus_3d_manual.updateIndexes();
688 trace.info() << "calculus_3d_manual=" << calculus_3d_manual << endl;
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});
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;
708 SCells3D cells_3d_manual;
710 std::set<Calculus3D::SCell> already_inserted;
712 for (int xx=0; xx<=4; xx++)
713 for (int yy=0; yy<=4; yy++)
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);
721 for (int xx=0; xx<=4; xx++)
722 for (int yy=0; yy<=4; yy++)
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);
730 for (int xx=0; xx<=4; xx++)
731 for (int yy=0; yy<=4; yy++)
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);
742 for (int xx=0; xx<=4; xx++)
743 for (int yy=0; yy<=4; yy++)
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);
754 for (int xx=0; xx<=4; xx++)
755 for (int yy=0; yy<=4; yy++)
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);
766 for (int xx=0; xx<=4; xx++)
767 for (int yy=0; yy<=4; yy++)
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);
778 for (int xx=0; xx<=4; xx++)
779 for (int yy=0; yy<=4; yy++)
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);
790 for (int xx=0; xx<=4; xx++)
791 for (int yy=0; yy<=4; yy++)
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);
802 for (int xx=0; xx<=4; xx++)
803 for (int yy=0; yy<=4; yy++)
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);
815 trace.beginBlock("checking operators");
816 trace.info() << calculus_3d_manual << endl;
817 trace.info() << cells_3d_manual.size() << endl;
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());
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;
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()
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()
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()
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()
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()
858 FATAL_ERROR( equal(primal_laplace_2d, primal_laplace_3d) );
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;
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()
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()
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()
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()
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()
887 FATAL_ERROR( equal(dual_laplace_2d, dual_laplace_3d) );
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());
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()
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()
905 trace.beginBlock("checking border");
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;
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");
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;
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});
936 trace.beginBlock("checking area");
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 );
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 );
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)
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.