48 trace.beginBlock(
"creating dec problem with neumann border condition");
57 const Calculus::KSpace kspace;
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)) );
76 for (SCells::const_iterator ci=cells.begin(), ce=cells.end(); ci!=ce; ci++)
calculus.insertSCell( *ci );
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;
93 board.
saveSVG(
"linear_structure_neumann_dirac.svg");
99 trace.beginBlock(
"solving problem with neumann border condition using sparse qr solver");
102 const Calculus::PrimalIdentity0 laplace =
calculus.laplace<
PRIMAL>();
104 trace.info() <<
"laplace = " << laplace << endl;
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;
116 Calculus::PrimalForm0 solved_solution = solver.
solve(dirac);
118 Calculus::PrimalForm0 solved_solution_ordered = reorder * solved_solution;
120 Calculus::PrimalForm0 analytic_solution(
calculus);
122 const Calculus::Index dirac_position = 17;
123 const Calculus::Index length = analytic_solution.length();
124 for (Calculus::Index kk=0; kk<length; kk++)
126 Calculus::Scalar alpha = 1. * (kk)/dirac_position * (kk+1.)/(dirac_position+1.);
127 if (kk>dirac_position)
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);
133 analytic_solution.myContainer(kk) = alpha;
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();
141 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
144 std::ofstream handle(
"linear_structure_neumann.dat");
145 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
147 handle << solved_solution_ordered.myContainer(kk) <<
" " << analytic_solution.myContainer(kk) << endl;
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 );
159 board << solved_solution;
160 board.
saveSVG(
"linear_structure_neumann_solution.svg");
164 Calculus::PrimalForm1 solved_solution_gradient =
calculus.derivative<0,
PRIMAL>() * solved_solution;
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");
177 trace.beginBlock(
"creating dec problem with dirichlet border condition");
185 dirac = Calculus::PrimalForm0::dirac(
calculus, dirac_cell);
187 trace.info() <<
"dirac_cell_index=" <<
calculus.getCellIndex(dirac_cell) << endl;
194 board.
saveSVG(
"linear_structure_dirichlet_dirac.svg");
200 trace.beginBlock(
"solving problem with dirichlet border condition using sparse qr solver");
203 const Calculus::PrimalIdentity0 laplace =
calculus.laplace<
PRIMAL>();
205 trace.info() <<
"laplace = " << laplace << endl;
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;
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;
222 Calculus::PrimalForm0 analytic_solution(
calculus);
224 const Calculus::Index dirac_position = 17;
225 const Calculus::Index length = analytic_solution.length();
226 for (Calculus::Index kk=0; kk<length; kk++)
228 Calculus::Scalar alpha = (kk+1.)/(dirac_position+1.);
229 if (kk>dirac_position)
231 alpha = 1. - (kk-dirac_position)/(1.*length-dirac_position);
233 analytic_solution.myContainer(kk) = alpha;
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();
241 trace.info() << solver.isValid() <<
" " << solver.myLinearAlgebraSolver.info() << endl;
244 std::ofstream handle(
"linear_structure_dirichlet.dat");
245 for (Calculus::Index kk=0; kk<analytic_solution.length(); kk++)
247 handle << solved_solution_ordered.myContainer(kk) <<
" " << analytic_solution.myContainer(kk) << endl;
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 );
259 board << solved_solution;
260 board.
saveSVG(
"linear_structure_dirichlet_solution.svg");
264 Calculus::PrimalForm1 solved_solution_gradient =
calculus.derivative<0,
PRIMAL>() * solved_solution;
269 board << solved_solution_gradient;
270 board <<
calculus.sharp(solved_solution_gradient);
271 board.saveSVG(
"linear_structure_dirichlet_solution_gradient.svg");
445 trace.beginBlock(
"testing 2d operators");
451 Calculus primal_calculus;
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)) );
463 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(
Point(1,2)) );
466 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(
Point(3,2), Calculus::KSpace::NEG) );
467 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(
Point(3,2), Calculus::KSpace::POS) );
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) );
475 primal_calculus.eraseCell( primal_calculus.myKSpace.uCell(
Point(1,2)) );
478 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(
Point(3,3)) );
479 primal_calculus.insertSCell( primal_calculus.myKSpace.sCell(
Point(3,5)) );
481 primal_calculus.updateIndexes();
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++)
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);
492 ASSERT( signed_cell == primal_calculus.getSCell(
dim,
PRIMAL, property.index) );
496 <<
" " << signed_cell
497 <<
" " <<
property.primal_size
498 <<
" " <<
property.dual_size
499 <<
" " <<
property.index
500 <<
" " << (
property.flipped ?
"negative" :
"positive")
507 Calculus dual_calculus;
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)) );
519 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(
Point(6,3)) );
522 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(
Point(8,3), Calculus::KSpace::NEG) );
523 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(
Point(8,3), Calculus::KSpace::POS) );
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) );
531 dual_calculus.eraseCell( dual_calculus.myKSpace.uCell(
Point(6,3)) );
534 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(
Point(8,4)) );
535 dual_calculus.insertSCell( dual_calculus.myKSpace.sCell(
Point(8,6)) );
537 dual_calculus.updateIndexes();
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++)
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);
548 ASSERT( signed_cell == dual_calculus.getSCell(
dim,
PRIMAL, property.index) );
552 <<
" " << signed_cell
553 <<
" " <<
property.primal_size
554 <<
" " <<
property.dual_size
555 <<
" " <<
property.index
556 <<
" " << (
property.flipped ?
"negative" :
"positive")
565 board << primal_calculus;
566 board << dual_calculus;
567 board.
saveSVG(
"operators_structure.svg");
570 trace.beginBlock(
"base operators");
572 const Calculus::PrimalDerivative0 primal_d0 = primal_calculus.derivative<0,
PRIMAL>();
573 const Calculus::DualDerivative0 dual_d0p = dual_calculus.derivative<0,
DUAL>();
578#if defined(TEST_HARDCODED_ORDER)
579 Eigen::MatrixXd d0_th(7, 6);
588 FATAL_ERROR( Eigen::MatrixXd(primal_d0.myContainer) == d0_th );
590 Eigen::MatrixXd d0p_th(7, 6);
599 FATAL_ERROR( Eigen::MatrixXd(dual_d0p.myContainer) == d0p_th );
603 const Calculus::PrimalDerivative1 primal_d1 = primal_calculus.derivative<1,
PRIMAL>();
604 const Calculus::DualDerivative1 dual_d1p = dual_calculus.derivative<1,
DUAL>();
609#if defined(TEST_HARDCODED_ORDER)
610 Eigen::MatrixXd d1_th(2, 7);
613 0, 0, -1, -1, 0, -1, -1;
614 FATAL_ERROR( Eigen::MatrixXd(primal_d1.myContainer) == d1_th );
616 Eigen::MatrixXd d1p_th(2, 7);
618 -1, -1, -1, 0, 0, -1, 0,
620 FATAL_ERROR( Eigen::MatrixXd(dual_d1p.myContainer) == d1p_th );
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) );
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>();
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) );
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>();
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) );
664 const Calculus::DualDerivative0 primal_d0p = primal_calculus.derivative<0,
DUAL>();
665 const Calculus::PrimalDerivative0 dual_d0 = dual_calculus.derivative<0,
PRIMAL>();
670#if defined(TEST_HARDCODED_ORDER)
671 Eigen::MatrixXd d0p_th_transpose(2, 7);
674 0, 0, -1, -1, 0, -1, -1;
675 FATAL_ERROR( Eigen::MatrixXd(primal_d0p.myContainer) == d0p_th_transpose.transpose() );
677 Eigen::MatrixXd minus_d0_th_transpose(2, 7);
678 minus_d0_th_transpose <<
679 -1, -1, -1, 0, 0, -1, 0,
681 FATAL_ERROR( Eigen::MatrixXd(dual_d0.myContainer) == -minus_d0_th_transpose.transpose() );
685 const Calculus::DualDerivative1 primal_d1p = primal_calculus.derivative<1,
DUAL>();
686 const Calculus::PrimalDerivative1 dual_d1 = dual_calculus.derivative<1,
PRIMAL>();
691#if defined(TEST_HARDCODED_ORDER)
692 Eigen::MatrixXd minus_d1p_th_transpose(7, 6);
693 minus_d1p_th_transpose <<
701 FATAL_ERROR( Eigen::MatrixXd(primal_d1p.myContainer) == -minus_d1p_th_transpose.transpose() );
703 Eigen::MatrixXd d1_th_transpose(7, 6);
712 FATAL_ERROR( Eigen::MatrixXd(dual_d1.myContainer) == d1_th_transpose.transpose() );
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>();
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) );
736 trace.beginBlock(
"laplace operators");
743 trace.beginBlock(
"sharp operators");
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);
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);
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");
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) );
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);
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);
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");
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) );
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);
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);
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");
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) );
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);
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);
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");
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) );
877 trace.beginBlock(
"flat operators");
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);
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);
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);
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");
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);
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");
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 );
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 );
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);
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);
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);
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");
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);
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");
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 );
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 );