31 #ifndef __PolynomialSet_H
32 #define __PolynomialSet_H 1
35 #include <boost/multi_array.hpp>
36 #include <boost/multi_array/extent_gen.hpp>
37 #include <boost/foreach.hpp>
38 #include <boost/optional.hpp>
39 #include <boost/mpl/min_max.hpp>
42 #include <feel/feelpoly/context.hpp>
48 #include <feel/feelpoly/tensorisedboundadapted.hpp>
74 template<
typename Poly,
template<u
int16_type>
class PolySetType = Scalar >
83 static const uint16_type nDim = Poly::nDim;
84 static const uint16_type nRealDim = Poly::nRealDim;
85 static const uint16_type nOrder = Poly::nOrder;
92 typedef PolynomialSet<Poly, PolySetType> self_type;
93 typedef boost::shared_ptr<self_type> self_ptrtype;
94 typedef typename Poly::value_type value_type;
95 typedef typename Poly::basis_type basis_type;
96 static const bool is_product = Poly::is_product;
99 typedef PolySetType<nRealDim> polyset_type;
100 static const bool is_tensor2 = polyset_type::is_tensor2;
101 static const bool is_vectorial = polyset_type::is_vectorial;
102 static const bool is_scalar = polyset_type::is_scalar;
103 static const uint16_type nComponents = polyset_type::nComponents;
104 static const uint16_type nComponents1 = polyset_type::nComponents1;
105 static const uint16_type nComponents2 = polyset_type::nComponents2;
106 static const uint16_type rank = polyset_type::rank;
108 typedef PolynomialSet<Poly, Scalar> component_type;
109 typedef Polynomial<Poly, PolySetType> polynomial_type;
111 typedef polynomial_type polynomial_view_type;
113 typedef typename Poly::convex_type convex_type;
114 typedef typename basis_type::matrix_type matrix_type;
115 typedef typename basis_type::points_type points_type;
118 typedef PolynomialSet<Poly, Vectorial> gradient_polynomialset_type;
120 BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, value_type>::value ) );
121 BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, typename points_type::value_type>::value ) );
137 PolynomialSet( Poly
const& p )
139 M_basis( p.
basis() ),
140 M_coeff( p.
coeff() ),
147 PolynomialSet( Poly
const& p, matrix_type
const& c,
bool __as_is =
false )
149 M_basis( p.
basis() ),
150 M_coeff( p.
coeff() ),
165 PolynomialSet( matrix_type
const& c,
bool __as_is =
false )
180 PolynomialSet( PolynomialSet
const & p )
182 M_basis( p.M_basis ),
183 M_coeff( p.M_coeff ),
191 virtual ~PolynomialSet()
200 self_type& operator=( self_type
const& pset )
204 M_basis = pset.M_basis;
205 M_coeff = pset.M_coeff;
218 BOOST_STATIC_ASSERT( is_vectorial );
219 FEELPP_ASSERT( i < nComponents )( i )( nComponents ).error (
"invalid component index" );
220 const int nrows = M_coeff.size1()/nComponents;
221 const int ncols = M_coeff.size2();
223 ublas::slice( nrows*i+i, nComponents, nrows/nComponents ),
224 ublas::slice( 0, 1, ncols ) ),
true );
277 static uint16_type numberOfComponents()
287 return M_coeff.size1()/nComponents;
295 return M_coeff.size2();
303 return ublas::norm_frobenius( M_coeff ) < Feel::type_traits<value_type>::epsilon();
326 std::string
name( std::string sep =
"." )
const
328 std::ostringstream os;
329 os << this->
familyName() << sep << nDim << sep << nOrder;
347 M_coeff = ublas::prod( __c, M_coeff );
357 M_coeff = ublas::prod( __c, polyset_type::toMatrix( M_coeff ) );
358 M_coeff = polyset_type::toType( M_coeff );
381 size_type new_dim_p = nComponents*list_p.size();
382 matrix_type
coeff( nComponents*nComponents*list_p.size(), M_coeff.size2() );
384 BOOST_FOREACH( uint16_type i, list_p )
386 for (
int c = 0; c < nComponents; ++c )
389 ublas::range( c*new_dim_p+nComponents*j, c*dim_p+nComponents*j+nComponents ),
390 ublas::range( 0, M_coeff.size2() ) ) =
392 ublas::range( c*dim_p+nComponents*i, c*dim_p+nComponents*i+nComponents ),
393 ublas::range( 0, M_coeff.size2() ) );
409 matrix_type
coeff( nComponents*nComponents*dim_p, M_coeff.size2() );
411 for (
int c = 0; c < nComponents; ++c )
415 ublas::range( c*nComponents*dim_p, ( c+1 )*nComponents*dim_p ),
416 ublas::range( 0, M_coeff.size2() ) ) =
418 ublas::range( nc, nc+nComponents*dim_p ),
419 ublas::range( 0, M_coeff.size2() ) );
435 uint16_type dim_p = dim_top-dim_bot;
436 matrix_type
coeff( nComponents*nComponents*dim_p,
439 for (
int c = 0; c < nComponents; ++c )
443 ublas::range( c*nComponents*dim_bot, ( c+1 )*nComponents*dim_top ),
444 ublas::range( 0, M_coeff.size2() ) ) =
446 ublas::range( nc+dim_bot, nc+nComponents*dim_top ),
447 ublas::range( 0, M_coeff.size2() ) );
462 matrix_type
coeff( nComponents, M_coeff.size2() );
466 ublas::range( 0, nComponents ),
467 ublas::range( 0, M_coeff.size2() ) ) =
470 ublas::range( nComponents*i, nComponents*( i+1 ) ),
471 ublas::range( 0, M_coeff.size2() ) );
483 template<
typename AE>
484 ublas::vector<value_type>
evaluate( uint16_type i, ublas::vector_expression<AE>
const& __pt )
const
486 return ublas::row( ublas::prod( M_coeff, M_basis( __pt ) ), i );
493 template<
typename AE>
494 ublas::matrix<value_type>
evaluate( ublas::vector_expression<AE>
const& __pt )
const
496 return ublas::prod( M_coeff, M_basis( __pt ) );
510 template<
typename AE>
511 matrix_type
evaluate( ublas::matrix_expression<AE>
const& __pts )
const
513 matrix_type m ( M_basis.evaluate( __pts ) );
514 FEELPP_ASSERT( M_coeff.size2() == m.size1() )( M_coeff.size2() )( m.size1() ).error(
"invalid size" );
515 return ublas::prod( M_coeff, m );
539 matrix_type
const&
d( uint16_type i )
const
541 return M_basis.d( i );
544 matrix_type
d( uint16_type i, uint16_type j )
const
546 return ublas::prod( M_basis.d( i ), M_basis.d( j ) );
556 return self_type( Poly(), ublas::prod( M_coeff, M_basis.d( l ) ),
true );
559 template<
typename AE>
560 ublas::vector<matrix_type>
derivate( ublas::matrix_expression<AE>
const& pts )
const
562 ublas::vector<matrix_type> der( M_basis.derivate( pts ) );
563 ublas::vector<matrix_type> res( nDim );
565 for ( uint16_type i = 0; i < nDim; ++i )
567 res[i].resize( M_coeff.size1(), pts().size2() );
568 ublas::axpy_prod( M_coeff, der[i], res[i] );
574 template<
typename AE>
575 matrix_type
derivate( uint16_type i, ublas::matrix_expression<AE>
const& pts )
const
577 ublas::vector<matrix_type> der( M_basis.derivate( pts ) );
578 matrix_type res( M_coeff.size1(), pts().size2() );
579 ublas::axpy_prod( M_coeff, der[i], res );
583 template<
typename AE>
584 matrix_type
derivate( uint16_type i, uint16_type j, ublas::matrix_expression<AE>
const& pts )
const
587 matrix_type eval( M_basis.evaluate( pts ) );
590 matrix_type p1 = ublas::prod( M_coeff, M_basis.d( i ) );
591 matrix_type p2 = ublas::prod( p1, M_basis.d( j ) );
592 return ublas::prod( p2, eval );
599 gradient_polynomialset_type
602 return gradient( mpl::int_<polyset_type::rank>() );
605 gradient_polynomialset_type
608 const int n1 = M_coeff.size1();
609 const int n2 = M_coeff.size2();
610 ublas::matrix<value_type> c ( nDim*nDim*n1, n2 );
613 for (
int i = 0; i <nDim; ++i )
616 ublas::slice( nDim*n1*i+i, nDim, n1 ),
617 ublas::slice( 0, 1, n2 ) ) = ublas::prod( M_coeff, M_basis.d( i ) );
620 return gradient_polynomialset_type( c,
true );
623 gradient_polynomialset_type
628 const int n1 = M_coeff.size1();
629 const int n2 = M_coeff.size2();
630 ublas::matrix<value_type> c ( nComponents*nComponents*n1, n2 );
633 std::cout <<
"M_coeff = " << M_coeff <<
"\n"
634 <<
"c = " << c <<
"\n";
637 for (
int i = 0; i <nDim; ++i )
640 std::cout <<
"i=" << i <<
"\n"
642 << ublas::prod( M_coeff, M_basis.d( i ) ) <<
"\n"
645 ublas::slice( nDim*n1*i+i, nDim, n1 ),
646 ublas::slice( 0, 1, n2 ) ) <<
"\n";
650 ublas::slice( nDim*n1*i+i, nDim, n1 ),
651 ublas::slice( 0, 1, n2 ) ) = ublas::prod( M_coeff, M_basis.d( i ) );
655 return gradient_polynomialset_type( c,
true );
666 return M_coeff.size1();
675 FEELPP_ASSERT( p.coeff().size2() ==
coeff().size2() )( p.coeff().size2() )(
coeff().size2() ).warn(
"invalid polynomial set" );
678 std::cout <<
"coeff = " << M_coeff <<
"\n"
679 <<
"p = " << p.coeff() <<
"\n";
689 matrix_type oldcoeff = M_coeff;
690 M_coeff.resize( M_coeff.size1() + p.coeff().size1(), p.coeff().size2(), false );
693 if ( oldcoeff.size1() > 0 )
696 ublas::range( 0, oldcoeff.size1() ),
697 ublas::range( 0, oldcoeff.size2() ) ) = oldcoeff;
701 ublas::range( oldcoeff.size1(), oldcoeff.size1()+p.coeff().size1() ),
702 ublas::range( 0, oldcoeff.size2() ) ) = p.coeff();
718 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
720 typedef boost::shared_ptr<reference_element_type> reference_element_ptrtype;
722 typedef typename reference_element_type::value_type value_type;
724 static const uint16_type nDim = reference_element_type::nDim;
725 static const uint16_type nComponents1 = reference_element_type::nComponents1;
726 static const uint16_type nComponents2 = reference_element_type::nComponents2;
728 static const uint16_type nComponents = reference_element_type::nComponents;
730 typedef typename reference_element_type::points_type matrix_node_t_type;
731 typedef ublas::matrix<value_type> matrix_type;
732 typedef Eigen::Matrix<value_type,nComponents1,1> id_type;
733 typedef Eigen::Matrix<value_type,nComponents1,nDim> g_type;
734 typedef Eigen::Matrix<value_type,nDim,nDim> h_type;
735 typedef boost::multi_array<id_type,2> functionvalue_type;
736 typedef boost::multi_array<g_type,2> grad_type;
737 typedef boost::multi_array<h_type,2> hessian_type;
745 PreCompute( matrix_node_t_type
const& __pts )
747 M_ref_ele( new reference_element_type() ),
753 init( M_ref_ele, __pts, mpl::int_<rank>() );
761 PreCompute( reference_element_ptrtype
const& __ref_ele,
762 matrix_node_t_type
const& __pts )
764 M_ref_ele( __ref_ele ),
770 init( M_ref_ele, __pts, mpl::int_<rank>() );
774 PreCompute( PreCompute
const& __pc )
776 M_ref_ele( __pc.M_ref_ele ),
777 M_nodes( __pc.M_nodes ),
779 M_grad( __pc.M_grad ),
780 M_hessian( __pc.M_hessian )
788 PreCompute&
operator=( PreCompute
const& __pc )
792 M_ref_ele = __pc.M_ref_ele;
793 M_nodes = __pc.M_nodes;
795 M_grad = __pc.M_grad;
796 M_hessian = __pc.M_hessian;
802 void update( matrix_node_t_type
const& __pts )
805 init( M_ref_ele, __pts, mpl::int_<rank>() );
811 uint16_type dim()
const
813 return reference_element_type::nDim;
820 uint16_type nComputedNodes()
const
822 return M_nodes.size2();
829 uint16_type nPoints()
const
831 return M_nodes.size2();
838 matrix_node_t_type
const& nodes()
const
847 ublas::matrix_column<matrix_node_t_type const> node( uint16_type __i )
const
849 return ublas::column( M_nodes, __i );
858 functionvalue_type
const& phi()
const
867 value_type phi( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
869 FEELPP_ASSERT( q < nComputedNodes() )( q )( nComputedNodes() ).error(
"invalid node index" );
870 FEELPP_ASSERT( i < M_ref_ele->
nbDof() )( i )( M_ref_ele->nbDof() ).error(
"invalid dof index" );
872 return M_phi[i][q]( c1,c2 );
875 grad_type
const& grad()
const
880 value_type grad(
size_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
882 return M_grad[i][q]( c1,c2 );
885 hessian_type
const&
hessian()
const
890 value_type
hessian(
size_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
892 return M_hessian[i][q]( c1,c2 );
898 init( reference_element_ptrtype
const& M_ref_ele,
899 matrix_node_t_type
const& __pts,
902 M_phi.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
903 M_grad.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
904 M_hessian.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
906 matrix_type phiv = M_ref_ele->evaluate( __pts );
907 ublas::vector<matrix_type> __grad( M_ref_ele->derivate( __pts ) );
908 matrix_type __hessian( M_ref_ele->gradient().gradient().evaluate( __pts ) );
910 typedef typename grad_type::index index;
911 const index I = M_ref_ele->nbDof();
912 const index Q = __pts.size2();
914 for ( index i = 0; i < I; ++i )
916 for ( index q = 0; q < Q; ++q )
917 M_phi[i][q]( 0,0 ) = phiv( i, q );
919 for ( index q = 0; q < Q; ++q )
920 for ( index j = 0; j < nDim; ++j )
921 M_grad[i][q]( 0,j ) = __grad[j]( i, q );
923 for ( index q = 0; q < Q; ++q )
924 for ( index j = 0; j < nDim; ++j )
925 for ( index k = j; k < nDim; ++k )
927 value_type t = __hessian( nDim*nDim*I*( nDim*k+j )+nDim*nDim*i+nDim*j+k, q );
928 M_hessian[i][q]( j,k ) = t;
929 M_hessian[i][q]( k,j ) = t;
935 init( reference_element_ptrtype
const& M_ref_ele,
936 matrix_node_t_type
const& __pts,
939 typedef typename grad_type::index index;
941 std::cout <<
"family = " << M_ref_ele->familyName() << std::endl;
942 std::cout <<
"is_product = " << reference_element_type::is_product << std::endl;
943 std::cout <<
"nbDof = " << M_ref_ele->nbDof() << std::endl;
945 const index I = ( reference_element_type::is_product?M_ref_ele->nbDof()/nRealDim/nRealDim:M_ref_ele->nbDof()/nRealDim );
949 const index Q = __pts.size2();
952 const int ncdof= ( reference_element_type::is_product?nRealDim:1 );
953 int nldof= ( reference_element_type::is_product?I*nRealDim:I );
955 M_phi.resize( boost::extents[nldof][__pts.size2()] );
956 M_grad.resize( boost::extents[nldof][__pts.size2()] );
958 matrix_type phiv = M_ref_ele->evaluate( __pts );
959 ublas::vector<matrix_type> __grad( M_ref_ele->derivate( __pts ) );
961 for ( index i = 0; i < I; ++i )
964 for ( index c1 = 0; c1 < ncdof; ++c1 )
966 for ( index q = 0; q < Q; ++q )
967 for ( index j = 0; j < nRealDim; ++j )
968 M_phi[I*c1+i][q]( j,0 ) = phiv( nldof*c1+nRealDim*i+j, q );
972 for ( index q = 0; q < Q; ++q )
973 for ( index j = 0; j < nRealDim; ++j )
974 for ( index l = 0; l < nDim; ++l )
977 M_grad[I*c1+i][q]( j,l ) = __grad[l]( nldof*c1+nRealDim*i+j, q );
986 reference_element_ptrtype M_ref_ele;
987 matrix_node_t_type M_nodes;
988 functionvalue_type M_phi;
990 hessian_type M_hessian;
994 typedef boost::shared_ptr<precompute_type> precompute_ptrtype;
997 preCompute( self_ptrtype p, points_type
const& P )
999 return precompute_ptrtype(
new PreCompute( p, P ) );
1002 typedef std::vector<std::map<typename convex_type::permutation_type, precompute_ptrtype> > faces_precompute_type;
1004 std::vector<std::map<typename convex_type::permutation_type, precompute_ptrtype> >
1005 preComputeOnFaces( self_ptrtype p, points_type
const& P )
1008 QuadMapped<pointset_type> qm;
1009 typedef typename QuadMapped<pointset_type>::permutation_type permutation_type;
1010 typename QuadMapped<pointset_type>::permutation_points_type ppts( qm( P ) );
1012 typedef typename convex_type::permutation_type permutation_type;
1013 std::vector<std::map<permutation_type, precompute_ptrtype> > geopc( convex_type::numTopologicalFaces );
1015 for ( uint16_type __f = 0; __f < convex_type::numTopologicalFaces; ++__f )
1017 for ( permutation_type __p( permutation_type::IDENTITY );
1018 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1029 template<
size_type context_v,
typename Basis_t,
typename Geo_t,
typename ElementType,
size_type context_g = context_v>
1033 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1034 static const int context = context_v;
1036 static const uint16_type PDim = ElementType::nDim;
1038 static const uint16_type NDim = ElementType::nRealDim;
1039 static const uint16_type nDof = Basis_t::nLocalDof;
1040 static const bool is_product = Basis_t::is_product;
1041 static const bool is_scalar = Basis_t::is_scalar;
1042 static const bool is_vectorial = Basis_t::is_vectorial;
1043 static const bool is_tensor2 = Basis_t::is_tensor2;
1044 static const uint16_type nComponents = Basis_t::nComponents;
1045 static const uint16_type nComponents1 = Basis_t::nComponents1;
1046 static const uint16_type nComponents2 = Basis_t::nComponents2;
1048 typedef typename Basis_t::polyset_type polyset_type;
1049 static const uint16_type rank = polyset_type::rank;
1051 typedef Basis_t reference_element_type;
1052 typedef boost::shared_ptr<Basis_t> reference_element_ptrtype;
1054 typedef typename reference_element_type::value_type value_type;
1056 typedef ElementType geometric_element_type;
1057 typedef typename Geo_t::template Context<context_g, ElementType> geometric_mapping_context_type;
1058 typedef boost::shared_ptr<geometric_mapping_context_type> geometric_mapping_context_ptrtype;
1060 typedef typename node<value_type>::type node_type;
1062 typedef ublas::matrix<value_type> matrix_type;
1063 typedef points_type matrix_node_t_type;
1065 typedef value_type dn_type;
1066 typedef value_type grad_type;
1067 typedef value_type div_type;
1070 typedef typename matrix_type::value_type phi_type;
1071 typedef typename matrix_type::value_type dphi_type;
1072 typedef typename matrix_type::value_type id_type;
1074 typedef Eigen::Matrix<value_type,nComponents1,1> id_type;
1075 typedef Eigen::Matrix<value_type,nComponents1,nDim> ref_grad_type;
1076 typedef Eigen::Matrix<value_type,nComponents1,NDim> grad_type;
1077 typedef Eigen::Matrix<value_type,NDim,NDim> hess_type;
1078 typedef Eigen::Matrix<value_type,1,1> div_type;
1079 typedef Eigen::Matrix<value_type,nComponents1,1> dn_type;
1080 typedef Eigen::Matrix<value_type,3,1> curl_type;
1081 typedef Eigen::Matrix<value_type,nComponents1,1> dx_type;
1082 typedef Eigen::Matrix<value_type,nComponents1,1> dy_type;
1083 typedef Eigen::Matrix<value_type,nComponents1,1> dz_type;
1086 typedef boost::multi_array<id_type,2> functionvalue_type;
1087 typedef boost::multi_array<g_type,2> grad_type;
1088 typedef boost::multi_array<h_type,2> hessian_type;
1089 typedef boost::multi_array<n_type,2> dn_type;
1090 typedef boost::multi_array<d_type,2> div_type;
1094 typedef geometric_mapping_context_type gmc_type;
1095 typedef Eigen::Matrix<value_type,Eigen::Dynamic, Eigen::Dynamic> matrix_eigen_type;
1096 typedef Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::PDim> matrix_eigen_NP_type;
1097 typedef Eigen::Matrix<value_type,gmc_type::PDim,gmc_type::NDim> matrix_eigen_PN_type;
1098 typedef Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::NDim> matrix_eigen_NN_type;
1099 typedef Eigen::Matrix<value_type,nComponents1,NDim> matrix_eigen_grad_type;
1100 typedef typename Eigen::Map<const Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::PDim,Eigen::RowMajor> > matrix_eigen_ublas_NP_type;
1101 typedef typename Eigen::Map<const Eigen::Matrix<value_type,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > matrix_eigen_ublas_type;
1104 template<u
int16_type TheRank = polyset_type::rank+2>
1107 static const uint16_type rank = TheRank;
1108 typedef boost::array<size_type,rank> index_type;
1109 typedef boost::detail::multi_array::extent_gen<rank> extents_type;
1110 template<u
int16_type>
friend class Rank;
1114 init( mpl::int_<rank>() );
1116 Index( Index
const& i )
1118 M_index( i.M_index ),
1119 M_extents( i.M_extents ),
1126 Index( Index<rank-1>
const& __index )
1128 M_extents( __index.extents()[RankUp<polyset_type>::type::nComponentsLast] )
1130 std::copy( __index.beginIndex(), __index.endIndex(), M_index.begin() );
1137 Index( Index<rank+1>
const& index )
1143 std::copy( index.M_extents.begin(),
1144 boost::prior( index.M_extents.end() ),
1145 M_extents.begin() );
1151 Index
const&
operator=( Index
const& i )
1155 M_index = i.M_index;
1156 M_extents = i.M_extents;
1162 typename index_type::iterator beginIndex()
1164 return M_index.begin();
1166 typename index_type::const_iterator beginIndex()
const
1168 return M_index.begin();
1170 typename index_type::iterator endIndex()
1172 return M_index.end();
1174 typename index_type::const_iterator endIndex()
const
1176 return M_index.end();
1179 template<
typename Tuple>
1180 void setIndex( Tuple
const& tu )
1182 setIndex( tu, mpl::int_<rank>() );
1185 void setIndex( uint16_type c,
size_type i )
1193 return index( mpl::int_<rank>() );
1198 return nComponents*nDof*M_index[1] + nComponents*M_index[0] + M_index[1];
1201 uint16_type component()
const
1206 Index<rank+1> rankUp()
const
1208 return Index<rank+1>( M_extents[RankUp<polyset_type>::type::nComponentsLast] );
1213 return index( mpl::int_<rank>() );
1216 extents_type extents()
const
1222 void init( mpl::int_<2> )
1224 M_extents = boost::extents[nDof][nComponents];
1227 void init( mpl::int_<3> )
1229 M_extents = boost::extents[nDof][nComponents][nComponents];
1234 void init( mpl::int_<4> )
1236 M_extents = boost::extents[nDof][nComponents][nComponents][nComponents];
1247 return nDof*M_index[1] + M_index[0];
1251 const size_type d = M_extents.ranges_[1].size();
1253 const size_type n = M_extents.ranges_[0].size();
1255 return ( d2*n*( d*M_index[1]+M_index[2] ) +
1257 d*M_index[1] + M_index[2] );
1259 return ( d2*n*M_index[1]*M_index[2] +
1260 d2*M_index[0]*M_index[2] +
1261 d*M_index[1]*M_index[2] +
1268 template<
typename Tuple>
1269 void setIndex( Tuple
const& tu, mpl::int_<2> )
1271 M_index[ 0 ] = boost::get<0>( tu );
1274 template<
typename Tuple>
1275 void setIndex( Tuple
const& tu, mpl::int_<3> )
1277 M_index[ 0 ] = boost::get<0>( tu );
1278 M_index[ 1 ] = boost::get<1>( tu );
1279 M_index[ 2 ] = boost::get<2>( tu );
1280 M_comp = M_index[ 1 ];
1282 template<
typename Tuple>
1283 void setIndex( Tuple
const& tu, mpl::int_<4> )
1285 M_index[ 0 ] = boost::get<0>( tu );
1286 M_index[ 1 ] = boost::get<1>( tu );
1287 M_index[ 2 ] = boost::get<2>( tu );
1293 extents_type M_extents;
1298 Context( reference_element_ptrtype
const& __RefEle,
1299 geometric_mapping_context_ptrtype
const& __gmc,
1300 precompute_ptrtype
const& __pc )
1303 M_npoints( __pc->nPoints() ),
1306 M_ref_ele( __RefEle ),
1309 M_phi( __pc->phi() ),
1310 M_gradphi( __pc->grad() ),
1317 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1319 const int ntdof = nDof*nComponents1;
1320 M_grad.resize( boost::extents[ntdof][M_npoints] );
1321 M_dx.resize( boost::extents[ntdof][M_npoints] );
1322 M_dy.resize( boost::extents[ntdof][M_npoints] );
1323 M_dz.resize( boost::extents[ntdof][M_npoints] );
1325 if ( vm::has_first_derivative_normal<context>::value )
1327 M_dn.resize( boost::extents[ntdof][M_npoints] );
1330 if ( vm::has_div<context>::value )
1332 M_div.resize( boost::extents[ntdof][M_npoints] );
1335 if ( vm::has_curl<context>::value )
1337 M_curl.resize( boost::extents[ntdof][M_npoints] );
1340 if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1342 M_hessian.resize( boost::extents[ntdof][M_npoints] );
1353 void transformationEquivalence( geometric_mapping_context_ptrtype
const& __gmc,
1369 void transformationEquivalence( geometric_mapping_context_ptrtype
const& __gmc,
1373 M_ref_ele->transform( __gmc, M_pc.get(), M_phi,
1374 M_gradphi, ( vm::has_div<context>::value || vm::has_curl<context>::value || vm::has_grad<context>::value || vm::has_first_derivative<context>::value ),
1375 M_hessphi, ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1379 void update( geometric_mapping_context_ptrtype
const& __gmc,
1380 precompute_ptrtype
const& __pc )
1385 if ( ( M_npoints != M_gmc->nPoints() ) ||
1386 ( M_npoints != __pc->nPoints() ) ||
1387 ( M_npoints != M_phi.shape()[1] ) )
1393 M_npoints = __pc->nPoints();
1395 const int ntdof = nDof*nComponents1;
1396 M_phi.resize( boost::extents[ntdof][M_npoints] );
1397 M_gradphi.resize( boost::extents[ntdof][M_npoints] );
1399 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1401 M_grad.resize( boost::extents[ntdof][M_npoints] );
1402 M_dx.resize( boost::extents[ntdof][M_npoints] );
1403 M_dy.resize( boost::extents[ntdof][M_npoints] );
1404 M_dz.resize( boost::extents[ntdof][M_npoints] );
1406 if ( vm::has_div<context>::value )
1408 M_div.resize( boost::extents[ntdof][M_npoints] );
1411 if ( vm::has_curl<context>::value )
1413 M_curl.resize( boost::extents[ntdof][M_npoints] );
1416 if ( vm::has_first_derivative_normal<context>::value )
1418 M_dn.resize( boost::extents[ntdof][M_npoints] );
1421 if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1423 M_hessian.resize( boost::extents[ntdof][M_npoints] );
1428 M_phi = M_pc.get()->phi();
1429 M_gradphi = M_pc.get()->grad();
1433 void update( geometric_mapping_context_ptrtype
const& __gmc )
1437 static const bool do_opt= ( nOrder<=1 ) && ( Geo_t::nOrder==1 ) && ( convex_type::is_simplex );
1438 transformationEquivalence( __gmc, mpl::bool_<Basis_t::isTransformationEquivalent>() );
1441 for (
int i = 0; i < M_gradphi.num_elements(); ++i )
1442 std::cout <<
"M_gradphi[" << i <<
"]=" << *( M_gradphi.data()+i ) <<
"\n";
1445 update( __gmc, mpl::int_<rank>(), mpl::bool_<do_opt>() );
1448 void update( geometric_mapping_context_ptrtype
const& __gmc, mpl::int_<0>, mpl::bool_<true> )
1454 geometric_mapping_context_type* thegmc = __gmc.get();
1456 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1458 const uint16_type Q = M_npoints;
1459 const uint16_type I = nDof;
1462 matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1463 matrix_eigen_grad_type grad_real = matrix_eigen_grad_type::Zero();
1465 for ( uint16_type i = 0; i < I; ++i )
1467 grad_real.noalias() = M_gradphi[i][0]*Bt.transpose();
1468 for ( uint16_type q = 0; q < Q; ++q )
1470 M_grad[i][q] = grad_real;
1471 M_dx[i][q] = M_grad[i][q].col( 0 );
1474 M_dy[i][q] = M_grad[i][q].col( 1 );
1477 M_dz[i][q] = M_grad[i][q].col( 2 );
1484 if ( vm::has_first_derivative_normal<context>::value )
1486 std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1487 const uint16_type I = M_ref_ele->nbDof()*nComponents;
1488 const uint16_type Q = nPoints();
1490 for (
int i = 0; i < I; ++i )
1492 for ( uint16_type q = 0; q < Q; ++q )
1494 for ( uint16_type l = 0; l < NDim; ++l )
1496 value_type n = thegmc->unitNormal( l, 0 );
1497 value_type gn = M_grad[i][0]( 0,l ) * n;
1498 M_dn[i][q]( 0,0 ) += gn;
1507 if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1510 const uint16_type Q = M_npoints;
1511 const uint16_type I = nDof;
1514 boost::multi_array<value_type,4>
const& B3 = thegmc->B3();
1516 std::fill( M_hessian.data(), M_hessian.data()+M_hessian.num_elements(), hess_type::Zero() );
1519 for ( uint16_type i = 0; i < I; ++i )
1521 for ( uint16_type q = 0; q < Q; ++q )
1524 for ( uint16_type l = 0; l < NDim; ++l )
1526 for ( uint16_type j = 0; j < NDim; ++j )
1528 for ( uint16_type p = 0; p < PDim; ++p )
1530 for ( uint16_type r = 0; r < PDim; ++r )
1533 value_type h = B3[l][j][p][r] * __pc->hessian( i, p, r, q );
1534 M_hessian[i][q]( j,l ) += h;
1544 void update( geometric_mapping_context_ptrtype
const& __gmc, mpl::int_<0>, mpl::bool_<false> )
1548 geometric_mapping_context_type* thegmc = __gmc.get();
1550 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1552 const uint16_type Q = M_npoints;
1553 const uint16_type I = nDof;
1555 for ( uint16_type i = 0; i < I; ++i )
1557 for ( uint16_type q = 0; q < Q; ++q )
1559 matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1560 M_grad[i][q] = M_gradphi[i][q]*Bt.transpose();
1561 M_dx[i][q] = M_grad[i][q].col( 0 );
1564 M_dy[i][q] = M_grad[i][q].col( 1 );
1567 M_dz[i][q] = M_grad[i][q].col( 2 );
1572 if ( vm::has_first_derivative_normal<context>::value )
1574 std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1575 const uint16_type I = M_ref_ele->nbDof()*nComponents;
1576 const uint16_type Q = nPoints();
1578 for (
int i = 0; i < I; ++i )
1580 for ( uint16_type q = 0; q < Q; ++q )
1582 for ( uint16_type l = 0; l < NDim; ++l )
1584 M_dn[i][q]( 0,0 ) += M_grad[i][q]( 0,l ) * thegmc->unitNormal( l, q );
1592 if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1595 const uint16_type Q = M_npoints;
1596 const uint16_type I = nDof;
1599 boost::multi_array<value_type,4>
const& B3 = thegmc->B3();
1601 std::fill( M_hessian.data(), M_hessian.data()+M_hessian.num_elements(), hess_type::Zero() );
1604 for ( uint16_type i = 0; i < I; ++i )
1606 for ( uint16_type q = 0; q < Q; ++q )
1608 for ( uint16_type l = 0; l < NDim; ++l )
1611 for ( uint16_type j = 0; j < NDim; ++j )
1613 for ( uint16_type p = 0; p < PDim; ++p )
1617 for ( uint16_type r = 0; r < PDim; ++r )
1620 value_type h = B3[l][j][p][r] * __pc->hessian( i, p, r, q );
1621 M_hessian[i][q]( l,j ) += h;
1631 void update( geometric_mapping_context_ptrtype
const& __gmc, mpl::int_<1>, mpl::bool_<false> )
1634 geometric_mapping_context_type* thegmc = __gmc.get();
1636 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1638 const uint16_type Q = M_npoints;
1639 const uint16_type I = nDof;
1641 typedef typename boost::multi_array<value_type,4>::index_range range;
1643 for ( uint16_type ii = 0; ii < I; ++ii )
1645 int ncomp= ( reference_element_type::is_product?nComponents1:1 );
1647 for ( uint16_type c = 0; c < ncomp; ++c )
1649 uint16_type i = I*c + ii;
1651 for ( uint16_type q = 0; q < Q; ++q )
1654 matrix_eigen_ublas_type Bt ( thegmc->B( q ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1656 M_grad[i][q] = M_gradphi[i][q]*Bt.transpose();
1657 M_dx[i][q] = M_grad[i][q].col( 0 );
1660 M_dy[i][q] = M_grad[i][q].col( 1 );
1663 M_dz[i][q] = M_grad[i][q].col( 2 );
1666 if ( vm::has_div<context>::value )
1668 M_div[i][q]( 0,0 ) = M_grad[i][q].trace();
1672 if ( vm::has_curl<context>::value )
1676 M_curl[i][q]( 0 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1677 M_curl[i][q]( 1 ) = M_curl[i][q]( 0 );
1678 M_curl[i][q]( 2 ) = M_curl[i][q]( 0 );
1681 else if ( NDim == 3 )
1683 M_curl[i][q]( 0 ) = M_grad[i][q]( 2,1 ) - M_grad[i][q]( 1,2 );
1684 M_curl[i][q]( 1 ) = M_grad[i][q]( 0,2 ) - M_grad[i][q]( 2,0 );
1685 M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1693 if ( vm::has_first_derivative_normal<context>::value )
1695 std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1696 const uint16_type I = nDof*nComponents1;
1697 const uint16_type Q = nPoints();
1699 for (
int i = 0; i < I; ++i )
1700 for ( uint16_type q = 0; q < Q; ++q )
1702 for ( uint16_type c1 = 0; c1 < NDim; ++c1 )
1704 for ( uint16_type l = 0; l < NDim; ++l )
1706 M_dn[i][q]( c1,0 ) += M_grad[i][q]( c1,l ) * thegmc->unitNormal( l, q );
1713 void update( geometric_mapping_context_ptrtype
const& __gmc, mpl::int_<1>, mpl::bool_<true> )
1716 geometric_mapping_context_type* thegmc = __gmc.get();
1718 if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1720 const uint16_type Q = M_npoints;
1721 const uint16_type I = nDof;
1723 typedef typename boost::multi_array<value_type,4>::index_range range;
1725 matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1726 matrix_eigen_grad_type grad_real = matrix_eigen_grad_type::Zero();
1729 for ( uint16_type ii = 0; ii < I; ++ii )
1731 int ncomp= ( reference_element_type::is_product?nComponents1:1 );
1733 for ( uint16_type c = 0; c < ncomp; ++c )
1736 uint16_type i = I*c + ii;
1738 grad_real.noalias() = M_gradphi[i][0]*Bt.transpose();
1740 for ( uint16_type q = 0; q < Q; ++q )
1742 M_grad[i][q] = grad_real;
1743 M_dx[i][q] = M_grad[i][q].col( 0 );
1746 M_dy[i][q] = M_grad[i][q].col( 1 );
1749 M_dz[i][q] = M_grad[i][q].col( 2 );
1752 if ( vm::has_div<context>::value )
1754 M_div[i][q]( 0,0 ) = M_grad[i][q].trace();
1758 if ( vm::has_curl<context>::value )
1763 M_curl[i][q]( 0 ) = 0;
1764 M_curl[i][q]( 1 ) = 0;
1766 M_curl[i][q]( 0 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1767 M_curl[i][q]( 1 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1769 M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1772 else if ( NDim == 3 )
1774 M_curl[i][q]( 0 ) = M_grad[i][q]( 2,1 ) - M_grad[i][q]( 1,2 );
1775 M_curl[i][q]( 1 ) = M_grad[i][q]( 0,2 ) - M_grad[i][q]( 2,0 );
1776 M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1784 if ( vm::has_first_derivative_normal<context>::value )
1786 std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1787 const uint16_type I = nDof*nComponents1;
1788 const uint16_type Q = nPoints();
1790 for (
int i = 0; i < I; ++i )
1791 for ( uint16_type q = 0; q < Q; ++q )
1793 for ( uint16_type c1 = 0; c1 < NDim; ++c1 )
1795 for ( uint16_type l = 0; l < NDim; ++l )
1797 M_dn[i][q]( c1,0 ) += M_grad[i][q]( c1,l ) * thegmc->unitNormal( l, q );
1810 uint16_type nPoints()
const
1818 geometric_mapping_context_ptrtype
const& gmContext()
const
1830 matrix_node_t_type
const& xRefs()
const
1832 return M_gmc->xRefs();
1836 precompute_ptrtype
const& pc()
const
1841 boost::multi_array<value_type,4>
const& id()
const
1846 value_type
const& id( uint32_type i,
1849 uint32_type q )
const
1851 return id( i, c1, c2, q, mpl::int_<rank>() );
1854 value_type
const& id( uint32_type i,
1858 mpl::int_<0> )
const
1860 return M_phi[i][q]( 0,0 );
1863 value_type
const& id( uint32_type i,
1867 mpl::int_<1> )
const
1869 detail::ignore_unused_variable_warning( c2 );
1870 return M_phi[i][q]( c1,0 );
1873 id_type
const& id( uint32_type i, uint32_type q )
const
1878 value_type
const& id( uint32_type i,
1882 mpl::int_<2> )
const
1884 return M_phi[i][q]( c1,c2 );
1887 value_type
d( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
1889 return M_grad[i][q]( c1,c2 );
1892 value_type
dx( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
1894 return dx( i, c1, c2, q, mpl::int_<rank>() );
1896 value_type
dx( uint32_type i, uint16_type , uint16_type , uint32_type q, mpl::int_<0> )
const
1898 BOOST_MPL_ASSERT_MSG( nDim >= 1, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<0> ) );
1899 return M_grad[i][q]( 0,0 );
1901 value_type
dx( uint32_type i, uint16_type c1, uint16_type , uint32_type q, mpl::int_<1> )
const
1903 BOOST_MPL_ASSERT_MSG( nDim >= 1, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<0> ) );
1904 return M_grad[i][q]( c1,0 );
1906 value_type
dy( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
1908 return dy( i, c1, c2, q, mpl::int_<rank>() );
1910 value_type
dy( uint32_type i, uint16_type , uint16_type , uint32_type q, mpl::int_<0> )
const
1912 BOOST_MPL_ASSERT_MSG( nDim >= 2, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<1> ) );
1913 return M_grad[i][q]( 0,1 );
1915 value_type
dy( uint32_type i, uint16_type c1, uint16_type , uint32_type q, mpl::int_<1> )
const
1917 BOOST_MPL_ASSERT_MSG( nDim >= 2, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<1> ) );
1918 return M_grad[i][q]( c1,1 );
1920 value_type
dz( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
1922 return dz( i, c1, c2, q, mpl::int_<rank>() );
1924 value_type
dz( uint32_type i, uint16_type , uint16_type , uint32_type q, mpl::int_<0> )
const
1926 BOOST_MPL_ASSERT_MSG( nDim >= 3, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<2> ) );
1927 return M_grad[i][q]( 0,2 );
1929 value_type
dz( uint32_type i, uint16_type c1, uint16_type , uint32_type q, mpl::int_<1> )
const
1931 BOOST_MPL_ASSERT_MSG( nDim >= 3, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<2> ) );
1932 return M_grad[i][q]( c1,2 );
1946 value_type
const& dn( uint32_type i,
1949 uint32_type q )
const
1951 return M_dn[i][q]( c1,c2 );
1954 dn_type
const& dn( uint16_type i, uint32_type q )
const
1959 value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
1961 return grad( i, c1, c2, q, mpl::int_<rank>() );
1963 value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
1965 detail::ignore_unused_variable_warning( c1 );
1966 detail::ignore_unused_variable_warning( c2 );
1967 return M_grad[i][q]( 0,c2 );
1969 value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
1971 return M_grad[i][q]( c1,c2 );
1974 grad_type
const& grad( uint16_type i, uint32_type q )
const
1976 return M_grad[i][q];
1978 dx_type
const&
dx( uint16_type i, uint32_type q )
const
1982 dy_type
const&
dy( uint16_type i, uint32_type q )
const
1986 dz_type
const&
dz( uint16_type i, uint32_type q )
const
1999 value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2001 detail::ignore_unused_variable_warning( c1 );
2002 detail::ignore_unused_variable_warning( c2 );
2003 return div( i, c1, c2, q, mpl::int_<rank>() );
2006 value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2008 detail::ignore_unused_variable_warning( i );
2009 detail::ignore_unused_variable_warning( q );
2010 detail::ignore_unused_variable_warning( c1 );
2011 detail::ignore_unused_variable_warning( c2 );
2015 FEELPP_ASSERT( 0 ).error(
"divergence of a scalar function is undefined." );
2019 value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
2021 detail::ignore_unused_variable_warning( c1 );
2022 detail::ignore_unused_variable_warning( c2 );
2023 return M_div[i][q]( 0,0 );
2025 value_type res = value_type( 0 );
2027 for (
int ii = 0; ii < nDim; ++ii )
2028 res += M_grad[i][ii]( ii,q );
2042 value_type
curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2044 detail::ignore_unused_variable_warning( c1 );
2045 detail::ignore_unused_variable_warning( c2 );
2046 return curl( i, c1, c2, q, mpl::int_<rank>() );
2049 value_type
curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2051 detail::ignore_unused_variable_warning( i );
2052 detail::ignore_unused_variable_warning( q );
2053 detail::ignore_unused_variable_warning( c1 );
2054 detail::ignore_unused_variable_warning( c2 );
2055 throw std::logic_error(
"invalid use of curl operator, field must be vectorial" );
2059 value_type
curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
2061 detail::ignore_unused_variable_warning( c2 );
2062 return M_curl[i][q]( c1 );
2064 value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2066 detail::ignore_unused_variable_warning( c1 );
2067 detail::ignore_unused_variable_warning( c2 );
2068 return curlx( i, c1, c2, q, mpl::int_<rank>() );
2070 value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2072 detail::ignore_unused_variable_warning(i);
2073 detail::ignore_unused_variable_warning(q);
2074 detail::ignore_unused_variable_warning(c1);
2075 detail::ignore_unused_variable_warning(c2);
2076 throw std::logic_error(
"invalid use of curl operator, field must be vectorial");
2079 value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
2081 detail::ignore_unused_variable_warning( c1 );
2082 detail::ignore_unused_variable_warning( c2 );
2083 return M_curl[i][q]( 0 );
2085 value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2087 detail::ignore_unused_variable_warning( c1 );
2088 detail::ignore_unused_variable_warning( c2 );
2089 return curly( i, c1, c2, q, mpl::int_<rank>() );
2091 value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2093 detail::ignore_unused_variable_warning(i);
2094 detail::ignore_unused_variable_warning(q);
2095 detail::ignore_unused_variable_warning(c1);
2096 detail::ignore_unused_variable_warning(c2);
2097 throw std::logic_error(
"invalid use of curl operator, field must be vectorial");
2100 value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
2102 detail::ignore_unused_variable_warning( c1 );
2103 detail::ignore_unused_variable_warning( c2 );
2104 return M_curl[i][q]( 1 );
2107 value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2109 detail::ignore_unused_variable_warning( c1 );
2110 detail::ignore_unused_variable_warning( c2 );
2111 return curlz( i, c1, c2, q, mpl::int_<rank>() );
2113 value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2115 detail::ignore_unused_variable_warning(i);
2116 detail::ignore_unused_variable_warning(q);
2117 detail::ignore_unused_variable_warning(c1);
2118 detail::ignore_unused_variable_warning(c2);
2119 throw std::logic_error(
"invalid use of curl operator, field must be vectorial");
2122 value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> )
const
2124 detail::ignore_unused_variable_warning( c1 );
2125 detail::ignore_unused_variable_warning( c2 );
2126 return M_curl[i][q]( 2 );
2129 value_type hess( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q )
const
2131 return hess( i, c1, c2, q, mpl::int_<rank>() );
2133 value_type hess( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> )
const
2135 return M_hessian[i][q]( c1,c2 );
2140 Context( Context
const& ) {}
2143 boost::optional<precompute_ptrtype> M_pc;
2144 uint16_type M_npoints;
2148 reference_element_ptrtype M_ref_ele;
2150 geometric_mapping_context_ptrtype M_gmc;
2152 boost::multi_array<id_type,2> M_phi;
2153 boost::multi_array<ref_grad_type,2> M_gradphi;
2154 boost::multi_array<hess_type,2> M_hessphi;
2155 boost::multi_array<dn_type,2> M_dn;
2156 boost::multi_array<grad_type,2> M_grad;
2157 boost::multi_array<dx_type,2> M_dx;
2158 boost::multi_array<dy_type,2> M_dy;
2159 boost::multi_array<dz_type,2> M_dz;
2160 boost::multi_array<div_type,2> M_div;
2161 boost::multi_array<curl_type,2> M_curl;
2162 boost::multi_array<hess_type,2> M_hessian;
2165 template<
size_type context_v,
size_type context_g,
typename BasisType,
typename GeoType,
typename ElementType>
2166 boost::shared_ptr<Context<context_v,BasisType, GeoType, ElementType> >
2167 context( boost::shared_ptr<BasisType> b, boost::shared_ptr<GeoType> gm, precompute_ptrtype& pc )
2169 return boost::shared_ptr<Context<context_v,BasisType, GeoType, ElementType> >(
2170 new Context<context_v, BasisType, GeoType, ElementType>( context_v,
2176 template<
int contextv,
int contextg,
typename BasisType,
typename GeoType,
typename ElementType>
2177 boost::shared_ptr<Context<contextv,BasisType, GeoType, ElementType> >
2178 ctx( boost::shared_ptr<BasisType>
const& b,
2179 boost::shared_ptr<
typename GeoType::template Context<contextg, ElementType> >
const& gm,
2180 precompute_ptrtype pc, ElementType& e )
2182 typedef Context<contextv,BasisType, GeoType, ElementType> ctx_type;
2183 return boost::shared_ptr<ctx_type>(
new ctx_type( b, gm, pc ) );
2186 template<
int contextv,
typename BasisType,
typename GeoType,
typename ElementType>
2187 boost::shared_ptr<Context<contextv,BasisType, GeoType, ElementType> >
2188 ctx( boost::shared_ptr<BasisType>
const& b,
2189 boost::shared_ptr<
typename GeoType::template Context<contextv, ElementType> >
const& gm,
2190 precompute_ptrtype pc, ElementType& e )
2192 typedef Context<contextv,BasisType, GeoType, ElementType> ctx_type;
2193 return boost::shared_ptr<ctx_type>(
new ctx_type( b, gm, pc ) );
2204 matrix_type M_coeff;
2205 std::string M_fname;
2208 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const bool PolynomialSet<Poly,PolySetType>::is_scalar;
2209 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const bool PolynomialSet<Poly,PolySetType>::is_vectorial;
2210 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const bool PolynomialSet<Poly,PolySetType>::is_tensor2;
2211 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const uint16_type PolynomialSet<Poly,PolySetType>::nComponents;
2212 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const uint16_type PolynomialSet<Poly,PolySetType>::nComponents1;
2213 template<
typename Poly,
template<u
int16_type>
class PolySetType>
const uint16_type PolynomialSet<Poly,PolySetType>::nComponents2;
void insert(PolynomialSet< Poly, PolySetType > const &p, bool erase=false)
Definition: polynomialset.hpp:673
self_type derivate(uint16_type l) const
Derivate with respect to the l-th direction.
Definition: polynomialset.hpp:554
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
basis_type const & basis() const
Definition: polynomialset.hpp:253
size_type polynomialDimension() const
Definition: polynomialset.hpp:285
Polynomial< Poly, Type > dx(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:63
PolynomialSet< Poly, Tensor2 > hessian(PolynomialSet< Poly, Scalar > const &p)
compute the gradient of a vectorial polynomial set
Definition: operations.hpp:231
polynomial class
Definition: polynomial.hpp:67
PolynomialSet< Poly, PolySetType > polynomials(std::vector< uint16_type > const &list_p) const
Definition: polynomialset.hpp:378
matrix_type const & coeff() const
Definition: polynomialset.hpp:245
matrix_type const & d(uint16_type i) const
derivatives of Dubiner polynomials the derivatives are computed at the nodes of the lattice ...
Definition: polynomialset.hpp:539
ublas::matrix< value_type > evaluate(ublas::vector_expression< AE > const &__pt) const
Definition: polynomialset.hpp:494
ublas::vector< value_type > evaluate(uint16_type i, ublas::vector_expression< AE > const &__pt) const
Definition: polynomialset.hpp:484
uint16_type degree() const
Definition: polynomialset.hpp:236
Polynomial< Poly, Type > dz(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:109
size_type polynomialDimensionPerComponent() const
Definition: polynomialset.hpp:293
Polynomial< Poly, PolySetType > polynomial(uint16_type i) const
Definition: polynomialset.hpp:459
component_type operator[](uint16_type i) const
extract the i-th component of a vectorial polynomial set
Definition: polynomialset.hpp:216
uint16_type nbDof() const
Definition: polynomialset.hpp:664
Polynomial< Poly, Vectorial > curl(Polynomial< Poly, Vectorial > const &p)
compute the curl of a vectorial polynomial p
Definition: operations.hpp:305
PolynomialSet< Poly, PolySetType > polynomialsUpToDimension(int dim_p) const
Definition: polynomialset.hpp:407
a Set of polynomials
Definition: fe.hpp:37
static bool isVectorial()
Definition: polynomialset.hpp:269
gradient_polynomialset_type gradient() const
Definition: polynomialset.hpp:600
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
PreCompute precompute_type
Definition: polynomialset.hpp:993
void setCoefficient(matrix_type const &__c, bool __as_is=false)
Definition: polynomialset.hpp:342
PolynomialSet< Poly, PolySetType > polynomialsRange(uint16_type dim_bot, uint16_type dim_top) const
Definition: polynomialset.hpp:433
Elements & operator=(Elements const &e)
Definition: elements.hpp:335
bool isZero() const
Definition: polynomialset.hpp:301
std::string name(std::string sep=".") const
Definition: polynomialset.hpp:326
matrix_type evaluate(ublas::matrix_expression< AE > const &__pts) const
Definition: polynomialset.hpp:511
Polynomial< Poly, Type > dy(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:85
virtual std::string familyName() const
Definition: polynomialset.hpp:310
static bool isScalar()
Definition: polynomialset.hpp:261