29 #ifndef __FEELPP_EVALUATORS_H
30 #define __FEELPP_EVALUATORS_H 1
32 #include <boost/timer.hpp>
33 #include <feel/feelcore/parameter.hpp>
34 #include <feel/feeldiscr/functionspace.hpp>
53 template<EvaluatorType iDim,
typename IteratorRange,
typename Pset,
typename ExprT>
63 static const size_type context = ExprT::context|vm::POINT;
65 typedef ExprT expression_type;
66 typedef typename expression_type::value_type value_type;
69 static const uint16_type imorder = 1;
70 static const bool imIsPoly =
true;
72 typedef typename boost::tuples::template element<0, IteratorRange>::type idim_type;
73 typedef typename boost::tuples::template element<1, IteratorRange>::type iterator_type;
74 typedef typename iterator_type::value_type mesh_element_type;
75 typedef IteratorRange range_iterator;
76 typedef typename mpl::if_<mpl::bool_<mesh_element_type::is_simplex>,
77 mpl::identity<typename Pset::template apply<mesh_element_type::nDim, value_type, Simplex>::type >,
78 mpl::identity<typename Pset::template apply<mesh_element_type::nDim, value_type, Hypercube>::type >
79 >::type::type pointset_type;
80 typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> element_type;
81 typedef Eigen::Matrix<value_type,mesh_element_type::nDim,Eigen::Dynamic> node_type;
82 typedef boost::tuple<element_type,node_type> eval_element_type;
89 Evaluator( IteratorRange
const& r,
91 expression_type
const& __expr,
92 GeomapStrategyType geomap_strategy )
97 M_geomap_strategy( geomap_strategy )
99 DVLOG(2) <<
"Evaluator constructor from expression\n";
103 Evaluator( Evaluator
const& __vfi )
105 M_range( __vfi.M_range ),
106 M_pset( __vfi.M_pset ),
107 M_expr( __vfi.M_expr ),
108 M_geomap_strategy( __vfi.M_geomap_strategy )
110 DVLOG(2) <<
"Evaluator copy constructor\n";
113 virtual ~Evaluator() {}
121 eval_element_type operator()()
const
123 return operator()( idim_type() );
138 expression_type
const& expression()
const
159 eval_element_type operator()( mpl::size_t<MESH_ELEMENTS> )
const;
160 eval_element_type operator()( mpl::size_t<MESH_FACES> )
const;
164 range_iterator M_range;
165 pointset_type M_pset;
166 expression_type
const& M_expr;
167 GeomapStrategyType M_geomap_strategy;
170 template<EvaluatorType iDim,
typename Iterator,
typename Pset,
typename ExprT>
171 typename Evaluator<iDim, Iterator, Pset, ExprT>::eval_element_type
172 Evaluator<iDim, Iterator, Pset, ExprT>::operator()( mpl::size_t<MESH_ELEMENTS> )
const
174 boost::timer __timer;
176 typedef typename mesh_element_type::gm_type gm_type;
177 typedef typename gm_type::template Context<context, mesh_element_type> gm_context_type;
178 typedef typename mesh_element_type::gm1_type gm1_type;
179 typedef typename gm1_type::template Context<context, mesh_element_type> gm1_context_type;
182 typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
183 typedef boost::shared_ptr<gm1_context_type> gm1_context_ptrtype;
184 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gm_context_ptrtype> > map_gmc_type;
185 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gm1_context_ptrtype> > map_gmc1_type;
188 typedef expression_type the_expression_type;
189 typedef typename boost::remove_reference<typename boost::remove_const<the_expression_type>::type >::type iso_expression_type;
190 typedef typename iso_expression_type::template tensor<map_gmc_type> t_expr_type;
191 typedef typename iso_expression_type::template tensor<map_gmc1_type> t_expr1_type;
192 typedef typename t_expr_type::value_type value_type;
198 typedef typename t_expr_type::shape shape;
201 iterator_type it, en;
202 boost::tie( boost::tuples::ignore, it, en ) = M_range;
204 int npoints = M_pset.points().size2();
205 element_type __v( M_pset.points().size2()*std::distance( it, en )*shape::M );
206 node_type __p( mesh_element_type::nDim, M_pset.points().size2()*std::distance( it, en ) );
214 return boost::make_tuple( __v, __p );
220 typename gm_type::precompute_ptrtype __geopc(
new typename gm_type::precompute_type( it->gm(),
222 typename gm1_type::precompute_ptrtype __geopc1(
new typename gm1_type::precompute_type( it->gm1(),
227 gm_context_ptrtype __c(
new gm_context_type( it->gm(),*it,__geopc ) );
228 gm1_context_ptrtype __c1(
new gm1_context_type( it->gm1(),*it,__geopc1 ) );
230 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
231 t_expr_type tensor_expr( M_expr, mapgmc );
233 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
235 t_expr1_type tensor_expr1( M_expr, mapgmc1 );
237 for (
int e = 0; it!=en ; ++it, ++e )
239 switch ( M_geomap_strategy )
241 case GeomapStrategyType::GEOMAP_HO:
244 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
245 tensor_expr.update( mapgmc );
247 for ( uint16_type p = 0; p < npoints; ++p )
249 for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
251 __p(c1, e*npoints+p) = __c->xReal(p)[c1];
254 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
256 __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr.evalq( c1, 0, p );
262 case GeomapStrategyType::GEOMAP_O1:
265 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
266 tensor_expr1.update( mapgmc1 );
268 for ( uint16_type p = 0; p < npoints; ++p )
270 for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
272 __p(c1, e*npoints+p) = __c1->xReal(p)[c1];
274 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
276 __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr1.evalq( c1, 0, p );
282 case GeomapStrategyType::GEOMAP_OPT:
284 if ( it->isOnBoundary() )
288 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
289 tensor_expr.update( mapgmc );
291 for ( uint16_type p = 0; p < npoints; ++p )
293 for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
295 __p(c1, e*npoints+p) = __c->xReal(p)[c1];
297 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
299 __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr.evalq( c1, 0, p );
307 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
308 tensor_expr1.update( mapgmc1 );
311 for ( uint16_type p = 0; p < npoints; ++p )
313 for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
315 __p(c1, e*npoints+p) = __c1->xReal(p)[c1];
317 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
319 __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr1.evalq( c1, 0, p );
327 return boost::make_tuple( __v, __p );
330 template<EvaluatorType iDim,
typename Iterator,
typename Pset,
typename ExprT>
331 typename Evaluator<iDim, Iterator, Pset, ExprT>::eval_element_type
332 Evaluator<iDim, Iterator, Pset, ExprT>::operator()( mpl::size_t<MESH_FACES> )
const
335 boost::timer __timer;
337 element_type __v( M_functionspace );
340 DVLOG(2) <<
"call project(MESH_FACES) " <<
"\n";
346 typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
347 typedef typename geoelement_type::face_type face_type;
350 typedef typename geoelement_type::gm_type gm_type;
351 typedef boost::shared_ptr<gm_type> gm_ptrtype;
352 typedef typename geoelement_type::gm1_type gm1_type;
353 typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
355 typedef typename gm_type::template Context<context, geoelement_type> gmc_type;
356 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
357 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc_ptrtype> > map_gmc_type;
358 typedef typename gm1_type::template Context<context, geoelement_type> gmc1_type;
359 typedef boost::shared_ptr<gmc1_type> gmc1_ptrtype;
360 typedef fusion::map<fusion::pair<vf::detail::gmc<0>, gmc1_ptrtype> > map_gmc1_type;
364 typedef typename element_type::functionspace_type::dof_type dof_type;
367 typedef typename element_type::functionspace_type::fe_type fe_type;
368 typedef typename fe_type::template Context< context, fe_type, gm_type, geoelement_type> fecontext_type;
369 typedef boost::shared_ptr<fecontext_type> fecontext_ptrtype;
370 typedef typename fe_type::template Context< context, fe_type, gm1_type, geoelement_type> fecontext1_type;
371 typedef boost::shared_ptr<fecontext1_type> fecontext1_ptrtype;
377 typedef expression_type the_expression_type;
378 typedef typename boost::remove_reference<typename boost::remove_const<the_expression_type>::type >::type iso_expression_type;
379 typedef typename iso_expression_type::template tensor<map_gmc_type> t_expr_type;
380 typedef typename iso_expression_type::template tensor<map_gmc1_type> t_expr1_type;
381 typedef typename t_expr_type::shape shape;
386 DVLOG(2) <<
"assembling Dirichlet conditions\n";
388 dof_type
const* __dof = __v.functionSpace()->dof().get();
390 fe_type
const* __fe = __v.functionSpace()->fe().get();
392 iterator_type __face_it, __face_en;
393 boost::tie( boost::tuples::ignore, __face_it, __face_en ) = M_range;
395 if ( __face_it == __face_en )
398 gm_ptrtype __gm(
new gm_type );
399 gm1_ptrtype __gm1(
new gm1_type );
407 typedef typename geoelement_type::permutation_type permutation_type;
408 typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
409 typedef typename gm_type::precompute_type geopc_type;
410 typedef typename gm1_type::precompute_ptrtype geopc1_ptrtype;
411 typedef typename gm1_type::precompute_type geopc1_type;
413 DVLOG(2) <<
"[integratoron] numTopologicalFaces = " << geoelement_type::numTopologicalFaces <<
"\n";
414 std::vector<std::map<permutation_type, geopc_ptrtype> > __geopc( geoelement_type::numTopologicalFaces );
415 std::vector<std::map<permutation_type, geopc1_ptrtype> > __geopc1( geoelement_type::numTopologicalFaces );
417 for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
419 for ( permutation_type __p( permutation_type::IDENTITY );
420 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
422 __geopc[__f][__p] = geopc_ptrtype(
new geopc_type( __gm, __fe->points( __f ) ) );
423 __geopc1[__f][__p] = geopc1_ptrtype(
new geopc1_type( __gm1, __fe->points( __f ) ) );
425 FEELPP_ASSERT( __geopc[__f][__p]->nPoints() ).error(
"invalid number of points for geopc" );
426 FEELPP_ASSERT( __geopc1[__f][__p]->nPoints() ).error(
"invalid number of points for geopc1" );
430 uint16_type __face_id = __face_it->pos_first();
431 gmc_ptrtype __c(
new gmc_type( __gm, __face_it->element( 0 ), __geopc, __face_id ) );
432 gmc1_ptrtype __c1(
new gmc1_type( __gm1, __face_it->element( 0 ), __geopc1, __face_id ) );
434 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
435 t_expr_type expr( basis_type::isomorphism( M_expr ), mapgmc );
436 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
437 t_expr1_type expr1( basis_type::isomorphism( M_expr ), mapgmc1 );
444 if ( !fe_type::is_modal )
445 nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
446 face_type::numEdges * fe_type::nDofPerEdge +
447 face_type::numFaces * fe_type::nDofPerFace );
450 nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
452 DVLOG(2) <<
"[projector::operator(MESH_FACES)] nbFaceDof = " << nbFaceDof <<
"\n";
454 std::vector<int> dofs;
455 std::vector<value_type> values;
457 for ( ; __face_it != __face_en; ++__face_it )
459 FEELPP_ASSERT( __face_it->isOnBoundary() && !__face_it->isConnectedTo1() )
460 ( __face_it->marker() )
461 ( __face_it->isOnBoundary() )
462 ( __face_it->ad_first() )
463 ( __face_it->pos_first() )
464 ( __face_it->ad_second() )
465 ( __face_it->pos_second() )
466 ( __face_it->id() ).warn(
"inconsistent data face" );
467 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id()
468 <<
" element id= " << __face_it->ad_first()
469 <<
" pos in elt= " << __face_it->pos_first()
470 <<
" marker: " << __face_it->marker() <<
"\n";
471 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id() <<
" real pts=" << __face_it->G() <<
"\n";
473 uint16_type __face_id = __face_it->pos_first();
475 std::pair<size_type,size_type> range_dof( std::make_pair( __v.start(),
476 __v.functionSpace()->nDof() ) );
477 DVLOG(2) <<
"[projector] dof start = " << range_dof.first <<
"\n";
478 DVLOG(2) <<
"[projector] dof range = " << range_dof.second <<
"\n";
480 switch ( M_geomap_strategy )
483 case GeomapStrategyType::GEOMAP_OPT:
484 case GeomapStrategyType::GEOMAP_HO:
486 __c->update( __face_it->element( 0 ), __face_id );
487 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id() <<
" ref pts=" << __c->xRefs() <<
"\n";
488 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id() <<
" real pts=" << __c->xReal() <<
"\n";
490 map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
492 expr.update( mapgmc );
494 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
495 for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
497 for ( uint16_type l = 0; l < nbFaceDof; ++l )
499 typename expression_type::value_type __value = expr.evalq( c1, c2, l );
502 boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
504 __v( thedof ) = __value;
510 case GeomapStrategyType::GEOMAP_O1:
512 __c1->update( __face_it->element( 0 ), __face_id );
513 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id() <<
" ref pts=" << __c1->xRefs() <<
"\n";
514 DVLOG(2) <<
"[projector] FACE_ID = " << __face_it->id() <<
" real pts=" << __c1->xReal() <<
"\n";
516 map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
518 expr1.update( mapgmc1 );
521 for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
522 for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
524 for ( uint16_type l = 0; l < nbFaceDof; ++l )
526 typename expression_type::value_type __value = expr1.evalq( c1, c2, l );
529 boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
531 __v( thedof ) = __value;
544 return boost::make_tuple( __v, __p );
552 template<
typename Args>
555 typedef typename clean_type<Args,tag::expr>::type _expr_type;
556 typedef typename clean_type<Args,tag::pset>::type _pset_type;
557 typedef typename clean_type<Args,tag::range>::type _range_type;
558 typedef details::Evaluator<EVAL_NODAL, _range_type, _pset_type, Expr<_expr_type> > eval_t;
559 typedef typename eval_t::mesh_element_type mesh_element_type;
560 typedef typename eval_t::eval_element_type element_type;
561 static const uint16_type nDim = mesh_element_type::nDim;
567 template<
typename IteratorRange,
typename Pset,
typename ExprT>
568 typename details::Evaluator<EVAL_NODAL, IteratorRange, Pset, Expr<ExprT> >::eval_element_type
569 evaluate_impl( IteratorRange
const& range_it,
571 Expr<ExprT>
const& __expr,
572 GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
574 typedef details::Evaluator<EVAL_NODAL, IteratorRange, Pset, Expr<ExprT> > proj_t;
575 proj_t p( range_it, pset, __expr, geomap );
589 BOOST_PARAMETER_FUNCTION(
590 (
typename vf::detail::evaluate<Args>::element_type ),
602 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
606 LOG(INFO) <<
"evaluate expression..." << std::endl;
607 return evaluate_impl( range, pset, expr, geomap );
608 LOG(INFO) <<
"evaluate expression done." << std::endl;
620 BOOST_PARAMETER_FUNCTION(
621 ( boost::tuple<
double, Eigen::Matrix<
double, vf::detail::evaluate<Args>::nDim,1> > ),
633 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
638 int proc_number = Environment::worldComm().globalRank();
640 LOG(INFO) <<
"evaluate expression..." << std::endl;
642 auto e = evaluate_impl( range, pset, expr, geomap );
644 double maxe = e.template get<0>().array().abs().maxCoeff(&index);
646 Eigen::Matrix<double, vf::detail::evaluate<Args>::nDim,1> n = e.template get<1>().col(index);
647 LOG(INFO) <<
"proc "<<proc_number<<
" index at which function (size: " << e.template get<0>().array().size() <<
") is maximal: "<< index <<
" coord = \n"<<n<<
"\n";
649 int world_size = Environment::worldComm().size();
650 std::vector<double> maxe_world( world_size );
651 mpi::all_gather( Environment::worldComm().globalComm(),
655 std::vector< Eigen::Matrix<double, vf::detail::evaluate<Args>::nDim,1> > n_world( world_size );
656 mpi::all_gather( Environment::worldComm().globalComm(),
660 auto it_max = std::max_element( maxe_world.begin() , maxe_world.end() );
661 int position = it_max - maxe_world.begin();
662 LOG(INFO)<<
"proc "<<proc_number<<
" : global max = "<<*it_max<<
" at position "<<position<<
" with coord : \n "<<n<<
"\n";
668 for(
int i = 0; i < e.template get<0>().size(); ++i )
670 if ( math::abs(e.template get<0>()(i)) > maxe2 )
672 maxe2 = math::abs(e.template get<0>()(i));
677 LOG_ASSERT( index2 == index ) <<
" index2 = " << index2 <<
" and index = " << index <<
"\n";
678 LOG(INFO) <<
"evaluate expression done." << std::endl;
680 return boost::make_tuple( *it_max, n_world[position] );
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319