35 #include <boost/lambda/if.hpp>
41 #include <feel/feelpoly/policy.hpp>
50 template<uint16_type Dim,
53 typename NormalizationPolicy,
55 template<
class>
class StoragePolicy>
59 template<uint16_type Dim,
62 typename NormalizationPolicy = Normalized<true>,
64 template<
class>
class StoragePolicy = StorageUBlas>
67 static const uint16_type nDim = Dim;
68 static const uint16_type nRealDim = RealDim;
69 static const uint16_type nOrder = Degree;
70 static const uint16_type nConvexOrderDiff = nDim+nOrder+1;
71 static const bool is_normalized = NormalizationPolicy::is_normalized;
82 template<u
int16_type order,
typename V = value_type>
85 typedef Simplex<nDim, order, nDim> type;
86 typedef Reference<Simplex<nDim, order, nDim>, nDim, order, nDim, V> reference_type;
89 template<
typename NewT>
90 struct ChangeValueType
92 typedef Dubiner<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> type;
93 typedef DubinerTraits<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> traits_type;
96 template<u
int16_type NewOrder>
99 typedef Dubiner<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> type;
100 typedef DubinerTraits<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> traits_type;
106 typedef typename Convex<nOrder>::type convex_type;
107 typedef typename Convex<nOrder>::reference_type reference_convex_type;
109 typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
110 typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
112 typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
113 mpl::identity<PointSetWarpBlend<diff_convex_type, nConvexOrderDiff, value_type> >,
114 mpl::identity<PointSetEquiSpaced<diff_convex_type, nConvexOrderDiff, value_type> > >::type::type diff_pointset_type;
119 typedef StoragePolicy<value_type> storage_policy;
120 typedef typename storage_policy::matrix_type matrix_type;
121 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
122 typedef typename storage_policy::matrix_node_type matrix_node_type;
123 typedef typename storage_policy::points_type points_type;
124 typedef typename storage_policy::node_type node_type;
127 template<
int D,
int O>
130 static const int Dim = D;
131 static const int Order = O;
161 template<uint16_type Dim,
164 typename NormalizationPolicy = Normalized<true>,
166 template<
class>
class StoragePolicy = StorageUBlas>
172 typedef DubinerTraits<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> traits_type;
174 static const uint16_type nDim = traits_type::nDim;
175 static const uint16_type nRealDim = traits_type::nRealDim;
176 static const uint16_type nOrder = traits_type::nOrder;
177 static const uint16_type nConvexOrderDiff = traits_type::nConvexOrderDiff;
178 static const bool is_normalized = traits_type::is_normalized;
179 static const bool isTransformationEquivalent =
true;
180 static const bool isContinuous =
false;
181 static const bool is_product =
true;
182 typedef Discontinuous continuity_type;
188 typedef Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> self_type;
194 typedef self_type basis_type;
196 typedef typename traits_type::value_type value_type;
201 typedef typename traits_type::convex_type convex_type;
202 typedef typename traits_type::reference_convex_type reference_convex_type;
204 typedef typename traits_type::diff_pointset_type diff_pointset_type;
209 typedef typename traits_type::storage_policy storage_policy;
210 typedef typename traits_type::matrix_type matrix_type;
211 typedef typename traits_type::vector_matrix_type vector_matrix_type;
212 typedef typename traits_type::matrix_node_type matrix_node_type;
213 typedef typename traits_type::points_type points_type;
214 typedef typename traits_type::node_type node_type;
225 M_pts( M_refconvex.makePoints( nDim, 0 ) )
227 this->initDerivation();
229 Dubiner( Dubiner
const &
d )
246 self_type
const& operator=( self_type
const&
d )
258 matrix_type operator()( node_type
const& pt )
const
260 points_type pts( pt.size(), 1 );
261 ublas::column( pts, 0 ) = pt;
265 matrix_type operator()( points_type
const& pts )
const
307 return is_normalized;
345 std::cout <<
"[Dubiner::coeff] coeff = "
346 << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
349 return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
358 static matrix_type
evaluate( points_type
const& __pts )
360 return evaluate( __pts, mpl::int_<nDim>() );
363 template<
typename AE>
364 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts )
366 return derivate( __pts, mpl::int_<nDim>() );
375 static matrix_type
const&
d( uint16_type i )
386 static matrix_type
const&
derivate( uint16_type i )
401 evaluate( points_type
const& __pts, mpl::int_<1> )
403 matrix_type m ( JacobiBatchEvaluation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts, 0 ) ) );
407 for ( uint16_type i = 0; i < m.size1(); ++i )
408 ublas::row( m, i ) *= math::sqrt( value_type( i )+0.5 );
418 template<
typename AE>
419 static vector_matrix_type
420 derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<1> )
422 FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).warn(
"invalid points" );
425 vector_matrix_type D( 1 );
426 D[0].resize( nOrder+1, __pts().size2() );
427 D[0] = JacobiBatchDerivation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts(),0 ) );
430 for ( uint16_type i = 0; i < nOrder+1; ++i )
431 ublas::row( D[0], i ) *= math::sqrt( value_type( i )+0.5 );
440 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<2> );
446 template<
typename AE>
447 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> );
453 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<3> );
459 template<
typename AE>
460 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> );
462 static void initDerivation();
464 reference_convex_type M_refconvex;
471 static bool _S_has_derivation;
477 static std::vector<matrix_type> _S_D;
481 template<uint16_type Dim,
484 typename NormalizationPolicy,
486 template<
class>
class StoragePolicy>
487 bool Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_has_derivation =
false;
489 template<uint16_type Dim,
492 typename NormalizationPolicy,
494 template<
class>
class StoragePolicy>
495 std::vector<typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type>
496 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_D;
498 template<uint16_type Dim,
501 typename NormalizationPolicy,
503 template<
class>
class StoragePolicy>
505 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::initDerivation()
508 typedef typename traits_type::convex_type convex_type;
509 typedef typename traits_type::reference_convex_type reference_convex_type;
511 typedef typename traits_type::diff_pointset_type diff_pointset_type;
513 typedef typename traits_type::storage_policy storage_policy;
514 typedef typename traits_type::matrix_type matrix_type;
515 typedef typename traits_type::vector_matrix_type vector_matrix_type;
516 typedef typename traits_type::matrix_node_type matrix_node_type;
517 typedef typename traits_type::points_type points_type;
518 typedef typename traits_type::node_type node_type;
521 if ( _S_has_derivation ==
false )
523 _S_has_derivation =
true;
525 reference_convex_type refconvex;
528 diff_pointset_type diff_pts( 1 );
529 matrix_type A( evaluate( diff_pts.points() ) );
532 matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
533 LU<matrix_type> lu( A );
534 matrix_type C = lu.solve( D );
536 vector_matrix_type d ( derivate( diff_pts.points() ) );
537 _S_D.resize( d.size() );
539 for (
size_type i = 0; i < d.size(); ++i )
541 _S_D[i] = ublas::prod( d[i], C );
542 glas::clean( _S_D[i] );
549 template<uint16_type Dim,
552 typename NormalizationPolicy,
554 template<
class>
class StoragePolicy>
555 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
558 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
560 details::etas<TRIANGLE, value_type> etas( __pts );
561 ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
562 ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
567 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
568 std::vector<matrix_type> bs( nOrder+1 );
570 for (
int i = 0; i < nOrder+1; ++i )
572 bs[ i ].resize( nOrder-i, eta2s.size() );
573 bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
577 details::scalings<nOrder, T> scalings( eta2s );
580 for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
582 for ( uint16_type i = 0; i < k+1; ++i,++cur )
584 uint16_type ii = k-i;
589 value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*( value_type( ii+jj )+1.0 ) );
591 for ( uint16_type l = 0; l < as.size2(); ++l )
592 res( cur, l ) = normalization*as( ii,l )*scalings()( ii,l )*bs[ii]( jj,l );
597 for ( uint16_type l = 0; l < as.size2(); ++l )
598 res( cur, l ) = as( ii,l )*scalings()( ii,l )*bs[ii]( jj,l );
606 template<uint16_type Dim,
609 typename NormalizationPolicy,
611 template<
class>
class StoragePolicy>
612 template<
typename AE>
613 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
614 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
616 vector_matrix_type res( 2 );
617 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
618 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
621 details::etas<TRIANGLE, value_type> etas( __pts );
622 ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
623 ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
629 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
630 matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
632 std::vector<matrix_type> bs( nOrder+1 );
633 std::vector<matrix_type> dbs( nOrder+1 );
635 for ( uint16_type i = 0; i < nOrder+1; ++i )
637 bs[ i ].resize( nOrder-i, eta2s.size() );
638 dbs[ i ].resize( nOrder-i, eta2s.size() );
639 bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
640 dbs[ i ] = dyna::JacobiBatchDerivation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
645 details::scalings<nOrder, T> scalings( eta2s );
647 ublas::vector<value_type> one( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
648 ublas::vector<value_type> tmp( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
651 for ( uint16_type k = 0, cur = 0; k < nOrder+1; ++k )
653 for ( uint16_type i = 0; i < k+1; ++i, ++cur )
655 uint16_type ii = k-i;
660 ublas::row( res[0], cur ) = ublas::element_prod( ublas::row( das, ii ),
661 ublas::row( bs[ii], jj ) );
664 ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
665 ublas::row( scalings(), ii-1 ) );
668 ublas::row( res[1], cur ) = ublas::element_prod( ublas::row( das, ii ),
669 ublas::row( bs[ii], jj ) );
670 ublas::row( res[1], cur ) = 0.5 * element_prod( ublas::row( res[1], cur ), ( one+eta1s ) );
673 ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
674 ublas::row( scalings(), ii-1 ) );
677 tmp = ublas::element_prod( ublas::row( scalings(), ii ),
678 ublas::row( dbs[ii], jj ) );
681 tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings(), ii-1 ),
682 ublas::row( bs[ii], jj ) );
685 ublas::row( res[1], cur ) += ublas::element_prod( ublas::row( as, ii ), tmp );
690 value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*( value_type( ii+jj )+1.0 ) );
691 ublas::row( res[0], cur ) *= normalization;
692 ublas::row( res[1], cur ) *= normalization;
700 template<uint16_type Dim,
703 typename NormalizationPolicy,
705 template<
class>
class StoragePolicy>
706 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
709 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
711 FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error(
"invalid space dimension" );
713 details::etas<TETRAHEDRON, value_type> etas( __pts );
714 ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
715 ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
716 ublas::vector<value_type> eta3s = ublas::row( etas(), 2 );
721 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
722 std::vector<matrix_type> bs( nOrder+1 );
723 ublas::matrix<matrix_type> cs( nOrder+1, nOrder+1 );
725 for (
int i = 0; i < nOrder+1; ++i )
727 bs[ i ].resize( nOrder-i, eta2s.size() );
728 bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
730 for (
int j = 0; j < nOrder+1-i; ++j )
732 cs( i, j ).resize( nOrder-i-j, eta3s.size() );
733 cs( i, j ) = dyna::JacobiBatchEvaluation( nOrder-i-j,
734 value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
739 details::scalings<nOrder, T> scalings2( eta2s );
740 details::scalings<nOrder, T> scalings3( eta3s );
742 for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
744 for ( uint16_type i = 0; i < k+1; ++i )
746 for ( uint16_type j = 0; j < k+1-i; ++j,++cur )
748 uint16_type ii = k-i-j;
755 value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*
756 ( value_type( ii+jj )+1.0 )*
757 ( value_type( ii+jj+kk )+1.5 ) );
759 for ( uint16_type l = 0; l < as.size2(); ++l )
760 res( cur, l ) = normalization*( as( ii,l )*
761 scalings2()( ii,l )*bs[ii]( jj,l )*
762 scalings3()( ii+jj,l )*cs( ii, jj )( kk,l ) );
767 for ( uint16_type l = 0; l < as.size2(); ++l )
768 res( cur, l ) = as( ii,l )*
769 scalings2()( ii,l )*bs[ii]( jj,l )*
770 scalings3()( ii+jj,l )*cs( ii, jj )( kk,l );
779 template<uint16_type Dim,
782 typename NormalizationPolicy,
784 template<
class>
class StoragePolicy>
785 template<
typename AE>
786 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
787 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
789 vector_matrix_type res( 3 );
790 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
791 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
792 res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
794 FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error(
"invalid space dimension" );
796 details::etas<TETRAHEDRON, value_type> etas( __pts );
797 ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
798 ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
799 ublas::vector<value_type> eta3s = ublas::row( etas(), 2 );
804 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
805 matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
806 std::vector<matrix_type> bs( nOrder+1 );
807 std::vector<matrix_type> dbs( nOrder+1 );
808 ublas::matrix<matrix_type> cs( nOrder+1, nOrder+1 );
809 ublas::matrix<matrix_type> dcs( nOrder+1, nOrder+1 );
811 for (
int i = 0; i < nOrder+1; ++i )
813 bs[ i ].resize( nOrder-i, eta2s.size() );
814 dbs[ i ].resize( nOrder-i, eta2s.size() );
815 bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
816 dbs[ i ] = dyna::JacobiBatchDerivation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
818 for (
int j = 0; j < nOrder+1-i; ++j )
820 cs( i, j ).resize( nOrder-i-j, eta3s.size() );
821 dcs( i, j ).resize( nOrder-i-j, eta3s.size() );
822 cs( i, j ) = dyna::JacobiBatchEvaluation( nOrder-i-j,
823 value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
824 dcs( i, j ) = dyna::JacobiBatchDerivation( nOrder-i-j,
825 value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
830 details::scalings<nOrder, T> scalings2( eta2s );
831 details::scalings<nOrder, T> scalings3( eta3s );
834 ublas::vector<value_type> one( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
835 ublas::vector<value_type> tmp( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
838 for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
840 for ( uint16_type i = 0; i < k+1; ++i )
842 for ( uint16_type j = 0; j < k+1-i; ++j,++cur )
844 uint16_type ii = k-i-j;
850 ublas::row( res[0], cur ) = ublas::element_prod( ublas::row( das, ii ),
851 ublas::row( bs[ii], jj ) );
852 ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
853 ublas::row( cs( ii, jj ), kk ) );
856 ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
857 ublas::row( scalings2(), ii-1 ) );
860 ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
861 ublas::row( scalings3(), ii+jj-1 ) );
864 ublas::row( res[1], cur ) = ublas::element_prod( ublas::row( das, ii ),
865 ublas::row( bs[ii], jj ) );
866 ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
867 ublas::row( cs( ii, jj ), kk ) );
868 ublas::row( res[1], cur ) = 0.5 * element_prod( ublas::row( res[1], cur ),
872 ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
873 ublas::row( scalings2(), ii-1 ) );
876 ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
877 ublas::row( scalings3(), ii+jj-1 ) );
880 tmp = ublas::element_prod( ublas::row( scalings2(), ii ),
881 ublas::row( dbs[ii], jj ) );
884 tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings2(), ii-1 ),
885 ublas::row( bs[ii], jj ) );
887 tmp = ublas::element_prod( tmp,
888 ublas::row( as, ii ) );
889 tmp = ublas::element_prod( tmp,
890 ublas::row( cs( ii, jj ), kk ) );
893 tmp = ublas::element_prod( tmp,
894 ublas::row( scalings3(), ii+jj-1 ) );
897 ublas::row( res[1], cur ) += tmp;
900 ublas::row( res[2], cur ) = ublas::element_prod( ublas::row( das, ii ),
901 ublas::row( bs[ii], jj ) );
902 ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
903 ublas::row( cs( ii, jj ), kk ) );
904 ublas::row( res[2], cur ) = 0.5 * element_prod( ublas::row( res[2], cur ),
908 ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
909 ublas::row( scalings2(), ii-1 ) );
912 ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
913 ublas::row( scalings3(), ii+jj-1 ) );
916 tmp = ublas::element_prod( ublas::row( scalings2(), ii ),
917 ublas::row( dbs[ii], jj ) );
920 tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings2(), ii-1 ),
921 ublas::row( bs[ii], jj ) );
923 tmp = ublas::element_prod( tmp,
924 ublas::row( as, ii ) );
925 tmp = ublas::element_prod( tmp,
926 ublas::row( cs( ii, jj ), kk ) );
927 tmp = 0.5 * element_prod( tmp, ( one+eta2s ) );
930 tmp = ublas::element_prod( tmp,
931 ublas::row( scalings3(), ii+jj-1 ) );
934 ublas::row( res[2], cur ) += tmp;
937 tmp = ublas::element_prod( ublas::row( scalings3(), ii+jj ),
938 ublas::row( dcs( ii, jj ), kk ) );
941 tmp -= 0.5*( ii+jj )*ublas::element_prod( ublas::row( cs( ii, jj ), kk ),
942 ublas::row( scalings3(), ii+jj-1 ) );
944 tmp = ublas::element_prod( tmp,
945 ublas::row( as, ii ) );
946 tmp = ublas::element_prod( tmp,
947 ublas::row( bs[ ii ], jj ) );
948 tmp = ublas::element_prod( tmp,
949 ublas::row( scalings2(), ii ) );
952 ublas::row( res[2], cur ) += tmp;
956 value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*
957 ( value_type( ii+jj )+1.0 )*
958 ( value_type( ii+jj+kk )+1.5 ) );
960 ublas::row( res[0], cur ) *= normalization;
961 ublas::row( res[1], cur ) *= normalization;
962 ublas::row( res[2], cur ) *= normalization;