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
bilinearform.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-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 __BilinearForm_H
31 #define __BilinearForm_H 1
32 
33 #include <Eigen/Eigen>
34 #include <Eigen/StdVector>
35 
36 
37 #include <set>
38 
39 #include <boost/parameter.hpp>
40 #include <boost/fusion/support/pair.hpp>
41 #include <boost/fusion/container.hpp>
42 #include <boost/fusion/sequence.hpp>
43 #include <boost/fusion/algorithm.hpp>
44 #include <boost/spirit/home/phoenix.hpp>
45 #include <boost/spirit/home/phoenix/core/argument.hpp>
46 #include <feel/feelcore/context.hpp>
47 #include <feel/feelalg/matrixvalue.hpp>
49 #include <feel/feelalg/backend.hpp>
50 #include <feel/feeldiscr/stencil.hpp>
51 #include <feel/feelvf/block.hpp>
52 #include <feel/feelvf/fec.hpp>
54 
55 
56 namespace Feel
57 {
58 namespace parameter = boost::parameter;
59 namespace fusion = boost::fusion;
60 namespace vf
61 {
62 
64 template<typename FE1,typename FE2,typename ElemContType> class BilinearForm;
65 namespace detail
66 {
67 
68 
69 template<typename BFType, typename ExprType>
70 struct BFAssign2
71 {
72  typedef typename BFType::matrix_type matrix_type;
73  BFAssign2( BFAssign2 const& lfa )
74  :
75  M_bf( lfa.M_bf ),
76  M_expr( lfa.M_expr ),
77  M_test_index( lfa.M_test_index )
78  {}
79  BFAssign2( BFType& lf, ExprType const& expr )
80  :
81  M_bf( lf ),
82  M_expr( expr ),
83  M_test_index( 0 )
84  {}
85  template<typename SpaceType>
86  void operator()( boost::shared_ptr<SpaceType> const& X ) const
87  {
88  DVLOG(2) << "[BFAssign2::operator()] start loop on trial spaces against test space index: " << M_test_index << "\n";
89 
90  if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
91  {
92  fusion::for_each( M_bf.trialSpace()->functionSpaces(),
93  make_bfassign1( M_bf, M_expr, X, M_test_index ) );
94  }
95 
96  DVLOG(2) << "[BFAssign2::operator()] stop loop on trial spaces against test space index: " << M_test_index << "\n";
97  ++M_test_index;
98 
99  }
100 private:
101  BFType& M_bf;
102  ExprType const& M_expr;
103  mutable size_type M_test_index;
104 };
105 template<typename BFType, typename ExprType>
106 BFAssign2<BFType,ExprType>
107 make_bfassign2( BFType& lf, ExprType const& expr )
108 {
109  return BFAssign2<BFType,ExprType>( lf, expr );
110 }
111 
112 template<typename BFType, typename ExprType, typename TestSpaceType>
113 struct BFAssign1
114 {
115  typedef typename BFType::matrix_type matrix_type;
116  BFAssign1( BFAssign1 const& lfa )
117  :
118  M_bf( lfa.M_bf ),
119  M_test( lfa.M_test ),
120  M_expr( lfa.M_expr ),
121  M_trial_index( lfa.M_trial_index )
122  {}
123  BFAssign1( BFType& lf,
124  ExprType const& expr,
125  boost::shared_ptr<TestSpaceType> const& Testh,
126  size_type test_index )
127  :
128  M_bf( lf ),
129  M_test( Testh ),
130  M_expr( expr ),
131  M_trial_index( 0 ),
132  M_test_index( test_index )
133  {}
134  template<typename SpaceType>
135  void operator()( boost::shared_ptr<SpaceType> const& trial ) const;
136 
137 private:
138  BFType& M_bf;
139  boost::shared_ptr<TestSpaceType> M_test;
140  ExprType const& M_expr;
141  mutable size_type M_trial_index;
142  size_type M_test_index;
143 };
144 template<typename BFType, typename ExprType, typename TestSpaceType>
145 BFAssign1<BFType,ExprType,TestSpaceType>
146 make_bfassign1( BFType& lf,
147  ExprType const& expr,
148  boost::shared_ptr<TestSpaceType> const& test_space,
149  size_type test_index )
150 {
151  return BFAssign1<BFType,ExprType,TestSpaceType>( lf, expr, test_space, test_index );
152 }
153 
154 
155 template<typename BFType, typename ExprType, typename TrialSpaceType>
156 struct BFAssign3
157 {
158  typedef typename BFType::matrix_type matrix_type;
159  BFAssign3( BFAssign3 const& lfa )
160  :
161  M_bf( lfa.M_bf ),
162  M_trial( lfa.M_trial ),
163  M_expr( lfa.M_expr ),
164  M_test_index( lfa.M_test_index )
165  {}
166  BFAssign3( BFType& lf,
167  ExprType const& expr,
168  boost::shared_ptr<TrialSpaceType> const& Trialh,
169  size_type trial_index )
170  :
171  M_bf( lf ),
172  M_trial( Trialh ),
173  M_expr( expr ),
174  M_trial_index( trial_index ),
175  M_test_index( 0 )
176  {}
177  template<typename SpaceType>
178  void operator()( boost::shared_ptr<SpaceType> const& test ) const;
179 
180 private:
181  BFType& M_bf;
182  boost::shared_ptr<TrialSpaceType> M_trial;
183  ExprType const& M_expr;
184  size_type M_trial_index;
185  mutable size_type M_test_index;
186 };
187 template<typename BFType, typename ExprType, typename TrialSpaceType>
188 BFAssign3<BFType,ExprType,TrialSpaceType>
189 make_bfassign3( BFType& lf,
190  ExprType const& expr,
191  boost::shared_ptr<TrialSpaceType> const& trial_space,
192  size_type trial_index )
193 {
194  return BFAssign3<BFType,ExprType,TrialSpaceType>( lf, expr, trial_space, trial_index );
195 }
196 
197 
205 template<typename FE1,
206  typename FE2,
207  typename ElemContType = VectorUblas<typename FE1::value_type> >
208 class BilinearForm
209 {
210 public:
211 
215  enum { nDim = FE1::nDim };
216 
217  typedef FE1 space_1_type;
218  typedef boost::shared_ptr<FE1> space_1_ptrtype;
219  typedef space_1_type test_space_type;
220  typedef boost::shared_ptr<space_1_type> test_space_ptrtype;
221 
222  typedef FE2 space_2_type;
223  typedef boost::shared_ptr<FE2> space_2_ptrtype;
224  typedef space_2_type trial_space_type;
225  typedef boost::shared_ptr<space_2_type> trial_space_ptrtype;
226 
227  typedef typename FE1::value_type value_type;
228  typedef typename space_1_type::template Element<value_type,ElemContType> element_1_type;
229 
230  typedef typename space_2_type::template Element<value_type,ElemContType> element_2_type;
231  typedef BilinearForm<FE1, FE2, ElemContType> self_type;
232 
233 #if 0
234  typedef typename space_1_type::component_fespace_type component_space_1_type;
235  typedef typename element_1_type::component_type component_1_type;
236  typedef typename space_2_type::component_fespace_type component_space_2_type;
237  typedef typename element_2_type::component_type component_2_type;
238 #endif // 0
239 
240  typedef typename space_1_type::mesh_type mesh_1_type;
241  typedef typename mesh_1_type::element_type mesh_element_1_type;
242  typedef typename mesh_element_1_type::permutation_type permutation_1_type;
243 
244  typedef typename space_2_type::mesh_type mesh_2_type;
245  typedef typename mesh_2_type::element_type mesh_element_2_type;
246  typedef typename mesh_element_2_type::permutation_type permutation_2_type;
247 
248  typedef typename space_1_type::fe_type fe_1_type;
249  typedef typename space_2_type::fe_type fe_2_type;
250 
251  typedef typename space_1_type::gm_type gm_1_type;
252  typedef typename space_1_type::gm_ptrtype gm_1_ptrtype;
253 
254  typedef typename space_1_type::gm1_type gm1_1_type;
255  typedef typename space_1_type::gm1_ptrtype gm1_1_ptrtype;
256 
257  typedef typename space_2_type::gm_type gm_2_type;
258  typedef typename space_2_type::gm_ptrtype gm_2_ptrtype;
259 
260  typedef typename space_2_type::gm1_type gm1_2_type;
261  typedef typename space_2_type::gm1_ptrtype gm1_2_ptrtype;
262 
263  //typedef ublas::compressed_matrix<value_type, ublas::row_major> csr_matrix_type;
264  typedef MatrixSparse<value_type> matrix_type;
265  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
266  static const bool is_row_major = true;//matrix_type::is_row_major;
267 
268  typedef typename mpl::if_<mpl::equal_to<mpl::bool_<is_row_major>, mpl::bool_<true> >,
269  mpl::identity<ublas::row_major>,
270  mpl::identity<ublas::column_major> >::type::type layout_type;
271 
272  template<int _N = 0>
273  struct test_precompute
274  {
275  //typedef typename space_1_type::basis_0_type::template precompute<_N>::type type;
276  typedef typename space_1_type::basis_0_type::PreCompute type;
277  typedef boost::shared_ptr<type> ptrtype;
278  };
279  template<int _N = 0>
280  struct trial_precompute
281  {
282  //typedef typename space_2_type::basis_0_type::template precompute<_N>::type type;
283  typedef typename space_2_type::basis_0_type::PreCompute type;
284  typedef boost::shared_ptr<type> ptrtype;
285  };
286 
288 
289 
294 
307  template<typename GeomapTestContext,
308  typename ExprT,
309  typename IM,
310  typename GeomapExprContext = GeomapTestContext,
311  typename GeomapTrialContext = GeomapTestContext
312  >
313  class Context //: public FormContextBase<GeomapTestContext,IM,GeomapExprContext>
314  {
315  typedef FormContextBase<GeomapTestContext,IM,GeomapExprContext> super;
316 
317  public:
318  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
319 
320 
321  typedef Context<GeomapTestContext,ExprT,IM,GeomapExprContext,GeomapTrialContext> form_context_type;
322  typedef BilinearForm<FE1, FE2, ElemContType> form_type;
323  typedef typename FE1::dof_type dof_1_type;
324  typedef typename FE2::dof_type dof_2_type;
325 
326  typedef typename form_type::value_type value_type;
327 
328 
329  typedef typename super::map_geometric_mapping_context_type map_test_geometric_mapping_context_type;
330  typedef typename super::geometric_mapping_context_ptrtype test_geometric_mapping_context_ptrtype;
331  typedef typename super::geometric_mapping_context_type test_geometric_mapping_context_type;
332  typedef typename super::geometric_mapping_type test_geometric_mapping_type;
333  typedef typename super::map_size test_map_size;
334  typedef typename super::gmc1 test_gmc1;
335  typedef typename super::left_gmc_ptrtype test_left_gmc_ptrtype;
336  typedef typename super::right_gmc_ptrtype test_right_gmc_ptrtype;
337  typedef typename super::map_left_gmc_type test_map_left_gmc_type;
338  typedef typename super::map_right_gmc_type test_map_right_gmc_type;
339 
340 
341  typedef FormContextBase<GeomapTrialContext,IM,GeomapExprContext> super2;
342 
343  typedef typename super2::map_geometric_mapping_context_type map_trial_geometric_mapping_context_type;
344  typedef typename super2::geometric_mapping_context_ptrtype trial_geometric_mapping_context_ptrtype;
345  typedef typename super2::geometric_mapping_context_type trial_geometric_mapping_context_type;
346  typedef typename super2::geometric_mapping_type trial_geometric_mapping_type;
347  typedef typename super2::map_size trial_map_size;
348  typedef typename super2::gmc1 trial_gmc1;
349  typedef typename super2::left_gmc_ptrtype trial_left_gmc_ptrtype;
350  typedef typename super2::right_gmc_ptrtype trial_right_gmc_ptrtype;
351  typedef typename super2::map_left_gmc_type trial_map_left_gmc_type;
352  typedef typename super2::map_right_gmc_type trial_map_right_gmc_type;
353 
354 
355 
356 
357 
358  static const uint16_type nDim = test_geometric_mapping_type::nDim;
359 
360  static const uint16_type nDimTest = test_space_type::mesh_type::nDim;
361  static const uint16_type nDimTrial = trial_space_type::mesh_type::nDim;
362  static const uint16_type nDimDiffBetweenTestTrial = ( nDimTest > nDimTrial )? nDimTest-nDimTrial : nDimTrial-nDimTest;
363 
364  typedef ExprT expression_type;
365 
366  typedef typename test_precompute<>::type test_precompute_type;
367  typedef typename test_precompute<>::ptrtype test_precompute_ptrtype;
368  typedef typename trial_precompute<>::type trial_precompute_type;
369  typedef typename trial_precompute<>::ptrtype trial_precompute_ptrtype;
370 
371  typedef typename FE2::fe_type trial_fe_type;
372  typedef boost::shared_ptr<trial_fe_type> trial_fe_ptrtype;
373  typedef typename trial_fe_type::template Context< trial_geometric_mapping_context_type::context,
374  trial_fe_type,
375  trial_geometric_mapping_type,
376  mesh_element_2_type> trial_fecontext_type;
377  typedef boost::shared_ptr<trial_fecontext_type> trial_fecontext_ptrtype;
378  typedef typename mpl::if_<mpl::equal_to<trial_map_size, mpl::int_<1> >,
379  mpl::identity<fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype> > >,
380  mpl::identity<fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype>,
381  fusion::pair<gmc<1>, trial_fecontext_ptrtype> > > >::type::type map_trial_fecontext_type;
382 
383  typedef fusion::map<fusion::pair<gmc<0>, trial_fecontext_ptrtype> > map_left_trial_fecontext_type;
384  typedef fusion::map<fusion::pair<trial_gmc1, trial_fecontext_ptrtype> > map_right_trial_fecontext_type;
385 
386  typedef typename FE1::fe_type test_fe_type;
387  typedef boost::shared_ptr<test_fe_type> test_fe_ptrtype;
388  typedef typename test_fe_type::template Context< test_geometric_mapping_context_type::context,
389  test_fe_type,
390  test_geometric_mapping_type,
391  mesh_element_1_type> test_fecontext_type;
392  typedef boost::shared_ptr<test_fecontext_type> test_fecontext_ptrtype;
393  typedef typename mpl::if_<mpl::equal_to<test_map_size, mpl::int_<1> >,
394  mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > >,
395  mpl::identity<fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype>,
396  fusion::pair<gmc<1>, test_fecontext_ptrtype> > > >::type::type map_test_fecontext_type;
397  typedef fusion::map<fusion::pair<gmc<0>, test_fecontext_ptrtype> > map_left_test_fecontext_type;
398  typedef fusion::map<fusion::pair<test_gmc1, test_fecontext_ptrtype> > map_right_test_fecontext_type;
399 
400  //--------------------------------------------------------------------------------------//
401 
402  typedef typename super::map_geometric_mapping_expr_context_type map_geometric_mapping_expr_context_type;
403 
404  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
405  map_left_test_fecontext_type,
406  map_left_trial_fecontext_type> eval00_expr_type;
407  typedef boost::shared_ptr<eval00_expr_type> eval00_expr_ptrtype;
408 
409  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
410  map_left_test_fecontext_type,
411  map_right_trial_fecontext_type> eval01_expr_type;
412  typedef boost::shared_ptr<eval01_expr_type> eval01_expr_ptrtype;
413 
414  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
415  map_right_test_fecontext_type,
416  map_left_trial_fecontext_type> eval10_expr_type;
417  typedef boost::shared_ptr<eval10_expr_type> eval10_expr_ptrtype;
418 
419  typedef typename ExprT::template tensor<map_geometric_mapping_expr_context_type,
420  map_right_test_fecontext_type,
421  map_right_trial_fecontext_type> eval11_expr_type;
422  typedef boost::shared_ptr<eval11_expr_type> eval11_expr_ptrtype;
423 
424 
425  static const int rep_shape = 4;//2+(eval_expr_type::shape::M-1>0)+(eval_expr_type::shape::N-1>0);
426  //typedef boost::multi_array<value_type, rep_shape> local_matrix_type;
427 
428  typedef typename FE1::dof_type test_dof_type;
429  typedef typename FE2::dof_type trial_dof_type;
430  static const int nDofPerElementTest = FE1::dof_type::nDofPerElement;
431  static const int nDofPerElementTrial = FE2::dof_type::nDofPerElement;
432  static const int nDofPerComponentTest = FE1::dof_type::fe_type::nLocalDof;
433  static const int nDofPerComponentTrial = FE2::dof_type::fe_type::nLocalDof;
434  static const int local_mat_traits = mpl::if_<mpl::equal_to<mpl::int_<nDofPerElementTrial>,mpl::int_<1> >,
435  mpl::int_<Eigen::ColMajor>,
436  mpl::int_<Eigen::RowMajor> >::type::value;
437  typedef Eigen::Matrix<value_type, nDofPerElementTest, nDofPerElementTrial,local_mat_traits> local_matrix_type;
438  typedef Eigen::Matrix<value_type, 2*nDofPerElementTest, 2*nDofPerElementTrial,Eigen::RowMajor> local2_matrix_type;
439  typedef Eigen::Matrix<value_type, nDofPerComponentTest, nDofPerComponentTrial,local_mat_traits> c_local_matrix_type;
440  typedef Eigen::Matrix<value_type, 2*nDofPerComponentTest, 2*nDofPerComponentTrial,Eigen::RowMajor> c_local2_matrix_type;
441  typedef Eigen::Matrix<int, nDofPerElementTest, 1> local_row_type;
442  typedef Eigen::Matrix<int, 2*nDofPerElementTest, 1> local2_row_type;
443  typedef Eigen::Matrix<int, nDofPerElementTrial, 1> local_col_type;
444  typedef Eigen::Matrix<int, 2*nDofPerElementTrial, 1> local2_col_type;
445 
446  typedef Eigen::Matrix<int, nDofPerComponentTest, 1> c_local_row_type;
447  typedef Eigen::Matrix<int, 2*nDofPerComponentTest, 1> c_local2_row_type;
448  typedef Eigen::Matrix<int, nDofPerComponentTrial, 1> c_local_col_type;
449  typedef Eigen::Matrix<int, 2*nDofPerComponentTrial, 1> c_local2_col_type;
450 
451 
452  public:
453 
454  template<typename Left, typename Right>
455  map_trial_fecontext_type getMap( Left left, Right right )
456  {
457  return getMap( left, right, boost::is_same<Left, map_trial_fecontext_type>() );
458  }
459  template<typename Left, typename Right>
460  map_trial_fecontext_type getMap( Left left, Right /*right*/, mpl::bool_<true> )
461  {
462  return left;
463  }
464  template<typename Left, typename Right>
465  map_trial_fecontext_type getMap( Left /*left*/, Right right, mpl::bool_<false> )
466  {
467  return map_trial_fecontext_type( right );
468  }
469 
470  template<typename Left, typename Right>
471  map_left_trial_fecontext_type getMapL( Left left, Right right )
472  {
473  return getMap( left, right, boost::is_same<Left, map_left_trial_fecontext_type>() );
474  }
475  template<typename Left, typename Right>
476  map_left_trial_fecontext_type getMapL( Left left, Right /*right*/, mpl::bool_<true> )
477  {
478  return left;
479  }
480  template<typename Left, typename Right>
481  map_left_trial_fecontext_type getMapL( Left /*left*/, Right right, mpl::bool_<false> )
482  {
483  return right;
484  }
485 
486  Context( form_type& __form,
487  map_test_geometric_mapping_context_type const& gmcTest,
488  map_trial_geometric_mapping_context_type const& _gmcTrial,
489  map_geometric_mapping_expr_context_type const & gmcExpr,
490  ExprT const& expr,
491  IM const& im );
492 
493  template<typename IMTest,typename IMTrial>
494  Context( form_type& __form,
495  map_test_geometric_mapping_context_type const& gmcTest,
496  map_trial_geometric_mapping_context_type const& _gmcTrial,
497  map_geometric_mapping_expr_context_type const & gmcExpr,
498  ExprT const& expr,
499  IM const& im,
500  IMTest const& imTest, IMTrial const& imTrial );
501 
502  template<typename IM2>
503  Context( form_type& __form,
504  map_test_geometric_mapping_context_type const& gmcTest,
505  map_trial_geometric_mapping_context_type const& _gmcTrial,
506  map_geometric_mapping_expr_context_type const & gmcExpr,
507  ExprT const& expr,
508  IM const& im,
509  IM2 const& im2 );
510 
511  template<typename IM2>
512  Context( form_type& __form,
513  map_test_geometric_mapping_context_type const& gmcTest,
514  map_trial_geometric_mapping_context_type const& _gmcTrial,
515  map_geometric_mapping_expr_context_type const & gmcExpr,
516  ExprT const& expr,
517  IM const& im,
518  IM2 const& im2,
519  mpl::int_<2> );
520 
521  size_type trialElementId( size_type trial_eid ) const
522  {
523  return trialElementId( trial_eid,mpl::int_<nDimDiffBetweenTestTrial>() );
524  }
525  size_type trialElementId( size_type trial_eid,mpl::int_<0> ) const
526  {
527  size_type idElem = trial_eid;
528  size_type domain_eid = idElem;
529  const bool test_related_to_trial = M_form.testSpace()->mesh()->isSubMeshFrom( M_form.trialSpace()->mesh() );
530  const bool trial_related_to_test = M_form.trialSpace()->mesh()->isSubMeshFrom( M_form.testSpace()->mesh() );
531  if ( test_related_to_trial )
532  {
533  domain_eid = M_form.testSpace()->mesh()->subMeshToMesh( idElem );
534  DVLOG(2) << "[test_related_to_trial] test element id: " << idElem << " trial element id : " << domain_eid << "\n";
535  }
536  if( trial_related_to_test )
537  {
538  domain_eid = M_form.trialSpace()->mesh()->meshToSubMesh( idElem );
539  DVLOG(2) << "[trial_related_to_test] test element id: " << idElem << " trial element id : " << domain_eid << "\n";
540  }
541  return domain_eid;
542  }
543  size_type trialElementId( size_type trial_eid,mpl::int_<1> ) const
544  {
545  size_type idElem = trial_eid;
546  size_type domain_eid = idElem;
547  const bool test_related_to_trial = M_form.testSpace()->mesh()->isSubMeshFrom( M_form.trialSpace()->mesh() );
548  const bool trial_related_to_test = M_form.trialSpace()->mesh()->isSubMeshFrom( M_form.testSpace()->mesh() );
549  if ( test_related_to_trial )
550  {
551  domain_eid = M_form.trialSpace()->mesh()->face( M_form.testSpace()->mesh()->subMeshToMesh( idElem )).element0().id();
552  DVLOG(2) << "[test_related_to_trial] test element id: " << idElem << " trial element id : " << domain_eid << "\n";
553  }
554  if( trial_related_to_test )
555  {
556  auto const& eltTest = M_form.testSpace()->mesh()->element(idElem);
557  std::set<size_type> idsFind;
558  for (uint16_type f=0;f< M_form.testSpace()->mesh()->numLocalFaces();++f)
559  {
560  const size_type idFind = M_form.trialSpace()->mesh()->meshToSubMesh( eltTest.face(f).id() );
561  if ( idFind != invalid_size_type_value ) idsFind.insert( idFind );
562  }
563  if ( idsFind.size()>1 ) std::cout << " TODO trialElementId " << std::endl;
564 
565  if ( idsFind.size()>0 )
566  domain_eid = *idsFind.begin();
567  else
568  domain_eid = invalid_size_type_value;
569 
570  DVLOG(2) << "[trial_related_to_test] test element id: " << idElem << " trial element id : " << domain_eid << "\n";
571  }
572  return domain_eid;
573  }
574  size_type trialElementId( typename mesh_1_type::element_iterator it ) const
575  {
576  return trialElementId( it->id() );
577  }
578  size_type trialElementId( typename mesh_1_type::element_type const& e ) const
579  {
580  return trialElementId( e.id() );
581  }
582  bool isZero( size_type i ) const
583  {
584  size_type domain_eid = trialElementId( i );
585  if ( domain_eid == invalid_size_type_value )
586  return true;
587  return false;
588  }
589  bool isZero( typename mesh_1_type::element_iterator it ) const
590  {
591  return this->isZero( it->id() );
592  }
593  bool isZero( typename mesh_1_type::element_type const& e ) const
594  {
595  return this->isZero( e.id() );
596  }
597 
598  void update( map_test_geometric_mapping_context_type const& gmcTest,
599  map_trial_geometric_mapping_context_type const& _gmcTrial,
600  map_geometric_mapping_expr_context_type const& gmcExpr );
601 
602  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
603  map_trial_geometric_mapping_context_type const& gmcTrial,
604  map_geometric_mapping_expr_context_type const& gmcExpr,
605  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad );
606 
607  void update( map_test_geometric_mapping_context_type const& gmcTest,
608  map_trial_geometric_mapping_context_type const& _gmcTrial,
609  map_geometric_mapping_expr_context_type const& gmcExpr,
610  mpl::int_<2> );
611 
612  void update( map_test_geometric_mapping_context_type const& gmcTest,
613  map_trial_geometric_mapping_context_type const& gmcTrial,
614  map_geometric_mapping_expr_context_type const& gmcExpr,
615  IM const& im )
616  {
617  M_integrator = im;
618  update( gmcTest, gmcTrial, gmcExpr );
619  }
620 
621  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
622  map_trial_geometric_mapping_context_type const& gmcTrial,
623  map_geometric_mapping_expr_context_type const& gmcExpr,
624  IM const& im,
625  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad )
626  {
627  M_integrator = im;
628  updateInCaseOfInterpolate( gmcTest, gmcTrial, gmcExpr, indexLocalToQuad );
629  }
630 
631  void update( map_test_geometric_mapping_context_type const& gmcTest,
632  map_trial_geometric_mapping_context_type const& gmcTrial,
633  map_geometric_mapping_expr_context_type const& gmcExpr,
634  IM const& im, mpl::int_<2> )
635  {
636  M_integrator = im;
637  update( gmcTest, gmcTrial, gmcExpr, mpl::int_<2>() );
638  }
639 
640  void integrate()
641  {
642  integrate( mpl::int_<fusion::result_of::size<GeomapTestContext>::type::value>() );
643  }
644  void integrateInCaseOfInterpolate( std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
645  bool isFirstExperience )
646  {
647  integrateInCaseOfInterpolate( mpl::int_<fusion::result_of::size<GeomapTestContext>::type::value>(),
648  indexLocalToQuad,
649  isFirstExperience );
650  }
651 
652 
653  void assemble()
654  {
655  assemble( fusion::at_key<gmc<0> >( M_test_gmc )->id() );
656  }
657  void assemble( size_type elt_0 );
658 
659  void assemble( mpl::int_<2> )
660  {
661  assemble( fusion::at_key<gmc<0> >( M_test_gmc )->id(),
662  fusion::at_key<test_gmc1>( M_test_gmc )->id() );
663  }
664  void assemble( size_type elt_0, size_type elt_1 );
665 
666  void assembleInCaseOfInterpolate();
667 
672  template<typename Pts>
673  void precomputeBasisAtPoints( Pts const& pts )
674  {
675  M_test_pc = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), pts ) );
676  M_trial_pc = trial_precompute_ptrtype( new trial_precompute_type( M_form.trialSpace()->fe(), pts ) );
677  }
678 
679  template<typename PtsTest,typename PtsTrial>
680  void precomputeBasisAtPoints( PtsTest const& pts1,PtsTrial const& pts2 )
681  {
682  M_test_pc = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), pts1 ) );
683  M_trial_pc = trial_precompute_ptrtype( new trial_precompute_type( M_form.trialSpace()->fe(), pts2 ) );
684  }
685 
686 
692  template<typename Pts>
693  void precomputeBasisAtPoints( uint16_type __f, permutation_1_type const& __p, Pts const& pts )
694  {
695  M_test_pc_face[__f][__p] = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), pts ) );
696  //FEELPP_ASSERT( M_test_pc_face.find(__f )->second )( __f ).error( "invalid test precompute type" );
697 
698  M_trial_pc_face[__f][__p] = trial_precompute_ptrtype( new trial_precompute_type( M_form.trialSpace()->fe(), pts ) );
699  //FEELPP_ASSERT( M_trial_pc_face.find(__f )->second )( __f ).error( "invalid trial precompute type" );
700  }
701 
702  template<typename PtsSet>
703  std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
704  precomputeTestBasisAtPoints( PtsSet const& pts )
705  {
706  typedef typename boost::is_same< permutation_1_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
707  return precomputeTestBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
708  }
709 
710  template<typename PtsSet>
711  std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
712  precomputeTestBasisAtPoints( PtsSet const& pts, mpl::bool_<false> )
713  {
714  std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> > testpc;
715  return testpc;
716  }
717 
718  template<typename PtsSet>
719  std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> >
720  precomputeTestBasisAtPoints( PtsSet const& pts, mpl::bool_<true> )
721  {
722  QuadMapped<PtsSet> qm;
723  typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
724  typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
725 
726  std::map<uint16_type, std::map<permutation_type,test_precompute_ptrtype> > testpc;
727 
728  for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
729  {
730  for ( permutation_type __p( permutation_type::IDENTITY );
731  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
732  {
733  testpc[__f][__p] = test_precompute_ptrtype( new test_precompute_type( M_form.testSpace()->fe(), ppts[__f].find( __p )->second ) );
734  }
735  }
736 
737  return testpc;
738  }
739 
740  template<typename PtsSet>
741  std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
742  precomputeTrialBasisAtPoints( PtsSet const& pts )
743  {
744  typedef typename boost::is_same< permutation_2_type, typename QuadMapped<PtsSet>::permutation_type>::type is_same_permuation_type;
745  return precomputeTrialBasisAtPoints( pts, mpl::bool_<is_same_permuation_type::value>() );
746  }
747 
748  template<typename PtsSet>
749  std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
750  precomputeTrialBasisAtPoints( PtsSet const& pts, mpl::bool_<false> )
751  {
752  std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> > trialpc;
753  return trialpc;
754  }
755 
756  template<typename PtsSet>
757  std::map<uint16_type, std::map<permutation_2_type,trial_precompute_ptrtype> >
758  precomputeTrialBasisAtPoints( PtsSet const& pts, mpl::bool_<true> )
759  {
760  QuadMapped<PtsSet> qm;
761  typedef typename QuadMapped<PtsSet>::permutation_type permutation_type;
762  typename QuadMapped<PtsSet>::permutation_points_type ppts( qm( pts ) );
763 
764  std::map<uint16_type, std::map<permutation_type,trial_precompute_ptrtype> > trialpc;
765 
766  for ( uint16_type __f = 0; __f < pts.nFaces(); ++__f )
767  {
768  for ( permutation_type __p( permutation_type::IDENTITY );
769  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
770  {
771  trialpc[__f][__p] = trial_precompute_ptrtype( new trial_precompute_type( M_form.trialSpace()->fe(), ppts[__f].find( __p )->second ) );
772  }
773  }
774 
775  return trialpc;
776  }
777 
785  test_precompute_ptrtype const& testPc( uint16_type __f,
786  permutation_1_type __p = permutation_1_type( permutation_1_type::NO_PERMUTATION ) ) const
787  {
788  if ( __f == invalid_uint16_type_value )
789  return M_test_pc;
790 
791  //FEELPP_ASSERT( M_test_pc_face.find(__f )->second )( __f ).error( "invalid test precompute type" );
792  return M_test_pc_face.find( __f )->second.find( __p )->second;
793  }
794 
801  trial_precompute_ptrtype const& trialPc( uint16_type __f,
802  permutation_2_type __p = permutation_2_type( permutation_2_type::NO_PERMUTATION ) ) const
803  {
804  if ( __f == invalid_uint16_type_value )
805  return M_trial_pc;
806 
807  //FEELPP_ASSERT( M_trial_pc_face.find(__f )->second )( __f ).error( "invalid trial precompute type" );
808  return M_trial_pc_face.find( __f )->second.find( __p )->second;
809  }
810 
811 
812  private:
813 
814  void update( map_test_geometric_mapping_context_type const& gmcTest,
815  map_trial_geometric_mapping_context_type const& gmcTrial,
816  map_geometric_mapping_expr_context_type const& gmcExpr,
817  mpl::bool_<false> );
818 
819  void update( map_test_geometric_mapping_context_type const& gmcTest,
820  map_trial_geometric_mapping_context_type const& gmcTrial,
821  map_geometric_mapping_expr_context_type const& gmcExpr,
822  mpl::bool_<true> );
823 
824  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
825  map_trial_geometric_mapping_context_type const& gmcTrial,
826  map_geometric_mapping_expr_context_type const& gmcExpr,
827  mpl::bool_<false> );
828 
829  void updateInCaseOfInterpolate( map_test_geometric_mapping_context_type const& gmcTest,
830  map_trial_geometric_mapping_context_type const& gmcTrial,
831  map_geometric_mapping_expr_context_type const& gmcExpr,
832  mpl::bool_<true> );
833 
834 
835  void integrate( mpl::int_<1> );
836 
837  void integrate( mpl::int_<2> );
838 
839  void integrateInCaseOfInterpolate( mpl::int_<1>,
840  std::vector<boost::tuple<size_type,size_type> > const& indexLocalToQuad,
841  bool isFirstExperience );
842  private:
843 
844  form_type& M_form;
845  const list_block_type& M_lb;
846  dof_1_type* M_test_dof;
847  dof_2_type* M_trial_dof;
848 
849  test_precompute_ptrtype M_test_pc;
850  std::map<uint16_type, std::map<permutation_1_type,test_precompute_ptrtype> > M_test_pc_face;
851  trial_precompute_ptrtype M_trial_pc;
852  std::map<uint16_type, std::map<permutation_2_type, trial_precompute_ptrtype> > M_trial_pc_face;
853 
854 
855  map_test_geometric_mapping_context_type M_test_gmc;
856  map_trial_geometric_mapping_context_type M_trial_gmc;
857 
858  map_test_fecontext_type M_test_fec;
859  map_left_test_fecontext_type M_test_fec0;
860  map_right_test_fecontext_type M_test_fec1;
861  map_trial_fecontext_type M_trial_fec;
862  map_left_trial_fecontext_type M_trial_fec0;
863  map_right_trial_fecontext_type M_trial_fec1;
864 
865 
866  local_matrix_type M_rep;
867  local2_matrix_type M_rep_2;
868  local_row_type M_local_rows;
869  local2_row_type M_local_rows_2;
870  local_col_type M_local_cols;
871  local2_col_type M_local_cols_2;
872  local_row_type M_local_rowsigns;
873  local2_row_type M_local_rowsigns_2;
874  local_col_type M_local_colsigns;
875  local2_col_type M_local_colsigns_2;
876 
877  c_local_matrix_type M_c_rep;
878  c_local2_matrix_type M_c_rep_2;
879  c_local_row_type M_c_local_rows;
880  c_local2_row_type M_c_local_rows_2;
881  c_local_col_type M_c_local_cols;
882  c_local2_col_type M_c_local_cols_2;
883  c_local_row_type M_c_local_rowsigns;
884  c_local2_row_type M_c_local_rowsigns_2;
885  c_local_col_type M_c_local_colsigns;
886  c_local2_col_type M_c_local_colsigns_2;
887 
888  eval00_expr_ptrtype M_eval_expr00;
889  eval01_expr_ptrtype M_eval_expr01;
890  eval10_expr_ptrtype M_eval_expr10;
891  eval11_expr_ptrtype M_eval_expr11;
892 
893  IM M_integrator;
894 
895  }; // Context
896 
898 
902 
906  BilinearForm( space_1_ptrtype const& __X1,
907  space_2_ptrtype const& __X2,
908  matrix_ptrtype& __M,
909  size_type rowstart = 0,
910  size_type colstart = 0,
911  bool build = true,
912  bool do_threshold = false,
913  value_type threshold = type_traits<value_type>::epsilon(),
914  size_type graph_hints = Pattern::COUPLED );
915 
916  BilinearForm( space_1_ptrtype const& __X1,
917  space_2_ptrtype const& __X2,
918  matrix_ptrtype& __M,
919  list_block_type const& __lb,
920  size_type rowstart = 0,
921  size_type colstart = 0,
922  bool do_threshold = false,
923  value_type threshold = type_traits<value_type>::epsilon(),
924  size_type graph_hints = Pattern::COUPLED );
925 
926  BilinearForm( BilinearForm const& __vf )
927  :
928  M_pattern( __vf.M_pattern ),
929  M_X1( __vf.M_X1 ),
930  M_X2( __vf.M_X2 ),
931  M_matrix( __vf.M_matrix ),
932  M_lb( __vf.M_lb ),
933  M_row_startInMatrix( __vf.M_row_startInMatrix ),
934  M_col_startInMatrix( __vf.M_col_startInMatrix ),
935  M_do_threshold( __vf.M_do_threshold ),
936  M_threshold( __vf.M_threshold )
937  {}
938 
939  ~BilinearForm()
940  {}
941 
942 
943 
945 
949 
953  BilinearForm&
954  operator=( BilinearForm const& form )
955  {
956  if ( this != &form )
957  {
958  M_pattern = form.M_pattern;
959  M_X1 = form.M_X1;
960  M_X2 = form.M_X2;
961  M_matrix = form.M_matrix;
962  M_row_startInMatrix = form.M_row_startInMatrix;
963  M_col_startInMatrix = form.M_col_startInMatrix;
964  M_lb = form.M_lb;
965  }
966 
967  return *this;
968  }
969 
973  template <class ExprT>
974  BilinearForm& operator=( Expr<ExprT> const& expr );
975 
979  template <class ExprT>
980  BilinearForm& operator+=( Expr<ExprT> const& expr );
981 
982  BilinearForm& operator+=( BilinearForm& a )
983  {
984  if ( this == &a )
985  return *this;
986 
987  M_matrix->addMatrix( 1.0, a.M_matrix );
988 
989  return *this;
990  }
991 
999  value_type operator()( element_1_type const& __v, element_2_type const& __u ) const
1000  {
1001  return M_matrix->energy( __v, __u );
1002  }
1003 
1007  // value_type& operator()( size_type i, size_type j ) { return M_matrix( i, j ); }
1009 
1013 
1017  size_type pattern() const
1018  {
1019  return M_pattern;
1020  }
1021 
1026  bool isPatternCoupled() const
1027  {
1028  Feel::Context ctx( M_pattern );
1029  return ctx.test( Pattern::COUPLED );
1030  }
1031 
1035  bool isPatternDefault() const
1036  {
1037  Feel::Context ctx( M_pattern );
1038  return ctx.test( Pattern::DEFAULT );
1039  }
1040 
1044  bool isPatternNeighbor() const
1045  {
1046  Feel::Context ctx( M_pattern );
1047  return ctx.test( Pattern::EXTENDED );
1048  }
1049  bool isPatternExtended() const
1050  {
1051  Feel::Context ctx( M_pattern );
1052  return ctx.test( Pattern::EXTENDED );
1053  }
1054 
1055  bool isPatternSymmetric() const
1056  {
1057  Feel::Context ctx( M_pattern );
1058  return ctx.test( Pattern::SYMMETRIC );
1059  }
1060 
1064  space_2_ptrtype const& trialSpace() const
1065  {
1066  return M_X2;
1067  }
1068 
1072  space_1_ptrtype const& testSpace() const
1073  {
1074  return M_X1;
1075  }
1076 
1080  gm_1_ptrtype const& gm() const
1081  {
1082  return M_X1->gm();
1083  }
1084 
1088  gm1_1_ptrtype const& gm1() const
1089  {
1090  return M_X1->gm1();
1091  }
1092 
1096  gm_2_ptrtype const& gmTrial() const
1097  {
1098  return M_X2->gm();
1099  }
1100 
1104  gm1_2_ptrtype const& gm1Trial() const
1105  {
1106  return M_X2->gm1();
1107  }
1108 
1109 
1113  matrix_type const& matrix() const
1114  {
1115  return *M_matrix;
1116  }
1117 
1118  matrix_type& matrix()
1119  {
1120  return *M_matrix;
1121  }
1122 
1123  matrix_ptrtype const& matrixPtr() const
1124  {
1125  return M_matrix;
1126  }
1127 
1128  matrix_ptrtype& matrixPtr()
1129  {
1130  return M_matrix;
1131  }
1132 
1133  list_block_type const& blockList() const
1134  {
1135  return M_lb;
1136  }
1137 
1138  size_type rowStartInMatrix() const
1139  {
1140  return M_row_startInMatrix;
1141  }
1142 
1143  size_type colStartInMatrix() const
1144  {
1145  return M_col_startInMatrix;
1146  }
1147 
1151  value_type threshold() const
1152  {
1153  return M_threshold;
1154  }
1155 
1159  bool doThreshold( value_type const& v ) const
1160  {
1161  return ( math::abs( v ) > M_threshold );
1162  }
1163 
1167  bool doThreshold() const
1168  {
1169  return M_do_threshold;
1170  }
1171 
1173 
1177 
1182  void setThreshold( value_type eps )
1183  {
1184  M_threshold = eps;
1185  }
1186 
1191  void setDoThreshold( bool do_threshold )
1192  {
1193  M_do_threshold = do_threshold;
1194  }
1195 
1197 
1201 
1202 
1203  // close matrix
1204  void close() { M_matrix->close(); }
1205 
1206 
1216  void zeroRows( std::vector<int> const& __dofs,
1217  Vector<value_type> const& __values,
1218  Vector<value_type>& rhs,
1219  Feel::Context const& on_context );
1220 
1225  void add( size_type i, size_type j, value_type const& v )
1226  {
1227  if ( M_do_threshold )
1228  {
1229  if ( doThreshold( v ) )
1230  M_matrix->add( i+this->rowStartInMatrix(),
1231  j+this->colStartInMatrix(),
1232  v );
1233  }
1234 
1235  else
1236  M_matrix->add( i+this->rowStartInMatrix(),
1237  j+this->colStartInMatrix(),
1238  v );
1239 
1240  }
1245  void addMatrix( int* rows, int nrows,
1246  int* cols, int ncols,
1247  value_type* data )
1248  {
1249  if ( this->rowStartInMatrix()!=0 )
1250  for ( int i=0; i<nrows; ++i )
1251  rows[i]+=this->rowStartInMatrix();
1252 
1253  if ( this->colStartInMatrix()!=0 )
1254  for ( int i=0; i<ncols; ++i )
1255  cols[i]+=this->colStartInMatrix();
1256 
1257  M_matrix->addMatrix( rows, nrows, cols, ncols, data );
1258  }
1259 
1260 
1261 
1266  void set( size_type i, size_type j, value_type const& v )
1267  {
1268  M_matrix->set( i, j, v );
1269  }
1270 
1271  void addToNOz( size_type i, size_type n )
1272  {
1273  M_n_oz[i] += n;
1274  }
1275  void addToNNz( size_type i, size_type n )
1276  {
1277  M_n_nz[i] += n;
1278  }
1279  size_type nOz( size_type i ) const
1280  {
1281  return M_n_oz[i];
1282  }
1283  size_type nNz( size_type i ) const
1284  {
1285  return M_n_nz[i];
1286  }
1287 
1288  BOOST_PARAMETER_MEMBER_FUNCTION( ( typename Backend<value_type>::solve_return_type ),
1289  solve,
1290  tag,
1291  ( required
1292  ( in_out( solution ),* )
1293  ( rhs, * ) )
1294  ( optional
1295  ( name, ( std::string ), "" )
1296  ( kind, ( BackendType ), BACKEND_PETSC )
1297  ( rebuild, ( bool ), false )
1298  ) )
1299  {
1300  return backend( _name=name, _kind=kind, _rebuild=rebuild )->solve( _matrix=this->matrixPtr(), _rhs=rhs.vectorPtr(),
1301  _solution=solution);
1302  }
1303 
1304  BOOST_PARAMETER_MEMBER_FUNCTION( ( typename Backend<value_type>::solve_return_type ),
1305  solveb,
1306  tag,
1307  ( required
1308  ( in_out( solution ),* )
1309  ( rhs, * )
1310  ( backend, *) )
1311  ( optional
1312  ( prec, ( preconditioner_ptrtype ),
1313  preconditioner( _prefix=backend->prefix(),
1314  _matrix=this->matrixPtr(),
1315  _pc=backend->pcEnumType()/*LU_PRECOND*/,
1316  _pcfactormatsolverpackage=backend->matSolverPackageEnumType(),
1317  _backend=backend ) )
1318  ) )
1319  {
1320  return backend->solve( _matrix=this->matrixPtr(), _rhs=rhs.vectorPtr(),
1321  _solution=solution, _prec = prec );
1322  }
1323 
1325 
1326 protected:
1327 
1328  BilinearForm();
1329 
1330 private:
1331 
1332 
1333  template <class ExprT> void assign( Expr<ExprT> const& __expr, bool init, mpl::bool_<false> );
1334  template <class ExprT> void assign( Expr<ExprT> const& __expr, bool init, mpl::bool_<true> );
1335  template <class ExprT> void assign( Expr<ExprT> const& __expr, mpl::bool_<true>, mpl::bool_<true> );
1336  template <class ExprT> void assign( Expr<ExprT> const& __expr, mpl::bool_<true>, mpl::bool_<false> );
1337  template <class ExprT> void assign( Expr<ExprT> const& __expr, mpl::bool_<true>, mpl::bool_<false>, mpl::bool_<true> );
1338  template <class ExprT> void assign( Expr<ExprT> const& __expr, mpl::bool_<true>, mpl::bool_<false>, mpl::bool_<false> );
1339 
1340 private:
1341 
1342  size_type M_pattern;
1343 
1344  space_1_ptrtype M_X1;
1345  space_2_ptrtype M_X2;
1346 
1347  matrix_ptrtype M_matrix;
1348 
1349  bool M_do_build;
1350 
1351  list_block_type M_lb;
1352  size_type M_row_startInMatrix,M_col_startInMatrix;
1353 
1354  bool M_do_threshold;
1355  value_type M_threshold;
1356 
1357  std::vector<size_type> M_n_nz;
1358  std::vector<size_type> M_n_oz;
1359 
1360 };
1361 template<typename FE1, typename FE2, typename ElemContType>
1362 BilinearForm<FE1, FE2, ElemContType>::BilinearForm( space_1_ptrtype const& Xh,
1363  space_2_ptrtype const& Yh,
1364  matrix_ptrtype& __M,
1365  size_type rowstart,
1366  size_type colstart,
1367  bool build,
1368  bool do_threshold,
1369  value_type threshold,
1370  size_type graph_hints )
1371  :
1372  M_pattern( graph_hints ),
1373  M_X1( Xh ),
1374  M_X2( Yh ),
1375  M_matrix( __M ),
1376  M_do_build( build ),
1377  M_lb(),
1378  M_row_startInMatrix( rowstart ),
1379  M_col_startInMatrix( colstart ),
1380  M_do_threshold( do_threshold ),
1381  M_threshold( threshold )
1382 {
1383  boost::timer tim;
1384  DVLOG(2) << "begin constructor with default listblock\n";
1385 
1386  if ( !M_matrix ) M_matrix = backend()->newMatrix( _test=M_X1, _trial=M_X2 );
1387  M_lb.push_back( Block ( 0, 0, 0, 0 ) );
1388 
1389  DVLOG(2) << " - form init in " << tim.elapsed() << "\n";
1390  DVLOG(2) << "begin constructor with default listblock done\n";
1391 }
1392 
1393 template<typename FE1, typename FE2, typename ElemContType>
1394 BilinearForm<FE1, FE2, ElemContType>::BilinearForm( space_1_ptrtype const& Xh,
1395  space_2_ptrtype const& Yh,
1396  matrix_ptrtype& __M,
1397  list_block_type const& __lb,
1398  size_type rowstart,
1399  size_type colstart,
1400  bool do_threshold ,
1401  value_type threshold,
1402  size_type graph_hints )
1403  :
1404  M_pattern( graph_hints ),
1405  M_X1( Xh ),
1406  M_X2( Yh ),
1407  M_matrix( __M ),
1408  M_do_build( false ),
1409  M_lb( __lb ),
1410  M_row_startInMatrix( rowstart ),
1411  M_col_startInMatrix( colstart ),
1412  M_do_threshold( do_threshold ),
1413  M_threshold( threshold )
1414 {
1415  if ( !M_matrix ) M_matrix = backend()->newMatrix( _test=M_X1, _trial=M_X2 );
1416 }
1417 
1418 template<typename FE1, typename FE2, typename ElemContType>
1419 template<typename ExprT>
1420 void
1421 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1422  bool init,
1423  mpl::bool_<false> )
1424 {
1425  if ( init )
1426  {
1427  DVLOG(2) << "[BilinearForm::assign<false>] start\n";
1428  typedef ublas::matrix_range<matrix_type> matrix_range_type;
1429  typename list_block_type::const_iterator __bit = M_lb.begin();
1430  typename list_block_type::const_iterator __ben = M_lb.end();
1431 
1432  for ( ; __bit != __ben; ++__bit )
1433  {
1434  size_type g_ic_start = M_row_startInMatrix + __bit->globalRowStart();
1435  size_type g_jc_start = M_col_startInMatrix + __bit->globalColumnStart();
1436  M_matrix->zero( g_ic_start, g_ic_start + M_X1->nDof(),
1437  g_jc_start, g_jc_start + M_X2->nDof() );
1438  }
1439  }
1440 
1441  __expr.assemble( M_X1, M_X2, *this );
1442 
1443  DVLOG(2) << "[BilinearForm::assign<false>] stop\n";
1444 
1445 }
1446 template<typename FE1, typename FE2, typename ElemContType>
1447 template<typename ExprT>
1448 void
1449 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1450  bool init,
1451  mpl::bool_<true> )
1452 {
1453  DVLOG(2) << "BilinearForm::assign() start loop on test spaces\n";
1454 
1455  if ( init ) M_matrix->zero();
1456 
1457  assign( __expr, mpl::bool_<true>(), mpl::bool_<( FE1::nSpaces > 1 && FE2::nSpaces > 1 )>() );
1458  DVLOG(2) << "BilinearForm::assign() stop loop on test spaces\n";
1459 
1460 }
1461 template<typename FE1, typename FE2, typename ElemContType>
1462 template<typename ExprT>
1463 void
1464 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1465  mpl::bool_<true>,
1466  mpl::bool_<true> )
1467 {
1468  fusion::for_each( M_X1->functionSpaces(), make_bfassign2( *this, __expr ) );
1469 }
1470 template<typename FE1, typename FE2, typename ElemContType>
1471 template<typename ExprT>
1472 void
1473 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1474  mpl::bool_<true>,
1475  mpl::bool_<false> )
1476 {
1477  assign( __expr, mpl::bool_<true>(), mpl::bool_<false>(), mpl::bool_<( FE1::nSpaces > 1 )>() );
1478 }
1479 
1480 template<typename FE1, typename FE2, typename ElemContType>
1481 template<typename ExprT>
1482 void
1483 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1484  mpl::bool_<true>,
1485  mpl::bool_<false>,
1486  mpl::bool_<true> )
1487 {
1488  fusion::for_each( M_X1->functionSpaces(), make_bfassign3( *this, __expr, M_X2, 0 ) );
1489 }
1490 template<typename FE1, typename FE2, typename ElemContType>
1491 template<typename ExprT>
1492 void
1493 BilinearForm<FE1, FE2, ElemContType>::assign( Expr<ExprT> const& __expr,
1494  mpl::bool_<true>,
1495  mpl::bool_<false>,
1496  mpl::bool_<false> )
1497 {
1498  fusion::for_each( M_X2->functionSpaces(), make_bfassign1( *this, __expr, M_X1, 0 ) );
1499 
1500 }
1501 template<typename FE1, typename FE2, typename ElemContType>
1502 template<typename ExprT>
1503 BilinearForm<FE1, FE2, ElemContType>&
1504 BilinearForm<FE1, FE2, ElemContType>::operator=( Expr<ExprT> const& __expr )
1505 {
1506  // loop(fusion::for_each) over sub-functionspaces in SpaceType
1507  // pass expression and initialize
1508  this->assign( __expr, true, mpl::bool_<mpl::or_< mpl::bool_< ( FE1::nSpaces > 1 )>,
1509  mpl::bool_< ( FE2::nSpaces > 1 )> >::type::value >() );
1510  return *this;
1511 }
1512 
1513 template<typename FE1, typename FE2, typename ElemContType>
1514 template<typename ExprT>
1515 BilinearForm<FE1, FE2, ElemContType>&
1516 BilinearForm<FE1, FE2, ElemContType>::operator+=( Expr<ExprT> const& __expr )
1517 {
1518  DVLOG(2) << "[BilinearForm::operator+=] start\n";
1519  this->assign( __expr, false, mpl::bool_<mpl::or_< mpl::bool_< ( FE1::nSpaces > 1 )>,
1520  mpl::bool_< ( FE2::nSpaces > 1 )> >::type::value >() );
1521  DVLOG(2) << "[BilinearForm::operator+=] stop\n";
1522  return *this;
1523 }
1524 
1525 
1526 template<typename FE1, typename FE2, typename ElemContType>
1527 void
1528 BilinearForm<FE1,FE2,ElemContType>::zeroRows( std::vector<int> const& __dofs,
1529  Vector<value_type> const&__values,
1530  Vector<value_type>& rhs,
1531  Feel::Context const& on_context )
1532 {
1533  M_matrix->zeroRows( __dofs, __values, rhs, on_context );
1534 }
1535 
1536 
1537 template<typename BFType, typename ExprType, typename TestSpaceType>
1538 template<typename SpaceType>
1539 void BFAssign1<BFType,ExprType,TestSpaceType>::operator()( boost::shared_ptr<SpaceType> const& trial ) const
1540 {
1541  if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
1542  {
1543  DVLOG(2) << "[BFAssign1::operator()] expression has test functions index "
1544  << M_test_index << " : "
1545  << ExprType::template HasTestFunction<typename TestSpaceType::reference_element_type>::result << " (0:no, 1:yes)\n";
1546  DVLOG(2) << "[BFAssign1::operator()] expression has trial functions index "
1547  << M_trial_index << " :"
1548  << ExprType::template HasTrialFunction<typename SpaceType::reference_element_type>::result << " (0:no, 1:yes)\n";
1549 
1550  if ( !ExprType::template HasTestFunction<typename TestSpaceType::reference_element_type>::result ||
1551  !ExprType::template HasTrialFunction<typename SpaceType::reference_element_type>::result )
1552  {
1553  ++M_trial_index;
1554  return;
1555  }
1556 
1557  DVLOG(2) << "[BFAssign1::operator()] terms found with "
1558  << "testindex: " << M_test_index << " trialindex: " << M_trial_index << "\n";
1559  typedef SpaceType trial_space_type;
1560  typedef TestSpaceType test_space_type;
1561 
1562  Feel::vf::list_block_type list_block;
1563 
1564  if ( M_bf.testSpace()->worldsComm()[M_test_index].globalSize()>1 )
1565  {
1566  //DVLOG(2) << "[BFAssign1::operator()] block: " << block << "\n";
1567  if (M_bf.testSpace()->hasEntriesForAllSpaces())
1568  list_block.push_back( Feel::vf::Block( 0, 0,
1569  M_bf.testSpace()->nLocalDofStart( M_test_index ),
1570  M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1571  else
1572  list_block.push_back( Feel::vf::Block( 0, 0, 0, M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1573 
1574  }
1575 
1576  else
1577  {
1578  //DVLOG(2) << "[BFAssign1::operator()] block: " << block << "\n";
1579  list_block.push_back( Feel::vf::Block( 0,0,
1580  M_bf.testSpace()->nDofStart( M_test_index ),
1581  M_bf.trialSpace()->nDofStart( M_trial_index ) ) );
1582  }
1583 
1584  typedef typename BFType::matrix_type matrix_type;
1585  typedef Feel::vf::detail::BilinearForm<test_space_type,
1586  trial_space_type,
1587  ublas::vector_range<ublas::vector<double> > > bf_type;
1588 
1589  bf_type bf( M_test,trial, M_bf.matrixPtr(), list_block, M_bf.rowStartInMatrix(), M_bf.colStartInMatrix(), M_bf.doThreshold(), M_bf.threshold(), M_bf.pattern() );
1590 
1591  bf += M_expr;
1592  }
1593 
1594  ++M_trial_index;
1595 }
1596 
1597 template<typename BFType, typename ExprType, typename TrialSpaceType>
1598 template<typename SpaceType>
1599 void BFAssign3<BFType,ExprType,TrialSpaceType>::operator()( boost::shared_ptr<SpaceType> const& test ) const
1600 {
1601  if ( M_bf.testSpace()->worldsComm()[M_test_index].isActive() )
1602  {
1603 
1604  DVLOG(2) << "[BFAssign3::operator()] expression has trial functions index "
1605  << M_test_index << " : "
1606  << ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result << " (0:no, 1:yes)\n";
1607  DVLOG(2) << "[BFAssign3::operator()] expression has test functions index "
1608  << M_trial_index << " :"
1609  << ExprType::template HasTrialFunction<typename TrialSpaceType::reference_element_type>::result << " (0:no, 1:yes)\n";
1610 
1611  if ( !ExprType::template HasTrialFunction<typename TrialSpaceType::reference_element_type>::result ||
1612  !ExprType::template HasTestFunction<typename SpaceType::reference_element_type>::result )
1613  {
1614  ++M_test_index;
1615  return;
1616  }
1617 
1618  DVLOG(2) << "[BFAssign3::operator()] terms found with "
1619  << "testindex: " << M_test_index << " trialindex: " << M_trial_index << "\n";
1620  typedef SpaceType test_space_type;
1621  typedef TrialSpaceType trial_space_type;
1622 
1623  //Feel::vf::Block block ( 0, 0,
1624  // M_bf.testSpace()->nDofStart( M_test_index ),
1625  // M_bf.trialSpace()->nDofStart( M_trial_index ) );
1626  Feel::vf::list_block_type list_block;
1627 
1628  // with mpi, dof start to 0 (thanks to the LocalToGlobal mapping).
1629  if ( M_bf.testSpace()->worldsComm()[M_test_index].globalSize()>1 )
1630  {
1631  //DVLOG(2) << "[BFAssign1::operator()] block: " << block << "\n";
1632  if (M_bf.testSpace()->hasEntriesForAllSpaces())
1633  list_block.push_back( Feel::vf::Block( 0, 0,
1634  M_bf.testSpace()->nLocalDofStart( M_test_index ),
1635  M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1636  else
1637  list_block.push_back( Feel::vf::Block( 0, 0, 0, M_bf.trialSpace()->nLocalDofStart( M_trial_index ) ) );
1638 
1639  }
1640 
1641  else
1642  {
1643  //DVLOG(2) << "[BFAssign1::operator()] block: " << block << "\n";
1644  list_block.push_back( Feel::vf::Block( 0,0,
1645  M_bf.testSpace()->nDofStart( M_test_index ),
1646  M_bf.trialSpace()->nDofStart( M_trial_index ) ) );
1647  }
1648 
1649  typedef typename BFType::matrix_type matrix_type;
1650  typedef Feel::vf::detail::BilinearForm<test_space_type,
1651  trial_space_type,
1652  ublas::vector_range<ublas::vector<double> > > bf_type;
1653 
1654  bf_type bf( test, M_trial, M_bf.matrixPtr(), list_block, M_bf.rowStartInMatrix(), M_bf.colStartInMatrix(), M_bf.doThreshold(), M_bf.threshold(), M_bf.pattern() );
1655 
1656  bf += M_expr;
1657  }
1658 
1659  ++M_test_index;
1660 }
1661 
1662 
1663 } // detail
1664 } // vf
1666 } // feel
1667 
1669 
1670 #endif /* __BilinearForm_H */

Generated on Fri Oct 25 2013 14:24:04 for Feel++ by doxygen 1.8.4