30 #ifndef __LinearForm_H
31 #define __LinearForm_H 1
33 #include <Eigen/Eigen>
34 #include <Eigen/StdVector>
36 #include <boost/fusion/support/pair.hpp>
37 #include <boost/fusion/container.hpp>
38 #include <boost/fusion/sequence.hpp>
39 #include <boost/fusion/algorithm.hpp>
40 #include <boost/multi_array.hpp>
42 #include <feel/feelalg/vectorvalue.hpp>
48 namespace fusion = boost::fusion;
73 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
82 enum { nDim = SpaceType::nDim };
84 typedef LinearForm<SpaceType, VectorType, ElemContType> self_type;
86 typedef SpaceType space_type;
87 typedef boost::shared_ptr<SpaceType> space_ptrtype;
89 typedef typename space_type::value_type value_type;
90 typedef typename space_type::real_type real_type;
91 typedef VectorType vector_type;
92 typedef boost::shared_ptr<VectorType> vector_ptrtype;
93 typedef typename space_type::template Element<value_type, ElemContType> element_type;
96 typedef typename space_type::component_fespace_type component_fespace_type;
97 typedef typename space_type::element_type::component_type component_type;
99 typedef typename space_type::mesh_type mesh_type;
100 typedef typename mesh_type::element_type mesh_test_element_type;
101 typedef typename mesh_type::element_type::permutation_type permutation_type;
102 typedef typename space_type::gm_type gm_type;
103 typedef typename space_type::gm_ptrtype gm_ptrtype;
104 typedef typename space_type::gm1_type gm1_type;
105 typedef typename space_type::gm1_ptrtype gm1_ptrtype;
107 typedef typename space_type::fe_type fe_type;
108 typedef typename space_type::basis_0_type::precompute_type test_precompute_type;
110 typedef boost::shared_ptr<test_precompute_type> test_precompute_ptrtype;
130 template<
typename GeomapContext,
133 typename GeomapExprContext = GeomapContext
137 typedef FormContextBase<GeomapContext, IM, GeomapExprContext> super;
140 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
142 typedef Context<GeomapContext,ExprT,IM,GeomapExprContext> form_context_type;
143 typedef LinearForm<SpaceType,VectorType, ElemContType> form_type;
144 typedef typename SpaceType::dof_type dof_type;
145 typedef typename form_type::value_type value_type;
149 typedef GeomapContext map_geometric_mapping_context_type;
150 typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type geometric_mapping_context_ptrtype;
151 typedef typename geometric_mapping_context_ptrtype::element_type geometric_mapping_context_type;
152 typedef typename geometric_mapping_context_type::gm_type geometric_mapping_type;
154 typedef typename super::map_geometric_mapping_context_type map_test_geometric_mapping_context_type;
155 typedef typename super::geometric_mapping_context_ptrtype test_geometric_mapping_context_ptrtype;
156 typedef typename super::geometric_mapping_context_type test_geometric_mapping_context_type;
157 typedef typename super::geometric_mapping_type test_geometric_mapping_type;
159 typedef typename super::map_geometric_mapping_context_type map_trial_geometric_mapping_context_type;
160 typedef typename super::geometric_mapping_context_ptrtype trial_geometric_mapping_context_ptrtype;
161 typedef typename super::geometric_mapping_context_type trial_geometric_mapping_context_type;
162 typedef typename super::geometric_mapping_type trial_geometric_mapping_type;
166 static const uint16_type nDim = test_geometric_mapping_type::nDim;
168 typedef ExprT expression_type;
170 typedef typename space_type::mesh_type mesh_type;
171 typedef typename mesh_type::element_type mesh_element_type;
172 typedef typename mesh_element_type::permutation_type permutation_type;
173 typedef typename space_type::fe_type test_fe_type;
174 typedef typename space_type::fe_type trial_fe_type;
175 typedef boost::shared_ptr<test_fe_type> test_fe_ptrtype;
176 typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
178 test_geometric_mapping_type,
179 mesh_element_type> test_fecontext_type;
180 typedef test_fecontext_type trial_fecontext_type;
183 typedef mpl::int_<fusion::result_of::template size<GeomapContext>::type::value> map_size;
185 typedef typename super::map_size map_size;
188 typedef boost::shared_ptr<test_fecontext_type> test_fecontext_ptrtype;
191 typedef typename mpl::if_<mpl::equal_to<map_size,mpl::int_<2> >,
195 typedef typename super::gmc1 gmc1;
199 typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type left_gmc_ptrtype;
200 typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type::element_type left_gmc_type;
201 typedef typename fusion::result_of::value_at_key<GeomapContext,gmc1 >::type right_gmc_ptrtype;
202 typedef typename fusion::result_of::value_at_key<GeomapContext,gmc1 >::type::element_type right_gmc_type;
204 typedef fusion::map<fusion::pair<gmc<0>, left_gmc_ptrtype> > map_left_gmc_type;
205 typedef fusion::map<fusion::pair<gmc<0>, right_gmc_ptrtype> > map_right_gmc_type;
207 typedef typename super::left_gmc_ptrtype left_gmc_ptrtype;
208 typedef typename super::left_gmc_type left_gmc_type;
209 typedef typename super::right_gmc_ptrtype right_gmc_ptrtype;
210 typedef typename super::right_gmc_type right_gmc_type;
212 typedef typename super::map_left_gmc_type map_left_gmc_type;
213 typedef typename super::map_right_gmc_type map_right_gmc_type;
216 typedef typename mpl::if_<mpl::equal_to<map_size, mpl::int_<1> >,
217 mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > >,
218 mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype>,
219 fusion::pair<gmc<1>, test_fecontext_ptrtype> > > >::type::type map_test_fecontext_type;
221 typedef fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > map_left_test_fecontext_type;
222 typedef fusion::map<fusion::pair<gmc<1>, test_fecontext_ptrtype> > map_right_test_fecontext_type;
224 typedef typename super::map_geometric_mapping_expr_context_type map_geometric_mapping_expr_context_type;
226 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type, map_left_test_fecontext_type> eval0_expr_type;
227 typedef boost::shared_ptr<eval0_expr_type> eval0_expr_ptrtype;
229 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type, map_right_test_fecontext_type> eval1_expr_type;
230 typedef boost::shared_ptr<eval1_expr_type> eval1_expr_ptrtype;
234 typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
236 test_geometric_mapping_type,
237 mesh_element_type>::template Index<> test_index_type;
244 static const int rep_shape = 2;
246 typedef typename space_type::dof_type test_dof_type;
247 static const int nDofPerElementTest = space_type::dof_type::nDofPerElement;
248 typedef Eigen::Matrix<value_type, nDofPerElementTest, 1> local_vector_type;
249 typedef Eigen::Matrix<value_type, 2*nDofPerElementTest, 1> local2_vector_type;
250 typedef Eigen::Matrix<int, nDofPerElementTest, 1> local_row_type;
251 typedef Eigen::Matrix<int, 2*nDofPerElementTest, 1> local2_row_type;
255 Context( form_type& __form,
256 map_test_geometric_mapping_context_type
const& _gmcTest,
257 map_trial_geometric_mapping_context_type
const & _gmcTrial,
258 map_geometric_mapping_expr_context_type
const& _gmcExpr,
262 template<
typename IM2>
263 Context( form_type& __form,
264 map_test_geometric_mapping_context_type
const& _gmcTest,
265 map_trial_geometric_mapping_context_type
const & _gmcTrial,
266 map_geometric_mapping_expr_context_type
const& _gmcExpr,
271 template<
typename IM2>
272 Context( form_type& __form,
273 map_test_geometric_mapping_context_type
const& _gmcTest,
274 map_trial_geometric_mapping_context_type
const & _gmcTrial,
275 map_geometric_mapping_expr_context_type
const& _gmcExpr,
285 bool isZero(
typename mesh_type::element_iterator it )
const
287 return this->isZero( it->id() );
289 bool isZero(
typename mesh_type::element_type
const& e )
const
291 return this->isZero( e.id() );
294 void update( map_test_geometric_mapping_context_type
const& _gmcTest,
295 map_trial_geometric_mapping_context_type
const & gmcTrial,
296 map_geometric_mapping_expr_context_type
const& _gmcExpr );
299 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
300 map_trial_geometric_mapping_context_type
const & gmcTrial,
301 map_geometric_mapping_expr_context_type
const& gmcExpr,
302 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad );
304 void update( map_test_geometric_mapping_context_type
const& _gmcTest,
305 map_trial_geometric_mapping_context_type
const & gmcTrial,
306 map_geometric_mapping_expr_context_type
const& _gmcExpr,
310 void update( map_test_geometric_mapping_context_type
const& _gmcTest,
311 map_trial_geometric_mapping_context_type
const & gmcTrial,
312 map_geometric_mapping_expr_context_type
const& _gmcExpr,
315 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
316 map_trial_geometric_mapping_context_type
const & gmcTrial,
317 map_geometric_mapping_expr_context_type
const& gmcExpr,
319 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad )
322 this->updateInCaseOfInterpolate( gmcTest, gmcTrial, gmcExpr, indexLocalToQuad );
326 void update( map_test_geometric_mapping_context_type
const& _gmcTest,
327 map_trial_geometric_mapping_context_type
const & gmcTrial,
328 map_geometric_mapping_expr_context_type
const& _gmcExpr,
329 IM
const& im, mpl::int_<2> );
334 integrate( mpl::int_<fusion::result_of::size<GeomapContext>::type::value>() );
337 void integrateInCaseOfInterpolate( std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad,
338 bool isFirstExperience )
340 integrateInCaseOfInterpolate( mpl::int_<fusion::result_of::size<GeomapContext>::type::value>(),
341 indexLocalToQuad, isFirstExperience );
347 assemble( M_gmc_left->id() );
351 void assemble( mpl::int_<2> )
353 assemble( M_gmc_left->id(), M_gmc_right->id() );
361 template<
typename Pts>
362 void precomputeBasisAtPoints( Pts
const& pts )
364 M_test_pc = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), pts ) );
372 template<
typename Pts>
373 void precomputeBasisAtPoints( uint16_type __f, permutation_type
const& __p, Pts
const& pts )
375 M_test_pc_face[__f][__p] = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), pts ) );
385 test_precompute_ptrtype
const& testPc( uint16_type __f,
386 permutation_type __p = permutation_type( permutation_type::NO_PERMUTATION ) )
const
392 return M_test_pc_face.find( __f )->second.find( __p )->second;
395 template<
typename PtsSet>
396 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
397 precomputeTestBasisAtPoints( PtsSet
const& pts )
399 typedef typename boost::is_same< permutation_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
400 return precomputeTestBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
404 template<
typename PtsSet>
405 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
406 precomputeTestBasisAtPoints( PtsSet
const& pts, mpl::bool_<false> )
408 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
412 template<
typename PtsSet>
413 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
414 precomputeTestBasisAtPoints( PtsSet
const& pts, mpl::bool_<true> )
416 QuadMapped<PtsSet> qm;
417 typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
418 typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
420 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
422 for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
424 for ( permutation_type __p( permutation_type::IDENTITY );
425 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
427 testpc[__f][__p] = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), ppts[__f].find( __p )->second ) );
441 void integrateInCaseOfInterpolate( mpl::int_<1>,
442 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad,
443 bool isFirstExperience );
448 dof_type* M_test_dof;
449 const list_block_type& M_lb;
451 test_precompute_ptrtype M_test_pc;
452 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > M_test_pc_face;
454 map_test_geometric_mapping_context_type M_gmc;
455 left_gmc_ptrtype M_gmc_left;
456 right_gmc_ptrtype M_gmc_right;
457 map_left_gmc_type M_left_map;
458 map_right_gmc_type M_right_map;
459 map_test_fecontext_type M_test_fec;
460 map_left_test_fecontext_type M_test_fec0;
461 map_right_test_fecontext_type M_test_fec1;
463 local_vector_type M_rep;
464 local2_vector_type M_rep_2;
465 local_row_type M_local_rows;
466 local2_row_type M_local_rows_2;
467 local_row_type M_local_rowsigns;
468 local2_row_type M_local_rowsigns_2;
470 eval0_expr_ptrtype M_eval0_expr;
471 eval1_expr_ptrtype M_eval1_expr;
484 LinearForm( LinearForm
const & __vf );
486 LinearForm( space_ptrtype
const& __X,
490 bool do_threshold =
false,
491 value_type threshold = type_traits<value_type>::epsilon() );
492 LinearForm( space_ptrtype
const& __X,
494 list_block_type
const& __lb,
497 bool do_threshold =
false,
498 value_type threshold = type_traits<value_type>::epsilon() );
517 template <
class ExprT>
518 LinearForm&
operator=( Expr<ExprT>
const& expr );
527 template <
class ExprT>
528 LinearForm& operator+=( Expr<ExprT>
const& expr );
535 LinearForm<component_fespace_type,
537 typename component_type::container_type>
538 operator()( component_type& __u )
540 typedef typename component_type::container_type container_type;
542 typedef LinearForm<component_fespace_type, vector_type,
543 container_type> vflf_type;
545 list_block_type __list_block;
546 __list_block.push_back( Block( 0, 0,
548 return vflf_type( __u, M_F, __list_block,
false );
556 value_type operator()(
size_type i )
const
567 value_type operator()( element_type
const& __v )
const
569 return M_F->dot( __v );
578 value_type operator()(
typename space_type::element_type
const& __v )
const
580 return M_F->dot( __v );
593 space_ptrtype
const& functionSpace()
const
601 space_ptrtype
const& testSpace()
const
609 gm_ptrtype
const& gm()
const
617 gm1_ptrtype
const& gm1()
const
628 test_precompute_ptrtype
const& testPc( uint16_type __f,
629 permutation_type __p = permutation_type( permutation_type::NO_PERMUTATION ) )
const
634 return M_test_pc_face.find( __f )->second.find( __p )->second;
637 vector_type& representation()
const
641 vector_ptrtype vectorPtr()
const
645 vector_ptrtype vectorPtr()
650 list_block_type
const& blockList()
const
657 return M_row_startInVector;
664 bool doThreshold( value_type
const& v )
const
666 return ( math::abs( v ) > M_threshold );
679 void setFunctionSpace( space_ptrtype
const& __X )
688 void setThreshold( value_type eps )
697 void setDoThreshold(
bool do_threshold )
699 M_do_threshold = do_threshold;
712 template<
typename Pts>
713 void precomputeBasisAtPoints( Pts
const& pts )
715 M_test_pc = test_precompute_ptrtype(
new test_precompute_type( functionSpace()->fe(), pts ) );
723 template<
typename Pts>
724 void precomputeBasisAtPoints( uint16_type __f, permutation_type __p, Pts
const& pts )
726 M_test_pc_face[__f][__p] = test_precompute_ptrtype(
new test_precompute_type( functionSpace()->fe(), pts ) );
733 void add(
size_type i, value_type
const& v )
735 if ( M_do_threshold )
737 if ( doThreshold( v ) )
738 M_F->add( i+this->rowStartInVector(), v );
742 M_F->add( i+this->rowStartInVector(), v );
749 void addVector(
int* i,
int n, value_type* v )
751 if ( this->rowStartInVector()!=0 )
752 for (
int k = 0; k< n ; ++k )
753 i[k]+=this->rowStartInVector();
755 M_F->addVector( i, n, v );
762 void set(
size_type i, value_type
const& v )
770 void zero() { M_F->zero(); }
772 LinearForm& operator+=( LinearForm& f )
788 template <
class ExprT>
void assign( Expr<ExprT>
const& expr,
bool init, mpl::bool_<false> );
789 template <
class ExprT>
void assign( Expr<ExprT>
const& expr,
bool init, mpl::bool_<true> );
796 list_block_type M_lb;
800 test_precompute_ptrtype M_test_pc;
802 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > M_test_pc_face;
805 value_type M_threshold;
808 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
809 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( LinearForm
const & __vf )
814 M_row_startInVector( __vf.M_row_startInVector ),
817 M_do_threshold( __vf.M_do_threshold ),
818 M_threshold( __vf.M_threshold )
821 DVLOG(2) <<
"LinearForm copy constructor\n";
822 DVLOG(2) <<
" n Dof : " << M_X->nDof() <<
"\n";
823 DVLOG(2) <<
" F size : " << M_F->size() <<
"\n";
824 DVLOG(2) <<
"block size : " << M_lb.size() <<
"\n";
827 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
828 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( space_ptrtype
const& __X,
833 value_type threshold )
838 M_row_startInVector( rowstart ),
839 M_do_threshold( do_threshold ),
840 M_threshold( threshold )
843 for ( uint16_type __i = 0; __i < M_X->qDim(); ++__i )
845 M_lb.push_back( Block( __i, 0,
846 __i*M_X->nDofPerComponent(),
848 DVLOG(2) <<
"[linearform::linearform] block: "
849 << Block( __i, 0, __i*M_X->nDofPerComponent(), 0 ) <<
"\n";
856 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
857 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( space_ptrtype
const& __X,
859 list_block_type
const& __lb,
863 value_type threshold )
868 M_row_startInVector( rowstart ),
869 M_do_threshold( do_threshold ),
870 M_threshold( threshold )
876 template<
typename LFType,
typename TheSpaceType,
typename ExprType>
879 LFAssign( LFAssign
const& lfa )
883 M_expr( lfa.M_expr ),
884 M_index( lfa.M_index ),
887 LFAssign( LFType& lf, TheSpaceType
const& Xh, ExprType
const& expr,
bool init )
895 template<
typename SpaceType>
896 void operator()( boost::shared_ptr<SpaceType>
const& X )
const
898 if ( M_lf.testSpace()->worldsComm()[M_index].isActive() )
900 DVLOG(2) <<
"expression has test functions ? :"
901 << ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result
904 if ( !ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result )
910 list_block_type __list_block;
912 if ( M_lf.testSpace()->worldsComm()[M_index].globalSize()>1 )
914 if (M_lf.testSpace()->hasEntriesForAllSpaces())
915 __list_block.push_back( Block( 0, 0, M_Xh->nLocalDofStart( M_index ), 0 ) );
917 __list_block.push_back( Block( 0, 0, 0, 0 ) );
920 __list_block.push_back( Block( 0, 0, M_Xh->nDofStart( M_index ), 0 ) );
922 LinearForm<SpaceType,typename LFType::vector_type, typename LFType::element_type> lf( X,
925 M_lf.rowStartInVector(),
951 if ( M_init ) M_lf.representation().zero();
958 TheSpaceType
const& M_Xh;
959 ExprType
const& M_expr;
966 template<
typename LFType,
typename TheSpaceType,
typename ExprType>
967 LFAssign<LFType,TheSpaceType,ExprType>
968 make_lfassign( LFType& lf, TheSpaceType
const& Xh, ExprType
const& expr,
bool init )
970 return LFAssign<LFType,TheSpaceType,ExprType>( lf, Xh, expr, init );
975 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
976 template<
typename ExprT>
978 LinearForm<SpaceType, VectorType, ElemContType>::assign( Expr<ExprT>
const& __expr,
bool init, mpl::bool_<true> )
980 fusion::for_each( M_X->functionSpaces(), make_lfassign( *
this, M_X, __expr, init ) );
982 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
983 template<
typename ExprT>
985 LinearForm<SpaceType, VectorType, ElemContType>::assign( Expr<ExprT>
const& __expr,
bool init, mpl::bool_<false> )
989 M_lb.push_back( Block( 0, 0, 0, 0 ) );
994 typename list_block_type::const_iterator __bit = M_lb.begin();
995 typename list_block_type::const_iterator __ben = M_lb.end();
997 for ( ; __bit != __ben; ++__bit )
999 DVLOG(2) <<
"LinearForm:: block: " << *__bit <<
"\n";
1000 size_type g_ic_start = __bit->globalRowStart();
1001 DVLOG(2) <<
"LinearForm:: g_ic_start: " << g_ic_start <<
"\n";
1003 M_F->zero( g_ic_start,g_ic_start + M_X->nDof() );
1007 __expr.assemble( M_X, *
this );
1012 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
1013 template<
typename ExprT>
1014 LinearForm<SpaceType, VectorType, ElemContType>&
1019 this->assign( __expr,
true, mpl::bool_<( SpaceType::nSpaces > 1 )>() );
1022 template<
typename SpaceType,
typename VectorType,
typename ElemContType>
1023 template<
typename ExprT>
1024 LinearForm<SpaceType, VectorType, ElemContType>&
1025 LinearForm<SpaceType, VectorType, ElemContType>::operator+=( Expr<ExprT>
const& __expr )
1027 this->assign( __expr,
false, mpl::bool_<( SpaceType::nSpaces > 1 )>() );
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Elements & operator=(Elements const &e)
Definition: elements.hpp:335