30 #ifndef __Projectors_H
31 #define __Projectors_H 1
33 #include <boost/timer.hpp>
34 #include <feel/feelcore/parameter.hpp>
35 #include <feel/feeldiscr/functionspace.hpp>
50 template<ProjectorType iDim,
typename FunctionSpaceType,
typename IteratorRange,
typename ExprT>
60 static const size_type context = ExprT::context|vm::POINT;
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;
70 static const uint16_type imorder = functionspace_type::basis_type::nOrder;
71 static const bool imIsPoly =
true;
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;
84 Projector( functionspace_ptrtype
const& __functionspace,
85 IteratorRange
const& r,
86 expression_type
const& __expr,
87 GeomapStrategyType geomap_strategy )
89 M_functionspace( __functionspace ),
92 M_geomap_strategy( geomap_strategy )
94 DVLOG(2) <<
"Projector constructor from expression\n";
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 )
105 DVLOG(2) <<
"Projector copy constructor\n";
116 element_type operator()(
const bool sum =
false )
const
118 return operator()( sum, idim_type() );
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;
160 functionspace_ptrtype
const& M_functionspace;
161 range_iterator M_range;
162 expression_type
const& M_expr;
163 GeomapStrategyType M_geomap_strategy;
166 template<ProjectorType iDim,
typename FunctionSpaceType,
typename Iterator,
typename ExprT>
170 boost::timer __timer;
172 element_type __v( M_functionspace );
173 FEELPP_ASSERT( __v.size() == M_functionspace->dof()->nDof() )( __v.size() )( M_functionspace->dof()->nDof() ).warn(
"invalid size" );
176 typedef typename functionspace_type::fe_type fe_type;
177 fe_type* __fe = __v.functionSpace()->fe().get();
179 typedef typename functionspace_type::mesh_type mesh_type;
180 typedef typename mesh_type::element_type element_type;
182 typedef typename element_type::gm_type gm_type;
184 typedef typename element_type::gm1_type gm1_type;
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;
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;
211 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( __v.functionSpace()->gm(),
213 typename gm1_type::precompute_ptrtype __geopc1(
new typename gm1_type::precompute_type( __v.mesh()->gm1(),
218 const uint16_type ndofv = functionspace_type::fe_type::nDof;
220 iterator_type it, en;
221 boost::tie( boost::tuples::ignore, it, en ) = M_range;
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 ) );
230 typedef typename t_expr_type::shape shape;
231 static const bool is_rank_ok = ( shape::M == FunctionSpaceType::nComponents1 &&
232 shape::N == FunctionSpaceType::nComponents2 );
234 BOOST_MPL_ASSERT_MSG( mpl::bool_<is_rank_ok>::value,
236 ( mpl::int_<shape::M>, mpl::int_<FunctionSpaceType::nComponents>, shape ) );
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 );
241 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
243 t_expr1_type tensor_expr1( basis_type::isomorphism( M_expr ), mapgmc1 );
245 std::vector<bool> points_done( __v.functionSpace()->dof()->nLocalDof()/__v.nComponents );
246 std::fill( points_done.begin(), points_done.end(),false );
248 for ( ; it!=en ; ++it )
250 switch ( M_geomap_strategy )
252 case GeomapStrategyType::GEOMAP_HO:
255 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
256 tensor_expr.update( mapgmc );
258 for ( uint16_type __j = 0; __j < ndofv; ++__j )
260 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
264 __v.plus_assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
267 __v.assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
273 case GeomapStrategyType::GEOMAP_O1:
276 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
277 tensor_expr1.update( mapgmc1 );
279 for ( uint16_type __j = 0; __j < ndofv; ++__j )
281 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
285 __v.plus_assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
288 __v.assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
294 case GeomapStrategyType::GEOMAP_OPT:
296 if ( it->isOnBoundary() )
300 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
301 tensor_expr.update( mapgmc );
303 for ( uint16_type __j = 0; __j < ndofv; ++__j )
305 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
309 __v.plus_assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
312 __v.assign( it->id(), __j, c1, tensor_expr.evalq( c1, 0, __j ) );
320 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
321 tensor_expr1.update( mapgmc1 );
323 for ( uint16_type __j = 0; __j < ndofv; ++__j )
325 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
329 __v.plus_assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
332 __v.assign( it->id(), __j, c1, tensor_expr1.evalq( c1, 0, __j ) );
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
348 boost::timer __timer;
350 element_type __v( M_functionspace );
353 DVLOG(2) <<
"call project(MESH_FACES) " <<
"\n";
359 typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
360 typedef typename geoelement_type::face_type face_type;
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;
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;
377 typedef typename element_type::functionspace_type::dof_type dof_type;
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;
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;
399 DVLOG(2) <<
"assembling Dirichlet conditions\n";
401 dof_type
const* __dof = __v.functionSpace()->dof().get();
403 fe_type
const* __fe = __v.functionSpace()->fe().get();
405 iterator_type __face_it, __face_en;
406 boost::tie( boost::tuples::ignore, __face_it, __face_en ) = M_range;
408 if ( __face_it == __face_en )
411 gm_ptrtype __gm(
new gm_type );
412 gm1_ptrtype __gm1(
new gm1_type );
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;
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 );
430 for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
432 for ( permutation_type __p( permutation_type::IDENTITY );
433 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
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 ) ) );
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" );
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 ) );
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 );
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 );
463 nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
465 DVLOG(2) <<
"[projector::operator(MESH_FACES)] nbFaceDof = " << nbFaceDof <<
"\n";
467 std::vector<int> dofs;
468 std::vector<value_type> values;
470 for ( ; __face_it != __face_en; ++__face_it )
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";
486 uint16_type __face_id = __face_it->pos_first();
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";
493 switch ( M_geomap_strategy )
496 case GeomapStrategyType::GEOMAP_OPT:
497 case GeomapStrategyType::GEOMAP_HO:
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";
503 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
505 expr.update( mapgmc );
507 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
508 for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
510 for ( uint16_type l = 0; l < nbFaceDof; ++l )
512 typename expression_type::value_type __value = expr.evalq( c1, c2, l );
515 boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
518 __v( thedof ) += __value;
521 __v( thedof ) = __value;
527 case GeomapStrategyType::GEOMAP_O1:
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";
533 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
535 expr1.update( mapgmc1 );
538 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
539 for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
541 for ( uint16_type l = 0; l < nbFaceDof; ++l )
543 typename expression_type::value_type __value = expr1.evalq( c1, c2, l );
546 boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
549 __v( thedof ) += __value;
552 __v( thedof ) = __value;
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
568 boost::timer __timer;
570 element_type __v( M_functionspace );
573 iterator_type pt_it, pt_en;
574 boost::tie( boost::tuples::ignore, pt_it, pt_en ) = M_range;
576 if ( pt_it == pt_en )
579 BOOST_FOREACH(
auto dof, M_functionspace->dof()->markerToDof( pt_it->marker() ) );
588 __v.add( dof.second, );
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 )
611 typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
612 proj_t p( __functionspace, range_it, __expr, geomap );
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 )
623 typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
624 proj_t p( __functionspace, range_it, __expr, geomap );
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 )
637 return project( __functionspace,
elements( __functionspace->mesh() ), __expr, geomap );
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 )
652 typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
653 proj_t p( __functionspace, range_it, __expr, geomap );
661 template<
typename FunctionSpaceType,
typename ExprT>
662 typename FunctionSpaceType::element_type
663 sum( boost::shared_ptr<FunctionSpaceType>
const& __functionspace, Expr<ExprT>
const& __expr )
665 return sum( __functionspace,
elements( __functionspace->mesh() ), __expr );
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 )
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 );
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 )
696 typedef details::Projector<NODAL, FunctionSpaceType, IteratorRange, Expr<ExprT> > proj_t;
697 proj_t p( __functionspace, range_it, __expr, geomap );
707 typedef typename S::element_type type;
716 template<
typename Args>
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;
742 BOOST_PARAMETER_FUNCTION(
743 (
typename vf::detail::project<Args>::element_type ),
749 ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<Feel::FunctionSpaceBase> > ) )
754 ( range, *,
elements( space->mesh() ) )
755 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
756 ( accumulate, *( boost::is_integral<mpl::_> ), false )
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;
765 typedef details::Projector<NODAL, _space_type, _range_type, Expr<_expr_type> > proj_t;
766 proj_t p( space, range, expr,geomap );
771 return project_impl( space, range, expr, geomap );