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
operatorlinearcomposite.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2013-04-26
7 
8  Copyright (C) 2013 Feel++ Consortium
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
30 #ifndef __OperatorLinearComposite_H
31 #define __OperatorLinearComposite_H 1
32 
33 #include <feel/feel.hpp>
34 #include <feel/feelcore/pslogger.hpp>
35 
36 namespace Feel
37 {
38 
39 template<class DomainSpace, class DualImageSpace>
40 class OperatorLinearComposite : public OperatorLinear<DomainSpace, DualImageSpace>
41 {
42 
43  typedef OperatorLinear<DomainSpace,DualImageSpace> super;
44  typedef boost::shared_ptr<super> super_ptrtype;
45 
46 public :
47 
48  typedef typename super::domain_space_type domain_space_type;
49  typedef typename super::dual_image_space_type dual_image_space_type;
50  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
51  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
52 
53  typedef typename super::backend_type backend_type;
54  typedef typename super::backend_ptrtype backend_ptrtype;
55 
56  typedef typename super::matrix_type matrix_type;
57  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
58 
59  typedef FsFunctionalLinear<DualImageSpace> image_element_type;
60  typedef typename super::domain_element_type domain_element_type;
61  typedef typename super::dual_image_element_type dual_image_element_type;
62 
63  typedef typename super::domain_element_range_type domain_element_range_type;
64  typedef typename super::domain_element_slice_type domain_element_slice_type;
65  typedef typename super::dual_image_element_range_type dual_image_element_range_type;
66  typedef typename super::dual_image_element_slice_type dual_image_element_slice_type;
67 
68  typedef OperatorLinearComposite<DomainSpace, DualImageSpace> this_type;
69  typedef boost::shared_ptr<this_type> this_ptrtype;
70 
71  OperatorLinearComposite (domain_space_ptrtype domainSpace,
72  dual_image_space_ptrtype dualImageSpace,
73  backend_ptrtype backend,
74  size_type pattern=Pattern::COUPLED)
75  :
76  super( domainSpace,dualImageSpace,backend , false ),
77  M_backend( backend ),
78  M_pattern( pattern )
79  {}
80 
81 
82  //if we have a list of operators
83  //i.e. \sum_{q=0}^Q Aq(.,.)
84  void addElement( super_ptrtype const& op )
85  {
86  int size = M_operators1.size();
87  M_operators1.insert( std::make_pair( size , op ) ) ;
88  }
89 
90 
91  void addElement( int q, super_ptrtype const& op )
92  {
93  M_operators1.insert( std::make_pair( q , op ) ) ;
94  }
95 
96  void addList( std::vector< super_ptrtype > const& vec )
97  {
98  int q_max=vec.size();
99  for(int q=0; q<q_max; q++)
100  this->addElement( q , vec[q] );
101  }
102 
103  void displayOperatorsNames()
104  {
105  int size1 = M_operators1.size();
106  int size2 = M_operators2.size();
107 
108  LOG( INFO ) << " the composite operator linear "<<this->name()<<" has following operators : ";
109 
110  if( size1 > 0 )
111  {
112  auto end = M_operators1.end();
113  for(auto it=M_operators1.begin(); it!=end; it++)
114  LOG(INFO)<<it->second->name();
115  }
116  else
117  {
118  auto end = M_operators2.end();
119  for(auto it=M_operators2.begin(); it!=end; it++)
120  LOG(INFO)<<it->second->name();
121  }
122  }
123 
124  int size()
125  {
126  int size1 = M_operators1.size();
127  int size2 = M_operators2.size();
128  return size1+size2;
129  }
130 
131  //if we have a list of list of operators
132  //i.e. \sum_{q=0}^Q \sum_{m=0}^M A_{qm}(.,.)
133  void addElement( boost::tuple<int,int> const& qm , super_ptrtype const& op )
134  {
135  M_operators2.insert( std::make_pair( qm , op ) );
136  }
137 
138  void addList( std::vector< std::vector< super_ptrtype > > const& vec )
139  {
140  int q_max = vec.size();
141  for(int q=0; q<q_max; q++)
142  {
143  int m_max = vec[q].size();
144  for(int m=0; m<m_max; m++)
145  {
146  auto tuple = boost::make_tuple( q , m );
147  this->addElement( tuple , vec[q][m] );
148  }
149  }
150  }
151 
152  void setScalars( std::vector<double> scalars )
153  {
154  M_scalars1 = scalars;
155  }
156 
157  void setScalars( std::vector< std::vector< double > > scalars )
158  {
159  M_scalars2 = scalars;
160  }
161 
162  void sumAllMatrices(matrix_ptrtype & matrix, bool use_scalar_one=false ) const
163  {
164  int size1 = M_operators1.size();
165  int size2 = M_operators2.size();
166  bool size_error=false;
167  if( size1 > 0 && size2 > 0 )
168  size_error=true;
169  if( (size1 + size2) == 0 )
170  size_error=true;
171 
172  FEELPP_ASSERT( !size_error )( size1 )( size2 ).error( "OperatorLinearComposite has no elements, or both maps have elements" );
173 
174  if( size1 > 0 )
175  sumAllMatrices1( matrix, use_scalar_one );
176  else
177  sumAllMatrices2( matrix, use_scalar_one );
178  }
179 
180  //return the sum of matrices given
181  //by all opertors in M_vectors
182  //arguments : a vector of scalars and a bool ( use scalar=1 if true )
183  void sumAllMatrices1( matrix_ptrtype & matrix, bool use_scalar_one=false ) const
184  {
185  int size1 = M_operators1.size();
186  int size2 = M_operators2.size();
187 
188  FEELPP_ASSERT( size1 > 0 )( size1 )( size2 ).error( "OperatorLinearComposite has no elements" );
189 
190  matrix->zero();
191  auto temp_matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
192 
193  auto end = M_operators1.end();
194  for(auto it=M_operators1.begin(); it!=end; ++it)
195  {
196  int position = it->first;
197  double scalar=1;
198  if( ! use_scalar_one )
199  scalar = M_scalars1[position];
200  it->second->matPtr(temp_matrix);
201  matrix->addMatrix( scalar , temp_matrix );
202  }
203 
204  }//sumAllMatrices
205 
206  //return the sum of matrices given
207  //by all opertors in M_vectors
208  //arguments : a vector of vector of scalars and a bool ( use scalar=1 if true )
209  void sumAllMatrices2( matrix_ptrtype & matrix, bool use_scalar_one=false ) const
210  {
211  int size1 = M_operators1.size();
212  int size2 = M_operators2.size();
213 
214  FEELPP_ASSERT( size2 > 0 )( size1 )( size2 ).error( "OperatorLinearComposite has no elements" );
215 
216  matrix->zero();
217  auto temp_matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
218 
219  auto end = M_operators2.end();
220  for(auto it=M_operators2.begin(); it!=end; ++it)
221  {
222  auto position = it->first;
223  int q = position.template get<0>();
224  int m = position.template get<1>();
225  double scalar=1;
226  if( ! use_scalar_one )
227  scalar = M_scalars2[q][m];
228  it->second->matPtr(temp_matrix);
229  matrix->addMatrix( scalar , temp_matrix );
230  }
231 
232  temp_matrix.reset();
233  }//sumAllMatrices
234 
235 
236  //Access to a specific matrix
237  void matrixPtr(int q , matrix_ptrtype& matrix)
238  {
239  int q_max = M_operators1.size();
240  FEELPP_ASSERT( q < q_max )( q_max )( q ).error( "OperatorLinearComposite has not enough elements" );
241  //fill matrix
242  M_operators1.template at(q)->matPtr(matrix);
243  }
244 
245  void matrixPtr(int q , matrix_ptrtype& matrix) const
246  {
247  int q_max = M_operators1.size();
248  FEELPP_ASSERT( q < q_max )( q_max )( q ).error( "OperatorLinearComposite has not enough elements" );
249  //fill matrix
250  M_operators1.template at(q)->matPtr(matrix);
251  }
252 
253  //Access to a specific matrix
254  void matrixPtr(int q, int m , matrix_ptrtype& matrix)
255  {
256  //fill matrix
257  auto tuple = boost::make_tuple(q,m);
258  M_operators2.template at(tuple)->matPtr(matrix);
259  }
260 
261  void matrixPtr(int q, int m , matrix_ptrtype& matrix) const
262  {
263  auto tuple = boost::make_tuple(q,m);
264  M_operators2.template at(tuple)->matPtr(matrix);
265  }
266 
267 
268  //Acces to operator
269  super_ptrtype & operatorlinear(int q)
270  {
271  int q_max = M_operators1.size();
272  FEELPP_ASSERT( q < q_max )( q_max )( q ).error( "OperatorLinearComposite has not enough elements" );
273  return M_operators1.template at(q);
274  }
275 
276  super_ptrtype & operatorlinear(int q, int m)
277  {
278  auto tuple = boost::make_tuple(q,m);
279  return M_operators2.template at(tuple);
280  }
281 
282 
283  virtual void
284  apply( const domain_element_type& de,
285  image_element_type& ie ) const
286  {
287 
288  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
289  sumAllMatrices( matrix, true );
290  matrix->close();
291 
292  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
293  *_v1 = de;_v1->close();
294  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
295  M_backend->prod( matrix, _v1, _v2 );
296  ie.container() = *_v2;
297  }
298 
299 
300  virtual double
301  energy( const typename domain_space_type::element_type & de,
302  const typename dual_image_space_type::element_type & ie ) const
303  {
304 
305  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
306  sumAllMatrices( matrix, true );
307  matrix->close();
308 
309  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
310  *_v1 = de;_v1->close();
311  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
312  *_v2 = ie;
313  vector_ptrtype _v3( M_backend->newVector( _test=ie.functionSpace() ) );
314  M_backend->prod( matrix, _v1, _v3 );
315  return inner_product( _v2, _v3 );
316  }
317 
318 
319  virtual void
320  apply( const typename domain_space_type::element_type & de,
321  typename dual_image_space_type::element_type & ie )
322  {
323  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
324  sumAllMatrices( matrix, true );
325  matrix->close();
326 
327  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
328  *_v1 = de;_v1->close();
329  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
330  M_backend->prod( matrix, _v1, _v2 );
331  ie.container() = *_v2;
332  }
333 
334 
335  virtual void
336  apply( const domain_element_range_type & de,
337  typename dual_image_space_type::element_type & ie )
338  {
339  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
340  sumAllMatrices(matrix, true );
341  matrix->close();
342 
343  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
344  *_v1 = de;_v1->close();
345  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
346  M_backend->prod( matrix, _v1, _v2 );
347  ie.container() = *_v2;
348  }
349 
350  virtual void
351  apply( const typename domain_space_type::element_type & de,
352  dual_image_element_range_type & ie )
353  {
354  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
355  sumAllMatrices( matrix, true );
356  matrix->close();
357 
358  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
359  *_v1 = de;_v1->close();
360  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
361  M_backend->prod( matrix, _v1, _v2 );
362  ie.container() = *_v2;
363  }
364 
365  virtual void
366  apply( const domain_element_range_type & de,
367  dual_image_element_range_type & ie )
368  {
369  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
370  sumAllMatrices( matrix, true );
371  matrix->close();
372 
373  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
374  *_v1 = de;_v1->close();
375  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
376  M_backend->prod( matrix, _v1, _v2 );
377  ie.container() = *_v2;
378  }
379 
380  virtual void
381  apply( const domain_element_slice_type & de,
382  typename dual_image_space_type::element_type & ie )
383  {
384 
385  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
386  sumAllMatrices( matrix, true );
387  matrix->close();
388 
389  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
390  *_v1 = de;_v1->close();
391  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
392  M_backend->prod( matrix, _v1, _v2 );
393  ie.container() = *_v2;
394  }
395 
396 
397  virtual void
398  apply( const typename domain_space_type::element_type & de,
399  dual_image_element_slice_type & ie )
400  {
401  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
402  sumAllMatrices( matrix, true );
403  matrix->close();
404 
405  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
406  *_v1 = de;_v1->close();
407  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
408  M_backend->prod( matrix, _v1, _v2 );
409  ie.container() = *_v2;
410  }
411 
412  virtual void
413  apply( /*const*/ domain_element_slice_type /*&*/ de,
414  dual_image_element_slice_type /*&*/ ie )
415  {
416  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
417  sumAllMatrices( matrix, true );
418  matrix->close();
419 
420  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
421  *_v1 = de;_v1->close();
422  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
423  M_backend->prod( matrix, _v1, _v2 );
424  ie.container() = *_v2;
425  }
426 
427 
428 
429  virtual void
430  apply( const domain_element_range_type & de,
431  dual_image_element_slice_type & ie )
432  {
433  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
434  sumAllMatrices( matrix, true );
435  matrix->close();
436 
437  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
438  *_v1 = de;_v1->close();
439  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
440  M_backend->prod( matrix, _v1, _v2 );
441  ie.container() = *_v2;
442  }
443 
444  virtual void
445  apply( const domain_element_slice_type & de,
446  dual_image_element_range_type & ie )
447  {
448  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
449  sumAllMatrices( matrix, true );
450  matrix->close();
451 
452  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
453  *_v1 = de;_v1->close();
454  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
455  M_backend->prod( matrix, _v1, _v2 );
456  ie.container() = *_v2;
457  }
458 
460  virtual void
461  applyInverse( domain_element_type& de,
462  const image_element_type& ie )
463  {
464  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
465  sumAllMatrices( matrix, true );
466  matrix->close();
467 
468  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
469  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
470  *_v2 = ie.container();
471  M_backend->solve( matrix, matrix, _v1, _v2 );
472  de = *_v1;
473  }
474 
475  this_type& operator=( this_type const& m )
476  {
477  M_backend = m.M_backend;
478  M_pattern = m.M_pattern;
479  M_operators1 = m.M_operators1;
480  M_operators2 = m.M_operators2;
481  M_scalars1 = m.M_scalars1;
482  M_scalars2 = m.M_scalars2;
483  return *this;
484  }
485 
486 
487 
488 private :
489 
490  std::map< int, super_ptrtype > M_operators1;
491  std::map< boost::tuple<int,int> , super_ptrtype > M_operators2;
492  backend_ptrtype M_backend;
493  size_type M_pattern;
494  std::vector< double > M_scalars1;
495  std::vector< std::vector<double> > M_scalars2;
496 };//class OperatorLinearComposite
497 
498 
499 namespace detail
500 {
501 
502 template<typename Args>
503 struct compute_opLinearComposite_return
504 {
505  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
506  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
507 
508  typedef OperatorLinearComposite<domain_space_type, image_space_type> type;
509  typedef boost::shared_ptr<OperatorLinearComposite<domain_space_type, image_space_type> > ptrtype;
510 };
511 }
512 
513 BOOST_PARAMETER_FUNCTION(
514  ( typename Feel::detail::compute_opLinearComposite_return<Args>::ptrtype ), // 1. return type
515  opLinearComposite, // 2. name of the function template
516  tag, // 3. namespace of tag types
517  ( required
518  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
519  ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
520  ) // required
521  ( optional
522  ( backend, *, Backend<typename Feel::detail::compute_opLinearComposite_return<Args>::domain_space_type::value_type>::build() )
523  ( pattern, *, (size_type)Pattern::COUPLED )
524  ) // optionnal
525 )
526 {
527 
528  Feel::detail::ignore_unused_variable_warning( args );
529  typedef typename Feel::detail::compute_opLinearComposite_return<Args>::type oplincomposite_type;
530  typedef typename Feel::detail::compute_opLinearComposite_return<Args>::ptrtype oplincomposite_ptrtype;
531  return oplincomposite_ptrtype ( new oplincomposite_type( domainSpace,imageSpace,backend,pattern) );
532 
533 } // opLinearComposite
534 
535 }//namespace Feel
536 #endif
537 

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