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
operatorlinear.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): Christoph Winkelmann <christoph.winkelmann@epfl.ch>
6  Date: 2006-11-16
7 
8  Copyright (C) 2006 EPFL
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 3.0 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 _OPERATORLINEAR_HPP_
31 #define _OPERATORLINEAR_HPP_
32 
34 #include <feel/feelvf/vf.hpp>
35 
36 namespace Feel
37 {
42 template<class DomainSpace, class DualImageSpace>
43 class OperatorLinear : public Operator<DomainSpace, DualImageSpace>
44 {
46 public:
47 
48  // -- TYPEDEFS --
51 
52  typedef typename super::domain_space_type domain_space_type;
53  typedef typename super::dual_image_space_type dual_image_space_type;
54  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
55  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
56  typedef typename domain_space_type::element_type domain_element_type;
57  typedef typename dual_image_space_type::element_type dual_image_element_type;
58 
59  typedef typename super::backend_type backend_type;
60  typedef typename super::backend_ptrtype backend_ptrtype;
62  typedef typename backend_type::vector_type vector_type;
63  typedef typename backend_type::vector_ptrtype vector_ptrtype;
64  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
65 
66  template<typename T, typename Storage>
67  struct domain_element: public super::domain_space_type::template Element<T,Storage> {};
68 
69  typedef typename domain_space_type::template Element<typename domain_space_type::value_type,
70  typename VectorUblas<typename domain_space_type::value_type>::range::type > domain_element_range_type;
71 
72  typedef typename domain_space_type::template Element<typename domain_space_type::value_type,
73  typename VectorUblas<typename domain_space_type::value_type>::slice::type > domain_element_slice_type;
74 
75  typedef typename dual_image_space_type::template Element<typename dual_image_space_type::value_type,
76  typename VectorUblas<typename dual_image_space_type::value_type>::range::type > dual_image_element_range_type;
77 
78  typedef typename dual_image_space_type::template Element<typename dual_image_space_type::value_type,
79  typename VectorUblas<typename dual_image_space_type::value_type>::slice::type > dual_image_element_slice_type;
80 
81  //template<typename T, typename Storage>
82  //struct dual_image_element: public super::dual_image_space_type::template Element<T,Storage> {};
83 
84  typedef FsFunctionalLinear<DualImageSpace> image_element_type;
85 
87  :
88  super_type(),
89  M_backend( backend_type::build( BACKEND_PETSC ) ),
90  M_matrix(),
91  M_pattern( Pattern::COUPLED ),
92  M_name("operatorlinear")
93  {}
94  OperatorLinear( OperatorLinear const& ol, bool deep_copy = false )
95  :
96  super_type( ol ),
97  M_backend( ol.M_backend ),
98  M_matrix(),
99  M_pattern( ol.M_pattern ),
100  M_name(ol.M_name )
101  {
102  if ( deep_copy )
103  {
104  M_matrix = new matrix_type( *ol.M_matrix );
105  }
106 
107  else
108  {
109  M_matrix = ol.M_matrix;
110  }
111  }
112  OperatorLinear( domain_space_ptrtype domainSpace,
113  dual_image_space_ptrtype dualImageSpace,
114  backend_ptrtype backend,
115  bool buildMatrix = true ,
116  size_type pattern=Pattern::COUPLED)
117  :
118  super_type( domainSpace, dualImageSpace ),
119  M_backend( backend ),
120  M_pattern( pattern ),
121  M_name( "operatorlinear" )
122  {
123  if ( buildMatrix ) M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace , _pattern=M_pattern );
124  }
125 
126  ~OperatorLinear() {}
127 
128  virtual void
129  init( domain_space_ptrtype domainSpace,
130  dual_image_space_ptrtype dualImageSpace,
131  backend_ptrtype backend,
132  bool buildMatrix = true ,
133  size_type pattern = Pattern::COUPLED )
134  {
135  this->setDomainSpace( domainSpace );
136  this->setDualImageSpace( dualImageSpace );
137  M_backend = backend;
138  M_pattern = pattern;
139  if ( buildMatrix ) M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace , _pattern=M_pattern );
140  }
141 
142  void setName( std::string name ) { M_name = name; }
143  std::string name() const { return M_name ; }
144 
145 
146  // apply the operator: ie := Op de
147  template<typename Storage>
148  void
149  apply( const domain_element<typename domain_element_type::value_type, Storage>& de,
150  image_element_type& ie ) const
151  {
152  if ( ! M_matrix->closed() )
153  {
154  M_matrix->close();
155  }
156 
157  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
158  *_v1 = de;_v1->close();
159  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
160  M_backend->prod( M_matrix, _v1, _v2 );
161  ie.container() = *_v2;
162  }
163 
164  virtual void
165  apply( const domain_element_type& de,
166  image_element_type& ie ) const
167  {
168  if ( ! M_matrix->closed() )
169  {
170  M_matrix->close();
171  }
172 
173  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
174  *_v1 = de;_v1->close();
175  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
176  M_backend->prod( M_matrix, _v1, _v2 );
177  ie.container() = *_v2;
178  }
179 
180  virtual double
181  energy( const typename domain_space_type::element_type & de,
182  const typename dual_image_space_type::element_type & ie ) const
183  {
184  if ( ! M_matrix->closed() )
185  {
186  M_matrix->close();
187  }
188 
189  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
190  *_v1 = de;_v1->close();
191  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
192  *_v2 = ie;
193  vector_ptrtype _v3( M_backend->newVector( _test=ie.functionSpace() ) );
194  M_backend->prod( M_matrix, _v1, _v3 );
195  return inner_product( _v2, _v3 );
196  }
197 
198 
199  virtual void
200  apply( const typename domain_space_type::element_type & de,
201  typename dual_image_space_type::element_type & ie )
202  {
203  if ( ! M_matrix->closed() )
204  {
205  M_matrix->close();
206  }
207 
208  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
209  *_v1 = de;_v1->close();
210  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
211  M_backend->prod( M_matrix, _v1, _v2 );
212  ie.container() = *_v2;
213  }
214 
215 
216  virtual void
217  //apply( const typename domain_space_type::template Element<typename domain_space_type::value_type,
218  // typename VectorUblas<typename domain_space_type::value_type>::range::type > & de,
219  apply( const domain_element_range_type & de,
220  typename dual_image_space_type::element_type & ie )
221  {
222  if ( ! M_matrix->closed() )
223  {
224  M_matrix->close();
225  }
226 
227  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
228  *_v1 = de;_v1->close();
229  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
230  M_backend->prod( M_matrix, _v1, _v2 );
231  ie.container() = *_v2;
232  }
233 
234  virtual void
235  apply( const typename domain_space_type::element_type & de,
236  dual_image_element_range_type & ie )
237  {
238  if ( ! M_matrix->closed() )
239  {
240  M_matrix->close();
241  }
242 
243  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
244  *_v1 = de;_v1->close();
245  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
246  M_backend->prod( M_matrix, _v1, _v2 );
247  ie.container() = *_v2;
248  }
249 
250 
251  virtual void
252  apply( const domain_element_range_type & de,
253  dual_image_element_range_type & ie )
254  {
255  if ( ! M_matrix->closed() )
256  {
257  M_matrix->close();
258  }
259 
260  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
261  *_v1 = de;_v1->close();
262  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
263  M_backend->prod( M_matrix, _v1, _v2 );
264  ie.container() = *_v2;
265  }
266 
267 
268  virtual void
269  apply( const domain_element_slice_type & de,
270  typename dual_image_space_type::element_type & ie )
271  {
272  if ( ! M_matrix->closed() )
273  {
274  M_matrix->close();
275  }
276 
277  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
278  *_v1 = de;_v1->close();
279  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
280  M_backend->prod( M_matrix, _v1, _v2 );
281  ie.container() = *_v2;
282  }
283 
284  virtual void
285  apply( const typename domain_space_type::element_type & de,
286  dual_image_element_slice_type & ie )
287  {
288  if ( ! M_matrix->closed() )
289  {
290  M_matrix->close();
291  }
292 
293  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
294  *_v1 = de;_v1->close();
295  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
296  M_backend->prod( M_matrix, _v1, _v2 );
297  ie.container() = *_v2;
298  }
299 
300 
301  virtual void
302  apply( /*const*/ domain_element_slice_type /*&*/ de,
303  dual_image_element_slice_type /*&*/ ie )
304  {
305  if ( ! M_matrix->closed() )
306  {
307  M_matrix->close();
308  }
309 
310  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
311  *_v1 = de;_v1->close();
312  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
313  M_backend->prod( M_matrix, _v1, _v2 );
314  ie.container() = *_v2;
315  }
316 #if 0
317  void
318  apply( const typename domain_space_type::element_type::component_type & de,
319  typename dual_image_space_type::element_type & ie )
320  {
321  if ( ! M_matrix->closed() )
322  {
323  M_matrix->close();
324  }
325 
326  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
327  *_v1 = de;_v1->close();
328  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
329  M_backend->prod( M_matrix, _v1, _v2 );
330  ie.container() = *_v2;
331  }
332 #endif
333 
334 
335 
336 
337  virtual void
338  apply( const domain_element_range_type & de,
339  dual_image_element_slice_type & ie )
340  {
341  if ( ! M_matrix->closed() )
342  {
343  M_matrix->close();
344  }
345 
346  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
347  *_v1 = de;_v1->close();
348  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
349  M_backend->prod( M_matrix, _v1, _v2 );
350  ie.container() = *_v2;
351  }
352 
353  virtual void
354  apply( const domain_element_slice_type & de,
355  dual_image_element_range_type & ie )
356  {
357  if ( ! M_matrix->closed() )
358  {
359  M_matrix->close();
360  }
361 
362  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
363  *_v1 = de;_v1->close();
364  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
365  M_backend->prod( M_matrix, _v1, _v2 );
366  ie.container() = *_v2;
367  }
368 
369 
370 
371 
372 #if 0
373  template<typename Storage1,typename Storage2>
374  void
375  apply( const domain_element<typename domain_element_type::value_type, Storage1>& de,
376  dual_image_element<typename dual_image_element_type::value_type, Storage2>& ie ) //const
377  {
378  if ( ! M_matrix->closed() )
379  {
380  M_matrix->close();
381  }
382 
383  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
384  *_v1 = de;_v1->close();
385  //vector_ptrtype _v2( M_backend->newVector( ie.space()->dof() ) );
386  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
387  M_backend->prod( M_matrix, _v1, _v2 );
388  ie.container() = *_v2;
389  }
390 #endif
391  template <typename T1 = typename domain_space_type::element_type,
392  typename T2 = typename dual_image_space_type::element_type >
393  T2
394  operator()( T1 & de )
395  {
396  T2 elt_image( this->dualImageSpace(),"oio" );
397  this->apply( de,elt_image );
398 
399  return elt_image;
400  }
401 
402  template <typename T1 = typename domain_space_type::element_type,
403  typename T2 = typename dual_image_space_type::element_type >
404  T2
405  operator()( boost::shared_ptr<T1> de )
406  {
407  T2 elt_image( this->dualImageSpace(),"oio" );
408  this->apply( *de,elt_image );
409 
410  return elt_image;
411  }
412 
413 
415  virtual void
416  applyInverse( domain_element_type& de,
417  const image_element_type& ie )
418  {
419  if ( ! M_matrix->closed() )
420  {
421  M_matrix->close();
422  }
423 
424  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
425  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
426  *_v2 = ie.container();
427  M_backend->solve( M_matrix, M_matrix, _v1, _v2 );
428  de = *_v1;
429 
430 #if 0
431 
432  if ( !M_backend->converged() )
433  {
434  std::cerr << "[OperatorLinear::applyInverse] "
435  << "solver failed to converge" << std::endl;
436  }
437 
438 #endif
439  }
440 
441 
442  // "improved version of applyInverse" :
443  // take directly the expression of the rhs
444  // return the result of inversion
445  template<typename RhsExpr>
446  domain_element_type
447  inv( RhsExpr const& rhs_expr )
448  {
449  if ( ! M_matrix->closed() )
450  {
451  M_matrix->close();
452  }
453 
454  domain_element_type de = this->domainSpace()->element();
455 
456  auto ie = M_backend->newVector( _test=this->dualImageSpace() );
457  form1( _test=this->dualImageSpace(), _vector=ie ) =
458  integrate( elements( this->domainSpace()->mesh() ),
459  rhs_expr * id( this->dualImageSpace()->element() ) );
460 
461  M_backend->solve( M_matrix, de, ie );
462 
463  return de;
464  }
465 
466 
467  // fill underlying matrix
468  template<class ExprT>
469  this_type& operator=( ExprT const& e )
470  {
471  // M_matrix->clear();
472  form2( _trial=this->domainSpace(),
473  _test=this->dualImageSpace(),
474  _matrix=M_matrix ) = e;
475  return *this;
476  }
477 
478  this_type& operator=( this_type const& m )
479  {
480  M_backend = m.M_backend;
481  M_matrix->zero();
482  M_matrix->addMatrix( 1.0, m.M_matrix );
483  M_pattern = m.M_pattern;
484 
485  return *this;
486  }
487 
488  // add to underlying matrix
489  template<class ExprT>
490  this_type& operator+=( ExprT const& e )
491  {
492  form2( _trial=this->domainSpace(),
493  _test=this->dualImageSpace(),
494  _matrix=M_matrix ) += e;
495  return *this;
496  }
497 
498  void close()
499  {
500  if ( ! M_matrix->closed() )
501  {
502  M_matrix->close();
503  }
504 
505  }
506  // retrieve underlying matrix
507  matrix_type& mat()
508  {
509  return *M_matrix;
510  }
511 
512  // retrieve underlying matrix
513  matrix_type const& mat() const
514  {
515  return *M_matrix;
516  }
517 
518  // retrieve underlying matrix
519  matrix_ptrtype const& matPtr() const
520  {
521  return M_matrix;
522  }
523 
524  // retrieve underlying matrix
525  matrix_ptrtype& matPtr()
526  {
527  return M_matrix;
528  }
529 
530  virtual void matPtr( matrix_ptrtype & matrix )
531  {
532  matrix = M_matrix;
533  }
534 
535  template<typename T>
536  OperatorLinear& add( T const& scalar, OperatorLinear const& ol )
537  {
538  this->close();
539  M_matrix->addMatrix( scalar, *ol.M_matrix );
540  return *this;
541  }
542 
543  backend_ptrtype& backend()
544  {
545  return M_backend;
546  }
547 
548  size_type pattern()
549  {
550  return M_pattern;
551  }
552 
553  template<typename T>
554  OperatorLinear& add( T const& scalar, boost::shared_ptr<OperatorLinear> ol )
555  {
556  this->close();
557  this->add( scalar, *ol );
558  return *this;
559  }
560 
561  OperatorLinear& operator+( boost::shared_ptr<OperatorLinear> ol )
562  {
563  this->close();
564  this->add( 1.0, *ol );
565  return *this;
566  }
567  OperatorLinear& operator+( OperatorLinear const& ol )
568  {
569  this->close();
570  this->add( 1.0, ol );
571  return *this;
572  }
573 
574 private:
575 
576  backend_ptrtype M_backend;
577  matrix_ptrtype M_matrix;
578  size_type M_pattern;
579  std::string M_name;
580 }; // class Operator
581 
582 
583 
584 template<typename Args>
585 struct compute_opLinear_return
586 {
587  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
588  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
589 
590  typedef boost::shared_ptr<OperatorLinear<domain_space_type, image_space_type> > type;
591 };
592 
593 BOOST_PARAMETER_FUNCTION(
594  ( typename compute_opLinear_return<Args>::type ), // 1. return type
595  opLinear, // 2. name of the function template
596  tag, // 3. namespace of tag types
597  ( required
598  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
599  ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
600  ) // required
601  ( optional
602  ( backend, *, Backend<typename compute_opLinear_return<Args>::domain_space_type::value_type>::build() )
603  ( pattern, *, (size_type)Pattern::COUPLED )
604  ) // optionnal
605 )
606 {
607  Feel::detail::ignore_unused_variable_warning( args );
608  typedef OperatorLinear<typename compute_opLinear_return<Args>::domain_space_type,
609  typename compute_opLinear_return<Args>::image_space_type> operatorlinear_type;
610 
611  boost::shared_ptr<operatorlinear_type> opI( new operatorlinear_type( domainSpace,imageSpace,backend,pattern ) );
612 
613  return opI;
614 
615 } // opLinear
616 
617 
618 
619 
620 
621 } // Feel
622 
623 #endif /* _OPERATORLINEAR_HPP_ */

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