File failed to load: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/config/TeX-MML-AM_CHTML/MathJax.js
DGtal 2.0.0
testLinearStructure.cpp File Reference
DGtal Test Files
#include "DGtal/math/linalg/EigenSupport.h"
#include "DGtal/dec/DiscreteExteriorCalculus.h"
#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
#include "DGtal/io/boards/Board2D.h"
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "boost/version.hpp"
Include dependency graph for testLinearStructure.cpp:

Go to the source code of this file.

Functions

void test_linear_structure ()
template<typename Operator>
void display_operator_info (const std::string &name, const Operator &op)
void test_linear_ring ()
void test_manual_operators_3d ()
void test_manual_operators_2d ()
int main ()

Detailed Description

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author
Pierre Gueth (pierr.nosp@m.e.gu.nosp@m.eth@l.nosp@m.iris.nosp@m..cnrs.nosp@m..fr ) Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
David Coeurjolly (david.nosp@m..coe.nosp@m.urjol.nosp@m.ly@l.nosp@m.iris..nosp@m.cnrs.nosp@m..fr )

This file is part of the DGtal library

Definition in file testLinearStructure.cpp.

Function Documentation

◆ display_operator_info()

template<typename Operator>
void display_operator_info ( const std::string & name,
const Operator & op )

Definition at line 280 of file testLinearStructure.cpp.

281{
282 trace.info() << name << endl << op << endl << Eigen::MatrixXd(op.myContainer) << endl;
283}
DGtal::trace
Trace trace

References DGtal::trace.

Referenced by test_linear_ring(), test_manual_operators_2d(), and test_manual_operators_3d().

◆ main()

int main ( void )

Definition at line 1009 of file testLinearStructure.cpp.

1010{
1011#if !defined(TEST_HARDCODED_ORDER)
1012 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1013#endif
1014 test_manual_operators_3d();
1015 test_manual_operators_2d();
1016 test_linear_ring();
1017 test_linear_structure();
1018#if !defined(TEST_HARDCODED_ORDER)
1019 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1020#endif
1021 return 0;
1022}
test_linear_structure
void test_linear_structure()
Definition testLinearStructure.cpp:46
test_linear_ring
void test_linear_ring()
Definition testLinearStructure.cpp:285
test_manual_operators_2d
void test_manual_operators_2d()
Definition testLinearStructure.cpp:443
test_manual_operators_3d
void test_manual_operators_3d()
Definition testLinearStructure.cpp:336

References test_linear_ring(), test_linear_structure(), test_manual_operators_2d(), test_manual_operators_3d(), and DGtal::trace.

◆ test_linear_ring()

void test_linear_ring ( )

Definition at line 285 of file testLinearStructure.cpp.

286{
287 trace.beginBlock("linear ring");
288
289 const Domain domain(Point(-5,-5), Point(5,5));
290
291 typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
292 Calculus calculus;
293 calculus.initKSpace<Domain>(domain);
294
295 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(-8,kk), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
296 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,10), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
297 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(10,kk)) );
298 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,-8)) );
299 calculus.updateIndexes();
300
301 {
302 trace.info() << calculus << endl;
303 Board2D board;
304 board << domain;
305 board << calculus;
306 board.saveSVG("ring_structure.svg");
307 }
308
309 const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
310 display_operator_info("d0", d0);
311
312 const Calculus::PrimalHodge0 h0 = calculus.hodge<0, PRIMAL>();
313 display_operator_info("h0", h0);
314
315 const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
316 display_operator_info("d0p", d0p);
317
318 const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
319 display_operator_info("h1", h1);
320
321 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
322 display_operator_info("laplace", laplace);
323
324 const auto laplace_size = calculus.kFormLength(0, PRIMAL);
325 const Eigen::MatrixXd laplace_dense(laplace.myContainer);
326
327 for (int ii=0; ii<laplace_size; ii++)
328 FATAL_ERROR( laplace_dense(ii,ii) == 2 );
329
330 FATAL_ERROR( laplace_dense.array().rowwise().sum().abs().sum() == 0 );
331 FATAL_ERROR( laplace_dense.transpose() == laplace_dense );
332
333 trace.endBlock();
334}
DGtal::Board2D
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition Board2D.h:71
DGtal::DiscreteExteriorCalculus
Aim: DiscreteExteriorCalculus represents a calculus in the dec package. This is the main structure in...
Definition DiscreteExteriorCalculus.h:98
LibBoard::Board::saveSVG
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition Board.cpp:1011
calculus
PolyCalculus * calculus
Definition dgtalCalculus-geodesic.cpp:70
DGtal::Z2i::Point
Space::Point Point
Definition StdDefs.h:95
DGtal::PRIMAL
@ PRIMAL
Definition Duality.h:61
DGtal::DUAL
@ DUAL
Definition Duality.h:62
display_operator_info
void display_operator_info(const std::string &name, const Operator &op)
Definition testLinearStructure.cpp:280
domain
Domain domain
Definition testProjection.cpp:88

References calculus, display_operator_info(), domain, DGtal::DUAL, DGtal::PRIMAL, LibBoard::Board::saveSVG(), and DGtal::trace.

Referenced by main().

◆ test_linear_structure()

void test_linear_structure ( )

[neumann-creation]

[neumann-creation]

[input-dirac]

[input-dirac]

[neumann-laplace-definition]

[neumann-laplace-definition]

[neumann-solve]

[neumann-solve]

[dirichlet-creation]

[dirichlet-creation]

[dirichlet-laplace-definition]

[dirichlet-laplace-definition]

[dirichlet-solve]

[dirichlet-solve]

Definition at line 46 of file testLinearStructure.cpp.

