29 #ifndef __Boundadapted_H
30 #define __Boundadapted_H 1
33 #include <boost/lambda/if.hpp>
34 #include <boost/numeric/ublas/banded.hpp>
44 template<
class Convex,u
int16_type O,
typename T>
class PointSetWarpBlend;
45 template<u
int16_type Dim,u
int16_type RealDim,u
int16_type Degree,
typename NormalizationPolicy,
typename T,
template<
class>
class StoragePolicy>
class Dubiner;
47 template<uint16_type Dim,
50 template<
class>
class StoragePolicy>
53 template<uint16_type Dim,
57 struct BoundaryAdaptedTraits
59 static const uint16_type nDim = Dim;
60 static const uint16_type nOrder = Degree;
61 static const uint16_type nConvexOrderDiff = nDim+nOrder+1;
62 static const bool is_normalized =
false;
67 template<u
int16_type order,
typename V = value_type>
74 template<
typename NewT>
75 struct ChangeValueType
78 typedef BoundaryAdaptedTraits<Dim, Degree, NewT, StoragePolicy> traits_type;
82 template<u
int16_type NewOrder>
85 typedef BoundaryAdapted<Dim, NewOrder, T, StoragePolicy> type;
86 typedef BoundaryAdaptedTraits<Dim, NewOrder, T, StoragePolicy> traits_type;
93 typedef typename Convex<nOrder>::type convex_type;
94 typedef typename Convex<nOrder>::reference_type reference_convex_type;
95 typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
96 typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
97 typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
98 mpl::identity<PointSetWarpBlend<diff_convex_type,nConvexOrderDiff,value_type> >,
99 mpl::identity<PointSetEquiSpaced<diff_convex_type, nConvexOrderDiff,value_type> > >::type::type diff_pointset_type;
101 static const uint16_type numVertices = reference_convex_type::numVertices;
102 static const uint16_type numFaces = reference_convex_type::numFaces;
108 typedef StoragePolicy<value_type> storage_policy;
109 typedef typename storage_policy::vector_type vector_type;
110 typedef typename storage_policy::matrix_type matrix_type;
111 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
112 typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
113 typedef typename storage_policy::matrix_node_type matrix_node_type;
114 typedef typename storage_policy::points_type points_type;
115 typedef typename storage_policy::node_type node_type;
118 template<
int D,
int O>
119 struct BoundaryAdaptedTag
121 static const int Dim = D;
122 static const int Order = O;
143 template<uint16_type Dim,
146 template<
class>
class StoragePolicy = StorageUBlas>
147 class BoundaryAdapted : Basis<BoundaryAdaptedTag<Dim, Degree>, T>
149 typedef Basis<BoundaryAdaptedTag<Dim, Degree>, T > super;
152 typedef BoundaryAdaptedTraits<Dim, Degree, T, StoragePolicy> traits_type;
154 static const uint16_type nDim = traits_type::nDim;
155 static const uint16_type nOrder = traits_type::nOrder;
156 static const uint16_type nConvexOrderDiff = traits_type::nConvexOrderDiff;
157 static const bool is_normalized = traits_type::is_normalized;
158 static const uint16_type numVertices = traits_type::numVertices;
159 static const uint16_type numFaces = traits_type::numFaces;
165 typedef BoundaryAdapted<Dim, Degree, T, StoragePolicy> self_type;
171 typedef self_type basis_type;
176 typedef T value_type;
177 typedef int16_type int_type;
178 typedef uint16_type uint_type;
183 typedef typename traits_type::convex_type convex_type;
184 typedef typename traits_type::reference_convex_type reference_convex_type;
185 typedef typename traits_type::diff_pointset_type diff_pointset_type;
190 typedef typename traits_type::storage_policy storage_policy;
191 typedef typename traits_type::vector_type vector_type;
192 typedef typename traits_type::matrix_type matrix_type;
193 typedef typename traits_type::vector_matrix_type vector_matrix_type;
194 typedef typename traits_type::vector_vector_matrix_type vector_vector_matrix_type;
195 typedef typename traits_type::matrix_node_type matrix_node_type;
196 typedef typename traits_type::points_type points_type;
197 typedef typename traits_type::node_type node_type;
199 typedef Principal<Degree, T,StoragePolicy> principal_type;
210 M_pts( nDim, numVertices ),
211 M_pts_face( reference_convex_type::numVertices )
213 PointSetEquiSpaced<convex_type,nOrder,value_type> pts;
216 M_pts = pts.pointsByEntity( 0 );
217 DVLOG(2) <<
"[boundaryadapted] pts= " << M_pts <<
"\n";
220 for ( uint16_type e = M_refconvex.entityRange( nDim-1 ).begin();
221 e < M_refconvex.entityRange( nDim-1 ).end();
224 M_pts_face[e] = pts.pointsBySubEntity( nDim-1, e, 1 );
225 DVLOG(2) <<
"[boundaryadapted] face " << e <<
" pts " << M_pts_face[e] <<
"\n";
229 BoundaryAdapted( BoundaryAdapted
const &
d )
232 M_refconvex( d.M_refconvex ),
234 M_pts_face( d.M_pts_face )
255 return M_pts_face[f];
262 self_type
const& operator=( self_type
const&
d )
266 M_refconvex = d.M_refconvex;
273 matrix_type operator()( node_type
const& pt )
const
275 points_type pts( pt.size(), 1 );
276 ublas::column( pts, 0 ) = pt;
280 matrix_type operator()( points_type
const& pts )
const
313 return "dubiner.boundary.adapted";
330 matrix_type coeff()
const
332 return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
341 static matrix_type
evaluate( points_type
const& __pts )
343 return evaluate( __pts, mpl::int_<nDim>() );
346 template<
typename AE>
347 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts )
349 return derivate( __pts, mpl::int_<nDim>() );
368 static evaluate( points_type
const& __pts, mpl::int_<1> )
370 matrix_type E = M_pfunc.
evaluate_1( ublas::row( __pts,0 ) );
372 D.resize( E.size1(), E.size2() );
374 ublas::row( D, 0 ) = ublas::row( E, 0 );
375 ublas::row( D, 1 ) = ublas::row( E, E.size2()-1 );
377 for (
unsigned int i=1; i < E.size2()-1 ; ++i )
379 ublas::row( D, i+1 ) = ublas::row( E, i );
388 template<
typename AE>
390 static derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<1> )
392 FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error(
"invalid points" );
394 vector_matrix_type E( 1 );
395 E[0].resize( nOrder+1, __pts().size2() );
396 E[0] = M_pfunc.derivate_1( ublas::row( __pts(),0 ) );
398 vector_matrix_type D( 1 );
399 D[0].resize( nOrder+1, __pts().size2() );
401 ublas::row( D[0], 0 ) = ublas::row( E[0], 0 );
402 ublas::row( D[0], 1 ) = ublas::row( E[0], E[0].size2()-1 );
404 for (
unsigned int i=1; i < E[0].size2()-1 ; ++i )
406 ublas::row( D[0], i+1 ) = ublas::row( E[0], i );
416 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<2> );
422 template<
typename AE>
423 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> );
429 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<3> );
435 template<
typename AE>
436 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> );
439 reference_convex_type M_refconvex;
441 std::vector<points_type> M_pts_face;
442 static principal_type M_pfunc;
445 template<uint16_type Dim,
448 template<
class>
class StoragePolicy>
449 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::principal_type
450 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::M_pfunc;
452 template<uint16_type Dim,
455 template<
class>
class StoragePolicy>
456 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
459 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
461 details::etas<TRIANGLE, value_type> etas( __pts );
462 vector_type eta1s = ublas::row( etas(), 0 );
463 vector_type eta2s = ublas::row( etas(), 1 );
465 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
466 vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
470 ublas::row( res,0 ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) );
471 ublas::row( res,1 ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) );
473 vector_type ones( ublas::scalar_vector<value_type>( __pts.size2(), value_type( 1.0 ) ) );
475 ublas::row( res,2 ) = ( ones + eta2s ) / 2.0;
484 for ( uint_type q = 1; q < nOrder; ++q )
487 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) );
491 for ( uint_type q = 1; q < nOrder; ++q )
494 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) );
497 ublas::row( res,G_i )*= value_type( -1.0 );
501 for ( uint_type p = 1; p < nOrder; ++p )
504 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) );
512 for ( uint_type p = 1; p < nOrder; ++p )
513 for ( uint_type q = 1; q < nOrder-p; ++q )
516 ublas::row( res,G_i ) = ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) );
525 template<uint16_type Dim,
528 template<
class>
class StoragePolicy>
529 template<
typename AE>
530 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
531 BoundaryAdapted<Dim, Degree, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
533 vector_matrix_type res( 2 );
534 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
535 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
537 details::etas<TRIANGLE, value_type> etas( __pts );
538 vector_type eta1s = ublas::row( etas(), 0 );
539 vector_type eta2s = ublas::row( etas(), 1 );
541 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
542 vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
544 matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
545 vector_matrix_type dpsi_2( M_pfunc.derivate_2( eta2s ) );
547 vector_type ones( ublas::scalar_vector<value_type>( __pts().size2(), value_type( 1.0 ) ) );
548 vector_type d1( value_type( 2.0 )*ublas::element_div( ones,ones-eta2s ) );
549 vector_type d2( ublas::element_div( ones+eta1s ,ones-eta2s ) );
553 ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) );
554 ublas::row( res[1],0 ) = ublas::element_prod( ublas::row( res[0],0 ) , d2 );
556 ublas::row( res[0],0 ) = ublas::element_prod( ublas::row( res[0],0 ) , d1 );
557 ublas::row( res[1],0 ) += ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) );
559 ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) );
560 ublas::row( res[1],1 ) = ublas::element_prod( ublas::row( res[0],1 ) , d2 );
562 ublas::row( res[0],1 ) = ublas::element_prod( ublas::row( res[0],1 ) , d1 );
563 ublas::row( res[1],1 ) += ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) );
565 ublas::row( res[0],2 ) = ublas::scalar_vector<value_type>( __pts().size2(),value_type( 0.0 ) );
566 ublas::row( res[1],2 ) = value_type( 0.5 ) * ones;
575 for ( uint_type q = 1; q < nOrder; ++q )
578 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) );
579 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
581 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
582 ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) );
586 for ( uint_type q = 1; q < nOrder; ++q )
589 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) );
590 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
592 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
593 ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) );
597 ublas::row( res[0],G_i )*= value_type( -1.0 );
598 ublas::row( res[1],G_i )*= value_type( -1.0 );
603 for ( uint_type p = 1; p < nOrder; ++p )
606 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) );
607 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
609 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
610 ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) );
615 for ( uint_type p = 1; p < nOrder; ++p )
616 for ( uint_type q = 1; q < nOrder-p; ++q )
619 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) );
620 ublas::row( res[1],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d2 );
622 ublas::row( res[0],G_i ) = ublas::element_prod( ublas::row( res[0],G_i ) , d1 );
623 ublas::row( res[1],G_i ) += ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) );
632 template<uint16_type Dim,
635 template<
class>
class StoragePolicy>
636 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::matrix_type
639 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
641 FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error(
"invalid space dimension" );
643 details::etas<TETRAHEDRON, value_type> etas( __pts );
644 vector_type eta1s = ublas::row( etas(), 0 );
645 vector_type eta2s = ublas::row( etas(), 1 );
646 vector_type eta3s = ublas::row( etas(), 2 );
648 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
649 vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
650 vector_vector_matrix_type psi_3( M_pfunc.evaluate_3( eta3s ) );
652 ublas::scalar_vector<value_type> ones( eta3s.size(), value_type( 1.0 ) );
656 ublas::row( res,0 ) =
658 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
659 ublas::row( psi_3[0][0],0 )
661 ublas::row( res,1 ) =
663 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
664 ublas::row( psi_3[nOrder][0],0 )
667 ublas::row( res,2 ) =
669 ublas::row( psi_2[0],nOrder ), ublas::row( psi_3[0][nOrder],0 )
672 ublas::row( res,3 ) = ( ones + eta3s ) / 2.0;
681 for ( uint_type q = 1; q < nOrder; ++q )
684 ublas::row( res,G_i ) =
686 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
687 ublas::row( psi_3[nOrder][q],0 )
692 for ( uint_type q = 1; q < nOrder; ++q )
695 ublas::row( res,G_i ) =
697 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
698 ublas::row( psi_3[0][q],0 )
705 ublas::row( res,G_i )*= value_type( -1.0 );
710 for ( uint_type p = 1; p < nOrder; ++p )
713 ublas::row( res,G_i ) =
715 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
716 ublas::row( psi_3[p][0],0 )
721 for ( uint_type r = 1; r < nOrder; ++r )
724 ublas::row( res,G_i ) =
726 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
727 ublas::row( psi_3[0][0],r )
732 for ( uint_type r = 1; r < nOrder; ++r )
735 ublas::row( res,G_i ) =
737 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
738 ublas::row( psi_3[nOrder][0],r )
743 for ( uint_type r = 1; r < nOrder; ++r )
746 ublas::row( res,G_i ) =
747 ublas::element_prod( ublas::row( psi_2[0],nOrder ),ublas::row( psi_3[0][nOrder],r ) );
757 for ( uint_type q = 1; q < nOrder; ++q )
758 for ( uint_type r = 1; r < nOrder-q; ++r )
761 ublas::row( res,G_i ) =
763 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
764 ublas::row( psi_3[nOrder][q],r )
769 for ( uint_type q = 1; q < nOrder; ++q )
770 for ( uint_type r = 1; r < nOrder-q; ++r )
773 ublas::row( res,G_i ) =
775 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
776 ublas::row( psi_3[0][q],r )
781 for ( uint_type p = 1; p < nOrder; ++p )
782 for ( uint_type r = 1; r < nOrder-p; ++r )
785 ublas::row( res,G_i ) =
787 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
788 ublas::row( psi_3[p][0],r )
794 for ( uint_type p = 1; p < nOrder; ++p )
795 for ( uint_type q = 1; q < nOrder-p; ++q )
798 ublas::row( res,G_i ) =
800 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
801 ublas::row( psi_3[p][q],0 )
810 for ( uint_type p = 1; p < nOrder-2; ++p )
811 for ( uint_type q = 1; q < nOrder-1-p; ++q )
812 for ( uint_type r = 1; r < nOrder-p-q; ++r )
815 ublas::row( res,G_i ) =
817 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
818 ublas::row( psi_3[p][q],r )
825 template<uint16_type Dim,
828 template<
class>
class StoragePolicy>
829 template<
typename AE>
830 typename BoundaryAdapted<Dim, Degree, T, StoragePolicy>::vector_matrix_type
833 vector_matrix_type res( 3 );
834 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
835 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
836 res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
838 FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error(
"invalid space dimension" );
840 details::etas<TETRAHEDRON, value_type> etas( __pts );
841 vector_type eta1s = ublas::row( etas(), 0 );
842 vector_type eta2s = ublas::row( etas(), 1 );
843 vector_type eta3s = ublas::row( etas(), 2 );
845 matrix_type psi_1( M_pfunc.evaluate_1( eta1s ) );
846 matrix_type dpsi_1( M_pfunc.derivate_1( eta1s ) );
848 vector_matrix_type psi_2( M_pfunc.evaluate_2( eta2s ) );
849 vector_matrix_type dpsi_2( M_pfunc.derivate_2( eta2s ) );
851 vector_vector_matrix_type psi_3( M_pfunc.evaluate_3( eta3s ) );
852 vector_vector_matrix_type dpsi_3( M_pfunc.derivate_3( eta3s ) );
854 ublas::scalar_vector<value_type> ones( __pts().size2(), value_type( 1.0 ) );
863 matrix_type Jac( 9, __pts().size2() );
867 ublas::row( Jac,0 ) =
868 ublas::element_div( ( -2.0 )*ones, ( ublas::row( __pts(), 1 ) + ublas::row( __pts(), 2 ) ) );
870 ublas::row( Jac,1 ) =
871 ublas::element_div( 2.0*( ones + ublas::row( __pts(), 0 ) ) , ( ublas::row( __pts(), 1 ) + ublas::row( __pts(), 2 ) ) );
873 ublas::row( Jac,1 ) =
874 ublas::element_div( ublas::row( Jac,1 ), ublas::row( __pts(), 1 )+ublas::row( __pts(), 2 ) );
876 ublas::row( Jac,2 ) =
881 ublas::row( Jac,3 ) =
882 ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
884 ublas::row( Jac,4 ) =
885 ublas::element_div( 2.0*ones, ( ones - ublas::row( __pts(), 2 ) ) );
887 ublas::row( Jac,5 ) =
888 ublas::element_div( 2.0*( ones + ublas::row( __pts(), 1 ) ) , ( ones - ublas::row( __pts(), 2 ) ) );
890 ublas::row( Jac,5 ) =
891 ublas::element_div( ublas::row( Jac,5 ), ones - ublas::row( __pts(), 2 ) );
895 ublas::row( Jac,6 ) =
898 ublas::row( Jac,7 ) =
901 ublas::row( Jac,8 ) = ones;
906 matrix_type d1 ( res[1].size1(), res[1].size2() );
907 matrix_type d2 ( res[1].size1(), res[1].size2() );
908 matrix_type d3 ( res[1].size1(), res[1].size2() );
914 ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) ),
915 ublas::row( psi_3[0][0],0 )
919 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) ),
920 ublas::row( psi_3[0][0],0 )
924 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
925 ublas::row( dpsi_3[0][0],0 )
930 ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
931 ublas::row( psi_3[nOrder][0],0 )
935 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) ),
936 ublas::row( psi_3[nOrder][0],0 )
940 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
941 ublas::row( dpsi_3[nOrder][0],0 )
944 ublas::row( d1,2 ) = ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
948 ublas::row( dpsi_2[0],nOrder ), ublas::row( psi_3[0][nOrder],0 )
953 ublas::row( psi_2[0],nOrder ), ublas::row( dpsi_3[0][nOrder],0 )
957 ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
960 ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
963 ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.5 ) );
971 for ( uint_type q = 1; q < nOrder; ++q )
974 ublas::row( d1,G_i ) =
976 ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
977 ublas::row( psi_3[nOrder][q],0 )
979 ublas::row( d2,G_i ) =
981 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) ),
982 ublas::row( psi_3[nOrder][q],0 )
984 ublas::row( d3,G_i ) =
986 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
987 ublas::row( dpsi_3[nOrder][q],0 )
992 for ( uint_type q = 1; q < nOrder; ++q )
995 ublas::row( d1,G_i ) =
997 ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) ),
998 ublas::row( psi_3[0][q],0 )
1000 ublas::row( d2,G_i ) =
1001 ublas::element_prod(
1002 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) ),
1003 ublas::row( psi_3[0][q],0 )
1005 ublas::row( d3,G_i ) =
1006 ublas::element_prod(
1007 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
1008 ublas::row( dpsi_3[0][q],0 )
1015 ublas::row( d1,G_i )*= value_type( -1.0 );
1016 ublas::row( d2,G_i )*= value_type( -1.0 );
1017 ublas::row( d3,G_i )*= value_type( -1.0 );
1022 for ( uint_type p = 1; p < nOrder; ++p )
1025 ublas::row( d1,G_i ) =
1026 ublas::element_prod(
1027 ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) ),
1028 ublas::row( psi_3[p][0],0 )
1030 ublas::row( d2,G_i ) =
1031 ublas::element_prod(
1032 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) ),
1033 ublas::row( psi_3[p][0],0 )
1035 ublas::row( d3,G_i ) =
1036 ublas::element_prod(
1037 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
1038 ublas::row( dpsi_3[p][0],0 )
1043 for ( uint_type r = 1; r < nOrder; ++r )
1046 ublas::row( d1,G_i ) =
1047 ublas::element_prod(
1048 ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],0 ) ),
1049 ublas::row( psi_3[0][0],r )
1051 ublas::row( d2,G_i ) =
1052 ublas::element_prod(
1053 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],0 ) ),
1054 ublas::row( psi_3[0][0],r )
1056 ublas::row( d3,G_i ) =
1057 ublas::element_prod(
1058 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],0 ) ),
1059 ublas::row( dpsi_3[0][0],r )
1064 for ( uint_type r = 1; r < nOrder; ++r )
1067 ublas::row( d1,G_i ) =
1068 ublas::element_prod(
1069 ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
1070 ublas::row( psi_3[nOrder][0],r )
1072 ublas::row( d2,G_i ) =
1073 ublas::element_prod(
1074 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],0 ) ),
1075 ublas::row( psi_3[nOrder][0],r )
1077 ublas::row( d3,G_i ) =
1078 ublas::element_prod(
1079 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],0 ) ),
1080 ublas::row( dpsi_3[nOrder][0],r )
1085 for ( uint_type r = 1; r < nOrder; ++r )
1088 ublas::row( d1,G_i ) = ublas::scalar_vector<value_type>( __pts().size2(), value_type( 0.0 ) );
1089 ublas::row( d2,G_i ) =
1090 ublas::element_prod( ublas::row( dpsi_2[0],nOrder ),ublas::row( psi_3[0][nOrder],r ) );
1091 ublas::row( d3,G_i ) =
1092 ublas::element_prod( ublas::row( psi_2[0],nOrder ),ublas::row( dpsi_3[0][nOrder],r ) );
1102 for ( uint_type q = 1; q < nOrder; ++q )
1103 for ( uint_type r = 1; r < nOrder-q; ++r )
1106 ublas::row( d1,G_i ) =
1107 ublas::element_prod(
1108 ublas::element_prod( ublas::row( dpsi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
1109 ublas::row( psi_3[nOrder][q],r )
1111 ublas::row( d2,G_i ) =
1112 ublas::element_prod(
1113 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( dpsi_2[nOrder],q ) ),
1114 ublas::row( psi_3[nOrder][q],r )
1116 ublas::row( d3,G_i ) =
1117 ublas::element_prod(
1118 ublas::element_prod( ublas::row( psi_1,nOrder ) , ublas::row( psi_2[nOrder],q ) ),
1119 ublas::row( dpsi_3[nOrder][q],r )
1124 for ( uint_type q = 1; q < nOrder; ++q )
1125 for ( uint_type r = 1; r < nOrder-q; ++r )
1128 ublas::row( d1,G_i ) =
1129 ublas::element_prod(
1130 ublas::element_prod( ublas::row( dpsi_1,0 ) , ublas::row( psi_2[0],q ) ),
1131 ublas::row( psi_3[0][q],r )
1133 ublas::row( d2,G_i ) =
1134 ublas::element_prod(
1135 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( dpsi_2[0],q ) ),
1136 ublas::row( psi_3[0][q],r )
1138 ublas::row( d3,G_i ) =
1139 ublas::element_prod(
1140 ublas::element_prod( ublas::row( psi_1,0 ) , ublas::row( psi_2[0],q ) ),
1141 ublas::row( dpsi_3[0][q],r )
1146 for ( uint_type p = 1; p < nOrder; ++p )
1147 for ( uint_type r = 1; r < nOrder-p; ++r )
1150 ublas::row( d1,G_i ) =
1151 ublas::element_prod(
1152 ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],0 ) ),
1153 ublas::row( psi_3[p][0],r )
1155 ublas::row( d2,G_i ) =
1156 ublas::element_prod(
1157 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],0 ) ),
1158 ublas::row( psi_3[p][0],r )
1160 ublas::row( d3,G_i ) =
1161 ublas::element_prod(
1162 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],0 ) ),
1163 ublas::row( dpsi_3[p][0],r )
1169 for ( uint_type p = 1; p < nOrder; ++p )
1170 for ( uint_type q = 1; q < nOrder-p; ++q )
1173 ublas::row( d1,G_i ) =
1174 ublas::element_prod(
1175 ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) ),
1176 ublas::row( psi_3[p][q],0 )
1178 ublas::row( d2,G_i ) =
1179 ublas::element_prod(
1180 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) ),
1181 ublas::row( psi_3[p][q],0 )
1183 ublas::row( d3,G_i ) =
1184 ublas::element_prod(
1185 ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ),
1186 ublas::row( dpsi_3[p][q],0 )
1195 for ( uint_type p = 1; p < nOrder-2; ++p )
1196 for ( uint_type q = 1; q < nOrder-1-p; ++q )
1197 for ( uint_type r = 1; r < nOrder-p-q; ++r )
1200 ublas::row( d1,G_i ) =
1201 ublas::element_prod( ublas::element_prod( ublas::row( dpsi_1,p ) , ublas::row( psi_2[p],q ) ), ublas::row( psi_3[p][q],r ) );
1202 ublas::row( d2,G_i ) =
1203 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( dpsi_2[p],q ) ), ublas::row( psi_3[p][q],r ) );
1204 ublas::row( d3,G_i ) =
1205 ublas::element_prod( ublas::element_prod( ublas::row( psi_1,p ) , ublas::row( psi_2[p],q ) ), ublas::row( dpsi_3[p][q],r ) );
1211 #if defined (FEELPP_HAS_QD_REAL)
1217 matrix_type j0( __pts().size2(), __pts().size2() );
1218 matrix_type j1( __pts().size2(), __pts().size2() );
1219 matrix_type j2( __pts().size2(), __pts().size2() );
1220 matrix_type j4( __pts().size2(), __pts().size2() );
1221 matrix_type j5( __pts().size2(), __pts().size2() );
1223 for ( uint_type i=0; i < __pts().size2(); ++i )
1224 for ( uint_type j=i; j < __pts().size2(); ++j )
1228 j0( i,i ) = ublas::row( Jac,0 )( i );
1229 j1( i,i ) = ublas::row( Jac,1 )( i );
1230 j2( i,i ) = ublas::row( Jac,2 )( i );
1231 j4( i,i ) = ublas::row( Jac,4 )( i );
1232 j5( i,i ) = ublas::row( Jac,5 )( i );
1237 j0( i,j ) = value_type( 0.0 );
1238 j1( i,j ) = value_type( 0.0 );
1239 j2( i,j ) = value_type( 0.0 );
1240 j4( i,j ) = value_type( 0.0 );
1241 j5( i,j ) = value_type( 0.0 );
1242 j0( j,i ) = value_type( 0.0 );
1243 j1( j,i ) = value_type( 0.0 );
1244 j2( j,i ) = value_type( 0.0 );
1245 j4( j,i ) = value_type( 0.0 );
1246 j5( j,i ) = value_type( 0.0 );
1251 ublas::diagonal_matrix<value_type> j0( __pts().size2() );
1252 ublas::diagonal_matrix<value_type> j1( __pts().size2() );
1253 ublas::diagonal_matrix<value_type> j2( __pts().size2() );
1254 ublas::diagonal_matrix<value_type> j4( __pts().size2() );
1255 ublas::diagonal_matrix<value_type> j5( __pts().size2() );
1257 for ( uint_type i=0; i < __pts().size2(); ++i )
1259 j0( i,i ) = ublas::row( Jac,0 )( i );
1260 j1( i,i ) = ublas::row( Jac,1 )( i );
1261 j2( i,i ) = ublas::row( Jac,2 )( i );
1262 j4( i,i ) = ublas::row( Jac,4 )( i );
1263 j5( i,i ) = ublas::row( Jac,5 )( i );
1268 res[0] = ublas::prod( d1, j0 );
1269 res[1] = ublas::prod( d1, j1 ) + ublas::prod( d2, j4 );
1270 res[2] = ublas::prod( d1, j2 ) + ublas::prod( d2, j5 ) + d3;
simplex of dimension Dim
Definition: simplex.hpp:73
std::string familyName() const
Definition: boundadapted.hpp:311
static matrix_type evaluate(points_type const &__pts)
Definition: boundadapted.hpp:341
matrix_type evaluate_1(vector_type const &__pts) const
Definition: principal.hpp:182
BoundaryAdapted Basis on simplex.
Definition: boundadapted.hpp:51
points_type const & points(int f) const
Definition: boundadapted.hpp:253
static matrix_type const & d(uint16_type i)
derivatives of Dubiner polynomials the derivatives are computed at the nodes of the lattice ...
Definition: basis.hpp:156
Dubiner polynomial orthonormal basis.
Definition: boundadapted.hpp:45
points_type const & points() const
Definition: boundadapted.hpp:245
Convex base class.
Definition: convex.hpp:43
Definition: feelpoly/policy.hpp:469
Reference convex.
Definition: refentity.hpp:58
self_type const & basis() const
Definition: boundadapted.hpp:303
uint_type degree() const
Definition: boundadapted.hpp:295