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
projector.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): Vincent Doyeux <vincent.doyeux@ujf-grenoble.fr>
6  Date: 2011-04-25
7 
8  Copyright (C) 2011 Université Joseph Fourier (Grenoble I)
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 */
29 #ifndef _PROJECTOR_HPP_
30 #define _PROJECTOR_HPP_
31 
32 #include <feel/feelcore/parameter.hpp>
34 #include <feel/feelvf/vf.hpp>
35 
36 namespace Feel
37 {
38  template<class DomainSpace, class DualImageSpace> class Projector;
39 
40 namespace detail
41 {
42 template<typename Args>
43 struct projector_args
44 {
45  typedef typename vf::detail::clean_type<Args,tag::domainSpace>::type::element_type domain_type;
46  typedef typename vf::detail::clean_type<Args,tag::imageSpace>::type::element_type image_type;
47  typedef boost::shared_ptr<Projector<domain_type,image_type> > return_type;
48 };
49 
50 template<typename Args>
51 struct lift_args
52 {
53  typedef typename vf::detail::clean_type<Args,tag::domainSpace>::type::element_type domain_type;
54  typedef boost::shared_ptr<Projector<domain_type,domain_type> > lift_return_type;
55 };
56 
57 } // detail
65 template<class DomainSpace, class DualImageSpace>
66 class Projector : public OperatorLinear<DomainSpace, DualImageSpace>
67 {
68  typedef Projector<DomainSpace,DualImageSpace> super;
69 
70 public :
71 
75 
76 
77  typedef OperatorLinear<DomainSpace, DualImageSpace> ol_type;
78 
79  typedef typename super::domain_space_type domain_space_type;
80  typedef typename super::dual_image_space_type dual_image_space_type;
81  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
82  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
83  typedef typename domain_space_type::element_type domain_element_type;
84  typedef typename dual_image_space_type::element_type dual_image_element_type;
85 
86  typedef typename super::backend_type backend_type;
87  typedef typename super::backend_ptrtype backend_ptrtype;
88  typedef typename backend_type::sparse_matrix_type matrix_type;
89  typedef typename backend_type::vector_type vector_type;
90  typedef typename backend_type::vector_ptrtype vector_ptrtype;
91  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
92 
93  typedef FsFunctionalLinear<DualImageSpace> image_element_type;
94 
96 
99 
100  Projector( domain_space_ptrtype domainSpace,
101  dual_image_space_ptrtype dualImageSpace,
102  backend_ptrtype backend = Backend<double>::build( BACKEND_PETSC ),
103  ProjectorType proj_type=L2,
104  double epsilon = 0.01,
105  double gamma = 20,
106  DirichletType dirichlet_type = WEAK
107  )
108  :
109  ol_type( domainSpace, dualImageSpace, backend ),
110  M_backend( backend ),
111  M_epsilon( epsilon ),
112  M_gamma( gamma ),
113  M_proj_type( proj_type ),
114  M_dir( dirichlet_type )
115 
116  {
117  M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace, _pattern=Pattern::EXTENDED ) ;
118  initMatrix();
119  }
120 
121  ~Projector() {}
123 
127  template<typename Args,typename IntEltsDefault>
128  struct integrate_type
129  {
130  typedef typename vf::detail::clean_type<Args,tag::expr>::type _expr_type;
131  typedef typename vf::detail::clean2_type<Args,tag::range,IntEltsDefault>::type _range_type;
132  typedef typename vf::detail::clean2_type<Args,tag::quad, _Q< vf::ExpressionOrder<_range_type,_expr_type>::value > >::type _quad_type;
133  typedef typename vf::detail::clean2_type<Args,tag::quad1, _Q< vf::ExpressionOrder<_range_type,_expr_type>::value_1 > >::type _quad1_type;
134  };
135 
136  template<typename Range, typename Expr>
137  void applyOn( Range range, Expr expr)
138  {
139  typedef typename boost::tuples::template element<0, Range>::type idim_type;
140  applyOn( range, expr, mpl::int_<idim_type::value>() );
141  }
142 
143 
144  BOOST_PARAMETER_MEMBER_FUNCTION( ( domain_element_type ),
145  project,
146  tag,
147  ( required
148  ( expr, * )
149  )
150  ( optional
151  ( range, *, elements( this->dualImageSpace()->mesh() ) )
152  ( quad, *, ( typename integrate_type<Args,decltype( elements( this->dualImageSpace()->mesh() ) )>::_quad_type() ) )
153  ( quad1, *, ( typename integrate_type<Args,decltype( elements( this->dualImageSpace()->mesh() ) )>::_quad1_type() ) )
154  ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
155  )
156  )
157  {
158  using namespace vf;
159 
160  auto sol = this->domainSpace()->element();
161 
162  ie = M_backend->newVector( this->dualImageSpace() );
163 
164  form1( _test=this->dualImageSpace(), _vector=ie, _init=true );
165 
166  if ( (M_proj_type != LIFT) )
167  {
168  form1( _test=this->dualImageSpace(), _vector=ie ) +=
169  integrate( _range=range, _expr=expr * id( this->dualImageSpace()->element() ),
170  _quad=quad, _quad1=quad1, _geomap=geomap );
171  }
172  else if ( ( M_proj_type == LIFT ) && ( M_dir == WEAK ) )
173  {
174  form1( _test=this->dualImageSpace(), _vector=ie ) +=
175  integrate( _range=range,
176  _expr=expr*( -grad( this->dualImageSpace()->element() )*vf::N() +
177  M_gamma / vf::hFace() *id( this->dualImageSpace()->element() ) ),
178  _quad=quad, _quad1=quad1, _geomap=geomap );
179  }
180 
181  //weak boundary conditions
182  if ( M_proj_type == DIFF )
183  {
184  form1( _test=this->dualImageSpace(), _vector=ie ) +=
185  integrate( _range=boundaryfaces( this->dualImageSpace()->mesh() ),
186  _expr=expr*M_epsilon*( -grad( this->domainSpace()->element() )*vf::N() +
187  M_gamma / vf::hFace() *id( this->dualImageSpace()->element() ) ),
188  _quad=quad, _quad1=quad1, _geomap=geomap );
189  }
190 
191  ie->close();
192 
193  M_matrixFull = M_backend->newMatrix( _trial=this->domainSpace(), _test=this->dualImageSpace(), _pattern=Pattern::EXTENDED );
194  auto bilinearForm = form2( _trial=this->domainSpace(), _test=this->dualImageSpace(), _matrix=M_matrixFull );
195 
196  if ( ( M_proj_type == LIFT ) && ( M_dir == WEAK ) )
197  {
198  bilinearForm +=
199  integrate( _range=range, _expr=
200  ( -trans( id( this->dualImageSpace()->element() ) )*gradt( this->domainSpace()->element() )*vf::N()
201  -trans( idt( this->domainSpace()->element() ) )* grad( this->dualImageSpace()->element() )*vf::N()
202  + M_gamma * trans( idt( this->domainSpace()->element() ) ) /*trial*/
203  *id( this->dualImageSpace()->element() ) / vf::hFace() /*test*/
204  ), _quad=quad, _quad1=quad1, _geomap=geomap );
205  }
206 
207  M_matrixFull->close();
208  M_matrixFull->addMatrix( 1., M_matrix );
209 
210  if ( ( M_proj_type == LIFT ) && ( M_dir == STRONG ) )
211  this->applyOn(range, expr);
212 
213  M_backend->solve( M_matrixFull, sol, ie );
214 
215  return sol;
216  }
217 
218  template<typename RhsExpr>
219  domain_element_type
220  operator()( RhsExpr const& rhs_expr )
221  {
222  return this->project( rhs_expr );
223  }
224 
225  template<typename RhsExpr>
226  void
227  operator()( domain_element_type& de, RhsExpr const& rhs_expr )
228  {
229  de = this->project( rhs_expr );
230  }
231 
232  domain_element_type
233  operator()( image_element_type const& ie )
234  {
235  domain_element_type de = this->domainSpace()->element();
236  M_backend->solve( M_matrix, de, ie );
237  return de ;
238  }
239 
240  void
241  operator()( domain_element_type &de, image_element_type const& ie )
242  {
243  M_backend->solve( M_matrix, de, ie );
244  }
245 
246  template<typename Range, typename Expr>
247  domain_element_type
248  operator()( Range const& range ,Expr const& expr )
249  {
250  return this->project( expr, range );
251  }
252 
253  void
254  apply( domain_element_type& de,
255  image_element_type const& ie )
256  {
257  M_backend->solve( M_matrix, de, ie );
258  }
259 
260  template<typename RhsExpr>
261  void
262  apply( domain_element_type& de,
263  RhsExpr const& rhs_expr )
264  {
265  de=this->project( rhs_expr );
266  }
267 
268 
269  template< typename Expr>
270  domain_element_type derivate( Expr expr )
271  {
272  de = this->domainSpace()->element();
273  ie = M_backend->newVector( this->dualImageSpace() );
274 
275  form1( _test=this->dualImageSpace(), _vector=ie, _init=true );
276 
277  form1( _test=this->dualImageSpace(), _vector=ie ) +=
278  integrate(_range = elements( this->dualImageSpace()->mesh() ),
279  _expr = - trace( expr * trans( grad( this->dualImageSpace()->element() ) ) ) );
280 
281  form1(_test=this->dualImageSpace(), _vector=ie) +=
282  integrate(_range = boundaryfaces( this->dualImageSpace()->mesh() ),
283  _expr = trans( id( this->dualImageSpace()->element() ) ) * expr * vf::N() );
284 
285  ie->close();
286 
287  M_backend->solve( M_matrix, de, ie );
288  return de;
289  }
290 
291 
293 
294 
295 private :
296 
297  void initMatrix()
298  {
299  using namespace vf;
300  auto a = form2 ( _trial=this->domainSpace(),
301  _test=this->dualImageSpace(),
302  _matrix=M_matrix,
303  _init=true );
304 
305  switch ( M_proj_type )
306  {
307  case L2:
308  {
309  a = integrate( elements( this->dualImageSpace()->mesh() ),
310  trans( idt( this->domainSpace()->element() ) ) /*trial*/
311  *id( this->dualImageSpace()->element() ) /*test*/
312  );
313  }
314  break;
315 
316  case H1:
317  {
318  a = integrate( elements( this->dualImageSpace()->mesh() ),
319  trans( idt( this->domainSpace()->element() ) ) /*trial*/
320  *id( this->dualImageSpace()->element() ) /*test*/
321  +
322  trace( gradt( this->domainSpace()->element() )
323  * trans( grad( this->dualImageSpace()->element() ) ) )
324  );
325  }
326  break;
327 
328  case DIFF:
329  {
330  a = integrate( elements( this->dualImageSpace()->mesh() ),
331  trans( idt( this->domainSpace()->element() ) ) /*trial*/
332  *id( this->dualImageSpace()->element() ) /*test*/
333  +
334  M_epsilon *
335  trace( gradt( this->domainSpace()->element() )
336  * trans( grad( this->dualImageSpace()->element() ) ) )
337  );
338  //weak boundary conditions
339  a += integrate( boundaryfaces( this->dualImageSpace()->mesh() ),
340  M_epsilon*( -trans( id( this->dualImageSpace()->element() ) )*gradt( this->domainSpace()->element() )*vf::N() ) );
341  a += integrate( boundaryfaces( this->dualImageSpace()->mesh() ),
342  M_epsilon*( -trans( idt( this->domainSpace()->element() ) )* grad( this->dualImageSpace()->element() )*vf::N() ) );
343  a += integrate( boundaryfaces( this->dualImageSpace()->mesh() ),
344  M_epsilon*( M_gamma * trans( idt( this->domainSpace()->element() ) ) /*trial*/
345  *id( this->dualImageSpace()->element() ) / vf::hFace() /*test*/
346  ) );
347  }
348  break;
349 
350  case HDIV:
351  {
352  a = integrate( elements( this->dualImageSpace()->mesh() ),
353  trans( idt( this->domainSpace()->element() ) ) /*trial*/
354  *id( this->dualImageSpace()->element() ) /*test*/
355  +
356  ( divt( this->domainSpace()->element() ) *
357  div( this->dualImageSpace()->element() ) )
358  );
359  }
360  break;
361 
362  case HCURL:
363  {
364  a = integrate( elements( this->dualImageSpace()->mesh() ),
365  trans( idt( this->domainSpace()->element() ) ) /*trial*/
366  *id( this->dualImageSpace()->element() ) /*test*/
367  +
368  // only for 2D, need to specialize this for 3D
369  curlzt( this->domainSpace()->element() )
370  * curlz( this->dualImageSpace()->element() )
371  );
372  }
373  break;
374 
375  case LIFT:
376  {
377  a = integrate( elements( this->dualImageSpace()->mesh() ),
378  trace( gradt( this->domainSpace()->element() )
379  * trans( grad( this->dualImageSpace()->element() ) ) )
380  );
381 
382  }
383  break;
384 
385  case CIP:
386  {
387  a = integrate( elements( this->dualImageSpace()->mesh() ),
388  trans( idt( this->domainSpace()->element() ) ) /*trial*/
389  *id( this->dualImageSpace()->element() ) /*test*/
390  );
391 
392  a += integrate( internalfaces( this->dualImageSpace()->mesh() ),
393  M_gamma * hFace() * hFace()
394  * trans(jumpt( gradt(this->domainSpace()->element()) ))
395  * jump( grad(this->dualImageSpace()->element()) )
396  );
397  }
398  break;
399 
400  case NODAL:
401  break;
402  }
403 
404  M_matrix->close();
405  }
406 
407  template<typename Range, typename Expr>
408  void applyOn( Range range, Expr expr, mpl::int_<MESH_ELEMENTS> ){}
409 
410  template<typename Range, typename Expr>
411  void applyOn( Range range, Expr expr, mpl::int_<MESH_FACES> )
412  {
413  form2 ( _trial=this->domainSpace(),
414  _test=this->dualImageSpace(),
415  _matrix=M_matrixFull ) += on( _range=range , _element=de, _rhs=ie, _expr=expr );
416  }
417 
418  backend_ptrtype M_backend;
419  const double M_epsilon;
420  const double M_gamma;
421  const ProjectorType M_proj_type;
422  DirichletType M_dir;
423  matrix_ptrtype M_matrix;
424  matrix_ptrtype M_matrixFull;
425  domain_element_type de;
426  vector_ptrtype ie;
427 
428 };//Projector
429 
430 #if 1
431 
439 template<typename TDomainSpace, typename TDualImageSpace>
440 boost::shared_ptr< Projector<TDomainSpace, TDualImageSpace> >
441 projector( boost::shared_ptr<TDomainSpace> const& domainspace,
442  boost::shared_ptr<TDualImageSpace> const& imagespace,
443  typename Projector<TDomainSpace, TDualImageSpace>::backend_ptrtype const& backend = Backend<double>::build( BACKEND_PETSC ),
444  ProjectorType proj_type=L2, double epsilon=0.01, double gamma = 20, DirichletType dirichlet_type = WEAK)
445 {
447  boost::shared_ptr<Proj_type> proj( new Proj_type( domainspace, imagespace, backend, proj_type, epsilon, gamma, dirichlet_type ) );
448  return proj;
449 }
450 
451 BOOST_PARAMETER_FUNCTION( ( typename Feel::detail::projector_args<Args>::return_type ),
452  opProjection,
453  tag,
454  ( required
455  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
456  ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
457  )
458  ( optional
459  ( type, (ProjectorType), L2 )
460  ( penaldir, *( boost::is_arithmetic<mpl::_> ), 20. )
461  ( backend, *, Backend<double>::build( BACKEND_PETSC ) )
462  ) )
463 {
464  return projector( domainSpace,imageSpace, backend, type, 0.01, penaldir );
465 }
466 
467 BOOST_PARAMETER_FUNCTION( ( typename Feel::detail::lift_args<Args>::lift_return_type ),
468  opLift,
469  tag,
470  ( required
471  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
472  )
473  ( optional
474  ( type, (DirichletType), WEAK )
475  ( penaldir, *( boost::is_arithmetic<mpl::_> ), 20. )
476  ( backend, *, Backend<double>::build( BACKEND_PETSC ) )
477  ) )
478 {
479  return projector( domainSpace, domainSpace, backend, ProjectorType::LIFT, 0.01 , penaldir, type );
480 }
481 
482 
483 #endif
484 
485 
486 
487 } //namespace Feel
488 
489 
490 #endif

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