47{
48 trace.beginBlock("creating dec problem with neumann border condition");
49
51 typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
52
53 SCells cells;
54
55 // fill cells container
56 {
57 const Calculus::KSpace kspace;
58
59 for (int kk=20; kk>0; kk--) cells.push_back( kspace.sCell(Point(0,kk), kk%2 != 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS) );
60 for (int kk=0; kk<10; kk++) cells.push_back( kspace.sCell(Point(kk,0)) );
61 for (int kk=0; kk<10; kk++) cells.push_back( kspace.sCell(Point(10,kk)) );
62 cells.push_back( kspace.sCell(Point(10,10)) );
63 cells.push_back( kspace.sCell(Point(9,10), Calculus::KSpace::NEG) );
64 for (int kk=10; kk<20; kk++) cells.push_back( kspace.sCell(Point(8,kk)) );
65 for (int kk=8; kk<12; kk++) cells.push_back( kspace.sCell(Point(kk,20)) );
66 for (int kk=20; kk>0; kk--) cells.push_back( kspace.sCell(Point(12,kk), kk%2 != 0 ? Calculus::KSpace::NEG : Calculus::KSpace::POS) );
67 cells.push_back( kspace.sCell(Point(12,0)) );
68 }
69
70 // fill calculus
71 const Domain domain(Point(-1,-1), Point(10,10));
72
73 // create DEC structure
74 Calculus calculus;
75 calculus.initKSpace<Domain>(domain);
76 for (SCells::const_iterator ci=cells.begin(), ce=cells.end(); ci!=ce; ci++) calculus.insertSCell( *ci );
77 calculus.updateIndexes();
79
80 trace.info() << calculus << endl;
81
83 Calculus::Cell dirac_cell = calculus.myKSpace.uCell(Point(10,4));
84 Calculus::PrimalForm0 dirac = Calculus::PrimalForm0::dirac(calculus, dirac_cell);
86 trace.info() << "dirac_cell_index=" << calculus.getCellIndex(dirac_cell) << endl;
87
88 {
89 Board2D board;
90 board << domain;
91 board << calculus;
92 board << dirac;
93 board.saveSVG("linear_structure_neumann_dirac.svg");
94 }
95
96 trace.endBlock();
97
98 {
99 trace.beginBlock("solving problem with neumann border condition using sparse qr solver");
100
102 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
104 trace.info() << "laplace = " << laplace << endl;
105
106 const Calculus::PrimalIdentity0 reorder = calculus.reorder<0, PRIMAL>(cells.begin(), cells.end());
107 const Calculus::PrimalIdentity0 laplace_ordered = reorder * laplace * reorder.transpose();
108 trace.info() << Eigen::MatrixXd(laplace_ordered.myContainer) << endl;
109
111 typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
112 typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;
113
114 Solver solver;
115 solver.compute(laplace);
116 Calculus::PrimalForm0 solved_solution = solver.solve(dirac);
118 Calculus::PrimalForm0 solved_solution_ordered = reorder * solved_solution;
119
120 Calculus::PrimalForm0 analytic_solution(calculus);
121 {
122 const Calculus::Index dirac_position = 17;
123 const Calculus::Index length = analytic_solution.length();
124 for (Calculus::Index kk=0; kk<length; kk++)
125 {
126 Calculus::Scalar alpha = 1. * (kk)/dirac_position * (kk+1.)/(dirac_position+1.);
127 if (kk>dirac_position)
128 {
129 alpha = 1. * (length-kk)/dirac_position * (length-kk-1.)/(dirac_position+1);
130 alpha -= 1. * (length-dirac_position)/dirac_position * (length-dirac_position-1.)/(dirac_position+1);
131 alpha += 1;
132 }
133 analytic_solution.myContainer(kk) = alpha;
134 }
135 }
136
137 const double shift = solved_solution_ordered.myContainer[0]-analytic_solution.myContainer[0];
138 for (Calculus::Index index=0; index<solved_solution_ordered.length(); index++) solved_solution_ordered.myContainer[index] -= shift;
139 solved_solution_ordered.myContainer /= solved_solution_ordered.myContainer.maxCoeff();
140
141 trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
142
143 {
144 std::ofstream handle("linear_structure_neumann.dat");
145 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
146 {
147 handle << solved_solution_ordered.myContainer(kk) << " " << analytic_solution.myContainer(kk) << endl;
148 }
149 }
150
151 const double error = (solved_solution_ordered-analytic_solution).myContainer.array().abs().maxCoeff();
152 trace.info() << "error=" << error << endl;
153 FATAL_ERROR( error < 1e-5 );
154
155 {
156 Board2D board;
157 board << domain;
158 board << calculus;
159 board << solved_solution;
160 board.saveSVG("linear_structure_neumann_solution.svg");
161 }
162
163 {
164 Calculus::PrimalForm1 solved_solution_gradient = calculus.derivative<0, PRIMAL>() * solved_solution;
165 Board2D board;
166 board << domain;
167 board << calculus;
168 board << solved_solution_gradient;
169 board << CustomStyle("VectorField", new VectorFieldStyle2D(1));
170 board << calculus.sharp(solved_solution_gradient);
171 board.saveSVG("linear_structure_neumann_solution_gradient.svg");
172 }
173
174 trace.endBlock();
175 }
176
177 trace.beginBlock("creating dec problem with dirichlet border condition");
178
180 calculus.insertSCell( calculus.myKSpace.sCell(Point(13,0)) );
181 calculus.insertSCell( calculus.myKSpace.sCell(Point(1,20), Calculus::KSpace::NEG) );
182 calculus.updateIndexes();
184
185 dirac = Calculus::PrimalForm0::dirac(calculus, dirac_cell);
186 trace.info() << calculus << endl;
187 trace.info() << "dirac_cell_index=" << calculus.getCellIndex(dirac_cell) << endl;
188
189 {
190 Board2D board;
191 board << domain;
192 board << calculus;
193 board << dirac;
194 board.saveSVG("linear_structure_dirichlet_dirac.svg");
195 }
196
197 trace.endBlock();
198
199 {
200 trace.beginBlock("solving problem with dirichlet border condition using sparse qr solver");
201
203 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
205 trace.info() << "laplace = " << laplace << endl;
206
207 const Calculus::PrimalIdentity0 reorder = calculus.reorder<0, PRIMAL>(cells.begin(), cells.end());
208 const Calculus::PrimalIdentity0 laplace_ordered = reorder * laplace * reorder.transpose();
209 trace.info() << Eigen::MatrixXd(laplace_ordered.myContainer) << endl;
210
212 typedef EigenLinearAlgebraBackend::SolverSparseQR LinearAlgebraSolver;
213 typedef DiscreteExteriorCalculusSolver<Calculus, LinearAlgebraSolver, 0, PRIMAL, 0, PRIMAL> Solver;
214
215 Solver solver;
216 solver.compute(laplace);
217 Calculus::PrimalForm0 solved_solution = solver.solve(dirac);
219 solved_solution.myContainer.array() /= solved_solution.myContainer.maxCoeff();
220 Calculus::PrimalForm0 solved_solution_ordered = reorder * solved_solution;
221
222 Calculus::PrimalForm0 analytic_solution(calculus);
223 {
224 const Calculus::Index dirac_position = 17;
225 const Calculus::Index length = analytic_solution.length();
226 for (Calculus::Index kk=0; kk<length; kk++)
227 {
228 Calculus::Scalar alpha = (kk+1.)/(dirac_position+1.);
229 if (kk>dirac_position)
230 {
231 alpha = 1. - (kk-dirac_position)/(1.*length-dirac_position);
232 }
233 analytic_solution.myContainer(kk) = alpha;
234 }
235 }
236
237 const double shift = solved_solution_ordered.myContainer[0]-analytic_solution.myContainer[0];
238 for (Calculus::Index index=0; index<solved_solution_ordered.length(); index++) solved_solution_ordered.myContainer[index] -= shift;
239 solved_solution_ordered.myContainer /= solved_solution_ordered.myContainer.maxCoeff();
240
241 trace.info() << solver.isValid() << " " << solver.myLinearAlgebraSolver.info() << endl;
242
243 {
244 std::ofstream handle("linear_structure_dirichlet.dat");
245 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
246 {
247 handle << solved_solution_ordered.myContainer(kk) << " " << analytic_solution.myContainer(kk) << endl;
248 }
249 }
250
251 const double error = (solved_solution_ordered-analytic_solution).myContainer.array().abs().maxCoeff();
252 trace.info() << "error=" << error << endl;
253 FATAL_ERROR( error < 1e-5 );
254
255 {
256 Board2D board;
257 board << domain;
258 board << calculus;
259 board << solved_solution;
260 board.saveSVG("linear_structure_dirichlet_solution.svg");
261 }
262
263 {
264 Calculus::PrimalForm1 solved_solution_gradient = calculus.derivative<0, PRIMAL>() * solved_solution;
265
266 Board2D board;
267 board << domain;
268 board << calculus;
269 board << solved_solution_gradient;
270 board << calculus.sharp(solved_solution_gradient);
271 board.saveSVG("linear_structure_dirichlet_solution_gradient.svg");
272 }
273
274 trace.endBlock();
275 }
276
277}
DGtal::DiscreteExteriorCalculusSolver
Aim: This wraps a linear algebra solver around a discrete exterior calculus.
Definition DiscreteExteriorCalculusSolver.h:70
DGtal::DiscreteExteriorCalculusSolver::solve
SolutionKForm solve(const InputKForm &input_kform) const
DGtal::DiscreteExteriorCalculusSolver::compute
DiscreteExteriorCalculusSolver & compute(const Operator &linear_operator)
DGtal::KForm::myContainer
Container myContainer
Definition KForm.h:131
DGtal::Z2i::SCells
KSpace::SCells SCells
Definition StdDefs.h:82
LibBoard::error
MessageStream error
DGtal::EigenLinearAlgebraBackend::SolverSparseQR
Eigen::SparseQR< SparseMatrix, Eigen::COLAMDOrdering< SparseMatrix::Index > > SolverSparseQR
Definition EigenSupport.h:110
index
unsigned int index(DGtal::uint32_t n, unsigned int b)
Definition testBits.cpp:44

