29 #ifndef __TensorisedBoundadapted_H
30 #define __TensorisedBoundadapted_H 1
33 #include <boost/lambda/if.hpp>
60 template<uint16_type Dim,
63 template<
class>
class StoragePolicy = StorageUBlas>
68 static const uint16_type nDim = Dim;
69 static const uint16_type nOrder = Degree;
70 static const uint16_type nConvexOrder = nDim+nOrder+1;
88 typedef int16_type int_type;
89 typedef uint16_type uint_type;
94 typedef Hypercube<nDim, nConvexOrder, nDim> convex_type;
96 typedef typename reference_convex_type::points_type points_type;
98 static const uint16_type numVertices = reference_convex_type::numVertices;
99 static const uint16_type numFaces = reference_convex_type::numFaces;
101 template<u
int16_type order>
104 typedef Hypercube<nDim, order, nDim> type;
110 typedef StoragePolicy<value_type> storage_policy;
111 typedef typename storage_policy::vector_type vector_type;
112 typedef typename storage_policy::matrix_type matrix_type;
113 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
114 typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
115 typedef typename storage_policy::matrix_node_type matrix_node_type;
116 typedef typename storage_policy::node_type node_type;
127 M_pts( nDim, numVertices ),
128 M_pts_face( numVertices )
130 PointSetEquiSpaced<convex_type, nOrder, value_type> pts;
133 M_pts = pts.pointsByEntity( 0 );
136 for ( uint16_type e = M_refconvex.entityRange( nDim-1 ).begin();
137 e < M_refconvex.entityRange( nDim-1 ).end();
140 M_pts_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
143 matrix_type A( ublas::trans(
evaluate( pts.points() ) ) );
144 matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
146 matrix_type C = lu.solve( D );
147 vector_matrix_type
d ( derivate( pts.points() ) );
148 M_D.resize( d.size() );
150 for (
size_type i = 0; i < d.size(); ++i )
152 M_D[i] = ublas::prod( d[i], C );
153 glas::clean( M_D[i] );
159 M_refconvex( d.M_refconvex ),
161 M_pts_face( d.M_pts_face )
185 return M_pts_face[f];
192 self_type
const& operator=( self_type
const& d )
203 matrix_type operator()( node_type
const& pt )
const
205 points_type pts( pt.size(), 1 );
206 ublas::column( pts, 0 ) = pt;
210 matrix_type operator()( points_type
const& pts )
const
243 return "tensorizedboundaryadapted";
261 matrix_type coeff()
const
263 return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
272 matrix_type
evaluate( points_type
const& __pts )
const
274 return evaluate( __pts, mpl::int_<nDim>() );
277 template<
typename AE>
278 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts )
const
280 return derivate( __pts, mpl::int_<nDim>() );
289 matrix_type
const&
d( uint16_type i )
320 evaluate( points_type
const& __pts, mpl::int_<1> )
const
322 return M_pfunc.evaluate_1( ublas::row( __pts,0 ) );
328 template<
typename AE>
330 derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<1> )
const
332 FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error(
"invalid points" );
334 vector_matrix_type D( 1 );
335 D[0].resize( nOrder+1, __pts().size2() );
336 D[0] = M_pfunc.derivate_1( ublas::row( __pts(),0 ) );
345 matrix_type
evaluate( points_type
const& __pts, mpl::int_<2> )
const;
351 template<
typename AE>
352 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
const;
358 matrix_type
evaluate( points_type
const& __pts, mpl::int_<3> )
const;
364 template<
typename AE>
365 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
const;
368 reference_convex_type M_refconvex;
370 std::vector<points_type> M_pts_face;
372 Principal<Degree, T, StoragePolicy> M_pfunc;
374 vector_matrix_type M_D;
377 template<uint16_type Dim,
380 template<
class>
class StoragePolicy>
381 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
384 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
386 vector_type eta1s = ublas::row( __pts, 0 );
387 vector_type eta2s = ublas::row( __pts, 1 );
389 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
390 matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
407 ublas::row( res,0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,0 ) );
409 ublas::row( res,1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,0 ) );
411 ublas::row( res,2 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,nOrder ) );
413 ublas::row( res,3 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,nOrder ) );
421 for ( uint_type p = 1; p < nOrder; ++p )
424 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) );
428 for ( uint_type q = 1; q < nOrder; ++q )
431 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,q ) );
434 ublas::row( res,G_i )*= value_type( -1.0 );
438 for ( uint_type p = 1; p < nOrder; ++p )
441 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) );
444 ublas::row( res,G_i )*= value_type( -1.0 );
448 for ( uint_type q = 1; q < nOrder; ++q )
451 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,q ) );
458 for ( uint_type p = 1; p < nOrder; ++p )
459 for ( uint_type q = 1; q < nOrder; ++q )
462 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) );
471 template<uint16_type Dim,
474 template<
class>
class StoragePolicy>
475 template<
typename AE>
476 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
477 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
const
479 vector_matrix_type res( 2 );
480 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
481 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
483 vector_type eta1s = ublas::row( __pts(), 0 );
484 vector_type eta2s = ublas::row( __pts(), 1 );
486 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
487 matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
489 matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
490 matrix_type dpsi_2( M_pfunc.derivate_1( eta2s ) );
495 ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,0 ) );
496 ublas::row( res[1],0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,0 ) );
498 ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,0 ) );
499 ublas::row( res[1],1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,0 ) );
501 ublas::row( res[0],2 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,nOrder ) );
502 ublas::row( res[1],2 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,nOrder ) );
504 ublas::row( res[0],3 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,nOrder ) );
505 ublas::row( res[1],3 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,nOrder ) );
513 for ( uint_type p = 1; p < nOrder; ++p )
516 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,0 ) );
517 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,0 ) );
521 for ( uint_type q = 1; q < nOrder; ++q )
524 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,q ) );
525 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,q ) );
529 for ( uint_type p = 1; p < nOrder; ++p )
532 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,nOrder ) );
533 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,nOrder ) );
537 for ( uint_type q = 1; q < nOrder; ++q )
540 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,q ) );
541 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,q ) );
548 for ( uint_type p = 1; p < nOrder; ++p )
549 for ( uint_type q = 1; q < nOrder; ++q )
552 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) );
553 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) );
561 template<uint16_type Dim,
564 template<
class>
class StoragePolicy>
565 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
568 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
570 ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
571 ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
572 ublas::vector<value_type> eta3s = ublas::row( __pts, 2 );
574 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
575 matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
576 matrix_type psi_3( M_pfunc.evaluate_1( eta3s ) );
598 ublas::row( res,0 ) =
599 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
601 ublas::row( res,1 ) =
602 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
604 ublas::row( res,2 ) =
605 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
607 ublas::row( res,3 ) =
608 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
610 ublas::row( res,4 ) =
611 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
613 ublas::row( res,5 ) =
614 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
616 ublas::row( res,6 ) =
617 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
619 ublas::row( res,7 ) =
620 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
630 for ( uint_type p = 1; p < nOrder; ++p )
633 ublas::row( res,G_i ) =
634 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
639 for ( uint_type p = 1; p < nOrder; ++p )
642 ublas::row( res,G_i ) =
643 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
648 for ( uint_type p = 1; p < nOrder; ++p )
651 ublas::row( res,G_i ) =
652 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
657 for ( uint_type p = 1; p < nOrder; ++p )
660 ublas::row( res,G_i ) =
661 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
666 for ( uint_type p = 1; p < nOrder; ++p )
669 ublas::row( res,G_i ) =
670 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
675 for ( uint_type p = 1; p < nOrder; ++p )
678 ublas::row( res,G_i ) =
679 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
684 for ( uint_type p = 1; p < nOrder; ++p )
687 ublas::row( res,G_i ) =
688 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
693 for ( uint_type p = 1; p < nOrder; ++p )
696 ublas::row( res,G_i ) =
697 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
702 for ( uint_type p = 1; p < nOrder; ++p )
705 ublas::row( res,G_i ) =
706 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
711 for ( uint_type p = 1; p < nOrder; ++p )
714 ublas::row( res,G_i ) =
715 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
720 for ( uint_type p = 1; p < nOrder; ++p )
723 ublas::row( res,G_i ) =
724 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
729 for ( uint_type p = 1; p < nOrder; ++p )
732 ublas::row( res,G_i ) =
733 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
744 for ( uint_type p = 1; p < nOrder; ++p )
745 for ( uint_type q = 1; q < nOrder; ++q )
748 ublas::row( res,G_i ) =
749 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,0 ) );
754 for ( uint_type p = 1; p < nOrder; ++p )
755 for ( uint_type q = 1; q < nOrder; ++q )
758 ublas::row( res,G_i ) =
759 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( psi_3,q ) );
764 for ( uint_type p = 1; p < nOrder; ++p )
765 for ( uint_type q = 1; q < nOrder; ++q )
768 ublas::row( res,G_i ) =
769 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
774 for ( uint_type p = 1; p < nOrder; ++p )
775 for ( uint_type q = 1; q < nOrder; ++q )
778 ublas::row( res,G_i ) =
779 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,q ) );
784 for ( uint_type p = 1; p < nOrder; ++p )
785 for ( uint_type q = 1; q < nOrder; ++q )
788 ublas::row( res,G_i ) =
789 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
794 for ( uint_type p = 1; p < nOrder; ++p )
795 for ( uint_type q = 1; q < nOrder; ++q )
798 ublas::row( res,G_i ) =
799 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,nOrder ) );
807 for ( uint_type p = 1; p < nOrder; ++p )
808 for ( uint_type q = 1; q < nOrder; ++q )
809 for ( uint_type r = 1; r < nOrder; ++r )
812 ublas::row( res,G_i ) =
813 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,r ) );
819 template<uint16_type Dim,
822 template<
class>
class StoragePolicy>
823 template<
typename AE>
824 typename TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
825 TensorisedBoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
const
827 vector_matrix_type res( 3 );
828 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
829 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
830 res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
833 vector_type eta1s = ublas::row( __pts(), 0 );
834 vector_type eta2s = ublas::row( __pts(), 1 );
835 vector_type eta3s = ublas::row( __pts(), 2 );
837 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
838 matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
840 matrix_type psi_2( M_pfunc.evaluate_1( eta2s ) );
841 matrix_type dpsi_2( M_pfunc.derivate_1( eta2s ) );
843 matrix_type psi_3( M_pfunc.evaluate_1( eta3s ) );
844 matrix_type dpsi_3( M_pfunc.derivate_1( eta3s ) );
848 ublas::row( res[0],0 ) =
849 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
850 ublas::row( res[1],0 ) =
851 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
852 ublas::row( res[2],0 ) =
853 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
855 ublas::row( res[0],1 ) =
856 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
857 ublas::row( res[1],1 ) =
858 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
859 ublas::row( res[2],1 ) =
860 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
862 ublas::row( res[0],2 ) =
863 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
864 ublas::row( res[1],2 ) =
865 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
866 ublas::row( res[2],2 ) =
867 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
869 ublas::row( res[0],3 ) =
870 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
871 ublas::row( res[1],3 ) =
872 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
873 ublas::row( res[2],3 ) =
874 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
876 ublas::row( res[0],4 ) =
877 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
878 ublas::row( res[1],4 ) =
879 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
880 ublas::row( res[2],4 ) =
881 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
883 ublas::row( res[0],5 ) =
884 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
885 ublas::row( res[1],5 ) =
886 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
887 ublas::row( res[2],5 ) =
888 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
890 ublas::row( res[0],6 ) =
891 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
892 ublas::row( res[1],6 ) =
893 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
894 ublas::row( res[2],6 ) =
895 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
897 ublas::row( res[0],7 ) =
898 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
899 ublas::row( res[1],7 ) =
900 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
901 ublas::row( res[2],7 ) =
902 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
912 for ( uint_type p = 1; p < nOrder; ++p )
915 ublas::row( res[0],G_i ) =
916 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,0 ) );
917 ublas::row( res[1],G_i ) =
918 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,0 ) );
919 ublas::row( res[2],G_i ) =
920 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,0 ) );
925 for ( uint_type p = 1; p < nOrder; ++p )
928 ublas::row( res[0],G_i ) =
929 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
930 ublas::row( res[1],G_i ) =
931 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,0 ) );
932 ublas::row( res[2],G_i ) =
933 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,0 ) );
938 for ( uint_type p = 1; p < nOrder; ++p )
941 ublas::row( res[0],G_i ) =
942 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,0 ) );
943 ublas::row( res[1],G_i ) =
944 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,0 ) );
945 ublas::row( res[2],G_i ) =
946 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,0 ) );
951 for ( uint_type p = 1; p < nOrder; ++p )
954 ublas::row( res[0],G_i ) =
955 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,0 ) );
956 ublas::row( res[1],G_i ) =
957 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,0 ) );
958 ublas::row( res[2],G_i ) =
959 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,0 ) );
964 for ( uint_type p = 1; p < nOrder; ++p )
967 ublas::row( res[0],G_i ) =
968 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
969 ublas::row( res[1],G_i ) =
970 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,p ) );
971 ublas::row( res[2],G_i ) =
972 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,p ) );
977 for ( uint_type p = 1; p < nOrder; ++p )
980 ublas::row( res[0],G_i ) =
981 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,p ) );
982 ublas::row( res[1],G_i ) =
983 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,p ) );
984 ublas::row( res[2],G_i ) =
985 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,p ) );
990 for ( uint_type p = 1; p < nOrder; ++p )
993 ublas::row( res[0],G_i ) =
994 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
995 ublas::row( res[1],G_i ) =
996 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,p ) );
997 ublas::row( res[2],G_i ) =
998 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,p ) );
1003 for ( uint_type p = 1; p < nOrder; ++p )
1006 ublas::row( res[0],G_i ) =
1007 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,p ) );
1008 ublas::row( res[1],G_i ) =
1009 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,p ) );
1010 ublas::row( res[2],G_i ) =
1011 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,p ) );
1016 for ( uint_type p = 1; p < nOrder; ++p )
1019 ublas::row( res[0],G_i ) =
1020 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( psi_3,nOrder ) );
1021 ublas::row( res[1],G_i ) =
1022 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,nOrder ) );
1023 ublas::row( res[2],G_i ) =
1024 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,nOrder ) );
1029 for ( uint_type p = 1; p < nOrder; ++p )
1032 ublas::row( res[0],G_i ) =
1033 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
1034 ublas::row( res[1],G_i ) =
1035 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,nOrder ) );
1036 ublas::row( res[2],G_i ) =
1037 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,nOrder ) );
1042 for ( uint_type p = 1; p < nOrder; ++p )
1045 ublas::row( res[0],G_i ) =
1046 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
1047 ublas::row( res[1],G_i ) =
1048 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,nOrder ) );
1049 ublas::row( res[2],G_i ) =
1050 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ), ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,nOrder ) );
1055 for ( uint_type p = 1; p < nOrder; ++p )
1058 ublas::row( res[0],G_i ) =
1059 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( psi_3,nOrder ) );
1060 ublas::row( res[1],G_i ) =
1061 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( dpsi_2,p ) ), ublas::row( psi_3,nOrder ) );
1062 ublas::row( res[2],G_i ) =
1063 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ), ublas::row( psi_2,p ) ), ublas::row( dpsi_3,nOrder ) );
1074 for ( uint_type p = 1; p < nOrder; ++p )
1075 for ( uint_type q = 1; q < nOrder; ++q )
1078 ublas::row( res[0],G_i ) =
1079 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,0 ) );
1080 ublas::row( res[1],G_i ) =
1081 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,0 ) );
1082 ublas::row( res[2],G_i ) =
1083 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,0 ) );
1088 for ( uint_type p = 1; p < nOrder; ++p )
1089 for ( uint_type q = 1; q < nOrder; ++q )
1092 ublas::row( res[0],G_i ) =
1093 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( psi_3,q ) );
1094 ublas::row( res[1],G_i ) =
1095 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,0 ) ), ublas::row( psi_3,q ) );
1096 ublas::row( res[2],G_i ) =
1097 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,0 ) ), ublas::row( dpsi_3,q ) );
1102 for ( uint_type p = 1; p < nOrder; ++p )
1103 for ( uint_type q = 1; q < nOrder; ++q )
1106 ublas::row( res[0],G_i ) =
1107 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
1108 ublas::row( res[1],G_i ) =
1109 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2,p ) ), ublas::row( psi_3,q ) );
1110 ublas::row( res[2],G_i ) =
1111 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2,p ) ), ublas::row( dpsi_3,q ) );
1116 for ( uint_type p = 1; p < nOrder; ++p )
1117 for ( uint_type q = 1; q < nOrder; ++q )
1120 ublas::row( res[0],G_i ) =
1121 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( psi_3,q ) );
1122 ublas::row( res[1],G_i ) =
1123 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,nOrder ) ), ublas::row( psi_3,q ) );
1124 ublas::row( res[2],G_i ) =
1125 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,nOrder ) ), ublas::row( dpsi_3,q ) );
1130 for ( uint_type p = 1; p < nOrder; ++p )
1131 for ( uint_type q = 1; q < nOrder; ++q )
1134 ublas::row( res[0],G_i ) =
1135 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( psi_3,q ) );
1136 ublas::row( res[1],G_i ) =
1137 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2,p ) ), ublas::row( psi_3,q ) );
1138 ublas::row( res[2],G_i ) =
1139 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2,p ) ), ublas::row( dpsi_3,q ) );
1144 for ( uint_type p = 1; p < nOrder; ++p )
1145 for ( uint_type q = 1; q < nOrder; ++q )
1148 ublas::row( res[0],G_i ) =
1149 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,nOrder ) );
1150 ublas::row( res[1],G_i ) =
1151 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,nOrder ) );
1152 ublas::row( res[2],G_i ) =
1153 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,nOrder ) );
1161 for ( uint_type p = 1; p < nOrder; ++p )
1162 for ( uint_type q = 1; q < nOrder; ++q )
1163 for ( uint_type r = 1; r < nOrder; ++r )
1166 ublas::row( res[0],G_i ) =
1167 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( psi_3,r ) );
1168 ublas::row( res[1],G_i ) =
1169 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2,q ) ), ublas::row( psi_3,r ) );
1170 ublas::row( res[2],G_i ) =
1171 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2,q ) ), ublas::row( dpsi_3,r ) );