Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
integrator.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-01-20
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2010 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __Integrators_H
31 #define __Integrators_H 1
32 #include <cxxabi.h>
33 #include <typeinfo>
34 #include <boost/timer.hpp>
35 #include <feel/feelcore/feel.hpp>
36 #include <feel/feelcore/parameter.hpp>
37 #include <feel/feelvf/block.hpp>
38 
39 
40 #include <Eigen/Eigen>
41 
42 
43 //#include <boost/numeric/bindings/traits/traits.hpp>
44 //#include <boost/numeric/bindings/traits/sparse_traits.hpp>
45 //#include <boost/numeric/bindings/traits/ublas_sparse.hpp>
46 
53 #if defined( FEELPP_HAS_GOOGLE_PROFILER_H )
54 #include <google/profiler.h>
55 #endif
56 
57 namespace Feel
58 {
59 namespace vf
60 {
61 
63 enum IntegratorType
64 {
65  INT_0D = -1,
66  INT_1D = 1,
67  INT_2D = 2,
68  INT_3D = 3,
69  INT_ELEMENTS = 4,
70  INT_ELEMENT_FACES = 5,
71  INT_INTERNAL_FACES = 6
73 };
74 
82 template<typename Elements, typename Im, typename Expr, typename Im2=Im>
83 class Integrator: public IntegratorBase
84 {
85 public:
86 
90 
91  static const size_type context = Expr::context;
92  static const bool is_terminal = false;
93 
94  static const uint16_type imorder = 0;
95  static const bool imIsPoly = true;
96 
97  template<typename Func>
98  struct HasTestFunction
99  {
100  static const bool result = Expr::template HasTestFunction<Func>::result;
101  };
102 
103  template<typename Func>
104  struct HasTrialFunction
105  {
106  static const bool result = Expr::template HasTrialFunction<Func>::result;
107  };
108 
109  static const size_type iDim = boost::tuples::template element<0, Elements>::type::value;
111 
115 
116  typedef Integrator<Elements, Im, Expr, Im2> self_type;
117 
118  typedef typename boost::tuples::template element<1, Elements>::type element_iterator;
119 
120  typedef Expr expression_type;
121  typedef typename expression_type::value_type expression_value_type;
122  typedef ublas::vector<int> vector_dof_type;
123 
124  struct eval
125  {
126  //
127  // some typedefs
128  //
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;
132 
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;
137 
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;
142 
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;
148  //typedef typename gm_type::template Context<expression_type::context, the_element_type, im_type::numPoints> gmc_type;
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;
153 #if 0
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;
156 #else
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;
161 #endif
162  //typedef typename eval_expr_type::value_type value_type;
163  //typedef typename strongest_numeric_type<typename Im::value_type, typename expression_type::value_type>::type value_type;
164  //typedef typename expression_type::value_type value_type;
165 
166  //
167  // Precompute some data in the reference element for
168  // geometric mapping and reference finite element
169  //
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;
173 
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; // should be the same as shape
177  //typedef typename shape_type::storage<value_type> storage_type;
178  /*
179  typedef mpl::if_<mpl::bool_<shape_type::is_scalar>,
180  mpl::identity<value_type>,
181  mpl::identity<ublas::vector<value_type,storage_type> > >::type::type ret_type;*/
182  //typedef ublas::matrix<value_type> ret_type;
183 #if 0
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;
188 #else
189  typedef Eigen::Matrix<expression_value_type, shape::M, shape::N> value_type;
190 #endif
191  typedef Eigen::Matrix<expression_value_type, shape::M, shape::N> matrix_type;
192  static value_type zero( mpl::bool_<false> )
193  {
194  return value_type::Zero();
195  }
196  static value_type zero( mpl::bool_<true> )
197  {
198  return 0;
199  }
200  static value_type zero()
201  {
202  return zero( boost::is_scalar<value_type>() );
203  }
204 
205  };
206 
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;
211  //typedef typename eval::value_type value_type;
212  typedef typename eval::matrix_type matrix_type;
213  typedef typename eval::matrix_type value_type;
215 
219 
220  Integrator( Elements const& elts, Im const& /*__im*/, expression_type const& __expr, GeomapStrategyType gt, Im2 const& /*__im2*/, bool use_tbb, int grainsize, std::string const& partitioner,
221  boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > qpl )
222  :
223  M_elts(),
224  M_eltbegin( elts.template get<1>() ),
225  M_eltend( elts.template get<2>() ),
226  M_im( ),
227  M_im2( ),
228  M_expr( __expr ),
229  M_gt( gt ),
230  M_use_tbb( use_tbb ),
231  M_grainsize( grainsize ),
232  M_partitioner( partitioner ),
233  M_QPL( qpl )
234  {
235  M_elts.push_back( elts );
236  DLOG(INFO) << "Integrator constructor from expression\n";
237  }
238 
239  Integrator( std::list<Elements> const& elts, Im const& /*__im*/, expression_type const& __expr,
240  GeomapStrategyType gt, Im2 const& /*__im2*/, bool use_tbb, int grainsize, std::string const& partitioner,
241  boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > qpl )
242  :
243  M_elts( elts ),
244  M_im( ),
245  M_im2( ),
246  M_expr( __expr ),
247  M_gt( gt ),
248  M_use_tbb( use_tbb ),
249  M_grainsize( grainsize ),
250  M_partitioner( partitioner ),
251  M_QPL( qpl )
252  {
253  DLOG(INFO) << "Integrator constructor from expression\n";
254  if ( elts.size() )
255  {
256  M_eltbegin = elts.begin()->template get<1>();
257  M_eltend = elts.begin()->template get<2>();
258  }
259  }
260 
261  Integrator( Integrator const& __vfi )
262  :
263  M_elts( __vfi.M_elts) ,
264  M_eltbegin( __vfi.M_eltbegin ),
265  M_eltend( __vfi.M_eltend ),
266  M_im( __vfi.M_im ),
267  M_im2( __vfi.M_im2 ),
268  M_expr( __vfi.M_expr ),
269  M_gt( __vfi.M_gt ),
270  M_use_tbb( __vfi.M_use_tbb ),
271  M_grainsize( __vfi.M_grainsize ),
272  M_partitioner( __vfi.M_partitioner ),
273  M_QPL( __vfi.M_QPL )
274  {
275  DLOG(INFO) << "Integrator copy constructor\n";
276  }
277 
278  virtual ~Integrator() {}
279 
281 
285 
286  template<typename TheExpr>
287  struct Lambda
288  {
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;
293  };
294 
295  template<typename ExprT>
296  typename Lambda<ExprT>::type
297  operator()( ExprT const& e )
298  {
299 #if 0
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 ) );
304 #endif
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;
311  quad_type quad;
312  quad1_type quad1;
313  //BOOST_STATIC_ASSERT( ( boost::is_same<expr_type,e_type> ) );
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 );
321  return i;
322  }
323 
325 
329 
336  im_type const& im() const
337  {
338  return M_im;
339  }
340 
347  im2_type const& im2() const
348  {
349  return M_im2;
350  }
351 
358  im_face_type im( uint16_type f ) const
359  {
360  return M_im.face( f );
361  }
362 
369  im2_face_type im2( uint16_type f ) const
370  {
371  return M_im2.face( f );
372  }
373 
380  expression_type const& expression() const
381  {
382  return M_expr;
383  }
384 
385 
391  GeomapStrategyType geomapIntegratorType() const
392  {
393  return M_gt;
394  }
395 
400  element_iterator beginElement() const
401  {
402  return M_eltbegin;
403  }
404 
409  element_iterator endElement() const
410  {
411  return M_eltend;
412  }
413 
420  bool isSymmetric() const
421  {
422  return M_expr.isSymmetric();
423  }
424 
426 
430  void setBeginElement( element_iterator it ) { M_eltbegin = it; }
431  void setEndElement( element_iterator en ) { M_eltend = en; }
432 
434 
438 
439 
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;
444 
445  template<typename Elem1, typename FormType>
446  void assemble( boost::shared_ptr<Elem1> const& __v,
447  FormType& __form ) const;
448 
449  //typename expression_type::template tensor<Geo_t>::value_type
450  template<typename P0hType>
451  typename P0hType::element_type broken( boost::shared_ptr<P0hType>& P0h ) const
452  {
453  auto p0 = broken( P0h, mpl::int_<iDim>() );
454  return p0;
455  }
456 #if 1
457  matrix_type
458  evaluate( bool parallel=true,
459  WorldComm const& worldcomm = Environment::worldComm() ) const
460 #else
461  //typename expression_type::template tensor<Geo_t>::value_type
462  BOOST_PARAMETER_MEMBER_FUNCTION( ( matrix_type ),
463  evaluate,
464  tag,
465  /*(required
466  //(h,*(double)))*/
467  ( optional
468  ( parallel,*( bool ), true ) ) )
469 #endif
470  {
471  typename eval::matrix_type loc = evaluate( mpl::int_<iDim>() );
472 
473  if ( !parallel )
474  return loc;
475 
476  else // parallel
477  {
478  typename eval::matrix_type glo( loc );
479  // maybe better to create anoter worldcomm which split the mesh worldcomm
480  // with only partition that contains at least one element (Vincent C.)
481  // and thus argument worldComm can be remove
482  // auto const& worldcomm = const_cast<MeshBase*>( this->beginElement()->mesh() )->worldComm();
483 
484  if ( worldcomm.localSize() > 1 )
485  {
486  mpi::all_reduce( worldcomm.localComm(),
487  loc,
488  glo,
489  [] ( matrix_type const& x, matrix_type const& y )
490  {
491  return x + y;
492  } );
493  }
494 
495  return glo;
496  }
497  }
498 
499 #if defined( FEELPP_HAS_TBB )
500  template<typename FormType, typename ExprType, typename IMType, typename EltType>
501  class Context
502  {
503  public:
504  //
505  // some typedefs
506  //
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;
514  //typedef vf::detail::FormContextBase<map_gmc_type,im_type> fcb_type;
515  typedef form_context_type fcb_type;
516  typedef boost::shared_ptr<fcb_type> fcb_ptrtype;
517 
518 
519  Context( FormType& _form,
520  ExprType const& _expr,
521  IMType const& _im,
522  EltType const& _elt )
523  :
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 ),
529  _expr,
530  _im ) )
531 
532  {
533  }
534  Context( Context const& c )
535  :
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 ) )
539  {
540  }
541  //std::vector<boost::reference_wrapper<const typename mesh_type::element_type> > _v;
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
544  {
545 #if 1
546  tbb::mutex m;
547  tbb::mutex::scoped_lock lock( m );
548  lock.release();
549  auto mapgmc = fusion::make_pair<vf::detail::gmc<0> >( M_c );
550 
551  for ( auto _elt = r.begin(); _elt != r.end(); ++_elt )
552  {
553  M_c->update( *_elt );
554  M_formc->update( mapgmc, mapgmc );
555  M_formc->integrate();
556  lock.acquire( m );
557  M_formc->assemble();
558  lock.release();
559  }
560 
561 #else
562 
563  for ( auto _elt = r.begin(); _elt != r.end(); ++_elt )
564  {
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();
569  }
570 
571 #endif
572  }
573  mutable gmpc_ptrtype M_geopc;
574  mutable gmc_ptrtype M_c;
575  mutable fcb_ptrtype M_formc;
576 
577  };
578 
579  template<typename ExprType, typename IMType, typename EltType>
580  class ContextEvaluate
581  {
582  public:
583  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
584  //
585  // some typedefs
586  //
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,
600  IMType const& _im,
601  EltType const& _elt )
602  :
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 ) ) ),
607  M_im( _im ),
608  M_ret( eval::matrix_type::Zero() )
609  {
610  }
611  ContextEvaluate( ContextEvaluate& c, tbb::split )
612  :
613  M_gm( new gm_type( *c.M_gm ) ),
614  //M_geopc( new gmpc_type( M_gm, c.M_im.points() ) ),
615  M_geopc( c.M_geopc ),
616  M_c( new gmc_type( M_gm, c.M_c->element(), M_geopc ) ),
617  M_expr( c.M_expr ),
618  M_im( c.M_im ),
619  M_ret( eval::matrix_type::Zero() )
620  {}
621 
622  ContextEvaluate( ContextEvaluate const& c )
623  :
624  M_gm( new gm_type( *c.M_gm ) ),
625  //M_geopc( new gmpc_type( M_gm, c.M_im.points() ) ),
626  M_geopc( c.M_geopc ),
627  M_c( new gmc_type( M_gm, c.M_c->element(), M_geopc ) ),
628  M_expr( c.M_expr ),
629  M_im( c.M_im ),
630  M_ret( c.M_ret )
631  {
632  }
633  //std::vector<boost::reference_wrapper<const typename mesh_type::element_type> > _v;
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 )
636  {
637 #if 1
638 
639  for ( auto _elt = r.begin(); _elt != r.end(); ++_elt )
640  {
641  M_c->update( *_elt );
642  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( M_c ) );
643 
644  M_expr.update( mapgmc );
645  M_im.update( *M_c );
646 
647 #if 1
648 
649  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
650  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
651  {
652  M_ret( c1,c2 ) += M_im( M_expr, c1, c2 );
653  }
654 
655 #endif
656  }
657 
658 #else
659 #if 0
660 
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 )
664  {
665  M_ret( c1,c2 ) += i*i; //M_im( M_expr, c1, c2 );
666  }
667 
668 #endif
669 #endif
670  }
671  void join( ContextEvaluate const& other )
672  {
673  M_ret += other.M_ret;
674  }
675 
676  value_type result() const
677  {
678  return M_ret;
679  }
680 
681  gm_ptrtype M_gm;
682  gmpc_ptrtype M_geopc;
683  gmc_ptrtype M_c;
684  eval_expr_type M_expr;
685  im_type M_im;
686  value_type M_ret;
687  };
688 #endif // FEELPP_HAS_TBB
689 
690 
691 private:
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;
696 
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;
701 
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;
710 
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;
719 
720 
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;
725 
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;
729 
730 private:
731 
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;
739  bool M_use_tbb;
740  int M_grainsize;
741  std::string M_partitioner;
742  mutable boost::shared_ptr<QuadPtLocalization<Elements, Im, Expr > > M_QPL;
743  // mutable boost::prof::basic_profiler<boost::prof::basic_profile_manager<std::string, double, boost::high_resolution_timer, boost::prof::empty_logging_policy, boost::prof::default_stats_policy<std::string, double> > > M_profile_local_assembly;
744 
745  // mutable boost::prof::basic_profiler<boost::prof::basic_profile_manager<std::string, double, boost::high_resolution_timer, boost::prof::empty_logging_policy, boost::prof::default_stats_policy<std::string, double> > > M_profile_global_assembly;
746 };
747 
748 template<typename Elements, typename Im, typename Expr, typename Im2>
749 template<typename Elem1, typename Elem2, typename FormType>
750 void
751 Integrator<Elements, Im, Expr, Im2>::assemble( boost::shared_ptr<Elem1> const& __u,
752  boost::shared_ptr<Elem2> const& __v,
753  FormType& __form ) const
754 {
755 #if 0
756  details::GlobalMatrixAssembler<iDim, self_type, Elem1, Elem2, FormType>( *this,
757  __u,
758  __v,
759  __form );
760 #endif
761 
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;
765 
766  element_iterator it, en;
767  bool findEltForInit = false;
768  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
769  {
770  it = lit->template get<1>();
771  en = lit->template get<2>();
772 
773  // check that we have elements to iterate over
774  if ( it == en )
775  continue;
776  else
777  findEltForInit = true;
778  }
779  if (!findEltForInit) return;
780 
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 )
785  {
786  DLOG(INFO) << "[integrator::assemble bilinear form] with_relation_mesh (same_mesh: " << same_mesh_type::value/*same_mesh*/ << " 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 );
788  }
789  else
790  {
791  DLOG(INFO) << "[integrator::assemble bilinear form] no_relation_mesh\n";
792  assemble( __form, mpl::int_<iDim>(), mpl::bool_<false>(), test_related_to_trial );
793  }
794 
795 }
796 
797 
798 template<typename Elements, typename Im, typename Expr, typename Im2>
799 template<typename Elem1, typename FormType>
800 void
801 Integrator<Elements, Im, Expr, Im2>::assemble( boost::shared_ptr<Elem1> const& __v,
802  FormType& __form ) const
803 {
804 #if 0
805  details::GlobalVectorAssembler<iDim, self_type, Elem1, FormType>( *this,
806  __v,
807  __form );
808 #endif
809 
810  typedef typename boost::is_same<typename eval::gmc_type::element_type,typename Elem1::mesh_type::element_type>::type same_mesh_type;
811 
812  element_iterator it, en;
813  bool findEltForInit = false;
814  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len && !findEltForInit; ++lit )
815  {
816  it = lit->template get<1>();
817  en = lit->template get<2>();
818 
819  // check that we have elements to iterate over
820  if ( it == en )
821  continue;
822  else
823  findEltForInit = true;
824  }
825  if (!findEltForInit) return;
826 
827 
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 );
830 
831  else
832  assemble( __form, mpl::int_<iDim>(), mpl::bool_<false>(),false );
833 
834  //assemble( __form, mpl::int_<iDim>(), mpl::bool_<true>() );
835 }
836 template<typename Elements, typename Im, typename Expr, typename Im2>
837 template<typename FormType>
838 void
839 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<true> , bool /*hasRelation*/ ) const
840 {
841  boost::timer __timer;
842  LOG(INFO) << "[integrator::assemble FormType& __form, mpl::int_<MESH_ELEMENTS> /**/, mpl::bool_<true>\n";
843 
844 #if defined(FEELPP_HAS_TBB)
845 
846  //std::cout << "Integrator Uses TBB: " << M_use_tbb << "\n";
847  if ( !M_use_tbb )
848 #else
849  if ( 1 )
850 #endif
851  {
852  //
853  // some typedefs
854  //
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;
861 
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;
866  //typedef vf::detail::FormContextBase<map_gmc_type,im_type> fcb_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;
871 
872  //
873  // Precompute some data in the reference element for
874  // geometric mapping and reference finite element
875  //
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() ) );
878 
879 
880  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
881  {
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";
885 
886  // check that we have elements to iterate over
887  if ( it == en )
888  continue;
889 
890  size_type idEltTestInit = it->id();
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 );
895 
896  auto const& eltTestInit = __form.testSpace()->mesh()->element( idEltTestInit );
897 
898 
899  gmc_ptrtype __c( new gmc_type( __form.gm(), eltTestInit, __geopc ) );
900  gmc1_ptrtype __c1( new gmc1_type( __form.gm1(), eltTestInit, __geopc1 ) );
901 
902 
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 ) );
905 
906  focb_ptrtype formc( new form_context_type( __form,
907  mapgmc,
908  mapgmc,
909  mapgmc,
910  this->expression(),
911  this->im() ) );
912  focb1_ptrtype formc1( new form1_context_type( __form,
913  mapgmc1,
914  mapgmc1,
915  mapgmc1,
916  this->expression(),
917  this->im2() ) );
918 
919  //int nelt = std::distance( this->beginElement(), this->endElement() );
920  boost::timer ti0,ti1, ti2, ti3;
921 
922  //double t0 = 0, t1 = 0,t2 = 0,t3 = 0;
923  //
924  // start the real intensive job:
925  // -# iterate over all elements to integrate over
926  // -# construct the associated geometric mapping with the reference element
927  // -# loop over quadrature loop and assemble the local matrix associated with the bilinear form
928  // -# assemble the local contribution in the global representation of the bilinear form
929  //
930  for ( ; it != en; ++it )
931  {
932  size_type idElt = it->id();
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 );
937 
938  auto const& eltTest = __form.testSpace()->mesh()->element( idElt );
939 
940  if ( formc->isZero( idElt ) )
941  continue;
942 
943  switch ( M_gt )
944  {
945  default:
946  case GeomapStrategyType::GEOMAP_HO:
947  {
948  //ti0.restart();
949  __c->update( eltTest );
950  //t0+=ti0.elapsed();
951 #if 0
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";
957 #endif
958  //ti1.restart();
959  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
960  formc->update( mapgmc,mapgmc,mapgmc );
961  //DLOG(INFO) << "update gmc : " << ti1.elapsed() << "\n";
962  //t1+=ti1.elapsed();
963 
964  //ti2.restart();
965  formc->integrate();
966  //DLOG(INFO) << "integrate : " << ti2.elapsed() << "\n";
967  //t2+=ti2.elapsed();
968 
969  //ti3.restart();
970  formc->assemble();
971  //DLOG(INFO) << "assemble : " << ti3.elapsed() << "\n";
972  //t3+=ti3.elapsed();
973  }
974  break;
975 
976  case GeomapStrategyType::GEOMAP_O1:
977  {
978  //ti0.restart();
979  __c1->update( eltTest );
980  //t0+=ti0.elapsed();
981 #if 0
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";
987 #endif
988  //ti1.restart();
989  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
990  formc1->update( mapgmc1,mapgmc1,mapgmc1 );
991  //DLOG(INFO) << "update gmc : " << ti1.elapsed() << "\n";
992  //t1+=ti1.elapsed();
993 
994  //ti2.restart();
995  formc1->integrate();
996  //DLOG(INFO) << "integrate : " << ti2.elapsed() << "\n";
997  //t2+=ti2.elapsed();
998 
999  //ti3.restart();
1000  formc1->assemble();
1001  //DLOG(INFO) << "assemble : " << ti3.elapsed() << "\n";
1002  //t3+=ti3.elapsed();
1003  }
1004  break;
1005 
1006  case GeomapStrategyType::GEOMAP_OPT:
1007  {
1008  if ( it->isOnBoundary() )
1009  {
1010  //ti0.restart();
1011  __c->update( eltTest );
1012  //t0+=ti0.elapsed();
1013 #if 0
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";
1019 #endif
1020  //ti1.restart();
1021  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
1022  formc->update( mapgmc,mapgmc,mapgmc );
1023 
1024  //DLOG(INFO) << "update gmc : " << ti1.elapsed() << "\n";
1025  //t1+=ti1.elapsed();
1026 
1027  //ti2.restart();
1028  formc->integrate();
1029  //DLOG(INFO) << "integrate : " << ti2.elapsed() << "\n";
1030  //t2+=ti2.elapsed();
1031 
1032  //ti3.restart();
1033  formc->assemble();
1034  //DLOG(INFO) << "assemble : " << ti3.elapsed() << "\n";
1035  //t3+=ti3.elapsed();
1036  }
1037 
1038  else
1039  {
1040  //ti0.restart();
1041  __c1->update( eltTest );
1042  //t0+=ti0.elapsed();
1043 #if 0
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";
1049 #endif
1050  //ti1.restart();
1051  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
1052  formc1->update( mapgmc1,mapgmc1,mapgmc1 );
1053 
1054  //DLOG(INFO) << "update gmc : " << ti1.elapsed() << "\n";
1055  //t1+=ti1.elapsed();
1056 
1057  //ti2.restart();
1058  formc1->integrate();
1059  //DLOG(INFO) << "integrate : " << ti2.elapsed() << "\n";
1060  //t2+=ti2.elapsed();
1061 
1062  //ti3.restart();
1063  formc1->assemble();
1064  //DLOG(INFO) << "assemble : " << ti3.elapsed() << "\n";
1065  //t3+=ti3.elapsed();
1066  }
1067  }
1068  break;
1069  }
1070  } // end loop on elements
1071  delete formc;
1072  delete formc1;
1073  }// end loop on list of elements
1074 #if 0
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";
1080 #endif
1081  DLOG(INFO) << "integrating over elements done in " << __timer.elapsed() << "s\n";
1082 
1083 
1084  }
1085 
1086 #if defined( FEELPP_HAS_TBB )
1087 
1088  else
1089  {
1090  element_iterator it = this->beginElement();
1091  element_iterator en = this->endElement();
1092 
1093  if ( it == en )
1094  return;
1095 
1096  std::vector<boost::reference_wrapper<const typename eval::element_type> > _v;
1097 
1098  for ( auto _it = it; _it != en; ++_it )
1099  _v.push_back( boost::cref( *_it ) );
1100 
1101  //tbb::blocked_range<decltype(_v.begin())> r( _v.begin(), _v.end(), M_grainsize );
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,
1104  this->expression(),
1105  this->im(),
1106  *it );
1107  tbb::parallel_for( r, thecontext );
1108  }
1109 
1110 #endif // FEELPP_HAS_TBB
1111 }
1112 
1113 template<typename Elements, typename Im, typename Expr, typename Im2>
1114 template<typename FormType>
1115 void
1116 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_ELEMENTS> , mpl::bool_<false> , bool hasRelation ) const
1117 {
1118  if ( hasRelation )
1119  {
1120  assembleWithRelationDifferentMeshType( __form,mpl::int_<MESH_ELEMENTS>() );
1121  }
1122  else
1123  {
1124  assembleInCaseOfInterpolate( __form,mpl::int_<MESH_ELEMENTS>() );
1125  }
1126 }
1127 
1128 namespace detail
1129 {
1130 template<typename SpaceType,typename ImType,typename GmcType,typename GmcExprType>
1131 boost::shared_ptr<GmcType>
1132 buildGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType> const& /*space*/,typename GmcType::gm_ptrtype /*gm*/,
1133  size_type /*idElt*/ ,boost::shared_ptr<GmcExprType> gmcExpr,mpl::int_<0> )
1134 {
1135  return gmcExpr;
1136 }
1137 
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> )
1142 {
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;
1149 
1150  ImType im;
1151  QuadMapped<ImType> qm;
1152  permutation_points_type ppts( qm( im ) );
1153  std::vector<std::map<permutation_type, geopc_ptrtype> > __geopc( im.nFaces() );
1154 
1155  for ( uint16_type __f = 0; __f < im.nFaces(); ++__f )
1156  {
1157  for ( permutation_type __p( permutation_type::IDENTITY );
1158  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1159  {
1160  __geopc[__f][__p] = geopc_ptrtype( new geopc_type( gm, ppts[__f].find( __p )->second ) );
1161  }
1162  }
1163 
1164  auto const& eltInit = space->mesh()->face( idElt );
1165  gmc_ptrtype gmc( new gmc_type( gm, eltInit.element0(), __geopc, eltInit.pos_first() /*face_id_in_elt_0*/ ) );
1166 
1167  return gmc;
1168 }
1169 template<typename SpaceType,typename ImType,typename GmcType,typename GmcExprType>
1170 void
1171 updateGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType> const& /*space*/,
1172  boost::shared_ptr<GmcType> /*gmc*/, boost::shared_ptr<GmcExprType> /*gmcExpr*/,
1173  size_type /*idElt*/, mpl::int_<0> )
1174 {
1175  // nothing to do!
1176 }
1177 template<typename SpaceType,typename ImType,typename GmcType,typename GmcExprType>
1178 void
1179 updateGmcWithRelationDifferentMeshType( boost::shared_ptr<SpaceType> const& space,
1180  boost::shared_ptr<GmcType> gmc, boost::shared_ptr<GmcExprType> gmcExpr,
1181  size_type idElt, mpl::int_<1> )
1182 {
1183  typedef typename QuadMapped<ImType>::permutation_type permutation_type;
1184 
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 )
1189  {
1190  gmc->update( theface.element0(), theface.pos_first(), __p );
1191 
1192  bool check=true;
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 );
1196 
1197  if (check) findPermutation=true;
1198  }
1199  CHECK(findPermutation) << "the permutation of quad point is not find\n";
1200 }
1201 
1202 } // namespace detail
1203 
1204 template<typename Elements, typename Im, typename Expr, typename Im2>
1205 template<typename FE1,typename FE2,typename ElemContType>
1206 void
1207 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1208  mpl::int_<MESH_ELEMENTS> ) const
1209 {
1210  LOG(INFO) << "[integrator::assembleWithRelationDifferentMeshType] vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS>\n";
1211 
1212  // typedef on integral mesh (expr) :
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;
1219 
1220  // typedef on form (trial and test):
1221  typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
1222  // test
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;
1228 
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;
1233  // trial
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;
1243 
1244  // typedef on formcontext
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;
1248 
1249  //static const bool gmTestIsGmExpr = boost::is_same<gm_expr_type,gm_formTest_type>::type::value;
1250  //static const bool gmTrialIsGmExpr = boost::is_same<gm_expr_type,gm_formTrial_type>::type::value;
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;
1256 
1257 
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;
1262 
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;
1267 
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;
1272 
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;
1277 
1278 
1279  //-----------------------------------------------//
1280 
1281  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
1282  {
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";
1286 
1287  // check that we have elements to iterate over
1288  if ( elt_it == elt_en )
1289  continue;
1290 
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 );
1302  //-----------------------------------------------//
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 ) );
1305  //-----------------------------------------------//
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>() );
1312  //-----------------------------------------------//
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 ) );
1316  //-----------------------------------------------//
1317  im_formtest_type imTest;
1318  im_formtrial_type imTrial;
1319  focb_ptrtype formc( new form_context_type( __form,
1320  mapgmcFormTest,
1321  mapgmcFormTrial,
1322  mapgmcExpr,
1323  this->expression(),
1324  this->im(),imTest,imTrial ) );
1325 
1326  //-----------------------------------------------//
1327 
1328  for ( ; elt_it != elt_en; ++elt_it )
1329  {
1330  const size_type idEltRange = elt_it->id();
1331  size_type idEltTest = idEltRange;
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 );
1336  size_type idEltTrial = 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 );
1341 
1342 #if 0
1343  if ( formc->isZero( eltTest.element0() /*idElt*/ ) )
1344  continue;
1345 #endif
1346 
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>() );
1354 
1355  formc->update( mapgmcFormTest,mapgmcFormTrial,mapgmcExpr );
1356  formc->integrate();
1357  formc->assemble();
1358 
1359  } // end loop on elements
1360  delete formc;
1361  //delete formc1;
1362  } // end loop on list of elements
1363 }
1364 
1365 template<typename Elements, typename Im, typename Expr, typename Im2>
1366 template<typename FE,typename VectorType,typename ElemContType>
1367 void
1368 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> ) const
1369 {
1370  CHECK ( false ) << "[assembleWithRelationDifferentMeshType<LinearForm,MESH_ELEMENTS>] : not implement\n";
1371 }
1372 
1373 template<typename Elements, typename Im, typename Expr, typename Im2>
1374 template<typename FE1,typename FE2,typename ElemContType>
1375 void
1376 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1377  mpl::int_<MESH_ELEMENTS> ) const
1378 {
1379  LOG(INFO) << "[integrator::assembleInCaseOfInterpolate] vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form, mpl::int_<MESH_ELEMENTS>\n";
1380 
1381  // typedef on integral mesh (expr) :
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;
1388 
1389  // typedef on form (trial and test):
1390  typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
1391  // test
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;
1399  // trial
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;
1407 
1408  // typedef on formcontext
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;
1412 
1413  //-----------------------------------------------//
1414 
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 )
1418  {
1419  elt_it = lit->template get<1>();
1420  elt_en = lit->template get<2>();
1421 
1422  // check that we have elements to iterate over
1423  if ( elt_it == elt_en )
1424  continue;
1425  else
1426  findEltForInit = true;
1427  }
1428  if (!findEltForInit) return;
1429 
1430  //-----------------------------------------------//
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 ) );
1434  //-----------------------------------------------//
1435  pc_formTest_ptrtype geopcFormTest( new pc_formTest_type( __form.gm(), __form.testSpace()->fe()->points()/* this->im().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 ) );
1438  //-----------------------------------------------//
1439  pc_formTrial_ptrtype geopcFormTrial( new pc_formTrial_type( __form.gmTrial(), __form.trialSpace()->fe()->points() /*this->im().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 ) );
1442  //-----------------------------------------------//
1443 
1444  focb_ptrtype formc( new form_context_type( __form,
1445  mapgmcFormTest,
1446  mapgmcFormTrial,
1447  mapgmcExpr,
1448  this->expression(),
1449  this->im() ) );
1450 
1451  //-----------------------------------------------//
1452 
1453  auto meshTrial = __form.trialSpace()->mesh();
1454  auto meshTest = __form.testSpace()->mesh();
1455 
1456  if (!M_QPL)
1457  {
1458  M_QPL.reset( new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
1459  M_QPL->update( meshTest,meshTrial );
1460  }
1461 
1462  if (!M_QPL->hasPrecompute())
1463  {
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 )
1467  {
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 )
1473  {
1474  auto const idEltTrial = map_it->first;
1475  auto const& eltTrial = meshTrial->element( idEltTrial );
1476  auto const& eltTest = meshTest->element( idEltTest );
1477 
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 )
1486  {
1487  geopcFormTest->update( gmcExpr_it->template get<2>() );
1488  geopcFormTrial->update( gmcExpr_it->template get<3>() );
1489 
1490  gmcFormTest->update( eltTest,geopcFormTest );
1491  gmcFormTrial->update( eltTrial,geopcFormTrial );
1492 
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>() ) );
1496 
1497  formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr, gmcExpr_it->template get<0>() );
1498 
1499  formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
1500  isFirstExperience = false;
1501  }
1502 
1503  formc->assembleInCaseOfInterpolate();
1504 
1505  }
1506  }
1507  } // if (!M_QPL->hasPrecomputeBF())
1508  else
1509  {
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)
1514  {
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)
1519  {
1520  map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>()/*gmcFormTest*/ ) );
1521  map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<3>()/*gmcFormTrial*/ ) );
1522  map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<1>() ) );
1523 
1524  formc->updateInCaseOfInterpolate( mapgmcFormTest, mapgmcFormTrial, mapgmcExpr, resQPLloc_it->template get<0>() );
1525 
1526  formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
1527  isFirstExperience = false;
1528  }
1529 
1530  formc->assembleInCaseOfInterpolate();
1531  }
1532 
1533  } //else !M_QPL->hasPrecomputeBF()
1534 
1535  delete formc;
1536 
1537 } // assembleInCaseOfInterpolate
1538 
1539 
1540 template<typename Elements, typename Im, typename Expr, typename Im2>
1541 template<typename FE,typename VectorType,typename ElemContType>
1542 void
1543 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_ELEMENTS> ) const
1544 {
1545 
1546  // typedef on integral mesh (expr) :
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;
1553 
1554  //typedef on form (trial and test):
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;
1563 
1564  // typedef on formcontext
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;
1568 
1569  //-----------------------------------------------//
1570 
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 )
1574  {
1575  elt_it = lit->template get<1>();
1576  elt_en = lit->template get<2>();
1577 
1578  // check that we have elements to iterate over
1579  if ( elt_it == elt_en )
1580  continue;
1581  else
1582  findEltForInit = true;
1583  }
1584  if (!findEltForInit) return;
1585 
1586  //-----------------------------------------------//
1587 
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 ) );
1591 
1592  //-----------------------------------------------//
1593 
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 ) );
1597 
1598  //-----------------------------------------------//
1599 
1600  focb_ptrtype formc( new form_context_type( __form,
1601  mapgmcForm,
1602  mapgmcForm,
1603  mapgmcExpr,
1604  this->expression(),
1605  this->im() ) );
1606 
1607  //-----------------------------------------------//
1608 
1609  auto meshTest = __form.testSpace()->mesh();
1610 
1611  if (!M_QPL)
1612  {
1613  M_QPL.reset( new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
1614  M_QPL->update( meshTest );
1615  }
1616 
1617  //-----------------------------------------------//
1618 
1619  if (!M_QPL->hasPrecompute())
1620  {
1621  auto res_it = M_QPL->resultLinear().begin();
1622  auto const res_en = M_QPL->resultLinear().end();
1623 
1624  for ( ; res_it != res_en ; ++res_it )
1625  {
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;
1634 
1635  for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
1636  {
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>() );
1642 
1643  formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
1644  isFirstExperience = false;
1645  }
1646 
1647  formc->assemble();
1648  }
1649  } //if (!M_QPL->hasPrecompute())
1650  else
1651  {
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)
1656  {
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)
1661  {
1662  map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>()/*gmcForm*/ ) );
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>() );
1665 
1666  formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
1667  isFirstExperience = false;
1668  }
1669  formc->assemble();
1670  }
1671  } //else
1672 
1673  delete formc;
1674 
1675 }
1676 
1677 
1678 template<typename Elements, typename Im, typename Expr, typename Im2>
1679 template<typename FormType>
1680 void
1681 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<true> , bool /*hasRelation*/ ) const
1682 {
1683  DLOG(INFO) << "integrating over "
1684  << std::distance( this->beginElement(), this->endElement() ) << " faces\n";
1685  boost::timer __timer;
1686 
1687  //
1688  // some typedefs
1689  //
1690  typedef typename eval::gm_type gm_type;
1691  typedef typename eval::gm1_type gm1_type;
1692  //typedef typename FormType::gm_type gm_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;
1699  //
1700  // Precompute some data in the reference element for
1701  // geometric mapping and reference finite element
1702  //
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;
1707  //typedef typename mpl::if_<mpl::equal_to<mpl::int_<FormType::nDim>, mpl::int_<2> >, mpl::identity<typename eval::element_type::edge_permutation_type>, mpl::identity<typename eval::element_type::face_permutation_type> >::type::type permutation_type;
1708 
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() ) );
1712 
1713 
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;
1718 
1719  //__typeof__(im(__face_id_in_elt_0 ) ) im_face ( im(__face_id_in_elt_0 ) );
1720  std::vector<face_im_type> face_ims( im().nFaces() );
1721  std::vector<face_im2_type> face_ims2( im2().nFaces() );
1722 
1723  for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
1724  {
1725  face_ims[__f] = this->im( __f );
1726  face_ims2[__f] = this->im2( __f );
1727 
1728  for ( permutation_type __p( permutation_type::IDENTITY );
1729  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1730  {
1731  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
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 ) );
1734  }
1735  }
1736 
1737  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
1738  {
1739  auto it = lit->template get<1>();
1740  auto en = lit->template get<2>();
1741 
1742  // check that we have elements to iterate over
1743  if ( (it == en) || (it->isConnectedTo0() == false) )
1744  continue;
1745 
1746  uint16_type __face_id_in_elt_0 = it->pos_first();
1747 
1748  // get the geometric mapping associated with element 0
1749  //DLOG(INFO) << "element " << it->element(0) << "face " << __face_id_in_elt_0 << " permutation " << it->element(0).permutation( __face_id_in_elt_0 ) << "\n";
1750  gm_ptrtype __gm = it->element( 0 ).gm();
1751  gm1_ptrtype __gm1 = it->element( 0 ).gm1();
1752  //DLOG(INFO) << "[integrator] evaluate(faces), gm is cached: " << __gm->isCached() << "\n";
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 ) );
1755 
1756 
1757 
1758  //
1759  // the case where the face is connected only to one element
1760  //
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;
1771 
1772  //
1773  // the case where the face is connected only to two elements
1774  //
1775  // get the geometric mapping associated with element 1
1776  gmc_ptrtype __c1;
1777  gmc1_ptrtype __c11;
1778 
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;
1787 
1788  bool isInitConnectionTo0=false;
1789  bool isInitConnectionTo1=false;
1790 
1791  // true if connected to another element, false otherwise
1792  if ( it->isConnectedTo1() )
1793  {
1794  uint16_type __face_id_in_elt_1 = it->pos_second();
1795 
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 ) );
1802 
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;
1806  }
1807 
1808  else
1809  {
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;
1813  }
1814 
1815  boost::timer ti0,ti1, ti2, ti3;
1816  //double t0 = 0, t1 = 0,t2 = 0,t3 = 0;
1817  DLOG(INFO) << "[Integrator::faces/forms] starting...\n";
1818 
1819  //
1820  // start the real intensive job:
1821  // -# iterate over all elements to integrate over
1822  // -# construct the associated geometric mapping with the reference element
1823  // -# loop over quadrature loop and assemble the local matrix associated with the bilinear form
1824  // -# assemble the local contribution in the global representation of the bilinear form
1825  //
1826  for ( ; it != en; ++it )
1827  {
1828  if ( it->isGhostFace() )
1829  {
1830  LOG(WARNING) << "face id : " << it->id() << " is a ghost face";
1831  continue;
1832  }
1833  if ( it->isConnectedTo1() )
1834  {
1835  if ( !isInitConnectionTo1 )
1836  {
1837  uint16_type __face_id_in_elt_1 = it->pos_second();
1838 
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 ) );
1845 
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;
1849  }
1850 
1851  switch ( M_gt )
1852  {
1853  default:
1854  case GeomapStrategyType::GEOMAP_HO:
1855  {
1856  FEELPP_ASSERT( it->isOnBoundary() == false )
1857  ( it->id() ).error( "face on boundary but connected on both sides" );
1858  //ti0.restart();
1859  // get the id of the face in each adjacent element
1860  uint16_type __face_id_in_elt_0 = it->pos_first();
1861  uint16_type __face_id_in_elt_1 = it->pos_second();
1862 
1863  __c0->update( it->element( 0 ), __face_id_in_elt_0 );
1864  __c1->update( it->element( 1 ), __face_id_in_elt_1 );
1865  //t0 += ti0.elapsed();
1866 
1867  //ti1.restart();
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>() );
1871  //t1 += ti1.elapsed();
1872 
1873  //ti2.restart();
1874  form2->integrate( );
1875  //t2 += ti2.elapsed();
1876 
1877  //ti3.restart();
1878  form2->assemble( it->element( 0 ).id(), it->element( 1 ).id() );
1879  //t3 += ti3.elapsed();
1880  }
1881  break;
1882 
1883  case GeomapStrategyType::GEOMAP_O1:
1884  case GeomapStrategyType::GEOMAP_OPT:
1885  {
1886  FEELPP_ASSERT( it->isOnBoundary() == false )
1887  ( it->id() ).error( "face on boundary but connected on both sides" );
1888  //ti0.restart();
1889  // get the id of the face in each adjacent element
1890  uint16_type __face_id_in_elt_0 = it->pos_first();
1891  uint16_type __face_id_in_elt_1 = it->pos_second();
1892 
1893  __c01->update( it->element( 0 ), __face_id_in_elt_0 );
1894  __c11->update( it->element( 1 ), __face_id_in_elt_1 );
1895  //t0 += ti0.elapsed();
1896 
1897  //ti1.restart();
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>() );
1901  //t1 += ti1.elapsed();
1902 
1903  //ti2.restart();
1904  form21->integrate( );
1905  //t2 += ti2.elapsed();
1906 
1907  //ti3.restart();
1908  form21->assemble( it->element( 0 ).id(), it->element( 1 ).id() );
1909  //t3 += ti3.elapsed();
1910  }
1911  break;
1912  }
1913  }
1914 
1915  else
1916  {
1917 #if 1
1918  uint16_type __face_id_in_elt_0 = it->pos_first();
1919 
1920  if ( !isInitConnectionTo0 )
1921  {
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;
1925  }
1926 
1927  //ti0.restart();
1928  __c0->update( it->element( 0 ),__face_id_in_elt_0 );
1929  //t0 += ti0.elapsed();
1930 
1931  FEELPP_ASSERT( __face_id_in_elt_0 == __c0->faceId() )
1932  ( __face_id_in_elt_0 )
1933  ( __c0->faceId() ).warn ( "invalid face id" );
1934 
1935  //ti1.restart();
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] );
1938  //t1 += ti1.elapsed();
1939 
1940  //ti2.restart();
1941  form->integrate( );
1942  //t2 += ti2.elapsed();
1943 
1944  //ti3.restart();
1945  form->assemble( it->element( 0 ).id() );
1946  //t3 += ti3.elapsed();
1947 #else
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] );
1950  form->integrate( );
1951 #endif
1952  } // end loop on elements
1953  }
1954  }// end loop on list of element
1955 
1956 #if 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";
1962 #endif
1963  DLOG(INFO) << "integrating over faces done in " << __timer.elapsed() << "s\n";
1964  //std::cout << "integrating over faces done in " << __timer.elapsed() << "s\n";
1965 }
1966 
1967 template<typename Elements, typename Im, typename Expr, typename Im2>
1968 template<typename FormType>
1969 void
1970 Integrator<Elements, Im, Expr, Im2>::assemble( FormType& __form, mpl::int_<MESH_FACES> , mpl::bool_<false> , bool hasRelation ) const
1971 {
1972  if ( false /*hasRelation*/ )
1973  {
1974  assembleWithRelationDifferentMeshType( __form,mpl::int_<MESH_FACES>() );
1975  }
1976  else
1977  {
1978  assembleInCaseOfInterpolate( __form,mpl::int_<MESH_FACES>() );
1979  }
1980 }
1981 
1982 template<typename Elements, typename Im, typename Expr, typename Im2>
1983 template<typename FE1,typename FE2,typename ElemContType>
1984 void
1985 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
1986  mpl::int_<MESH_FACES> ) const
1987 {
1988  CHECK ( false ) << "[assembleWithRelationDifferentMeshType<BilinearForm,MESH_FACES>] : not implement\n";
1989 }
1990 
1991 template<typename Elements, typename Im, typename Expr, typename Im2>
1992 template<typename FE,typename VectorType,typename ElemContType>
1993 void
1994 Integrator<Elements, Im, Expr, Im2>::assembleWithRelationDifferentMeshType(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> ) const
1995 {
1996  CHECK ( false ) << "[assembleWithRelationDifferentMeshType<LinearForm,MESH_FACES>] : not implement\n";
1997 }
1998 
1999 template<typename Elements, typename Im, typename Expr, typename Im2>
2000 template<typename FE1,typename FE2,typename ElemContType>
2001 void
2002 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::BilinearForm<FE1,FE2,ElemContType>& __form,
2003  mpl::int_<MESH_FACES> ) const
2004 {
2005 
2006  // typedef on integral mesh (expr) :
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;
2013 
2014  //typedef on form (trial and test):
2015  typedef vf::detail::BilinearForm<FE1,FE2,ElemContType> FormType;
2016  /*typedef typename FormType::gm_type gm_form_type;
2017  typedef typename FormType::mesh_element_type geoelement_form_type;
2018  typedef typename gm_form_type::template Context<expression_type::context|vm::POINT,geoelement_form_type> gmc_form_type;
2019  typedef boost::shared_ptr<gmc_form_type> gmc_form_ptrtype;
2020  typedef typename gm_form_type::precompute_type pc_form_type;
2021  typedef typename gm_form_type::precompute_ptrtype pc_form_ptrtype;
2022  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_form_ptrtype> > map_gmc_form_type;*/
2023 
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;
2031  // trial
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;
2039 
2040 
2041 
2042  // typedef on formcontext
2043  typedef typename im_type::face_quadrature_type face_im_type;
2044 
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;
2048 
2049 
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 )
2053  {
2054  elt_it = lit->template get<1>();
2055  elt_en = lit->template get<2>();
2056 
2057  // check that we have elements to iterate over
2058  if ( elt_it == elt_en )
2059  continue;
2060  else
2061  findEltForInit = true;
2062  }
2063  if (!findEltForInit) return;
2064 
2065 
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() ) );
2069 
2070  std::vector<std::map<permutation_type, pc_expr_ptrtype> > __geopcExpr( im().nFaces() );
2071  std::vector<face_im_type> face_ims( im().nFaces() );
2072 
2073  for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2074  {
2075  face_ims[__f] = this->im( __f );
2076 
2077  for ( permutation_type __p( permutation_type::IDENTITY );
2078  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2079  {
2080  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
2081  __geopcExpr[__f][__p] = pc_expr_ptrtype( new pc_expr_type( elt_it->element( 0 ).gm(), ppts[__f].find( __p )->second ) );
2082  }
2083  }
2084 
2085 
2086  uint16_type __face_id_in_elt_0 = elt_it->pos_first();
2087 
2088  gmc_expr_ptrtype gmcExpr( new gmc_expr_type( elt_it->element( 0 ).gm(),
2089  elt_it->element( 0 ),
2090  __geopcExpr,
2091  __face_id_in_elt_0 ) );
2092 
2093 
2094  map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
2095 
2096  //-----------------------------------------------//
2097 
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 ) );
2101 
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 ) );
2105 
2106  //-----------------------------------------------//
2107 
2108  focb_ptrtype formc( new form_context_type( __form,
2109  mapgmcFormTest,
2110  mapgmcFormTrial,
2111  mapgmcExpr,
2112  this->expression(),
2113  face_ims[__face_id_in_elt_0],
2114  this->im() ) );
2115 
2116  //-----------------------------------------------//
2117 
2118  auto meshTrial = __form.trialSpace()->mesh();
2119  auto meshTest = __form.testSpace()->mesh();
2120 
2121  if (!M_QPL)
2122  {
2123  M_QPL.reset( new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
2124  M_QPL->update( meshTest,meshTrial );
2125  }
2126 
2127  //-----------------------------------------------//
2128 
2129  if (!M_QPL->hasPrecompute())
2130  {
2131  auto res_it = M_QPL->result().begin();
2132  auto res_en = M_QPL->result().end();
2133 
2134  for ( ; res_it != res_en ; ++res_it )
2135  {
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();
2140 
2141  for ( ; map_it != map_en ; ++map_it )
2142  {
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>();
2149 
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;
2154 
2155  for ( ; gmcExpr_it != gmcExpr_en ; ++gmcExpr_it )
2156  {
2157  geopcFormTest->update( gmcExpr_it->template get<2>() );
2158  geopcFormTrial->update( gmcExpr_it->template get<3>() );
2159 
2160  gmcFormTest->update( eltTest,geopcFormTest );
2161  gmcFormTrial->update( eltTrial,geopcFormTrial );
2162  //std::cout << "eltTest.id() " << eltTest.id() << " eltTest.G() " << eltTest.G()
2163  // << " eltTrial.id() "<< eltTrial.id() << " eltTrial.G() " << eltTrial.G() << std::endl;
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>() );
2169 
2170  formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
2171  isFirstExperience = false;
2172  }
2173 
2174  formc->assembleInCaseOfInterpolate();
2175  }
2176 
2177  }
2178  } // if (!M_QPL->hasPrecompute())
2179  else
2180  {
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)
2185  {
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)
2190  {
2191  map_gmc_formTest_type mapgmcFormTest( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>()/*gmcFormTest*/ ) );
2192  map_gmc_formTrial_type mapgmcFormTrial( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<3>()/*gmcFormTrial*/ ) );
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>() );
2196 
2197  formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
2198  isFirstExperience = false;
2199  }
2200 
2201  formc->assembleInCaseOfInterpolate();
2202  }
2203  } //else
2204 
2205  delete formc;
2206 }
2207 
2208 template<typename Elements, typename Im, typename Expr, typename Im2>
2209 template<typename FE,typename VectorType,typename ElemContType>
2210 void
2211 Integrator<Elements, Im, Expr, Im2>::assembleInCaseOfInterpolate(vf::detail::LinearForm<FE,VectorType,ElemContType>& __form, mpl::int_<MESH_FACES> ) const
2212 {
2213  // typedef on integral mesh (expr) :
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;
2220 
2221  //typedef on form (test):
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;
2230 
2231  // typedef on formcontext
2232  typedef typename im_type::face_quadrature_type face_im_type;
2233 
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;
2237 
2238 
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 )
2242  {
2243  elt_it = lit->template get<1>();
2244  elt_en = lit->template get<2>();
2245 
2246  // check that we have elements to iterate over
2247  if ( elt_it == elt_en )
2248  continue;
2249  else
2250  findEltForInit = true;
2251  }
2252  if (!findEltForInit) return;
2253 
2254  //-----------------------------------------------//
2255 
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() ) );
2259 
2260  std::vector<std::map<permutation_type, pc_expr_ptrtype> > __geopcExpr( im().nFaces() );
2261  std::vector<face_im_type> face_ims( im().nFaces() );
2262 
2263  for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2264  {
2265  face_ims[__f] = this->im( __f );
2266 
2267  for ( permutation_type __p( permutation_type::IDENTITY );
2268  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2269  {
2270  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
2271  __geopcExpr[__f][__p] = pc_expr_ptrtype( new pc_expr_type( elt_it->element( 0 ).gm(), ppts[__f].find( __p )->second ) );
2272  }
2273  }
2274 
2275 
2276  uint16_type __face_id_in_elt_0 = elt_it->pos_first();
2277 
2278  gmc_expr_ptrtype gmcExpr( new gmc_expr_type( elt_it->element( 0 ).gm(),
2279  elt_it->element( 0 ),
2280  __geopcExpr,
2281  __face_id_in_elt_0 ) );
2282 
2283  map_gmc_expr_type mapgmcExpr( fusion::make_pair<vf::detail::gmc<0> >( gmcExpr ) );
2284 
2285  //-----------------------------------------------//
2286 
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 ) );
2290 
2291  //-----------------------------------------------//
2292 
2293  focb_ptrtype formc( new form_context_type( __form,
2294  mapgmcForm,
2295  mapgmcForm,
2296  mapgmcExpr,
2297  this->expression(),
2298  face_ims[__face_id_in_elt_0],
2299  this->im() ) );
2300 
2301  //-----------------------------------------------//
2302 
2303  auto meshTest = __form.testSpace()->mesh();
2304 
2305  if (!M_QPL)
2306  {
2307  M_QPL.reset( new QuadPtLocalization<Elements, Im, Expr >( M_elts ) );
2308  M_QPL->update( meshTest );
2309  }
2310 
2311  //-----------------------------------------------//
2312 
2313  if (!M_QPL->hasPrecompute())
2314  {
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 )
2318  {
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>();
2323 
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 )
2329  {
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>() );
2336 
2337  formc->integrateInCaseOfInterpolate( gmcExpr_it->template get<0>(),isFirstExperience );
2338  isFirstExperience = false;
2339  }
2340 
2341  formc->assemble();
2342  }
2343  }
2344  else
2345  {
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)
2350  {
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)
2355  {
2356  map_gmc_form_type mapgmcForm( fusion::make_pair<vf::detail::gmc<0> >( resQPLloc_it->template get<2>()/*gmcForm*/ ) );
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>() );
2360 
2361  formc->integrateInCaseOfInterpolate( resQPLloc_it->template get<0>(),isFirstExperience );
2362  isFirstExperience = false;
2363  }
2364  formc->assemble();
2365  }
2366  } // else
2367 
2368  delete formc;
2369 }
2370 
2371 
2372 
2373 
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
2377 {
2378  DLOG(INFO) << "integrating over "
2379  << std::distance( this->beginElement(), this->endElement() ) << " elements\n";
2380  boost::timer __timer;
2381 
2382 #if defined(FEELPP_HAS_TBB)
2383 
2384  //std::cout << "Integrator Uses TBB: " << M_use_tbb << "\n";
2385  if ( !M_use_tbb )
2386 #else
2387  if ( 1 )
2388 #endif
2389  {
2390  //
2391  // some typedefs
2392  //
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;
2400 
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;
2405 
2406  //typedef typename eval_expr_type::value_type value_type;
2407  //typedef typename Im::value_type value_type;
2408  typename eval::matrix_type res( eval::matrix_type::Zero() );
2409 
2410  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2411  {
2412  auto it = lit->template get<1>();
2413  auto en = lit->template get<2>();
2414 
2415  // make sure that we have elements to iterate over (return 0
2416  // otherwise)
2417  if ( it == en )
2418  continue;
2419 
2420  //std::cout << "0" << std::endl;
2421 
2422  //
2423  // Precompute some data in the reference element for
2424  // geometric mapping and reference finite element
2425  //
2426  // warning this is not efficient here, we want to use the geometric mapping
2427  // from the elements in order to take advantage of the cache if possible
2428  // this change hsa been made in order to circumvent a bug which is not yet found
2429  //#warning INEFFICIENT CODE HERE : TO DEBUG
2430  //gm_ptrtype gm( new gm_type) ;//it->gm();
2431  gm_ptrtype gm( it->gm() );
2432  //std::cout << "0.5" << std::endl;
2433  gm1_ptrtype gm1( new gm1_type ); //it->gm1();
2434  //std::cout << "0.6: " << gm1.use_count() << " " << gm.use_count() << std::endl;
2435  //DDLOG(INFO) << "[integrator] evaluate(elements), gm is cached: " << gm->isCached() << "\n";
2436  typename eval::gmpc_ptrtype __geopc( new typename eval::gmpc_type( gm,
2437  this->im().points() ) );
2438  //std::cout << "1" << std::endl;
2439  typename eval::gmpc1_ptrtype __geopc1( new typename eval::gmpc1_type( gm1,
2440  this->im().points() ) );
2441 
2442  //std::cout << "2" << std::endl;
2443  //it = this->beginElement();
2444 
2445  // wait for all the guys
2446 #ifdef FEELPP_HAS_MPI
2447  auto const& worldComm = const_cast<MeshBase*>( it->mesh() )->worldComm();
2448 #if 0
2449  if ( worldComm.localSize() > 1 )
2450  {
2451  worldComm.localComm().barrier();
2452  }
2453 #endif
2454 #endif
2455 
2456  // possibly high order
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 ) );
2460  //std::cout << "3" << std::endl;
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;
2464  //std::cout << "4" << std::endl;
2465 
2466  // order 1
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 ) );
2470  //std::cout << "5" << std::endl;
2471  typedef typename expression_type::template tensor<map_gmc1_type> eval_expr1_type;
2472  eval_expr1_type expr1( expression(), mapgmc1 );
2473 
2474  //std::cout << "6" << std::endl;
2475 
2476 
2477  //value_type res1 = 0;
2478  for ( ; it != en; ++it )
2479  {
2480  switch ( M_gt )
2481  {
2482  default:
2483  case GeomapStrategyType::GEOMAP_HO :
2484  {
2485  __c->update( *it );
2486  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2487  expr.update( mapgmc );
2488  const gmc_type& gmc = *__c;
2489 
2490  M_im.update( gmc );
2491 
2492 
2493  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2494  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2495  {
2496  res( c1,c2 ) += M_im( expr, c1, c2 );
2497  }
2498  }
2499  break;
2500 
2501  case GeomapStrategyType::GEOMAP_O1:
2502  {
2503  //DDLOG(INFO) << "geomap o1" << "\n";
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;
2508 
2509  M_im.update( gmc );
2510 
2511 
2512  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2513  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2514  {
2515  res( c1,c2 ) += M_im( expr1, c1, c2 );
2516  }
2517  }
2518  break;
2519 
2520  case GeomapStrategyType::GEOMAP_OPT:
2521  {
2522  //DDLOG(INFO) << "geomap opt" << "\n";
2523  if ( it->isOnBoundary() )
2524  {
2525  //DDLOG(INFO) << "boundary element using ho" << "\n";
2526  __c->update( *it );
2527  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2528  expr.update( mapgmc );
2529  const gmc_type& gmc = *__c;
2530 
2531  M_im.update( gmc );
2532 
2533 
2534  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2535  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2536  {
2537  res( c1,c2 ) += M_im( expr, c1, c2 );
2538  }
2539  }
2540 
2541  else
2542  {
2543  //DDLOG(INFO) << "interior element using order 1" << "\n";
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;
2548 
2549  M_im.update( gmc );
2550 
2551  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2552  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2553  {
2554  res( c1,c2 ) += M_im( expr1, c1, c2 );
2555  }
2556  }
2557  }
2558 
2559  //break;
2560  }
2561  }
2562  }
2563  DLOG(INFO) << "integrating over elements done in " << __timer.elapsed() << "s\n";
2564  return res;
2565  }
2566 
2567  else
2568  {
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;
2573 
2574  if ( it == en )
2575  return eval::zero();
2576 
2577  std::vector<boost::reference_wrapper<const typename eval::element_type> > _v;
2578 
2579  for ( auto _it = it; _it != en; ++_it )
2580  _v.push_back( boost::cref( *_it ) );
2581 
2582  tbb::blocked_range<decltype( _v.begin() )> r( _v.begin(), _v.end(), M_grainsize );
2583  context_type thecontext( this->expression(), this->im(), *it );
2584 
2585  if ( M_partitioner == "auto" )
2586  tbb::parallel_reduce( r, thecontext );
2587 
2588  else if ( M_partitioner == "simple" )
2589  tbb::parallel_reduce( r, thecontext, tbb::simple_partitioner() );
2590 
2591  //else if ( M_partitioner == "affinity" )
2592  //tbb::parallel_reduce( r, thecontext, tbb::affinity_partitioner() );
2593  return thecontext.result();
2594 #endif // FEELPP_HAS_TBB
2595  }
2596 
2597 }
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
2601 {
2602  DLOG(INFO) << "integrating over "
2603  << std::distance( this->beginElement(), this->endElement() ) << "faces\n";
2604  boost::timer __timer;
2605 
2606  //
2607  // some typedefs
2608  //
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;
2616  //typedef typename eval_expr_type::value_type value_type;
2617  //typedef typename Im::value_type value_type;
2618 
2619  //BOOST_MPL_ASSERT_MSG( the_element_type::nDim > 1, INVALID_DIM, (mpl::int_<the_element_type::nDim>, mpl::int_<the_face_element_type::nDim>, mpl::identity<the_face_element_type>, mpl::identity<the_element_type> ) );;
2620 
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;
2624 
2625  //
2626  // Precompute some data in the reference element for
2627  // geometric mapping and reference finite element
2628  //
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;
2633 
2634  std::vector<im_face_type> __integrators;
2635 
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() );
2639 
2640 
2641  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2642  {
2643  auto it = lit->template get<1>();
2644  auto en = lit->template get<2>();
2645 
2646  // make sure that we have elements to iterate over (return 0
2647  // otherwise)
2648  if ( (it == en) || (it->isConnectedTo0()==false) )
2649  continue;
2650  //return typename eval::matrix_type( eval::matrix_type::Zero() );
2651 
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();
2656 
2657  //DDLOG(INFO) << "[integrator] evaluate(faces), gm is cached: " << gm->isCached() << "\n";
2658  for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2659  {
2660  __integrators.push_back( im( __f ) );
2661 
2662  for ( permutation_type __p( permutation_type::IDENTITY );
2663  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2664  {
2665  //FEELPP_ASSERT( ppts[__f][__p]->size2() != 0 ).warn( "invalid quadrature type" );
2666  __geopc[__f][__p] = pc_ptrtype( new pc_type( gm, ppts[__f].find( __p )->second ) );
2667  }
2668  }
2669 
2670  uint16_type __face_id_in_elt_0 = it->pos_first();
2671 
2672  // get the geometric mapping associated with element 0
2673  gmc_ptrtype __c0( new gmc_type( gm, it->element( 0 ), __geopc, __face_id_in_elt_0 ) );
2674 
2675  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2676 
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 ) );
2681  expr->init( im() );
2682 
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;
2687 
2688  // true if connected to another element, false otherwise
2689  bool isConnectedTo1 = it->isConnectedTo1();
2690 
2691  // get the geometric mapping associated with element 1
2692  gmc_ptrtype __c1;
2693 
2694  //value_type res = 0;
2695  //value_type res1 = 0;
2696  if ( isConnectedTo1 )
2697  {
2698  uint16_type __face_id_in_elt_1 = it->pos_second();
2699 
2700  __c1 = gmc_ptrtype( new gmc_type( gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
2701 
2702  map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
2703  fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
2704 
2705  expr2 = eval2_expr_ptrtype( new eval2_expr_type( expression(), mapgmc ) );
2706  expr2->init( im() );
2707  }
2708 
2709  //
2710  // start the real intensive job:
2711  // -# iterate over all elements to integrate over
2712  // -# construct the associated geometric mapping with the reference element
2713  // -# loop over quadrature loop and assemble the local matrix associated with the bilinear form
2714  // -# assemble the local contribution in the global representation of the bilinear form
2715  //
2716  for ( ; it != en; ++it )
2717  {
2718  if ( it->isGhostFace() )
2719  {
2720  LOG(WARNING) << "face id : " << it->id() << " is a ghost face";
2721  continue;
2722  }
2723  if ( it->isConnectedTo1() )
2724  {
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();
2729 
2730  __c0->update( it->element( 0 ), __face_id_in_elt_0 );
2731  __c1->update( it->element( 1 ), __face_id_in_elt_1 );
2732 
2733 #if 0
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";
2750 #endif
2751 
2752  __typeof__( im( __face_id_in_elt_0 ) ) im_face ( im( __face_id_in_elt_0 ) );
2753  //std::cout << "pts = " << im_face.points() << "\n";
2754  map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
2755  fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
2756 
2757  expr2->update( mapgmc, __face_id_in_elt_0 );
2758  const gmc_type& gmc = *__c0;
2759 
2760  __integrators[__face_id_in_elt_0].update( gmc );
2761 
2762  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2763  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2764  {
2765  res( c1,c2 ) += __integrators[__face_id_in_elt_0]( *expr2, c1, c2 );
2766  }
2767  }
2768 
2769  else
2770  {
2771  //LOG_IF( !it->isConnectedTo0(), WARN ) << "integration invalid boundary face";
2772  if ( !it->isConnectedTo0() || it->pos_first() == invalid_uint16_type_value )
2773  continue;
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 );
2778  //expr->update( mapgmc );
2779  const gmc_type& gmc = *__c0;
2780 
2781  __integrators[__face_id_in_elt_0].update( gmc );
2782 
2783  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2784  for ( uint16_type c2 = 0; c2 < eval::shape::N; ++c2 )
2785  {
2786  res( c1,c2 ) += __integrators[__face_id_in_elt_0]( *expr, c1, c2 );
2787  }
2788  } // !isConnectedTo1
2789  } // for loop on face
2790  }
2791  //std::cout << "res=" << res << "\n";
2792  //std::cout << "res1=" << res1 << "\n";
2793  DLOG(INFO) << "integrating over faces done in " << __timer.elapsed() << "s\n";
2794  return res;
2795 }
2796 
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
2800 {
2801  DLOG(INFO) << "integrating over "
2802  << std::distance( this->beginElement(), this->endElement() ) << " points\n";
2803 
2804  // first loop on the points, then retrieve the elements to which they belong
2805  // and evaluate the integrand expression and accumulate it
2806 
2807 
2808 }
2809 
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
2814 {
2815  DLOG(INFO) << "integrating over "
2816  << std::distance( this->beginElement(), this->endElement() ) << " elements\n";
2817  boost::timer __timer;
2818 
2819  //
2820  // some typedefs
2821  //
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;
2829  //typedef typename eval::gm_type gmc_type;
2830  //typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
2831 
2832  //typedef typename eval_expr_type::value_type value_type;
2833  //typedef typename Im::value_type value_type;
2834 
2835  auto p0 = P0h->element( "p0" );
2836  // set to 0 first
2837  p0.zero();
2838 
2839  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2840  {
2841  auto it = lit->template get<1>();
2842  auto en = lit->template get<2>();
2843 
2844 
2845  // make sure that we have elements to iterate over (return 0
2846  // otherwise)
2847  if ( it == en )
2848  continue;
2849  //return p0;
2850 
2851  //
2852  // Precompute some data in the reference element for
2853  // geometric mapping and reference finite element
2854  //
2855  gm_ptrtype gm = it->gm();
2856  //DDLOG(INFO) << "[integrator] evaluate(elements), gm is cached: " << gm->isCached() << "\n";
2857  typename eval::gmpc_ptrtype __geopc( new typename eval::gmpc_type( gm,
2858  this->im().points() ) );
2859 
2860 
2861  it = this->beginElement();
2862  // wait for all the guys
2863 #ifdef FEELPP_HAS_MPI
2864  auto const& worldComm = const_cast<MeshBase*>( it->mesh() )->worldComm();
2865 
2866 #if 0
2867  if ( worldComm.size() > 1 )
2868  {
2869  worldComm.barrier();
2870  }
2871 #endif
2872 #endif
2873 
2874 
2875 
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 ) );
2879 
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;
2883 
2884  //value_type res1 = 0;
2885  for ( ; it != en; ++it )
2886  {
2887  boost::timer ti;
2888  __c->update( *it );
2889  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
2890  expr.update( mapgmc );
2891  const gmc_type& gmc = *__c;
2892 
2893  M_im.update( gmc );
2894 
2895 
2896  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
2897  {
2898  size_type i= P0h->dof()->localToGlobal( it->id(), 0, c1 ).index();
2899  double v = M_im( expr, c1, 0 );
2900  p0.set( i, v );
2901  }
2902  }
2903  }
2904  //std::cout << "res=" << res << "\n";
2905  //std::cout << "res1=" << res1 << "\n";
2906  DLOG(INFO) << "integrating over elements done in " << __timer.elapsed() << "s\n";
2907 
2908  return p0;
2909 }
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
2914 {
2915  DLOG(INFO) << "integrating over "
2916  << std::distance( this->beginElement(), this->endElement() ) << "faces\n";
2917  boost::timer __timer;
2918 
2919  //
2920  // some typedefs
2921  //
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;
2929  //typedef typename eval_expr_type::value_type value_type;
2930  //typedef typename Im::value_type value_type;
2931 
2932  //BOOST_MPL_ASSERT_MSG( the_element_type::nDim > 1, INVALID_DIM, (mpl::int_<the_element_type::nDim>, mpl::int_<the_face_element_type::nDim>, mpl::identity<the_face_element_type>, mpl::identity<the_element_type> ) );;
2933 
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;
2937 
2938  //
2939  // Precompute some data in the reference element for
2940  // geometric mapping and reference finite element
2941  //
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;
2946 
2947  std::vector<im_face_type> __integrators;
2948 
2949  auto p0 = P0h->element( "p0" );
2950  // set to 0 first
2951  p0.zero();
2952 
2953  for( auto lit = M_elts.begin(), len = M_elts.end(); lit != len; ++lit )
2954  {
2955  auto it = lit->template get<1>();
2956  auto en = lit->template get<2>();
2957 
2958  // make sure that we have elements to iterate over (return 0
2959  // otherwise)
2960  if ( it == en )
2961  continue;
2962  //return p0;
2963 
2964  gm_ptrtype gm = it->element( 0 ).gm();
2965 
2966  //DDLOG(INFO) << "[integrator] evaluate(faces), gm is cached: " << gm->isCached() << "\n";
2967  for ( uint16_type __f = 0; __f < im().nFaces(); ++__f )
2968  {
2969  __integrators.push_back( im( __f ) );
2970 
2971  for ( permutation_type __p( permutation_type::IDENTITY );
2972  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
2973  {
2974  //FEELPP_ASSERT( ppts[__f][__p]->size2() != 0 ).warn( "invalid quadrature type" );
2975  __geopc[__f][__p] = pc_ptrtype( new pc_type( gm, ppts[__f].find( __p )->second ) );
2976  }
2977  }
2978 
2979 
2980  uint16_type __face_id_in_elt_0 = it->pos_first();
2981 
2982  // get the geometric mapping associated with element 0
2983  gmc_ptrtype __c0( new gmc_type( gm, it->element( 0 ), __geopc, __face_id_in_elt_0 ) );
2984 
2985  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
2986 
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 ) );
2991  expr->init( im() );
2992 
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;
2997 
2998  // true if connected to another element, false otherwise
2999  bool isConnectedTo1 = it->isConnectedTo1();
3000 
3001  // get the geometric mapping associated with element 1
3002  gmc_ptrtype __c1;
3003 
3004  //value_type res = 0;
3005  //value_type res1 = 0;
3006  if ( isConnectedTo1 )
3007  {
3008  uint16_type __face_id_in_elt_1 = it->pos_second();
3009 
3010  __c1 = gmc_ptrtype( new gmc_type( gm, it->element( 1 ), __geopc, __face_id_in_elt_1 ) );
3011 
3012  map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
3013  fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
3014 
3015  expr2 = eval2_expr_ptrtype( new eval2_expr_type( expression(), mapgmc ) );
3016  expr2->init( im() );
3017  }
3018 
3019  //
3020  // start the real intensive job:
3021  // -# iterate over all elements to integrate over
3022  // -# construct the associated geometric mapping with the reference element
3023  // -# loop over quadrature loop and assemble the local matrix associated with the bilinear form
3024  // -# assemble the local contribution in the global representation of the bilinear form
3025  //
3026  for ( ; it != en; ++it )
3027  {
3028 
3029  if ( it->isConnectedTo1() )
3030  {
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();
3035 
3036  __c0->update( it->element( 0 ), __face_id_in_elt_0 );
3037  __c1->update( it->element( 1 ), __face_id_in_elt_1 );
3038 
3039 #if 0
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";
3056 #endif
3057 
3058  __typeof__( im( __face_id_in_elt_0 ) ) im_face ( im( __face_id_in_elt_0 ) );
3059  //std::cout << "pts = " << im_face.points() << "\n";
3060  map2_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c0 ),
3061  fusion::make_pair<vf::detail::gmc<1> >( __c1 ) );
3062 
3063  expr2->update( mapgmc, __face_id_in_elt_0 );
3064  const gmc_type& gmc = *__c0;
3065 
3066  __integrators[__face_id_in_elt_0].update( gmc );
3067 
3068  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
3069  {
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 );
3073  p0.add( i0, v );
3074  p0.add( i1, v );
3075  }
3076  }
3077 
3078  else
3079  {
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 );
3084  //expr->update( mapgmc );
3085  const gmc_type& gmc = *__c0;
3086 
3087  __integrators[__face_id_in_elt_0].update( gmc );
3088 
3089  for ( uint16_type c1 = 0; c1 < eval::shape::M; ++c1 )
3090  {
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 );
3093  p0.add( i0, v );
3094  }
3095  } // !isConnectedTo1
3096  } // for loop on face
3097  }
3098  //std::cout << "res=" << res << "\n";
3099  //std::cout << "res1=" << res1 << "\n";
3100  DLOG(INFO) << "integrating over faces done in " << __timer.elapsed() << "s\n";
3101  return p0;
3102 }
3104 
3105 #if 0
3106 
3110 template<typename Elts, typename Im, typename ExprT>
3111 Expr<Integrator<Elts, Im, ExprT, Im> >
3112 integrate( Elts const& elts,
3113  Im const& im,
3114  ExprT const& expr,
3115  GeomapStrategyType gt = GeomapStrategyType::GEOMAP_HO )
3116 {
3117  typedef Integrator<Elts, Im, ExprT, Im> expr_t;
3118  return Expr<expr_t>( expr_t( elts, im, expr, gt, im ) );
3119 }
3120 #endif
3121 
3122 //Macro which get the good integration order
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 \
3129 
3130 
3131 #if 0
3132 
3136 template<typename Elts, typename ExprT>
3137 Expr<Integrator<Elts, _Q< ExpressionOrder<Elts,ExprT>::value >, ExprT, _Q< ExpressionOrder<Elts,ExprT>::value > > >
3138 integrate( Elts const& elts,
3139  ExprT const& expr,
3140  GeomapStrategyType gt = GeomapStrategyType::GEOMAP_HO )
3141 {
3142 
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 );
3146 
3147 }
3148 #endif
3149 
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,
3157  Im const& im,
3158  ExprT const& expr,
3159  GeomapStrategyType const& gt,
3160  Im2 const& im2,
3161  bool use_tbb,
3162  int grainsize,
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 > >() )
3166 {
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 ) );
3169 }
3170 
3171 
3173 namespace detail
3174 {
3175 template<typename TheArgs, typename Tag>
3176 struct clean_type
3177 {
3178  typedef typename boost::remove_pointer<
3179  typename boost::remove_const<
3180  typename boost::remove_reference<
3181  typename parameter::binding<TheArgs, Tag>::type
3182  >::type
3183  >::type
3184  >::type type;
3185 };
3186 template<typename TheArgs, typename Tag, typename Default>
3187 struct clean2_type
3188 {
3189  typedef typename boost::remove_pointer<
3190  typename boost::remove_const<
3191  typename boost::remove_reference<
3192  typename parameter::binding<TheArgs, Tag, Default>::type
3193  >::type
3194  >::type
3195  >::type type;
3196 };
3197 
3198 template<typename Args>
3199 struct integrate_type
3200 {
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;
3203  //typedef _Q< ExpressionOrder<_range_type,_expr_type>::value > the_quad_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;
3207 
3208  typedef boost::shared_ptr<QuadPtLocalization<_range_type,_quad_type,_expr_type > > _quadptloc_ptrtype;
3209 };
3210 }
3212 
3223 BOOST_PARAMETER_FUNCTION(
3224  ( typename vf::detail::integrate_type<Args>::expr_type ), // return type
3225  integrate, // 2. function name
3226 
3227  tag, // 3. namespace of tag types
3228 
3229  ( required
3230  ( range, * )
3231  ( expr, * )
3232  ) // 4. one required parameter, and
3233 
3234  ( optional
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() )
3243  )
3244 )
3245 {
3246 
3247  auto ret = integrate_impl( range, quad, expr, geomap, quad1, use_tbb, grainsize, partitioner, quadptloc );
3248 
3249  if ( verbose )
3250  {
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";
3255  //std::cout << " -- integrate: geomap = " << geomap << "\n";
3256  }
3257 
3258  return ret;
3259 }
3260 
3261 BOOST_PARAMETER_FUNCTION(
3262  ( double ), // return type
3263  normL2, // 2. function name
3264 
3265  tag, // 3. namespace of tag types
3266 
3267  ( required
3268  ( range, * )
3269  ( expr, * )
3270  ) // 4. one required parameter, and
3271 
3272  ( optional
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 )
3280  )
3281 )
3282 {
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 );
3287 }
3288 
3289 BOOST_PARAMETER_FUNCTION(
3290  ( double ), // return type
3291  normL2Squared, // 2. function name
3292 
3293  tag, // 3. namespace of tag types
3294 
3295  ( required
3296  ( range, * )
3297  ( expr, * )
3298  ) // 4. one required parameter, and
3299 
3300  ( optional
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 )
3308  )
3309 )
3310 {
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 );
3314 }
3315 
3316 BOOST_PARAMETER_FUNCTION(
3317  ( double ), // return type
3318  normH1, // 2. function name
3319 
3320  tag, // 3. namespace of tag types
3321 
3322  ( required
3323  ( range, * )
3324  ( expr, * )
3325  ( grad_expr, *)
3326  ) // 4. one required parameter, and
3327 
3328  ( optional
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 )
3336  )
3337 )
3338 {
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 );
3346 }
3347 
3348 BOOST_PARAMETER_FUNCTION(
3349  ( double ), // return type
3350  normSemiH1, // 2. function name
3351 
3352  tag, // 3. namespace of tag types
3353 
3354  ( required
3355  ( range, * )
3356  ( grad_expr, *)
3357  ) // 4. one required parameter, and
3358 
3359  ( optional
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 )
3367  )
3368 )
3369 {
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 );
3374 }
3375 
3376 BOOST_PARAMETER_FUNCTION(
3377  ( typename vf::detail::integrate_type<Args>::expr_type::expression_type::matrix_type ), // return type
3378  mean, // 2. function name
3379 
3380  tag, // 3. namespace of tag types
3381 
3382  ( required
3383  ( range, * )
3384  ( expr, * )
3385  ) // 4. one required parameter, and
3386 
3387  ( optional
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 )
3395  )
3396 )
3397 {
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";
3407  return eint/meas;
3408 }
3409 
3410 } // vf
3411 
3412 
3413 } // feel
3415 
3416 
3417 #endif /* __Integrator_H */
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

Generated on Sun Dec 22 2013 13:11:05 for Feel++ by doxygen 1.8.5