References calculus, DGtal::DiscreteExteriorCalculusSolver< TCalculus, TLinearAlgebraSolver, order_in, duality_in, order_out, duality_out >::compute(), domain, index(), DGtal::KForm< TCalculus, order, duality >::myContainer, DGtal::PRIMAL, LibBoard::Board::saveSVG(), DGtal::DiscreteExteriorCalculusSolver< TCalculus, TLinearAlgebraSolver, order_in, duality_in, order_out, duality_out >::solve(), and DGtal::trace.

Referenced by main().

◆ test_manual_operators_2d()

void test_manual_operators_2d ( )

Definition at line 443 of file testLinearStructure.cpp.

444{
445 trace.beginBlock("testing 2d operators");
446
447 const Domain domain(Point(0,0), Point(5,4));
448
449 typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;
450
451 Calculus primal_calculus;
452 primal_calculus.initKSpace<Domain>(domain);
453
454 { // filling primal calculus
455 // 0-cells
456 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,2)) );
457 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,2)) );
458 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,4)) );
459 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,4)) );
460 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,6)) );
461 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,6)) );
462
463 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(1,2)) ); // insert cell
464
465 // 1-cells
466 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,2), Calculus::KSpace::NEG) ); // insert negative cell
467 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,2), Calculus::KSpace::POS) ); // then reinserting negative cell in structure
468 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,3), Calculus::KSpace::NEG) );
469 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,3), Calculus::KSpace::POS) );
470 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,4), Calculus::KSpace::NEG) );
471 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(2,5), Calculus::KSpace::POS) );
472 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(4,5), Calculus::KSpace::NEG) );
473 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,6), Calculus::KSpace::POS) );
474
475 primal_calculus.eraseCell( primal_calculus.myKSpace.uCell(Point(1,2)) ); // then remove it
476
477 // 2-cells
478 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,3)) );
479 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(Point(3,5)) );
480
481 primal_calculus.updateIndexes();
482
483 trace.beginBlock("primal calculus");
484 trace.info() << primal_calculus << endl;
485 for (Calculus::ConstIterator iter_property=primal_calculus.begin(), iter_property_end=primal_calculus.end(); iter_property!=iter_property_end; iter_property++)
486 {
487 const Calculus::Cell cell = iter_property->first;
488 const Calculus::Property property = iter_property->second;
489 const Dimension dim = primal_calculus.myKSpace.uDim(cell);
490 const Calculus::SCell signed_cell = primal_calculus.myKSpace.signs(cell, property.flipped ? Calculus::KSpace::NEG : Calculus::KSpace::POS);
491
492 ASSERT( signed_cell == primal_calculus.getSCell(dim, PRIMAL, property.index) );
493
494 trace.info() << cell
495 << " " << dim
496 << " " << signed_cell
497 << " " << property.primal_size
498 << " " << property.dual_size
499 << " " << property.index
500 << " " << (property.flipped ? "negative" : "positive")
501 << endl;
502 }
503 trace.endBlock();
504 }
505
506
507 Calculus dual_calculus;
508 dual_calculus.initKSpace<Domain>(domain);
509
510 { // filling dual calculus
511 // 2-cells
512 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,3)) );
513 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,3)) );
514 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,5)) );
515 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,5)) );
516 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,7)) );
517 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,7)) );
518
519 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(6,3)) ); // insert cell
520
521 // 1-cells
522 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,3), Calculus::KSpace::NEG) ); // insert negative cell
523 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,3), Calculus::KSpace::POS) ); // then reinserting negative cell in structure
524 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,4), Calculus::KSpace::POS) );
525 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,4), Calculus::KSpace::NEG) );
526 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,5), Calculus::KSpace::NEG) );
527 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(7,6), Calculus::KSpace::NEG) );
528 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(9,6), Calculus::KSpace::POS) );
529 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,7), Calculus::KSpace::POS) );
530
531 dual_calculus.eraseCell( dual_calculus.myKSpace.uCell(Point(6,3)) ); // then remove it
532
533 // 0-cells
534 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,4)) );
535 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(Point(8,6)) );
536
537 dual_calculus.updateIndexes();
538
539 trace.beginBlock("dual calculus");
540 trace.info() << dual_calculus << endl;
541 for (Calculus::ConstIterator iter_property=dual_calculus.begin(), iter_property_end=dual_calculus.end(); iter_property!=iter_property_end; iter_property++)
542 {
543 const Calculus::Cell cell = iter_property->first;
544 const Calculus::Property property = iter_property->second;
545 const Dimension dim = dual_calculus.myKSpace.uDim(cell);
546 const Calculus::SCell signed_cell = dual_calculus.myKSpace.signs(cell, property.flipped ? Calculus::KSpace::NEG : Calculus::KSpace::POS);
547
548 ASSERT( signed_cell == dual_calculus.getSCell(dim, PRIMAL, property.index) );
549
550 trace.info() << cell
551 << " " << dim
552 << " " << signed_cell
553 << " " << property.primal_size
554 << " " << property.dual_size
555 << " " << property.index
556 << " " << (property.flipped ? "negative" : "positive")
557 << endl;
558 }
559 trace.endBlock();
560 }
561
562 {
563 Board2D board;
564 board << domain;
565 board << primal_calculus;
566 board << dual_calculus;
567 board.saveSVG("operators_structure.svg");
568 }
569
570 trace.beginBlock("base operators");
571
572 const Calculus::PrimalDerivative0 primal_d0 = primal_calculus.derivative<0, PRIMAL>();
573 const Calculus::DualDerivative0 dual_d0p = dual_calculus.derivative<0, DUAL>();
574 {
575 display_operator_info("primal d0", primal_d0);
576 display_operator_info("dual d0p", dual_d0p);
577
578#if defined(TEST_HARDCODED_ORDER)
579 Eigen::MatrixXd d0_th(7, 6);
580 d0_th <<
581 -1, 1, 0, 0, 0, 0,
582 1, 0, 0, -1, 0, 0,
583 0, 0, 1, 0, 0, -1,
584 0, 0, 0, 0, -1, 1,
585 0, -1, 1, 0, 0, 0,
586 0, 0, -1, 1, 0, 0,
587 0, 0, 0, -1, 1, 0;
588 FATAL_ERROR( Eigen::MatrixXd(primal_d0.myContainer) == d0_th );
589
590 Eigen::MatrixXd d0p_th(7, 6);
591 d0p_th <<
592 1, -1, 0, 0, 0, 0,
593 0, 0, -1, 1, 0, 0,
594 0, 1, 0, -1, 0, 0,
595 0, 0, 1, 0, -1, 0,
596 0, 0, 0, 0, 1, -1,
597 -1, 0, 1, 0, 0, 0,
598 0, 0, 0, -1, 0, 1;
599 FATAL_ERROR( Eigen::MatrixXd(dual_d0p.myContainer) == d0p_th );
600#endif
601 }
602
603 const Calculus::PrimalDerivative1 primal_d1 = primal_calculus.derivative<1, PRIMAL>();
604 const Calculus::DualDerivative1 dual_d1p = dual_calculus.derivative<1, DUAL>();
605 {
606 display_operator_info("primal d1", primal_d1);
607 display_operator_info("dual d1p", dual_d1p);
608
609#if defined(TEST_HARDCODED_ORDER)
610 Eigen::MatrixXd d1_th(2, 7);
611 d1_th <<
612 1, 1, 0, 0, 1, 1, 0,
613 0, 0, -1, -1, 0, -1, -1;
614 FATAL_ERROR( Eigen::MatrixXd(primal_d1.myContainer) == d1_th );
615
616 Eigen::MatrixXd d1p_th(2, 7);
617 d1p_th <<
618 -1, -1, -1, 0, 0, -1, 0,
619 0, 1, 0, 1, 1, 0, 1;
620 FATAL_ERROR( Eigen::MatrixXd(dual_d1p.myContainer) == d1p_th );
621#endif
622 }
623
624 {
625 display_operator_info("primal d1*d0", primal_d1*primal_d0);
626 display_operator_info("dual d1p*d0p", dual_d1p*dual_d0p);
627
628 FATAL_ERROR( Eigen::MatrixXd((primal_d1*primal_d0).myContainer) == Eigen::MatrixXd::Zero(2,6) );
629 FATAL_ERROR( Eigen::MatrixXd((dual_d1p*dual_d0p).myContainer) == Eigen::MatrixXd::Zero(2,6) );
630 }
631
632 const Calculus::PrimalHodge0 primal_h0 = primal_calculus.hodge<0, PRIMAL>();
633 const Calculus::DualHodge0 dual_h0p = dual_calculus.hodge<0, DUAL>();
634 const Calculus::DualHodge2 primal_h2p = primal_calculus.hodge<2, DUAL>();
635 const Calculus::PrimalHodge2 dual_h2 = dual_calculus.hodge<2, PRIMAL>();
636 {
637 display_operator_info("primal h0", primal_h0);
638 display_operator_info("dual h0p", dual_h0p);
639 display_operator_info("primal h2p", primal_h2p);
640 display_operator_info("dual h2", dual_h2);
641
642 FATAL_ERROR( Eigen::MatrixXd(primal_h0.myContainer) == Eigen::MatrixXd::Identity(6,6) );
643 FATAL_ERROR( Eigen::MatrixXd(dual_h0p.myContainer) == Eigen::MatrixXd::Identity(6,6) );
644 FATAL_ERROR( Eigen::MatrixXd(primal_h2p.myContainer) == Eigen::MatrixXd::Identity(6,6) );
645 FATAL_ERROR( Eigen::MatrixXd(dual_h2.myContainer) == Eigen::MatrixXd::Identity(6,6) );
646 }
647
648 const Calculus::PrimalHodge2 primal_h2 = primal_calculus.hodge<2, PRIMAL>();
649 const Calculus::DualHodge2 dual_h2p = dual_calculus.hodge<2, DUAL>();
650 const Calculus::DualHodge0 primal_h0p = primal_calculus.hodge<0, DUAL>();
651 const Calculus::PrimalHodge0 dual_h0 = dual_calculus.hodge<0, PRIMAL>();
652 {
653 display_operator_info("primal h2", primal_h2);
654 display_operator_info("dual h2p", dual_h2p);
655 display_operator_info("primal h0p", primal_h0p);
656 display_operator_info("dual h0", dual_h0);
657
658 FATAL_ERROR( Eigen::MatrixXd(primal_h2.myContainer) == Eigen::MatrixXd::Identity(2,2) );
659 FATAL_ERROR( Eigen::MatrixXd(dual_h2p.myContainer) == Eigen::MatrixXd::Identity(2,2) );
660 FATAL_ERROR( Eigen::MatrixXd(primal_h0p.myContainer) == Eigen::MatrixXd::Identity(2,2) );
661 FATAL_ERROR( Eigen::MatrixXd(dual_h0.myContainer) == Eigen::MatrixXd::Identity(2,2) );
662 }
663
664 const Calculus::DualDerivative0 primal_d0p = primal_calculus.derivative<0, DUAL>();
665 const Calculus::PrimalDerivative0 dual_d0 = dual_calculus.derivative<0, PRIMAL>();
666 {
667 display_operator_info("primal d0p", primal_d0p);
668 display_operator_info("dual d0", dual_d0);
669
670#if defined(TEST_HARDCODED_ORDER)
671 Eigen::MatrixXd d0p_th_transpose(2, 7);
672 d0p_th_transpose <<
673 1, 1, 0, 0, 1, 1, 0,
674 0, 0, -1, -1, 0, -1, -1;
675 FATAL_ERROR( Eigen::MatrixXd(primal_d0p.myContainer) == d0p_th_transpose.transpose() );
676
677 Eigen::MatrixXd minus_d0_th_transpose(2, 7);
678 minus_d0_th_transpose <<
679 -1, -1, -1, 0, 0, -1, 0,
680 0, 1, 0, 1, 1, 0, 1;
681 FATAL_ERROR( Eigen::MatrixXd(dual_d0.myContainer) == -minus_d0_th_transpose.transpose() );
682#endif
683 }
684
685 const Calculus::DualDerivative1 primal_d1p = primal_calculus.derivative<1, DUAL>();
686 const Calculus::PrimalDerivative1 dual_d1 = dual_calculus.derivative<1, PRIMAL>();
687 {
688 display_operator_info("primal d1p", primal_d1p);
689 display_operator_info("dual d1", dual_d1);
690
691#if defined(TEST_HARDCODED_ORDER)
692 Eigen::MatrixXd minus_d1p_th_transpose(7, 6);
693 minus_d1p_th_transpose <<
694 -1, 1, 0, 0, 0, 0,
695 1, 0, 0, -1, 0, 0,
696 0, 0, 1, 0, 0, -1,
697 0, 0, 0, 0, -1, 1,
698 0, -1, 1, 0, 0, 0,
699 0, 0, -1, 1, 0, 0,
700 0, 0, 0, -1, 1, 0;
701 FATAL_ERROR( Eigen::MatrixXd(primal_d1p.myContainer) == -minus_d1p_th_transpose.transpose() );
702
703 Eigen::MatrixXd d1_th_transpose(7, 6);
704 d1_th_transpose <<
705 1, -1, 0, 0, 0, 0,
706 0, 0, -1, 1, 0, 0,
707 0, 1, 0, -1, 0, 0,
708 0, 0, 1, 0, -1, 0,
709 0, 0, 0, 0, 1, -1,
710 -1, 0, 1, 0, 0, 0,
711 0, 0, 0, -1, 0, 1;
712 FATAL_ERROR( Eigen::MatrixXd(dual_d1.myContainer) == d1_th_transpose.transpose() );
713#endif
714 }
715
716 const Calculus::PrimalHodge1 primal_h1 = primal_calculus.hodge<1, PRIMAL>();
717 const Calculus::DualHodge1 dual_h1p = dual_calculus.hodge<1, DUAL>();
718 const Calculus::DualHodge1 primal_h1p = primal_calculus.hodge<1, DUAL>();
719 const Calculus::PrimalHodge1 dual_h1 = dual_calculus.hodge<1, PRIMAL>();
720 {
721 display_operator_info("primal h1", primal_h1);
722 display_operator_info("dual h1p", dual_h1p);
723 display_operator_info("primal h1p", primal_h1p);
724 display_operator_info("dual h1", dual_h1);
725
726 FATAL_ERROR( Eigen::MatrixXd(primal_h1.myContainer) == Eigen::MatrixXd::Identity(7,7) );
727 FATAL_ERROR( Eigen::MatrixXd(dual_h1p.myContainer) == -Eigen::MatrixXd::Identity(7,7) );
728 FATAL_ERROR( Eigen::MatrixXd((primal_h1p*primal_h1).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
729 FATAL_ERROR( Eigen::MatrixXd((dual_h1*dual_h1p).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
730 FATAL_ERROR( Eigen::MatrixXd((primal_h1*primal_h1p).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
731 FATAL_ERROR( Eigen::MatrixXd((dual_h1p*dual_h1).myContainer) == -Eigen::MatrixXd::Identity(7,7) );
732 }
733
734 trace.endBlock();
735
736 trace.beginBlock("laplace operators");
737 display_operator_info("primal laplace", primal_calculus.laplace<PRIMAL>());
738 display_operator_info("dual laplacep", dual_calculus.laplace<DUAL>());
739 display_operator_info("primal laplacep", primal_calculus.laplace<DUAL>());
740 display_operator_info("dual laplace", dual_calculus.laplace<PRIMAL>());
741 trace.endBlock();
742
743 trace.beginBlock("sharp operators");
744
745 { // primal sharp
746 display_operator_info("primal #x", primal_calculus.sharpDirectional<PRIMAL>(0));
747 display_operator_info("dual #xp", dual_calculus.sharpDirectional<DUAL>(0));
748 display_operator_info("primal #y", primal_calculus.sharpDirectional<PRIMAL>(1));
749 display_operator_info("dual #yp", dual_calculus.sharpDirectional<DUAL>(1));
750
751 {
752 Calculus::PrimalForm1::Container dx_container(7);
753 dx_container << -1, 0, 0, -1, 0, 1, 0;
754 const Calculus::PrimalForm1 primal_dx(primal_calculus, dx_container);
755 const Calculus::PrimalVectorField primal_dx_field = primal_calculus.sharp(primal_dx);
756
757 Calculus::PrimalForm1::Container dxp_container(7);
758 dxp_container << 1, -1, 0, 0, 1, 0, 0;
759 const Calculus::DualForm1 dual_dx(dual_calculus, dxp_container);
760 const Calculus::DualVectorField dual_dx_field = dual_calculus.sharp(dual_dx);
761
762 {
763 Board2D board;
764 board << domain;
765 board << primal_calculus;
766 board << primal_dx << primal_dx_field;
767 board << dual_calculus;
768 board << dual_dx << dual_dx_field;
769 board.saveSVG("operators_sharp_dx_primal.svg");
770 }
771
772#if defined(TEST_HARDCODED_ORDER)
773 FATAL_ERROR( primal_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(6) );
774 FATAL_ERROR( primal_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(6) );
775 FATAL_ERROR( dual_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(6) );
776 FATAL_ERROR( dual_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(6) );
777#endif
778 }
779
780 {
781 Calculus::PrimalForm1::Container dy_container(7);
782 dy_container << 0, 1, 1, 0, -1, 0, -1;
783 const Calculus::PrimalForm1 primal_dy(primal_calculus, dy_container);
784 const Calculus::PrimalVectorField primal_dy_field = primal_calculus.sharp(primal_dy);
785
786 Calculus::PrimalForm1::Container dyp_container(7);
787 dyp_container << 0, 0, 1, 1, 0, -1, -1;
788 const Calculus::DualForm1 dual_dy(dual_calculus, dyp_container);
789 const Calculus::DualVectorField dual_dy_field = dual_calculus.sharp(dual_dy);
790
791 {
792 Board2D board;
793 board << domain;
794 board << primal_calculus;
795 board << primal_dy << primal_dy_field;
796 board << dual_calculus;
797 board << dual_dy << dual_dy_field;
798 board.saveSVG("operators_sharp_dy_primal.svg");
799 }
800
801#if defined(TEST_HARDCODED_ORDER)
802 FATAL_ERROR( primal_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(6) );
803 FATAL_ERROR( primal_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(6) );
804 FATAL_ERROR( dual_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(6) );
805 FATAL_ERROR( dual_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(6) );
806#endif
807 }
808 }
809
810 { // dual sharp
811 display_operator_info("primal #xp", primal_calculus.sharpDirectional<DUAL>(0));
812 display_operator_info("dual #x", dual_calculus.sharpDirectional<PRIMAL>(0));
813 display_operator_info("primal #yp", primal_calculus.sharpDirectional<DUAL>(1));
814 display_operator_info("dual #y", dual_calculus.sharpDirectional<PRIMAL>(1));
815
816 {
817 Calculus::DualForm1::Container dx_container(7);
818 dx_container << 0, -1, -1, 0, 1, 0, 1;
819 const Calculus::DualForm1 primal_dx(primal_calculus, dx_container);
820 const Calculus::DualVectorField primal_dx_field = primal_calculus.sharp(primal_dx);
821
822 Calculus::DualForm1::Container dxp_container(7);
823 dxp_container << 0, 0, 1, 1, 0, -1, -1;
824 const Calculus::PrimalForm1 dual_dx(dual_calculus, dxp_container);
825 const Calculus::PrimalVectorField dual_dx_field = dual_calculus.sharp(dual_dx);
826
827 {
828 Board2D board;
829 board << domain;
830 board << primal_calculus;
831 board << primal_dx << primal_dx_field;
832 board << dual_calculus;
833 board << dual_dx << dual_dx_field;
834 board.saveSVG("operators_sharp_dx_dual.svg");
835 }
836
837#if defined(TEST_HARDCODED_ORDER)
838 FATAL_ERROR( primal_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(2) );
839 FATAL_ERROR( primal_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(2) );
840 FATAL_ERROR( dual_dx_field.myCoordinates.col(0) == Eigen::VectorXd::Ones(2) );
841 FATAL_ERROR( dual_dx_field.myCoordinates.col(1) == Eigen::VectorXd::Zero(2) );
842#endif
843 }
844
845 {
846 Calculus::DualForm1::Container dy_container(7);
847 dy_container << -1, 0, 0, -1, 0 , 1, 0;
848 const Calculus::DualForm1 primal_dy(primal_calculus, dy_container);
849 const Calculus::DualVectorField primal_dy_field = primal_calculus.sharp(primal_dy);
850
851 Calculus::DualForm1::Container dyp_container(7);
852 dyp_container << -1, 1, 0, 0, -1, 0, 0;
853 const Calculus::PrimalForm1 dual_dy(dual_calculus, dyp_container);
854 const Calculus::PrimalVectorField dual_dy_field = dual_calculus.sharp(dual_dy);
855
856 {
857 Board2D board;
858 board << domain;
859 board << primal_calculus;
860 board << primal_dy << primal_dy_field;
861 board << dual_calculus;
862 board << dual_dy << dual_dy_field;
863 board.saveSVG("operators_sharp_dy_dual.svg");
864 }
865
866#if defined(TEST_HARDCODED_ORDER)
867 FATAL_ERROR( primal_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(2) );
868 FATAL_ERROR( primal_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(2) );
869 FATAL_ERROR( dual_dy_field.myCoordinates.col(0) == Eigen::VectorXd::Zero(2) );
870 FATAL_ERROR( dual_dy_field.myCoordinates.col(1) == Eigen::VectorXd::Ones(2) );
871#endif
872 }
873 }
874
875 trace.endBlock();
876
877 trace.beginBlock("flat operators");
878
879 { // primal flat
880 display_operator_info("primal bx", primal_calculus.flatDirectional<PRIMAL>(0));
881 display_operator_info("dual bxp", dual_calculus.flatDirectional<DUAL>(0));
882 display_operator_info("primal by", primal_calculus.flatDirectional<PRIMAL>(1));
883 display_operator_info("dual byp", dual_calculus.flatDirectional<DUAL>(1));
884
885 Calculus::PrimalVectorField::Coordinates dx_coords(6,2);
886 dx_coords.col(0) = Eigen::VectorXd::Ones(6);
887 dx_coords.col(1) = Eigen::VectorXd::Zero(6);
888
889 Calculus::PrimalVectorField::Coordinates dy_coords(6,2);
890 dy_coords.col(0) = Eigen::VectorXd::Zero(6);
891 dy_coords.col(1) = Eigen::VectorXd::Ones(6);
892
893 const Calculus::PrimalVectorField primal_dx_field(primal_calculus, dx_coords);
894 const Calculus::PrimalForm1 primal_dx = primal_calculus.flat(primal_dx_field);
895 const Calculus::DualVectorField dual_dx_field(dual_calculus, dx_coords);
896 const Calculus::DualForm1 dual_dx = dual_calculus.flat(dual_dx_field);
897
898 {
899 Board2D board;
900 board << domain;
901 board << primal_calculus;
902 board << primal_dx << primal_dx_field;
903 board << dual_calculus;
904 board << dual_dx << dual_dx_field;
905 board.saveSVG("operators_flat_dx_primal.svg");
906 }
907
908 const Calculus::PrimalVectorField primal_dy_field(primal_calculus, dy_coords);
909 const Calculus::PrimalForm1 primal_dy = primal_calculus.flat(primal_dy_field);
910 const Calculus::DualVectorField dual_dy_field(dual_calculus, dy_coords);
911 const Calculus::DualForm1 dual_dy = dual_calculus.flat(dual_dy_field);
912
913 {
914 Board2D board;
915 board << domain;
916 board << primal_calculus;
917 board << primal_dy << primal_dy_field;
918 board << dual_calculus;
919 board << dual_dy << dual_dy_field;
920 board.saveSVG("operators_flat_dy_primal.svg");
921 }
922
923#if defined(TEST_HARDCODED_ORDER)
924 Calculus::PrimalForm1::Container dx_container(7);
925 dx_container << -1, 0, 0, -1, 0, 1, 0;
926 Calculus::PrimalForm1::Container dxp_container(7);
927 dxp_container << 1, -1, 0, 0, 1, 0, 0;
928 FATAL_ERROR( primal_dx.myContainer == dx_container );
929 FATAL_ERROR( dual_dx.myContainer == dxp_container );
930
931 Calculus::PrimalForm1::Container dy_container(7);
932 dy_container << 0, 1, 1, 0, -1, 0, -1;
933 Calculus::PrimalForm1::Container dyp_container(7);
934 dyp_container << 0, 0, 1, 1, 0, -1, -1;
935 FATAL_ERROR( primal_dy.myContainer == dy_container );
936 FATAL_ERROR( dual_dy.myContainer == dyp_container );
937#endif
938 }
939
940 { // dual flat
941 display_operator_info("primal bxp", primal_calculus.flatDirectional<DUAL>(0));
942 display_operator_info("dual bx", dual_calculus.flatDirectional<PRIMAL>(0));
943 display_operator_info("primal byp", primal_calculus.flatDirectional<DUAL>(1));
944 display_operator_info("dual by", dual_calculus.flatDirectional<PRIMAL>(1));
945
946 Calculus::PrimalVectorField::Coordinates dx_coords(2,2);
947 dx_coords.col(0) = Eigen::VectorXd::Ones(2);
948 dx_coords.col(1) = Eigen::VectorXd::Zero(2);
949
950 Calculus::PrimalVectorField::Coordinates dy_coords(2,2);
951 dy_coords.col(0) = Eigen::VectorXd::Zero(2);
952 dy_coords.col(1) = Eigen::VectorXd::Ones(2);
953
954 const Calculus::DualVectorField primal_dx_field(primal_calculus, dx_coords);
955 const Calculus::DualForm1 primal_dx = primal_calculus.flat(primal_dx_field);
956 const Calculus::PrimalVectorField dual_dx_field(dual_calculus, dx_coords);
957 const Calculus::PrimalForm1 dual_dx = dual_calculus.flat(dual_dx_field);
958
959 {
960 Board2D board;
961 board << domain;
962 board << primal_calculus;
963 board << primal_dx << primal_dx_field;
964 board << dual_calculus;
965 board << dual_dx << dual_dx_field;
966 board.saveSVG("operators_flat_dx_dual.svg");
967 }
968
969 const Calculus::DualVectorField primal_dy_field(primal_calculus, dy_coords);
970 const Calculus::DualForm1 primal_dy = primal_calculus.flat(primal_dy_field);
971 const Calculus::PrimalVectorField dual_dy_field(dual_calculus, dy_coords);
972 const Calculus::PrimalForm1 dual_dy = dual_calculus.flat(dual_dy_field);
973
974 {
975 Board2D board;
976 board << domain;
977 board << primal_calculus;
978 board << primal_dy << primal_dy_field;
979 board << dual_calculus;
980 board << dual_dy << dual_dy_field;
981 board.saveSVG("operators_flat_dy_dual.svg");
982 }
983
984#if defined(TEST_HARDCODED_ORDER)
985 Calculus::PrimalForm1::Container dx_container(7);
986 dx_container << 0, -1, -1, 0, 1, 0, 1;
987 Calculus::PrimalForm1::Container dxp_container(7);
988 dxp_container << 0, 0, 1, 1, 0, -1, -1;
989 FATAL_ERROR( primal_dx.myContainer == dx_container );
990 FATAL_ERROR( dual_dx.myContainer == dxp_container );
991
992 Calculus::PrimalForm1::Container dy_container(7);
993 dy_container << -1, 0, 0, -1, 0, 1, 0;
994 Calculus::PrimalForm1::Container dyp_container(7);
995 dyp_container << -1, 1, 0, 0, -1, 0, 0;
996 FATAL_ERROR( primal_dy.myContainer == dy_container );
997 FATAL_ERROR( dual_dy.myContainer == dyp_container );
998#endif
999 }
1000
1001 trace.endBlock();
1002
1003
1004 trace.endBlock();
1005}
DGtal::DiscreteExteriorCalculus::initKSpace
void initKSpace(ConstAlias< TDomain > domain)
DGtal::Dimension
DGtal::uint32_t Dimension
Definition Common.h:119

References dim, display_operator_info(), domain, DGtal::DUAL, DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger >::initKSpace(), DGtal::PRIMAL, LibBoard::Board::saveSVG(), and DGtal::trace.

Referenced by main().

◆ test_manual_operators_3d()

void test_manual_operators_3d ( )

Definition at line 336 of file testLinearStructure.cpp.

337{
338 trace.beginBlock("testing 3d operators");
339
340 const Z3i::Domain domain(Z3i::Point(0,0,0), Z3i::Point(1,1,1));
341
342 typedef DiscreteExteriorCalculus<3, 3, EigenLinearAlgebraBackend> Calculus;
343
344 Calculus calculus;
345 calculus.initKSpace<Z3i::Domain>(domain);
346
347 { // filling primal calculus
348 // 0-cells
349 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,0,0)) );
350 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,0,0)) );
351
352 // 1-cells
353 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,1,0), Calculus::KSpace::POS) );
354 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,0,1), Calculus::KSpace::POS) );
355 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,0,0), Calculus::KSpace::POS) );
356 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,1,0), Calculus::KSpace::NEG) );
357 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,0,1), Calculus::KSpace::NEG) );
358
359 // 2-cells
360 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(0,1,1), Calculus::KSpace::NEG) );
361 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,0,1), Calculus::KSpace::POS) ); //FIXME strange cell orientation
362 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(2,1,1), Calculus::KSpace::POS) );
363 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,1,0), Calculus::KSpace::NEG) );
364
365 // 3-cells
366 calculus.insertSCell( calculus.myKSpace.sCell(Z3i::Point(1,1,1)) );
367
368 calculus.updateIndexes();
369
370 trace.info() << calculus << endl;
371 }
372
373 trace.beginBlock("base operators");
374
375 {
376 const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
377 const Calculus::DualDerivative2 d2p = calculus.derivative<2, DUAL>();
378 display_operator_info("d0", d0);
379 display_operator_info("d2p", d2p);
380
381#if defined(TEST_HARDCODED_ORDER)
382 Eigen::MatrixXd d0_th(5, 2);
383 d0_th <<
384 -1, 0,
385 -1, 0,
386 -1, 1,
387 0, 1,
388 0, 1;
389
390 FATAL_ERROR( Eigen::MatrixXd(d0.myContainer) == d0_th );
391#endif
392 FATAL_ERROR( Eigen::MatrixXd(d2p.transpose().myContainer) == Eigen::MatrixXd(d0.myContainer) );
393 }
394
395 {
396 const Calculus::PrimalDerivative1 d1 = calculus.derivative<1, PRIMAL>();
397 const Calculus::DualDerivative1 d1p = calculus.derivative<1, DUAL>();
398 display_operator_info("d1", d1);
399 display_operator_info("d1p", d1p);
400
401#if defined(TEST_HARDCODED_ORDER)
402 Eigen::MatrixXd d1_th(4, 5);
403 d1_th <<
404 0, -1, 1, 0, -1,
405 1, 0, -1, 1, 0,
406 0, 0, 0, -1, 1,
407 -1, 1, 0, 0, 0;
408
409 FATAL_ERROR( Eigen::MatrixXd(d1.myContainer) == d1_th );
410#endif
411 FATAL_ERROR( Eigen::MatrixXd(d1p.transpose().myContainer) == Eigen::MatrixXd(d1.myContainer) );
412 }
413
414 {
415 const Calculus::PrimalDerivative2 d2 = calculus.derivative<2, PRIMAL>();
416 const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
417 display_operator_info("d2", d2);
418 display_operator_info("d0p", d0p);
419
420#if defined(TEST_HARDCODED_ORDER)
421 Eigen::MatrixXd d2_th(1, 4);
422 d2_th << -1, -1, -1, -1;
423
424 FATAL_ERROR( Eigen::MatrixXd(d2.myContainer) == d2_th );
425#endif
426 FATAL_ERROR( Eigen::MatrixXd(d0p.transpose().myContainer) == Eigen::MatrixXd(d2.myContainer) );
427 }
428
429 trace.endBlock();
430
431 /*
432 typedef PolyscopeViewer<Z3i::Space, Z3i::KSpace> Viewer;
433 Viewer* viewer = new Viewer(calculus.myKSpace);
434 (*viewer) << DGtal::Color(0,0,0);
435 (*viewer) << domain;
436 (*viewer) << calculus;
437 viewer->show();
438 */
439
440 trace.endBlock();
441}
DGtal::Z3i::Domain
HyperRectDomain< Space > Domain
Definition StdDefs.h:172
DGtal::Z3i::Point
Space::Point Point
Definition StdDefs.h:168

References calculus, display_operator_info(), domain, DGtal::DUAL, DGtal::PRIMAL, and DGtal::trace.

Referenced by main().