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
evaluator.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: 2012-05-27
7 
8  Copyright (C) 2012 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 __FEELPP_EVALUATORS_H
30 #define __FEELPP_EVALUATORS_H 1
31 
32 #include <boost/timer.hpp>
33 #include <feel/feelcore/parameter.hpp>
34 #include <feel/feeldiscr/functionspace.hpp>
35 
36 namespace Feel
37 {
38 namespace vf
39 {
41 enum EvaluatorType
42 {
43  EVAL_NODAL = 0
44 };
45 namespace details
46 {
53 template<EvaluatorType iDim, typename IteratorRange, typename Pset, typename ExprT>
54 class Evaluator
55 {
56 public:
57 
58 
62 
63  static const size_type context = ExprT::context|vm::POINT;
64 
65  typedef ExprT expression_type;
66  typedef typename expression_type::value_type value_type;
67 
68 
69  static const uint16_type imorder = 1;
70  static const bool imIsPoly = true;
71 
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;
84 
88 
89  Evaluator( IteratorRange const& r,
90  Pset const& pset,
91  expression_type const& __expr,
92  GeomapStrategyType geomap_strategy )
93  :
94  M_range( r ),
95  M_pset(),
96  M_expr( __expr ),
97  M_geomap_strategy( geomap_strategy )
98  {
99  DVLOG(2) << "Evaluator constructor from expression\n";
100  }
101 
102 
103  Evaluator( Evaluator const& __vfi )
104  :
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 )
109  {
110  DVLOG(2) << "Evaluator copy constructor\n";
111  }
112 
113  virtual ~Evaluator() {}
114 
116 
120 
121  eval_element_type operator()() const
122  {
123  return operator()( idim_type() );
124  }
125 
127 
131 
138  expression_type const& expression() const
139  {
140  return M_expr;
141  }
142 
144 
148 
150 
154 
156 
157 private:
158 
159  eval_element_type operator()( mpl::size_t<MESH_ELEMENTS> ) const;
160  eval_element_type operator()( mpl::size_t<MESH_FACES> ) const;
161 
162 private:
163 
164  range_iterator M_range;
165  pointset_type M_pset;
166  expression_type const& M_expr;
167  GeomapStrategyType M_geomap_strategy;
168 };
169 
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
173 {
174  boost::timer __timer;
175 
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;
180 
181 
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;
186  //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;
187  //typedef decltype( basis_type::isomorphism( M_expr ) ) the_expression_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;
193 
194  // we should manipulate the same type of functions on the left and
195  // on the right
196  //BOOST_STATIC_ASSERT(( boost::is_same<return_value_type, typename functionspace_type::return_value_type>::value ));
197 
198  typedef typename t_expr_type::shape shape;
199 
200 
201  iterator_type it, en;
202  boost::tie( boost::tuples::ignore, it, en ) = M_range;
203 
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 ) );
207  __v.setZero();
208  __p.setZero();
209 
210 
211 
212  // return if no elements
213  if ( it == en )
214  return boost::make_tuple( __v, __p );
215 
216  //
217  // Precompute some data in the reference element for
218  // geometric mapping and reference finite element
219  //
220  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( it->gm(),
221  M_pset.points() ) );
222  typename gm1_type::precompute_ptrtype __geopc1( new typename gm1_type::precompute_type( it->gm1(),
223  M_pset.points() ) );
224 
225 
226 
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 ) );
229 
230  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
231  t_expr_type tensor_expr( M_expr, mapgmc );
232 
233  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
234 
235  t_expr1_type tensor_expr1( M_expr, mapgmc1 );
236 
237  for ( int e = 0; it!=en ; ++it, ++e )
238  {
239  switch ( M_geomap_strategy )
240  {
241  case GeomapStrategyType::GEOMAP_HO:
242  {
243  __c->update( *it );
244  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
245  tensor_expr.update( mapgmc );
246 
247  for ( uint16_type p = 0; p < npoints; ++p )
248  {
249  for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
250  {
251  __p(c1, e*npoints+p) = __c->xReal(p)[c1];
252  }
253 
254  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
255  {
256  __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr.evalq( c1, 0, p );
257  }
258  }
259  }
260  break;
261 
262  case GeomapStrategyType::GEOMAP_O1:
263  {
264  __c1->update( *it );
265  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
266  tensor_expr1.update( mapgmc1 );
267 
268  for ( uint16_type p = 0; p < npoints; ++p )
269  {
270  for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
271  {
272  __p(c1, e*npoints+p) = __c1->xReal(p)[c1];
273  }
274  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
275  {
276  __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr1.evalq( c1, 0, p );
277  }
278  }
279  }
280  break;
281 
282  case GeomapStrategyType::GEOMAP_OPT:
283  {
284  if ( it->isOnBoundary() )
285  {
286  // HO if on boundary
287  __c->update( *it );
288  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
289  tensor_expr.update( mapgmc );
290 
291  for ( uint16_type p = 0; p < npoints; ++p )
292  {
293  for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
294  {
295  __p(c1, e*npoints+p) = __c->xReal(p)[c1];
296  }
297  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
298  {
299  __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr.evalq( c1, 0, p );
300  }
301  }
302  }
303 
304  else
305  {
306  __c1->update( *it );
307  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
308  tensor_expr1.update( mapgmc1 );
309 
310 
311  for ( uint16_type p = 0; p < npoints; ++p )
312  {
313  for ( uint16_type c1 = 0; c1 < mesh_element_type::nDim; ++c1 )
314  {
315  __p(c1, e*npoints+p) = __c1->xReal(p)[c1];
316  }
317  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
318  {
319  __v( e*npoints*shape::M+shape::M*p+c1) = tensor_expr1.evalq( c1, 0, p );
320  }
321  }
322  }
323  }
324  break;
325  }
326  }
327  return boost::make_tuple( __v, __p );
328 }
329 
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
333 {
334 #if 0
335  boost::timer __timer;
336 
337  element_type __v( M_functionspace );
338  __v.setZero();
339 
340  DVLOG(2) << "call project(MESH_FACES) " << "\n";
341  //
342  // a few typedefs
343  //
344 
345  // mesh element
346  typedef typename element_type::functionspace_type::mesh_type::element_type geoelement_type;
347  typedef typename geoelement_type::face_type face_type;
348 
349  // geometric mapping context
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;
354 
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;
361 
362 
363  // dof
364  typedef typename element_type::functionspace_type::dof_type dof_type;
365 
366  // basis
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;
372  //typedef fusion::map<fusion::pair<vf::detail::gmc<0>, fecontext_ptrtype> > map_gmc_type;
373 
374  // expression
375  //typedef typename expression_type::template tensor<map_gmc_type,fecontext_type> t_expr_type;
376  //typedef decltype( basis_type::isomorphism( M_expr ) ) the_expression_type;
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;
382 
383  //
384  // start
385  //
386  DVLOG(2) << "assembling Dirichlet conditions\n";
387 
388  dof_type const* __dof = __v.functionSpace()->dof().get();
389 
390  fe_type const* __fe = __v.functionSpace()->fe().get();
391 
392  iterator_type __face_it, __face_en;
393  boost::tie( boost::tuples::ignore, __face_it, __face_en ) = M_range;
394 
395  if ( __face_it == __face_en )
396  return __v;
397 
398  gm_ptrtype __gm( new gm_type );
399  gm1_ptrtype __gm1( new gm1_type );
400 
401 
402 
403  //
404  // Precompute some data in the reference element for
405  // geometric mapping and reference finite element
406  //
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;
412 
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 );
416 
417  for ( uint16_type __f = 0; __f < geoelement_type::numTopologicalFaces; ++__f )
418  {
419  for ( permutation_type __p( permutation_type::IDENTITY );
420  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
421  {
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 ) ) );
424  //DVLOG(2) << "[geopc] FACE_ID = " << __f << " ref pts=" << __fe->dual().points( __f ) << "\n";
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" );
427  }
428  }
429 
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 ) );
433 
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 );
438 
439 
440 
441 
443 
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 );
448 
449  else
450  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
451 
452  DVLOG(2) << "[projector::operator(MESH_FACES)] nbFaceDof = " << nbFaceDof << "\n";
453 
454  std::vector<int> dofs;
455  std::vector<value_type> values;
456 
457  for ( ; __face_it != __face_en; ++__face_it )
458  {
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";
472 
473  uint16_type __face_id = __face_it->pos_first();
474 
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";
479 
480  switch ( M_geomap_strategy )
481  {
482  default:
483  case GeomapStrategyType::GEOMAP_OPT:
484  case GeomapStrategyType::GEOMAP_HO:
485  {
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";
489 
490  map_gmc_type mapgmc( fusion::make_pair<vf::detail::gmc<0> >( __c ) );
491 
492  expr.update( mapgmc );
493 
494  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
495  for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
496  {
497  for ( uint16_type l = 0; l < nbFaceDof; ++l )
498  {
499  typename expression_type::value_type __value = expr.evalq( c1, c2, l );
500 
501  size_type thedof = __v.start() +
502  boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
503 
504  __v( thedof ) = __value;
505  }
506  }
507  }
508  break;
509 
510  case GeomapStrategyType::GEOMAP_O1:
511  {
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";
515 
516  map_gmc1_type mapgmc1( fusion::make_pair<vf::detail::gmc<0> >( __c1 ) );
517 
518  expr1.update( mapgmc1 );
519 
520 
521  for ( uint16_type c1 = 0; c1 < shape::M; ++c1 )
522  for ( uint16_type c2 = 0; c2 < shape::N; ++c2 )
523  {
524  for ( uint16_type l = 0; l < nbFaceDof; ++l )
525  {
526  typename expression_type::value_type __value = expr1.evalq( c1, c2, l );
527 
528  size_type thedof = __v.start() +
529  boost::get<0>( __dof->faceLocalToGlobal( __face_it->id(), l, c1 ) );
530 
531  __v( thedof ) = __value;
532  }
533  }
534  }
535  break;
536  }
537 
538  } // face_it
539 #else
540  element_type __v;
541  node_type __p;
542 #endif
543 
544  return boost::make_tuple( __v, __p );
545 }
546 }
548 
550 namespace detail
551 {
552 template<typename Args>
553 struct evaluate
554 {
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;
562 };
563 }
565 
566 
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,
570  Pset const& pset,
571  Expr<ExprT> const& __expr,
572  GeomapStrategyType geomap = GeomapStrategyType::GEOMAP_HO )
573 {
574  typedef details::Evaluator<EVAL_NODAL, IteratorRange, Pset, Expr<ExprT> > proj_t;
575  proj_t p( range_it, pset, __expr, geomap );
576  return p();
577 }
578 
579 
589 BOOST_PARAMETER_FUNCTION(
590  ( typename vf::detail::evaluate<Args>::element_type ), // return type
591  evaluate, // 2. function name
592 
593  tag, // 3. namespace of tag types
594 
595  ( required
596  ( range, * )
597  ( pset, * )
598  ( expr, * )
599  ) // 4. one required parameter, and
600 
601  ( optional
602  ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
603  )
604 )
605 {
606  LOG(INFO) << "evaluate expression..." << std::endl;
607  return evaluate_impl( range, pset, expr, geomap );
608  LOG(INFO) << "evaluate expression done." << std::endl;
609 }
610 
620 BOOST_PARAMETER_FUNCTION(
621  ( boost::tuple<double, Eigen::Matrix<double, vf::detail::evaluate<Args>::nDim,1> > ), // return type
622  normLinf, // 2. function name
623 
624  tag, // 3. namespace of tag types
625 
626  ( required
627  ( range, * )
628  ( pset, * )
629  ( expr, * )
630  ) // 4. one required parameter, and
631 
632  ( optional
633  ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
634  )
635 )
636 {
637 
638  int proc_number = Environment::worldComm().globalRank();
639 
640  LOG(INFO) << "evaluate expression..." << std::endl;
641 
642  auto e = evaluate_impl( range, pset, expr, geomap );
643  int index;
644  double maxe = e.template get<0>().array().abs().maxCoeff(&index);
645 
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";
648 
649  int world_size = Environment::worldComm().size();
650  std::vector<double> maxe_world( world_size );
651  mpi::all_gather( Environment::worldComm().globalComm(),
652  maxe,
653  maxe_world );
654 
655  std::vector< Eigen::Matrix<double, vf::detail::evaluate<Args>::nDim,1> > n_world( world_size );
656  mpi::all_gather( Environment::worldComm().globalComm(),
657  n,
658  n_world );
659 
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";
663  //for(int i=0;i<vf::detail::evaluate<Args>::nDim;i++) LOG(INFO) << n_world[i](i) <<" - ";
664  //LOG(INFO)<<"\n";
665 
666  int index2=0;
667  double maxe2 = 0;
668  for( int i = 0; i < e.template get<0>().size(); ++i )
669  {
670  if ( math::abs(e.template get<0>()(i)) > maxe2 )
671  {
672  maxe2 = math::abs(e.template get<0>()(i));
673  index2 = i;
674  }
675  }
676 
677  LOG_ASSERT( index2 == index ) << " index2 = " << index2 << " and index = " << index << "\n";
678  LOG(INFO) << "evaluate expression done." << std::endl;
679 
680  return boost::make_tuple( *it_max, n_world[position] );
681 }
682 
683 } // vf
684 } // feel
685 
686 
687 #endif /* __FEELPP_EVALUATORS_H */
688 
689 

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