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
linearform.hpp
Go to the documentation of this file.
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-18
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2011 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 __LinearForm_H
31 #define __LinearForm_H 1
32 
33 #include <Eigen/Eigen>
34 #include <Eigen/StdVector>
35 
36 #include <boost/fusion/support/pair.hpp>
37 #include <boost/fusion/container.hpp>
38 #include <boost/fusion/sequence.hpp>
39 #include <boost/fusion/algorithm.hpp>
40 #include <boost/multi_array.hpp>
41 #include <feel/feelvf/block.hpp>
42 #include <feel/feelalg/vectorvalue.hpp>
43 #include <feel/feelvf/fec.hpp>
45 
46 namespace Feel
47 {
48 namespace fusion = boost::fusion;
49 namespace vf
50 {
52 namespace detail
53 {
73 template<typename SpaceType,typename VectorType,typename ElemContType>
74 class LinearForm
75 {
76 public:
77 
78 
82  enum { nDim = SpaceType::nDim };
83 
84  typedef LinearForm<SpaceType, VectorType, ElemContType> self_type;
85 
86  typedef SpaceType space_type;
87  typedef boost::shared_ptr<SpaceType> space_ptrtype;
88 
89  typedef typename space_type::value_type value_type;
90  typedef typename space_type::real_type real_type;
91  typedef VectorType vector_type;
92  typedef boost::shared_ptr<VectorType> vector_ptrtype;
93  typedef typename space_type::template Element<value_type, ElemContType> element_type;
94 
95 #if 0
96  typedef typename space_type::component_fespace_type component_fespace_type;
97  typedef typename space_type::element_type::component_type component_type;
98 #endif
99  typedef typename space_type::mesh_type mesh_type;
100  typedef typename mesh_type::element_type mesh_test_element_type;
101  typedef typename mesh_type::element_type::permutation_type permutation_type;
102  typedef typename space_type::gm_type gm_type;
103  typedef typename space_type::gm_ptrtype gm_ptrtype;
104  typedef typename space_type::gm1_type gm1_type;
105  typedef typename space_type::gm1_ptrtype gm1_ptrtype;
106 
107  typedef typename space_type::fe_type fe_type;
108  typedef typename space_type::basis_0_type::precompute_type test_precompute_type;
109 
110  typedef boost::shared_ptr<test_precompute_type> test_precompute_ptrtype;
112 
117 
130  template<typename GeomapContext,
131  typename ExprT,
132  typename IM,
133  typename GeomapExprContext = GeomapContext
134  >
135  class Context //: public FormContextBase<GeomapContext, IM>
136  {
137  typedef FormContextBase<GeomapContext, IM, GeomapExprContext> super;
138 
139  public:
140  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
141 
142  typedef Context<GeomapContext,ExprT,IM,GeomapExprContext> form_context_type;
143  typedef LinearForm<SpaceType,VectorType, ElemContType> form_type;
144  typedef typename SpaceType::dof_type dof_type;
145  typedef typename form_type::value_type value_type;
146 
147 
148 #if 0
149  typedef GeomapContext map_geometric_mapping_context_type;
150  typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type geometric_mapping_context_ptrtype;
151  typedef typename geometric_mapping_context_ptrtype::element_type geometric_mapping_context_type;
152  typedef typename geometric_mapping_context_type::gm_type geometric_mapping_type;
153 #else
154  typedef typename super::map_geometric_mapping_context_type map_test_geometric_mapping_context_type;
155  typedef typename super::geometric_mapping_context_ptrtype test_geometric_mapping_context_ptrtype;
156  typedef typename super::geometric_mapping_context_type test_geometric_mapping_context_type;
157  typedef typename super::geometric_mapping_type test_geometric_mapping_type;
158 
159  typedef typename super::map_geometric_mapping_context_type map_trial_geometric_mapping_context_type;
160  typedef typename super::geometric_mapping_context_ptrtype trial_geometric_mapping_context_ptrtype;
161  typedef typename super::geometric_mapping_context_type trial_geometric_mapping_context_type;
162  typedef typename super::geometric_mapping_type trial_geometric_mapping_type;
163 
164 #endif
165 
166  static const uint16_type nDim = test_geometric_mapping_type::nDim;
167 
168  typedef ExprT expression_type;
169 
170  typedef typename space_type::mesh_type mesh_type;
171  typedef typename mesh_type::element_type mesh_element_type;
172  typedef typename mesh_element_type::permutation_type permutation_type;
173  typedef typename space_type::fe_type test_fe_type;
174  typedef typename space_type::fe_type trial_fe_type;
175  typedef boost::shared_ptr<test_fe_type> test_fe_ptrtype;
176  typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
177  test_fe_type,
178  test_geometric_mapping_type,
179  mesh_element_type> test_fecontext_type;
180  typedef test_fecontext_type trial_fecontext_type;
181 
182 #if 0
183  typedef mpl::int_<fusion::result_of::template size<GeomapContext>::type::value> map_size;
184 #else
185  typedef typename super::map_size map_size;
186 #endif
187 
188  typedef boost::shared_ptr<test_fecontext_type> test_fecontext_ptrtype;
189 
190 #if 0
191  typedef typename mpl::if_<mpl::equal_to<map_size,mpl::int_<2> >,
192  gmc<1>,
193  gmc<0> >::type gmc1;
194 #else
195  typedef typename super::gmc1 gmc1;
196 #endif
197 
198 #if 0
199  typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type left_gmc_ptrtype;
200  typedef typename fusion::result_of::value_at_key<GeomapContext,gmc<0> >::type::element_type left_gmc_type;
201  typedef typename fusion::result_of::value_at_key<GeomapContext,gmc1 >::type right_gmc_ptrtype;
202  typedef typename fusion::result_of::value_at_key<GeomapContext,gmc1 >::type::element_type right_gmc_type;
203 
204  typedef fusion::map<fusion::pair<gmc<0>, left_gmc_ptrtype> > map_left_gmc_type;
205  typedef fusion::map<fusion::pair<gmc<0>, right_gmc_ptrtype> > map_right_gmc_type;
206 #else
207  typedef typename super::left_gmc_ptrtype left_gmc_ptrtype;
208  typedef typename super::left_gmc_type left_gmc_type;
209  typedef typename super::right_gmc_ptrtype right_gmc_ptrtype;
210  typedef typename super::right_gmc_type right_gmc_type;
211 
212  typedef typename super::map_left_gmc_type map_left_gmc_type;
213  typedef typename super::map_right_gmc_type map_right_gmc_type;
214 #endif
215 
216  typedef typename mpl::if_<mpl::equal_to<map_size, mpl::int_<1> >,
217  mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > >,
218  mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype>,
219  fusion::pair<gmc<1>, test_fecontext_ptrtype> > > >::type::type map_test_fecontext_type;
220 
221  typedef fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > map_left_test_fecontext_type;
222  typedef fusion::map<fusion::pair<gmc<1>, test_fecontext_ptrtype> > map_right_test_fecontext_type;
223 
224  typedef typename super::map_geometric_mapping_expr_context_type map_geometric_mapping_expr_context_type;
225 
226  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type, map_left_test_fecontext_type> eval0_expr_type;
227  typedef boost::shared_ptr<eval0_expr_type> eval0_expr_ptrtype;
228 
229  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type, map_right_test_fecontext_type> eval1_expr_type;
230  typedef boost::shared_ptr<eval1_expr_type> eval1_expr_ptrtype;
231  //typedef typename ExprT::template tensor<map_right_gmc_type, map_test_fecontext_type> eval1_expr_type;
232 
233 
234  typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
235  test_fe_type,
236  test_geometric_mapping_type,
237  mesh_element_type>::template Index<> test_index_type;
238 
239 
240  //typedef typename ExprT::template tensor<map_geometric_mapping_context_type, map0_test_fecontext_type> eval0_expr_type;
241  //typedef typename ExprT::template tensor<map_geometric_mapping_context_type, map1_test_fecontext_type> eval1_expr_type;
242 
243  //typedef ublas::vector<value_type> local_vector_type;
244  static const int rep_shape = 2;//1+(eval_expr_type::shape::M-1>0)+(eval_expr_type::shape::N-1>0);
245  //typedef boost::multi_array<value_type, rep_shape> local_vector_type;
246  typedef typename space_type::dof_type test_dof_type;
247  static const int nDofPerElementTest = space_type::dof_type::nDofPerElement;
248  typedef Eigen::Matrix<value_type, nDofPerElementTest, 1> local_vector_type;
249  typedef Eigen::Matrix<value_type, 2*nDofPerElementTest, 1> local2_vector_type;
250  typedef Eigen::Matrix<int, nDofPerElementTest, 1> local_row_type;
251  typedef Eigen::Matrix<int, 2*nDofPerElementTest, 1> local2_row_type;
252 
253  public:
254 
255  Context( form_type& __form,
256  map_test_geometric_mapping_context_type const& _gmcTest,
257  map_trial_geometric_mapping_context_type const & _gmcTrial, //useless(fix compilation)
258  map_geometric_mapping_expr_context_type const& _gmcExpr,
259  ExprT const& expr,
260  IM const& im );
261 
262  template<typename IM2>
263  Context( form_type& __form,
264  map_test_geometric_mapping_context_type const& _gmcTest,
265  map_trial_geometric_mapping_context_type const & _gmcTrial, //useless(fix compilation)
266  map_geometric_mapping_expr_context_type const& _gmcExpr,
267  ExprT const& expr,
268  IM const& im,
269  IM2 const& im2 );
270 
271  template<typename IM2>
272  Context( form_type& __form,
273  map_test_geometric_mapping_context_type const& _gmcTest,
274  map_trial_geometric_mapping_context_type const & _gmcTrial, //useless(fix compilation)
275  map_geometric_mapping_expr_context_type const& _gmcExpr,
276  ExprT const& expr,
277  IM const& im,
278  IM2 const& im2,
279  mpl::int_<2> );
280 
281  bool isZero( size_type i ) const
282  {
283  return false;
284  }
285  bool isZero( typename mesh_type::element_iterator it ) const
286  {
287  return this->isZero( it->id() );
288  }
289  bool isZero( typename mesh_type::element_type const& e ) const
290  {
291  return this->isZero( e.id() );
292  }
293 
294  void update( map_test_geometric_mapping_context_type const& _gmcTest,
295  map_trial_geometric_mapping_context_type const & gmcTrial,
296  map_geometric_mapping_expr_context_type const& _gmcExpr );
297 
298 
299  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
300  map_trial_geometric_mapping_context_type const & gmcTrial,
301  map_geometric_mapping_expr_context_type const& gmcExpr,
302  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad );
303 
304  void update( map_test_geometric_mapping_context_type const& _gmcTest,
305  map_trial_geometric_mapping_context_type const & gmcTrial,
306  map_geometric_mapping_expr_context_type const& _gmcExpr,
307  mpl::int_<2> );
308 
309 
310  void update( map_test_geometric_mapping_context_type const& _gmcTest,
311  map_trial_geometric_mapping_context_type const & gmcTrial,
312  map_geometric_mapping_expr_context_type const& _gmcExpr,
313  IM const& im );
314 
315  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
316  map_trial_geometric_mapping_context_type const & gmcTrial,
317  map_geometric_mapping_expr_context_type const& gmcExpr,
318  IM const& im,
319  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad )
320  {
321  M_integrator = im;
322  this->updateInCaseOfInterpolate( gmcTest, gmcTrial, gmcExpr, indexLocalToQuad );
323  }
324 
325 
326  void update( map_test_geometric_mapping_context_type const& _gmcTest,
327  map_trial_geometric_mapping_context_type const & gmcTrial,
328  map_geometric_mapping_expr_context_type const& _gmcExpr,
329  IM const& im, mpl::int_<2> );
330 
331 
332  void integrate()
333  {
334  integrate( mpl::int_<fusion::result_of::size<GeomapContext>::type::value>() );
335  }
336 
337  void integrateInCaseOfInterpolate( std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
338  bool isFirstExperience )
339  {
340  integrateInCaseOfInterpolate( mpl::int_<fusion::result_of::size<GeomapContext>::type::value>(),
341  indexLocalToQuad, isFirstExperience );
342  }
343 
344 
345  void assemble()
346  {
347  assemble( M_gmc_left->id() );
348  }
349  void assemble( size_type elt_0 );
350 
351  void assemble( mpl::int_<2> )
352  {
353  assemble( M_gmc_left->id(), M_gmc_right->id() );
354  }
355  void assemble( size_type elt_0, size_type elt_1 );
356 
361  template<typename Pts>
362  void precomputeBasisAtPoints( Pts const& pts )
363  {
364  M_test_pc = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), pts ) );
365  }
366 
372  template<typename Pts>
373  void precomputeBasisAtPoints( uint16_type __f, permutation_type const& __p, Pts const& pts )
374  {
375  M_test_pc_face[__f][__p] = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), pts ) );
376  //FEELPP_ASSERT( M_test_pc_face.find(__f )->second )( __f ).error( "invalid test precompute type" );
377  }
385  test_precompute_ptrtype const& testPc( uint16_type __f,
386  permutation_type __p = permutation_type( permutation_type::NO_PERMUTATION ) ) const
387  {
388  if ( __f == invalid_uint16_type_value )
389  return M_test_pc;
390 
391  //FEELPP_ASSERT( M_test_pc_face.find(__f )->second )( __f ).error( "invalid test precompute type" );
392  return M_test_pc_face.find( __f )->second.find( __p )->second;
393  }
394 
395  template<typename PtsSet>
396  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
397  precomputeTestBasisAtPoints( PtsSet const& pts )
398  {
399  typedef typename boost::is_same< permutation_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
400  return precomputeTestBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
401  }
402 
403 
404  template<typename PtsSet>
405  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
406  precomputeTestBasisAtPoints( PtsSet const& pts, mpl::bool_<false> )
407  {
408  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
409  return testpc;
410  }
411 
412  template<typename PtsSet>
413  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> >
414  precomputeTestBasisAtPoints( PtsSet const& pts, mpl::bool_<true> )
415  {
416  QuadMapped<PtsSet> qm;
417  typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
418  typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
419 
420  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
421 
422  for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
423  {
424  for ( permutation_type __p( permutation_type::IDENTITY );
425  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
426  {
427  testpc[__f][__p] = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), ppts[__f].find( __p )->second ) );
428  }
429  }
430 
431  return testpc;
432  }
433 
434  private:
435 
436 
437  void integrate( mpl::int_<1> );
438 
439  void integrate( mpl::int_<2> );
440 
441  void integrateInCaseOfInterpolate( mpl::int_<1>,
442  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
443  bool isFirstExperience );
444 
445  private:
446 
447  form_type& M_form;
448  dof_type* M_test_dof;
449  const list_block_type& M_lb;
450 
451  test_precompute_ptrtype M_test_pc;
452  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > M_test_pc_face;
453 
454  map_test_geometric_mapping_context_type M_gmc;
455  left_gmc_ptrtype M_gmc_left;
456  right_gmc_ptrtype M_gmc_right;
457  map_left_gmc_type M_left_map;
458  map_right_gmc_type M_right_map;
459  map_test_fecontext_type M_test_fec;
460  map_left_test_fecontext_type M_test_fec0;
461  map_right_test_fecontext_type M_test_fec1;
462 
463  local_vector_type M_rep;
464  local2_vector_type M_rep_2;
465  local_row_type M_local_rows;
466  local2_row_type M_local_rows_2;
467  local_row_type M_local_rowsigns;
468  local2_row_type M_local_rowsigns_2;
469 
470  eval0_expr_ptrtype M_eval0_expr;
471  eval1_expr_ptrtype M_eval1_expr;
472 
473  IM M_integrator;
474 
475  }; // Context
476 
477 
479 
483 
484  LinearForm( LinearForm const & __vf );
485 
486  LinearForm( space_ptrtype const& __X,
487  vector_ptrtype __F,
488  size_type rowstart = 0,
489  bool init = true,
490  bool do_threshold = false,
491  value_type threshold = type_traits<value_type>::epsilon() );
492  LinearForm( space_ptrtype const& __X,
493  vector_ptrtype __F,
494  list_block_type const& __lb,
495  size_type rowstart = 0,
496  bool init = true,
497  bool do_threshold = false,
498  value_type threshold = type_traits<value_type>::epsilon() );
499 
500  ~LinearForm()
501  {}
502 
503 
504 
506 
510 
517  template <class ExprT>
518  LinearForm& operator=( Expr<ExprT> const& expr );
519 
527  template <class ExprT>
528  LinearForm& operator+=( Expr<ExprT> const& expr );
529 
530 #if 0
531 
535  LinearForm<component_fespace_type,
536  vector_type,
537  typename component_type::container_type>
538  operator()( component_type& __u )
539  {
540  typedef typename component_type::container_type container_type;
541 
542  typedef LinearForm<component_fespace_type, vector_type,
543  container_type> vflf_type;
544 
545  list_block_type __list_block;
546  __list_block.push_back( Block( 0, 0,
547  __u.start(), 0 ) );
548  return vflf_type( __u, M_F, __list_block, false );
549 
550  }
551 #endif
552 
556  value_type operator()( size_type i ) const
557  {
558  return M_F( i );
559  }
560 
567  value_type operator()( element_type const& __v ) const
568  {
569  return M_F->dot( __v );
570  }
571 
578  value_type operator()( typename space_type::element_type const& __v ) const
579  {
580  return M_F->dot( __v );
581  }
582 
583 
585 
589 
593  space_ptrtype const& functionSpace() const
594  {
595  return M_X;
596  }
597 
601  space_ptrtype const& testSpace() const
602  {
603  return M_X;
604  }
605 
609  gm_ptrtype const& gm() const
610  {
611  return M_X->gm();
612  }
613 
617  gm1_ptrtype const& gm1() const
618  {
619  return M_X->gm1();
620  }
621 
628  test_precompute_ptrtype const& testPc( uint16_type __f,
629  permutation_type __p = permutation_type( permutation_type::NO_PERMUTATION ) ) const
630  {
631  if ( __f == invalid_uint16_type_value )
632  return M_test_pc;
633 
634  return M_test_pc_face.find( __f )->second.find( __p )->second;
635  }
636 
637  vector_type& representation() const
638  {
639  return *M_F;
640  }
641  vector_ptrtype vectorPtr() const
642  {
643  return M_F;
644  }
645  vector_ptrtype vectorPtr()
646  {
647  return M_F;
648  }
649 
650  list_block_type const& blockList() const
651  {
652  return M_lb;
653  }
654 
655  size_type rowStartInVector() const
656  {
657  return M_row_startInVector;
658  }
659 
660 
664  bool doThreshold( value_type const& v ) const
665  {
666  return ( math::abs( v ) > M_threshold );
667  }
668 
670 
674 
679  void setFunctionSpace( space_ptrtype const& __X )
680  {
681  M_X = __X;
682  }
683 
688  void setThreshold( value_type eps )
689  {
690  M_threshold = eps;
691  }
692 
697  void setDoThreshold( bool do_threshold )
698  {
699  M_do_threshold = do_threshold;
700  }
701 
703 
707 
712  template<typename Pts>
713  void precomputeBasisAtPoints( Pts const& pts )
714  {
715  M_test_pc = test_precompute_ptrtype( new test_precompute_type( functionSpace()->fe(), pts ) );
716  }
717 
723  template<typename Pts>
724  void precomputeBasisAtPoints( uint16_type __f, permutation_type __p, Pts const& pts )
725  {
726  M_test_pc_face[__f][__p] = test_precompute_ptrtype( new test_precompute_type( functionSpace()->fe(), pts ) );
727  }
728 
733  void add( size_type i, value_type const& v )
734  {
735  if ( M_do_threshold )
736  {
737  if ( doThreshold( v ) )
738  M_F->add( i+this->rowStartInVector(), v );
739  }
740 
741  else
742  M_F->add( i+this->rowStartInVector(), v );
743  }
744 
749  void addVector( int* i, int n, value_type* v )
750  {
751  if ( this->rowStartInVector()!=0 )
752  for ( int k = 0; k< n ; ++k )
753  i[k]+=this->rowStartInVector();
754 
755  M_F->addVector( i, n, v );
756  }
757 
762  void set( size_type i, value_type const& v )
763  {
764  M_F->set( i, v );
765  }
766 
770  void zero() { M_F->zero(); }
771 
772  LinearForm& operator+=( LinearForm& f )
773  {
774  if ( this == &f )
775  return *this;
776 
777  *M_F += *f.M_F;
778 
779  return *this;
780  }
782 
783 protected:
784 
785  LinearForm();
786 
787 private:
788  template <class ExprT> void assign( Expr<ExprT> const& expr, bool init, mpl::bool_<false> );
789  template <class ExprT> void assign( Expr<ExprT> const& expr, bool init, mpl::bool_<true> );
790 private:
791 
792  space_ptrtype M_X;
793 
794  vector_ptrtype M_F;
795 
796  list_block_type M_lb;
797 
798  size_type M_row_startInVector;
799 
800  test_precompute_ptrtype M_test_pc;
801 
802  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > M_test_pc_face;
803 
804  bool M_do_threshold;
805  value_type M_threshold;
806 };
807 
808 template<typename SpaceType, typename VectorType, typename ElemContType>
809 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( LinearForm const & __vf )
810  :
811  M_X( __vf.M_X ),
812  M_F( __vf.M_F ),
813  M_lb( __vf.M_lb ),
814  M_row_startInVector( __vf.M_row_startInVector ),
815  M_test_pc(),
816  M_test_pc_face(),
817  M_do_threshold( __vf.M_do_threshold ),
818  M_threshold( __vf.M_threshold )
819 
820 {
821  DVLOG(2) << "LinearForm copy constructor\n";
822  DVLOG(2) << " n Dof : " << M_X->nDof() << "\n";
823  DVLOG(2) << " F size : " << M_F->size() << "\n";
824  DVLOG(2) << "block size : " << M_lb.size() << "\n";
825 }
826 
827 template<typename SpaceType, typename VectorType, typename ElemContType>
828 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( space_ptrtype const& __X,
829  vector_ptrtype __F,
830  size_type rowstart,
831  bool init,
832  bool do_threshold,
833  value_type threshold )
834  :
835  M_X( __X ),
836  M_F( __F ),
837  M_lb(),
838  M_row_startInVector( rowstart ),
839  M_do_threshold( do_threshold ),
840  M_threshold( threshold )
841 {
842 
843  for ( uint16_type __i = 0; __i < M_X->qDim(); ++__i )
844  {
845  M_lb.push_back( Block( __i, 0,
846  __i*M_X->nDofPerComponent(),
847  0 ) );
848  DVLOG(2) << "[linearform::linearform] block: "
849  << Block( __i, 0, __i*M_X->nDofPerComponent(), 0 ) << "\n";
850  }
851 
852  if ( init )
853  M_F->zero();
854 }
855 
856 template<typename SpaceType, typename VectorType, typename ElemContType>
857 LinearForm<SpaceType, VectorType, ElemContType>::LinearForm( space_ptrtype const& __X,
858  vector_ptrtype __F,
859  list_block_type const& __lb,
860  size_type rowstart,
861  bool init,
862  bool do_threshold,
863  value_type threshold )
864  :
865  M_X( __X ),
866  M_F( __F ),
867  M_lb( __lb ),
868  M_row_startInVector( rowstart ),
869  M_do_threshold( do_threshold ),
870  M_threshold( threshold )
871 {
872  if ( init )
873  M_F->zero();
874 }
875 
876 template<typename LFType, typename TheSpaceType, typename ExprType>
877 struct LFAssign
878 {
879  LFAssign( LFAssign const& lfa )
880  :
881  M_lf( lfa.M_lf ),
882  M_Xh( lfa.M_Xh ),
883  M_expr( lfa.M_expr ),
884  M_index( lfa.M_index ),
885  M_init( lfa.M_init )
886  {}
887  LFAssign( LFType& lf, TheSpaceType const& Xh, ExprType const& expr, bool init )
888  :
889  M_lf( lf ),
890  M_Xh( Xh ),
891  M_expr( expr ),
892  M_index( 0 ),
893  M_init( init )
894  {}
895  template<typename SpaceType>
896  void operator()( boost::shared_ptr<SpaceType> const& X ) const
897  {
898  if ( M_lf.testSpace()->worldsComm()[M_index].isActive() )
899  {
900  DVLOG(2) << "expression has test functions ? :"
901  << ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result
902  << "\n";
903 
904  if ( !ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result )
905  {
906  ++M_index;
907  return;
908  }
909 
910  list_block_type __list_block;
911 
912  if ( M_lf.testSpace()->worldsComm()[M_index].globalSize()>1 )
913  {
914  if (M_lf.testSpace()->hasEntriesForAllSpaces())
915  __list_block.push_back( Block( 0, 0, M_Xh->nLocalDofStart( M_index ), 0 ) );
916  else
917  __list_block.push_back( Block( 0, 0, 0, 0 ) );
918  }
919  else
920  __list_block.push_back( Block( 0, 0, M_Xh->nDofStart( M_index ), 0 ) );
921 
922  LinearForm<SpaceType,typename LFType::vector_type, typename LFType::element_type> lf( X,
923  M_lf.vectorPtr(),
924  __list_block,
925  M_lf.rowStartInVector(),
926  false );
927 
928  //
929  // in composite integration, make sure that if M_init is \p
930  // true for the first space, it is set to \p false for the
931  // next spaces otherwise it will erase/clear to 0 the previous
932  // assemblies
933  //
934  if ( M_init )
935  {
936  // assembly
937  lf = M_expr;
938 
939  // make sure we won't erase the assembly we just did
940  M_init = false;
941  }
942 
943  else
944  {
945  lf += M_expr;
946  }
947  }
948 
949  else // not active : there is the init case with a close in zero
950  {
951  if ( M_init ) M_lf.representation().zero();
952  }
953 
954  ++M_index;
955  }
956 private:
957  LFType& M_lf;
958  TheSpaceType const& M_Xh;
959  ExprType const& M_expr;
960  mutable size_type M_index;
961  mutable bool M_init;
962 
963 };
964 
965 // implementation
966 template<typename LFType, typename TheSpaceType, typename ExprType>
967 LFAssign<LFType,TheSpaceType,ExprType>
968 make_lfassign( LFType& lf, TheSpaceType const& Xh, ExprType const& expr, bool init )
969 {
970  return LFAssign<LFType,TheSpaceType,ExprType>( lf, Xh, expr, init );
971 }
972 
973 
974 
975 template<typename SpaceType, typename VectorType, typename ElemContType>
976 template<typename ExprT>
977 void
978 LinearForm<SpaceType, VectorType, ElemContType>::assign( Expr<ExprT> const& __expr, bool init, mpl::bool_<true> )
979 {
980  fusion::for_each( M_X->functionSpaces(), make_lfassign( *this, M_X, __expr, init ) );
981 }
982 template<typename SpaceType, typename VectorType, typename ElemContType>
983 template<typename ExprT>
984 void
985 LinearForm<SpaceType, VectorType, ElemContType>::assign( Expr<ExprT> const& __expr, bool init, mpl::bool_<false> )
986 {
987  if ( M_lb.empty() )
988  {
989  M_lb.push_back( Block( 0, 0, 0, 0 ) );
990  }
991 
992  if ( init )
993  {
994  typename list_block_type::const_iterator __bit = M_lb.begin();
995  typename list_block_type::const_iterator __ben = M_lb.end();
996 
997  for ( ; __bit != __ben; ++__bit )
998  {
999  DVLOG(2) << "LinearForm:: block: " << *__bit << "\n";
1000  size_type g_ic_start = __bit->globalRowStart();
1001  DVLOG(2) << "LinearForm:: g_ic_start: " << g_ic_start << "\n";
1002 
1003  M_F->zero( g_ic_start,g_ic_start + M_X->nDof() );
1004  }
1005  }
1006 
1007  __expr.assemble( M_X, *this );
1008  //vector_range_type r( M_F, ublas::range( M_v.start(),M_v.start()+ M_v.size() ) );
1009  //std::cout << "r = " << r << "\n";
1010  //std::cout << "after F= " << M_F << "\n";
1011 }
1012 template<typename SpaceType, typename VectorType, typename ElemContType>
1013 template<typename ExprT>
1014 LinearForm<SpaceType, VectorType, ElemContType>&
1016 {
1017  // loop(fusion::for_each) over sub-functionspaces in SpaceType
1018  // pass expression and initialize
1019  this->assign( __expr, true, mpl::bool_<( SpaceType::nSpaces > 1 )>() );
1020  return *this;
1021 }
1022 template<typename SpaceType, typename VectorType, typename ElemContType>
1023 template<typename ExprT>
1024 LinearForm<SpaceType, VectorType, ElemContType>&
1025 LinearForm<SpaceType, VectorType, ElemContType>::operator+=( Expr<ExprT> const& __expr )
1026 {
1027  this->assign( __expr, false, mpl::bool_<( SpaceType::nSpaces > 1 )>() );
1028  return *this;
1029 }
1030 
1031 
1032 } // detail
1034 } // vf
1035 } // feel
1036 
1038 #endif /* __LinearForm_H */

Generated on Sun Oct 20 2013 08:25:01 for Feel++ by doxygen 1.8.4