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
operatorlinearfree.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 -*-
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 __OperatorLinearFree_H
31 #define __OperatorLinearFree_H 1
32 
33 #include <feel/feel.hpp>
34 #include <boost/ref.hpp>
35 #include <boost/next_prior.hpp>
36 #include <boost/type_traits.hpp>
37 #include <boost/tuple/tuple.hpp>
38 #if BOOST_VERSION >= 104700
39 #include <boost/math/special_functions/nonfinite_num_facets.hpp>
40 #endif
41 #include <feel/feelvf/vf.hpp>
42 namespace Feel
43 {
44 
50 template<class DomainSpace, class DualImageSpace, typename ExprType>
51 class OperatorLinearFree : public OperatorLinear<DomainSpace, DualImageSpace>
52 {
54 public:
55 
56  typedef ExprType expr_type;
57 
58  typedef typename super::domain_space_type domain_space_type;
59  typedef typename super::dual_image_space_type dual_image_space_type;
60  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
61  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
62 
63  typedef typename super::backend_type backend_type;
64  typedef typename super::backend_ptrtype backend_ptrtype;
65 
66  typedef typename super::matrix_type matrix_type;
67  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
68 
69  typedef FsFunctionalLinear<DualImageSpace> image_element_type;
70  typedef typename super::domain_element_type domain_element_type;
71  typedef typename super::dual_image_element_type dual_image_element_type;
72 
73  typedef typename super::domain_element_range_type domain_element_range_type;
74  typedef typename super::domain_element_slice_type domain_element_slice_type;
75  typedef typename super::dual_image_element_range_type dual_image_element_range_type;
76  typedef typename super::dual_image_element_slice_type dual_image_element_slice_type;
77 
79  typedef boost::shared_ptr<this_type> this_ptrtype;
80 
82  :
83  super()
84  {}
85 
86  OperatorLinearFree (domain_space_ptrtype domainSpace,
87  dual_image_space_ptrtype dualImageSpace,
88  backend_ptrtype backend,
89  expr_type expr,
90  size_type pattern=Pattern::COUPLED )
91  :
92  super( domainSpace , dualImageSpace , backend , false , pattern ),
93  M_backend( backend ),
94  M_expr( expr ),
95  M_pattern( pattern )
96  {}
97 
98 
99  virtual void
100  init( domain_space_ptrtype domainSpace,
101  dual_image_space_ptrtype dualImageSpace,
102  backend_ptrtype backend,
103  bool buildMatrix = false ,
104  size_type pattern = Pattern::COUPLED )
105  {
106  this->setDomainSpace( domainSpace );
107  this->setDualImageSpace( dualImageSpace );
108  M_backend = backend;
109  M_pattern = pattern;
110  }
111 
112 
113  expr_type expr()
114  {
115  return M_expr;
116  }
117 
118  virtual void
119  apply( const domain_element_type& de,
120  image_element_type& ie ) const
121  {
122 
123  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
124  _trial=this->domainSpace(),
125  _pattern=M_pattern );
126 
127 
128  form2( _trial=this->domainSpace(),
129  _test=this->dualImageSpace(),
130  _matrix=matrix,
131  _pattern=M_pattern) = M_expr;
132 
133  matrix->close();
134 
135  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
136  *_v1 = de;_v1->close();
137  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
138  M_backend->prod( matrix, _v1, _v2 );
139  ie.container() = *_v2;
140  }
141 
142 
143  virtual double
144  energy( const typename domain_space_type::element_type & de,
145  const typename dual_image_space_type::element_type & ie ) const
146  {
147 
148  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
149  _trial=this->domainSpace(),
150  _pattern=M_pattern );
151 
152  form2( _trial=this->domainSpace(),
153  _test=this->dualImageSpace(),
154  _matrix=matrix,
155  _pattern=M_pattern) = M_expr;
156 
157  matrix->close();
158 
159  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
160  *_v1 = de;_v1->close();
161  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
162  *_v2 = ie;
163  vector_ptrtype _v3( M_backend->newVector( _test=ie.functionSpace() ) );
164  M_backend->prod( matrix, _v1, _v3 );
165  return inner_product( _v2, _v3 );
166  }
167 
168 
169  virtual void
170  apply( const typename domain_space_type::element_type & de,
171  typename dual_image_space_type::element_type & ie )
172  {
173  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
174  _trial=this->domainSpace(),
175  _pattern=M_pattern );
176 
177  form2( _trial=this->domainSpace(),
178  _test=this->dualImageSpace(),
179  _matrix=matrix,
180  _pattern=M_pattern) = M_expr;
181 
182  matrix->close();
183 
184  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
185  *_v1 = de;_v1->close();
186  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
187  M_backend->prod( matrix, _v1, _v2 );
188  ie.container() = *_v2;
189  }
190 
191 
192  virtual void
193  apply( const domain_element_range_type & de,
194  typename dual_image_space_type::element_type & ie )
195  {
196  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
197  _trial=this->domainSpace(),
198  _pattern=M_pattern );
199 
200  form2( _trial=this->domainSpace(),
201  _test=this->dualImageSpace(),
202  _matrix=matrix,
203  _pattern=M_pattern) = M_expr;
204 
205  matrix->close();
206 
207  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
208  *_v1 = de;_v1->close();
209  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
210  M_backend->prod( matrix, _v1, _v2 );
211  ie.container() = *_v2;
212  }
213 
214  virtual void
215  apply( const typename domain_space_type::element_type & de,
216  dual_image_element_range_type & ie )
217  {
218 
219  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
220  _trial=this->domainSpace(),
221  _pattern=M_pattern );
222 
223  form2( _trial=this->domainSpace(),
224  _test=this->dualImageSpace(),
225  _matrix=matrix,
226  _pattern=M_pattern) = M_expr;
227 
228  matrix->close();
229 
230  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
231  *_v1 = de;_v1->close();
232  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
233  M_backend->prod( matrix, _v1, _v2 );
234  ie.container() = *_v2;
235  }
236 
237  virtual void
238  apply( const domain_element_range_type & de,
239  dual_image_element_range_type & ie )
240  {
241 
242  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
243  _trial=this->domainSpace(),
244  _pattern=M_pattern );
245 
246  form2( _trial=this->domainSpace(),
247  _test=this->dualImageSpace(),
248  _matrix=matrix,
249  _pattern=M_pattern) = M_expr;
250 
251  matrix->close();
252 
253  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
254  *_v1 = de;_v1->close();
255  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
256  M_backend->prod( matrix, _v1, _v2 );
257  ie.container() = *_v2;
258  }
259 
260  virtual void
261  apply( const domain_element_slice_type & de,
262  typename dual_image_space_type::element_type & ie )
263  {
264  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
265  _trial=this->domainSpace(),
266  _pattern=M_pattern );
267 
268  form2( _trial=this->domainSpace(),
269  _test=this->dualImageSpace(),
270  _matrix=matrix,
271  _pattern=M_pattern) = M_expr;
272 
273  matrix->close();
274 
275  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
276  *_v1 = de;_v1->close();
277  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
278  M_backend->prod( matrix, _v1, _v2 );
279  ie.container() = *_v2;
280  }
281 
282 
283  virtual void
284  apply( const typename domain_space_type::element_type & de,
285  dual_image_element_slice_type & ie )
286  {
287  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
288  _trial=this->domainSpace(),
289  _pattern=M_pattern );
290 
291  form2( _trial=this->domainSpace(),
292  _test=this->dualImageSpace(),
293  _matrix=matrix,
294  _pattern=M_pattern) = M_expr;
295 
296  matrix->close();
297 
298  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
299  *_v1 = de;_v1->close();
300  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
301  M_backend->prod( matrix, _v1, _v2 );
302  ie.container() = *_v2;
303  }
304 
305  virtual void
306  apply( /*const*/ domain_element_slice_type /*&*/ de,
307  dual_image_element_slice_type /*&*/ ie )
308  {
309  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
310  _trial=this->domainSpace(),
311  _pattern=M_pattern );
312 
313  form2( _trial=this->domainSpace(),
314  _test=this->dualImageSpace(),
315  _matrix=matrix,
316  _pattern=M_pattern) = M_expr;
317 
318  matrix->close();
319 
320  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
321  *_v1 = de;_v1->close();
322  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
323  M_backend->prod( matrix, _v1, _v2 );
324  ie.container() = *_v2;
325  }
326 
327 
328 
329  virtual void
330  apply( const domain_element_range_type & de,
331  dual_image_element_slice_type & ie )
332  {
333  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
334  _trial=this->domainSpace(),
335  _pattern=M_pattern );
336 
337  form2( _trial=this->domainSpace(),
338  _test=this->dualImageSpace(),
339  _matrix=matrix,
340  _pattern=M_pattern) = M_expr;
341 
342  matrix->close();
343 
344  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
345  *_v1 = de;_v1->close();
346  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
347  M_backend->prod( matrix, _v1, _v2 );
348  ie.container() = *_v2;
349  }
350 
351  virtual void
352  apply( const domain_element_slice_type & de,
353  dual_image_element_range_type & ie )
354  {
355  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
356  _trial=this->domainSpace(),
357  _pattern=M_pattern );
358 
359  form2( _trial=this->domainSpace(),
360  _test=this->dualImageSpace(),
361  _matrix=matrix,
362  _pattern=M_pattern) = M_expr;
363 
364  matrix->close();
365 
366  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
367  *_v1 = de;_v1->close();
368  vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
369  M_backend->prod( matrix, _v1, _v2 );
370  ie.container() = *_v2;
371  }
372 
374  virtual void
375  applyInverse( domain_element_type& de,
376  const image_element_type& ie )
377  {
378  auto matrix = M_backend->newMatrix( _test=this->dualImageSpace() ,
379  _trial=this->domainSpace(),
380  _pattern=M_pattern );
381 
382  form2( _trial=this->domainSpace(),
383  _test=this->dualImageSpace(),
384  _matrix=matrix,
385  _pattern=M_pattern) = M_expr;
386 
387  matrix->close();
388 
389  vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
390  vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
391  *_v2 = ie.container();
392  M_backend->solve( matrix, matrix, _v1, _v2 );
393  de = *_v1;
394 
395  }
396 
397 
398  this_type& operator=( this_type const& m )
399  {
400  M_backend = m.M_backend;
401  M_expr = m.M_expr;
402  M_pattern = m.M_pattern;
403  return *this;
404  }
405 
406 
407  // retrieve underlying matrix
408  virtual void matPtr( matrix_ptrtype & matrix )
409  {
410  form2( _trial=this->domainSpace(),
411  _test=this->dualImageSpace(),
412  _matrix=matrix,
413  _pattern=M_pattern) = M_expr;
414 
415  matrix->close();
416  }
417 
418 private:
419  backend_ptrtype M_backend;
420  expr_type M_expr;
421  size_type M_pattern;
422 };//OperatorLinearFree
423 
424 
425 namespace detail
426 {
427 
428 template<typename Args>
429 struct compute_opLinearFree_return
430 {
431  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
432  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
433  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::expr>::type>::type expr_type;
434 
435  typedef OperatorLinearFree<domain_space_type, image_space_type, expr_type> type;
436  typedef boost::shared_ptr<OperatorLinearFree<domain_space_type, image_space_type, expr_type> > ptrtype;
437 };
438 }
439 
440 BOOST_PARAMETER_FUNCTION(
441  ( typename Feel::detail::compute_opLinearFree_return<Args>::ptrtype ), // 1. return type
442  opLinearFree, // 2. name of the function template
443  tag, // 3. namespace of tag types
444  ( required
445  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
446  ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
447  ( expr , * )
448  ) // required
449  ( optional
450  ( backend, *, Backend<typename Feel::detail::compute_opLinearFree_return<Args>::domain_space_type::value_type>::build() )
451  ( pattern, *, (size_type)Pattern::COUPLED )
452  ) // optionnal
453 )
454 {
455 
456  Feel::detail::ignore_unused_variable_warning( args );
457  typedef typename Feel::detail::compute_opLinearFree_return<Args>::type oplinfree_type;
458  typedef typename Feel::detail::compute_opLinearFree_return<Args>::ptrtype oplinfree_ptrtype;
459  return oplinfree_ptrtype ( new oplinfree_type( domainSpace,imageSpace,backend,expr ) );
460 
461 } // opLinearFree
462 
463 }//namespace Feel
464 #endif
Linear Operator between function spaces, represented by a matrix.
Definition: operatorlinear.hpp:43
base class for all linear algebra backends
Definition: backend.hpp:133
type_traits< T >::real_type inner_product(Vector< T > const &v1, Vector< T > const &v2)
Definition: vector.hpp:579
virtual void applyInverse(domain_element_type &de, const image_element_type &ie)
apply the inverse of the operator:
Definition: operatorlinearfree.hpp:375
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Definition: matrixsparse.hpp:76
Operator between function spaces.
Definition: operator.hpp:44
Linear Operator Free between function spaces, represented by a matrix but this matrix is not automati...
Definition: operatorlinearfree.hpp:51

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