30 #ifndef __BilinearForm_H
31 #define __BilinearForm_H 1
33 #include <Eigen/Eigen>
34 #include <Eigen/StdVector>
39 #include <boost/parameter.hpp>
40 #include <boost/fusion/support/pair.hpp>
41 #include <boost/fusion/container.hpp>
42 #include <boost/fusion/sequence.hpp>
43 #include <boost/fusion/algorithm.hpp>
44 #include <boost/spirit/home/phoenix.hpp>
45 #include <boost/spirit/home/phoenix/core/argument.hpp>
46 #include <feel/feelcore/context.hpp>
47 #include <feel/feelalg/matrixvalue.hpp>
50 #include <feel/feeldiscr/stencil.hpp>
58 namespace parameter = boost::parameter;
59 namespace fusion = boost::fusion;
64 template<
typename FE1,
typename FE2,
typename ElemContType>
class BilinearForm;
69 template<
typename BFType,
typename ExprType>
72 typedef typename BFType::matrix_type matrix_type;
73 BFAssign2( BFAssign2
const& lfa )
77 M_test_index( lfa.M_test_index )
79 BFAssign2( BFType& lf, ExprType
const& expr )
85 template<
typename SpaceType>
86 void operator()( boost::shared_ptr<SpaceType>
const& X )
const
88 DVLOG(2) <<
"[BFAssign2::operator()] start loop on trial spaces against test space index: " << M_test_index <<
"\n";
90 if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
92 fusion::for_each( M_bf.trialSpace()->functionSpaces(),
93 make_bfassign1( M_bf, M_expr, X, M_test_index ) );
96 DVLOG(2) <<
"[BFAssign2::operator()] stop loop on trial spaces against test space index: " << M_test_index <<
"\n";
102 ExprType
const& M_expr;
105 template<
typename BFType,
typename ExprType>
106 BFAssign2<BFType,ExprType>
107 make_bfassign2( BFType& lf, ExprType
const& expr )
109 return BFAssign2<BFType,ExprType>( lf, expr );
112 template<
typename BFType,
typename ExprType,
typename TestSpaceType>
115 typedef typename BFType::matrix_type matrix_type;
116 BFAssign1( BFAssign1
const& lfa )
119 M_test( lfa.M_test ),
120 M_expr( lfa.M_expr ),
121 M_trial_index( lfa.M_trial_index )
123 BFAssign1( BFType& lf,
124 ExprType
const& expr,
125 boost::shared_ptr<TestSpaceType>
const& Testh,
132 M_test_index( test_index )
134 template<
typename SpaceType>
135 void operator()( boost::shared_ptr<SpaceType>
const& trial )
const;
139 boost::shared_ptr<TestSpaceType> M_test;
140 ExprType
const& M_expr;
144 template<
typename BFType,
typename ExprType,
typename TestSpaceType>
145 BFAssign1<BFType,ExprType,TestSpaceType>
146 make_bfassign1( BFType& lf,
147 ExprType
const& expr,
148 boost::shared_ptr<TestSpaceType>
const& test_space,
151 return BFAssign1<BFType,ExprType,TestSpaceType>( lf, expr, test_space, test_index );
155 template<
typename BFType,
typename ExprType,
typename TrialSpaceType>
158 typedef typename BFType::matrix_type matrix_type;
159 BFAssign3( BFAssign3
const& lfa )
162 M_trial( lfa.M_trial ),
163 M_expr( lfa.M_expr ),
164 M_test_index( lfa.M_test_index )
166 BFAssign3( BFType& lf,
167 ExprType
const& expr,
168 boost::shared_ptr<TrialSpaceType>
const& Trialh,
174 M_trial_index( trial_index ),
177 template<
typename SpaceType>
178 void operator()( boost::shared_ptr<SpaceType>
const& test )
const;
182 boost::shared_ptr<TrialSpaceType> M_trial;
183 ExprType
const& M_expr;
187 template<
typename BFType,
typename ExprType,
typename TrialSpaceType>
188 BFAssign3<BFType,ExprType,TrialSpaceType>
189 make_bfassign3( BFType& lf,
190 ExprType
const& expr,
191 boost::shared_ptr<TrialSpaceType>
const& trial_space,
194 return BFAssign3<BFType,ExprType,TrialSpaceType>( lf, expr, trial_space, trial_index );
205 template<
typename FE1,
207 typename ElemContType = VectorUblas<typename FE1::value_type> >
215 enum { nDim = FE1::nDim };
217 typedef FE1 space_1_type;
218 typedef boost::shared_ptr<FE1> space_1_ptrtype;
219 typedef space_1_type test_space_type;
220 typedef boost::shared_ptr<space_1_type> test_space_ptrtype;
222 typedef FE2 space_2_type;
223 typedef boost::shared_ptr<FE2> space_2_ptrtype;
224 typedef space_2_type trial_space_type;
225 typedef boost::shared_ptr<space_2_type> trial_space_ptrtype;
227 typedef typename FE1::value_type value_type;
228 typedef typename space_1_type::template Element<value_type,ElemContType> element_1_type;
230 typedef typename space_2_type::template Element<value_type,ElemContType> element_2_type;
231 typedef BilinearForm<FE1, FE2, ElemContType> self_type;
234 typedef typename space_1_type::component_fespace_type component_space_1_type;
235 typedef typename element_1_type::component_type component_1_type;
236 typedef typename space_2_type::component_fespace_type component_space_2_type;
237 typedef typename element_2_type::component_type component_2_type;
240 typedef typename space_1_type::mesh_type mesh_1_type;
241 typedef typename mesh_1_type::element_type mesh_element_1_type;
242 typedef typename mesh_element_1_type::permutation_type permutation_1_type;
244 typedef typename space_2_type::mesh_type mesh_2_type;
245 typedef typename mesh_2_type::element_type mesh_element_2_type;
246 typedef typename mesh_element_2_type::permutation_type permutation_2_type;
248 typedef typename space_1_type::fe_type fe_1_type;
249 typedef typename space_2_type::fe_type fe_2_type;
251 typedef typename space_1_type::gm_type gm_1_type;
252 typedef typename space_1_type::gm_ptrtype gm_1_ptrtype;
254 typedef typename space_1_type::gm1_type gm1_1_type;
255 typedef typename space_1_type::gm1_ptrtype gm1_1_ptrtype;
257 typedef typename space_2_type::gm_type gm_2_type;
258 typedef typename space_2_type::gm_ptrtype gm_2_ptrtype;
260 typedef typename space_2_type::gm1_type gm1_2_type;
261 typedef typename space_2_type::gm1_ptrtype gm1_2_ptrtype;
264 typedef MatrixSparse<value_type> matrix_type;
265 typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
266 static const bool is_row_major =
true;
268 typedef typename mpl::if_<mpl::equal_to<mpl::bool_<is_row_major>, mpl::bool_<true> >,
269 mpl::identity<ublas::row_major>,
270 mpl::identity<ublas::column_major> >::type::type layout_type;
273 struct test_precompute
276 typedef typename space_1_type::basis_0_type::PreCompute type;
277 typedef boost::shared_ptr<type> ptrtype;
280 struct trial_precompute
283 typedef typename space_2_type::basis_0_type::PreCompute type;
284 typedef boost::shared_ptr<type> ptrtype;
307 template<
typename GeomapTestContext,
310 typename GeomapExprContext = GeomapTestContext,
311 typename GeomapTrialContext = GeomapTestContext
315 typedef FormContextBase<GeomapTestContext,IM,GeomapExprContext> super;
318 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
321 typedef Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext> form_context_type;
322 typedef BilinearForm<FE1, FE2, ElemContType> form_type;
323 typedef typename FE1::dof_type dof_1_type;
324 typedef typename FE2::dof_type dof_2_type;
326 typedef typename form_type::value_type value_type;
329 typedef typename super::map_geometric_mapping_context_type map_test_geometric_mapping_context_type;
330 typedef typename super::geometric_mapping_context_ptrtype test_geometric_mapping_context_ptrtype;
331 typedef typename super::geometric_mapping_context_type test_geometric_mapping_context_type;
332 typedef typename super::geometric_mapping_type test_geometric_mapping_type;
333 typedef typename super::map_size test_map_size;
334 typedef typename super::gmc1 test_gmc1;
335 typedef typename super::left_gmc_ptrtype test_left_gmc_ptrtype;
336 typedef typename super::right_gmc_ptrtype test_right_gmc_ptrtype;
337 typedef typename super::map_left_gmc_type test_map_left_gmc_type;
338 typedef typename super::map_right_gmc_type test_map_right_gmc_type;
341 typedef FormContextBase<GeomapTrialContext,IM,GeomapExprContext> super2;
343 typedef typename super2::map_geometric_mapping_context_type map_trial_geometric_mapping_context_type;
344 typedef typename super2::geometric_mapping_context_ptrtype trial_geometric_mapping_context_ptrtype;
345 typedef typename super2::geometric_mapping_context_type trial_geometric_mapping_context_type;
346 typedef typename super2::geometric_mapping_type trial_geometric_mapping_type;
347 typedef typename super2::map_size trial_map_size;
348 typedef typename super2::gmc1 trial_gmc1;
349 typedef typename super2::left_gmc_ptrtype trial_left_gmc_ptrtype;
350 typedef typename super2::right_gmc_ptrtype trial_right_gmc_ptrtype;
351 typedef typename super2::map_left_gmc_type trial_map_left_gmc_type;
352 typedef typename super2::map_right_gmc_type trial_map_right_gmc_type;
358 static const uint16_type nDim = test_geometric_mapping_type::nDim;
360 static const uint16_type nDimTest = test_space_type::mesh_type::nDim;
361 static const uint16_type nDimTrial = trial_space_type::mesh_type::nDim;
362 static const uint16_type nDimDiffBetweenTestTrial = ( nDimTest > nDimTrial )? nDimTest-nDimTrial : nDimTrial-nDimTest;
364 typedef ExprT expression_type;
366 typedef typename test_precompute<>::type test_precompute_type;
367 typedef typename test_precompute<>::ptrtype test_precompute_ptrtype;
368 typedef typename trial_precompute<>::type trial_precompute_type;
369 typedef typename trial_precompute<>::ptrtype trial_precompute_ptrtype;
371 typedef typename FE2::fe_type trial_fe_type;
372 typedef boost::shared_ptr<trial_fe_type> trial_fe_ptrtype;
373 typedef typename trial_fe_type::template Context< trial_geometric_mapping_context_type::context,
375 trial_geometric_mapping_type,
376 mesh_element_2_type> trial_fecontext_type;
377 typedef boost::shared_ptr<trial_fecontext_type> trial_fecontext_ptrtype;
378 typedef typename mpl::if_<mpl::equal_to<trial_map_size, mpl::int_<1> >,
379 mpl::identity<fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype> > >,
380 mpl::identity<fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype>,
381 fusion::pair<gmc<1>, trial_fecontext_ptrtype> > > >::type::type map_trial_fecontext_type;
383 typedef fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype> > map_left_trial_fecontext_type;
384 typedef fusion::map<fusion::pair<trial_gmc1, trial_fecontext_ptrtype> > map_right_trial_fecontext_type;
386 typedef typename FE1::fe_type test_fe_type;
387 typedef boost::shared_ptr<test_fe_type> test_fe_ptrtype;
388 typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
390 test_geometric_mapping_type,
391 mesh_element_1_type> test_fecontext_type;
392 typedef boost::shared_ptr<test_fecontext_type> test_fecontext_ptrtype;
393 typedef typename mpl::if_<mpl::equal_to<test_map_size, mpl::int_<1> >,
394 mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > >,
395 mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype>,
396 fusion::pair<gmc<1>, test_fecontext_ptrtype> > > >::type::type map_test_fecontext_type;
397 typedef fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > map_left_test_fecontext_type;
398 typedef fusion::map<fusion::pair<test_gmc1, test_fecontext_ptrtype> > map_right_test_fecontext_type;
402 typedef typename super::map_geometric_mapping_expr_context_type map_geometric_mapping_expr_context_type;
404 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
405 map_left_test_fecontext_type,
406 map_left_trial_fecontext_type> eval00_expr_type;
407 typedef boost::shared_ptr<eval00_expr_type> eval00_expr_ptrtype;
409 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
410 map_left_test_fecontext_type,
411 map_right_trial_fecontext_type> eval01_expr_type;
412 typedef boost::shared_ptr<eval01_expr_type> eval01_expr_ptrtype;
414 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
415 map_right_test_fecontext_type,
416 map_left_trial_fecontext_type> eval10_expr_type;
417 typedef boost::shared_ptr<eval10_expr_type> eval10_expr_ptrtype;
419 typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
420 map_right_test_fecontext_type,
421 map_right_trial_fecontext_type> eval11_expr_type;
422 typedef boost::shared_ptr<eval11_expr_type> eval11_expr_ptrtype;
425 static const int rep_shape = 4;
428 typedef typename FE1::dof_type test_dof_type;
429 typedef typename FE2::dof_type trial_dof_type;
430 static const int nDofPerElementTest = FE1::dof_type::nDofPerElement;
431 static const int nDofPerElementTrial = FE2::dof_type::nDofPerElement;
432 static const int nDofPerComponentTest = FE1::dof_type::fe_type::nLocalDof;
433 static const int nDofPerComponentTrial = FE2::dof_type::fe_type::nLocalDof;
434 static const int local_mat_traits = mpl::if_<mpl::equal_to<mpl::int_<nDofPerElementTrial>,mpl::int_<1> >,
435 mpl::int_<Eigen::ColMajor>,
436 mpl::int_<Eigen::RowMajor> >::type::value;
437 typedef Eigen::Matrix<value_type, nDofPerElementTest, nDofPerElementTrial,local_mat_traits> local_matrix_type;
438 typedef Eigen::Matrix<value_type, 2*nDofPerElementTest, 2*nDofPerElementTrial,Eigen::RowMajor> local2_matrix_type;
439 typedef Eigen::Matrix<value_type, nDofPerComponentTest, nDofPerComponentTrial,local_mat_traits> c_local_matrix_type;
440 typedef Eigen::Matrix<value_type, 2*nDofPerComponentTest, 2*nDofPerComponentTrial,Eigen::RowMajor> c_local2_matrix_type;
441 typedef Eigen::Matrix<int, nDofPerElementTest, 1> local_row_type;
442 typedef Eigen::Matrix<int, 2*nDofPerElementTest, 1> local2_row_type;
443 typedef Eigen::Matrix<int, nDofPerElementTrial, 1> local_col_type;
444 typedef Eigen::Matrix<int, 2*nDofPerElementTrial, 1> local2_col_type;
446 typedef Eigen::Matrix<int, nDofPerComponentTest, 1> c_local_row_type;
447 typedef Eigen::Matrix<int, 2*nDofPerComponentTest, 1> c_local2_row_type;
448 typedef Eigen::Matrix<int, nDofPerComponentTrial, 1> c_local_col_type;
449 typedef Eigen::Matrix<int, 2*nDofPerComponentTrial, 1> c_local2_col_type;
454 template<
typename Left,
typename Right>
455 map_trial_fecontext_type getMap( Left left, Right right )
457 return getMap( left, right, boost::is_same<Left, map_trial_fecontext_type>() );
459 template<
typename Left,
typename Right>
460 map_trial_fecontext_type getMap( Left left, Right , mpl::bool_<true> )
464 template<
typename Left,
typename Right>
465 map_trial_fecontext_type getMap( Left , Right right, mpl::bool_<false> )
467 return map_trial_fecontext_type( right );
470 template<
typename Left,
typename Right>
471 map_left_trial_fecontext_type getMapL( Left left, Right right )
473 return getMap( left, right, boost::is_same<Left, map_left_trial_fecontext_type>() );
475 template<
typename Left,
typename Right>
476 map_left_trial_fecontext_type getMapL( Left left, Right , mpl::bool_<true> )
480 template<
typename Left,
typename Right>
481 map_left_trial_fecontext_type getMapL( Left , Right right, mpl::bool_<false> )
486 Context( form_type& __form,
487 map_test_geometric_mapping_context_type
const& gmcTest,
488 map_trial_geometric_mapping_context_type
const& _gmcTrial,
489 map_geometric_mapping_expr_context_type
const & gmcExpr,
493 template<
typename IMTest,
typename IMTrial>
494 Context( form_type& __form,
495 map_test_geometric_mapping_context_type
const& gmcTest,
496 map_trial_geometric_mapping_context_type
const& _gmcTrial,
497 map_geometric_mapping_expr_context_type
const & gmcExpr,
500 IMTest
const& imTest, IMTrial
const& imTrial );
502 template<
typename IM2>
503 Context( form_type& __form,
504 map_test_geometric_mapping_context_type
const& gmcTest,
505 map_trial_geometric_mapping_context_type
const& _gmcTrial,
506 map_geometric_mapping_expr_context_type
const & gmcExpr,
511 template<
typename IM2>
512 Context( form_type& __form,
513 map_test_geometric_mapping_context_type
const& gmcTest,
514 map_trial_geometric_mapping_context_type
const& _gmcTrial,
515 map_geometric_mapping_expr_context_type
const & gmcExpr,
523 return trialElementId( trial_eid,mpl::int_<nDimDiffBetweenTestTrial>() );
529 const bool test_related_to_trial = M_form.testSpace()->mesh()->isSubMeshFrom( M_form.trialSpace()->mesh() );
530 const bool trial_related_to_test = M_form.trialSpace()->mesh()->isSubMeshFrom( M_form.testSpace()->mesh() );
531 if ( test_related_to_trial )
533 domain_eid = M_form.testSpace()->mesh()->subMeshToMesh( idElem );
534 DVLOG(2) <<
"[test_related_to_trial] test element id: " << idElem <<
" trial element id : " << domain_eid <<
"\n";
536 if( trial_related_to_test )
538 domain_eid = M_form.trialSpace()->mesh()->meshToSubMesh( idElem );
539 DVLOG(2) <<
"[trial_related_to_test] test element id: " << idElem <<
" trial element id : " << domain_eid <<
"\n";
547 const bool test_related_to_trial = M_form.testSpace()->mesh()->isSubMeshFrom( M_form.trialSpace()->mesh() );
548 const bool trial_related_to_test = M_form.trialSpace()->mesh()->isSubMeshFrom( M_form.testSpace()->mesh() );
549 if ( test_related_to_trial )
551 domain_eid = M_form.trialSpace()->mesh()->face( M_form.testSpace()->mesh()->subMeshToMesh( idElem )).element0().id();
552 DVLOG(2) <<
"[test_related_to_trial] test element id: " << idElem <<
" trial element id : " << domain_eid <<
"\n";
554 if( trial_related_to_test )
556 auto const& eltTest = M_form.testSpace()->mesh()->element(idElem);
557 std::set<size_type> idsFind;
558 for (uint16_type f=0;f< M_form.testSpace()->mesh()->numLocalFaces();++f)
560 const size_type idFind = M_form.trialSpace()->mesh()->meshToSubMesh( eltTest.face(f).id() );
563 if ( idsFind.size()>1 ) std::cout <<
" TODO trialElementId " << std::endl;
565 if ( idsFind.size()>0 )
566 domain_eid = *idsFind.begin();
570 DVLOG(2) <<
"[trial_related_to_test] test element id: " << idElem <<
" trial element id : " << domain_eid <<
"\n";
574 size_type trialElementId(
typename mesh_1_type::element_iterator it )
const
576 return trialElementId( it->id() );
578 size_type trialElementId(
typename mesh_1_type::element_type
const& e )
const
580 return trialElementId( e.id() );
584 size_type domain_eid = trialElementId( i );
589 bool isZero(
typename mesh_1_type::element_iterator it )
const
591 return this->isZero( it->id() );
593 bool isZero(
typename mesh_1_type::element_type
const& e )
const
595 return this->isZero( e.id() );
598 void update( map_test_geometric_mapping_context_type
const& gmcTest,
599 map_trial_geometric_mapping_context_type
const& _gmcTrial,
600 map_geometric_mapping_expr_context_type
const& gmcExpr );
602 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
603 map_trial_geometric_mapping_context_type
const& gmcTrial,
604 map_geometric_mapping_expr_context_type
const& gmcExpr,
605 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad );
607 void update( map_test_geometric_mapping_context_type
const& gmcTest,
608 map_trial_geometric_mapping_context_type
const& _gmcTrial,
609 map_geometric_mapping_expr_context_type
const& gmcExpr,
612 void update( map_test_geometric_mapping_context_type
const& gmcTest,
613 map_trial_geometric_mapping_context_type
const& gmcTrial,
614 map_geometric_mapping_expr_context_type
const& gmcExpr,
618 update( gmcTest, gmcTrial, gmcExpr );
621 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
622 map_trial_geometric_mapping_context_type
const& gmcTrial,
623 map_geometric_mapping_expr_context_type
const& gmcExpr,
625 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad )
628 updateInCaseOfInterpolate( gmcTest, gmcTrial, gmcExpr, indexLocalToQuad );
631 void update( map_test_geometric_mapping_context_type
const& gmcTest,
632 map_trial_geometric_mapping_context_type
const& gmcTrial,
633 map_geometric_mapping_expr_context_type
const& gmcExpr,
634 IM
const& im, mpl::int_<2> )
637 update( gmcTest, gmcTrial, gmcExpr, mpl::int_<2>() );
642 integrate( mpl::int_<fusion::result_of::size<GeomapTestContext>::type::value>() );
644 void integrateInCaseOfInterpolate( std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad,
645 bool isFirstExperience )
647 integrateInCaseOfInterpolate( mpl::int_<fusion::result_of::size<GeomapTestContext>::type::value>(),
655 assemble( fusion::at_key<gmc<0> >( M_test_gmc )->
id() );
659 void assemble( mpl::int_<2> )
661 assemble( fusion::at_key<gmc<0> >( M_test_gmc )->
id(),
662 fusion::at_key<test_gmc1>( M_test_gmc )->
id() );
666 void assembleInCaseOfInterpolate();
672 template<
typename Pts>
673 void precomputeBasisAtPoints( Pts
const& pts )
675 M_test_pc = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), pts ) );
676 M_trial_pc = trial_precompute_ptrtype(
new trial_precompute_type( M_form.trialSpace()->fe(), pts ) );
679 template<
typename PtsTest,
typename PtsTrial>
680 void precomputeBasisAtPoints( PtsTest
const& pts1,PtsTrial
const& pts2 )
682 M_test_pc = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), pts1 ) );
683 M_trial_pc = trial_precompute_ptrtype(
new trial_precompute_type( M_form.trialSpace()->fe(), pts2 ) );
692 template<
typename Pts>
693 void precomputeBasisAtPoints( uint16_type __f, permutation_1_type
const& __p, Pts
const& pts )
695 M_test_pc_face[__f][__p] = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), pts ) );
698 M_trial_pc_face[__f][__p] = trial_precompute_ptrtype(
new trial_precompute_type( M_form.trialSpace()->fe(), pts ) );
702 template<
typename PtsSet>
703 std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
704 precomputeTestBasisAtPoints( PtsSet
const& pts )
706 typedef typename boost::is_same< permutation_1_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
707 return precomputeTestBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
710 template<
typename PtsSet>
711 std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
712 precomputeTestBasisAtPoints( PtsSet
const& pts, mpl::bool_<false> )
714 std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> > testpc;
718 template<
typename PtsSet>
719 std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
720 precomputeTestBasisAtPoints( PtsSet
const& pts, mpl::bool_<true> )
722 QuadMapped<PtsSet> qm;
723 typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
724 typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
726 std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
728 for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
730 for ( permutation_type __p( permutation_type::IDENTITY );
731 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
733 testpc[__f][__p] = test_precompute_ptrtype(
new test_precompute_type( M_form.testSpace()->fe(), ppts[__f].find( __p )->second ) );
740 template<
typename PtsSet>
741 std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
742 precomputeTrialBasisAtPoints( PtsSet
const& pts )
744 typedef typename boost::is_same< permutation_2_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
745 return precomputeTrialBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
748 template<
typename PtsSet>
749 std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
750 precomputeTrialBasisAtPoints( PtsSet
const& pts, mpl::bool_<false> )
752 std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> > trialpc;
756 template<
typename PtsSet>
757 std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
758 precomputeTrialBasisAtPoints( PtsSet
const& pts, mpl::bool_<true> )
760 QuadMapped<PtsSet> qm;
761 typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
762 typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
764 std::map<uint16_type, std::map<permutation_type,trial_precompute_ptrtype> > trialpc;
766 for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
768 for ( permutation_type __p( permutation_type::IDENTITY );
769 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
771 trialpc[__f][__p] = trial_precompute_ptrtype(
new trial_precompute_type( M_form.trialSpace()->fe(), ppts[__f].find( __p )->second ) );
785 test_precompute_ptrtype
const& testPc( uint16_type __f,
786 permutation_1_type __p = permutation_1_type( permutation_1_type::NO_PERMUTATION ) )
const
792 return M_test_pc_face.find( __f )->second.find( __p )->second;
801 trial_precompute_ptrtype
const& trialPc( uint16_type __f,
802 permutation_2_type __p = permutation_2_type( permutation_2_type::NO_PERMUTATION ) )
const
808 return M_trial_pc_face.find( __f )->second.find( __p )->second;
814 void update( map_test_geometric_mapping_context_type
const& gmcTest,
815 map_trial_geometric_mapping_context_type
const& gmcTrial,
816 map_geometric_mapping_expr_context_type
const& gmcExpr,
819 void update( map_test_geometric_mapping_context_type
const& gmcTest,
820 map_trial_geometric_mapping_context_type
const& gmcTrial,
821 map_geometric_mapping_expr_context_type
const& gmcExpr,
824 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
825 map_trial_geometric_mapping_context_type
const& gmcTrial,
826 map_geometric_mapping_expr_context_type
const& gmcExpr,
829 void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type
const& gmcTest,
830 map_trial_geometric_mapping_context_type
const& gmcTrial,
831 map_geometric_mapping_expr_context_type
const& gmcExpr,
839 void integrateInCaseOfInterpolate( mpl::int_<1>,
840 std::vector<boost::tuple<size_type,size_type> >
const& indexLocalToQuad,
841 bool isFirstExperience );
845 const list_block_type& M_lb;
846 dof_1_type* M_test_dof;
847 dof_2_type* M_trial_dof;
849 test_precompute_ptrtype M_test_pc;
850 std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> > M_test_pc_face;
851 trial_precompute_ptrtype M_trial_pc;
852 std::map<uint16_type, std::map<permutation_2_type, trial_precompute_ptrtype> > M_trial_pc_face;
855 map_test_geometric_mapping_context_type M_test_gmc;
856 map_trial_geometric_mapping_context_type M_trial_gmc;
858 map_test_fecontext_type M_test_fec;
859 map_left_test_fecontext_type M_test_fec0;
860 map_right_test_fecontext_type M_test_fec1;
861 map_trial_fecontext_type M_trial_fec;
862 map_left_trial_fecontext_type M_trial_fec0;
863 map_right_trial_fecontext_type M_trial_fec1;
866 local_matrix_type M_rep;
867 local2_matrix_type M_rep_2;
868 local_row_type M_local_rows;
869 local2_row_type M_local_rows_2;
870 local_col_type M_local_cols;
871 local2_col_type M_local_cols_2;
872 local_row_type M_local_rowsigns;
873 local2_row_type M_local_rowsigns_2;
874 local_col_type M_local_colsigns;
875 local2_col_type M_local_colsigns_2;
877 c_local_matrix_type M_c_rep;
878 c_local2_matrix_type M_c_rep_2;
879 c_local_row_type M_c_local_rows;
880 c_local2_row_type M_c_local_rows_2;
881 c_local_col_type M_c_local_cols;
882 c_local2_col_type M_c_local_cols_2;
883 c_local_row_type M_c_local_rowsigns;
884 c_local2_row_type M_c_local_rowsigns_2;
885 c_local_col_type M_c_local_colsigns;
886 c_local2_col_type M_c_local_colsigns_2;
888 eval00_expr_ptrtype M_eval_expr00;
889 eval01_expr_ptrtype M_eval_expr01;
890 eval10_expr_ptrtype M_eval_expr10;
891 eval11_expr_ptrtype M_eval_expr11;
906 BilinearForm( space_1_ptrtype
const& __X1,
907 space_2_ptrtype
const& __X2,
912 bool do_threshold =
false,
913 value_type threshold = type_traits<value_type>::epsilon(),
914 size_type graph_hints = Pattern::COUPLED );
916 BilinearForm( space_1_ptrtype
const& __X1,
917 space_2_ptrtype
const& __X2,
919 list_block_type
const& __lb,
922 bool do_threshold =
false,
923 value_type threshold = type_traits<value_type>::epsilon(),
924 size_type graph_hints = Pattern::COUPLED );
926 BilinearForm( BilinearForm
const& __vf )
928 M_pattern( __vf.M_pattern ),
931 M_matrix( __vf.M_matrix ),
933 M_row_startInMatrix( __vf.M_row_startInMatrix ),
934 M_col_startInMatrix( __vf.M_col_startInMatrix ),
935 M_do_threshold( __vf.M_do_threshold ),
936 M_threshold( __vf.M_threshold )
958 M_pattern = form.M_pattern;
961 M_matrix = form.M_matrix;
962 M_row_startInMatrix = form.M_row_startInMatrix;
963 M_col_startInMatrix = form.M_col_startInMatrix;
973 template <
class ExprT>
974 BilinearForm&
operator=( Expr<ExprT>
const& expr );
979 template <
class ExprT>
980 BilinearForm& operator+=( Expr<ExprT>
const& expr );
982 BilinearForm& operator+=( BilinearForm& a )
987 M_matrix->addMatrix( 1.0, a.M_matrix );
999 value_type operator()( element_1_type
const& __v, element_2_type
const& __u )
const
1001 return M_matrix->energy( __v, __u );
1026 bool isPatternCoupled()
const
1029 return ctx.test( Pattern::COUPLED );
1035 bool isPatternDefault()
const
1038 return ctx.test( Pattern::DEFAULT );
1044 bool isPatternNeighbor()
const
1047 return ctx.test( Pattern::EXTENDED );
1049 bool isPatternExtended()
const
1052 return ctx.test( Pattern::EXTENDED );
1055 bool isPatternSymmetric()
const
1058 return ctx.test( Pattern::SYMMETRIC );
1064 space_2_ptrtype
const& trialSpace()
const
1072 space_1_ptrtype
const& testSpace()
const
1080 gm_1_ptrtype
const& gm()
const
1088 gm1_1_ptrtype
const& gm1()
const
1096 gm_2_ptrtype
const& gmTrial()
const
1104 gm1_2_ptrtype
const& gm1Trial()
const
1113 matrix_type
const& matrix()
const
1118 matrix_type& matrix()
1123 matrix_ptrtype
const& matrixPtr()
const
1128 matrix_ptrtype& matrixPtr()
1133 list_block_type
const& blockList()
const
1140 return M_row_startInMatrix;
1145 return M_col_startInMatrix;
1151 value_type threshold()
const
1159 bool doThreshold( value_type
const& v )
const
1161 return ( math::abs( v ) > M_threshold );
1167 bool doThreshold()
const
1169 return M_do_threshold;
1182 void setThreshold( value_type eps )
1191 void setDoThreshold(
bool do_threshold )
1193 M_do_threshold = do_threshold;
1204 void close() { M_matrix->close(); }
1216 void zeroRows( std::vector<int>
const& __dofs,
1217 Vector<value_type>
const& __values,
1218 Vector<value_type>& rhs,
1227 if ( M_do_threshold )
1229 if ( doThreshold( v ) )
1230 M_matrix->add( i+this->rowStartInMatrix(),
1231 j+this->colStartInMatrix(),
1236 M_matrix->add( i+this->rowStartInMatrix(),
1237 j+this->colStartInMatrix(),
1245 void addMatrix(
int* rows,
int nrows,
1246 int* cols,
int ncols,
1249 if ( this->rowStartInMatrix()!=0 )
1250 for (
int i=0; i<nrows; ++i )
1251 rows[i]+=this->rowStartInMatrix();
1253 if ( this->colStartInMatrix()!=0 )
1254 for (
int i=0; i<ncols; ++i )
1255 cols[i]+=this->colStartInMatrix();
1257 M_matrix->addMatrix( rows, nrows, cols, ncols, data );
1268 M_matrix->set( i, j, v );
1288 BOOST_PARAMETER_MEMBER_FUNCTION( (
typename Backend<value_type>::solve_return_type ),
1292 ( in_out( solution ),* )
1295 ( name, ( std::string ),
"" )
1297 ( rebuild, (
bool ),
false )
1300 return backend( _name=name, _kind=kind, _rebuild=rebuild )->solve( _matrix=this->matrixPtr(), _rhs=rhs.vectorPtr(),
1301 _solution=solution);
1304 BOOST_PARAMETER_MEMBER_FUNCTION( (
typename Backend<value_type>::solve_return_type ),
1308 ( in_out( solution ),* )
1312 ( prec, ( preconditioner_ptrtype ),
1313 preconditioner( _prefix=backend->prefix(),
1314 _matrix=this->matrixPtr(),
1315 _pc=backend->pcEnumType(),
1316 _pcfactormatsolverpackage=backend->matSolverPackageEnumType(),
1317 _backend=backend ) )
1320 return backend->solve( _matrix=this->matrixPtr(), _rhs=rhs.vectorPtr(),
1321 _solution=solution, _prec = prec );
1333 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr,
bool init, mpl::bool_<false> );
1334 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr,
bool init, mpl::bool_<true> );
1335 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr, mpl::bool_<true>, mpl::bool_<true> );
1336 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr, mpl::bool_<true>, mpl::bool_<false> );
1337 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr, mpl::bool_<true>, mpl::bool_<false>, mpl::bool_<true> );
1338 template <
class ExprT>
void assign( Expr<ExprT>
const& __expr, mpl::bool_<true>, mpl::bool_<false>, mpl::bool_<false> );
1344 space_1_ptrtype M_X1;
1345 space_2_ptrtype M_X2;
1347 matrix_ptrtype M_matrix;
1351 list_block_type M_lb;
1352 size_type M_row_startInMatrix,M_col_startInMatrix;
1354 bool M_do_threshold;
1355 value_type M_threshold;
1357 std::vector<size_type> M_n_nz;
1358 std::vector<size_type> M_n_oz;
1361 template<
typename FE1,
typename FE2,
typename ElemContType>
1362 BilinearForm<FE1, FE2, ElemContType>::BilinearForm( space_1_ptrtype
const& Xh,
1363 space_2_ptrtype
const& Yh,
1364 matrix_ptrtype& __M,
1369 value_type threshold,
1372 M_pattern( graph_hints ),
1376 M_do_build( build ),
1378 M_row_startInMatrix( rowstart ),
1379 M_col_startInMatrix( colstart ),
1380 M_do_threshold( do_threshold ),
1381 M_threshold( threshold )
1384 DVLOG(2) <<
"begin constructor with default listblock\n";
1386 if ( !M_matrix ) M_matrix = backend()->newMatrix( _test=M_X1, _trial=M_X2 );
1387 M_lb.push_back( Block ( 0, 0, 0, 0 ) );
1389 DVLOG(2) <<
" - form init in " << tim.elapsed() <<
"\n";
1390 DVLOG(2) <<
"begin constructor with default listblock done\n";
1393 template<
typename FE1,
typename FE2,
typename ElemContType>
1394 BilinearForm<FE1, FE2, ElemContType>::BilinearForm( space_1_ptrtype
const& Xh,
1395 space_2_ptrtype
const& Yh,
1396 matrix_ptrtype& __M,
1397 list_block_type
const& __lb,
1401 value_type threshold,
1404 M_pattern( graph_hints ),
1408 M_do_build( false ),
1410 M_row_startInMatrix( rowstart ),
1411 M_col_startInMatrix( colstart ),
1412 M_do_threshold( do_threshold ),
1413 M_threshold( threshold )
1415 if ( !M_matrix ) M_matrix = backend()->newMatrix( _test=M_X1, _trial=M_X2 );
1418 template<
typename FE1,
typename FE2,
typename ElemContType>
1419 template<
typename ExprT>
1421 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1427 DVLOG(2) <<
"[BilinearForm::assign<false>] start\n";
1428 typedef ublas::matrix_range<matrix_type> matrix_range_type;
1429 typename list_block_type::const_iterator __bit = M_lb.begin();
1430 typename list_block_type::const_iterator __ben = M_lb.end();
1432 for ( ; __bit != __ben; ++__bit )
1434 size_type g_ic_start = M_row_startInMatrix + __bit->globalRowStart();
1435 size_type g_jc_start = M_col_startInMatrix + __bit->globalColumnStart();
1436 M_matrix->zero( g_ic_start, g_ic_start + M_X1->nDof(),
1437 g_jc_start, g_jc_start + M_X2->nDof() );
1441 __expr.assemble( M_X1, M_X2, *
this );
1443 DVLOG(2) <<
"[BilinearForm::assign<false>] stop\n";
1446 template<
typename FE1,
typename FE2,
typename ElemContType>
1447 template<
typename ExprT>
1449 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1453 DVLOG(2) <<
"BilinearForm::assign() start loop on test spaces\n";
1455 if ( init ) M_matrix->zero();
1457 assign( __expr, mpl::bool_<true>(), mpl::bool_<( FE1::nSpaces > 1 && FE2::nSpaces > 1 )>() );
1458 DVLOG(2) <<
"BilinearForm::assign() stop loop on test spaces\n";
1461 template<
typename FE1,
typename FE2,
typename ElemContType>
1462 template<
typename ExprT>
1464 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1468 fusion::for_each( M_X1->functionSpaces(), make_bfassign2( *
this, __expr ) );
1470 template<
typename FE1,
typename FE2,
typename ElemContType>
1471 template<
typename ExprT>
1473 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1477 assign( __expr, mpl::bool_<true>(), mpl::bool_<false>(), mpl::bool_<( FE1::nSpaces > 1 )>() );
1480 template<
typename FE1,
typename FE2,
typename ElemContType>
1481 template<
typename ExprT>
1483 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1488 fusion::for_each( M_X1->functionSpaces(), make_bfassign3( *
this, __expr, M_X2, 0 ) );
1490 template<
typename FE1,
typename FE2,
typename ElemContType>
1491 template<
typename ExprT>
1493 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT>
const& __expr,
1498 fusion::for_each( M_X2->functionSpaces(), make_bfassign1( *
this, __expr, M_X1, 0 ) );
1501 template<
typename FE1,
typename FE2,
typename ElemContType>
1502 template<
typename ExprT>
1503 BilinearForm<FE1, FE2, ElemContType>&
1508 this->assign( __expr,
true, mpl::bool_<mpl::or_< mpl::bool_< ( FE1::nSpaces > 1 )>,
1509 mpl::bool_< ( FE2::nSpaces > 1 )> >::type::value >() );
1513 template<
typename FE1,
typename FE2,
typename ElemContType>
1514 template<
typename ExprT>
1515 BilinearForm<FE1, FE2, ElemContType>&
1516 BilinearForm<FE1, FE2, ElemContType>::operator+=( Expr<ExprT>
const& __expr )
1518 DVLOG(2) <<
"[BilinearForm::operator+=] start\n";
1519 this->assign( __expr,
false, mpl::bool_<mpl::or_< mpl::bool_< ( FE1::nSpaces > 1 )>,
1520 mpl::bool_< ( FE2::nSpaces > 1 )> >::type::value >() );
1521 DVLOG(2) <<
"[BilinearForm::operator+=] stop\n";
1526 template<
typename FE1,
typename FE2,
typename ElemContType>
1528 BilinearForm<FE1,FE2,ElemContType>::zeroRows( std::vector<int>
const& __dofs,
1529 Vector<value_type>
const&__values,
1530 Vector<value_type>& rhs,
1533 M_matrix->zeroRows( __dofs, __values, rhs, on_context );
1537 template<
typename BFType,
typename ExprType,
typename TestSpaceType>
1538 template<
typename SpaceType>
1539 void BFAssign1<BFType,ExprType,TestSpaceType>::operator()( boost::shared_ptr<SpaceType>
const& trial )
const
1541 if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
1543 DVLOG(2) <<
"[BFAssign1::operator()] expression has test functions index "
1544 << M_test_index <<
" : "
1545 << ExprType::template HasTestFunction<typename TestSpaceType::reference_element_type>::result <<
" (0:no, 1:yes)\n";
1546 DVLOG(2) <<
"[BFAssign1::operator()] expression has trial functions index "
1547 << M_trial_index <<
" :"
1548 << ExprType::template HasTrialFunction<typename SpaceType::reference_element_type>::result <<
" (0:no, 1:yes)\n";
1550 if ( !ExprType::template HasTestFunction<typename TestSpaceType::reference_element_type>::result ||
1551 !ExprType::template HasTrialFunction<typename SpaceType::reference_element_type>::result )
1557 DVLOG(2) <<
"[BFAssign1::operator()] terms found with "
1558 <<
"testindex: " << M_test_index <<
" trialindex: " << M_trial_index <<
"\n";
1559 typedef SpaceType trial_space_type;
1560 typedef TestSpaceType test_space_type;
1562 Feel::vf::list_block_type list_block;
1564 if ( M_bf.testSpace()->worldsComm()[M_test_index].globalSize()>1 )
1567 if (M_bf.testSpace()->hasEntriesForAllSpaces())
1568 list_block.push_back( Feel::vf::Block( 0, 0,
1569 M_bf.testSpace()->nLocalDofStart( M_test_index ),
1570 M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1572 list_block.push_back( Feel::vf::Block( 0, 0, 0, M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1579 list_block.push_back( Feel::vf::Block( 0,0,
1580 M_bf.testSpace()->nDofStart( M_test_index ),
1581 M_bf.trialSpace()->nDofStart( M_trial_index ) ) );
1584 typedef typename BFType::matrix_type matrix_type;
1585 typedef Feel::vf::detail::BilinearForm<test_space_type,
1587 ublas::vector_range<ublas::vector<double> > > bf_type;
1589 bf_type bf( M_test,trial, M_bf.matrixPtr(), list_block, M_bf.rowStartInMatrix(), M_bf.colStartInMatrix(), M_bf.doThreshold(), M_bf.threshold(), M_bf.pattern() );
1597 template<
typename BFType,
typename ExprType,
typename TrialSpaceType>
1598 template<
typename SpaceType>
1599 void BFAssign3<BFType,ExprType,TrialSpaceType>::operator()( boost::shared_ptr<SpaceType>
const& test )
const
1601 if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
1604 DVLOG(2) <<
"[BFAssign3::operator()] expression has trial functions index "
1605 << M_test_index <<
" : "
1606 << ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result <<
" (0:no, 1:yes)\n";
1607 DVLOG(2) <<
"[BFAssign3::operator()] expression has test functions index "
1608 << M_trial_index <<
" :"
1609 << ExprType::template HasTrialFunction<typename TrialSpaceType::reference_element_type>::result <<
" (0:no, 1:yes)\n";
1611 if ( !ExprType::template HasTrialFunction<typename TrialSpaceType::reference_element_type>::result ||
1612 !ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result )
1618 DVLOG(2) <<
"[BFAssign3::operator()] terms found with "
1619 <<
"testindex: " << M_test_index <<
" trialindex: " << M_trial_index <<
"\n";
1620 typedef SpaceType test_space_type;
1621 typedef TrialSpaceType trial_space_type;
1626 Feel::vf::list_block_type list_block;
1629 if ( M_bf.testSpace()->worldsComm()[M_test_index].globalSize()>1 )
1632 if (M_bf.testSpace()->hasEntriesForAllSpaces())
1633 list_block.push_back( Feel::vf::Block( 0, 0,
1634 M_bf.testSpace()->nLocalDofStart( M_test_index ),
1635 M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1637 list_block.push_back( Feel::vf::Block( 0, 0, 0, M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1644 list_block.push_back( Feel::vf::Block( 0,0,
1645 M_bf.testSpace()->nDofStart( M_test_index ),
1646 M_bf.trialSpace()->nDofStart( M_trial_index ) ) );
1649 typedef typename BFType::matrix_type matrix_type;
1650 typedef Feel::vf::detail::BilinearForm<test_space_type,
1652 ublas::vector_range<ublas::vector<double> > > bf_type;
1654 bf_type bf( test, M_trial, M_bf.matrixPtr(), list_block, M_bf.rowStartInMatrix(), M_bf.colStartInMatrix(), M_bf.doThreshold(), M_bf.threshold(), M_bf.pattern() );