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
projectors.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-05-31
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006-2011 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 __Projectors_H
31 #define __Projectors_H 1
32 
33 #include <boost/timer.hpp>
34 #include <feel/feelcore/parameter.hpp>
35 #include <feel/feeldiscr/functionspace.hpp>
36 
37 namespace Feel
38 {
39 namespace vf
40 {
41 namespace details
42 {
50 template<ProjectorType iDim, typename FunctionSpaceType, typename IteratorRange, typename ExprT>
51 class Projector
52 {
53 public:
54 
55 
59 
60  static const size_type context = ExprT::context|vm::POINT;
61 
62  typedef FunctionSpaceType functionspace_type;
63  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
64  typedef typename functionspace_type::element_type element_type;
65  typedef typename functionspace_type::basis_type basis_type;
66  typedef ExprT expression_type;
67  typedef typename expression_type::value_type value_type;
68 
69 
70  static const uint16_type imorder = functionspace_type::basis_type::nOrder;
71  static const bool imIsPoly = true;
72 
73  typedef typename boost::tuples::template element<0, IteratorRange>::type idim_type;
74  typedef typename boost::tuples::template element<1, IteratorRange>::type iterator_type;
75  typedef IteratorRange range_iterator;
76 
77 
79 
83 
84  Projector( functionspace_ptrtype const& __functionspace,
85  IteratorRange const& r,
86  expression_type const& __expr,
87  GeomapStrategyType geomap_strategy )
88  :
89  M_functionspace( __functionspace ),
90  M_range( r ),
91  M_expr( __expr ),
92  M_geomap_strategy( geomap_strategy )
93  {
94  DVLOG(2) << "Projector constructor from expression\n";
95  }
96 
97 
98  Projector( Projector const& __vfi )
99  :
100  M_functionspace( __vfi.M_functionspace ),
101  M_range( __vfi.M_range ),
102  M_expr( __vfi.M_expr ),
103  M_geomap_strategy( __vfi.M_geomap_strategy )
104  {
105  DVLOG(2) << "Projector copy constructor\n";
106  }
107 
108  virtual ~Projector() {}
109 
111 
115 
116  element_type operator()( const bool sum = false ) const
117  {
118  return operator()( sum, idim_type() );
119  }
120 
122 
126 
133  expression_type const& expression() const
134  {
135  return M_expr;
136  }
137 
139 
143 
145 
149 
151 
152 private:
153 
154  element_type operator()( const bool sum, mpl::size_t<MESH_ELEMENTS> ) const;
155  element_type operator()( const bool sum, mpl::size_t<MESH_FACES> ) const;
156  element_type operator()( const bool sum, mpl::size_t<MESH_POINTS> ) const;
157 
158 private:
159 
160  functionspace_ptrtype const& M_functionspace;
161  range_iterator M_range;
162  expression_type const& M_expr;
163  GeomapStrategyType M_geomap_strategy;
164 };
165 
166 template<ProjectorType iDim, typename FunctionSpaceType, typename Iterator, typename ExprT>
168 Projector<iDim, FunctionSpaceType, Iterator, ExprT>::operator()( const bool sum, mpl::size_t<MESH_ELEMENTS> ) const
169 {
170  boost::timer __timer;
171 
172  element_type __v( M_functionspace );
173  FEELPP_ASSERT( __v.size() == M_functionspace->dof()->nDof() )( __v.size() )( M_functionspace->dof()->nDof() ).warn( "invalid size" );
174  __v.setZero();
175 
176  typedef typename functionspace_type::fe_type fe_type;
177  fe_type* __fe = __v.functionSpace()->fe().get();
178 
179  typedef typename functionspace_type::mesh_type mesh_type;
180  typedef typename mesh_type::element_type element_type;
181 
182  typedef typename element_type::gm_type gm_type;
183  typedef typename gm_type::template Context<context, element_type> gm_context_type;
184  typedef typename element_type::gm1_type gm1_type;
185  typedef typename gm1_type::template Context<context, element_type> gm1_context_type;
186 
187 
188  typedef typename fe_type::template Context<context, fe_type, gm_type, element_type> fecontext_type;
189  typedef typename fe_type::template Context<context, fe_type, gm1_type, element_type> fecontext1_type;
190 
191  typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
192  typedef boost::shared_ptr<gm1_context_type> gm1_context_ptrtype;
193  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gm_context_ptrtype> > map_gmc_type;
194  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gm1_context_ptrtype> > map_gmc1_type;
195  //typedef typename expression_type::template tensor<map_gmc_type,fusion::map<fusion::pair<vf::detail::gmc<0>,boost::shared_ptr<fecontext_type> > > > t_expr_type;
196  //typedef decltype( basis_type::isomorphism( M_expr ) ) the_expression_type;
197  typedef expression_type the_expression_type;
198  typedef typename boost::remove_reference<typename boost::remove_const<the_expression_type>::type >::type iso_expression_type;
199  typedef typename iso_expression_type::template tensor<map_gmc_type> t_expr_type;
200  typedef typename iso_expression_type::template tensor<map_gmc1_type> t_expr1_type;
201  typedef typename t_expr_type::value_type value_type;
202 
203  // we should manipulate the same type of functions on the left and
204  // on the right
205  //BOOST_STATIC_ASSERT(( boost::is_same<return_value_type, typename functionspace_type::return_value_type>::value ));
206 
207  //
208  // Precompute some data in the reference element for
209  // geometric mapping and reference finite element
210  //
211  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( __v.functionSpace()->gm(),
212  __fe->points() ) );
213  typename gm1_type::precompute_ptrtype __geopc1( new typename gm1_type::precompute_type( __v.mesh()->gm1(),
214  __fe->points() ) );
215 
216 
217 
218  const uint16_type ndofv = functionspace_type::fe_type::nDof;
219 
220  iterator_type it, en;
221  boost::tie( boost::tuples::ignore, it, en ) = M_range;
222 
223  // return if no elements
224  if ( it == en )
225  return __v;
226 
227  gm_context_ptrtype __c( new gm_context_type( __v.functionSpace()->gm(),*it,__geopc ) );
228  gm1_context_ptrtype __c1( new gm1_context_type( __v.mesh()->gm1(),*it,__geopc1 ) );
229 
230  typedef typename t_expr_type::shape shape;
231  static const bool is_rank_ok = ( shape::M == FunctionSpaceType::nComponents1 &&
232  shape::N == FunctionSpaceType::nComponents2 );
233 
234  BOOST_MPL_ASSERT_MSG( mpl::bool_<is_rank_ok>::value,
235  INVALID_TENSOR_RANK,
236  ( mpl::int_<shape::M>, mpl::int_<FunctionSpaceType::nComponents>, shape ) );
237 
238  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
239  t_expr_type tensor_expr( basis_type::isomorphism( M_expr ), mapgmc );
240 
241  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
242 
243  t_expr1_type tensor_expr1( basis_type::isomorphism( M_expr ), mapgmc1 );
244 
245  std::vector<bool> points_done( __v.functionSpace()->dof()->nLocalDof()/__v.nComponents );
246  std::fill( points_done.begin(), points_done.end(),false );
247 
248  for ( ; it!=en ; ++it )
249  {
250  switch ( M_geomap_strategy )
251  {
252  case GeomapStrategyType::GEOMAP_HO:
253  {
254  __c->update( *it );
255  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
256  tensor_expr.update( mapgmc );
257 
258  for ( uint16_type __j = 0; __j < ndofv; ++__j )
259  {
260  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
261  //for ( uint16_type c2 = 0; c2 < shape::N;++c2 )
262  {
263  if ( sum )
264  __v.plus_assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
265 
266  else
267  __v.assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
268  }
269  }
270  }
271  break;
272 
273  case GeomapStrategyType::GEOMAP_O1:
274  {
275  __c1->update( *it );
276  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
277  tensor_expr1.update( mapgmc1 );
278 
279  for ( uint16_type __j = 0; __j < ndofv; ++__j )
280  {
281  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
282  //for ( uint16_type c2 = 0; c2 < shape::N;++c2 )
283  {
284  if ( sum )
285  __v.plus_assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
286 
287  else
288  __v.assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
289  }
290  }
291  }
292  break;
293 
294  case GeomapStrategyType::GEOMAP_OPT:
295  {
296  if ( it->isOnBoundary() )
297  {
298  // HO if on boundary
299  __c->update( *it );
300  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
301  tensor_expr.update( mapgmc );
302 
303  for ( uint16_type __j = 0; __j < ndofv; ++__j )
304  {
305  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
306  //for ( uint16_type c2 = 0; c2 < shape::N;++c2 )
307  {
308  if ( sum )
309  __v.plus_assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
310 
311  else
312  __v.assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
313  }
314  }
315  }
316 
317  else
318  {
319  __c1->update( *it );
320  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
321  tensor_expr1.update( mapgmc1 );
322 
323  for ( uint16_type __j = 0; __j < ndofv; ++__j )
324  {
325  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
326  //for ( uint16_type c2 = 0; c2 < shape::N;++c2 )
327  {
328  if ( sum )
329  __v.plus_assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
330 
331  else
332  __v.assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
333  }
334  }
335  }
336  }
337  break;
338  }
339  }
340 
341  return __v;
342 }
343 
344 template<ProjectorType iDim, typename FunctionSpaceType, typename Iterator, typename ExprT>
345 typename Projector<iDim, FunctionSpaceType, Iterator, ExprT>::element_type
346 Projector<iDim, FunctionSpaceType, Iterator, ExprT>::operator()( const bool sum, mpl::size_t<MESH_FACES> ) const
347 {
348  boost::timer __timer;
349 
350  element_type __v( M_functionspace );
351  __v.setZero();
352 
353  DVLOG(2) << "call project(MESH_FACES) " << "\n";
354  //
355  // a few typedefs
356  //
357 
358  // mesh element
359  typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
360  typedef typename geoelement_type::face_type face_type;
361 
362  // geometric mapping context
363  typedef typename geoelement_type::gm_type gm_type;
364  typedef boost::shared_ptr<gm_type> gm_ptrtype;
365  typedef typename geoelement_type::gm1_type gm1_type;
366  typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
367 
368  typedef typename gm_type::template Context<context, geoelement_type> gmc_type;
369  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
370  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
371  typedef typename gm1_type::template Context<context, geoelement_type> gmc1_type;
372  typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
373  typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
374 
375 
376  // dof
377  typedef typename element_type::functionspace_type::dof_type dof_type;
378 
379  // basis
380  typedef typename element_type::functionspace_type::fe_type fe_type;
381  typedef typename fe_type::template Context< context, fe_type, gm_type, geoelement_type> fecontext_type;
382  typedef boost::shared_ptr<fecontext_type> fecontext_ptrtype;
383  typedef typename fe_type::template Context< context, fe_type, gm1_type, geoelement_type> fecontext1_type;
384  typedef boost::shared_ptr<fecontext1_type> fecontext1_ptrtype;
385  //typedef fusion::map<fusion::pair<vf::detail::gmc<0>, fecontext_ptrtype> > map_gmc_type;
386 
387  // expression
388  //typedef typename expression_type::template tensor<map_gmc_type,fecontext_type> t_expr_type;
389  //typedef decltype( basis_type::isomorphism( M_expr ) ) the_expression_type;
390  typedef expression_type the_expression_type;
391  typedef typename boost::remove_reference<typename boost::remove_const<the_expression_type>::type >::type iso_expression_type;
392  typedef typename iso_expression_type::template tensor<map_gmc_type> t_expr_type;
393  typedef typename iso_expression_type::template tensor<map_gmc1_type> t_expr1_type;
394  typedef typename t_expr_type::shape shape;
395 
396  //
397  // start
398  //
399  DVLOG(2) << "assembling Dirichlet conditions\n";
400 
401  dof_type const* __dof = __v.functionSpace()->dof().get();
402 
403  fe_type const* __fe = __v.functionSpace()->fe().get();
404 
405  iterator_type __face_it, __face_en;
406  boost::tie( boost::tuples::ignore, __face_it, __face_en ) = M_range;
407 
408  if ( __face_it == __face_en )
409  return __v;
410 
411  gm_ptrtype __gm( new gm_type );
412  gm1_ptrtype __gm1( new gm1_type );
413 
414 
415 
416  //
417  // Precompute some data in the reference element for
418  // geometric mapping and reference finite element
419  //
420  typedef typename geoelement_type::permutation_type permutation_type;
421  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
422  typedef typename gm_type::precompute_type geopc_type;
423  typedef typename gm1_type::precompute_ptrtype geopc1_ptrtype;
424  typedef typename gm1_type::precompute_type geopc1_type;
425 
426  DVLOG(2) << "[integratoron] numTopologicalFaces = " << geoelement_type::numTopologicalFaces << "\n";
427  std::vector<std::map<permutation_type, geopc_ptrtype> > __geopc( geoelement_type::numTopologicalFaces );
428  std::vector<std::map<permutation_type, geopc1_ptrtype> > __geopc1( geoelement_type::numTopologicalFaces );
429 
430  for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
431  {
432  for ( permutation_type __p( permutation_type::IDENTITY );
433  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
434  {
435  __geopc[__f][__p] = geopc_ptrtype( new geopc_type( __gm, __fe->points( __f ) ) );
436  __geopc1[__f][__p] = geopc1_ptrtype( new geopc1_type( __gm1, __fe->points( __f ) ) );
437  //DVLOG(2) << "[geopc] FACE_ID = " << __f << " ref pts=" << __fe->dual().points( __f ) << "\n";
438  FEELPP_ASSERT( __geopc[__f][__p]->nPoints() ).error( "invalid number of points for geopc" );
439  FEELPP_ASSERT( __geopc1[__f][__p]->nPoints() ).error( "invalid number of points for geopc1" );
440  }
441  }
442 
443  uint16_type __face_id = __face_it->pos_first();
444  gmc_ptrtype __c( new gmc_type( __gm, __face_it->element( 0 ), __geopc, __face_id ) );
445  gmc1_ptrtype __c1( new gmc1_type( __gm1, __face_it->element( 0 ), __geopc1, __face_id ) );
446 
447  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
448  t_expr_type expr( basis_type::isomorphism( M_expr ), mapgmc );
449  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
450  t_expr1_type expr1( basis_type::isomorphism( M_expr ), mapgmc1 );
451 
452 
453 
454 
456 
457  if ( !fe_type::is_modal )
458  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
459  face_type::numEdges * fe_type::nDofPerEdge +
460  face_type::numFaces * fe_type::nDofPerFace );
461 
462  else
463  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
464 
465  DVLOG(2) << "[projector::operator(MESH_FACES)] nbFaceDof = " << nbFaceDof << "\n";
466 
467  std::vector<int> dofs;
468  std::vector<value_type> values;
469 
470  for ( ; __face_it != __face_en; ++__face_it )
471  {
472  FEELPP_ASSERT( __face_it->isOnBoundary() && !__face_it->isConnectedTo1() )
473  ( __face_it->marker() )
474  ( __face_it->isOnBoundary() )
475  ( __face_it->ad_first() )
476  ( __face_it->pos_first() )
477  ( __face_it->ad_second() )
478  ( __face_it->pos_second() )
479  ( __face_it->id() ).warn( "inconsistent data face" );
480  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id()
481  << " element id= " << __face_it->ad_first()
482  << " pos in elt= " << __face_it->pos_first()
483  << " marker: " << __face_it->marker() << "\n";
484  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id() << " real pts=" << __face_it->G() << "\n";
485 
486  uint16_type __face_id = __face_it->pos_first();
487 
488  std::pair<size_type,size_type> range_dof( std::make_pair( __v.start(),
489  __v.functionSpace()->nDof() ) );
490  DVLOG(2) << "[projector] dof start = " << range_dof.first << "\n";
491  DVLOG(2) << "[projector] dof range = " << range_dof.second << "\n";
492 
493  switch ( M_geomap_strategy )
494  {
495  default:
496  case GeomapStrategyType::GEOMAP_OPT:
497  case GeomapStrategyType::GEOMAP_HO:
498  {
499  __c->update( __face_it->element( 0 ), __face_id );
500  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id() << " ref pts=" << __c->xRefs() << "\n";
501  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id() << " real pts=" << __c->xReal() << "\n";
502 
503  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
504 
505  expr.update( mapgmc );
506 
507  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
508  for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
509  {
510  for ( uint16_type l = 0; l < nbFaceDof; ++l )
511  {
512  typename expression_type::value_type __value = expr.evalq( c1, c2, l );
513 
514  size_type thedof = __v.start() +
515  boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
516 
517  if ( sum )
518  __v( thedof ) += __value;
519 
520  else
521  __v( thedof ) = __value;
522  }
523  }
524  }
525  break;
526 
527  case GeomapStrategyType::GEOMAP_O1:
528  {
529  __c1->update( __face_it->element( 0 ), __face_id );
530  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id() << " ref pts=" << __c1->xRefs() << "\n";
531  DVLOG(2) << "[projector] FACE_ID = " << __face_it->id() << " real pts=" << __c1->xReal() << "\n";
532 
533  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
534 
535  expr1.update( mapgmc1 );
536 
537 
538  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
539  for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
540  {
541  for ( uint16_type l = 0; l < nbFaceDof; ++l )
542  {
543  typename expression_type::value_type __value = expr1.evalq( c1, c2, l );
544 
545  size_type thedof = __v.start() +
546  boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
547 
548  if ( sum )
549  __v( thedof ) += __value;
550 
551  else
552  __v( thedof ) = __value;
553  }
554  }
555  }
556  break;
557  }
558 
559  } // face_it
560 
561  return __v;
562 }
563 
564 template<ProjectorType iDim, typename FunctionSpaceType, typename Iterator, typename ExprT>
565 typename Projector<iDim, FunctionSpaceType, Iterator, ExprT>::element_type
566 Projector<iDim, FunctionSpaceType, Iterator, ExprT>::operator()( const bool sum, mpl::size_t<MESH_POINTS> ) const
567 {
568  boost::timer __timer;
569 
570  element_type __v( M_functionspace );
571  __v.setZero();
572 
573  iterator_type pt_it, pt_en;
574  boost::tie( boost::tuples::ignore, pt_it, pt_en ) = M_range;
575 
576  if ( pt_it == pt_en )
577  return __v;
578 #if 0
579  BOOST_FOREACH( auto dof, M_functionspace->dof()->markerToDof( pt_it->marker() ) );
580  {
581 
582  // get the first element to which the point/dof belong and then build
583  // the proper geomap context in order to evaluate the expression at the
584  // point
585 
586 #if 0
587  if ( sum )
588  __v.add( dof.second, );
589 #endif
590  }
591 #endif
592  return __v;
593 }
594 
595 } // detail
597 
598 
604 template<typename FunctionSpaceType, typename IteratorRange, typename ExprT>
605 typename FunctionSpaceType::element_type
606 project( boost::shared_ptr<FunctionSpaceType> const& __functionspace,
607  IteratorRange const& range_it,
608  Expr<ExprT> const& __expr,
609  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
610 {
611  typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
612  proj_t p( __functionspace, range_it, __expr, geomap );
613  return p();
614 }
615 
616 template<typename FunctionSpaceType, typename IteratorRange, typename ExprT>
617 typename FunctionSpaceType::element_type
618 project_impl( boost::shared_ptr<FunctionSpaceType> const& __functionspace,
619  IteratorRange const& range_it,
620  Expr<ExprT> const& __expr,
621  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
622 {
623  typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
624  proj_t p( __functionspace, range_it, __expr, geomap );
625  return p();
626 }
627 
632 template<typename FunctionSpaceType, typename ExprT>
633 typename FunctionSpaceType::element_type
634 project( boost::shared_ptr<FunctionSpaceType> const& __functionspace, Expr<ExprT> const& __expr,
635  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
636 {
637  return project( __functionspace, elements( __functionspace->mesh() ), __expr, geomap );
638 }
639 
645 template<typename FunctionSpaceType, typename IteratorRange, typename ExprT>
646 typename FunctionSpaceType::element_type
647 sum( boost::shared_ptr<FunctionSpaceType> const& __functionspace,
648  IteratorRange const& range_it,
649  Expr<ExprT> const& __expr,
650  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
651 {
652  typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
653  proj_t p( __functionspace, range_it, __expr, geomap );
654  return p( true );
655 }
656 
661 template<typename FunctionSpaceType, typename ExprT>
662 typename FunctionSpaceType::element_type
663 sum( boost::shared_ptr<FunctionSpaceType> const& __functionspace, Expr<ExprT> const& __expr )
664 {
665  return sum( __functionspace, elements( __functionspace->mesh() ), __expr );
666 }
667 
673 template<typename FunctionSpaceType, typename ExprT>
674 typename FunctionSpaceType::element_type
675 project( FunctionSpaceType const& __functionspace, Expr<ExprT> const& __expr,
676  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
677 {
678  typedef __typeof__( __functionspace->mesh()->elementsRange() ) IteratorRange;
679  typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
680  proj_t p( __functionspace, __functionspace->mesh()->elementsRange(), __expr, geomap );
681  return p();
682 }
683 
689 template<typename FunctionSpaceType, typename IteratorRange, typename ExprT>
690 typename FunctionSpaceType::element_type
691 project( FunctionSpaceType const& __functionspace,
692  IteratorRange const& range_it,
693  Expr<ExprT> const& __expr,
694  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
695 {
696  typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
697  proj_t p( __functionspace, range_it, __expr, geomap );
698  return p();
699 }
700 
702 namespace detail
703 {
704 template<typename S>
705 struct space_ptr
706 {
707  typedef typename S::element_type type;
708 };
709 
710 template<typename S>
711 struct space_value
712 {
713  typedef S type;
714 };
715 
716 template<typename Args>
717 struct project
718 {
719  typedef typename clean_type<Args,tag::space>::type the_space_type;
720  typedef typename mpl::if_<is_shared_ptr<the_space_type>,
721  mpl::identity<space_ptr<the_space_type> >,
722  mpl::identity<space_value<the_space_type> > >::type::type space_type;
723  typedef typename space_type::type _space_type;
724  typedef boost::shared_ptr<_space_type> _space_ptrtype;
725  typedef typename _space_type::element_type element_type;
726  //typedef typename clean_type<Args,tag::expr>::type _expr_type;
727  //typedef typename clean_type<Args,tag::range>::type _range_type;
728 };
729 }
731 
742 BOOST_PARAMETER_FUNCTION(
743  ( typename vf::detail::project<Args>::element_type ), // return type
744  project, // 2. function name
745 
746  tag, // 3. namespace of tag types
747 
748  ( required
749  ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<Feel::FunctionSpaceBase> > ) )
750  ( expr, * )
751  ) // 4. one required parameter, and
752 
753  ( optional
754  ( range, *, elements( space->mesh() ) )
755  ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
756  ( accumulate, *( boost::is_integral<mpl::_> ), false )
757  )
758 )
759 {
760 #if 0
761  typedef typename vf::detail::project<Args>::_space_type _space_type;
762  typedef typename vf::detail::project<Args>::_range_type _range_type;
763  typedef typename vf::detail::project<Args>::_expr_type _expr_type;
764 
765  typedef details::Projector<NODAL, _space_type, _range_type, Expr<_expr_type> > proj_t;
766  proj_t p( space, range, expr,geomap );
767  return p( sum );
768 #else
769  // if ( accumulate )
770  //return sum( space, range, expr, geomap );
771  return project_impl( space, range, expr, geomap );
772 #endif
773 }
774 
775 } // vf
776 } // feel
777 
778 
779 #endif /* __Projectors_H */

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