30 #ifndef __Integrators_H
31 #define __Integrators_H 1
34 #include <boost/timer.hpp>
35 #include <feel/feelcore/feel.hpp>
36 #include <feel/feelcore/parameter.hpp>
40 #include <Eigen/Eigen>
53 #if defined( FEELPP_HAS_GOOGLE_PROFILER_H )
54 #include <google/profiler.h>
70 INT_ELEMENT_FACES = 5,
71 INT_INTERNAL_FACES = 6
82 template<
typename Elements,
typename Im,
typename Expr,
typename Im2=Im>
83 class Integrator:
public IntegratorBase
91 static const size_type context = Expr::context;
92 static const bool is_terminal =
false;
94 static const uint16_type imorder = 0;
95 static const bool imIsPoly =
true;
97 template<
typename Func>
98 struct HasTestFunction
100 static const bool result = Expr::template HasTestFunction<Func>::result;
103 template<
typename Func>
104 struct HasTrialFunction
106 static const bool result = Expr::template HasTrialFunction<Func>::result;
109 static const size_type iDim = boost::tuples::template element<0, Elements>::type::value;
116 typedef Integrator<Elements, Im, Expr, Im2> self_type;
118 typedef typename boost::tuples::template element<1, Elements>::type element_iterator;
120 typedef Expr expression_type;
121 typedef typename expression_type::value_type expression_value_type;
122 typedef ublas::vector<int> vector_dof_type;
129 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
130 typedef typename boost::remove_const<const_t>::type the_face_element_type;
131 typedef typename the_face_element_type::super2::template Element<the_face_element_type>::type the_element_type;
133 typedef typename mpl::if_<mpl::bool_<the_element_type::is_simplex>,
134 mpl::identity<typename Im::template apply<the_element_type::nDim, expression_value_type, Simplex>::type >,
135 mpl::identity<typename Im::template apply<the_element_type::nDim, expression_value_type, Hypercube>::type >
136 >::type::type im_type;
138 typedef typename mpl::if_<mpl::bool_<the_element_type::is_simplex>,
139 mpl::identity<typename Im2::template apply<the_element_type::nDim, expression_value_type, Simplex>::type >,
140 mpl::identity<typename Im2::template apply<the_element_type::nDim, expression_value_type, Hypercube>::type >
141 >::type::type im2_type;
143 typedef the_element_type element_type;
144 typedef typename the_element_type::gm_type gm_type;
145 typedef boost::shared_ptr<gm_type> gm_ptrtype;
146 typedef typename the_element_type::gm1_type gm1_type;
147 typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
149 typedef typename gm_type::template Context<expression_type::context|vm::JACOBIAN, the_element_type> gmc_type;
150 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
151 typedef typename gm1_type::template Context<expression_type::context|vm::JACOBIAN, the_element_type> gmc1_type;
152 typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
154 typedef typename gm_type::template precompute<im_type::numPoints>::type gmpc_type;
155 typedef typename gm_type::template precompute<im_type::numPoints>::ptrtype gmpc_ptrtype;
157 typedef typename gm_type::PreCompute gmpc_type;
158 typedef boost::shared_ptr<gmpc_type> gmpc_ptrtype;
159 typedef typename gm1_type::PreCompute gmpc1_type;
160 typedef boost::shared_ptr<gmpc1_type> gmpc1_ptrtype;
170 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
171 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
172 typedef typename eval_expr_type::shape shape;
174 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
175 typedef typename expression_type::template tensor<map_gmc1_type> eval_expr1_type;
176 typedef typename eval_expr1_type::shape shape1;
184 typedef typename mpl::if_<mpl::and_<mpl::equal_to<mpl::int_<shape::M>, mpl::int_<1> >,
185 mpl::equal_to<mpl::int_<shape::N>, mpl::int_<1> > >,
186 mpl::identity<expression_value_type>,
187 mpl::identity<Eigen::Matrix<expression_value_type, shape::M, shape::N> > >::type::type value_type;
189 typedef Eigen::Matrix<expression_value_type, shape::M, shape::N> value_type;
191 typedef Eigen::Matrix<expression_value_type, shape::M, shape::N> matrix_type;
192 static value_type zero( mpl::bool_<false> )
194 return value_type::Zero();
196 static value_type zero( mpl::bool_<true> )
200 static value_type zero()
202 return zero( boost::is_scalar<value_type>() );
207 typedef typename eval::im_type im_type;
208 typedef typename eval::im2_type im2_type;
209 typedef typename im_type::face_quadrature_type im_face_type;
210 typedef typename im2_type::face_quadrature_type im2_face_type;
212 typedef typename eval::matrix_type matrix_type;
213 typedef typename eval::matrix_type value_type;
220 Integrator( Elements
const& elts, Im
const& , expression_type
const& __expr, GeomapStrategyType gt, Im2
const& ,
bool use_tbb,
int grainsize, std::string
const& partitioner,
221 boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > qpl )
224 M_eltbegin( elts.template get<1>() ),
225 M_eltend( elts.template get<2>() ),
230 M_use_tbb( use_tbb ),
231 M_grainsize( grainsize ),
232 M_partitioner( partitioner ),
235 M_elts.push_back( elts );
236 DLOG(INFO) <<
"Integrator constructor from expression\n";
239 Integrator( std::list<Elements>
const& elts, Im
const& , expression_type
const& __expr,
240 GeomapStrategyType gt, Im2
const& ,
bool use_tbb,
int grainsize, std::string
const& partitioner,
241 boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > qpl )
248 M_use_tbb( use_tbb ),
249 M_grainsize( grainsize ),
250 M_partitioner( partitioner ),
253 DLOG(INFO) <<
"Integrator constructor from expression\n";
256 M_eltbegin = elts.begin()->template get<1>();
257 M_eltend = elts.begin()->template get<2>();
261 Integrator( Integrator
const& __vfi )
263 M_elts( __vfi.M_elts) ,
264 M_eltbegin( __vfi.M_eltbegin ),
265 M_eltend( __vfi.M_eltend ),
267 M_im2( __vfi.M_im2 ),
268 M_expr( __vfi.M_expr ),
270 M_use_tbb( __vfi.M_use_tbb ),
271 M_grainsize( __vfi.M_grainsize ),
272 M_partitioner( __vfi.M_partitioner ),
275 DLOG(INFO) <<
"Integrator copy constructor\n";
278 virtual ~Integrator() {}
286 template<
typename TheExpr>
289 typedef typename expression_type::template Lambda<TheExpr>::type expr_type;
290 typedef _Q< ExpressionOrder<Elements,expr_type>::value > quad_type;
291 typedef _Q< ExpressionOrder<Elements,expr_type>::value_1 > quad1_type;
292 typedef Integrator<Elements, quad_type, expr_type, quad1_type> type;
295 template<
typename ExprT>
296 typename Lambda<ExprT>::type
297 operator()( ExprT
const& e )
300 typedef decltype(expr(M_expr(e))) t1;
301 typedef typename Lambda<ExprT>::expr_type t2;
302 BOOST_MPL_ASSERT_MSG( (boost::is_same<t1,t2>), INVALID_TYPE_IN_META_EXPRESSION,
303 (decltype(expr(M_expr(e))), typename Lambda<ExprT>::expr_type ) );
305 auto new_expr = M_expr(e);
306 typedef decltype(new_expr) expr_type;
307 typedef typename Lambda<ExprT>::expr_type e_type;
308 typedef _Q< ExpressionOrder<Elements,e_type>::value > quad_type;
309 typedef _Q< ExpressionOrder<Elements,e_type>::value_1 > quad1_type;
310 typedef boost::shared_ptr<QuadPtLocalization<Elements,quad_type,expr_type > > quadptloc_ptrtype;
314 auto i = Integrator<Elements, quad_type, expr_type, quad1_type>( M_elts, quad, new_expr, M_gt, quad1, M_use_tbb, M_grainsize, M_partitioner, quadptloc_ptrtype() );
315 DLOG(INFO) << " -- M_elts size=" << M_elts.size() << "\n";
316 DLOG(INFO) << " -- nelts=" << std::distance( M_eltbegin, M_eltend ) << "\n";
317 DLOG(INFO) << " --
integrate: quad = " << i.im().nPoints() << "\n";
318 DLOG(INFO) << " --
integrate: quad1 = " << i.im2().nPoints() << "\n";
319 i.setBeginElement( M_eltbegin );
320 i.setEndElement( M_eltend );
336 im_type const& im()
const
347 im2_type
const& im2()
const
358 im_face_type im( uint16_type f )
const
360 return M_im.face( f );
369 im2_face_type im2( uint16_type f )
const
371 return M_im2.face( f );
380 expression_type
const& expression()
const
391 GeomapStrategyType geomapIntegratorType()
const
400 element_iterator beginElement()
const
409 element_iterator endElement()
const
420 bool isSymmetric()
const
422 return M_expr.isSymmetric();
430 void setBeginElement( element_iterator it ) { M_eltbegin = it; }
431 void setEndElement( element_iterator en ) { M_eltend = en; }
440 template<
typename Elem1,
typename Elem2,
typename FormType>
441 void assemble( boost::shared_ptr<Elem1>
const& __u,
442 boost::shared_ptr<Elem2>
const& __v,
443 FormType& __form )
const;
445 template<
typename Elem1,
typename FormType>
446 void assemble( boost::shared_ptr<Elem1>
const& __v,
447 FormType& __form )
const;
450 template<
typename P0hType>
451 typename P0hType::element_type broken( boost::shared_ptr<P0hType>& P0h )
const
453 auto p0 = broken( P0h, mpl::int_<iDim>() );
458 evaluate(
bool parallel=
true,
459 WorldComm
const& worldcomm = Environment::worldComm() ) const
462 BOOST_PARAMETER_MEMBER_FUNCTION( ( matrix_type ),
468 ( parallel,*(
bool ),
true ) ) )
471 typename eval::matrix_type loc = evaluate( mpl::int_<iDim>() );
478 typename eval::matrix_type glo( loc );
484 if ( worldcomm.localSize() > 1 )
486 mpi::all_reduce( worldcomm.localComm(),
489 [] ( matrix_type
const& x, matrix_type
const& y )
499 #if defined( FEELPP_HAS_TBB )
500 template<
typename FormType,
typename ExprType,
typename IMType,
typename EltType>
507 typedef typename eval::gm_type gm_type;
508 typedef typename eval::gmc_type gmc_type;
509 typedef typename eval::gmpc_type gmpc_type;
510 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
511 typedef boost::shared_ptr<gmpc_type> gmpc_ptrtype;
512 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
513 typedef typename FormType::template Context<map_gmc_type, expression_type, im_type> form_context_type;
515 typedef form_context_type fcb_type;
516 typedef boost::shared_ptr<fcb_type> fcb_ptrtype;
519 Context( FormType& _form,
520 ExprType
const& _expr,
522 EltType
const& _elt )
524 M_geopc( new gmpc_type( _form.gm(), _im.
points() ) ),
525 M_c( new gmc_type( _form.gm(), _elt, M_geopc ) ),
526 M_formc( new form_context_type( _form,
527 fusion::make_pair<vf::detail::gmc<0> >( M_c ),
528 fusion::make_pair<vf::detail::gmc<0> >( M_c ),
534 Context( Context
const& c )
536 M_geopc( new gmpc_type( *c.M_geopc ) ),
537 M_c( new gmc_type( *c.M_c ) ),
538 M_formc( new form_context_type( *c.M_formc ) )
542 typedef typename std::vector<boost::reference_wrapper<const typename eval::element_type> >::iterator elt_iterator;
543 void operator() (
const tbb::blocked_range<elt_iterator>& r )
const
547 tbb::mutex::scoped_lock lock( m );
549 auto mapgmc = fusion::make_pair<vf::detail::gmc<0> >( M_c );
551 for (
auto _elt = r.begin(); _elt != r.end(); ++_elt )
553 M_c->update( *_elt );
554 M_formc->update( mapgmc, mapgmc );
555 M_formc->integrate();
563 for (
auto _elt = r.begin(); _elt != r.end(); ++_elt )
565 M_c->update( *_elt );
566 M_formc->update( fusion::make_pair<vf::detail::gmc<0> >( M_c ),
567 fusion::make_pair<vf::detail::gmc<0> >( M_c ) );
568 M_formc->integrate();
573 mutable gmpc_ptrtype M_geopc;
574 mutable gmc_ptrtype M_c;
575 mutable fcb_ptrtype M_formc;
579 template<
typename ExprType,
typename IMType,
typename EltType>
580 class ContextEvaluate
583 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
587 typedef ExprType expression_type;
588 typedef typename eval::gm_type gm_type;
589 typedef boost::shared_ptr<gm_type> gm_ptrtype;
590 typedef typename eval::gmc_type gmc_type;
591 typedef typename eval::gmpc_type gmpc_type;
592 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
593 typedef boost::shared_ptr<gmpc_type> gmpc_ptrtype;
594 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
595 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
596 typedef typename eval_expr_type::shape shape;
597 typedef IMType im_type;
598 typedef typename eval::matrix_type value_type;
599 ContextEvaluate( ExprType
const& _expr,
601 EltType
const& _elt )
603 M_gm( new gm_type( *_elt.gm() ) ),
604 M_geopc( new gmpc_type( M_gm, _im.
points() ) ),
605 M_c( new gmc_type( M_gm, _elt, M_geopc ) ),
606 M_expr( _expr, map_gmc_type( fusion::make_pair<vf::detail::gmc<0> >( M_c ) ) ),
608 M_ret( eval::matrix_type::Zero() )
611 ContextEvaluate( ContextEvaluate& c, tbb::split )
613 M_gm( new gm_type( *c.M_gm ) ),
615 M_geopc( c.M_geopc ),
616 M_c( new gmc_type( M_gm, c.M_c->element(), M_geopc ) ),
619 M_ret( eval::matrix_type::Zero() )
622 ContextEvaluate( ContextEvaluate
const& c )
624 M_gm( new gm_type( *c.M_gm ) ),
626 M_geopc( c.M_geopc ),
627 M_c( new gmc_type( M_gm, c.M_c->element(), M_geopc ) ),
634 typedef typename std::vector<boost::reference_wrapper<const typename eval::element_type> >::iterator elt_iterator;
635 void operator() (
const tbb::blocked_range<elt_iterator>& r )
639 for (
auto _elt = r.begin(); _elt != r.end(); ++_elt )
641 M_c->update( *_elt );
642 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( M_c ) );
644 M_expr.update( mapgmc );
649 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
650 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
652 M_ret( c1,c2 ) += M_im( M_expr, c1, c2 );
661 for (
int i = 0; i < 10000; ++i )
662 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
663 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
665 M_ret( c1,c2 ) += i*i;
671 void join( ContextEvaluate
const& other )
673 M_ret += other.M_ret;
676 value_type result()
const
682 gmpc_ptrtype M_geopc;
684 eval_expr_type M_expr;
688 #endif // FEELPP_HAS_TBB
692 template<
typename FormType>
693 void assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<true> ,
bool hasRelation )
const;
694 template<
typename FormType>
695 void assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<false> ,
bool hasRelation )
const;
697 template<
typename FormType>
698 void assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<true> ,
bool hasRelation )
const;
699 template<
typename FormType>
700 void assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<false> ,
bool hasRelation )
const;
702 template<
typename FE1,
typename FE2,
typename ElemContType>
703 void assembleWithRelationDifferentMeshType( vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const;
704 template<
typename FE,
typename VectorType,
typename ElemContType>
705 void assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const;
706 template<
typename FE1,
typename FE2,
typename ElemContType>
707 void assembleWithRelationDifferentMeshType( vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_FACES> )
const;
708 template<
typename FE,
typename VectorType,
typename ElemContType>
709 void assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> )
const;
711 template<
typename FE1,
typename FE2,
typename ElemContType>
712 void assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const;
713 template<
typename FE1,
typename FE2,
typename ElemContType>
714 void assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_FACES> )
const;
715 template<
typename FE,
typename VectorType,
typename ElemContType>
716 void assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const;
717 template<
typename FE,
typename VectorType,
typename ElemContType>
718 void assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> )
const;
721 template<
typename P0hType>
722 typename P0hType::element_type broken( boost::shared_ptr<P0hType>& P0h, mpl::int_<MESH_ELEMENTS> )
const;
723 template<
typename P0hType>
724 typename P0hType::element_type broken( boost::shared_ptr<P0hType>& P0h, mpl::int_<MESH_FACES> )
const;
726 typename eval::matrix_type evaluate( mpl::int_<MESH_ELEMENTS> )
const;
727 typename eval::matrix_type evaluate( mpl::int_<MESH_FACES> )
const;
728 typename eval::matrix_type evaluate( mpl::int_<MESH_POINTS> )
const;
732 std::list<Elements> M_elts;
733 element_iterator M_eltbegin;
734 element_iterator M_eltend;
735 mutable im_type M_im;
736 mutable im2_type M_im2;
737 expression_type M_expr;
738 GeomapStrategyType M_gt;
741 std::string M_partitioner;
742 mutable boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > M_QPL;
748 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
749 template<
typename Elem1,
typename Elem2,
typename FormType>
751 Integrator<Elements, Im, Expr, Im2>::assemble( boost::shared_ptr<Elem1>
const& __u,
752 boost::shared_ptr<Elem2>
const& __v,
753 FormType& __form )
const
756 details::GlobalMatrixAssembler<iDim, self_type, Elem1, Elem2, FormType>( *
this,
762 typedef typename boost::is_same<typename eval::gmc_type::element_type,typename Elem1::mesh_type::element_type>::type same1_mesh_type;
763 typedef typename boost::is_same<typename eval::gmc_type::element_type,typename Elem2::mesh_type::element_type>::type same2_mesh_type;
764 typedef typename boost::mpl::and_< same1_mesh_type,same2_mesh_type>::type same_mesh_type;
766 element_iterator it, en;
767 bool findEltForInit =
false;
768 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
770 it = lit->template get<1>();
771 en = lit->template get<2>();
777 findEltForInit =
true;
779 if (!findEltForInit)
return;
781 bool same_mesh = (
dynamic_cast<void*
>(
const_cast<MeshBase*
>( it->mesh() ) ) ==
dynamic_cast<void*
>( __u->mesh().get() ) &&
782 dynamic_cast<void*>( const_cast<MeshBase*>( it->mesh() ) ) ==
dynamic_cast<void*
>( __v->mesh().get() ) );
783 const bool test_related_to_trial = __v->mesh()->isRelatedTo( __u->mesh() ) && ( __u->mesh()->isSameMesh( it->mesh() ) || __v->mesh()->isSameMesh( it->mesh() ) );
784 if ( same_mesh || test_related_to_trial )
786 DLOG(INFO) <<
"[integrator::assemble bilinear form] with_relation_mesh (same_mesh: " << same_mesh_type::value <<
" test_related_to_trial: " << test_related_to_trial <<
")\n";
787 assemble( __form, mpl::int_<iDim>(), mpl::bool_<same_mesh_type::value>(), test_related_to_trial );
791 DLOG(INFO) <<
"[integrator::assemble bilinear form] no_relation_mesh\n";
792 assemble( __form, mpl::int_<iDim>(), mpl::bool_<false>(), test_related_to_trial );
798 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
799 template<
typename Elem1,
typename FormType>
801 Integrator<Elements, Im, Expr, Im2>::assemble( boost::shared_ptr<Elem1>
const& __v,
802 FormType& __form )
const
805 details::GlobalVectorAssembler<iDim, self_type, Elem1, FormType>( *
this,
810 typedef typename boost::is_same<typename eval::gmc_type::element_type,typename Elem1::mesh_type::element_type>::type same_mesh_type;
812 element_iterator it, en;
813 bool findEltForInit =
false;
814 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
816 it = lit->template get<1>();
817 en = lit->template get<2>();
823 findEltForInit =
true;
825 if (!findEltForInit)
return;
828 if ( dynamic_cast<void*>( const_cast<MeshBase*>( it->mesh() ) ) ==
dynamic_cast<void*
>( __v->mesh().get() ) )
829 assemble( __form, mpl::int_<iDim>(), mpl::bool_<same_mesh_type::value>(),
true );
832 assemble( __form, mpl::int_<iDim>(), mpl::bool_<false>(),
false );
836 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
837 template<
typename FormType>
839 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<true> ,
bool )
const
841 boost::timer __timer;
842 LOG(INFO) <<
"[integrator::assemble FormType& __form, mpl::int_<MESH_ELEMENTS> /**/, mpl::bool_<true>\n";
844 #if defined(FEELPP_HAS_TBB)
855 typedef typename eval::gm_type gm_type;
856 typedef typename eval::gm1_type gm1_type;
857 typedef typename eval::gmc_type gmc_type;
858 typedef typename eval::gmc1_type gmc1_type;
859 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
860 typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
862 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
863 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
864 typedef typename FormType::template Context<map_gmc_type, expression_type, im_type> form_context_type;
865 typedef typename FormType::template Context<map_gmc1_type, expression_type, im2_type> form1_context_type;
867 typedef form_context_type fcb_type;
868 typedef form1_context_type fcb1_type;
869 typedef fcb_type* focb_ptrtype;
870 typedef fcb1_type* focb1_ptrtype;
876 typename eval::gmpc_ptrtype __geopc(
new typename eval::gmpc_type( __form.gm(), this->im().points() ) );
877 typename eval::gmpc1_ptrtype __geopc1(
new typename eval::gmpc1_type( __form.gm1(), this->im2().points() ) );
880 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
882 element_iterator it = lit->template get<1>();
883 element_iterator en = lit->template get<2>();
884 DLOG(INFO) <<
"integrating over " << std::distance( it, en ) <<
" elements\n";
891 if ( it->mesh()->isSubMeshFrom( __form.testSpace()->mesh() ) )
892 idEltTestInit = it->mesh()->subMeshToMesh( idEltTestInit );
893 else if ( __form.testSpace()->mesh()->isSubMeshFrom( it->mesh() ) )
894 idEltTestInit = __form.testSpace()->mesh()->meshToSubMesh( idEltTestInit );
896 auto const& eltTestInit = __form.testSpace()->mesh()->element( idEltTestInit );
899 gmc_ptrtype __c(
new gmc_type( __form.gm(), eltTestInit, __geopc ) );
900 gmc1_ptrtype __c1(
new gmc1_type( __form.gm1(), eltTestInit, __geopc1 ) );
903 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
904 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
906 focb_ptrtype formc(
new form_context_type( __form,
912 focb1_ptrtype formc1(
new form1_context_type( __form,
920 boost::timer ti0,ti1, ti2, ti3;
930 for ( ; it != en; ++it )
933 if ( it->mesh()->isSubMeshFrom( __form.testSpace()->mesh() ) )
934 idElt = it->mesh()->subMeshToMesh( idElt );
935 else if ( __form.testSpace()->mesh()->isSubMeshFrom( it->mesh() ) )
936 idElt = __form.testSpace()->mesh()->meshToSubMesh( idElt );
938 auto const& eltTest = __form.testSpace()->mesh()->element( idElt );
940 if ( formc->isZero( idElt ) )
946 case GeomapStrategyType::GEOMAP_HO:
949 __c->update( eltTest );
952 std::cout <<
"Element: " << it->id() <<
"\n"
953 <<
" o - points : " << it->G() <<
"\n"
954 <<
" o - quadrature :\n"
955 <<
" ref : " << this->im().points() <<
"\n"
956 <<
" real : " << __c->xReal() <<
"\n";
959 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
960 formc->update( mapgmc,mapgmc,mapgmc );
976 case GeomapStrategyType::GEOMAP_O1:
979 __c1->update( eltTest );
982 DLOG(INFO) <<
"Element: " << it->id() <<
"\n"
983 <<
" o - points : " << it->G() <<
"\n"
984 <<
" o - quadrature :\n"
985 <<
" ref : " << this->im().points() <<
"\n"
986 <<
" real : " << __c->xReal() <<
"\n";
989 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
990 formc1->update( mapgmc1,mapgmc1,mapgmc1 );
1006 case GeomapStrategyType::GEOMAP_OPT:
1008 if ( it->isOnBoundary() )
1011 __c->update( eltTest );
1014 DLOG(INFO) <<
"Element: " << it->id() <<
"\n"
1015 <<
" o - points : " << it->G() <<
"\n"
1016 <<
" o - quadrature :\n"
1017 <<
" ref : " << this->im().points() <<
"\n"
1018 <<
" real : " << __c->xReal() <<
"\n";
1021 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
1022 formc->update( mapgmc,mapgmc,mapgmc );
1041 __c1->update( eltTest );
1044 DLOG(INFO) <<
"Element: " << it->id() <<
"\n"
1045 <<
" o - points : " << it->G() <<
"\n"
1046 <<
" o - quadrature :\n"
1047 <<
" ref : " << this->im().points() <<
"\n"
1048 <<
" real : " << __c->xReal() <<
"\n";
1051 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
1052 formc1->update( mapgmc1,mapgmc1,mapgmc1 );
1058 formc1->integrate();
1075 DLOG(INFO) <<
"[elements] Overall geometric mapping update time : " << ( t0+t1+t2 ) <<
" per element:" << ( t0+t1+t2 )/std::distance( this->beginElement(), this->endElement() ) <<
"\n";
1076 DLOG(INFO) <<
"[elements] Overall geometric mapping update time : " << t0 <<
"\n";
1077 DLOG(INFO) <<
"[elements] Overall form update time : " << t1 <<
"\n";
1078 DLOG(INFO) <<
"[elements] Overall local assembly time : " << t2 <<
"\n";
1079 DLOG(INFO) <<
"[elements] Overall global assembly time : " << t3 <<
"\n";
1081 DLOG(INFO) <<
"integrating over elements done in " << __timer.elapsed() <<
"s\n";
1086 #if defined( FEELPP_HAS_TBB )
1090 element_iterator it = this->beginElement();
1091 element_iterator en = this->endElement();
1096 std::vector<boost::reference_wrapper<const typename eval::element_type> > _v;
1098 for (
auto _it = it; _it != en; ++_it )
1099 _v.push_back( boost::cref( *_it ) );
1102 tbb::blocked_range<decltype( _v.begin() )> r( _v.begin(), _v.end(), std::distance( it, en ) );
1103 Context<FormType,expression_type, im_type, typename eval::the_element_type> thecontext ( __form,
1107 tbb::parallel_for( r, thecontext );
1110 #endif // FEELPP_HAS_TBB
1113 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1114 template<
typename FormType>
1116 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<false> ,
bool hasRelation )
const
1120 assembleWithRelationDifferentMeshType( __form,mpl::int_<MESH_ELEMENTS>() );
1124 assembleInCaseOfInterpolate( __form,mpl::int_<MESH_ELEMENTS>() );
1130 template<
typename SpaceType,
typename ImType,
typename GmcType,
typename GmcExprType>
1131 boost::shared_ptr<GmcType>
1132 buildGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType>
const& ,
typename GmcType::gm_ptrtype ,
1133 size_type ,boost::shared_ptr<GmcExprType> gmcExpr,mpl::int_<0> )
1138 template<
typename SpaceType,
typename ImType,
typename GmcType,
typename GmcExprType>
1139 boost::shared_ptr<GmcType>
1140 buildGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType>
const& space,
typename GmcType::gm_ptrtype gm,
1141 size_type idElt,boost::shared_ptr<GmcExprType> ,mpl::int_<1> )
1143 typedef typename SpaceType::gm_type::precompute_type geopc_type;
1144 typedef typename SpaceType::gm_type::precompute_ptrtype geopc_ptrtype;
1145 typedef typename QuadMapped<ImType>::permutation_type permutation_type;
1146 typedef typename QuadMapped<ImType>::permutation_points_type permutation_points_type;
1147 typedef GmcType gmc_type;
1148 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1151 QuadMapped<ImType> qm;
1152 permutation_points_type ppts( qm( im ) );
1153 std::vector<std::map<permutation_type, geopc_ptrtype> > __geopc( im.nFaces() );
1155 for ( uint16_type __f = 0; __f < im.nFaces(); ++__f )
1157 for ( permutation_type __p( permutation_type::IDENTITY );
1158 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1160 __geopc[__f][__p] = geopc_ptrtype(
new geopc_type( gm, ppts[__f].find( __p )->second ) );
1164 auto const& eltInit = space->mesh()->face( idElt );
1165 gmc_ptrtype gmc(
new gmc_type( gm, eltInit.element0(), __geopc, eltInit.pos_first() ) );
1169 template<
typename SpaceType,
typename ImType,
typename GmcType,
typename GmcExprType>
1171 updateGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType>
const& ,
1172 boost::shared_ptr<GmcType> , boost::shared_ptr<GmcExprType> ,
1177 template<
typename SpaceType,
typename ImType,
typename GmcType,
typename GmcExprType>
1179 updateGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType>
const& space,
1180 boost::shared_ptr<GmcType> gmc, boost::shared_ptr<GmcExprType> gmcExpr,
1183 typedef typename QuadMapped<ImType>::permutation_type permutation_type;
1185 auto const& theface = space->mesh()->face( idElt );
1186 bool findPermutation=
false;
1187 for ( permutation_type __p( permutation_type::IDENTITY );
1188 __p < permutation_type( permutation_type::N_PERMUTATIONS ) && !findPermutation; ++__p )
1190 gmc->update( theface.element0(), theface.pos_first(), __p );
1193 for ( uint16_type i=0;i<gmc->nPoints() && check;++i )
1194 for (uint16_type d=0;d<GmcType::NDim;++d)
1195 check = check && ( std::abs(gmc->xReal(i)[d]-gmcExpr->xReal(i)[d])<1e-8 );
1197 if (check) findPermutation=
true;
1199 CHECK(findPermutation) <<
"the permutation of quad point is not find\n";
1204 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1205 template<
typename FE1,
typename FE2,
typename ElemContType>
1207 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1208 mpl::int_<MESH_ELEMENTS> )
const
1210 LOG(INFO) <<
"[integrator::assembleWithRelationDifferentMeshType] vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS>\n";
1213 typedef typename eval::gm_type gm_expr_type;
1214 typedef typename gm_expr_type::template Context<expression_type::context|vm::POINT|vm::JACOBIAN, typename eval::element_type> gmc_expr_type;
1215 typedef boost::shared_ptr<gmc_expr_type> gmc_expr_ptrtype;
1216 typedef typename gm_expr_type::precompute_type pc_expr_type;
1217 typedef typename gm_expr_type::precompute_ptrtype pc_expr_ptrtype;
1218 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_expr_ptrtype> > map_gmc_expr_type;
1221 typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
1223 typedef typename FormType::gm_1_type gm_formTest_type;
1224 typedef typename FormType::gm1_1_type gm1_formTest_type;
1225 typedef typename FormType::mesh_element_1_type geoelement_formTest_type;
1226 typedef typename gm_formTest_type::template Context<expression_type::context|vm::POINT,geoelement_formTest_type> gmc_formTest_type;
1227 typedef typename gm1_formTest_type::template Context<expression_type::context|vm::POINT,geoelement_formTest_type> gmc1_formTest_type;
1229 typedef boost::shared_ptr<gmc_formTest_type> gmc_formTest_ptrtype;
1230 typedef typename gm_formTest_type::precompute_type pc_formTest_type;
1231 typedef typename gm_formTest_type::precompute_ptrtype pc_formTest_ptrtype;
1232 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTest_ptrtype> > map_gmc_formTest_type;
1234 typedef typename FormType::gm_2_type gm_formTrial_type;
1235 typedef typename FormType::gm1_2_type gm1_formTrial_type;
1236 typedef typename FormType::mesh_element_2_type geoelement_formTrial_type;
1237 typedef typename gm_formTrial_type::template Context<expression_type::context|vm::POINT,geoelement_formTrial_type> gmc_formTrial_type;
1238 typedef typename gm1_formTrial_type::template Context<expression_type::context|vm::POINT,geoelement_formTrial_type> gmc1_formTrial_type;
1239 typedef boost::shared_ptr<gmc_formTrial_type> gmc_formTrial_ptrtype;
1240 typedef typename gm_formTrial_type::precompute_type pc_formTrial_type;
1241 typedef typename gm_formTrial_type::precompute_ptrtype pc_formTrial_ptrtype;
1242 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTrial_ptrtype> > map_gmc_formTrial_type;
1245 typedef typename FormType::template Context<map_gmc_formTest_type, expression_type, im_type,map_gmc_expr_type, map_gmc_formTrial_type> form_context_type;
1246 typedef form_context_type fcb_type;
1247 typedef fcb_type* focb_ptrtype;
1251 static const uint16_type nDimTest = gm_formTest_type::nDim;
1252 static const uint16_type nDimTrial = gm_formTrial_type::nDim;
1253 static const uint16_type nDimRange = gm_expr_type::nDim;
1254 static const uint16_type gmTestRangeRelation = ( nDimTest > nDimRange )? nDimTest-nDimRange : nDimRange-nDimTest;
1255 static const uint16_type gmTrialRangeRelation = ( nDimTrial > nDimRange )? nDimTrial-nDimRange : nDimRange-nDimTrial;
1258 typedef typename mpl::if_<mpl::bool_<eval::the_element_type::is_simplex>,
1259 mpl::identity<typename Im::template apply<geoelement_formTest_type::nDim, expression_value_type, Simplex>::type >,
1260 mpl::identity<typename Im::template apply<geoelement_formTest_type::nDim, expression_value_type, Hypercube>::type >
1261 >::type::type im_formtest_type;
1263 typedef typename mpl::if_<mpl::bool_<eval::the_element_type::is_simplex>,
1264 mpl::identity<typename Im2::template apply<geoelement_formTest_type::nDim, expression_value_type, Simplex>::type >,
1265 mpl::identity<typename Im2::template apply<geoelement_formTest_type::nDim, expression_value_type, Hypercube>::type >
1266 >::type::type im1_formtest_type;
1268 typedef typename mpl::if_<mpl::bool_<eval::the_element_type::is_simplex>,
1269 mpl::identity<typename Im::template apply<geoelement_formTrial_type::nDim, expression_value_type, Simplex>::type >,
1270 mpl::identity<typename Im::template apply<geoelement_formTrial_type::nDim, expression_value_type, Hypercube>::type >
1271 >::type::type im_formtrial_type;
1273 typedef typename mpl::if_<mpl::bool_<eval::the_element_type::is_simplex>,
1274 mpl::identity<typename Im2::template apply<geoelement_formTrial_type::nDim, expression_value_type, Simplex>::type >,
1275 mpl::identity<typename Im2::template apply<geoelement_formTrial_type::nDim, expression_value_type, Hypercube>::type >
1276 >::type::type im1_formtrial_type;
1281 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
1283 auto elt_it = lit->template get<1>();
1284 auto const elt_en = lit->template get<2>();
1285 DLOG(INFO) <<
"integrating over " << std::distance( elt_it, elt_en ) <<
" elements\n";
1288 if ( elt_it == elt_en )
1291 const size_type idEltRangeInit = elt_it->id();
1292 size_type idEltTestInit = idEltRangeInit;
1293 if ( elt_it->mesh()->isSubMeshFrom( __form.testSpace()->mesh() ) )
1294 idEltTestInit = elt_it->mesh()->subMeshToMesh( idEltRangeInit );
1295 else if ( __form.testSpace()->mesh()->isSubMeshFrom( elt_it->mesh() ) )
1296 idEltTestInit = __form.testSpace()->mesh()->meshToSubMesh( idEltRangeInit );
1297 size_type idEltTrialInit = idEltRangeInit;
1298 if ( elt_it->mesh()->isSubMeshFrom( __form.trialSpace()->mesh() ) )
1299 idEltTrialInit = elt_it->mesh()->subMeshToMesh( idEltRangeInit );
1300 else if ( __form.trialSpace()->mesh()->isSubMeshFrom( elt_it->mesh() ) )
1301 idEltTrialInit = __form.trialSpace()->mesh()->meshToSubMesh( idEltRangeInit );
1303 pc_expr_ptrtype geopcExpr(
new pc_expr_type( elt_it->gm(), this->im().points() ) );
1304 gmc_expr_ptrtype gmcExpr(
new gmc_expr_type( elt_it->gm(),*elt_it, geopcExpr ) );
1306 auto gmcFormTest = detail::buildGmcWithRelationDifferentMeshType<
typename FormType::space_1_type,im_formtest_type,
1307 gmc_formTest_type,gmc_expr_type>( __form.testSpace(), __form.testSpace()->gm(),
1308 idEltTestInit, gmcExpr, mpl::int_<gmTestRangeRelation>() );
1309 auto gmcFormTrial = detail::buildGmcWithRelationDifferentMeshType<
typename FormType::space_2_type,im_formtrial_type,
1310 gmc_formTrial_type,gmc_expr_type>( __form.trialSpace(), __form.trialSpace()->gm(),
1311 idEltTrialInit, gmcExpr, mpl::int_<gmTrialRangeRelation>() );
1313 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
1314 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTest ) );
1315 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTrial ) );
1317 im_formtest_type imTest;
1318 im_formtrial_type imTrial;
1319 focb_ptrtype formc(
new form_context_type( __form,
1324 this->im(),imTest,imTrial ) );
1328 for ( ; elt_it != elt_en; ++elt_it )
1330 const size_type idEltRange = elt_it->id();
1332 if ( elt_it->mesh()->isSubMeshFrom( __form.testSpace()->mesh() ) )
1333 idEltTest = elt_it->mesh()->subMeshToMesh( idEltRange );
1334 else if ( __form.testSpace()->mesh()->isSubMeshFrom( elt_it->mesh() ) )
1335 idEltTest = __form.testSpace()->mesh()->meshToSubMesh( idEltRange );
1337 if ( elt_it->mesh()->isSubMeshFrom( __form.trialSpace()->mesh() ) )
1338 idEltTrial = elt_it->mesh()->subMeshToMesh( idEltRange );
1339 else if ( __form.trialSpace()->mesh()->isSubMeshFrom( elt_it->mesh() ) )
1340 idEltTrial = __form.trialSpace()->mesh()->meshToSubMesh( idEltRange );
1343 if ( formc->isZero( eltTest.element0() ) )
1347 gmcExpr->update(*elt_it);
1348 detail::updateGmcWithRelationDifferentMeshType<
typename FormType::space_1_type,im_formtest_type,
1349 gmc_formTest_type,gmc_expr_type>(__form.testSpace(), gmcFormTest, gmcExpr,
1350 idEltTest, mpl::int_<gmTestRangeRelation>() );
1351 detail::updateGmcWithRelationDifferentMeshType<
typename FormType::space_2_type,im_formtrial_type,
1352 gmc_formTrial_type,gmc_expr_type>(__form.trialSpace(), gmcFormTrial, gmcExpr,
1353 idEltTrial, mpl::int_<gmTrialRangeRelation>() );
1355 formc->update( mapgmcFormTest,mapgmcFormTrial,mapgmcExpr );
1365 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1366 template<
typename FE,
typename VectorType,
typename ElemContType>
1368 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const
1370 CHECK (
false ) <<
"[assembleWithRelationDifferentMeshType<LinearForm,MESH_ELEMENTS>] : not implement\n";
1373 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1374 template<
typename FE1,
typename FE2,
typename ElemContType>
1376 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1377 mpl::int_<MESH_ELEMENTS> )
const
1379 LOG(INFO) <<
"[integrator::assembleInCaseOfInterpolate] vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS>\n";
1382 typedef typename eval::gm_type gm_expr_type;
1383 typedef typename gm_expr_type::template Context<expression_type::context|vm::POINT|vm::JACOBIAN, typename eval::element_type> gmc_expr_type;
1384 typedef boost::shared_ptr<gmc_expr_type> gmc_expr_ptrtype;
1385 typedef typename gm_expr_type::precompute_type pc_expr_type;
1386 typedef typename gm_expr_type::precompute_ptrtype pc_expr_ptrtype;
1387 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_expr_ptrtype> > map_gmc_expr_type;
1390 typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
1392 typedef typename FormType::gm_1_type gm_formTest_type;
1393 typedef typename FormType::mesh_element_1_type geoelement_formTest_type;
1394 typedef typename gm_formTest_type::template Context<expression_type::context|vm::POINT,geoelement_formTest_type> gmc_formTest_type;
1395 typedef boost::shared_ptr<gmc_formTest_type> gmc_formTest_ptrtype;
1396 typedef typename gm_formTest_type::precompute_type pc_formTest_type;
1397 typedef typename gm_formTest_type::precompute_ptrtype pc_formTest_ptrtype;
1398 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTest_ptrtype> > map_gmc_formTest_type;
1400 typedef typename FormType::gm_2_type gm_formTrial_type;
1401 typedef typename FormType::mesh_element_2_type geoelement_formTrial_type;
1402 typedef typename gm_formTrial_type::template Context<expression_type::context|vm::POINT,geoelement_formTrial_type> gmc_formTrial_type;
1403 typedef boost::shared_ptr<gmc_formTrial_type> gmc_formTrial_ptrtype;
1404 typedef typename gm_formTrial_type::precompute_type pc_formTrial_type;
1405 typedef typename gm_formTrial_type::precompute_ptrtype pc_formTrial_ptrtype;
1406 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTrial_ptrtype> > map_gmc_formTrial_type;
1409 typedef typename FormType::template Context<map_gmc_formTest_type, expression_type, im_type,map_gmc_expr_type, map_gmc_formTrial_type> form_context_type;
1410 typedef form_context_type fcb_type;
1411 typedef fcb_type* focb_ptrtype;
1415 element_iterator elt_it, elt_en;
1416 bool findEltForInit =
false;
1417 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
1419 elt_it = lit->template get<1>();
1420 elt_en = lit->template get<2>();
1423 if ( elt_it == elt_en )
1426 findEltForInit =
true;
1428 if (!findEltForInit)
return;
1431 pc_expr_ptrtype geopcExpr(
new pc_expr_type( elt_it->gm(), this->im().points() ) );
1432 gmc_expr_ptrtype gmcExpr(
new gmc_expr_type( elt_it->gm(),*elt_it, geopcExpr ) );
1433 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
1435 pc_formTest_ptrtype geopcFormTest(
new pc_formTest_type( __form.gm(), __form.testSpace()->fe()->points() ) );
1436 gmc_formTest_ptrtype gmcFormTest(
new gmc_formTest_type( __form.gm(), __form.testSpace()->mesh()->element( 0 ), geopcFormTest ) );
1437 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTest ) );
1439 pc_formTrial_ptrtype geopcFormTrial(
new pc_formTrial_type( __form.gmTrial(), __form.trialSpace()->fe()->points() ) );
1440 gmc_formTrial_ptrtype gmcFormTrial(
new gmc_formTrial_type( __form.gmTrial(), __form.trialSpace()->mesh()->element( 0 ), geopcFormTrial ) );
1441 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTrial ) );
1444 focb_ptrtype formc(
new form_context_type( __form,
1453 auto meshTrial = __form.trialSpace()->mesh();
1454 auto meshTest = __form.testSpace()->mesh();
1458 M_QPL.reset(
new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
1459 M_QPL->update( meshTest,meshTrial );
1462 if (!M_QPL->hasPrecompute())
1464 auto res_it = M_QPL->result().begin();
1465 auto const res_en = M_QPL->result().end();
1466 for ( ; res_it != res_en ; ++res_it )
1468 auto const idEltTest = res_it->template get<0>();
1469 auto const& map = res_it->template get<1>();
1470 auto map_it = map.begin();
1471 auto const map_en = map.end();
1472 for ( ; map_it != map_en ; ++map_it )
1474 auto const idEltTrial = map_it->first;
1475 auto const& eltTrial = meshTrial->element( idEltTrial );
1476 auto const& eltTest = meshTest->element( idEltTest );
1478 auto const& ptRefTest = map_it->second.template get<1>();
1479 auto const& ptRefTrial = map_it->second.template get<2>();
1480 auto const& themapQuad = map_it->second.template get<0>();
1481 auto vec_gmcExpr = M_QPL->getUsableDataInFormContext( themapQuad,ptRefTest,ptRefTrial );
1482 auto gmcExpr_it = vec_gmcExpr.begin();
1483 auto const gmcExpr_en = vec_gmcExpr.end();
1484 bool isFirstExperience =
true;
1485 for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
1487 geopcFormTest->update( gmcExpr_it->template get<2>() );
1488 geopcFormTrial->update( gmcExpr_it->template get<3>() );
1490 gmcFormTest->update( eltTest,geopcFormTest );
1491 gmcFormTrial->update( eltTrial,geopcFormTrial );
1493 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTest ) );
1494 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTrial ) );
1495 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr_it->template get<1>() ) );
1497 formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr, gmcExpr_it->template get<0>() );
1499 formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
1500 isFirstExperience =
false;
1503 formc->assembleInCaseOfInterpolate();
1510 auto const& resQPL = M_QPL->getPrecompute(__form);
1511 auto resQPL_it = resQPL.begin();
1512 auto const resQPL_en = resQPL.end();
1513 for ( ; resQPL_it != resQPL_en ; ++resQPL_it)
1515 auto resQPLloc_it = resQPL_it->begin();
1516 auto const resQPLloc_en = resQPL_it->end();
1517 bool isFirstExperience =
true;
1518 for ( ; resQPLloc_it != resQPLloc_en ; ++resQPLloc_it)
1520 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>() ) );
1521 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<3>() ) );
1522 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<1>() ) );
1524 formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr, resQPLloc_it->template get<0>() );
1526 formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
1527 isFirstExperience =
false;
1530 formc->assembleInCaseOfInterpolate();
1540 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1541 template<
typename FE,
typename VectorType,
typename ElemContType>
1543 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> )
const
1547 typedef typename eval::gm_type gm_expr_type;
1548 typedef typename gm_expr_type::template Context<expression_type::context|vm::POINT|vm::JACOBIAN, typename eval::element_type> gmc_expr_type;
1549 typedef boost::shared_ptr<gmc_expr_type> gmc_expr_ptrtype;
1550 typedef typename gm_expr_type::precompute_type pc_expr_type;
1551 typedef typename gm_expr_type::precompute_ptrtype pc_expr_ptrtype;
1552 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_expr_ptrtype> > map_gmc_expr_type;
1555 typedef vf::detail::LinearForm<FE,VectorType,ElemContType> FormType;
1556 typedef typename FormType::gm_type gm_form_type;
1557 typedef typename FormType::mesh_test_element_type geoelement_form_type;
1558 typedef typename gm_form_type::template Context<expression_type::context|vm::POINT|vm::JACOBIAN,geoelement_form_type> gmc_form_type;
1559 typedef boost::shared_ptr<gmc_form_type> gmc_form_ptrtype;
1560 typedef typename gm_form_type::precompute_type pc_form_type;
1561 typedef typename gm_form_type::precompute_ptrtype pc_form_ptrtype;
1562 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_form_ptrtype> > map_gmc_form_type;
1565 typedef typename FormType::template Context<map_gmc_form_type, expression_type, im_type,map_gmc_expr_type> form_context_type;
1566 typedef form_context_type fcb_type;
1567 typedef fcb_type* focb_ptrtype;
1571 element_iterator elt_it, elt_en;
1572 bool findEltForInit =
false;
1573 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
1575 elt_it = lit->template get<1>();
1576 elt_en = lit->template get<2>();
1579 if ( elt_it == elt_en )
1582 findEltForInit =
true;
1584 if (!findEltForInit)
return;
1588 pc_expr_ptrtype geopcExpr(
new pc_expr_type( elt_it->gm(), this->im().points() ) );
1589 gmc_expr_ptrtype gmcExpr(
new gmc_expr_type( elt_it->gm(),*elt_it, geopcExpr ) );
1590 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
1594 pc_form_ptrtype geopcForm(
new pc_form_type( __form.gm(), __form.testSpace()->fe()->points() ) );
1595 gmc_form_ptrtype gmcForm(
new gmc_form_type( __form.gm(), __form.testSpace()->mesh()->element( 0 ), geopcForm ) );
1596 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( gmcForm ) );
1600 focb_ptrtype formc(
new form_context_type( __form,
1609 auto meshTest = __form.testSpace()->mesh();
1613 M_QPL.reset(
new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
1614 M_QPL->update( meshTest );
1619 if (!M_QPL->hasPrecompute())
1621 auto res_it = M_QPL->resultLinear().begin();
1622 auto const res_en = M_QPL->resultLinear().end();
1624 for ( ; res_it != res_en ; ++res_it )
1626 auto const idEltTest = res_it->template get<0>();
1627 auto const& eltTest = meshTest->element( idEltTest );
1628 auto const& ptRefTest = res_it->template get<2>();
1629 auto const& themapQuad = res_it->template get<1>();
1630 auto vec_gmcExpr = M_QPL->getUsableDataInFormContext( themapQuad,ptRefTest );
1631 auto gmcExpr_it = vec_gmcExpr.begin();
1632 auto const gmcExpr_en = vec_gmcExpr.end();
1633 bool isFirstExperience =
true;
1635 for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
1637 geopcForm->update( gmcExpr_it->template get<2>() );
1638 gmcForm->update( eltTest,geopcForm );
1639 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( gmcForm ) );
1640 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr_it->template get<1>() ) );
1641 formc->updateInCaseOfInterpolate( mapgmcForm, mapgmcForm, mapgmcExpr,gmcExpr_it->template get<0>() );
1643 formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
1644 isFirstExperience =
false;
1652 auto const& resQPL = M_QPL->getPrecompute(__form);
1653 auto resQPL_it = resQPL.begin();
1654 auto const resQPL_en = resQPL.end();
1655 for ( ; resQPL_it != resQPL_en ; ++resQPL_it)
1657 auto resQPLloc_it = resQPL_it->begin();
1658 auto const resQPLloc_en = resQPL_it->end();
1659 bool isFirstExperience =
true;
1660 for ( ; resQPLloc_it != resQPLloc_en ; ++resQPLloc_it)
1662 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>() ) );
1663 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<1>() ) );
1664 formc->updateInCaseOfInterpolate( mapgmcForm, mapgmcForm, mapgmcExpr,resQPLloc_it->template get<0>() );
1666 formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
1667 isFirstExperience =
false;
1678 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1679 template<
typename FormType>
1681 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<true> ,
bool )
const
1683 DLOG(INFO) <<
"integrating over "
1684 << std::distance( this->beginElement(), this->endElement() ) <<
" faces\n";
1685 boost::timer __timer;
1690 typedef typename eval::gm_type gm_type;
1691 typedef typename eval::gm1_type gm1_type;
1693 typedef boost::shared_ptr<gm_type> gm_ptrtype;
1694 typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
1695 typedef typename gm_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, typename eval::element_type> gmc_type;
1696 typedef typename gm1_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, typename eval::element_type> gmc1_type;
1697 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1698 typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
1703 typedef typename eval::gmpc_type pc_type;
1704 typedef typename eval::gmpc1_type pc1_type;
1705 typedef typename eval::gmpc_ptrtype pc_ptrtype;
1706 typedef typename eval::gmpc1_ptrtype pc1_ptrtype;
1709 QuadMapped<im_type> qm;
1710 typedef typename QuadMapped<im_type>::permutation_type permutation_type;
1711 typename QuadMapped<im_type>::permutation_points_type ppts( qm( im() ) );
1714 std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( im().nFaces() );
1715 std::vector<std::map<permutation_type, pc1_ptrtype> > __geopc1( im().nFaces() );
1716 typedef typename im_type::face_quadrature_type face_im_type;
1717 typedef typename im2_type::face_quadrature_type face_im2_type;
1720 std::vector<face_im_type> face_ims( im().nFaces() );
1721 std::vector<face_im2_type> face_ims2( im2().nFaces() );
1723 for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
1725 face_ims[__f] = this->im( __f );
1726 face_ims2[__f] = this->im2( __f );
1728 for ( permutation_type __p( permutation_type::IDENTITY );
1729 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1732 __geopc[__f][__p] = pc_ptrtype(
new pc_type( __form.gm(), ppts[__f].find( __p )->second ) );
1733 __geopc1[__f][__p] = pc1_ptrtype(
new pc1_type( __form.gm1(), ppts[__f].find( __p )->second ) );
1737 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
1739 auto it = lit->template get<1>();
1740 auto en = lit->template get<2>();
1743 if ( (it == en) || (it->isConnectedTo0() ==
false) )
1746 uint16_type __face_id_in_elt_0 = it->pos_first();
1750 gm_ptrtype __gm = it->element( 0 ).gm();
1751 gm1_ptrtype __gm1 = it->element( 0 ).gm1();
1753 gmc_ptrtype __c0(
new gmc_type( __gm, it->element( 0 ), __geopc, __face_id_in_elt_0 ) );
1754 gmc1_ptrtype __c01(
new gmc1_type( __gm1, it->element( 0 ), __geopc1, __face_id_in_elt_0 ) );
1761 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
1762 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
1763 typedef typename FormType::template Context<map_gmc_type, expression_type, face_im_type> form_context_type;
1764 typedef typename FormType::template Context<map_gmc1_type, expression_type, face_im2_type> form1_context_type;
1765 typedef boost::shared_ptr<form_context_type> form_context_ptrtype;
1766 typedef boost::shared_ptr<form1_context_type> form1_context_ptrtype;
1767 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
1768 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c01 ) );
1769 form_context_ptrtype form;
1770 form1_context_ptrtype form1;
1779 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype>, fusion::pair<vf::detail::gmc<1>, gmc_ptrtype> > map2_gmc_type;
1780 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype>, fusion::pair<vf::detail::gmc<1>, gmc1_ptrtype> > map21_gmc_type;
1781 typedef typename FormType::template Context<map2_gmc_type, expression_type, face_im_type> form2_context_type;
1782 typedef typename FormType::template Context<map21_gmc_type, expression_type, face_im2_type> form21_context_type;
1783 typedef boost::shared_ptr<form2_context_type> form2_context_ptrtype;
1784 typedef boost::shared_ptr<form21_context_type> form21_context_ptrtype;
1785 form2_context_ptrtype form2;
1786 form21_context_ptrtype form21;
1788 bool isInitConnectionTo0=
false;
1789 bool isInitConnectionTo1=
false;
1792 if ( it->isConnectedTo1() )
1794 uint16_type __face_id_in_elt_1 = it->pos_second();
1796 __c1 = gmc_ptrtype(
new gmc_type( __gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
1797 __c11 = gmc1_ptrtype(
new gmc1_type( __gm1, it->element( 1 ), __geopc1, __face_id_in_elt_1 ) );
1798 map2_gmc_type mapgmc2( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
1799 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
1800 map21_gmc_type mapgmc21( fusion::make_pair<vf::detail::gmc<0> >( __c01 ),
1801 fusion::make_pair<vf::detail::gmc<1> >( __c11 ) );
1803 form2 = form2_context_ptrtype(
new form2_context_type( __form, mapgmc2, mapgmc2, mapgmc2, expression(), face_ims[__face_id_in_elt_0], this->im(), mpl::int_<2>() ) );
1804 form21 = form21_context_ptrtype(
new form21_context_type( __form, mapgmc21, mapgmc21, mapgmc21, expression(), face_ims2[__face_id_in_elt_0], this->im2(), mpl::int_<2>() ) );
1805 isInitConnectionTo1=
true;
1810 form = form_context_ptrtype(
new form_context_type( __form, mapgmc, mapgmc, mapgmc, expression(), face_ims[__face_id_in_elt_0], this->im() ) );
1811 form1 = form1_context_ptrtype(
new form1_context_type( __form, mapgmc1, mapgmc1, mapgmc1, expression(), face_ims2[__face_id_in_elt_0], this->im2() ) );
1812 isInitConnectionTo0=
true;
1815 boost::timer ti0,ti1, ti2, ti3;
1817 DLOG(INFO) <<
"[Integrator::faces/forms] starting...\n";
1826 for ( ; it != en; ++it )
1828 if ( it->isGhostFace() )
1830 LOG(WARNING) <<
"face id : " << it->id() <<
" is a ghost face";
1833 if ( it->isConnectedTo1() )
1835 if ( !isInitConnectionTo1 )
1837 uint16_type __face_id_in_elt_1 = it->pos_second();
1839 __c1 = gmc_ptrtype(
new gmc_type( __gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
1840 __c11 = gmc1_ptrtype(
new gmc1_type( __gm1, it->element( 1 ), __geopc1, __face_id_in_elt_1 ) );
1841 map2_gmc_type mapgmc2( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
1842 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
1843 map21_gmc_type mapgmc21( fusion::make_pair<vf::detail::gmc<0> >( __c01 ),
1844 fusion::make_pair<vf::detail::gmc<1> >( __c11 ) );
1846 form2 = form2_context_ptrtype(
new form2_context_type( __form, mapgmc2, mapgmc2, mapgmc2, expression(), face_ims[__face_id_in_elt_0], this->im(), mpl::int_<2>() ) );
1847 form21 = form21_context_ptrtype(
new form21_context_type( __form, mapgmc21, mapgmc21, mapgmc21, expression(), face_ims2[__face_id_in_elt_0], this->im2(), mpl::int_<2>() ) );
1848 isInitConnectionTo1=
true;
1854 case GeomapStrategyType::GEOMAP_HO:
1856 FEELPP_ASSERT( it->isOnBoundary() == false )
1857 ( it->id() ).error(
"face on boundary but connected on both sides" );
1860 uint16_type __face_id_in_elt_0 = it->pos_first();
1861 uint16_type __face_id_in_elt_1 = it->pos_second();
1863 __c0->update( it->element( 0 ), __face_id_in_elt_0 );
1864 __c1->update( it->element( 1 ), __face_id_in_elt_1 );
1868 map2_gmc_type mapgmc2 = map2_gmc_type( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
1869 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
1870 form2->update( mapgmc2, mapgmc2, mapgmc2, face_ims[__face_id_in_elt_0], mpl::int_<2>() );
1874 form2->integrate( );
1878 form2->assemble( it->element( 0 ).id(), it->element( 1 ).id() );
1883 case GeomapStrategyType::GEOMAP_O1:
1884 case GeomapStrategyType::GEOMAP_OPT:
1886 FEELPP_ASSERT( it->isOnBoundary() == false )
1887 ( it->id() ).error(
"face on boundary but connected on both sides" );
1890 uint16_type __face_id_in_elt_0 = it->pos_first();
1891 uint16_type __face_id_in_elt_1 = it->pos_second();
1893 __c01->update( it->element( 0 ), __face_id_in_elt_0 );
1894 __c11->update( it->element( 1 ), __face_id_in_elt_1 );
1898 map21_gmc_type mapgmc21 = map21_gmc_type( fusion::make_pair<vf::detail::gmc<0> >( __c01 ),
1899 fusion::make_pair<vf::detail::gmc<1> >( __c11 ) );
1900 form21->update( mapgmc21, mapgmc21, mapgmc21, face_ims2[__face_id_in_elt_0], mpl::int_<2>() );
1904 form21->integrate( );
1908 form21->assemble( it->element( 0 ).id(), it->element( 1 ).id() );
1918 uint16_type __face_id_in_elt_0 = it->pos_first();
1920 if ( !isInitConnectionTo0 )
1922 form = form_context_ptrtype(
new form_context_type( __form, mapgmc, mapgmc, mapgmc, expression(), face_ims[__face_id_in_elt_0], this->im() ) );
1923 form1 = form1_context_ptrtype(
new form1_context_type( __form, mapgmc1, mapgmc1, mapgmc1, expression(), face_ims2[__face_id_in_elt_0], this->im2() ) );
1924 isInitConnectionTo0=
true;
1928 __c0->update( it->element( 0 ),__face_id_in_elt_0 );
1931 FEELPP_ASSERT( __face_id_in_elt_0 == __c0->faceId() )
1932 ( __face_id_in_elt_0 )
1933 ( __c0->faceId() ).warn (
"invalid face id" );
1936 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
1937 form->update( mapgmc, mapgmc, mapgmc, face_ims[__face_id_in_elt_0] );
1945 form->assemble( it->element( 0 ).id() );
1948 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
1949 form->update( mapgmc, mapgmc, mapgmc, face_ims[__face_id_in_elt_0] );
1957 DLOG(INFO) <<
"[faces] Overall integration time : " << ( t0+t1+t2+t3 ) <<
" per element:" << ( t0+t1+t2+t3 )/std::distance( this->beginElement(), this->endElement() ) <<
"for " << std::distance( this->beginElement(), this->endElement() ) <<
"elements\n";
1958 DLOG(INFO) <<
"[faces] Overall geometric mapping update time : " << t0 <<
"\n";
1959 DLOG(INFO) <<
"[faces] Overall form update time : " << t1 <<
"\n";
1960 DLOG(INFO) <<
"[faces] Overall local assembly time : " << t2 <<
"\n";
1961 DLOG(INFO) <<
"[faces] Overall global assembly time : " << t3 <<
"\n";
1963 DLOG(INFO) <<
"integrating over faces done in " << __timer.elapsed() <<
"s\n";
1967 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1968 template<
typename FormType>
1970 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<false> ,
bool hasRelation )
const
1974 assembleWithRelationDifferentMeshType( __form,mpl::int_<MESH_FACES>() );
1978 assembleInCaseOfInterpolate( __form,mpl::int_<MESH_FACES>() );
1982 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1983 template<
typename FE1,
typename FE2,
typename ElemContType>
1985 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1986 mpl::int_<MESH_FACES> )
const
1988 CHECK (
false ) <<
"[assembleWithRelationDifferentMeshType<BilinearForm,MESH_FACES>] : not implement\n";
1991 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
1992 template<
typename FE,
typename VectorType,
typename ElemContType>
1994 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> )
const
1996 CHECK (
false ) <<
"[assembleWithRelationDifferentMeshType<LinearForm,MESH_FACES>] : not implement\n";
1999 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2000 template<
typename FE1,
typename FE2,
typename ElemContType>
2002 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
2003 mpl::int_<MESH_FACES> )
const
2007 typedef typename eval::gm_type gm_expr_type;
2008 typedef typename gm_expr_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, typename eval::element_type> gmc_expr_type;
2009 typedef boost::shared_ptr<gmc_expr_type> gmc_expr_ptrtype;
2010 typedef typename gm_expr_type::precompute_type pc_expr_type;
2011 typedef typename gm_expr_type::precompute_ptrtype pc_expr_ptrtype;
2012 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_expr_ptrtype> > map_gmc_expr_type;
2015 typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
2024 typedef typename FormType::gm_1_type gm_formTest_type;
2025 typedef typename FormType::mesh_element_1_type geoelement_formTest_type;
2026 typedef typename gm_formTest_type::template Context<expression_type::context|vm::POINT,geoelement_formTest_type> gmc_formTest_type;
2027 typedef boost::shared_ptr<gmc_formTest_type> gmc_formTest_ptrtype;
2028 typedef typename gm_formTest_type::precompute_type pc_formTest_type;
2029 typedef typename gm_formTest_type::precompute_ptrtype pc_formTest_ptrtype;
2030 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTest_ptrtype> > map_gmc_formTest_type;
2032 typedef typename FormType::gm_2_type gm_formTrial_type;
2033 typedef typename FormType::mesh_element_2_type geoelement_formTrial_type;
2034 typedef typename gm_formTrial_type::template Context<expression_type::context|vm::POINT,geoelement_formTrial_type> gmc_formTrial_type;
2035 typedef boost::shared_ptr<gmc_formTrial_type> gmc_formTrial_ptrtype;
2036 typedef typename gm_formTrial_type::precompute_type pc_formTrial_type;
2037 typedef typename gm_formTrial_type::precompute_ptrtype pc_formTrial_ptrtype;
2038 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_formTrial_ptrtype> > map_gmc_formTrial_type;
2043 typedef typename im_type::face_quadrature_type face_im_type;
2045 typedef typename FormType::template Context<map_gmc_formTest_type, expression_type, face_im_type,map_gmc_expr_type, map_gmc_formTrial_type> form_context_type;
2046 typedef form_context_type fcb_type;
2047 typedef fcb_type* focb_ptrtype;
2050 element_iterator elt_it, elt_en;
2051 bool findEltForInit =
false;
2052 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
2054 elt_it = lit->template get<1>();
2055 elt_en = lit->template get<2>();
2058 if ( elt_it == elt_en )
2061 findEltForInit =
true;
2063 if (!findEltForInit)
return;
2066 QuadMapped<im_type> qm;
2067 typedef typename QuadMapped<im_type>::permutation_type permutation_type;
2068 typename QuadMapped<im_type>::permutation_points_type ppts( qm( im() ) );
2070 std::vector<std::map<permutation_type, pc_expr_ptrtype> > __geopcExpr( im().nFaces() );
2071 std::vector<face_im_type> face_ims( im().nFaces() );
2073 for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2075 face_ims[__f] = this->im( __f );
2077 for ( permutation_type __p( permutation_type::IDENTITY );
2078 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2081 __geopcExpr[__f][__p] = pc_expr_ptrtype(
new pc_expr_type( elt_it->element( 0 ).gm(), ppts[__f].find( __p )->second ) );
2086 uint16_type __face_id_in_elt_0 = elt_it->pos_first();
2088 gmc_expr_ptrtype gmcExpr(
new gmc_expr_type( elt_it->element( 0 ).gm(),
2089 elt_it->element( 0 ),
2091 __face_id_in_elt_0 ) );
2094 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
2098 pc_formTest_ptrtype geopcFormTest(
new pc_formTest_type( __form.gm(), __form.testSpace()->fe()->points() ) );
2099 gmc_formTest_ptrtype gmcFormTest(
new gmc_formTest_type( __form.gm(), __form.testSpace()->mesh()->element( 0 ), geopcFormTest ) );
2100 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTest ) );
2102 pc_formTrial_ptrtype geopcFormTrial(
new pc_formTrial_type( __form.gmTrial(), __form.trialSpace()->fe()->points() ) );
2103 gmc_formTrial_ptrtype gmcFormTrial(
new gmc_formTrial_type( __form.gmTrial(), __form.trialSpace()->mesh()->element( 0 ), geopcFormTrial ) );
2104 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTrial ) );
2108 focb_ptrtype formc(
new form_context_type( __form,
2113 face_ims[__face_id_in_elt_0],
2118 auto meshTrial = __form.trialSpace()->mesh();
2119 auto meshTest = __form.testSpace()->mesh();
2123 M_QPL.reset(
new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
2124 M_QPL->update( meshTest,meshTrial );
2129 if (!M_QPL->hasPrecompute())
2131 auto res_it = M_QPL->result().begin();
2132 auto res_en = M_QPL->result().end();
2134 for ( ; res_it != res_en ; ++res_it )
2136 auto const idEltTest = res_it->template get<0>();
2137 auto const& map = res_it->template get<1>();
2138 auto map_it = map.begin();
2139 auto const map_en = map.end();
2141 for ( ; map_it != map_en ; ++map_it )
2143 auto const idEltTrial = map_it->first;
2144 auto const& eltTrial = meshTrial->element( idEltTrial );
2145 auto const& eltTest = meshTest->element( idEltTest );
2146 auto const& ptRefTest = map_it->second.template get<1>();
2147 auto const& ptRefTrial = map_it->second.template get<2>();
2148 auto const& themapQuad = map_it->second.template get<0>();
2150 auto vec_gmcExpr = M_QPL->getUsableDataInFormContext( themapQuad,ptRefTest,ptRefTrial );
2151 auto gmcExpr_it = vec_gmcExpr.begin();
2152 auto const gmcExpr_en = vec_gmcExpr.end();
2153 bool isFirstExperience =
true;
2155 for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
2157 geopcFormTest->update( gmcExpr_it->template get<2>() );
2158 geopcFormTrial->update( gmcExpr_it->template get<3>() );
2160 gmcFormTest->update( eltTest,geopcFormTest );
2161 gmcFormTrial->update( eltTrial,geopcFormTrial );
2164 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTest ) );
2165 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( gmcFormTrial ) );
2166 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr_it->template get<1>() ) );
2167 __face_id_in_elt_0 = gmcExpr_it->template get<1>()->faceId();
2168 formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr,face_ims[__face_id_in_elt_0],gmcExpr_it->template get<0>() );
2170 formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
2171 isFirstExperience =
false;
2174 formc->assembleInCaseOfInterpolate();
2181 auto const& resQPL = M_QPL->getPrecompute(__form);
2182 auto resQPL_it = resQPL.begin();
2183 auto const resQPL_en = resQPL.end();
2184 for ( ; resQPL_it != resQPL_en ; ++resQPL_it)
2186 auto resQPLloc_it = resQPL_it->begin();
2187 auto const resQPLloc_en = resQPL_it->end();
2188 bool isFirstExperience =
true;
2189 for ( ; resQPLloc_it != resQPLloc_en ; ++resQPLloc_it)
2191 map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>() ) );
2192 map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<3>() ) );
2193 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<1>() ) );
2194 __face_id_in_elt_0 = resQPLloc_it->template get<1>()->faceId();
2195 formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr,face_ims[__face_id_in_elt_0],resQPLloc_it->template get<0>() );
2197 formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
2198 isFirstExperience =
false;
2201 formc->assembleInCaseOfInterpolate();
2208 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2209 template<
typename FE,
typename VectorType,
typename ElemContType>
2211 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> )
const
2214 typedef typename eval::gm_type gm_expr_type;
2215 typedef typename gm_expr_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, typename eval::element_type> gmc_expr_type;
2216 typedef boost::shared_ptr<gmc_expr_type> gmc_expr_ptrtype;
2217 typedef typename gm_expr_type::precompute_type pc_expr_type;
2218 typedef typename gm_expr_type::precompute_ptrtype pc_expr_ptrtype;
2219 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_expr_ptrtype> > map_gmc_expr_type;
2222 typedef vf::detail::LinearForm<FE,VectorType,ElemContType> FormType;
2223 typedef typename FormType::gm_type gm_form_type;
2224 typedef typename FormType::mesh_test_element_type geoelement_form_type;
2225 typedef typename gm_form_type::template Context<expression_type::context|vm::POINT,geoelement_form_type> gmc_form_type;
2226 typedef boost::shared_ptr<gmc_form_type> gmc_form_ptrtype;
2227 typedef typename gm_form_type::precompute_type pc_form_type;
2228 typedef typename gm_form_type::precompute_ptrtype pc_form_ptrtype;
2229 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_form_ptrtype> > map_gmc_form_type;
2232 typedef typename im_type::face_quadrature_type face_im_type;
2234 typedef typename FormType::template Context<map_gmc_form_type, expression_type, face_im_type,map_gmc_expr_type> form_context_type;
2235 typedef form_context_type fcb_type;
2236 typedef fcb_type* focb_ptrtype;
2239 element_iterator elt_it, elt_en;
2240 bool findEltForInit =
false;
2241 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
2243 elt_it = lit->template get<1>();
2244 elt_en = lit->template get<2>();
2247 if ( elt_it == elt_en )
2250 findEltForInit =
true;
2252 if (!findEltForInit)
return;
2256 QuadMapped<im_type> qm;
2257 typedef typename QuadMapped<im_type>::permutation_type permutation_type;
2258 typename QuadMapped<im_type>::permutation_points_type ppts( qm( im() ) );
2260 std::vector<std::map<permutation_type, pc_expr_ptrtype> > __geopcExpr( im().nFaces() );
2261 std::vector<face_im_type> face_ims( im().nFaces() );
2263 for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2265 face_ims[__f] = this->im( __f );
2267 for ( permutation_type __p( permutation_type::IDENTITY );
2268 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2271 __geopcExpr[__f][__p] = pc_expr_ptrtype(
new pc_expr_type( elt_it->element( 0 ).gm(), ppts[__f].find( __p )->second ) );
2276 uint16_type __face_id_in_elt_0 = elt_it->pos_first();
2278 gmc_expr_ptrtype gmcExpr(
new gmc_expr_type( elt_it->element( 0 ).gm(),
2279 elt_it->element( 0 ),
2281 __face_id_in_elt_0 ) );
2283 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
2287 pc_form_ptrtype geopcForm(
new pc_form_type( __form.gm(), __form.testSpace()->fe()->points() ) );
2288 gmc_form_ptrtype gmcForm(
new gmc_form_type( __form.gm(), __form.testSpace()->mesh()->element( 0 ), geopcForm ) );
2289 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( gmcForm ) );
2293 focb_ptrtype formc(
new form_context_type( __form,
2298 face_ims[__face_id_in_elt_0],
2303 auto meshTest = __form.testSpace()->mesh();
2307 M_QPL.reset(
new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
2308 M_QPL->update( meshTest );
2313 if (!M_QPL->hasPrecompute())
2315 auto res_it = M_QPL->resultLinear().begin();
2316 auto const res_en = M_QPL->resultLinear().end();
2317 for ( ; res_it != res_en ; ++res_it )
2319 auto const idEltTest = res_it->template get<0>();
2320 auto const& eltTest = meshTest->element( idEltTest );
2321 auto const& ptRefTest = res_it->template get<2>();
2322 auto const& themapQuad = res_it->template get<1>();
2324 auto vec_gmcExpr = M_QPL->getUsableDataInFormContext( themapQuad,ptRefTest );
2325 auto gmcExpr_it = vec_gmcExpr.begin();
2326 auto const gmcExpr_en = vec_gmcExpr.end();
2327 bool isFirstExperience =
true;
2328 for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
2330 geopcForm->update( gmcExpr_it->template get<2>() );
2331 gmcForm->update( eltTest,geopcForm );
2332 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( gmcForm ) );
2333 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr_it->template get<1>() ) );
2334 __face_id_in_elt_0 = gmcExpr_it->template get<1>()->faceId();
2335 formc->updateInCaseOfInterpolate( mapgmcForm, mapgmcForm, mapgmcExpr,face_ims[__face_id_in_elt_0],gmcExpr_it->template get<0>() );
2337 formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
2338 isFirstExperience =
false;
2346 auto const& resQPL = M_QPL->getPrecompute(__form);
2347 auto resQPL_it = resQPL.begin();
2348 auto const resQPL_en = resQPL.end();
2349 for ( ; resQPL_it != resQPL_en ; ++resQPL_it)
2351 auto resQPLloc_it = resQPL_it->begin();
2352 auto const resQPLloc_en = resQPL_it->end();
2353 bool isFirstExperience =
true;
2354 for ( ; resQPLloc_it != resQPLloc_en ; ++resQPLloc_it)
2356 map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>() ) );
2357 map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<1>() ) );
2358 __face_id_in_elt_0 = resQPLloc_it->template get<1>()->faceId();
2359 formc->updateInCaseOfInterpolate( mapgmcForm, mapgmcForm, mapgmcExpr,face_ims[__face_id_in_elt_0],resQPLloc_it->template get<0>() );
2361 formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
2362 isFirstExperience =
false;
2374 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2375 typename Integrator<Elements, Im, Expr, Im2>::eval::matrix_type
2376 Integrator<Elements, Im, Expr, Im2>::evaluate( mpl::int_<MESH_ELEMENTS> )
const
2378 DLOG(INFO) <<
"integrating over "
2379 << std::distance( this->beginElement(), this->endElement() ) <<
" elements\n";
2380 boost::timer __timer;
2382 #if defined(FEELPP_HAS_TBB)
2393 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
2394 typedef typename boost::remove_const<const_t>::type the_element_type;
2395 typedef the_element_type element_type;
2396 typedef typename the_element_type::gm_type gm_type;
2397 typedef boost::shared_ptr<gm_type> gm_ptrtype;
2398 typedef typename eval::gmc_type gmc_type;
2399 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2401 typedef typename the_element_type::gm1_type gm1_type;
2402 typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
2403 typedef typename eval::gmc1_type gmc1_type;
2404 typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
2408 typename eval::matrix_type res( eval::matrix_type::Zero() );
2410 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2412 auto it = lit->template get<1>();
2413 auto en = lit->template get<2>();
2431 gm_ptrtype gm( it->gm() );
2433 gm1_ptrtype gm1(
new gm1_type );
2436 typename eval::gmpc_ptrtype __geopc(
new typename eval::gmpc_type( gm,
2439 typename eval::gmpc1_ptrtype __geopc1(
new typename eval::gmpc1_type( gm1,
2446 #ifdef FEELPP_HAS_MPI
2447 auto const& worldComm =
const_cast<MeshBase*
>( it->mesh() )->worldComm();
2449 if ( worldComm.localSize() > 1 )
2451 worldComm.localComm().barrier();
2457 gmc_ptrtype __c(
new gmc_type( gm, *it, __geopc ) );
2458 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2459 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2461 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
2462 eval_expr_type expr( expression(), mapgmc );
2463 typedef typename eval_expr_type::shape shape;
2467 gmc1_ptrtype __c1(
new gmc1_type( gm1, *it, __geopc1 ) );
2468 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
2469 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
2471 typedef typename expression_type::template tensor<map_gmc1_type> eval_expr1_type;
2472 eval_expr1_type expr1( expression(), mapgmc1 );
2478 for ( ; it != en; ++it )
2483 case GeomapStrategyType::GEOMAP_HO :
2486 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2487 expr.update( mapgmc );
2488 const gmc_type& gmc = *__c;
2493 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2494 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2496 res( c1,c2 ) += M_im( expr, c1, c2 );
2501 case GeomapStrategyType::GEOMAP_O1:
2504 __c1->update( *it );
2505 map_gmc1_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
2506 expr1.update( mapgmc );
2507 const gmc1_type& gmc = *__c1;
2512 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2513 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2515 res( c1,c2 ) += M_im( expr1, c1, c2 );
2520 case GeomapStrategyType::GEOMAP_OPT:
2523 if ( it->isOnBoundary() )
2527 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2528 expr.update( mapgmc );
2529 const gmc_type& gmc = *__c;
2534 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2535 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2537 res( c1,c2 ) += M_im( expr, c1, c2 );
2544 __c1->update( *it );
2545 map_gmc1_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
2546 expr1.update( mapgmc );
2547 const gmc1_type& gmc = *__c1;
2551 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2552 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2554 res( c1,c2 ) += M_im( expr1, c1, c2 );
2563 DLOG(INFO) <<
"integrating over elements done in " << __timer.elapsed() <<
"s\n";
2569 #if defined(FEELPP_HAS_TBB)
2570 element_iterator it = this->beginElement();
2571 element_iterator en = this->endElement();
2572 typedef ContextEvaluate<expression_type, im_type, typename eval::the_element_type> context_type;
2575 return eval::zero();
2577 std::vector<boost::reference_wrapper<const typename eval::element_type> > _v;
2579 for (
auto _it = it; _it != en; ++_it )
2580 _v.push_back( boost::cref( *_it ) );
2582 tbb::blocked_range<decltype( _v.begin() )> r( _v.begin(), _v.end(), M_grainsize );
2583 context_type thecontext( this->expression(), this->im(), *it );
2585 if ( M_partitioner ==
"auto" )
2586 tbb::parallel_reduce( r, thecontext );
2588 else if ( M_partitioner ==
"simple" )
2589 tbb::parallel_reduce( r, thecontext, tbb::simple_partitioner() );
2593 return thecontext.result();
2594 #endif // FEELPP_HAS_TBB
2598 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2599 typename Integrator<Elements, Im, Expr, Im2>::eval::matrix_type
2600 Integrator<Elements, Im, Expr, Im2>::evaluate( mpl::int_<MESH_FACES> )
const
2602 DLOG(INFO) <<
"integrating over "
2603 << std::distance( this->beginElement(), this->endElement() ) <<
"faces\n";
2604 boost::timer __timer;
2609 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
2610 typedef typename boost::remove_const<const_t>::type the_face_element_type;
2611 typedef typename the_face_element_type::super2::entity_type the_element_type;
2612 typedef typename the_element_type::gm_type gm_type;
2613 typedef boost::shared_ptr<gm_type> gm_ptrtype;
2614 typedef typename gm_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, the_element_type> gmc_type;
2615 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2621 QuadMapped<im_type> qm;
2622 typename QuadMapped<im_type>::permutation_points_type ppts( qm( im() ) );
2623 typedef typename QuadMapped<im_type>::permutation_type permutation_type;
2629 typedef typename eval::gmpc_type pc_type;
2630 typedef typename eval::gmpc_ptrtype pc_ptrtype;
2631 std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( im().nFaces() );
2632 typedef typename im_type::face_quadrature_type face_im_type;
2634 std::vector<im_face_type> __integrators;
2636 typename eval::matrix_type res( eval::matrix_type::Zero() );
2637 typename eval::matrix_type res0( eval::matrix_type::Zero() );
2638 typename eval::matrix_type res1( eval::matrix_type::Zero() );
2641 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2643 auto it = lit->template get<1>();
2644 auto en = lit->template get<2>();
2648 if ( (it == en) || (it->isConnectedTo0()==
false) )
2652 CHECK( it->isConnectedTo0() ) <<
"invalid face with id=" << it->id();
2653 CHECK( it->element(0).gm() ) <<
"invalid geometric transformation assocated to face id="
2654 << it->id() <<
" and element id " << it->element(0).id();
2655 gm_ptrtype gm = it->element( 0 ).gm();
2658 for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2660 __integrators.push_back( im( __f ) );
2662 for ( permutation_type __p( permutation_type::IDENTITY );
2663 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2666 __geopc[__f][__p] = pc_ptrtype(
new pc_type( gm, ppts[__f].find( __p )->second ) );
2670 uint16_type __face_id_in_elt_0 = it->pos_first();
2673 gmc_ptrtype __c0(
new gmc_type( gm, it->element( 0 ), __geopc, __face_id_in_elt_0 ) );
2675 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2677 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
2678 typedef boost::shared_ptr<eval_expr_type> eval_expr_ptrtype;
2679 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
2680 eval_expr_ptrtype expr(
new eval_expr_type( expression(), mapgmc ) );
2683 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype>, fusion::pair<vf::detail::gmc<1>, gmc_ptrtype> > map2_gmc_type;
2684 typedef typename expression_type::template tensor<map2_gmc_type> eval2_expr_type;
2685 typedef boost::shared_ptr<eval2_expr_type> eval2_expr_ptrtype;
2686 eval2_expr_ptrtype expr2;
2689 bool isConnectedTo1 = it->isConnectedTo1();
2696 if ( isConnectedTo1 )
2698 uint16_type __face_id_in_elt_1 = it->pos_second();
2700 __c1 = gmc_ptrtype(
new gmc_type( gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
2702 map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
2703 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
2705 expr2 = eval2_expr_ptrtype(
new eval2_expr_type( expression(), mapgmc ) );
2706 expr2->init( im() );
2716 for ( ; it != en; ++it )
2718 if ( it->isGhostFace() )
2720 LOG(WARNING) <<
"face id : " << it->id() <<
" is a ghost face";
2723 if ( it->isConnectedTo1() )
2725 FEELPP_ASSERT( it->isOnBoundary() == false )
2726 ( it->id() ).error(
"face on boundary but connected on both sides" );
2727 uint16_type __face_id_in_elt_0 = it->pos_first();
2728 uint16_type __face_id_in_elt_1 = it->pos_second();
2730 __c0->update( it->element( 0 ), __face_id_in_elt_0 );
2731 __c1->update( it->element( 1 ), __face_id_in_elt_1 );
2734 std::cout <<
"face " << it->id() <<
"\n"
2735 <<
" id in elt = " << __face_id_in_elt_1 <<
"\n"
2736 <<
" elt 0 : " << it->element( 0 ).id() <<
"\n"
2737 <<
" elt 0 G: " << it->element( 0 ).G() <<
"\n"
2738 <<
" node elt 0 0 :" << it->element( 0 ).point( it->element( 0 ).fToP( __face_id_in_elt_0, 0 ) ).node() <<
"\n"
2739 <<
" node elt 0 1 :" << it->element( 0 ).point( it->element( 0 ).fToP( __face_id_in_elt_0, 1 ) ).node() <<
"\n"
2740 <<
" ref nodes 0 :" << __c0->xRefs() <<
"\n"
2741 <<
" real nodes 0: " << __c0->xReal() <<
"\n";
2742 std::cout <<
"face " << it->id() <<
"\n"
2743 <<
" id in elt = " << __face_id_in_elt_1 <<
"\n"
2744 <<
" elt 1 : " << it->element( 1 ).id() <<
"\n"
2745 <<
" elt 1 G: " << it->element( 1 ).G() <<
"\n"
2746 <<
" node elt 1 0 :" << it->element( 1 ).point( it->element( 1 ).fToP( __face_id_in_elt_1, 1 ) ).node() <<
"\n"
2747 <<
" node elt 1 1 :" << it->element( 1 ).point( it->element( 1 ).fToP( __face_id_in_elt_1, 0 ) ).node() <<
"\n"
2748 <<
" ref nodes 1 :" << __c1->xRefs() <<
"\n"
2749 <<
" real nodes 1:" << __c1->xReal() <<
"\n";
2752 __typeof__( im( __face_id_in_elt_0 ) ) im_face ( im( __face_id_in_elt_0 ) );
2754 map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
2755 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
2757 expr2->update( mapgmc, __face_id_in_elt_0 );
2758 const gmc_type& gmc = *__c0;
2760 __integrators[__face_id_in_elt_0].update( gmc );
2762 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2763 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2765 res( c1,c2 ) += __integrators[__face_id_in_elt_0]( *expr2, c1, c2 );
2774 uint16_type __face_id_in_elt_0 = it->pos_first();
2775 __c0->update( it->element( 0 ), __face_id_in_elt_0 );
2776 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
2777 expr->update( mapgmc, __face_id_in_elt_0 );
2779 const gmc_type& gmc = *__c0;
2781 __integrators[__face_id_in_elt_0].update( gmc );
2783 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2784 for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2786 res( c1,c2 ) += __integrators[__face_id_in_elt_0]( *expr, c1, c2 );
2793 DLOG(INFO) <<
"integrating over faces done in " << __timer.elapsed() <<
"s\n";
2797 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2798 typename Integrator<Elements, Im, Expr, Im2>::eval::matrix_type
2799 Integrator<Elements, Im, Expr, Im2>::evaluate( mpl::int_<MESH_POINTS> )
const
2801 DLOG(INFO) <<
"integrating over "
2802 << std::distance( this->beginElement(), this->endElement() ) <<
" points\n";
2810 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2811 template<
typename P0hType>
2812 typename P0hType::element_type
2813 Integrator<Elements, Im, Expr, Im2>::broken( boost::shared_ptr<P0hType>& P0h, mpl::int_<MESH_ELEMENTS> )
const
2815 DLOG(INFO) <<
"integrating over "
2816 << std::distance( this->beginElement(), this->endElement() ) <<
" elements\n";
2817 boost::timer __timer;
2822 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
2823 typedef typename boost::remove_const<const_t>::type the_element_type;
2824 typedef the_element_type element_type;
2825 typedef typename the_element_type::gm_type gm_type;
2826 typedef boost::shared_ptr<gm_type> gm_ptrtype;
2827 typedef typename gm_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::POINT, the_element_type> gmc_type;
2828 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2835 auto p0 = P0h->element(
"p0" );
2839 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2841 auto it = lit->template get<1>();
2842 auto en = lit->template get<2>();
2855 gm_ptrtype gm = it->gm();
2857 typename eval::gmpc_ptrtype __geopc(
new typename eval::gmpc_type( gm,
2861 it = this->beginElement();
2863 #ifdef FEELPP_HAS_MPI
2864 auto const& worldComm =
const_cast<MeshBase*
>( it->mesh() )->worldComm();
2867 if ( worldComm.size() > 1 )
2869 worldComm.barrier();
2876 gmc_ptrtype __c(
new gmc_type( gm, *it, __geopc ) );
2877 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2878 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2880 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
2881 eval_expr_type expr( expression(), mapgmc );
2882 typedef typename eval_expr_type::shape shape;
2885 for ( ; it != en; ++it )
2889 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2890 expr.update( mapgmc );
2891 const gmc_type& gmc = *__c;
2896 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2898 size_type i= P0h->dof()->localToGlobal( it->id(), 0, c1 ).index();
2899 double v = M_im( expr, c1, 0 );
2906 DLOG(INFO) <<
"integrating over elements done in " << __timer.elapsed() <<
"s\n";
2910 template<
typename Elements,
typename Im,
typename Expr,
typename Im2>
2911 template<
typename P0hType>
2912 typename P0hType::element_type
2913 Integrator<Elements, Im, Expr, Im2>::broken( boost::shared_ptr<P0hType>& P0h, mpl::int_<MESH_FACES> )
const
2915 DLOG(INFO) <<
"integrating over "
2916 << std::distance( this->beginElement(), this->endElement() ) <<
"faces\n";
2917 boost::timer __timer;
2922 typedef typename boost::remove_reference<typename element_iterator::reference>::type const_t;
2923 typedef typename boost::remove_const<const_t>::type the_face_element_type;
2924 typedef typename the_face_element_type::super2::entity_type the_element_type;
2925 typedef typename the_element_type::gm_type gm_type;
2926 typedef boost::shared_ptr<gm_type> gm_ptrtype;
2927 typedef typename gm_type::template Context<expression_type::context|vm::JACOBIAN|vm::KB|vm::NORMAL|vm::POINT, the_element_type> gmc_type;
2928 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2934 QuadMapped<im_type> qm;
2935 typename QuadMapped<im_type>::permutation_points_type ppts( qm( im() ) );
2936 typedef typename QuadMapped<im_type>::permutation_type permutation_type;
2942 typedef typename eval::gmpc_type pc_type;
2943 typedef typename eval::gmpc_ptrtype pc_ptrtype;
2944 std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( im().nFaces() );
2945 typedef typename im_type::face_quadrature_type face_im_type;
2947 std::vector<im_face_type> __integrators;
2949 auto p0 = P0h->element(
"p0" );
2953 for(
auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2955 auto it = lit->template get<1>();
2956 auto en = lit->template get<2>();
2964 gm_ptrtype gm = it->element( 0 ).gm();
2967 for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2969 __integrators.push_back( im( __f ) );
2971 for ( permutation_type __p( permutation_type::IDENTITY );
2972 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2975 __geopc[__f][__p] = pc_ptrtype(
new pc_type( gm, ppts[__f].find( __p )->second ) );
2980 uint16_type __face_id_in_elt_0 = it->pos_first();
2983 gmc_ptrtype __c0(
new gmc_type( gm, it->element( 0 ), __geopc, __face_id_in_elt_0 ) );
2985 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2987 typedef typename expression_type::template tensor<map_gmc_type> eval_expr_type;
2988 typedef boost::shared_ptr<eval_expr_type> eval_expr_ptrtype;
2989 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
2990 eval_expr_ptrtype expr(
new eval_expr_type( expression(), mapgmc ) );
2993 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype>, fusion::pair<vf::detail::gmc<1>, gmc_ptrtype> > map2_gmc_type;
2994 typedef typename expression_type::template tensor<map2_gmc_type> eval2_expr_type;
2995 typedef boost::shared_ptr<eval2_expr_type> eval2_expr_ptrtype;
2996 eval2_expr_ptrtype expr2;
2999 bool isConnectedTo1 = it->isConnectedTo1();
3006 if ( isConnectedTo1 )
3008 uint16_type __face_id_in_elt_1 = it->pos_second();
3010 __c1 = gmc_ptrtype(
new gmc_type( gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
3012 map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
3013 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
3015 expr2 = eval2_expr_ptrtype(
new eval2_expr_type( expression(), mapgmc ) );
3016 expr2->init( im() );
3026 for ( ; it != en; ++it )
3029 if ( it->isConnectedTo1() )
3031 FEELPP_ASSERT( it->isOnBoundary() == false )
3032 ( it->id() ).error(
"face on boundary but connected on both sides" );
3033 uint16_type __face_id_in_elt_0 = it->pos_first();
3034 uint16_type __face_id_in_elt_1 = it->pos_second();
3036 __c0->update( it->element( 0 ), __face_id_in_elt_0 );
3037 __c1->update( it->element( 1 ), __face_id_in_elt_1 );
3040 std::cout <<
"face " << it->id() <<
"\n"
3041 <<
" id in elt = " << __face_id_in_elt_1 <<
"\n"
3042 <<
" elt 0 : " << it->element( 0 ).id() <<
"\n"
3043 <<
" elt 0 G: " << it->element( 0 ).G() <<
"\n"
3044 <<
" node elt 0 0 :" << it->element( 0 ).point( it->element( 0 ).fToP( __face_id_in_elt_0, 0 ) ).node() <<
"\n"
3045 <<
" node elt 0 1 :" << it->element( 0 ).point( it->element( 0 ).fToP( __face_id_in_elt_0, 1 ) ).node() <<
"\n"
3046 <<
" ref nodes 0 :" << __c0->xRefs() <<
"\n"
3047 <<
" real nodes 0: " << __c0->xReal() <<
"\n";
3048 std::cout <<
"face " << it->id() <<
"\n"
3049 <<
" id in elt = " << __face_id_in_elt_1 <<
"\n"
3050 <<
" elt 1 : " << it->element( 1 ).id() <<
"\n"
3051 <<
" elt 1 G: " << it->element( 1 ).G() <<
"\n"
3052 <<
" node elt 1 0 :" << it->element( 1 ).point( it->element( 1 ).fToP( __face_id_in_elt_1, 1 ) ).node() <<
"\n"
3053 <<
" node elt 1 1 :" << it->element( 1 ).point( it->element( 1 ).fToP( __face_id_in_elt_1, 0 ) ).node() <<
"\n"
3054 <<
" ref nodes 1 :" << __c1->xRefs() <<
"\n"
3055 <<
" real nodes 1:" << __c1->xReal() <<
"\n";
3058 __typeof__( im( __face_id_in_elt_0 ) ) im_face ( im( __face_id_in_elt_0 ) );
3060 map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
3061 fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
3063 expr2->update( mapgmc, __face_id_in_elt_0 );
3064 const gmc_type& gmc = *__c0;
3066 __integrators[__face_id_in_elt_0].update( gmc );
3068 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
3070 size_type i0 = P0h->dof()->localToGlobal( it->element( 0 ), 0, c1 ).index();
3071 size_type i1 = P0h->dof()->localToGlobal( it->element( 1 ), 0, c1 ).index();
3072 double v = __integrators[__face_id_in_elt_0]( *expr2, c1, 0 );
3080 uint16_type __face_id_in_elt_0 = it->pos_first();
3081 __c0->update( it->element( 0 ), __face_id_in_elt_0 );
3082 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ) );
3083 expr->update( mapgmc, __face_id_in_elt_0 );
3085 const gmc_type& gmc = *__c0;
3087 __integrators[__face_id_in_elt_0].update( gmc );
3089 for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
3091 size_type i0 = P0h->dof()->localToGlobal( it->element( 0 ), 0, c1 ).index();
3092 double v = __integrators[__face_id_in_elt_0]( *expr, c1, 0 );
3100 DLOG(INFO) <<
"integrating over faces done in " << __timer.elapsed() <<
"s\n";
3110 template<
typename Elts,
typename Im,
typename ExprT>
3111 Expr<Integrator<Elts, Im, ExprT, Im> >
3115 GeomapStrategyType gt = GeomapStrategyType::GEOMAP_HO )
3117 typedef Integrator<Elts, Im, ExprT, Im> expr_t;
3118 return Expr<expr_t>( expr_t( elts, im, expr, gt, im ) );
3123 # define VF_VALUE_OF_IM(O) \
3124 boost::mpl::if_< boost::mpl::bool_< O::imIsPoly > , \
3125 typename boost::mpl::if_< boost::mpl::greater< boost::mpl::int_<O::imorder>, boost::mpl::int_<19> > , \
3126 boost::mpl::int_<19>, \
3127 boost::mpl::int_<O::imorder> >::type >::type , \
3128 boost::mpl::int_<10> >::type::value \
3136 template<
typename Elts,
typename ExprT>
3137 Expr<Integrator<Elts, _Q< ExpressionOrder<Elts,ExprT>::value >, ExprT, _Q< ExpressionOrder<Elts,ExprT>::value > > >
3140 GeomapStrategyType gt = GeomapStrategyType::GEOMAP_HO )
3143 DLOG(INFO) <<
"[integrate] order to integrate = " << ExpressionOrder<Elts,ExprT>::value <<
"\n";
3144 _Q< ExpressionOrder<Elts,ExprT>::value > quad;
3145 return integrate_impl( elts, quad, expr, gt, quad );
3154 template<
typename Elts,
typename Im,
typename ExprT,
typename Im2 = Im>
3155 Expr<Integrator<typename Feel::detail::quadptlocrangetype<Elts>::type, Im, ExprT, Im2> >
3156 integrate_impl( Elts
const& elts,
3159 GeomapStrategyType
const& gt,
3163 std::string
const& partitioner,
3164 boost::shared_ptr<QuadPtLocalization<
typename Feel::detail::quadptlocrangetype<Elts>::type, Im, ExprT > > quadptloc
3165 = boost::shared_ptr<QuadPtLocalization<
typename Feel::detail::quadptlocrangetype<Elts>::type, Im, ExprT > >() )
3167 typedef Integrator<typename Feel::detail::quadptlocrangetype<Elts>::type, Im, ExprT, Im2> expr_t;
3168 return Expr<expr_t>( expr_t( elts, im, expr, gt, im2, use_tbb, grainsize, partitioner, quadptloc ) );
3175 template<
typename TheArgs,
typename Tag>
3178 typedef typename boost::remove_pointer<
3179 typename boost::remove_const<
3180 typename boost::remove_reference<
3181 typename parameter::binding<TheArgs, Tag>::type
3186 template<
typename TheArgs,
typename Tag,
typename Default>
3189 typedef typename boost::remove_pointer<
3190 typename boost::remove_const<
3191 typename boost::remove_reference<
3192 typename parameter::binding<TheArgs, Tag, Default>::type
3198 template<
typename Args>
3199 struct integrate_type
3201 typedef typename clean_type<Args,tag::expr>::type _expr_type;
3202 typedef typename Feel::detail::quadptlocrangetype<typename clean_type<Args,tag::range>::type>::type _range_type;
3204 typedef typename clean2_type<Args,tag::quad, _Q< ExpressionOrder<_range_type,_expr_type>::value > >::type _quad_type;
3205 typedef typename clean2_type<Args,tag::quad1, _Q< ExpressionOrder<_range_type,_expr_type>::value_1 > >::type _quad1_type;
3206 typedef Expr<Integrator<_range_type, _quad_type, _expr_type, _quad1_type> > expr_type;
3208 typedef boost::shared_ptr<QuadPtLocalization<_range_type,_quad_type,_expr_type > > _quadptloc_ptrtype;
3223 BOOST_PARAMETER_FUNCTION(
3224 (
typename vf::detail::integrate_type<Args>::expr_type ),
3235 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3236 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3237 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3238 ( use_tbb, (
bool ),
false )
3239 ( grainsize, (
int ), 100 )
3240 ( partitioner, *,
"auto" )
3241 ( verbose, (
bool ),
false )
3242 ( quadptloc, *,
typename vf::detail::integrate_type<Args>::_quadptloc_ptrtype() )
3247 auto ret = integrate_impl( range, quad, expr, geomap, quad1, use_tbb, grainsize, partitioner, quadptloc );
3251 std::cout <<
" -- integrate: size(range) = " << std::distance( ret.expression().beginElement(),
3252 ret.expression().endElement() ) <<
"\n";
3253 std::cout <<
" -- integrate: quad = " << ret.expression().im().nPoints() <<
"\n";
3254 std::cout <<
" -- integrate: quad1 = " << ret.expression().im2().nPoints() <<
"\n";
3261 BOOST_PARAMETER_FUNCTION(
3273 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3274 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3275 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3276 ( use_tbb, (
bool ),
false )
3277 ( grainsize, (
int ), 100 )
3278 ( partitioner, *,
"auto" )
3279 ( verbose, (
bool ),
false )
3283 double a =
integrate( _range=range, _expr=inner(expr,expr), _quad=quad, _geomap=geomap,
3284 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3285 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3286 return math::sqrt( a );
3289 BOOST_PARAMETER_FUNCTION(
3301 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3302 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3303 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3304 ( use_tbb, (
bool ),
false )
3305 ( grainsize, (
int ), 100 )
3306 ( partitioner, *,
"auto" )
3307 ( verbose, (
bool ),
false )
3311 return integrate( _range=range, _expr=inner(expr,expr), _quad=quad, _geomap=geomap,
3312 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3313 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3316 BOOST_PARAMETER_FUNCTION(
3329 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3330 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3331 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3332 ( use_tbb, (
bool ),
false )
3333 ( grainsize, (
int ), 100 )
3334 ( partitioner, *,
"auto" )
3335 ( verbose, (
bool ),
false )
3339 double a =
integrate( _range=range, _expr=inner(expr,expr), _quad=quad, _geomap=geomap,
3340 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3341 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3342 double b =
integrate( _range=range, _expr=grad_expr*trans(grad_expr), _quad=quad, _geomap=geomap,
3343 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3344 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3345 return math::sqrt( a + b );
3348 BOOST_PARAMETER_FUNCTION(
3360 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3361 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3362 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3363 ( use_tbb, (
bool ),
false )
3364 ( grainsize, (
int ), 100 )
3365 ( partitioner, *,
"auto" )
3366 ( verbose, (
bool ),
false )
3370 double a =
integrate( _range=range, _expr=grad_expr*trans(grad_expr), _quad=quad, _geomap=geomap,
3371 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3372 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3373 return math::sqrt( a );
3376 BOOST_PARAMETER_FUNCTION(
3377 (
typename vf::detail::integrate_type<Args>::expr_type::expression_type::matrix_type ),
3388 ( quad, *,
typename vf::detail::integrate_type<Args>::_quad_type() )
3389 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
3390 ( quad1, *,
typename vf::detail::integrate_type<Args>::_quad1_type() )
3391 ( use_tbb, (
bool ),
false )
3392 ( grainsize, (
int ), 100 )
3393 ( partitioner, *,
"auto" )
3394 ( verbose, (
bool ),
false )
3398 double meas =
integrate( _range=range, _expr=cst(1.0), _quad=quad, _quad1=quad1, _geomap=geomap,
3399 _use_tbb=use_tbb, _grainsize=grainsize,
3400 _partitioner=partitioner, _verbose=verbose ).evaluate()( 0, 0 );
3401 DLOG(INFO) <<
"[mean] measure = " << meas <<
"\n";
3402 CHECK( math::abs(meas) > 1e-13 ) <<
"Invalid domain measure : " << meas <<
", domain range: " <<
nelements( range ) <<
"\n";
3403 auto eint =
integrate( _range=range, _expr=expr, _quad=quad, _geomap=geomap,
3404 _quad1=quad1, _use_tbb=use_tbb, _grainsize=grainsize,
3405 _partitioner=partitioner, _verbose=verbose ).evaluate();
3406 DLOG(INFO) <<
"[ein] eint = " << eint <<
"\n";
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
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
size_type nelements(boost::tuple< MT, Iterator, Iterator > const &its, bool global=false)
Definition: filters.hpp:1375
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319