DGtal 1.3.0
Loading...
Searching...
No Matches
Functions
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"

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 281 of file testLinearStructure.cpp.

282{
283 trace.info() << name << endl << op << endl << Eigen::MatrixXd(op.myContainer) << endl;
284}
DGtal::Trace::info
std::ostream & info()
DGtal::trace
Trace trace
Definition: Common.h:154

References DGtal::Trace::info(), and DGtal::trace.

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

◆ main()

int main ( void  )

Definition at line 1016 of file testLinearStructure.cpp.

1017{
1018#if !defined(TEST_HARDCODED_ORDER)
1019 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1020#endif
1021 test_manual_operators_3d();
1022 test_manual_operators_2d();
1023 test_linear_ring();
1024 test_linear_structure();
1025#if !defined(TEST_HARDCODED_ORDER)
1026 trace.warning() << "hardcoded order tests are NOT performed" << endl;
1027#endif
1028 return 0;
1029}
DGtal::Trace::warning
std::ostream & warning()
test_linear_structure
void test_linear_structure()
Definition: testLinearStructure.cpp:46
test_linear_ring
void test_linear_ring()
Definition: testLinearStructure.cpp:286
test_manual_operators_2d
void test_manual_operators_2d()
Definition: testLinearStructure.cpp:450
test_manual_operators_3d
void test_manual_operators_3d()
Definition: testLinearStructure.cpp:337

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

◆ test_linear_ring()

void test_linear_ring ( )

Definition at line 286 of file testLinearStructure.cpp.

287{
288 trace.beginBlock("linear ring");
289
290 const Domain domain(Point(-5,-5), Point(5,5));
291
292 typedef DiscreteExteriorCalculus<1, 2, EigenLinearAlgebraBackend> Calculus;
293 Calculus calculus;
294 calculus.initKSpace<Domain>(domain);
295
296 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(-8,kk), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
297 for (int kk=-8; kk<10; kk++) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,10), kk%2 == 0 ? Calculus::KSpace::POS : Calculus::KSpace::NEG) );
298 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(10,kk)) );
299 for (int kk=10; kk>-8; kk--) calculus.insertSCell( calculus.myKSpace.sCell(Point(kk,-8)) );
300 calculus.updateIndexes();
301
302 {
303 trace.info() << calculus << endl;
304 Board2D board;
305 board << domain;
306 board << calculus;
307 board.saveSVG("ring_structure.svg");
308 }
309
310 const Calculus::PrimalDerivative0 d0 = calculus.derivative<0, PRIMAL>();
311 display_operator_info("d0", d0);
312
313 const Calculus::PrimalHodge0 h0 = calculus.hodge<0, PRIMAL>();
314 display_operator_info("h0", h0);
315
316 const Calculus::DualDerivative0 d0p = calculus.derivative<0, DUAL>();
317 display_operator_info("d0p", d0p);
318
319 const Calculus::PrimalHodge1 h1 = calculus.hodge<1, PRIMAL>();
320 display_operator_info("h1", h1);
321
322 const Calculus::PrimalIdentity0 laplace = calculus.laplace<PRIMAL>();
323 display_operator_info("laplace", laplace);
324
325 const int laplace_size = calculus.kFormLength(0, PRIMAL);
326 const Eigen::MatrixXd laplace_dense(laplace.myContainer);
327
328 for (int ii=0; ii<laplace_size; ii++)
329 FATAL_ERROR( laplace_dense(ii,ii) == 2 );
330
331 FATAL_ERROR( laplace_dense.array().rowwise().sum().abs().sum() == 0 );
332 FATAL_ERROR( laplace_dense.transpose() == laplace_dense );
333
334 trace.endBlock();
335}
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
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::endBlock
double endBlock()
LibBoard::Board::saveSVG
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
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:281
domain
Domain domain
Definition: testProjection.cpp:88

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

References DGtal::Trace::beginBlock(), DGtal::DiscreteExteriorCalculusSolver< TCalculus, TLinearAlgebraSolver, order_in, duality_in, order_out, duality_out >::compute(), domain, DGtal::Trace::endBlock(), DGtal::Trace::info(), DGtal::KForm< TCalculus, order, duality >::myContainer, DGtal::PRIMAL, LibBoard::Board::saveSVG(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::sharp(), 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 450 of file testLinearStructure.cpp.

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

References DGtal::Trace::beginBlock(), dim, display_operator_info(), domain, DGtal::DUAL, DGtal::Trace::endBlock(), DGtal::Trace::info(), 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 337 of file testLinearStructure.cpp.

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

References DGtal::Trace::beginBlock(), display_operator_info(), domain, DGtal::DUAL, DGtal::Trace::endBlock(), DGtal::Trace::info(), DGtal::PRIMAL, and DGtal::trace.

Referenced by main().

Generated on Sun Nov 27 2022 15:15:59 for DGtal by doxygen 1.9.5