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
functionspace.hpp
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: 2004-11-22
7 
8  Copyright (C) 2004 EPFL
9  Copyright (C) 2006-2012 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 __FunctionSpace_H
31 #define __FunctionSpace_H 1
32 
33 #include <boost/static_assert.hpp>
34 
35 #include <boost/mpl/at.hpp>
36 #include <boost/mpl/vector.hpp>
37 #include <boost/mpl/transform.hpp>
38 #include <boost/fusion/support/pair.hpp>
39 #include <boost/fusion/support/is_sequence.hpp>
40 #include <boost/fusion/sequence.hpp>
41 #include <boost/fusion/container/vector.hpp>
42 #include <boost/fusion/algorithm.hpp>
43 #include <boost/fusion/adapted/mpl.hpp>
44 
45 #include <boost/serialization/vector.hpp>
46 #include <boost/serialization/array.hpp>
47 #include <boost/serialization/version.hpp>
48 
49 #include <boost/shared_ptr.hpp>
50 #include <boost/function.hpp>
51 #include <boost/numeric/ublas/vector.hpp>
52 #include <boost/numeric/ublas/vector_proxy.hpp>
53 // #include <boost/numeric/ublas/vector_serialize.hpp>
54 #include <boost/numeric/ublas/io.hpp>
55 #include <boost/lambda/lambda.hpp>
56 #include <boost/lambda/bind.hpp>
57 #include <boost/optional.hpp>
58 #include <boost/preprocessor/control/if.hpp>
59 
60 #include <boost/smart_ptr/enable_shared_from_this.hpp>
61 
62 
63 //#include<boost/filesystem.hpp>
64 
65 
66 #include <stdexcept>
67 #include <sstream>
68 #include <limits>
69 
70 #include <feel/feelcore/parameter.hpp>
72 
74 #include <feel/feelalg/glas.hpp>
76 
77 #include <feel/feelmesh/regiontree.hpp>
78 #include <feel/feelpoly/geomap.hpp>
79 
80 
81 #include <feel/feeldiscr/mesh.hpp>
86 #include <feel/feeldiscr/parameter.hpp>
87 #include <feel/feeldiscr/bases.hpp>
90 
92 
93 namespace Feel
94 {
95 namespace fusion = boost::fusion;
96 namespace parameter = boost::parameter;
97 
98 namespace detail
99 {
100 
101 template<typename T>
102 struct vector_plus
103 {
104  std::vector<T> operator()( std::vector<T> const& v1, std::vector<T> const& v2 ) const
105  {
106  FEELPP_ASSERT( v1.size() == v2.size() )( v1.size() )( v2.size() ).error( "invalid vector size for vector_plus<>" );
107  std::vector<T> res( v1.size() );
108 
109  for ( size_type i = 0; i < v1.size(); ++i )
110  res[i]=v1[i]+v2[i];
111 
112  return res;
113  }
114 };
115 template<typename T,int M, int N>
116 struct ID
117 {
118  friend class boost::serialization::access;
119  typedef T value_type;
120  typedef Eigen::Matrix<value_type,M,N> m_type;
121  typedef boost::multi_array<m_type,1> array_type;
122  typedef typename array_type::index_range range;
123 
124  struct result
125  {
126 
127  typedef typename array_type::template array_view<1>::type type;
128  };
129 
130  ID()
131  :
132  M_id()
133  {}
134 
135  ID( ID const& id )
136  :
137  M_id( id.M_id )
138  {}
139  ID( array_type const& id )
140  :
141  M_id( id )
142  {}
143 
144 
145  template<typename Elem, typename ContextType>
146  ID( Elem const& elem, ContextType const & context )
147  :
148  M_id( elem.idExtents( context ) )
149  {
150  elem.id_( context, M_id );
151  }
152 
153  ID& operator=( ID const& id )
154  {
155  if ( this != &id )
156  {
157  //VLOG(1) "[ID] extent = " <<
158 
159  M_id = id.M_id;
160  }
161 
162  return *this;
163 
164  }
165  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q ) const
166  {
167  return M_id[q]( c1,c2 );
168  }
169 
170  template <typename ExtentList>
171  void resize( const ExtentList& extents )
172  {
173  M_id.resize( extents );
174  }
175 
176  array_type M_id;
177 
178  template<class Archive>
179  void save( Archive & ar, const unsigned int /*version*/ ) const
180  {
181  size_type e1 = M_id.shape()[0];
182  DVLOG(2) << "saving in archive e1= " << e1 << "\n";
183  ar & e1;
184  size_type e2 = M_id.shape()[1];
185  DVLOG(2) << "saving in archive e2= " << e2 << "\n";
186  ar & e2;
187  size_type e3 = M_id.shape()[2];
188  DVLOG(2) << "saving in archive e3= " << e3 << "\n";
189  ar & e3;
190  DVLOG(2) << "saving in archive array of size = " << M_id.num_elements() << "\n";
191  ar & boost::serialization::make_array( M_id.data(), M_id.num_elements() );
192  DVLOG(2) << "saving in archive done\n";
193  }
194  template<class Archive>
195  void load( Archive & ar, const unsigned int /*version*/ )
196  {
197  size_type e1, e2, e3;
198  ar & e1;
199  DVLOG(2) << "loading from archive e1= " << e1 << "\n";
200  ar & e2;
201  DVLOG(2) << "loading from archive e2= " << e2 << "\n";
202  ar & e3;
203  DVLOG(2) << "loading from archive e3= " << e3 << "\n";
204  M_id.resize( boost::extents[e1] );
205  DVLOG(2) << "loading from archive array of size = " << M_id.num_elements() << "\n";
206  ar & boost::serialization::make_array( M_id.data(), M_id.num_elements() );
207  DVLOG(2) << "loading from archive done\n";
208  DVLOG(2) << "creating view interpolation context done\n";
209  }
210  BOOST_SERIALIZATION_SPLIT_MEMBER()
211 };
212 template<typename T,int M,int N>
213 struct DD
214 {
215  typedef T value_type;
216  friend class boost::serialization::access;
217  typedef Eigen::Matrix<value_type,M,N> m_type;
218  typedef boost::multi_array<m_type,1> array_type;
219  typedef typename array_type::index_range range;
220  struct result
221  {
222 
223  typedef array_type type;
224  };
225 
226  DD()
227  :
228  M_grad()
229  {}
230 
231  template<typename Elem, typename ContextType>
232  DD( Elem const& elem, ContextType const & context )
233  :
234  M_grad( elem.gradExtents( context ) )
235  {
236  elem.grad_( context, M_grad );
237  }
238 
239  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q ) const
240  {
241  return M_grad[q]( c1,c2 );
242  }
243  array_type M_grad;
244 
245  template<class Archive>
246  void save( Archive & ar, const unsigned int /*version*/ ) const
247  {
248  size_type e1 = M_grad.shape()[0];
249  DVLOG(2) << "saving in archive e1= " << e1 << "\n";
250  ar & e1;
251  size_type e2 = M_grad.shape()[1];
252  DVLOG(2) << "saving in archive e2= " << e2 << "\n";
253  ar & e2;
254  size_type e3 = M_grad.shape()[2];
255  DVLOG(2) << "saving in archive e3= " << e3 << "\n";
256  ar & e3;
257  DVLOG(2) << "saving in archive array of size = " << M_grad.num_elements() << "\n";
258  ar & boost::serialization::make_array( M_grad.data(), M_grad.num_elements() );
259  DVLOG(2) << "saving in archive done\n";
260  }
261  template<class Archive>
262  void load( Archive & ar, const unsigned int /*version*/ )
263  {
264  size_type e1, e2, e3;
265  ar & e1;
266  DVLOG(2) << "loading from archive e1= " << e1 << "\n";
267  ar & e2;
268  DVLOG(2) << "loading from archive e2= " << e2 << "\n";
269  ar & e3;
270  DVLOG(2) << "loading from archive e3= " << e3 << "\n";
271  M_grad.resize( boost::extents[e1] );
272  DVLOG(2) << "loading from archive array of size = " << M_grad.num_elements() << "\n";
273  ar & boost::serialization::make_array( M_grad.data(), M_grad.num_elements() );
274  DVLOG(2) << "loading from archive done\n";
275  DVLOG(2) << "creating view interpolation context done\n";
276  }
277  BOOST_SERIALIZATION_SPLIT_MEMBER()
278 
279 };
280 template<typename T,int N, int M, int P>
281 struct D
282 {
283  typedef T value_type;
284  typedef Eigen::Matrix<value_type,M,P> m_type;
285  typedef boost::multi_array<m_type,1> array_type;
286  typedef typename array_type::index_range range;
287 
288  struct result
289  {
290 
291  typedef typename array_type::template array_view<1>::type type;
292  };
293 
294  D()
295  :
296  M_grad()
297  {}
298 
299  template<typename Elem, typename ContextType>
300  D( Elem const& elem, ContextType const & context )
301  :
302  M_grad( elem.dExtents( context ) )
303  {
304  elem.d_( N, context, M_grad );
305  }
306 
307  value_type operator()( uint16_type c1, uint16_type /*c2*/, uint16_type q ) const
308  {
309  return M_grad[q]( c1,0 );
310  }
311  array_type M_grad;
312 };
313 
314 template<typename T>
315 struct Div
316 {
317  typedef T value_type;
318  typedef Eigen::Matrix<value_type,1,1> m_type;
319  typedef boost::multi_array<m_type,1> array_type;
320  typedef typename array_type::index_range range;
321  struct result
322  {
323 
324  typedef typename array_type::template array_view<1>::type type;
325  };
326 
327  Div()
328  :
329  M_div()
330  {}
331 
332  template<typename Elem, typename ContextType>
333  Div( Elem const& elem, ContextType const & context )
334  :
335  M_div( elem.divExtents( context ) )
336  {
337  elem.div_( context, M_div );
338 #if 0
339  uint16_type nComponents1 = elem.nComponents1;
340  std::fill( M_div.data(), M_div.data()+M_div.num_elements(), value_type( 0 ) );
341 
342  M_grad( elem.div_( context, pc, M_div ) ),
343 
344 
345  const uint16_type nq = context.xRefs().size2();
346 
347  for ( int c1 = 0; c1 < nComponents1; ++c1 )
348  for ( uint16_type q = 0; q < nq ; ++q )
349  {
350  M_div[q]( 0,0 ) += M_grad[q]( c1,c1 );
351  }
352 
353 #endif
354  }
355 
356  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q ) const
357  {
358  return M_div[q]( c1,c2 );
359  }
360  array_type M_div;
361 };
362 template<typename T, int N, int D>
363 struct Curl
364 {
365  typedef T value_type;
366  typedef Eigen::Matrix<value_type,D,1> m_type;
367  typedef boost::multi_array<m_type,1> array_type;
368  typedef typename array_type::index_range range;
369  struct result
370  {
371 
372  typedef typename array_type::template array_view<1>::type type;
373  };
374 
375  Curl()
376  :
377  M_curl()
378  {}
379 
380  template<typename Elem, typename ContextType>
381  Curl( Elem const& elem, ContextType const & context )
382  :
383  M_curl( elem.curlExtents( context ) )
384  {
385  init( elem, context, boost::is_same<mpl::int_<N>, mpl::int_<-1> >() );
386  }
387  template<typename Elem, typename ContextType>
388  void
389  init( Elem const& elem, ContextType const& context, mpl::bool_<true> )
390  {
391  elem.curl_( context, M_curl );
392  }
393  template<typename Elem, typename ContextType>
394  void
395  init( Elem const& elem, ContextType const& context, mpl::bool_<false> )
396  {
397  elem.curl_( context, M_curl, N );
398  }
399 
400  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q ) const
401  {
402  return this->operator()( c1, c2, q, mpl::int_<N>() );
403  }
404  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<-1> ) const
405  {
406  return M_curl[q]( c1,c2 );
407  }
408  value_type operator()( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q, mpl::int_<0> ) const
409  {
410  return M_curl[q]( 0,0 );
411  }
412  value_type operator()( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q, mpl::int_<1> ) const
413  {
414  return M_curl[q]( 0,0 );
415  }
416  value_type operator()( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q, mpl::int_<2> ) const
417  {
418  return M_curl[q]( 0,0 );
419  }
420  array_type M_curl;
421 };
422 
423 template<typename T,int M, int N>
424 struct H
425 {
426  friend class boost::serialization::access;
427  typedef T value_type;
428  typedef Eigen::Matrix<value_type,M,N> m_type;
429  typedef boost::multi_array<m_type,1> array_type;
430  typedef typename array_type::index_range range;
431 
432  struct result
433  {
434 
435  typedef typename array_type::template array_view<1>::type type;
436  };
437 
438  H()
439  :
440  M_hess()
441  {}
442 
443  template<typename Elem, typename ContextType>
444  H( Elem const& elem, ContextType const & context )
445  :
446  M_hess( elem.hessExtents( context ) )
447  {
448  elem.hess_( context, M_hess );
449  }
450 
451  value_type operator()( uint16_type c1, uint16_type c2, uint16_type q ) const
452  {
453  return M_hess[q]( c1,c2 );
454  }
455  array_type M_hess;
456 };
457 
458 template<typename MeshPtrType, typename PeriodicityType = Periodicity<NoPeriodicity> >
459 struct InitializeSpace
460 {
461  InitializeSpace( MeshPtrType const& mesh,
462  PeriodicityType const& periodicity,
463  std::vector<Dof> const& dofindices,
464  std::vector<WorldComm> const & worldsComm )
465  :
466  M_cursor( 0 ),
467  M_worldsComm( worldsComm ),
468  M_mesh( mesh ),
469  M_dofindices( dofindices ),
470  M_periodicity( periodicity )
471 
472  {}
473  template <typename T>
474  void operator()( boost::shared_ptr<T> & x ) const
475  {
476  operator()( x, is_shared_ptr<MeshPtrType>() );
477  }
478  template <typename T>
479  void operator()( boost::shared_ptr<T> & x, mpl::bool_<true> ) const
480  {
481  auto p = *fusion::find<typename T::periodicity_0_type>(M_periodicity);
482  x = boost::shared_ptr<T>( new T( M_mesh, M_dofindices, p,
483  std::vector<WorldComm>( 1,M_worldsComm[M_cursor] ) ) );
484  FEELPP_ASSERT( x ).error( "invalid function space" );
485 
486  ++M_cursor;// warning M_cursor < nb color
487  }
488  template <typename T>
489  void operator()( boost::shared_ptr<T> & x, mpl::bool_<false> ) const
490  {
491  auto p = *fusion::find<typename T::periodicity_0_type>(M_periodicity);
492 
493  // look for T::mesh_ptrtype in MeshPtrType
494  auto m = *fusion::find<typename T::mesh_ptrtype>(M_mesh);
495  x = boost::shared_ptr<T>( new T( m, M_dofindices, p,
496  std::vector<WorldComm>( 1,M_worldsComm[M_cursor] ) ) );
497  FEELPP_ASSERT( x ).error( "invalid function space" );
498 
499  ++M_cursor;// warning M_cursor < nb color
500  }
501  mutable uint16_type M_cursor;
502  std::vector<WorldComm> M_worldsComm;
503  MeshPtrType M_mesh;
504  std::vector<Dof> const& M_dofindices;
505  PeriodicityType M_periodicity;
506 };
507 template<typename DofType>
508 struct updateDataMapProcess
509 {
510  updateDataMapProcess( std::vector<WorldComm> const & worldsComm,
511  WorldComm const& worldCommFusion,
512  uint16_type lastCursor )
513  :
514  M_cursor( 0 ),
515  M_start_index( 0 ),
516  M_lastCursor( lastCursor ),
517  M_worldsComm( worldsComm ),
518  M_dm( new DofType( worldCommFusion ) ),
519  M_dmOnOff( new DofType( worldCommFusion ) )
520  {}
521 
522  template <typename T>
523  void operator()( boost::shared_ptr<T> & x ) const
524  {
525 
526  if ( M_worldsComm[M_cursor].isActive() )
527  {
528  size_type nLocWithGhost=x->nLocalDofWithGhost();
529  size_type nLocWithoutGhost=x->nLocalDofWithoutGhost();
530  M_dm->setFirstDof( M_dm->worldComm().globalRank(), x->dof()->firstDof() );
531  M_dm->setLastDof( M_dm->worldComm().globalRank(), x->dof()->lastDof() );
532  M_dm->setFirstDofGlobalCluster( M_dm->worldComm().globalRank(), M_start_index + x->dof()->firstDofGlobalCluster() );
533  M_dm->setLastDofGlobalCluster( M_dm->worldComm().globalRank(), M_start_index + x->dof()->lastDofGlobalCluster() );
534  M_dm->setNLocalDofWithoutGhost( M_dm->worldComm().globalRank(), x->dof()->nLocalDofWithoutGhost() );
535  M_dm->setNLocalDofWithGhost( M_dm->worldComm().globalRank(), x->dof()->nLocalDofWithGhost() );
536 
537  M_dm->resizeMapGlobalProcessToGlobalCluster( nLocWithGhost );
538  M_dm->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhost );
539 
540  for ( size_type i=0; i<nLocWithGhost; ++i )
541  M_dm->setMapGlobalProcessToGlobalCluster( i, M_start_index + x->dof()->mapGlobalProcessToGlobalCluster( i ) );
542 
543  for ( size_type i=0; i<nLocWithoutGhost; ++i )
544  M_dm->setMapGlobalClusterToGlobalProcess( i, x->dof()->mapGlobalClusterToGlobalProcess( i ) );
545  }
546 
547 
548  if ( M_cursor==0 )
549  {
550  M_dmOnOff->setFirstDof( M_dmOnOff->worldComm().globalRank(), x->dof()->firstDof() );
551  M_dmOnOff->setFirstDofGlobalCluster( M_dmOnOff->worldComm().globalRank(),
552  M_start_index + x->dof()->firstDofGlobalCluster() );
553  M_dmOnOff->setNLocalDofWithoutGhost( M_dmOnOff->worldComm().globalRank(),
554  0 );
555  M_dmOnOff->setNLocalDofWithGhost( M_dmOnOff->worldComm().globalRank(),
556  0 );
557  }
558 
559  if ( M_cursor==M_lastCursor )
560  {
561  M_dmOnOff->setLastDof( M_dmOnOff->worldComm().globalRank(),
562  M_start_index + x->dof()->lastDof() );
563  M_dmOnOff->setLastDofGlobalCluster( M_dmOnOff->worldComm().globalRank(),
564  M_start_index + x->dof()->lastDofGlobalCluster() );
565  }
566 
567  // update nLoc
568  size_type nLocWithoutGhostOnOff= M_dmOnOff->nLocalDofWithoutGhost() + x->dof()->nLocalDofWithoutGhost();
569  size_type nLocWithGhostOnOff= M_dmOnOff->nLocalDofWithGhost() + x->dof()->nLocalDofWithGhost();
570 
571  M_dmOnOff->setNLocalDofWithoutGhost( M_dmOnOff->worldComm().globalRank(),
572  nLocWithoutGhostOnOff );
573  M_dmOnOff->setNLocalDofWithGhost( M_dmOnOff->worldComm().globalRank(),
574  nLocWithGhostOnOff );
575 
576  // update map
577  M_dmOnOff->resizeMapGlobalProcessToGlobalCluster( nLocWithGhostOnOff );
578  M_dmOnOff->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhostOnOff );
579 
580  size_type startGlobClusterDof = M_dmOnOff->nLocalDofWithoutGhost() - x->dof()->nLocalDofWithoutGhost();
581  size_type startGlobProcessDof = M_dmOnOff->nLocalDofWithGhost() - x->dof()->nLocalDofWithGhost();
582 
583  for ( size_type i=0; i<x->dof()->nLocalDofWithGhost(); ++i )
584  {
585  M_dmOnOff->setMapGlobalProcessToGlobalCluster( startGlobProcessDof + i, M_start_index + x->dof()->mapGlobalProcessToGlobalCluster( i ) );
586  }
587 
588  for ( size_type i=0; i<x->dof()->nLocalDofWithoutGhost(); ++i )
589  {
590  M_dmOnOff->setMapGlobalClusterToGlobalProcess( startGlobClusterDof + i, x->dof()->mapGlobalClusterToGlobalProcess( i ) );
591  }
592 
593 
594  M_start_index+=x->nDof();
595 
596  ++M_cursor;// warning M_cursor < nb color
597  }
598 
599  boost::shared_ptr<DofType> dataMap() const
600  {
601  return M_dm;
602  }
603  boost::shared_ptr<DofType> dataMapOnOff() const
604  {
605  return M_dmOnOff;
606  }
607 
608  mutable uint16_type M_cursor;
609  mutable size_type M_start_index;
610  uint16_type M_lastCursor;
611  std::vector<WorldComm> M_worldsComm;
612  mutable boost::shared_ptr<DofType> M_dm;
613  mutable boost::shared_ptr<DofType> M_dmOnOff;
614 }; // updateDataMapProcess
615 
616 
617 
618 template<typename DofType>
619 struct updateDataMapProcessStandard
620 {
621  updateDataMapProcessStandard( std::vector<WorldComm> const & worldsComm,
622  WorldComm const& worldCommFusion,
623  uint16_type lastCursor,
624  std::vector<size_type> const& startDofGlobalCluster,
625  std::vector<size_type> nLocWithoutGhost, std::vector<size_type> nLocWithGhost)
626  :
627  M_cursor( 0 ),
628  M_start_index( worldCommFusion.globalSize(), 0 ),// M_startDofGlobalCluster[worldCommFusion.globalRank()] ),
629  M_startIndexWithGhost( 0 ),
630  M_lastCursor( lastCursor ),
631  M_worldsComm( worldsComm ),
632  M_dm( new DofType( worldCommFusion ) ),
633  M_startDofGlobalCluster(startDofGlobalCluster),
634  M_nLocWithoutGhost(nLocWithoutGhost),
635  M_nLocWithGhost(nLocWithGhost)
636  {
637 
638  const int myrank = M_dm->worldComm().globalRank();
639  const int worldsize = M_dm->worldComm().globalSize();
640  for (int proc = 0 ; proc < worldsize ; ++proc)
641  {
642  M_dm->setNLocalDofWithoutGhost( proc, nLocWithoutGhost[proc] );
643  M_dm->setNLocalDofWithGhost( proc, nLocWithGhost[proc] );
644  M_start_index[proc] = M_startDofGlobalCluster[proc];// nLocWithoutGhost;
645  }
646 
647  M_dm->resizeMapGlobalProcessToGlobalCluster( nLocWithGhost[myrank] );
648  M_dm->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhost[myrank] );
649  }
650 
651  template <typename T>
652  void operator()( boost::shared_ptr<T> & x ) const
653  {
654  const int myrank = M_dm->worldComm().globalRank();
655  const int worldsize = M_dm->worldComm().globalSize();
656  const size_type nLocWithGhost=x->dof()->nLocalDofWithGhost(myrank);
657  const size_type nLocWithoutGhost=x->dof()->nLocalDofWithoutGhost(myrank);
658 
659  for (int proc = 0 ; proc < worldsize ; ++proc)
660  {
661  if ( M_cursor==0 )
662  {
663  M_dm->setFirstDof( proc, x->dof()->firstDof(proc) );
664  M_dm->setFirstDofGlobalCluster( proc, M_startDofGlobalCluster[proc] );
665  M_dm->setLastDof( proc, x->dof()->lastDof(proc) );
666  if (x->dof()->nLocalDofWithoutGhost(proc) == 0)
667  M_dm->setLastDofGlobalCluster( proc, M_startDofGlobalCluster[proc] );
668  else
669  M_dm->setLastDofGlobalCluster( proc, M_startDofGlobalCluster[proc] + x->dof()->nLocalDofWithoutGhost(proc) - 1 );
670  }
671  else
672  {
673  M_dm->setLastDof( proc, M_dm->lastDof(proc) + x->dof()->nLocalDofWithGhost(proc) );
674  M_dm->setLastDofGlobalCluster( proc, M_dm->lastDofGlobalCluster(proc) + x->dof()->nLocalDofWithoutGhost(proc) );
675  }
676  }
677 
678  const size_type firstDofGC = M_dm->firstDofGlobalCluster();
679  const size_type firstDofGCSubSpace = x->dof()->firstDofGlobalCluster();
680 
681  for (size_type gdof = x->dof()->firstDof(myrank) ; gdof < x->dof()->nLocalDofWithGhost(myrank) ; ++gdof )
682  {
683  const size_type localDof = M_startIndexWithGhost+gdof;//nLocalDofStart+gdof;
684  const size_type gdofGC = x->dof()->mapGlobalProcessToGlobalCluster(gdof);
685 
686  if ( x->dof()->dofGlobalClusterIsOnProc( gdofGC ) )
687  {
688  const size_type globalDof = M_start_index[myrank] + ( gdofGC-firstDofGCSubSpace );
689  M_dm->setMapGlobalProcessToGlobalCluster( localDof, globalDof );
690  M_dm->setMapGlobalClusterToGlobalProcess( globalDof-firstDofGC ,localDof );
691  }
692  else
693  {
694  const int realproc = x->dof()->procOnGlobalCluster(gdofGC);
695  const size_type globalDof = M_start_index[realproc] + (gdofGC- x->dof()->firstDofGlobalCluster(realproc));
696  M_dm->setMapGlobalProcessToGlobalCluster( localDof, globalDof );
697  }
698  }
699 
700  for (int proc = 0 ; proc < worldsize ; ++proc)
701  M_start_index[proc] += x->dof()->nLocalDofWithoutGhost(proc);
702  M_startIndexWithGhost+=nLocWithGhost;
703 
704  ++M_cursor;
705  }
706 
707  boost::shared_ptr<DofType> dataMap() const
708  {
709  return M_dm;
710  }
711  boost::shared_ptr<DofType> dataMapOnOff() const
712  {
713  return M_dm;
714  }
715 
716  mutable uint16_type M_cursor;
717  mutable std::vector<size_type> M_start_index;
718  mutable size_type M_startIndexWithGhost;
719  uint16_type M_lastCursor;
720  std::vector<WorldComm> M_worldsComm;
721  mutable boost::shared_ptr<DofType> M_dm;
722  std::vector<size_type> M_startDofGlobalCluster;
723  std::vector<size_type> M_nLocWithoutGhost, M_nLocWithGhost;
724 }; // updateDataMapProcessStandard
725 
726 
727 
728 
729 
730 struct NbDof
731 {
732  typedef size_type result_type;
733  NbDof( size_type start = 0, size_type size = invalid_size_type_value )
734  :
735  M_cursor( start ),
736  M_finish( size )
737  {}
738  template<typename Sig>
739  struct result;
740 
741  template<typename T, typename S>
742 #if BOOST_VERSION < 104200
743  struct result<NbDof( T,S )>
744 #else
745  struct result<NbDof( S,T )>
746 #endif
747 :
748  boost::remove_reference<S>
749  {};
750  template <typename T>
751  size_type
752  operator()( T const& x, size_type s ) const
753  {
754  size_type ret = s;
755 
756  if ( !x )
757  return ret;
758 
759  if ( M_cursor < M_finish )
760  ret += x->nDof();
761 
762  ++M_cursor;
763  return ret;
764  }
765 
766  template <typename T>
767  size_type
768  operator()( size_type s, T const& x ) const
769  {
770  return this->operator()( x, s );
771  }
772 private:
773  mutable size_type M_cursor;
774  size_type M_finish;
775 };
776 
777 #if 0
778 struct NLocalDof
779 {
780  NLocalDof( size_type start = 0, size_type size = invalid_size_type_value )
781  :
782  M_cursor( start ),
783  M_finish( size )
784  {}
785  template<typename Sig>
786  struct result;
787 
788  template<typename T, typename S>
789 #if BOOST_VERSION < 104200
790  struct result<NLocalDof( T,S )>
791 #else
792  struct result<NLocalDof( S,T )>
793 #endif
794 :
795  boost::remove_reference<S>
796  {};
797  template <typename T>
798  size_type
799  operator()( T const& x, size_type s ) const
800  {
801  size_type ret = s;
802 
803  if ( M_cursor < M_finish )
804  ret += x->nLocalDof();
805 
806  ++M_cursor;
807  return ret;
808  }
809  template <typename T>
810  size_type
811  operator()( size_type s, T const& x ) const
812  {
813  return this->operator()( x, s );
814  }
815 private:
816  mutable size_type M_cursor;
817  size_type M_finish;
818 };
819 #else // MPI
820 template< typename IsWithGhostType>
821 struct NLocalDof
822 {
823 
824  NLocalDof( std::vector<WorldComm> const & worldsComm = std::vector<WorldComm>( 1,Environment::worldComm() ),
825  bool useOffSubSpace = false,
826  size_type start = 0, size_type size = invalid_size_type_value )
827  :
828  M_cursor( start ),
829  M_finish( size ),
830  M_worldsComm( worldsComm ),
831  M_useOffSubSpace( useOffSubSpace )
832  {}
833  template<typename Sig>
834  struct result;
835 
836  template<typename T, typename S>
837 #if BOOST_VERSION < 104200
838  struct result<NLocalDof( T,S )>
839 #else
840  struct result<NLocalDof( S,T )>
841 #endif
842 :
843  boost::remove_reference<S>
844  {};
845 
846  template <typename T>
847  size_type
848  nLocalDof( T const& x, mpl::bool_<true> ) const
849  {
850  return x->nLocalDofWithGhost();
851  }
852 
853  template <typename T>
854  size_type
855  nLocalDof( T const& x, mpl::bool_<false> ) const
856  {
857  return x->nLocalDofWithoutGhost();
858  }
859 
860  template <typename T>
861  size_type
862  operator()( T const& x, size_type s ) const
863  {
864  size_type ret = s;
865 
866  if ( M_cursor < M_finish )
867  {
868  if ( M_useOffSubSpace )
869  {
870  ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
871  }
872 
873  else
874  {
875  if ( M_worldsComm[M_cursor].isActive() )
876  ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
877  }
878  }
879 
880  ++M_cursor;
881  return ret;
882  }
883 
884  template <typename T>
885  size_type
886  operator()( size_type s, T const& x ) const
887  {
888  return this->operator()( x, s );
889  }
890 private:
891  mutable size_type M_cursor;
892  size_type M_finish;
893  std::vector<WorldComm> const& M_worldsComm;
894  bool M_useOffSubSpace;
895 };
896 #endif // end MPI
897 
898 
899 template< typename IsWithGhostType>
900 struct NLocalDofOnProc
901 {
902 
903  NLocalDofOnProc( const int proc,
904  std::vector<WorldComm> const & worldsComm = std::vector<WorldComm>( 1,Environment::worldComm() ),
905  bool useOffSubSpace = false,
906  size_type start = 0, size_type size = invalid_size_type_value )
907  :
908  M_proc(proc),
909  M_cursor( start ),
910  M_finish( size ),
911  M_worldsComm( worldsComm ),
912  M_useOffSubSpace( useOffSubSpace )
913  {}
914 
915  template<typename Sig>
916  struct result;
917 
918  template<typename T, typename S>
919 #if BOOST_VERSION < 104200
920  struct result<NLocalDofOnProc( T,S )>
921 #else
922  struct result<NLocalDofOnProc( S,T )>
923 #endif
924 :
925  boost::remove_reference<S>
926  {};
927 
928  template <typename T>
929  size_type
930  nLocalDof( T const& x, mpl::bool_<true> ) const
931  {
932  return x->nLocalDofWithGhostOnProc(M_proc);
933  }
934 
935  template <typename T>
936  size_type
937  nLocalDof( T const& x, mpl::bool_<false> ) const
938  {
939  return x->nLocalDofWithoutGhostOnProc(M_proc);
940  }
941 
942  template <typename T>
943  size_type
944  operator()( T const& x, size_type s ) const
945  {
946  size_type ret = s;
947 
948  if ( M_cursor < M_finish )
949  {
950  if ( M_useOffSubSpace )
951  {
952  ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
953  }
954 
955  else
956  {
957  if ( M_worldsComm[M_cursor].isActive() )
958  ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
959  }
960  }
961 
962  ++M_cursor;
963  return ret;
964  }
965 
966  template <typename T>
967  size_type
968  operator()( size_type s, T const& x ) const
969  {
970  return this->operator()( x, s );
971  }
972 private:
973  int M_proc;
974  mutable size_type M_cursor;
975  size_type M_finish;
976  std::vector<WorldComm> const& M_worldsComm;
977  bool M_useOffSubSpace;
978 }; // NLocalDofOnProc
979 
980 
981 template<int i,typename SpaceCompositeType>
982 struct InitializeContainersOff
983 {
984  InitializeContainersOff( boost::shared_ptr<SpaceCompositeType> const& _space )
985  :
986  M_cursor( 0 ),
987  M_space( _space )
988  {}
989  template <typename T>
990  void operator()( boost::shared_ptr<T> & x ) const
991  {
992  if ( M_cursor==i && !x )
993  x = boost::shared_ptr<T>( new T( M_space->template functionSpace<i>()->dof() ) );
994 
995  ++M_cursor;// warning M_cursor < nb color
996  }
997  mutable uint16_type M_cursor;
998  boost::shared_ptr<SpaceCompositeType> M_space;
999 };
1000 
1001 
1002 template<int i,typename SpaceCompositeType>
1003 struct SendContainersOn
1004 {
1005  SendContainersOn( boost::shared_ptr<SpaceCompositeType> const& _space,
1006  std::vector<double> const& _dataToSend )
1007  :
1008  M_cursor( 0 ),
1009  M_space( _space ),
1010  M_dataToSend( _dataToSend )
1011  {}
1012  template <typename T>
1013  void operator()( boost::shared_ptr<T> & x ) const
1014  {
1015  if ( M_cursor!=i )
1016  {
1017  int locRank=M_space->worldComm().localRank();
1018  int globRank=M_space->worldComm().localColorToGlobalRank( M_cursor,locRank );
1019  int tag = 0;
1020  //std::cout << "\n I am proc " << M_space->worldComm().globalRank()
1021  // << " I send to proc " << globRank << std::endl;
1022  M_space->worldComm().globalComm().send( globRank,tag,M_dataToSend );
1023  }
1024 
1025  ++M_cursor;// warning M_cursor < nb color
1026  }
1027  mutable uint16_type M_cursor;
1028  boost::shared_ptr<SpaceCompositeType> M_space;
1029  std::vector<double> M_dataToSend;
1030 };
1031 
1032 
1033 template<int i,typename SpaceCompositeType>
1034 struct RecvContainersOff
1035 {
1036  RecvContainersOff( boost::shared_ptr<SpaceCompositeType> const& _space )
1037  :
1038  M_cursor( 0 ),
1039  M_space( _space )
1040  {}
1041  template <typename T>
1042  void operator()( boost::shared_ptr<T> & x ) const
1043  {
1044  if ( M_cursor==i )
1045  {
1046  std::vector<double> dataToRecv( M_space->template functionSpace<i>()->nLocalDof() );
1047  int locRank=M_space->worldComm().localRank();
1048  int globRank=M_space->worldComm().localColorToGlobalRank( i,locRank );
1049  int tag = 0;//locRank;
1050  //std::cout << "\n I am proc " << M_space->worldComm().globalRank()
1051  // << " I recv to proc " << globRank << std::endl;
1052  M_space->worldComm().globalComm().recv( globRank,tag,dataToRecv );
1053  std::copy( dataToRecv.begin(), dataToRecv.end(), x->begin() );
1054  }
1055 
1056  ++M_cursor;// warning M_cursor < nb color
1057  }
1058  mutable uint16_type M_cursor;
1059  boost::shared_ptr<SpaceCompositeType> M_space;
1060 };
1061 
1062 
1063 
1064 template< typename map_type >
1065 struct searchIndicesBySpace
1066 {
1067  searchIndicesBySpace()
1068  {}
1069 
1070  searchIndicesBySpace( map_type& /*u*/ )
1071  {}
1072 
1073  template<typename T>
1074  searchIndicesBySpace( T const& fspace, map_type& u )
1075  {
1076  u = getIndicesFromSpace( fspace,u );
1077  }
1078 
1079  template<typename Sig>
1080  struct result;
1081 
1082  template<typename T, typename M>
1083 #if BOOST_VERSION < 104200
1084  struct result<searchIndicesBySpace( T,M )>
1085 #else
1086  struct result<searchIndicesBySpace( M,T )>
1087 #endif
1088 :
1089  boost::remove_reference<M>
1090  {};
1091 
1092  template < typename T >
1093  map_type getIndicesFromSpace( T const& fspace, map_type t ) const
1094  {
1095  if ( fspace->mesh()->numElements() == 0 )
1096  return t;
1097 
1098  size_type nProc = fspace->dof()->nProcessors();
1099 
1100  //search for the biggest index already in t; this will give the shift for the dofs
1101  std::vector< size_type > max_per_space;
1102 
1103  for ( size_type j=0; j<t.size(); j++ )
1104  {
1105  size_type _end = t[j].size();
1106 
1107  if ( _end )
1108  max_per_space.push_back( t[j][_end-1] );
1109  }
1110 
1111  //from all max indices found, determine the biggest
1112  size_type max_index = 0;
1113 
1114  if ( t.size() )
1115  max_index = *max_element( max_per_space.begin(), max_per_space.end() ) + 1;
1116 
1117  //std::cout << "maximum index " << max_index << "\n";
1118 
1119  //loop in all processors
1120  for ( size_type i=0; i<nProc; i++ )
1121  {
1122  /*
1123  std::cout << "Processor " << i << " has dofs"
1124  << " from " << fspace->dof()->firstDof(i)
1125  << " to " << fspace->dof()->lastDof(i) << "\n";
1126  */
1127 
1128  size_type _first = fspace->dof()->firstDof( i );
1129  size_type _last = fspace->dof()->lastDof( i );
1130 
1131  //the dofs numbering for the current space start at max_index+1
1132  for ( size_type j=_first; j<=_last; j++ )
1133  t[i].push_back( max_index + j );
1134  }
1135 
1136  return t;
1137  }
1138  template <typename T>
1139  map_type
1140 #if BOOST_VERSION < 104200
1141  operator()( T const& fspace, map_type t ) const
1142 #else
1143  operator()( map_type t, T const& fspace ) const
1144 #endif
1145  {
1146  return getIndicesFromSpace( fspace,t );
1147  }
1148 };
1149 
1150 // get start for each proc ->( proc0 : 0 ), (proc1 : sumdofproc0 ), (proc2 : sumdofproc0+sumdofproc1 ) ....
1151 struct computeStartOfFieldSplit
1152 {
1153  typedef boost::tuple< uint16_type , size_type > result_type;
1154 
1155  template<typename T>
1156  result_type operator()( result_type const & previousRes, T const& t )
1157  {
1158  auto cptSpaces = previousRes.get<0>();
1159  auto start = previousRes.get<1>();
1160 
1161  for (int proc=0;proc<t->dof()->worldComm().globalSize();++proc)
1162  {
1163  if (proc < t->dof()->worldComm().globalRank())
1164  start+=t->dof()->nLocalDofWithoutGhost(proc);
1165  }
1166  return boost::make_tuple( ++cptSpaces, start );
1167  }
1168 
1169 };
1170 
1171 // compute split
1172 struct computeNDofForEachSpace
1173 {
1174 
1175  computeNDofForEachSpace(size_type startSplit)
1176  :
1177  M_startSplit(startSplit)
1178  {}
1179 
1180  typedef boost::tuple< uint16_type, size_type, std::vector<std::vector<size_type> > > result_type;
1181 
1182  template<typename T>
1183  result_type operator()( result_type const & previousRes, T const& t )
1184  {
1185  const size_type nDofWithoutGhost = t->dof()->nLocalDofWithoutGhost();
1186  const size_type nDofWithGhost = t->dof()->nLocalDofWithGhost();
1187  const size_type firstDof = t->dof()->firstDofGlobalCluster();
1188  uint16_type cptSpaces = previousRes.get<0>();
1189  const size_type start = previousRes.get<1>();
1190  auto is = previousRes.get<2>();
1191 
1192  //std::cout << "compo " << cptSpaces << " start " << start <<" split nDofWithout " << nDof << " with ghost "<< nDofWithGhost<< " M_startSplit " << M_startSplit <<std::endl;
1193 
1194  //is.push_back( std::vector<size_type>( nDofWithoutGhost ) );
1195  is[cptSpaces].resize(nDofWithoutGhost);
1196 
1197  for ( size_type i=0; i<nDofWithGhost; ++i )
1198  {
1199  if ( t->dof()->dofGlobalProcessIsGhost(i) ) continue;
1200  const size_type globalDof = t->dof()->mapGlobalProcessToGlobalCluster(i);
1201  is[cptSpaces][globalDof - firstDof ] = M_startSplit + start + (globalDof - firstDof);
1202  }
1203 
1204  return boost::make_tuple( ++cptSpaces, ( start+nDofWithoutGhost ), is );
1205  }
1206 
1207  size_type M_startSplit;
1208 };
1209 
1210 struct rebuildDofPointsTool
1211 {
1212 
1213  template <typename T>
1214  void operator()( boost::shared_ptr<T> & x ) const
1215  {
1216  x->dof()->rebuildDofPoints( *x->mesh() );
1217  }
1218 };
1219 
1220 struct BasisName
1221 {
1222  typedef std::string result_type;
1223 
1224  template<typename T>
1225  result_type operator()( result_type const & previousRes, T const& t )
1226  {
1227  std::ostringstream os;
1228 
1229  if ( previousRes.size() )
1230  os << previousRes << "_" << t->basis()->familyName();
1231 
1232  else
1233  os << t->basis()->familyName();
1234 
1235  return os.str();
1236  }
1237 };
1238 
1239 struct BasisOrder
1240 {
1241  typedef std::vector<int> result_type;
1242 
1243  template<typename T>
1244  result_type operator()( result_type const & previousRes, T const& t )
1245  {
1246  std::vector<int> res( previousRes );
1247  res.push_back( t->nSubFunctionSpace() );
1248  return res;
1249  }
1250 };
1251 
1252 template<typename ElementType>
1253 struct CreateElementVector
1254 {
1255 public:
1256  template<typename Sig>
1257  struct result;
1258 
1259  template<typename Lhs, typename Rhs>
1260  struct result<CreateElementVector( Lhs,Rhs )>
1261  {
1262  typedef typename boost::remove_const<typename boost::remove_reference<Lhs>::type>::type lhs_noref_type;
1263  typedef typename boost::remove_const<typename boost::remove_reference<Rhs>::type>::type::element_type rhs_noref_type;
1264  typedef typename boost::remove_const<typename boost::remove_reference<ElementType>::type>::type ElementType_noref_type;
1265  typedef typename boost::fusion::result_of::size<lhs_noref_type>::type index;
1266  typedef typename ElementType_noref_type::template sub_element<index::value>::type elt_type;
1267  BOOST_MPL_ASSERT( ( boost::is_same<typename elt_type::functionspace_type,rhs_noref_type> ) );
1268  typedef typename boost::fusion::result_of::make_vector<elt_type>::type v_elt_type;
1269 
1270  typedef typename boost::fusion::result_of::push_back<lhs_noref_type, elt_type>::type ptype;
1271  typedef typename boost::fusion::result_of::as_vector<ptype>::type type;
1272  };
1273  CreateElementVector( ElementType const& e ) : M_e( e ), M_names() {}
1274  CreateElementVector( ElementType const& e, std::vector<std::string> const& names ) : M_e( e ), M_names( names ) {}
1275  ElementType const& M_e;
1276  std::vector<std::string> M_names;
1277 
1278  template<typename Lhs, typename Rhs>
1279  typename result<CreateElementVector( Lhs,Rhs )>::type
1280  operator()( Lhs const& lhs, Rhs const& rhs ) const
1281  {
1282  typedef typename boost::remove_const<typename boost::remove_reference<Lhs>::type>::type lhs_noref_type;
1283  typedef typename boost::remove_const<typename boost::remove_reference<Rhs>::type>::type rhs_noref_type;
1284  typedef typename boost::remove_const<typename boost::remove_reference<ElementType>::type>::type ElementType_noref_type;
1285  typedef typename boost::fusion::result_of::size<lhs_noref_type>::type index;
1286  typename ElementType_noref_type::template sub_element<index::value>::type elt = M_e.template element<index::value>();
1287  static const int s = mpl::size<typename ElementType::functionspace_type::bases_list>::type::value;
1288  BOOST_STATIC_ASSERT( (boost::is_same<decltype(elt), typename ElementType::template sub_element<index::value>::type>::value ) );
1289  if ( !M_names.empty() && M_names.size() > index::value )
1290  {
1291 
1292  FEELPP_ASSERT( M_names.size() == s )
1293  ( M_names.size() )( s ).error( "incompatible number of function names and functions");
1294  elt.setName( M_names[index::value] );
1295  }
1296  else if ( ( M_names.size() == 1 ) && s > 1 )
1297  {
1298  elt.setName( (boost::format( "%1%-%2%" ) % M_names[0] % index::value ).str() );
1299  }
1300  else
1301  {
1302  elt.setName( (boost::format( "%1%-%2%" ) % M_e.name() % index::value ).str() );
1303  }
1304  return boost::fusion::as_vector( boost::fusion::push_back( lhs, elt ) );
1305  }
1306 };
1307 
1308 } // detail
1309 
1310 enum ComponentType
1311 {
1312  NO_COMPONENT = -1,
1313  X = 0,
1314  Y,
1315  Z,
1316  NX,
1317  NY,
1318  NZ,
1319  TX,
1320  TY,
1321  TZ
1322 };
1323 enum FunctionSpaceType
1324 {
1325  SCALAR = 0,
1326  VECTORIAL = 1
1327 };
1328 
1329 template<uint16_type PN,
1330  uint16_type GN = 1>
1331 struct Order
1332 {
1333  static const uint16_type PolynomialOrder = PN;
1334  static const uint16_type GeometricOrder = GN;
1335 
1336  static const bool is_isoparametric = ( PN == GN );
1337  static const bool is_subparametric = ( PN > GN );
1338  static const bool is_surparametric = ( PN < GN );
1339 };
1340 
1341 typedef parameter::parameters<
1342 // parameter::required<tag::mesh_type, mpl::or_<boost::is_base_and_derived<MeshBase,_> >, mpl::or_<fusion::traits::is_sequence<_>, mpl::is_sequence<_> > >
1343 parameter::required<tag::mesh_type, boost::is_base_and_derived<MeshBase,_> >
1344 #if 0
1345 , parameter::optional<parameter::deduced<tag::bases_list>, mpl::or_<boost::is_base_and_derived<detail::bases_base,_>,
1346 mpl::or_<fusion::traits::is_sequence<_>,
1347 mpl::is_sequence<_> > > >
1348 #else
1349 , parameter::optional<parameter::deduced<tag::bases_list>, boost::is_base_and_derived<detail::bases_base,_> >
1350 //, parameter::optional<parameter::deduced<tag::bases_list>, fusion::traits::is_sequence<_> >
1351 #endif
1352 , parameter::optional<parameter::deduced<tag::value_type>, boost::is_floating_point<_> >
1353 , parameter::optional<parameter::deduced<tag::periodicity_type>, boost::is_base_and_derived<detail::periodicity_base,_> >
1354 > functionspace_signature;
1355 
1356 
1367 //template<typename MeshType, typename Basis_t, typename T_type = double, typename PeriodicityType = NoPeriodicity>
1368 template<
1369 typename A0,
1370  typename A1 = parameter::void_,
1371  typename A2 = parameter::void_,
1372  typename A3 = parameter::void_,
1373  typename A4 = parameter::void_>
1375  :
1376 public FunctionSpaceBase,
1377 public boost::enable_shared_from_this<FunctionSpace<A0,A1,A2,A3,A4> >
1378 {
1379 public:
1380  typedef typename functionspace_signature::bind<A0,A1,A2,A3,A4>::type args;
1381 
1382  typedef typename parameter::binding<args, tag::mesh_type>::type meshes_list;
1383  typedef typename parameter::binding<args, tag::value_type, double>::type value_type;
1384  typedef typename parameter::binding<args, tag::periodicity_type, Periodicity<NoPeriodicity> >::type periodicity_type;
1385  typedef typename parameter::binding<args, tag::bases_list, detail::bases<Lagrange<1,Scalar> > >::type bases_list;
1386 
1387  BOOST_MPL_ASSERT_NOT( ( boost::is_same<mpl::at<bases_list,mpl::int_<0> >, mpl::void_> ) );
1388 
1389 public:
1390 
1391  template<typename ThePeriodicityType, int pos>
1392  struct GetPeriodicity
1393  {
1394 #if 0
1395  typedef typename boost::remove_reference<periodicity_type>::type periodicity_list_noref;
1396  typedef typename fusion::result_of::at_c<periodicity_list_noref, pos>::type _type;
1397  typedef typename boost::remove_reference<_type>::type type;
1398 #else
1399  typedef typename mpl::if_<mpl::equal_to<fusion::result_of::size<ThePeriodicityType>,mpl::int_<1> >,
1400  mpl::identity<fusion::result_of::at_c<ThePeriodicityType,0> >,
1401  mpl::identity<fusion::result_of::at_c<ThePeriodicityType,pos> > >::type::type::type _type;
1402  typedef typename boost::remove_reference<_type>::type type;
1403 
1404 #endif
1405  };
1406  template<typename BasisType>
1407  struct ChangeMesh
1408  {
1409  typedef typename boost::remove_reference<meshes_list>::type meshes_list_noref;
1410  typedef typename boost::remove_reference<bases_list>::type bases_list_noref;
1411  typedef typename fusion::result_of::distance<typename fusion::result_of::begin<meshes_list_noref>::type,
1412  typename fusion::result_of::find<bases_list_noref,BasisType>::type>::type pos;
1413  typedef typename fusion::result_of::at_c<meshes_list_noref, pos::value >::type _mesh_type;
1414 
1416  detail::bases<BasisType>,value_type,
1417  Periodicity<typename GetPeriodicity<periodicity_type,pos::value>::type > > _type;
1418  typedef boost::shared_ptr<_type> type;
1419 
1420 
1421  };
1422  template<typename BasisType>
1423  struct ChangeBasis
1424  {
1425  //typedef typename mpl::if_<mpl::and_<boost::is_base_of<MeshBase, meshes_list >, boost::is_base_of<Feel::detail::periodic_base, periodicity_type > >,
1426  typedef typename boost::remove_reference<bases_list>::type bases_list_noref;
1427  typedef typename fusion::result_of::distance<typename fusion::result_of::begin<bases_list_noref>::type,
1428  typename fusion::result_of::find<bases_list_noref,BasisType>::type>::type pos;
1429 
1430  typedef typename mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1431  mpl::identity<mpl::identity<boost::shared_ptr<FunctionSpace<meshes_list,detail::bases<BasisType>,value_type, Periodicity<typename GetPeriodicity<periodicity_type,pos::value>::type> > > > >,
1432  mpl::identity<ChangeMesh<BasisType> > >::type::type::type type;
1433 
1434 //mpl::identity<typename mpl::transform<meshes_list, ChangeMesh<mpl::_1,BasisType>, mpl::back_inserter<fusion::vector<> > >::type > >::type::type type;
1435  };
1436  typedef typename mpl::transform<bases_list, ChangeBasis<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type functionspace_vector_type;
1437 
1438  template<typename BasisType>
1439  struct ChangeBasisToComponentBasis
1440  {
1441  typedef typename BasisType::component_basis_type component_basis_type;
1442  //typedef typename mpl::if_<mpl::and_<boost::is_base_of<MeshBase, meshes_list >, boost::is_base_of<Feel::detail::periodic_base, periodicity_type > >,
1443  typedef typename mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1444  mpl::identity<mpl::identity<boost::shared_ptr<FunctionSpace<meshes_list,detail::bases<component_basis_type>,value_type, periodicity_type> > > >,
1445  mpl::identity<ChangeMesh<component_basis_type> > >::type::type::type type;
1446  };
1447 
1448  typedef typename mpl::transform<bases_list,
1449  ChangeBasisToComponentBasis<mpl::_1>,
1450  mpl::back_inserter<fusion::vector<> > >::type component_functionspace_vector_type;
1451 
1452  template<typename BasisType>
1453  struct GetComponentBasis
1454  {
1455  typedef typename BasisType::component_basis_type type;
1456  };
1457 
1458  typedef detail::bases2<typename mpl::transform<bases_list,
1459  GetComponentBasis<mpl::_1>,
1460  mpl::back_inserter<fusion::vector<> > >::type> component_basis_vector_type;
1461 
1462 
1463 
1464 public:
1465 
1469  static const bool is_composite = ( mpl::size<bases_list>::type::value > 1 );
1470 
1471  template<typename MeshListType,int N>
1472  struct GetMesh
1473  {
1474  typedef typename mpl::if_<mpl::or_<boost::is_base_of<MeshBase, MeshListType >,
1475  is_shared_ptr<MeshListType> >,
1476  mpl::identity<mpl::identity<MeshListType> >,
1477  mpl::identity<mpl::at_c<MeshListType,N> > >::type::type::type type;
1478  };
1479  // mesh
1480  typedef meshes_list MeshesListType;
1481  typedef typename GetMesh<meshes_list,0>::type mesh_0_type;
1482  typedef typename mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1483  mpl::identity<meshes_list>,
1484  mpl::identity<mesh_0_type> >::type::type mesh_type;
1485 
1486  template<typename MeshType>
1487  struct ChangeToMeshPtr
1488  {
1489  typedef boost::shared_ptr<MeshType> type;
1490  };
1491 
1492  typedef typename mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1493  mpl::identity<boost::shared_ptr<mesh_type> >,
1494  mpl::identity<typename mpl::transform<meshes_list, ChangeToMeshPtr<mpl::_1>, mpl::back_inserter<meshes<> > >::type > >::type::type mesh_ptrtype;
1495  typedef typename mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1496  mpl::identity<typename mesh_type::element_type>,
1497  mpl::identity<typename mesh_0_type::element_type> >::type::type convex_type;
1498 
1499  template<typename BasisType>
1500  struct GetNComponents
1501  {
1502  //typedef mpl::int_<BasisType::template apply<mesh_type::nDim,value_type,typename mesh_type::element_type>::type::nComponents> type;
1503  typedef mpl::int_<BasisType::template apply<mesh_type::nDim,
1504  mesh_type::nRealDim,
1505  value_type,
1506  typename mesh_type::element_type>::type::nComponents> type;
1507  };
1508  struct nodim { static const int nDim = -1; static const int nRealDim = -1; };
1509  static const uint16_type nDim = mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1510  mpl::identity<meshes_list >,
1511  mpl::identity<nodim> >::type::type::nDim;
1512  static const uint16_type nRealDim = mpl::if_<boost::is_base_of<MeshBase, meshes_list >,
1513  mpl::identity<meshes_list>,
1514  mpl::identity<nodim> >::type::type::nRealDim;
1515 
1516 
1517 
1518  //typedef typename mpl::at_c<bases_list,0>::type::template apply<mesh_type::nDim,value_type,typename mesh_type::element_type>::type basis_0_type;
1519  typedef typename mpl::at_c<bases_list,0>::type::template apply<GetMesh<meshes_list,0>::type::nDim,
1520  GetMesh<meshes_list,0>::type::nRealDim,
1521  value_type,
1522  typename GetMesh<meshes_list,0>::type::element_type>::type basis_0_type;
1523 
1524  static const uint16_type rank = ( is_composite? invalid_uint16_type_value : basis_0_type::rank );
1525  static const bool is_scalar = ( is_composite? false : basis_0_type::is_scalar );
1526  static const bool is_vectorial = ( is_composite? false : basis_0_type::is_vectorial );
1527  static const bool is_tensor2 = ( is_composite? false : basis_0_type::is_tensor2 );
1528  static const bool is_continuous = ( is_composite? false : basis_0_type::isContinuous );
1529  static const bool is_modal = ( is_composite? false : basis_0_type::is_modal );
1530  static const uint16_type nComponents1 = ( is_composite? invalid_uint16_type_value : basis_0_type::nComponents1 );
1531  static const uint16_type nComponents2 = ( is_composite? invalid_uint16_type_value : basis_0_type::nComponents2 );
1532  static const bool is_product = ( is_composite? invalid_uint16_type_value : basis_0_type::is_product );
1533  typedef typename basis_0_type::continuity_type continuity_type;
1534 
1535  static const uint16_type nComponents = mpl::transform<bases_list,
1536  GetNComponents<mpl::_1>,
1537  mpl::inserter<mpl::int_<0>,mpl::plus<mpl::_,mpl::_> > >::type::value;
1538  static const uint16_type N_COMPONENTS = nComponents;
1539  static const uint16_type nSpaces = mpl::size<bases_list>::type::value;
1540 
1541  typedef typename GetPeriodicity<periodicity_type,0>::type periodicity_0_type;
1542  static const bool is_periodic = periodicity_0_type::is_periodic;
1543 
1545 
1549 
1550  typedef typename ublas::type_traits<value_type>::real_type real_type;
1551  typedef typename node<value_type>::type node_type;
1552 
1553 
1555  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
1556  typedef boost::shared_ptr<functionspace_type> pointer_type;
1557 
1559  typedef boost::shared_ptr<component_functionspace_type> component_functionspace_ptrtype;
1560 
1561 
1562  // basis
1563  typedef bases_list BasisType;
1564  template<int N>
1565  struct Basis
1566  {
1567  typedef typename mpl::at_c<bases_list,N>::type type;
1568  };
1569  typedef typename mpl::if_<mpl::bool_<is_composite>,
1570  mpl::identity<bases_list>,
1571  mpl::identity<basis_0_type> >::type::type basis_type;
1572 
1573 
1574  typedef boost::shared_ptr<basis_type> basis_ptrtype;
1575  typedef basis_type reference_element_type;
1576  typedef boost::shared_ptr<reference_element_type> reference_element_ptrtype;
1577  typedef reference_element_type fe_type;
1578  typedef reference_element_ptrtype fe_ptrtype;
1579  typedef typename mpl::if_<mpl::bool_<is_composite>,
1580  mpl::identity<boost::none_t>,
1581  mpl::identity<typename basis_0_type::PreCompute> >::type pc_type;
1582  typedef boost::shared_ptr<pc_type> pc_ptrtype;
1583 
1584  // component basis
1585 #if 0
1586  typedef typename mpl::if_<mpl::bool_<is_composite>,
1587  mpl::identity<component_basis_vector_type>,
1588  typename mpl::at_c<component_basis_vector_type,0>::type::template apply<nDim,value_type,convex_type> >::type::type component_basis_type;
1589 #else
1590  typedef typename mpl::if_<mpl::bool_<is_composite>,
1591  mpl::identity<component_basis_vector_type>,
1592  typename mpl::at_c<component_basis_vector_type,0>::type::template apply<nDim,
1593  nRealDim,
1594  value_type,
1595  convex_type> >::type::type component_basis_type;
1596 #endif
1597  typedef boost::shared_ptr<component_basis_type> component_basis_ptrtype;
1598 
1599  // trace space
1600  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1601  mpl::identity<typename mesh_type::trace_mesh_type>,
1602  mpl::identity<mpl::void_> >::type::type trace_mesh_type;
1603  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1604  mpl::identity<typename mesh_type::trace_mesh_ptrtype>,
1605  mpl::identity<mpl::void_> >::type::type trace_mesh_ptrtype;
1606  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1607  mpl::identity<FunctionSpace<trace_mesh_type, bases_list> >,
1608  mpl::identity<mpl::void_> >::type::type trace_functionspace_type;
1609  typedef typename boost::shared_ptr<trace_functionspace_type> trace_functionspace_ptrtype;
1610 
1611  // wirebasket
1612  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1613  mpl::identity<typename mesh_type::trace_trace_mesh_type>,
1614  mpl::identity<mpl::void_> >::type::type trace_trace_mesh_type;
1615  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1616  mpl::identity<typename mesh_type::trace_trace_mesh_ptrtype>,
1617  mpl::identity<mpl::void_> >::type::type trace_trace_mesh_ptrtype;
1618  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1619  mpl::identity<FunctionSpace<trace_trace_mesh_type, bases_list> >,
1620  mpl::identity<mpl::void_> >::type::type trace_trace_functionspace_type;
1621  typedef typename boost::shared_ptr<trace_trace_functionspace_type> trace_trace_functionspace_ptrtype;
1622 
1623 #if 0
1624  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<1> >,
1625  mpl::identity<typename trace_functionspace_type::element_type>,
1626  mpl::identity<mpl::void_> >::type::type trace_element_type;
1627 #endif
1628 
1629  // geomap
1630  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename mesh_type::gm_type>, mpl::identity<mpl::void_> >::type::type gm_type;
1631  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename mesh_type::gm1_type>, mpl::identity<mpl::void_> >::type::type gm1_type;
1632  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename mesh_type::element_type>, mpl::identity<mpl::void_> >::type::type geoelement_type;
1633  typedef boost::shared_ptr<gm_type> gm_ptrtype;
1634  typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
1635  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename gm_type::template Context<vm::POINT|vm::JACOBIAN|vm::HESSIAN|vm::KB, geoelement_type> >,
1636  mpl::identity<mpl::void_> >::type::type gmc_type;
1637  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1638  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename gm_type::precompute_ptrtype>, mpl::identity<mpl::void_> >::type::type geopc_ptrtype;
1639  typedef typename mpl::if_<mpl::greater<mpl::int_<nDim>, mpl::int_<0> >,mpl::identity<typename gm_type::precompute_type>, mpl::identity<mpl::void_> >::type::type geopc_type;
1640 
1641  // dof
1642  typedef typename mpl::if_<mpl::bool_<is_composite>,
1643  mpl::identity<DofComposite>,
1644  mpl::identity<DofTable<mesh_type, basis_type, periodicity_0_type> > >::type::type dof_type;
1645 
1646  typedef boost::shared_ptr<dof_type> dof_ptrtype;
1647 
1648  // return types
1649  //typedef typename bases_list::polyset_type return_value_type;
1650 
1651  typedef typename mpl::if_<mpl::equal_to<mpl::int_<N_COMPONENTS>, mpl::int_<1> >,
1652  mpl::identity<value_type>,
1653  mpl::identity<node_type> >::type::type return_type;
1654 
1655  typedef boost::function<return_type ( node_type const& )> function_type;
1657 
1658  template<int i>
1659  struct sub_functionspace
1660  {
1661  typedef typename mpl::at_c<functionspace_vector_type,i>::type type;
1662  };
1663 
1668 
1669 
1673  template<typename T = double, typename Cont = VectorUblas<T> >
1674  class Element
1675  :
1676  public Cont,boost::addable<Element<T,Cont> >, boost::subtractable<Element<T,Cont> >
1677  {
1678  public:
1679  typedef T value_type;
1680 
1681  template<typename BasisType>
1682  struct ChangeElement
1683  {
1684  typedef T value_type;
1685  BOOST_MPL_ASSERT_NOT( ( boost::is_same<BasisType,mpl::void_> ) );
1686  typedef typename ChangeBasis<BasisType>::type::element_type fs_type;
1687  typedef typename fs_type::template Element<value_type, typename VectorUblas<T>::range::type > element_type;
1688  typedef element_type type;
1689  };
1690  typedef typename mpl::transform<bases_list, ChangeElement<mpl::_1>, mpl::back_inserter<mpl::vector<> > >::type element_vector_type;
1691 
1692  //typedef typename fusion::result_of::accumulate<bases_list, fusion::vector<>, ChangeElement<> >
1693  typedef typename VectorUblas<T>::range::type ct_type;
1694 
1699  template<typename BasisType>
1701  {
1702  BOOST_MPL_ASSERT_NOT( ( boost::is_same<BasisType,mpl::void_> ) );
1703  typedef boost::shared_ptr<Cont> type;
1704  };
1705  typedef typename mpl::transform<bases_list, AddOffContainer<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type container_vector_type;
1706 
1707 
1708 
1710  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
1711 
1712  static const uint16_type nDim = mesh_type::nDim;
1713  static const uint16_type nRealDim = mesh_type::nRealDim;
1714  static const bool is_composite = functionspace_type::is_composite;
1715  static const bool is_scalar = functionspace_type::is_scalar;
1716  static const bool is_vectorial = functionspace_type::is_vectorial;
1717  static const bool is_tensor2 = functionspace_type::is_tensor2;
1718  static const bool is_continuous = functionspace_type::is_continuous;
1719  static const uint16_type nComponents1 = functionspace_type::nComponents1;
1720  static const uint16_type nComponents2 = functionspace_type::nComponents2;
1721  static const uint16_type nComponents = functionspace_type::nComponents;
1722  static const uint16_type nSpaces = functionspace_type::nSpaces;
1723 
1727 
1728 
1729  typedef typename ublas::type_traits<value_type>::real_type real_type;
1730  typedef T element_type;
1731 
1732  typedef Cont super;
1733  typedef Cont container_type;
1734  typedef container_type vector_temporary_type;
1735 
1736  typedef typename mpl::if_<mpl::bool_<is_composite>,
1737  mpl::identity<boost::none_t>,
1738  mpl::identity<typename basis_0_type::polyset_type> >::type::type polyset_type;
1739 
1740  typedef typename mpl::if_<mpl::bool_<is_composite>,
1741  mpl::identity<boost::none_t>,
1742  mpl::identity<typename basis_0_type::PreCompute> >::type::type pc_type;
1743  typedef boost::shared_ptr<pc_type> pc_ptrtype;
1744  //typedef typename basis_type::polyset_type return_value_type;
1745  typedef typename functionspace_type::return_type return_type;
1746 
1747  typedef typename matrix_node<value_type>::type matrix_node_type;
1748 
1749  typedef typename mpl::if_<mpl::bool_<is_composite>,
1750  mpl::identity<boost::none_t>,
1751  mpl::identity<typename basis_0_type::polynomial_type> >::type::type polynomial_view_type;
1752 
1753  typedef Element<T,Cont> this_type;
1754  template<int i>
1755  struct sub_element
1756  {
1757  typedef typename mpl::at_c<element_vector_type,i>::type type;
1758  };
1760  typedef typename functionspace_type::component_functionspace_ptrtype component_functionspace_ptrtype;
1761  typedef typename component_functionspace_type::template Element<T,typename VectorUblas<value_type>::slice::type> component_type;
1762 
1766  typedef typename mesh_type::element_type geoelement_type;
1767  typedef typename functionspace_type::gm_type gm_type;
1768  typedef boost::shared_ptr<gm_type> gm_ptrtype;
1769  typedef typename gm_type::template Context<vm::POINT|vm::JACOBIAN|vm::HESSIAN|vm::KB, geoelement_type> gmc_type;
1770  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1771  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
1772  typedef typename gm_type::precompute_type geopc_type;
1773 
1775 
1779 
1780  Element();
1781 
1782  Element( Element const& __e );
1783 
1784  friend class FunctionSpace<A0,A1,A2,A3,A4>;
1785 
1786  Element( functionspace_ptrtype const& __functionspace,
1787  std::string const& __name = "unknown",
1788  size_type __start = 0,
1789  ComponentType __ct = NO_COMPONENT );
1790 
1791  Element( functionspace_ptrtype const& __functionspace,
1792  container_type const& __c,
1793  std::string const& __name = "unknown",
1794  size_type __start = 0,
1795  ComponentType __ct = NO_COMPONENT );
1796 
1797  ~Element();
1798 
1799  void initFromSpace( functionspace_ptrtype const& __functionspace,
1800  container_type const& __c );
1801 
1803 
1807 
1808  Element& operator=( Element const& __e );
1809 #if 0
1810  template<typename ExprT>
1811  Element& operator=( vf::Expr<ExprT> const& __expr )
1812  {
1813  *this = project( this->functionspace(), elements( this->functionspace()->mesh() ), __expr );
1814  return *this;
1815  }
1816 #endif
1817  template<typename VectorExpr>
1818  Element& operator=( VectorExpr const& __v );
1819 
1820 
1824  super const& container() const
1825  {
1826  return *this;
1827  }
1828 
1832  super & container()
1833  {
1834  return *this;
1835  }
1836 
1843  template<ComponentType THECOMP>
1844  component_type
1845  comp()
1846  {
1847  //return comp( typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1848  return this->template comp<THECOMP>( typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1849  }
1850  template<ComponentType THECOMP>
1851  component_type
1852  comp( mpl::bool_<true> )
1853  {
1854  BOOST_STATIC_ASSERT( THECOMP >= X && THECOMP < ( ComponentType )N_COMPONENTS );
1855  auto s = ublas::slice( THECOMP, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1856  std::string __name = this->name() + "_" + componentToString( THECOMP );
1857  return component_type( compSpace(),
1858  typename component_type::container_type( this->vec().data().expression(), s ),
1859  __name,
1860  start()+THECOMP,
1861  THECOMP );
1862  }
1863  template<ComponentType THECOMP>
1864  component_type
1865  comp( mpl::bool_<false> )
1866  {
1867  BOOST_STATIC_ASSERT( THECOMP >= X && THECOMP < ( ComponentType )N_COMPONENTS );
1868  auto s = ublas::slice( THECOMP, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1869  std::string __name = this->name() + "_" + componentToString( THECOMP );
1870  return component_type( compSpace(),
1871  typename component_type::container_type( ( VectorUblas<value_type>& )*this, s ),
1872  __name,
1873  start()+THECOMP,
1874  THECOMP );
1875  }
1876 
1884  component_type
1885  comp( ComponentType i ) const
1886  {
1887  //return comp( i, mpl::bool_<boost::is_same<>is_composite>() );
1888  return comp( i, typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1889  }
1890  component_type
1891  comp( ComponentType i, mpl::bool_<true> ) const
1892  {
1893  FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1894  auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1895  std::string __name = this->name() + "_" + componentToString( i );
1896 
1897  component_type c( compSpace(),
1898  typename component_type::container_type( this->vec().data().expression(), s ),
1899  __name,
1900  start()+i,
1901  i );
1902  return c;
1903  }
1904  component_type
1905  comp( ComponentType i, mpl::bool_<false> ) const
1906  {
1907  FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1908  auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1909  std::string __name = this->name() + "_" + componentToString( i );
1910 
1911  component_type c( compSpace(),
1912  typename component_type::container_type( ( VectorUblas<value_type>& )*this, s ),
1913  //typename component_type::container_type( this->data().expression(), r ),
1914  __name,
1915  start()+i,
1916  i );
1917  return c;
1918  }
1919 
1926  component_type
1927  comp( ComponentType i )
1928  {
1929  //return comp( i, mpl::bool_<is_composite>() );
1930  return comp( i, typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1931  }
1932  component_type
1933  comp( ComponentType i, mpl::bool_<true> )
1934  {
1935  FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1936  auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1937 
1938  std::string __name = this->name() + "_" + componentToString( i );
1939  component_type c( compSpace(),
1940  typename component_type::container_type( this->vec().data().expression(), s ),
1941  __name,
1942  start()+i,
1943  i );
1944  return c;
1945  }
1946  component_type
1947  comp( ComponentType i, mpl::bool_<false> )
1948  {
1949  FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1950  auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1951 
1952  std::string __name = this->name() + "_" + componentToString( i );
1953  component_type c( compSpace(),
1954  typename component_type::container_type( ( VectorUblas<value_type>& )*this, s ),
1955  __name,
1956  start()+i,
1957  i );
1958  return c;
1959  }
1960 
1961  value_type localToGlobal( size_type e, size_type l, int c ) const
1962  {
1963  size_type index=start()+boost::get<0>( M_functionspace->dof()->localToGlobal( e, l, c ) );
1964  return super::operator()( index );
1965  }
1966 #if 0
1967  value_type& localToGlobal( size_type e, size_type l, int c )
1968  {
1969  size_type index=start()+boost::get<0>( M_functionspace->dof()->localToGlobal( e, l, c ) );
1970  return super::operator()( index );
1971  }
1972 #endif
1973  value_type operator()( size_t i ) const
1974  {
1975  return super::operator()( i );
1976  }
1977  value_type& operator()( size_t i )
1978  {
1979  return super::operator()( i );
1980  }
1981  Element& operator+=( Element const& _e )
1982  {
1983  for ( int i=0; i < _e.nLocalDof(); ++i )
1984  this->operator()( i ) += _e( i );
1985 
1986  return *this;
1987  }
1988  Element& operator-=( Element const& _e )
1989  {
1990  for ( size_type i=0; i < _e.nLocalDof(); ++i )
1991  this->operator()( i ) -= _e( i );
1992 
1993  return *this;
1994  }
1998  void operator()( MESH_CHANGES mesh_changes )
1999  {
2000  DVLOG(2) << "Update element after a change in the mesh\n";
2001  }
2002 
2003  template<typename AE>
2004  container_type& assign( const ublas::vector_expression<AE> &ae )
2005  {
2006  return super::assign( ae );
2007  }
2008  void assign( size_type ie, uint16_type il, uint16_type c, value_type const& __v )
2009  {
2010  size_type index=start()+ M_functionspace->dof()->localToGlobal( ie, il, c ).index();
2011  this->operator[]( index ) = __v;
2012  }
2013  void plus_assign( size_type ie, uint16_type il, uint16_type c, value_type const& __v )
2014  {
2015  size_type index=start()+ M_functionspace->dof()->localToGlobal( ie, il, c ).index();
2016  this->operator[]( index ) += __v;
2017  }
2018 
2020 
2024 
2025  typedef boost::multi_array<value_type,3> array_type;
2026  typedef Eigen::Matrix<value_type,nComponents1,1> _id_type;
2027  typedef Eigen::Matrix<value_type,nComponents1,nRealDim> _grad_type;
2028  typedef Eigen::Matrix<value_type,nRealDim,nRealDim> _hess_type;
2029  typedef Eigen::Matrix<value_type,1,1> _div_type;
2030  typedef Eigen::Matrix<value_type,nRealDim,1> _curl_type;
2031  typedef boost::multi_array<_id_type,1> id_array_type;
2032  typedef boost::multi_array<_grad_type,1> grad_array_type;
2033  typedef boost::multi_array<_hess_type,1> hess_array_type;
2034  typedef boost::multi_array<_div_type,1> div_array_type;
2035  typedef boost::multi_array<_curl_type,1> curl_array_type;
2036  typedef boost::multi_array<_div_type,1> comp_curl_array_type;
2037 
2041  DataMap const& map() const
2042  {
2043  return M_functionspace->map();
2044  }
2045 
2049  mesh_ptrtype mesh()
2050  {
2051  return M_functionspace->mesh();
2052  }
2053 
2057  mesh_ptrtype mesh() const
2058  {
2059  return M_functionspace->mesh();
2060  }
2061 
2065  Element<T,Cont> sqrt() const
2066  {
2067  Element<T,Cont> _e( M_functionspace );
2068 
2069  for ( int i=0; i < _e.nLocalDof(); ++i ) _e( i ) = math::sqrt( this->operator()( i ) );
2070 
2071  return _e;
2072  }
2073 
2077  Element<T,Cont> pow( int n ) const
2078  {
2079  Element<T,Cont> _e( M_functionspace );
2080 
2081  for ( int i=0; i < _e.nLocalDof(); ++i ) _e( i ) = math::pow( this->operator()( i ),n );
2082 
2083  return _e;
2084  }
2085 
2086  // Only works for scalar fields
2087  template < typename p0_space_type >
2088  typename p0_space_type::element_type extremeValue( boost::shared_ptr<p0_space_type> const& P0h, std::string extreme )
2089  {
2090  // check if the mesh coming from P0h and the class elements is the same
2091  FEELPP_ASSERT( P0h->mesh() == this->mesh() ).error( "mesh is not the same" );
2092  FEELPP_ASSERT( is_scalar ).error( "only works for scalar fields" );
2093 
2094  typename p0_space_type::element_type p0Element( P0h );
2095 
2096  for ( auto elt_it = P0h->mesh()->beginElement(); elt_it != P0h->mesh()->endElement(); ++elt_it )
2097  {
2098  size_type eid = elt_it->id();
2099 
2100  size_type dofp0 = boost::get<0>( P0h->dof()->localToGlobal( eid, 0, 0 ) );
2101  std::vector<value_type> values ( functionspace_type::fe_type::nLocalDof );
2102 
2103  size_type dofpn = 0;
2104 
2105  for ( uint16_type local_id=0; local_id < functionspace_type::fe_type::nLocalDof; ++local_id )
2106  {
2107  dofpn = boost::get<0>( this->functionSpace()->dof()->localToGlobal( eid, local_id, 0 ) );
2108  values[local_id] = this->operator()( dofpn );
2109  }
2110 
2111  if ( extreme == "max" )
2112  p0Element.assign( eid, 0, 0, *( std::max_element( values.begin(), values.end() ) ) );
2113 
2114  else
2115  p0Element.assign( eid, 0, 0, *( std::min_element( values.begin(), values.end() ) ) );
2116  }
2117 
2118  return p0Element;
2119  }
2120 
2121  value_type max() const
2122  {
2123  return super::max();
2124  }
2125 
2126  template < typename p0_space_type >
2127  typename p0_space_type::element_type max( boost::shared_ptr<p0_space_type> const& P0h )
2128  {
2129  return this->extremeValue( P0h, "max" );
2130  }
2131 
2132  value_type min() const
2133  {
2134  return super::min();
2135  }
2136 
2137  template < typename p0_space_type >
2138  typename p0_space_type::element_type min( boost::shared_ptr<p0_space_type> const& P0h )
2139  {
2140  return this->extremeValue( P0h, "min" );
2141  }
2142 
2144 
2145 
2149  typedef detail::ID<value_type,nComponents1,nComponents2> id_type;
2150 
2151 
2156  template<typename ContextType>
2157  boost::array<typename array_type::index, 1>
2158  idExtents( ContextType const & context ) const
2159  {
2160  boost::array<typename array_type::index, 1> shape;
2161  shape[0] = context.xRefs().size2();
2162  //shape[1] = nComponents1;
2163  //shape[2] = nComponents2;
2164  return shape;
2165  }
2166 
2167  template<typename Context_t>
2168  id_type
2169  id( Context_t const & context ) const
2170  {
2171  return id_type( *this, context );
2172  }
2173 
2185  template<typename Context_t>
2186  void
2187  id_( Context_t const & context, id_array_type& v ) const;
2188 
2189  template<typename Context_t>
2190  void
2191  id( Context_t const & context, id_array_type& v ) const
2192  {
2193  id_( context, v );
2194  }
2195 
2196 
2197  void
2198  idInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const;
2199 
2200  /*
2201  * Get the reals points matrix in a context
2202  * 1 : Element
2203  *
2204  * \todo store a geometric mapping context to evaluate the real points
2205  * from a set of point in the referene element, should probably done in
2206  * the real element (geond)
2207  */
2208  template<typename Context_t>
2209  matrix_node_type
2210  ptsInContext( Context_t const & context, mpl::int_<1> ) const
2211  {
2212  //new context for evaluate the points
2213  typedef typename Context_t::gm_type::template Context< Context_t::context|vm::POINT, typename Context_t::element_type> gmc_interp_type;
2214  typedef boost::shared_ptr<gmc_interp_type> gmc_interp_ptrtype;
2215 
2216  gmc_interp_ptrtype __c_interp( new gmc_interp_type( context.geometricMapping(), context.element_c(), context.pc() ) );
2217 
2218  return __c_interp->xReal();
2219  }
2220 
2221  /*
2222  * Get the real point matrix in a context
2223  * 2 : Face
2224  * \todo see above
2225  */
2226  template<typename Context_t>
2227  matrix_node_type
2228  ptsInContext( Context_t const & context, mpl::int_<2> ) const
2229  {
2230  //new context for the interpolation
2231  typedef typename Context_t::gm_type::template Context< Context_t::context|vm::POINT, typename Context_t::element_type> gmc_interp_type;
2232  typedef boost::shared_ptr<gmc_interp_type> gmc_interp_ptrtype;
2233 
2234  typedef typename Context_t::gm_type::template Context<Context_t::context,typename Context_t::element_type>::permutation_type permutation_type;
2235  typedef typename Context_t::gm_type::template Context<Context_t::context,typename Context_t::element_type>::precompute_ptrtype precompute_ptrtype;
2236 
2237  //not good because ?
2238  //gmc_interp_ptrtype __c_interp( new gmc_interp_type( context.geometricMapping(), context.element_c(),context.pcFaces(), context.faceId()) );
2239  //good with this
2240  std::vector<std::map<permutation_type, precompute_ptrtype> > __geo_pcfaces = context.pcFaces();
2241  gmc_interp_ptrtype __c_interp( new gmc_interp_type( context.geometricMapping(), context.element_c(), __geo_pcfaces , context.faceId() ) );
2242 
2243  return __c_interp->xReal();
2244  }
2245 
2250  gmc_ptrtype geomapPtr( geoelement_type const& elt ) const
2251  {
2252  gm_ptrtype __gm = functionSpace()->gm();
2253  geopc_ptrtype __geopc( new geopc_type( __gm, elt.G() ) );
2254  return gmc_ptrtype( new gmc_type( __gm, elt, __geopc ) );
2255  }
2260  pc_ptrtype pcPtr( geoelement_type const& elt ) const
2261  {
2262  return pc_ptrtype( new pc_type( functionSpace()->fe(), elt.G() ) );
2263  }
2264 
2265  id_type operator()( Eigen::Matrix<value_type,nDim,1> const& __x, bool extrapolate = false ) const
2266  {
2267  node_type n( nDim );
2268  for(int i = 0; i < nDim; ++i ) n[i]=__x[i];
2269  return operator()( n, extrapolate );
2270  }
2271 
2277  id_type operator()( node_type const& __x, bool extrapolate = false ) const;
2278 
2280 
2282 
2283 
2284  typedef detail::DD<value_type, nComponents1, nRealDim> grad_type;
2285  typedef detail::D<value_type,0,nComponents1, 1> dx_type;
2286  typedef detail::D<value_type,1,nComponents1, 1> dy_type;
2287  typedef detail::D<value_type,2,nComponents1, 1> dz_type;
2288 
2289  template<typename ContextType>
2290  boost::array<typename array_type::index, 1>
2291  gradExtents( ContextType const & context ) const
2292  {
2293  boost::array<typename array_type::index, 1> shape;
2294  shape[0] = context.xRefs().size2();
2295  //shape[1] = nComponents1;
2296  //shape[2] = nRealDim;
2297 
2298  return shape;
2299  }
2300  template<typename ContextType>
2301  grad_type
2302  grad( ContextType const & context ) const
2303  {
2304  return grad_type( *this, context );
2305  }
2306  template<typename ElementType>
2307  void
2308  grad( ElementType const& elt, grad_array_type& v )
2309  {
2310  gmc_ptrtype gmc( geomapPtr( elt ) );
2311  v.resize( gradExtents( *gmc ) );
2312  grad_( *gmc, *pcPtr( elt ), v );
2313  }
2314 
2329  template<typename ContextType>
2330  void grad_( ContextType const & context, grad_array_type& v ) const;
2331 
2332  template<typename ContextType>
2333  void grad( ContextType const & context, grad_array_type& v ) const
2334  {
2335  grad_( context, v );
2336  }
2337 
2338  void
2339  gradInterpolate( matrix_node_type __ptsReal, grad_array_type& v ) const;
2340 
2346  grad_type grad( node_type const& __x ) const;
2347 
2348  template<typename ContextType>
2349  dx_type
2350  dx( ContextType const & context ) const
2351  {
2352  return dx_type( *this, context );
2353  }
2354 
2355  template<typename ContextType>
2356  dy_type
2357  dy( ContextType const & context ) const
2358  {
2359  return dy_type( *this, context );
2360  }
2361  template<typename ContextType>
2362  dz_type
2363  dz( ContextType const & context ) const
2364  {
2365  return dz_type( *this, context );
2366  }
2367 
2368  template<typename ContextType>
2369  boost::array<typename array_type::index, 1>
2370  dxExtents( ContextType const & context ) const
2371  {
2372  boost::array<typename array_type::index, 1> shape;
2373  shape[0] = context.xRefs().size2();
2374  //shape[1] = nComponents1;
2375  //shape[2] = nComponents2;
2376  return shape;
2377  }
2378  template<typename ContextType>
2379  boost::array<typename array_type::index, 1>
2380  dyExtents( ContextType const & context ) const
2381  {
2382  return dxExtents( context );
2383  }
2384  template<typename ContextType>
2385  boost::array<typename array_type::index, 1>
2386  dzExtents( ContextType const & context ) const
2387  {
2388  return dxExtents( context );
2389  }
2390  template<typename ContextType>
2391  boost::array<typename array_type::index, 1>
2392  dExtents( ContextType const & context ) const
2393  {
2394  return dxExtents( context );
2395  }
2396 
2397  template<typename ContextType>
2398  void d_( int N, ContextType const & context, id_array_type& v ) const;
2399 
2400  template<typename ContextType>
2401  void dx( ContextType const & context, id_array_type& v ) const
2402  {
2403  d_( 0, context, v );
2404  }
2405  template<typename ContextType>
2406  void dy( ContextType const & context, id_array_type& v ) const
2407  {
2408  d_( 1, context, v );
2409  }
2410  template<typename ContextType>
2411  void dz( ContextType const & context, id_array_type& v ) const
2412  {
2413  d_( 2, context, v );
2414  }
2415 
2416  void
2417  dxInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const;
2418 
2419  void
2420  dyInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const;
2421 
2422  void
2423  dzInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const;
2424 
2425 
2427 
2428  typedef detail::Div<value_type> div_type;
2429 
2430  template<typename ContextType>
2431  boost::array<typename array_type::index, 1>
2432  divExtents( ContextType const & context ) const
2433  {
2434  boost::array<typename array_type::index, 1> shape;
2435  shape[0] = context.xRefs().size2();
2436  //shape[1] = 1;
2437  //shape[2] = 1;
2438  return shape;
2439  }
2440  template<typename ContextType>
2441  div_type
2442  div( ContextType const & context ) const
2443  {
2444  //BOOST_STATIC_ASSERT( (rank == 1) || (rank==2) );
2445  return div_type( *this, context );
2446  }
2447 
2448 
2449  template<typename ContextType>
2450  void
2451  div( ContextType const & context, div_array_type& v ) const
2452  {
2453  //BOOST_STATIC_ASSERT( (rank == 1) || (rank==2) );
2454  div_( context, v );
2455  }
2456  template<typename ContextType>
2457  void div_( ContextType const & context, div_array_type& v ) const;
2458 
2459  void
2460  divInterpolate( matrix_node_type __ptsReal, div_array_type& v ) const;
2461 
2462  typedef detail::Curl<value_type,-1,nRealDim> curl_type;
2463  typedef detail::Curl<value_type,0,1> curlx_type;
2464  typedef detail::Curl<value_type,1,1> curly_type;
2465  typedef detail::Curl<value_type,2,1> curlz_type;
2466 
2467  template<typename ContextType>
2468  curl_type
2469  curl( ContextType const & context ) const
2470  {
2471  return curl_type( *this, context );
2472  }
2473 
2474  template<typename ContextType>
2475  curlx_type
2476  curlx( ContextType const & context ) const
2477  {
2478  return curlx_type( *this, context );
2479  }
2480 
2481  template<typename ContextType>
2482  curly_type
2483  curly( ContextType const & context ) const
2484  {
2485  return curly_type( *this, context );
2486  }
2487 
2488  template<typename ContextType>
2489  curlz_type
2490  curlz( ContextType const & context ) const
2491  {
2492  return curlz_type( *this, context );
2493  }
2494 
2495  template<typename ContextType>
2496  boost::array<typename array_type::index, 1>
2497  curlExtents( ContextType const & context ) const
2498  {
2499  //BOOST_MPL_ASSERT_MSG( ( rank == 1 ), INVALID_RANK_FOR_CURL, (mpl::int_<rank>) );
2500 
2501  boost::array<typename array_type::index, 1> shape;
2502  shape[0] = context.xRefs().size2();
2503  //shape[1] = nComponents1;
2504  //shape[2] = 1;
2505  return shape;
2506  }
2507  template<typename ContextType>
2508  void
2509  curl( ContextType const & context, curl_array_type& v ) const
2510  {
2511  //BOOST_MPL_ASSERT_MSG( ( rank == 1 ), INVALID_RANK_FOR_CURL, (mpl::int_<rank>) );
2512 
2513  curl_( context, v );
2514  }
2515 
2516  void
2517  curlInterpolate( matrix_node_type __ptsReal, curl_array_type& v ) const;
2518 
2519  template<typename ContextType>
2520  void curl_( ContextType const & context, curl_array_type& v ) const;
2521 
2522  template<typename ContextType>
2523  void curl_( ContextType const & context, comp_curl_array_type& v, int comp ) const;
2524 
2525 
2526  template<typename ContextType>
2527  boost::array<typename array_type::index, 1>
2528  curlxExtents( ContextType const & context ) const
2529  {
2530  boost::array<typename array_type::index, 1> shape;
2531  shape[0] = context.xRefs().size2();
2532  //shape[1] = nComponents1;
2533  //shape[2] = nComponents2;
2534  return shape;
2535  }
2536  template<typename ContextType>
2537  boost::array<typename array_type::index, 1>
2538  curlyExtents( ContextType const & context ) const
2539  {
2540  return curlxExtents( context );
2541  }
2542  template<typename ContextType>
2543  boost::array<typename array_type::index, 1>
2544  curlzExtents( ContextType const & context ) const
2545  {
2546  return curlxExtents( context );
2547  }
2548 
2549  template<typename ContextType>
2550  void
2551  curlx( ContextType const & context, comp_curl_array_type& v ) const
2552  {
2553  //BOOST_STATIC_ASSERT( rank == 1 );
2554  curl_( context, v, 0 );
2555  }
2556  template<typename ContextType>
2557  void
2558  curly( ContextType const & context, comp_curl_array_type& v ) const
2559  {
2560  //BOOST_STATIC_ASSERT( rank == 1 );
2561  curl_( context, v, 1 );
2562  }
2563  template<typename ContextType>
2564  void
2565  curlz( ContextType const & context, comp_curl_array_type& v ) const
2566  {
2567  //BOOST_STATIC_ASSERT( rank == 1 );
2568  curl_( context, v, 2 );
2569  }
2570 
2571  void
2572  curlxInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const;
2573 
2574  void
2575  curlyInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const;
2576 
2577  void
2578  curlzInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const;
2579 
2580  template<typename ContextType>
2581  boost::array<typename array_type::index, 1>
2582  hessExtents( ContextType const & context ) const
2583  {
2584  BOOST_STATIC_ASSERT( rank == 0 );
2585 
2586  boost::array<typename array_type::index, 1> shape;
2587  shape[0] = context.xRefs().size2();
2588  //shape[1] = nRealDim;
2589  //shape[2] = nRealDim;
2590  return shape;
2591  }
2605  template<typename ContextType>
2606  void
2607  hess_( ContextType const & context, hess_array_type& v ) const;
2608 
2609  template<typename ContextType>
2610  void
2611  hess_( ContextType const & context, hess_array_type& v, mpl::int_<0> ) const;
2612 
2613  void
2614  hessInterpolate( matrix_node_type __ptsReal, hess_array_type& v ) const;
2615 
2616  typedef detail::H<value_type,nRealDim,nRealDim> hess_type;
2617 
2618  template<typename ContextType>
2619  hess_type
2620  hess( ContextType const & context ) const
2621  {
2622  return hess_type( *this, context );
2623  }
2624 
2625  template<typename ContextType>
2626  void
2627  hess( ContextType const & context, hess_array_type& v ) const
2628  {
2629  hess_( context, v );
2630  }
2631  template<typename ElementType>
2632  void
2633  hess( ElementType const& elt, hess_array_type& v )
2634  {
2635  gmc_ptrtype gmc( geomapPtr( elt ) );
2636  v.resize( hessExtents( *gmc ) );
2637  hess_( *gmc, *pcPtr( elt ), v );
2638  }
2639 
2643  functionspace_ptrtype const& functionSpace() const
2644  {
2645  return M_functionspace;
2646  }
2647 
2651  //functionspace_vector_type const&
2652  //functionSpaces() const { return M_functionspace->functionSpaces(); }
2653 
2654  functionspace_vector_type const&
2655  functionSpaces() const
2656  {
2657  return M_functionspace->functionSpaces();
2658  }
2659 
2663  template<int i>
2664  typename mpl::at_c<functionspace_vector_type,i>::type
2665  functionSpace() const
2666  {
2667  return M_functionspace->template functionSpace<i>();
2668  }
2669 
2670  Eigen::Matrix<value_type, Eigen::Dynamic, 1>
2671  extractValuesWithMarker( std::string const& m )
2672  {
2673  auto r = functionSpace()->dof()->markerToDof( m );
2674  Eigen::Matrix<value_type, Eigen::Dynamic, 1> res( std::distance( r.first, r.second ) );
2675  size_type i = 0;
2676  for( auto it = r.first, en = r.second; it != en; ++it, ++i )
2677  res( i ) = this->operator[]( it->second );
2678  return res;
2679  }
2680  Eigen::Matrix<value_type, Eigen::Dynamic, 1>
2681  extractValuesWithoutMarker( std::string const& m )
2682  {
2683  auto r1 = functionSpace()->dof()->markerToDofLessThan( m );
2684  auto r2 = functionSpace()->dof()->markerToDofGreaterThan( m );
2685  size_type s = std::distance( r1.first, r1.second ) + std::distance( r2.first, r2.second );
2686  Eigen::Matrix<value_type, Eigen::Dynamic, 1> res( s );
2687  size_type i = 0;
2688  for( auto it = r1.first, en = r1.second; it != en; ++it, ++i )
2689  res( i ) = this->operator[]( it->second );
2690  for( auto it = r2.first, en = r2.second; it != en; ++it, ++i )
2691  res( i ) = this->operator[]( it->second );
2692  return res;
2693  }
2694 
2695  template<int i>
2696  typename mpl::at_c<element_vector_type,i>::type
2697  element( std::string const& name ="u", bool updateOffViews=true )
2698  {
2699  if ( this->worldComm().globalSize()>1 ) this->worldComm().globalComm().barrier();
2700 
2701  size_type nbdof_start = fusion::accumulate( this->functionSpaces(),
2702  size_type( 0 ),
2703  detail::NLocalDof<mpl::bool_<true> >( this->worldsComm(), false, 0, i ) );
2704 
2705  typename mpl::at_c<functionspace_vector_type,i>::type space( M_functionspace->template functionSpace<i>() );
2706  DVLOG(2) << "Element <" << i << ">::start : "<< nbdof_start << "\n";
2707  DVLOG(2) << "Element <" << i << ">::size : "<< space->nDof()<< "\n";
2708  DVLOG(2) << "Element <" << i << ">::local size : "<< space->nLocalDof()<< "\n";
2709  DVLOG(2) << "Element <" << -1 << ">::size : "<< this->size() << "\n";
2710 
2711  if ( this->functionSpace()->template functionSpace<i>()->worldComm().isActive() )
2712  {
2713  ct_type ct( *this, ublas::range( nbdof_start, nbdof_start+space->nLocalDof() ),
2714  M_functionspace->template functionSpace<i>()->dof() );
2715 
2716 #if defined(FEELPP_ENABLE_MPI_MODE)
2717 
2718  // update M_containersOffProcess<i> : send
2719  if ( this->worldComm().globalSize()>1 && updateOffViews && !this->functionSpace()->hasEntriesForAllSpaces() )
2720  {
2721  std::vector<double> dataToSend( space->nLocalDof() );
2722  std::copy( ct.begin(), ct.end(), dataToSend.begin() );
2723 
2724  if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2725 
2726  fusion::for_each( *M_containersOffProcess, detail::SendContainersOn<i,functionspace_type>( this->functionSpace(), dataToSend ) );
2727  }
2728 
2729 #endif
2730  DVLOG(2) << "Element <" << i << ">::range.size : "<< ct.size()<< "\n";
2731  DVLOG(2) << "Element <" << i << ">::range.start : "<< ct.start()<< "\n";
2732  return typename mpl::at_c<element_vector_type,i>::type( space, ct, name );
2733  }
2734 
2735  else
2736  {
2737  // initialize if not the case
2738  if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2739 
2740  fusion::for_each( *M_containersOffProcess, detail::InitializeContainersOff<i,functionspace_type>( this->functionSpace() ) );
2741 
2742 #if defined(FEELPP_ENABLE_MPI_MODE)
2743 
2744  // update M_containersOffProcess<i> : recv
2745  if ( this->worldComm().globalSize()>1 && updateOffViews && !this->functionSpace()->hasEntriesForAllSpaces() )
2746  {
2747  fusion::for_each( *M_containersOffProcess, detail::RecvContainersOff<i,functionspace_type>( this->functionSpace() ) );
2748  }
2749 
2750 #endif
2751  // build a subrange view identical
2752  ct_type ct( *fusion::at_c<i>( *M_containersOffProcess ),
2753  ublas::range( 0, space->nLocalDof() ),
2754  M_functionspace->template functionSpace<i>()->dof() );
2755 
2756  DVLOG(2) << "Element <" << i << ">::range.size : "<< ct.size()<< "\n";
2757  DVLOG(2) << "Element <" << i << ">::range.start : "<< ct.start()<< "\n";
2758 
2759  return typename mpl::at_c<element_vector_type,i>::type( space, ct, name );
2760  }
2761  }
2762 
2763  template<int i>
2764  typename mpl::at_c<element_vector_type,i>::type
2765  element( std::string const& name ="u", bool updateOffViews=true ) const
2766  {
2767  size_type nbdof_start = fusion::accumulate( M_functionspace->functionSpaces(),
2768  size_type( 0 ),
2769  detail::NLocalDof<mpl::bool_<true> >( this->worldsComm(), false, 0, i ) );
2770  typename mpl::at_c<functionspace_vector_type,i>::type space( M_functionspace->template functionSpace<i>() );
2771 
2772  DVLOG(2) << "Element <" << i << ">::start : "<< nbdof_start << "\n";
2773  DVLOG(2) << "Element <" << i << ">::size : "<< space->nDof()<< "\n";
2774  DVLOG(2) << "Element <" << i << ">::local size : "<< space->nLocalDof()<< "\n";
2775  DVLOG(2) << "Element <" << -1 << ">::size : "<< this->size() << "\n";
2776 
2777  if ( this->functionSpace()->worldsComm()[i].isActive() )
2778  {
2779  ct_type ct( const_cast<VectorUblas<value_type>&>( dynamic_cast<VectorUblas<value_type> const&>( *this ) ),
2780  ublas::range( nbdof_start, nbdof_start+space->nLocalDof() ),
2781  M_functionspace->template functionSpace<i>()->dof() );
2782 
2783 #if defined(FEELPP_ENABLE_MPI_MODE)
2784 
2785  // update M_containersOffProcess<i> : send
2786  if ( this->worldComm().globalSize()>1 && updateOffViews && !this->functionSpace()->hasEntriesForAllSpaces() )
2787  {
2788  std::vector<double> dataToSend( space->nLocalDof() );
2789  std::copy( ct.begin(), ct.end(), dataToSend.begin() );
2790 
2791  if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2792 
2793  fusion::for_each( *M_containersOffProcess, detail::SendContainersOn<i,functionspace_type>( this->functionSpace(), dataToSend ) );
2794  }
2795 
2796 #endif
2797 
2798  DVLOG(2) << "Element <" << i << ">::range.size : "<< ct.size()<< "\n";
2799  DVLOG(2) << "Element <" << i << ">::range.start : "<< ct.start()<< "\n";
2800  return typename mpl::at_c<element_vector_type,i>::type( space, ct, name );
2801  }
2802 
2803  else
2804  {
2805 #if defined(FEELPP_ENABLE_MPI_MODE)
2806 
2807  // update M_containersOffProcess<i> : recv
2808  if ( this->worldComm().globalSize()>1 && updateOffViews && !this->functionSpace()->hasEntriesForAllSpaces() )
2809  {
2810  fusion::for_each( *M_containersOffProcess, detail::RecvContainersOff<i,functionspace_type>( this->functionSpace() ) );
2811  }
2812 
2813 #endif
2814 
2815  // build a subrange view identical
2816  ct_type ct( *fusion::at_c<i>( *M_containersOffProcess ),
2817  ublas::range( 0, space->nLocalDof() ),
2818  M_functionspace->template functionSpace<i>()->dof() );
2819 
2820  DVLOG(2) << "Element <" << i << ">::range.size : "<< ct.size()<< "\n";
2821  DVLOG(2) << "Element <" << i << ">::range.start : "<< ct.start()<< "\n";
2822  return typename mpl::at_c<element_vector_type,i>::type( space, ct, name );
2823  }
2824 
2825  }
2830  component_functionspace_ptrtype const& compSpace() const
2831  {
2832  return M_functionspace->compSpace();
2833  }
2834 
2838  std::vector<WorldComm> const& worldsComm() const
2839  {
2840  return M_functionspace->worldsComm();
2841  }
2842 
2846  WorldComm const& worldComm() const
2847  {
2848  return M_functionspace->worldComm();
2849  }
2850 
2854  size_type nDof() const
2855  {
2856  return M_functionspace->nDof();
2857  }
2858 
2862  size_type nLocalDof() const
2863  {
2864  return M_functionspace->nLocalDof();
2865  }
2866 
2870  size_type nDofPerComponent() const
2871  {
2872  return M_functionspace->nDofPerComponent();
2873  }
2874 
2878  std::string const& name() const
2879  {
2880  return M_name;
2881  }
2882 
2883  size_type start() const
2884  {
2885  return M_start;
2886  }
2887 
2888  bool isAComponent() const
2889  {
2890  return M_ct >= X && M_ct <= Z;
2891  }
2892 
2893  ComponentType component() const
2894  {
2895  if ( M_ct < X || M_ct > Z )
2896  {
2897  std::ostringstream __str;
2898  __str << "invalid component extraction (should be 0 or 1 or 2): " << M_ct;
2899  throw std::logic_error( __str.str() );
2900  }
2901 
2902  return M_ct;
2903  }
2904 
2905  std::string componentToString( ComponentType ct ) const
2906  {
2907  if ( ct < X || ct > Z )
2908  {
2909  std::ostringstream __str;
2910  __str << "invalid component extraction (should be 0 or 1 or 2): " << ct;
2911  throw std::logic_error( __str.str() );
2912  }
2913 
2914  switch ( ct )
2915  {
2916  case X:
2917  return "X";
2918 
2919  case Y:
2920  return "Y";
2921 
2922  case Z:
2923  return "Z";
2924 
2925  default:
2926  std::ostringstream __str;
2927  __str << "invalid component extraction (should be 0 or 1 or 2): " << ct;
2928  throw std::logic_error( __str.str() );
2929  }
2930  }
2931 
2933 
2937 
2942  void setName( std::string const & __name )
2943  {
2944  M_name = __name;
2945  }
2946 
2947  void setFunctionSpace( functionspace_ptrtype space )
2948  {
2949  M_functionspace = space;
2950  super::init( M_functionspace->dof() );
2951  //super::init( M_functionspace->nDof(), M_functionspace->nLocalDof() );
2952  }
2953 
2954  BOOST_PARAMETER_MEMBER_FUNCTION( ( void ),
2955  save,
2956  tag,
2957  ( required
2958  ( path,* ) )
2959  ( optional
2960  ( type,( std::string ),std::string( "binary" ) )
2961  ( suffix,( std::string ),std::string( "" ) )
2962  ( sep,( std::string ),std::string( "" ) )
2963  ) )
2964  {
2965  Feel::detail::ignore_unused_variable_warning( args );
2966 
2967  if ( !fs::exists( fs::path( path ) ) )
2968  {
2969  fs::create_directories( fs::path( path ) );
2970  }
2971 
2972  std::ostringstream os1;
2973  os1 << M_name << sep << suffix << "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank() << ".fdb";
2974  fs::path p = fs::path( path ) / os1.str();
2975  LOG(INFO) << "saving " << p << "\n";
2976  fs::ofstream ofs( p );
2977 
2978  if ( type == "binary" )
2979  {
2980  boost::archive::binary_oarchive oa( ofs );
2981  oa << *this;
2982  }
2983 
2984  else if ( type == "text" )
2985  {
2986  boost::archive::text_oarchive oa( ofs );
2987  oa << *this;
2988  }
2989 
2990  else if ( type == "xml" )
2991  {
2992  //boost::archive::xml_oarchive oa(ofs);
2993  //oa << *this;
2994  }
2995  }
2996  BOOST_PARAMETER_MEMBER_FUNCTION(
2997  ( bool ),
2998  load,
2999  tag,
3000  ( required
3001  ( path,* ) )
3002  ( optional
3003  ( type,( std::string ),std::string( "binary" ) )
3004  ( suffix,( std::string ),std::string( "" ) )
3005  ( sep,( std::string ),std::string( "" ) )
3006  )
3007  )
3008  {
3009  Feel::detail::ignore_unused_variable_warning( args );
3010  std::ostringstream os1;
3011  os1 << M_name << sep << suffix << "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank() << ".fdb";
3012  fs::path p = fs::path( path ) / os1.str();
3013  fs::path partial_path = fs::path(path);
3014 
3015  fs::path full_path_dir_sol(fs::current_path());
3016  full_path_dir_sol = full_path_dir_sol/partial_path;
3017  //std::cout << " In load the first full path is " << p << std::endl;
3018  if ( !fs::exists( p ) )
3019  {
3020  std::ostringstream os2;
3021  os2 << M_name << sep << suffix<< "-" << this->worldComm().globalSize() << "." << this->worldComm().globalRank();
3022  p = fs::path( path ) / os2.str();
3023 
3024  if ( !fs::exists( p ) )
3025  {
3026  LOG(WARNING) << "[load] :" << full_path_dir_sol << " FILE : " << os1.str() << " OR " << os2.str() << " DO NOT EXIST" << std::endl ;
3027  //std::cerr << "ATTENTION : p does not exist
3028  return 0;
3029  }
3030  }
3031  LOG(INFO) << p << " exists, is is a regular file : " << fs::is_regular_file( p ) << "\n";
3032  if ( !fs::is_regular_file( p ) )
3033  {
3034 
3035  LOG(WARNING) << "[load] : " << full_path_dir_sol << p << " is not a regular_file !" << std::endl;
3036  return 0;
3037  }
3038 
3039  fs::ifstream ifs( p );
3040 
3041  if ( type == "binary" )
3042  {
3043  boost::archive::binary_iarchive ia( ifs );
3044  ia >> *this;
3045  }
3046 
3047  else if ( type == "text" )
3048  {
3049  boost::archive::text_iarchive ia( ifs );
3050  ia >> *this;
3051  }
3052 
3053  else if ( type == "xml" )
3054  {
3055  //boost::archive::xml_iarchive ia(ifs);
3056  //ia >> *this;
3057  return false;
3058  }
3059  return true;
3060  }
3061 
3062  void printMatlab( std::string fname, bool gmsh = false ) const
3063  {
3064  container_type m( *this );
3065  if ( gmsh )
3066  {
3067  auto relation = this->functionSpace()->dof()->pointIdToDofRelation();
3068  for( size_type i = 0; i < this->localSize(); ++i )
3069  {
3070  m[relation.first[i]-1] = this->operator[]( i );
3071  }
3072  }
3073  m.printMatlab( fname );
3074  }
3075 
3077  private:
3078 
3079  friend class boost::serialization::access;
3080 
3081  template<class Archive>
3082  void serialize( Archive & ar, const unsigned int version )
3083  {
3084  //ar & BOOST_SERIALIZATION_NVP( boost::serialization::base_object<super>(*this) );
3085  ar & boost::serialization::make_nvp( "name", M_name );
3086  DVLOG(2) << "got name " << M_name << "\n";
3087 
3088  if ( Archive::is_saving::value )
3089  {
3090  //std::cout << "saving in version " << version << "\n";
3091  size_type s = this->functionSpace()->nLocalDofWithGhost();
3092  ar & boost::serialization::make_nvp( "size", s );
3093 
3094  std::vector<int> no = M_functionspace->basisOrder();
3095  //for( int i = 0; i < no.size(); ++i ) std::cout << no[i] << std::endl;
3096  std::string family = M_functionspace->basisName();
3097  //std::cout << "family name = " << family << std::endl;
3098 
3099  ar & boost::serialization::make_nvp( "order", no );
3100  //std::cout << "saving order done" << std::endl;
3101  ar & boost::serialization::make_nvp( "family", family );
3102  //std::cout << "saving family done" << std::endl;
3103 
3104  typename container_type::const_iterator it = this->begin();
3105  typename container_type::const_iterator en = this->end();
3106 
3107  for ( size_type i = 0; it != en; ++it, ++i )
3108  {
3109  T value = this->operator[]( i );
3110  std::ostringstream v_str;
3111  v_str << "value_" << i;
3112 
3113  // DVLOG(2) << "save value " << value << " at " << v_str.str() << "\n";
3114 
3115  ar & boost::serialization::make_nvp( v_str.str().c_str(), value );
3116  }
3117  }
3118 
3119  if ( Archive::is_loading::value )
3120  {
3121  //std::cout << "loading in version " << version << "\n";
3122 
3123  size_type s( 0 );
3124  ar & boost::serialization::make_nvp( "size", s );
3125 
3126  // verify number of degree of freedom
3127  DVLOG(2) << "loading ublas::vector of size " << s << "\n";
3128 
3129  if ( s != this->functionSpace()->nLocalDofWithGhost() )
3130  throw std::logic_error( ( boost::format( "load function: invalid number of degrees of freedom, read %1% but has %2%" ) % s % this->functionSpace()->nLocalDofWithGhost() ).str() );
3131 
3132 
3133  std::vector<int> order;
3134  std::string family;
3135  ar & boost::serialization::make_nvp( "order", order );
3136  //for( int i = 0; i < order.size(); ++i ) std::cout << order[i] << std::endl;
3137 
3138  ar & boost::serialization::make_nvp( "family", family );
3139  //std::cout << "family name = " << family << std::endl;
3140 #if 0
3141  auto orders = M_functionspace->basisOrder();
3142 
3143  if ( order != orders )
3144  throw std::logic_error( ( boost::format( "load function: invalid polynomial order, read %1% but has %2%" ) % order % orders ).str() );
3145 
3146  std::string bname = M_functionspace->basisName();
3147 
3148  if ( family != bname )
3149  throw std::logic_error( ( boost::format( "load function: invalid polynomial family, read %1% but has %2%" ) % family % bname ).str() );
3150 
3151 #endif
3152 
3153  for ( size_type i = 0; i < s ; ++i )
3154  {
3155  value_type value( 0 );
3156  std::ostringstream v_str;
3157  v_str << "value_" << i;
3158  // DVLOG(2) << "load value at " << v_str.str() << "\n";
3159  ar & boost::serialization::make_nvp( v_str.str().c_str(), value );
3160  // DVLOG(2) << "got value " << value << " at index " << i << "\n";
3161  this->operator[]( i ) = value;
3162  }
3163  }
3164 
3165 
3166  }
3167 
3168  private:
3169 
3173  functionspace_ptrtype M_functionspace;
3174 
3175  std::string M_name;
3176 
3177  size_type M_start;
3178 
3180  ComponentType M_ct;
3181 
3182  // only init in // with composite case : ex p = U.element<1>()
3183  mutable boost::optional<container_vector_type> M_containersOffProcess;
3184 
3185  }; // Element
3186 
3188 
3191  typedef Element<value_type> element_type;
3192  typedef Element<value_type> real_element_type;
3193  typedef boost::shared_ptr<element_type> element_ptrtype;
3194 
3195  typedef std::map< size_type, std::vector< size_type > > proc_dist_map_type;
3196 
3198 
3202 
3208  FunctionSpace( mesh_ptrtype const& mesh,
3209  size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3210  periodicity_type periodicity = periodicity_type(),
3211  std::vector<WorldComm> const& _worldsComm = Environment::worldsComm(nSpaces) )
3212  :
3213  M_worldsComm( _worldsComm ),
3214  M_worldComm( new WorldComm( _worldsComm[0] ) )
3215  {
3216  this->init( mesh, mesh_components, periodicity );
3217  }
3218 
3219  FunctionSpace( mesh_ptrtype const& mesh,
3220  std::vector<Dof > const& dofindices,
3221  periodicity_type periodicity = periodicity_type(),
3222  std::vector<WorldComm> const& _worldsComm = Environment::worldsComm(nSpaces) )
3223  :
3224  M_worldsComm( _worldsComm ),
3225  M_worldComm( new WorldComm( _worldsComm[0] ) )
3226  {
3227  this->init( mesh, 0, dofindices, periodicity );
3228  }
3229 
3234 #if 0 // ambiguous call with new below
3235  static pointer_type New( mesh_ptrtype const& __m, size_type mesh_components = MESH_RENUMBER | MESH_CHECK )
3236  {
3237  return pointer_type( new functionspace_type( __m, mesh_components ) );
3238  }
3239 #endif // 0
3240 
3241  static pointer_type New( mesh_ptrtype const& __m, std::vector<Dof > const& dofindices )
3242  {
3243  return pointer_type( new functionspace_type( __m, dofindices ) );
3244  }
3245 #if !defined( FEELPP_ENABLE_MPI_MODE)
3246  BOOST_PARAMETER_MEMBER_FUNCTION( ( pointer_type ),
3247  static New,
3248  tag,
3249  ( required
3250  ( mesh,* )
3251  )
3252  ( optional
3253  ( components, ( size_type ), MESH_RENUMBER | MESH_CHECK )
3254  ( periodicity,*,periodicity_type() )
3255  )
3256  )
3257  {
3258  return NewImpl( mesh, components, periodicity );
3259  }
3260  static pointer_type NewImpl( mesh_ptrtype const& __m,
3261  size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3262  periodicity_type periodicity = periodicity_type() )
3263  {
3264  return pointer_type( new functionspace_type( __m, mesh_components, periodicity ) );
3265  }
3266 #else
3267  BOOST_PARAMETER_MEMBER_FUNCTION( ( pointer_type ),
3268  static New,
3269  tag,
3270  ( required
3271  ( mesh,* )
3272  )
3273  ( optional
3274  ( worldscomm, *, Environment::worldsComm(nSpaces) )
3275  ( components, ( size_type ), MESH_RENUMBER | MESH_CHECK )
3276  ( periodicity,*,periodicity_type() )
3277  )
3278  )
3279  {
3280  return NewImpl( mesh, worldscomm, components, periodicity );
3281  }
3282 
3283  static pointer_type NewImpl( mesh_ptrtype const& __m,
3284  std::vector<WorldComm> const& worldsComm = Environment::worldsComm(nSpaces),
3285  size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3286  periodicity_type periodicity = periodicity_type() )
3287  {
3288 
3289  return pointer_type( new functionspace_type( __m, mesh_components, periodicity, worldsComm ) );
3290  }
3291 
3292 #endif
3293 
3296  void init( mesh_ptrtype const& mesh,
3297  size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3298  periodicity_type periodicity = periodicity_type() )
3299  {
3300  Context ctx( mesh_components );
3301  DVLOG(2) << "component MESH_RENUMBER: " << ctx.test( MESH_RENUMBER ) << "\n";
3302  DVLOG(2) << "component MESH_UPDATE_EDGES: " << ctx.test( MESH_UPDATE_EDGES ) << "\n";
3303  DVLOG(2) << "component MESH_UPDATE_FACES: " << ctx.test( MESH_UPDATE_FACES ) << "\n";
3304  DVLOG(2) << "component MESH_PARTITION: " << ctx.test( MESH_PARTITION ) << "\n";
3305 
3306  this->init( mesh, mesh_components, periodicity, std::vector<Dof >(), mpl::bool_<is_composite>() );
3307  //mesh->addObserver( *this );
3308  }
3309 
3310  void init( mesh_ptrtype const& mesh,
3311  size_type mesh_components,
3312  std::vector<Dof > const& dofindices,
3313  periodicity_type periodicity = periodicity_type() )
3314  {
3315 
3316  Context ctx( mesh_components );
3317  DVLOG(2) << "component MESH_RENUMBER: " << ctx.test( MESH_RENUMBER ) << "\n";
3318  DVLOG(2) << "component MESH_UPDATE_EDGES: " << ctx.test( MESH_UPDATE_EDGES ) << "\n";
3319  DVLOG(2) << "component MESH_UPDATE_FACES: " << ctx.test( MESH_UPDATE_FACES ) << "\n";
3320  DVLOG(2) << "component MESH_PARTITION: " << ctx.test( MESH_PARTITION ) << "\n";
3321 
3322  this->init( mesh, mesh_components, periodicity, dofindices, mpl::bool_<is_composite>() );
3323  //mesh->addObserver( *this );
3324  }
3325 
3328  {
3329  M_dof.reset();
3330  M_dofOnOff.reset();
3331  M_ref_fe.reset();
3332  }
3333 
3334 
3335  void setWorldsComm( std::vector<WorldComm> const& _worldsComm )
3336  {
3337  M_worldsComm=_worldsComm;
3338  };
3339  void setWorldComm( WorldComm const& _worldComm )
3340  {
3341  M_worldComm.reset( new WorldComm( _worldComm ) );
3342  };
3343 
3344  std::vector<WorldComm> const& worldsComm() const
3345  {
3346  return M_worldsComm;
3347  };
3348  WorldComm const& worldComm() const
3349  {
3350  return *M_worldComm;
3351  };
3352 
3353  bool hasEntriesForAllSpaces()
3354  {
3355  return (this->template mesh<0>()->worldComm().localSize() == this->template mesh<0>()->worldComm().globalSize() );
3356  }
3358 
3362 
3366  void operator()( MESH_CHANGES mesh_changes )
3367  {
3368  DVLOG(2) << "Update function space after a change in the mesh\n";
3369 
3370  }
3371 
3372 
3374 
3378 
3379 
3380  /*
3381  * Get the real point matrix in a context
3382  * 2 : Face
3383  * \todo see above
3384  */
3385  template<typename Context_t>
3387  ptsInContext( Context_t const & context, mpl::int_<2> ) const
3388  {
3389  //new context for the interpolation
3390  typedef typename Context_t::gm_type::template Context< Context_t::context|vm::POINT, typename Context_t::element_type> gmc_interp_type;
3391  typedef boost::shared_ptr<gmc_interp_type> gmc_interp_ptrtype;
3392 
3393  typedef typename Context_t::gm_type::template Context<Context_t::context,typename Context_t::element_type>::permutation_type permutation_type;
3394  typedef typename Context_t::gm_type::template Context<Context_t::context,typename Context_t::element_type>::precompute_ptrtype precompute_ptrtype;
3395 
3396  //not good because ?
3397  //gmc_interp_ptrtype __c_interp( new gmc_interp_type( context.geometricMapping(), context.element_c(),context.pcFaces(), context.faceId()) );
3398  //good with this
3399  std::vector<std::map<permutation_type, precompute_ptrtype> > __geo_pcfaces = context.pcFaces();
3400  gmc_interp_ptrtype __c_interp( new gmc_interp_type( context.geometricMapping(), context.element_c(), __geo_pcfaces , context.faceId() ) );
3401 
3402  return __c_interp->xReal();
3403  }
3404 
3405 
3409  uint16_type qDim() const
3410  {
3411  return N_COMPONENTS;
3412  }
3413 
3417  size_type nDof() const
3418  {
3419  return this->nDof( mpl::bool_<is_composite>() );
3420  }
3421 
3426  {
3427  return this->nLocalDof( mpl::bool_<is_composite>() );
3428  }
3429 
3430  size_type nLocalDofWithGhost() const
3431  {
3432  return this->nLocalDofWithGhost( mpl::bool_<is_composite>() );
3433  }
3434 
3435  size_type nLocalDofWithoutGhost() const
3436  {
3437  return this->nLocalDofWithoutGhost( mpl::bool_<is_composite>() );
3438  }
3439 
3440  size_type nLocalDofWithGhostOnProc( const int proc ) const
3441  {
3442  return this->nLocalDofWithGhostOnProc( proc, mpl::bool_<is_composite>() );
3443  }
3444 
3445  size_type nLocalDofWithoutGhostOnProc(const int proc) const
3446  {
3447  return this->nLocalDofWithoutGhostOnProc( proc, mpl::bool_<is_composite>() );
3448  }
3449 
3453  proc_dist_map_type getProcDistMap() const
3454  {
3455  return procDistMap;
3456  }
3457 
3461  size_type nDofStart( size_type i = /*invalid_size_type_value*/0 ) const
3462  {
3463  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NbDof( 0, i ) );
3464  return start;
3465  }
3466 
3467  size_type nLocalDofStart( size_type i = 0 ) const
3468  {
3469  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NLocalDof<mpl::bool_<true> >( this->worldsComm(),true,0,i ) );
3470  return start;
3471  }
3472 
3473  size_type nLocalDofWithGhostStart( size_type i = 0 ) const
3474  {
3475  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NLocalDof<mpl::bool_<true> >( this->worldsComm(),true,0,i ) );
3476  return start;
3477  }
3478 
3479  size_type nLocalDofWithoutGhostStart( size_type i = 0 ) const
3480  {
3481  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NLocalDof<mpl::bool_<false> >( this->worldsComm(),true,0,i ) );
3482  return start;
3483  }
3484 
3485  size_type nLocalDofWithGhostOnProcStart( const int proc, size_type i = 0 ) const
3486  {
3487  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<true> >( proc, this->worldsComm(),true,0,i ) );
3488  return start;
3489  }
3490 
3491  size_type nLocalDofWithoutGhostOnProcStart( const int proc, size_type i = 0 ) const
3492  {
3493  size_type start = fusion::accumulate( this->functionSpaces(), size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<false> >( proc, this->worldsComm(),true,0,i ) );
3494  return start;
3495  }
3496 
3497  uint16_type nSubFunctionSpace() const
3498  {
3499  return nSpaces;
3500  }
3501 
3506  {
3507  return this->nDof()/qDim();
3508  }
3509 
3510 
3514  mesh_ptrtype& mesh()
3515  {
3516  return M_mesh;
3517  }
3518 
3522  mesh_ptrtype const& mesh() const
3523  {
3524  return M_mesh;
3525  }
3526 
3527  template<int i>
3528  typename GetMesh<mesh_ptrtype,i>::type
3529  mesh() const
3530  {
3531  return mesh<i>( is_shared_ptr<mesh_ptrtype>() );
3532  }
3533  template<int i>
3534  typename GetMesh<mesh_ptrtype,i>::type
3535  mesh( mpl::bool_<true> ) const
3536  {
3537  return M_mesh;
3538  }
3539  template<int i>
3540  typename GetMesh<mesh_ptrtype,i>::type
3541  mesh( mpl::bool_<false> ) const
3542  {
3543  return fusion::at_c<i>(M_mesh);
3544  }
3545 
3549  basis_ptrtype const& basis() const
3550  {
3551  return M_ref_fe;
3552  }
3553 
3560  std::string basisName() const
3561  {
3562  return basisName( mpl::bool_<( nSpaces>1 )>() );
3563  }
3564  std::string basisName( mpl::bool_<true> ) const
3565  {
3566  return fusion::accumulate( this->functionSpaces(), std::string(), detail::BasisName() );
3567  }
3568  std::string basisName( mpl::bool_<false> ) const
3569  {
3570  return this->basis()->familyName();
3571  }
3572 
3579  std::vector<int> basisOrder() const
3580  {
3581  return basisOrder( mpl::bool_<( nSpaces>1 )>() );
3582  }
3583  std::vector<int> basisOrder( mpl::bool_<true> ) const
3584  {
3585  return fusion::accumulate( this->functionSpaces(), std::vector<int>(), detail::BasisOrder() );
3586  }
3587  std::vector<int> basisOrder( mpl::bool_<false> ) const
3588  {
3589  std::vector<int> o( 1 );
3590  o[0]=basis_type::nOrder;
3591  return o;
3592  }
3593 
3597  reference_element_ptrtype const& fe() const
3598  {
3599  return M_ref_fe;
3600  }
3601 
3605  gm_ptrtype const& gm() const
3606  {
3607  return M_mesh->gm();
3608  }
3609 
3613  gm1_ptrtype const& gm1() const
3614  {
3615  return M_mesh->gm1();
3616  }
3617 
3621  DataMap const& map() const
3622  {
3623  return *M_dof;
3624  }
3625 
3629  DataMap const& mapOn() const
3630  {
3631  return *M_dof;
3632  }
3633 
3637  DataMap const& mapOnOff() const
3638  {
3639  return *M_dofOnOff;
3640  }
3641 
3645  dof_ptrtype dof()
3646  {
3647  return M_dof;
3648  }
3649 
3653  dof_ptrtype dofOn()
3654  {
3655  return M_dof;
3656  }
3657 
3661  dof_ptrtype dofOnOff()
3662  {
3663  return M_dofOnOff;
3664  }
3665 
3669  dof_ptrtype const& dof() const
3670  {
3671  return M_dof;
3672  }
3673 
3677  dof_ptrtype const& dofOn() const
3678  {
3679  return M_dof;
3680  }
3681 
3685  dof_ptrtype const& dofOnOff() const
3686  {
3687  return M_dofOnOff;
3688  }
3689 
3693  //functionspace_vector_type const&
3694  //functionSpaces() const { return M_functionspaces; }
3695 
3696  functionspace_vector_type const&
3698  {
3699  return M_functionspaces;
3700  }
3701 
3702  std::vector<std::vector<size_type> > dofIndexSplit()
3703  {
3704  if ( nSpaces > 1 )
3705  {
3706  uint16_type cptSpaces=0;
3707  size_type start=0;
3708  std::vector<std::vector<size_type> > is(nSpaces);
3709  auto initial = boost::make_tuple( cptSpaces,start,is );
3710 
3711  // get start for each proc ->( proc0 : 0 ), (proc1 : sumdofproc0 ), (proc2 : sumdofproc0+sumdofproc1 ) ....
3712  auto startSplit = boost::fusion::fold( functionSpaces(), boost::make_tuple(0,0), detail::computeStartOfFieldSplit() ).template get<1>();
3713  // compute split
3714  auto result = boost::fusion::fold( functionSpaces(), initial, detail::computeNDofForEachSpace(startSplit) );
3715  is = result.template get<2>();
3716 
3717 
3718 #if 0
3719  for ( int proc = 0; proc<this->worldComm().globalSize(); ++proc )
3720  {
3721  this->worldComm().globalComm().barrier();
3722  if ( proc==this->worldComm().globalRank() )
3723  {
3724  std::cout << "proc " << proc << "\n";
3725  std::cout << "split size=" << result.template get<2>().size() << " nspace=" << nSpaces << "\n";
3726  std::cout << "split:\n";
3727 
3728  //std::cout << "\n\n";
3729 
3730 
3731  for ( int s= 0; s < nSpaces; ++s )
3732  {
3733  std::cout << "space: " << is[s].size() << "\n";
3734 
3735  for ( int i = 0; i < is[s].size(); ++i )
3736  {
3737  std::cout << is[s][i] << " ";
3738  }
3739 
3740  std::cout << "\n\n";
3741  }
3742  }
3743  }
3744 
3745 #endif
3746  return is;
3747  }
3748 
3749  std::vector<std::vector<size_type> > is;
3750  is.push_back( std::vector<size_type>( nLocalDof() ) );
3751  int index = 0;
3752  //for( int& i : is[0] ) { i = index++; }
3753  BOOST_FOREACH( auto& i, is[0] )
3754  {
3755  i = index++;
3756  }
3757  return is;
3758 
3759  }
3763  element_type
3764  element( std::string const& name = "u" )
3765  {
3766  element_type u( this->shared_from_this(), name );
3767  u.zero();
3768  return u;
3769  }
3770 
3774  element_ptrtype
3775  elementPtr( std::string const& name = "u" )
3776  {
3777  element_ptrtype u( new element_type( this->shared_from_this(), name ) );
3778  u->zero();
3779  return u;
3780  }
3781 
3785  template<int i>
3786  typename mpl::at_c<functionspace_vector_type,i>::type
3788  {
3789  return fusion::at_c<i>( M_functionspaces );
3790  }
3791 
3796  component_functionspace_ptrtype const& compSpace() const
3797  {
3798  return M_comp_space;
3799  }
3800 
3804  trace_functionspace_ptrtype
3805  trace() const
3806  {
3807  //return trace( mpl::greater<mpl::int_<nDim>,mpl::int_<1> >::type() )
3808  return trace_functionspace_type::New( mesh()->trace( boundaryfaces( mesh() ) ) );
3809  }
3810  template<typename RangeT>
3811  trace_functionspace_ptrtype
3812  trace( RangeT range ) const
3813  {
3814  return trace_functionspace_type::New( mesh()->trace( range ) );
3815  }
3816 
3817 
3821  trace_trace_functionspace_ptrtype
3822  wireBasket() const
3823  {
3824  //return trace( mpl::greater<mpl::int_<nDim>,mpl::int_<1> >::type() )
3825  return trace_trace_functionspace_type::New( mesh()->wireBasket( markededges( mesh(),"WireBasket" ) ) );
3826  }
3827  template<typename RangeT>
3828  trace_trace_functionspace_ptrtype
3829  wireBasket( RangeT range ) const
3830  {
3831  return trace_trace_functionspace_type::New( mesh()->wireBasket( range ) );
3832  }
3833 
3834 
3838  bool hasRegionTree() const
3839  {
3840  return M_rt;
3841  }
3842 
3846  region_tree_ptrtype const& regionTree() const;
3847 
3851  void updateRegionTree() const;
3852 
3863  bool findPoint( node_type const& pt, size_type &cv, node_type& ptr ) const;
3864 
3868  void rebuildDofPoints() { rebuildDofPoints(mpl::bool_<is_composite>()); }
3869 
3871 
3875 
3876 
3878 
3882 
3883  void printInfo() const
3884  {
3885  LOG(INFO) << " number of components : " << qDim() << "\n";
3886  LOG(INFO) << " n Global Dof : " << nDof() << "\n";
3887  LOG(INFO) << " n Local Dof : " << nLocalDof() << "\n";
3888  }
3889 
3891 
3892 
3893  FunctionSpace( FunctionSpace const& __fe )
3894  :
3895  M_worldsComm( __fe.M_worldsComm ),
3896  M_worldComm( __fe.M_worldComm ),
3897  M_mesh( __fe.M_mesh ),
3898  M_ref_fe( __fe.M_ref_fe ),
3899  M_comp_space( __fe.M_comp_space ),
3900  M_dof( __fe.M_dof ),
3901  M_dofOnOff( __fe.M_dofOnOff ),
3902  M_rt( __fe.M_rt )
3903  {
3904  DVLOG(2) << "copying FunctionSpace\n";
3905  }
3906 
3907 protected:
3908 
3909 private:
3910 
3911  void init( mesh_ptrtype const& mesh,
3912  size_type mesh_components,
3913  periodicity_type const& periodicity,
3914  std::vector<Dof > const& dofindices,
3915  mpl::bool_<false> );
3916  void init( mesh_ptrtype const& mesh,
3917  size_type mesh_components,
3918  periodicity_type const& periodicity,
3919  std::vector<Dof > const& dofindices,
3920  mpl::bool_<true> );
3921 
3922  size_type nDof( mpl::bool_<false> ) const;
3923  size_type nDof( mpl::bool_<true> ) const;
3924 
3925  size_type nLocalDof( mpl::bool_<false> ) const;
3926  size_type nLocalDof( mpl::bool_<true> ) const;
3927 
3928  size_type nLocalDofWithGhost( mpl::bool_<false> ) const;
3929  size_type nLocalDofWithGhost( mpl::bool_<true> ) const;
3930  size_type nLocalDofWithoutGhost( mpl::bool_<false> ) const;
3931  size_type nLocalDofWithoutGhost( mpl::bool_<true> ) const;
3932 
3933  size_type nLocalDofWithGhostOnProc( const int proc, mpl::bool_<false> ) const;
3934  size_type nLocalDofWithGhostOnProc( const int proc, mpl::bool_<true> ) const;
3935  size_type nLocalDofWithoutGhostOnProc( const int proc, mpl::bool_<false> ) const;
3936  size_type nLocalDofWithoutGhostOnProc( const int proc, mpl::bool_<true> ) const;
3937 
3938  void rebuildDofPoints( mpl::bool_<false> );
3939  void rebuildDofPoints( mpl::bool_<true> );
3940 
3941  friend class ComponentSpace;
3942  class ComponentSpace
3943  {
3944  public:
3945  typedef FunctionSpace<A0,A1,A2,A3,A4> functionspace_type;
3946  typedef functionspace_type* functionspace_ptrtype;
3947  typedef functionspace_type const* functionspace_cptrtype;
3948 
3949  typedef typename FunctionSpace<A0,A1,A2,A3,A4>::component_functionspace_type component_functionspace_type;
3950  typedef typename FunctionSpace<A0,A1,A2,A3,A4>::component_functionspace_ptrtype component_functionspace_ptrtype;
3951  typedef component_functionspace_type const* component_functionspace_cptrtype;
3952 
3953  ComponentSpace( FunctionSpace<A0,A1,A2,A3,A4> * __functionspace,
3954  mesh_ptrtype __m )
3955  :
3956  M_functionspace( __functionspace ),
3957  M_mesh( __m )
3958  {}
3959  component_functionspace_ptrtype operator()()
3960  {
3961  return operator()( mpl::bool_<functionspace_type::is_scalar>() );
3962  }
3963  component_functionspace_ptrtype operator()( mpl::bool_<true> )
3964  {
3965  FEELPP_ASSERT( 0 ).error( "invalid call for component space extraction" );
3966  //return M_functionspace;
3967  return component_functionspace_ptrtype();
3968  }
3969  component_functionspace_ptrtype operator()( mpl::bool_<false> )
3970  {
3971  return component_functionspace_type::New( M_mesh );
3972  }
3973 
3974  private:
3975 
3976  FunctionSpace<A0,A1,A2,A3,A4> * M_functionspace;
3977  mesh_ptrtype M_mesh;
3978  };
3979 
3980 #if 0
3981  template<typename ElementRange, typename OnExpr>
3982  void
3983  on( ElementRange const& _range,
3984  OnExpr const& _expr )
3985  {
3986  M_constraints.push_back( vf::project( this->shared_from_this(), _range, _expr ) );
3987  }
3988 #endif // 0
3989 
3990 protected:
3991 
3992  //friend class FunctionSpace<mesh_type, typename bases_list::component_basis_type, value_type>;
3993  //friend class FunctionSpace<mesh_type, bases_list, value_type>;
3994 
3995  std::vector<WorldComm> M_worldsComm;
3996  boost::shared_ptr<WorldComm> M_worldComm;
3997 
3998  // finite element mesh
3999  mesh_ptrtype M_mesh;
4000 
4002  reference_element_ptrtype M_ref_fe;
4003 
4005  component_functionspace_ptrtype M_comp_space;
4006 
4008  dof_ptrtype M_dof;
4009 
4011  dof_ptrtype M_dofOnOff;
4012 
4014  mutable boost::optional<region_tree_ptrtype> M_rt;
4015 
4016  //mutable boost::prof::basic_profiler<boost::prof::basic_profile_manager<std::string, real_type, boost::high_resolution_timer, boost::prof::empty_logging_policy, boost::prof::default_stats_policy<std::string, real_type> > > M_prof_find_points;
4017 private:
4018 
4020  FunctionSpace();
4021 
4022  functionspace_vector_type M_functionspaces;
4023 
4024  proc_dist_map_type procDistMap;
4025 #if 0
4026  std::list<Element> M_constraints;
4027 #endif
4028 }; // FunctionSpace
4029 
4030 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4031 template<typename T, typename Cont>
4032 const bool FunctionSpace<A0,A1,A2,A3, A4>::Element<T,Cont>::is_scalar;
4033 
4034 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4035 template<typename T, typename Cont>
4036 const bool FunctionSpace<A0,A1,A2,A3,A4>::Element<T,Cont>::is_vectorial;
4037 
4038 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4039 template<typename T, typename Cont>
4040 const bool FunctionSpace<A0,A1,A2,A3,A4>::Element<T,Cont>::is_tensor2;
4041 
4042 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4043 template<typename T, typename Cont>
4044 const uint16_type FunctionSpace<A0,A1,A2,A3,A4>::Element<T,Cont>::nComponents;
4045 
4046 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4047 void
4048 FunctionSpace<A0, A1, A2, A3, A4>::init( mesh_ptrtype const& __m,
4049  size_type mesh_components,
4050  periodicity_type const& periodicity,
4051  std::vector<Dof> const& dofindices,
4052  mpl::bool_<false> )
4053 {
4054  DVLOG(2) << "calling init(<space>) begin\n";
4055  DVLOG(2) << "calling init(<space>) is_periodic: " << is_periodic << "\n";
4056  M_mesh = __m;
4057 
4058 
4059  if ( basis_type::nDofPerEdge || nDim >= 3 )
4060  mesh_components |= MESH_UPDATE_EDGES;
4061 
4062  /*
4063  * update faces info in mesh only if dofs exists on faces or the
4064  * expansion is continuous between elements. This case handles strong
4065  * Dirichlet imposition
4066  */
4067  if ( basis_type::nDofPerFace || is_continuous || nDim >= 3 )
4068  mesh_components |= MESH_UPDATE_FACES;
4069 
4070  M_mesh->components().set( mesh_components );
4071 
4072  M_mesh->updateForUse();
4073 
4074  if ( is_periodic )
4075  {
4076  M_mesh->removeFacesFromBoundary( { periodicity.tag1(), periodicity.tag2() } );
4077  }
4078 
4079  M_ref_fe = basis_ptrtype( new basis_type );
4080 
4081  M_dof = dof_ptrtype( new dof_type( M_ref_fe, fusion::at_c<0>(periodicity), this->worldsComm()[0] ) );
4082 
4083  DVLOG(2) << "[functionspace] Dof indices is empty ? " << dofindices.empty() << "\n";
4084  M_dof->setDofIndices( dofindices );
4085  DVLOG(2) << "[functionspace] is_periodic = " << is_periodic << "\n";
4086 
4087  M_dof->build( M_mesh );
4088 
4089  M_dofOnOff = M_dof;
4090 
4091  if ( is_vectorial )
4092  {
4093  // Warning: this works regarding the communicator . for the component space
4094  // it will use in mixed spaces only numberofSudomains/numberofspace processors
4095  //
4096  M_comp_space = component_functionspace_ptrtype( new component_functionspace_type( M_mesh,
4097  MESH_COMPONENTS_DEFAULTS,
4098  periodicity,
4099  std::vector<WorldComm>( 1,this->worldsComm()[0] ) ) );
4100  }
4101 
4102  DVLOG(2) << "nb dim : " << qDim() << "\n";
4103  DVLOG(2) << "nb dof : " << nDof() << "\n";
4104  DVLOG(2) << "nb dof per component: " << nDofPerComponent() << "\n";
4105 
4106  if ( is_vectorial )
4107  {
4108  DVLOG(2) << "component space :: nb dim : " << M_comp_space->qDim() << "\n";
4109  DVLOG(2) << "component space :: nb dof : " << M_comp_space->nDof() << "\n";
4110  DVLOG(2) << "component space :: nb dof per component: " << M_comp_space->nDofPerComponent() << "\n";
4111  }
4112 
4113  //detail::searchIndicesBySpace<proc_dist_map_type>( this, procDistMap);
4114 
4115  DVLOG(2) << "calling init(<space>) end\n";
4116 
4117 }
4118 
4119 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4120 void
4121 FunctionSpace<A0, A1, A2, A3, A4>::init( mesh_ptrtype const& __m,
4122  size_type mesh_components,
4123  periodicity_type const& periodicity,
4124  std::vector<Dof> const& dofindices,
4125  mpl::bool_<true> )
4126 {
4127  DVLOG(2) << "calling init(<composite>) begin\n";
4128  M_mesh = __m;
4129 
4130  // todo : check worldsComm size and M_functionspaces are the same!
4131  fusion::for_each( M_functionspaces, detail::InitializeSpace<mesh_ptrtype,periodicity_type>( __m, periodicity,
4132  dofindices, this->worldsComm() ) );
4133 
4134 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
4135  M_dof = dof_ptrtype( new dof_type( this->nDof(), this->nLocalDof() ) );
4136  DVLOG(2) << "calling nDof(<composite>)" << this->nDof() << "\n";
4137  DVLOG(2) << "calling init(<composite>) end\n";
4138 
4139  proc_dist_map_type emptyMap;
4140  procDistMap = fusion::accumulate( M_functionspaces,
4141  emptyMap,
4142  detail::searchIndicesBySpace<proc_dist_map_type>() );
4143  M_dofOnOff = M_dof;
4144 #else // new version with MPI
4145 
4146  if ( this->worldComm().globalSize()>1 )
4147  {
4148  if ( this->hasEntriesForAllSpaces() )
4149  {
4150  // construction with same partionment for all subspaces
4151  // and each processors has entries for all subspaces
4152  DVLOG(2) << "init(<composite>) type hasEntriesForAllSpaces\n";
4153  // build usefull data for detail::updateDataMapProcessStandard
4154  const int worldsize = this->worldComm().globalSize();
4155  std::vector<size_type> startDofGlobalCluster(worldsize);
4156  std::vector<size_type> nLocalDofWithoutGhostWorld(worldsize),nLocalDofWithGhostWorld(worldsize);
4157  for (int proc = 0; proc<worldsize ; ++proc )
4158  {
4159  nLocalDofWithoutGhostWorld[proc] = this->nLocalDofWithoutGhostOnProc(proc);
4160  nLocalDofWithGhostWorld[proc] = this->nLocalDofWithGhostOnProc(proc);
4161  }
4162 
4163  startDofGlobalCluster[0]=0;
4164  for (int p=1;p<worldsize;++p)
4165  {
4166  startDofGlobalCluster[p] = startDofGlobalCluster[p-1] + nLocalDofWithoutGhostWorld[p-1];
4167  }
4168 
4169  // build datamap
4170  auto dofInitTool=detail::updateDataMapProcessStandard<dof_type>( this->worldsComm(), this->worldComm(),
4171  this->nSubFunctionSpace()-1,
4172  startDofGlobalCluster,
4173  nLocalDofWithoutGhostWorld,
4174  nLocalDofWithGhostWorld
4175  );
4176  fusion::for_each( M_functionspaces, dofInitTool );
4177  // finish update datamap
4178  M_dof = dofInitTool.dataMap();
4179  M_dof->setNDof( this->nDof() );
4180  M_dofOnOff = M_dof;
4181  }
4182  else
4183  {
4184  // construction with same partionment for all subspaces
4185  // and one processor has entries for only one subspace
4186  DVLOG(2) << "init(<composite>) type Not hasEntriesForAllSpaces\n";
4187 
4188  // build the WorldComm associated to mix space
4189  WorldComm mixSpaceWorldComm = this->worldsComm()[0];
4190 
4191  if ( this->worldsComm().size()>1 )
4192  for ( int i=1; i<( int )this->worldsComm().size(); ++i )
4193  {
4194  mixSpaceWorldComm = mixSpaceWorldComm + this->worldsComm()[i];
4195  }
4196 
4197  this->setWorldComm( mixSpaceWorldComm );
4198  //mixSpaceWorldComm.showMe();
4199 
4200  // update DofTable for the mixedSpace (we have 2 dofTables : On and OnOff)
4201  auto dofInitTool=detail::updateDataMapProcess<dof_type>( this->worldsComm(), mixSpaceWorldComm, this->nSubFunctionSpace()-1 );
4202  fusion::for_each( M_functionspaces, dofInitTool );
4203  // finish update datamap
4204  M_dof = dofInitTool.dataMap();
4205  M_dof->setNDof( this->nDof() );
4206  M_dof->updateDataInWorld();
4207  M_dofOnOff = dofInitTool.dataMapOnOff();
4208  M_dofOnOff->setNDof( this->nDof() );
4209  M_dofOnOff->updateDataInWorld();
4210  }
4211  }
4212 
4213  else // sequential
4214  {
4215  // update DofTable for the mixedSpace (here On is not build properly but OnOff yes and On=OnOff, see detail::updateDataMapProcess)
4216  auto dofInitTool=detail::updateDataMapProcess<dof_type>( this->worldsComm(), this->worldComm(), this->nSubFunctionSpace()-1 );
4217  fusion::for_each( M_functionspaces, dofInitTool );
4218  M_dof = dofInitTool.dataMapOnOff();
4219  M_dof->setNDof( this->nDof() );
4220  M_dofOnOff = dofInitTool.dataMapOnOff();
4221  M_dofOnOff->setNDof( this->nDof() );
4222  }
4223 
4224 #endif
4225 }
4226 
4227 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4228 size_type
4229 FunctionSpace<A0, A1, A2, A3, A4>::nDof( mpl::bool_<true> ) const
4230 {
4231  DVLOG(2) << "calling nDof(<composite>) begin\n";
4232  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NbDof() );
4233  DVLOG(2) << "calling nDof(<composite>) end\n";
4234  return ndof;
4235 }
4236 
4237 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4238 size_type
4239 FunctionSpace<A0, A1, A2, A3, A4>::nDof( mpl::bool_<false> ) const
4240 {
4241  return M_dof->nDof();
4242 }
4243 
4244 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4245 size_type
4246 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDof( mpl::bool_<true> ) const
4247 {
4248  DVLOG(2) << "calling nLocalDof(<composite>) begin\n";
4249  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NLocalDof<mpl::bool_<true> >( this->worldsComm() ) );
4250  DVLOG(2) << "calling nLocalDof(<composite>) end\n";
4251  return ndof;
4252 }
4253 
4254 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4255 size_type
4256 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDof( mpl::bool_<false> ) const
4257 {
4258  //return M_dof->nLocalDof();
4259  return M_dof->nLocalDofWithGhost();
4260 }
4261 
4262 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4263 size_type
4264 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhost( mpl::bool_<true> ) const
4265 {
4266  DVLOG(2) << "calling nLocalDof(<composite>) begin\n";
4267  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NLocalDof<mpl::bool_<true> >( this->worldsComm() ) );
4268  DVLOG(2) << "calling nLocalDof(<composite>) end\n";
4269  return ndof;
4270 }
4271 
4272 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4273 size_type
4274 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhost( mpl::bool_<false> ) const
4275 {
4276  return M_dof->nLocalDofWithGhost();
4277 }
4278 
4279 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4280 size_type
4281 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhostOnProc( const int proc, mpl::bool_<true> ) const
4282 {
4283  DVLOG(2) << "calling nLocalDof(<composite>) begin\n";
4284  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<true> >( proc, this->worldsComm() ) );
4285  DVLOG(2) << "calling nLocalDof(<composite>) end\n";
4286  return ndof;
4287 }
4288 
4289 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4290 size_type
4291 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhostOnProc( const int proc, mpl::bool_<false> ) const
4292 {
4293  return M_dof->nLocalDofWithGhost(proc);
4294 }
4295 
4296 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4297 size_type
4298 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhost( mpl::bool_<true> ) const
4299 {
4300  DVLOG(2) << "calling nLocalDof(<composite>) begin\n";
4301  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NLocalDof<mpl::bool_<false> >( this->worldsComm() ) );
4302  DVLOG(2) << "calling nLocalDof(<composite>) end\n";
4303  return ndof;
4304 }
4305 
4306 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4307 size_type
4308 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhost( mpl::bool_<false> ) const
4309 {
4310  return M_dof->nLocalDofWithoutGhost();
4311 }
4312 
4313 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4314 size_type
4315 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhostOnProc( const int proc, mpl::bool_<true> ) const
4316 {
4317  DVLOG(2) << "calling nLocalDof(<composite>) begin\n";
4318  size_type ndof = fusion::accumulate( M_functionspaces, size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<false> >( proc, this->worldsComm() ) );
4319  DVLOG(2) << "calling nLocalDof(<composite>) end\n";
4320  return ndof;
4321 }
4322 
4323 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4324 size_type
4325 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhostOnProc( const int proc, mpl::bool_<false> ) const
4326 {
4327  return M_dof->nLocalDofWithoutGhost(proc);
4328 }
4329 
4330 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4331 void
4333 {
4334  M_dof->rebuildDofPoints( *M_mesh );
4335 }
4336 
4337 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4338 void
4340 {
4341  fusion::for_each( M_functionspaces, detail::rebuildDofPointsTool() );
4342 }
4343 
4344 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4345 void
4347 {
4348  scalar_type EPS=1E-13;
4349 
4351 
4352  __rt->clear();
4353  BoundingBox<> __bb( M_mesh->gm()->isLinear() );
4354 
4355  typedef typename mesh_type::element_iterator mesh_element_iterator;
4356  mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().rank() );
4357  mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().rank() );
4358 
4359  for ( size_type __i = 0; it != en; ++__i, ++it )
4360  {
4361  __bb.make( it->G() );
4362 
4363  for ( unsigned k=0; k < __bb.min.size(); ++k )
4364  {
4365  __bb.min[k]-=EPS;
4366  __bb.max[k]+=EPS;
4367  }
4368 
4369  __rt->addBox( __bb.min, __bb.max, it->id() );
4370  }
4371 
4372  //__rt->dump();
4373  M_rt = __rt;
4374 }
4375 
4376 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4377 region_tree_ptrtype const&
4379 {
4380  if ( !M_rt )
4381  {
4382  scalar_type EPS=1E-13;
4383 
4385 
4386  __rt->clear();
4387  BoundingBox<> __bb( M_mesh->gm()->isLinear() );
4388 
4389  typedef typename mesh_type::element_iterator mesh_element_iterator;
4390  mesh_element_iterator it = M_mesh->beginElementWithProcessId( M_mesh->comm().rank() );
4391  mesh_element_iterator en = M_mesh->endElementWithProcessId( M_mesh->comm().rank() );
4392 
4393  for ( size_type __i = 0; it != en; ++__i, ++it )
4394  {
4395  __bb.make( it->G() );
4396 
4397  for ( unsigned k=0; k < __bb.min.size(); ++k )
4398  {
4399  __bb.min[k]-=EPS;
4400  __bb.max[k]+=EPS;
4401  }
4402 
4403  __rt->addBox( __bb.min, __bb.max, it->id() );
4404  }
4405 
4406  //__rt->dump();
4407  M_rt = __rt;
4408  }
4409 
4410  return M_rt.get();
4411 }
4412 
4413 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4414 bool
4415 FunctionSpace<A0, A1, A2, A3, A4>::findPoint( node_type const& pt,size_type &cv , node_type &ptr ) const
4416 {
4417  if ( !hasRegionTree() )
4418  regionTree();
4419 
4420  //M_prof_find_points.resume();
4421 
4422  region_tree_type* __rt = M_rt.get().get();
4423 
4424  region_tree_type::pbox_set_type boxlst;
4425  __rt->findBoxesAtPoint( pt, boxlst );
4426 
4427  typedef typename gm_type::Inverse inv_trans_type;
4428  typename gm_type::reference_convex_type refelem;
4429 
4430  std::pair<size_type,value_type> closest = std::make_pair( invalid_size_type_value, -1e30 );
4431 
4432  region_tree_type::pbox_set_type::const_iterator it = boxlst.begin();
4433  region_tree_type::pbox_set_type::const_iterator ite = boxlst.end();
4434 
4435  for ( ; it != ite; ++it )
4436  {
4437  inv_trans_type __git( M_mesh->gm(), M_mesh->element( ( *it )->id ), this->worldComm().subWorldCommSeq() );
4438 
4439  size_type cv_stored = ( *it )->id;
4440 
4441 
4442  DVLOG(2) << "[FunctionSpace::findPoint] id : " << cv_stored << "\n";
4443 
4444  __git.setXReal( pt );
4445  ptr = __git.xRef();
4446 
4447 
4448  bool isin;
4449  value_type dmin;
4450  boost::tie( isin, dmin ) = refelem.isIn( ptr );
4451  DVLOG(2) << "[FunctionSpace::findPoint] isin: " << isin << " dmin: " << dmin << "\n";
4452 
4453  closest = ( dmin > closest.second )?std::make_pair( cv_stored, dmin ):closest;
4454 
4455  if ( isin )
4456  {
4457  DVLOG(2) << "[FunctionSpace::findPoint] id of the convex where " << pt << " belongs : " << cv_stored << "\n";
4458  DVLOG(2) << "[FunctionSpace::findPoint] ref coordinate: " << ptr << "\n";
4459  cv = ( *it )->id;
4460  //M_prof_find_points.pause();
4461  return true;
4462  }
4463  }
4464 
4465  cv=closest.first;
4466  //M_prof_find_points.pause();
4467  return false;
4468 }
4469 
4470 //
4471 // Element implementation
4472 //
4473 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4474 template<typename Y, typename Cont>
4476  :
4477  super(),
4478  M_start( 0 ),
4479  M_ct( NO_COMPONENT ),
4480  M_containersOffProcess( boost::none )
4481 {}
4482 
4483 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4484 template<typename Y, typename Cont>
4485 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::Element( Element const& __e )
4486  :
4487  super( __e ),
4488  M_functionspace( __e.M_functionspace ),
4489  M_name( __e.M_name ),
4490  M_start( __e.M_start ),
4491  M_ct( __e.M_ct ),
4492  M_containersOffProcess( __e.M_containersOffProcess )
4493 {
4494  DVLOG(2) << "Element<copy>::range::start = " << this->start() << "\n";
4495  DVLOG(2) << "Element<copy>::range::size = " << this->size() << "\n";
4496 
4497 }
4498 
4499 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4500 template<typename Y, typename Cont>
4501 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::Element( functionspace_ptrtype const& __functionspace,
4502  std::string const& __name,
4503  size_type __start,
4504  ComponentType __ct )
4505  :
4506  super( __functionspace->dof() ),
4507  M_functionspace( __functionspace ),
4508  M_name( __name ),
4509  M_start( __start ),
4510  M_ct( __ct ),
4511  M_containersOffProcess( boost::none )
4512 {
4513  DVLOG(2) << "Element::start = " << this->start() << "\n";
4514  DVLOG(2) << "Element::size = " << this->size() << "\n";
4515  DVLOG(2) << "Element::ndof = " << this->nDof() << "\n";
4516  DVLOG(2) << "Element::nlocaldof = " << this->nLocalDof() << "\n";
4517 }
4518 
4519 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4520 template<typename Y, typename Cont>
4521 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::Element( functionspace_ptrtype const& __functionspace,
4522  container_type const& __c,
4523  std::string const& __name,
4524  size_type __start,
4525  ComponentType __ct )
4526  :
4527  super( __c ),
4528  M_functionspace( __functionspace ),
4529  M_name( __name ),
4530  M_start( __start ),
4531  M_ct( __ct ),
4532  M_containersOffProcess( boost::none )
4533 {
4534  DVLOG(2) << "Element<range>::range::start = " << __c.start() << "\n";
4535  DVLOG(2) << "Element<range>::range::size = " << __c.size() << "\n";
4536  DVLOG(2) << "Element<range>::start = " << this->start() << "\n";
4537  DVLOG(2) << "Element<range>::size = " << this->size() << "\n";
4538  DVLOG(2) << "Element<range>::ndof = " << this->nDof() << "\n";
4539  DVLOG(2) << "Element<range>::nlocaldof = " << this->nLocalDof() << "\n";
4540  M_start = __c.start();
4541 }
4542 
4543 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4544 template<typename Y, typename Cont>
4545 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::~Element()
4546 {}
4547 
4548 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4549 template<typename Y, typename Cont>
4550 void
4551 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::initFromSpace( functionspace_ptrtype const& __functionspace,
4552  container_type const& __c )
4553 {
4554  M_functionspace = __functionspace;
4555  ( container_type )*this = __c;
4556 }
4557 
4558 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4559 template<typename Y, typename Cont>
4560 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>&
4562 {
4563  if ( this != &__e )
4564  {
4565  M_functionspace = __e.M_functionspace;
4566 
4567  if ( __e.M_name != "unknown" )
4568  M_name = __e.M_name;
4569 
4570  M_start = __e.M_start;
4571  M_ct = __e.M_ct;
4572  M_containersOffProcess = __e.M_containersOffProcess;
4573  this->resize( M_functionspace->nLocalDof() );
4574  super::operator=( __e );
4575  this->outdateGlobalValues();
4576  }
4577 
4578  return *this;
4579 }
4580 
4581 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4582 template<typename Y, typename Cont>
4583 template<typename VectorExpr>
4584 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>&
4586 {
4587  if ( __v.size() != this->size() )
4588  {
4589  std::ostringstream __err;
4590  __err << "Invalid vector size this->size()=" << this->size()
4591  << " and v.size()=" << __v.size();
4592  throw std::logic_error( __err.str() );
4593  }
4594 
4595  this->outdateGlobalValues();
4596  super::operator=( __v );
4597  return *this;
4598 }
4599 
4600 
4601 //
4602 // Interpolation tools
4603 //
4604 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4605 template<typename Y, typename Cont>
4606 typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::id_type
4607 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::operator()( node_type const& __x, bool extrapolate ) const
4608 {
4609  this->updateGlobalValues();
4610 
4611  node_type __x_ref;
4612  size_type __cv_id;
4613  int rank = functionSpace()->mesh()->comm().rank();
4614  int nprocs = functionSpace()->mesh()->comm().size();
4615  std::vector<int> found_pt( nprocs, 0 );
4616  std::vector<int> global_found_pt( nprocs, 0 );
4617 
4618  if ( functionSpace()->findPoint( __x, __cv_id, __x_ref ) || extrapolate )
4619  {
4620 #if !defined( NDEBUG )
4621  DVLOG(2) << "Point " << __x << " is in element " << __cv_id << " pt_ref=" << __x_ref << "\n";
4622 #endif
4623  gm_ptrtype __gm = functionSpace()->gm();
4624  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
4625  typedef typename gm_type::precompute_type geopc_type;
4626  typename matrix_node<value_type>::type pts( __x_ref.size(), 1 );
4627  ublas::column( pts, 0 ) = __x_ref;
4628  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
4629 
4630 
4631  typedef typename gm_type::template Context<vm::POINT|vm::GRAD|vm::KB|vm::JACOBIAN, geoelement_type> gmc_type;
4632  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
4633  gmc_ptrtype __c( new gmc_type( __gm,
4634  functionSpace()->mesh()->element( __cv_id ),
4635  __geopc ) );
4636  pc_ptrtype pc( new pc_type( this->functionSpace()->fe(), pts ) );
4637 
4638  typedef typename mesh_type::element_type geoelement_type;
4639  typedef typename functionspace_type::fe_type fe_type;
4640  typedef typename fe_type::template Context<vm::POINT|vm::GRAD|vm::KB|vm::JACOBIAN, fe_type, gm_type, geoelement_type> fectx_type;
4641  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
4642  fectx_ptrtype fectx( new fectx_type( this->functionSpace()->fe(),
4643  __c,
4644  pc ) );
4645  found_pt[ rank ] = 1;
4646 
4647 #if defined(FEELPP_HAS_MPI)
4648 
4649  if ( nprocs > 1 )
4650  {
4651  //mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, std::plus<std::vector<int> >() );
4652  mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4653  }
4654 
4655 #else
4656  global_found_pt[ 0 ] = found_pt[ 0 ];
4657 #endif /* FEELPP_HAS_MPI */
4658 
4659  id_type __id( this->id( *fectx ) );
4660 
4661  //DVLOG(2) << "[interpolation] id = " << __id << "\n";
4662 #if defined(FEELPP_HAS_MPI)
4663  DVLOG(2) << "sending interpolation context to all processors from " << functionSpace()->mesh()->comm().rank() << "\n";
4664 
4665  if ( functionSpace()->mesh()->comm().size() > 1 )
4666  {
4667  mpi::broadcast( functionSpace()->mesh()->comm(), __id, functionSpace()->mesh()->comm().rank() );
4668  }
4669 
4670  //DVLOG(2) << "[interpolation] after broadcast id = " << __id << "\n";
4671 #endif /* FEELPP_HAS_MPI */
4672  return __id;
4673  }
4674 
4675  else
4676  {
4677 #if defined(FEELPP_HAS_MPI)
4678 
4679  if ( functionSpace()->mesh()->comm().size() > 1 )
4680  {
4681  //mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, std::plus<std::vector<int> >() );
4682  mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4683  }
4684 
4685 #endif /* FEELPP_HAS_MPI */
4686  bool found = false;
4687  size_type i = 0;
4688 
4689  for ( ; i < global_found_pt.size(); ++i )
4690  if ( global_found_pt[i] != 0 )
4691  {
4692  DVLOG(2) << "processor " << i << " has the point " << __x << "\n";
4693  found = true;
4694  break;
4695  }
4696 
4697  id_type __id;
4698 
4699  if ( found )
4700  {
4701  DVLOG(2) << "receiving interpolation context from processor " << i << "\n";
4702 #if defined(FEELPP_HAS_MPI)
4703 
4704  if ( functionSpace()->mesh()->comm().size() > 1 )
4705  mpi::broadcast( functionSpace()->mesh()->comm(), __id, i );
4706 
4707 #endif /* FEELPP_HAS_MPI */
4708  }
4709 
4710  else
4711  {
4712  LOG(WARNING) << "no processor seems to have the point " << __x << "\n";
4713  }
4714 
4715  return __id;
4716  }
4717 
4718 } // operator()
4719 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4720 template<typename Y, typename Cont>
4721 template<typename Context_t>
4722 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
4723 void
4724 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::id_( Context_t const & context, id_array_type& v ) const
4725 {
4726  if ( !this->areGlobalValuesUpdated() )
4727  this->updateGlobalValues();
4728 
4729  size_type elt_id = context.eId();
4730  if ( context.gmContext()->element().mesh()->isSubMeshFrom( this->mesh() ) )
4731  elt_id = context.gmContext()->element().mesh()->subMeshToMesh( context.eId() );
4732  if ( context.gmContext()->element().mesh()->isParentMeshOf( this->mesh() ) )
4733  elt_id = this->mesh()->meshToSubMesh( context.eId() );
4734  if ( elt_id == invalid_size_type_value )
4735  return;
4736 
4737  const uint16_type nq = context.xRefs().size2();
4738 
4739  //array_type v( boost::extents[nComponents1][nComponents2][context.xRefs().size2()] );
4740  for ( int l = 0; l < basis_type::nDof; ++l )
4741  {
4742  const int ncdof = is_product?nComponents1:1;
4743 
4744  for ( typename array_type::index c1 = 0; c1 < ncdof; ++c1 )
4745  {
4746  typename array_type::index ldof = basis_type::nDof*c1+l;
4747  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( elt_id, l, c1 ) );
4748  //std::cout << "ldof = " << ldof << "\n";
4749  //std::cout << "gdof = " << gdof << "\n";
4750 
4751  FEELPP_ASSERT( gdof < this->size() )
4752  ( context.eId() )
4753  ( l )( c1 )( ldof )( gdof )
4754  ( this->size() )( this->localSize() )
4755  ( this->firstLocalIndex() )( this->lastLocalIndex() )
4756  .warn( "FunctionSpace::Element invalid access index" );
4757 
4758  value_type v_ = this->globalValue( gdof );
4759 
4760  //std::cout << "v_ =" << v_ << "\n";
4761  //for( typename array_type::index c2 = 0; c2 < nComponents2; ++c2 )
4762  for ( uint16_type q = 0; q < nq; ++q )
4763  {
4764  for ( typename array_type::index i = 0; i < nComponents1; ++i )
4765  //for( typename array_type::index j = 0; j < nComponents2; ++j )
4766  {
4767  v[q]( i,0 ) += v_*context.id( ldof, i, 0, q );
4768  //v[q](i,0) += v_*context.gmc()->J(*)*context.pc()->phi( ldof, i, 0, q );
4769  }
4770  }
4771  }
4772  }
4773 
4774  //return v;
4775 }
4776 
4777 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4778 template<typename Y, typename Cont>
4779 void
4780 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::idInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const
4781 {
4782 
4783  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
4784  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
4785  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
4786 
4787  // create analysys map : id -> List of pt
4788  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
4789  __loc->run_analysis( __ptsReal,invalid_size_type_value );
4790  analysis_iterator_type it = __loc->result_analysis_begin();
4791  analysis_iterator_type it_end = __loc->result_analysis_end();
4792  analysis_output_iterator_type itL,itL_end;
4793 
4794  //geomap
4795  gm_ptrtype __gm = this->functionSpace()->gm();
4796 
4797  //if analysis map is empty : no interpolation
4798  if ( it==it_end ) return;
4799 
4800  //alocate a point matrix
4801  size_type nbPtsElt=it->second.size();
4802  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
4803  matrix_node_type pts( nbCoord, nbPtsElt );
4804 
4805  //init the geomap context and precompute basis function
4806  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
4807  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
4808 
4809  gmc_ptrtype __c( new gmc_type( __gm,
4810  this->functionSpace()->mesh()->element( it->first ),
4811  __geopc ) );
4812 
4813  typedef typename mesh_type::element_type geoelement_type;
4814  typedef typename functionspace_type::fe_type fe_type;
4815  typedef typename fe_type::template Context<vm::JACOBIAN|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
4816  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
4817  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
4818  __c,
4819  __pc ) );
4820 
4821  for ( ; it!=it_end; ++it )
4822  {
4823  nbPtsElt = it->second.size();
4824 
4825  //iterate in the list pt for a element
4826  itL=it->second.begin();
4827  itL_end=it->second.end();
4828 
4829  //compute a point matrix with the list of point
4830  pts= matrix_node_type( nbCoord, nbPtsElt );
4831 
4832  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
4833  ublas::column( pts, i ) = boost::get<1>( *itL );
4834 
4835  //update geomap context
4836  __geopc->update( pts );
4837 
4838  //update precompute of basis functions
4839  __pc->update( pts );
4840 
4841  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
4842  __ctx->update( __c, __pc );
4843 
4844  //evaluate element for these points
4845  id_type __id( this->id( *__ctx ) );
4846 
4847  //update the output data
4848  itL=it->second.begin();
4849 
4850  for ( uint k=0; k<nbPtsElt; ++k,++itL )
4851  {
4852  for ( typename array_type::index i = 0; i < nComponents1; ++i )
4853  {
4854  v[boost::get<0>( *itL )]( i,0 ) = __id( i,0,k );
4855  }
4856  }
4857  }
4858 
4859 }
4860 
4861 //
4862 // Grad
4863 //
4864 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4865 template<typename Y, typename Cont>
4866 typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::grad_type
4867 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::grad( node_type const& __x ) const
4868 {
4869  this->updateGlobalValues();
4870 
4871  node_type __x_ref;
4872  size_type __cv_id;
4873  std::vector<int> found_pt( functionSpace()->mesh()->comm().size(), 0 );
4874  std::vector<int> global_found_pt( functionSpace()->mesh()->comm().size(), 0 );
4875 
4876  if ( functionSpace()->findPoint( __x, __cv_id, __x_ref ) )
4877  {
4878 #if !defined( NDEBUG )
4879  DVLOG(2) << "Point " << __x << " is in element " << __cv_id << " pt_ref=" << __x_ref << "\n";
4880 #endif
4881  gm_ptrtype __gm = functionSpace()->gm();
4882  typedef typename gm_type::precompute_ptrtype geopc_ptrtype;
4883  typedef typename gm_type::precompute_type geopc_type;
4884  typename matrix_node<value_type>::type pts( __x_ref.size(), 1 );
4885  ublas::column( pts, 0 ) = __x_ref;
4886  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
4887 
4888 
4889  typedef typename gm_type::template Context<vm::POINT|vm::JACOBIAN|vm::KB, geoelement_type> gmc_type;
4890  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
4891  gmc_ptrtype __c( new gmc_type( __gm,
4892  functionSpace()->mesh()->element( __cv_id ),
4893  __geopc ) );
4894  pc_ptrtype pc( new pc_type( this->functionSpace()->fe(), pts ) );
4895  typedef typename mesh_type::element_type geoelement_type;
4896  typedef typename functionspace_type::fe_type fe_type;
4897  typedef typename fe_type::template Context<vm::GRAD, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
4898  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
4899  fectx_ptrtype fectx( new fectx_type( this->functionSpace()->fe(),
4900  __c,
4901  pc ) );
4902 
4903  found_pt[ functionSpace()->mesh()->comm().rank() ] = 1;
4904 
4905 #if defined(FEELPP_HAS_MPI)
4906 
4907  if ( functionSpace()->mesh()->comm().size() > 1 )
4908  {
4909  //mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, std::plus<std::vector<int> >() );
4910  mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4911  }
4912 
4913 #else
4914  global_found_pt[ 0 ] = found_pt[ 0 ];
4915 #endif /* FEELPP_HAS_MPI */
4916 
4917  grad_type g_( this->grad( *fectx ) );
4918  //DVLOG(2) << "[interpolation] id = " << v << "\n";
4919 #if defined(FEELPP_HAS_MPI)
4920  DVLOG(2) << "sending interpolation context to all processors from " << functionSpace()->mesh()->comm().rank() << "\n";
4921 
4922  if ( functionSpace()->mesh()->comm().size() > 1 )
4923  {
4924  mpi::broadcast( functionSpace()->mesh()->comm(), g_, functionSpace()->mesh()->comm().rank() );
4925  }
4926 
4927  //DVLOG(2) << "[interpolation] after broadcast g_ = " << g_ << "\n";
4928 #endif /* FEELPP_HAS_MPI */
4929  return g_;
4930  }
4931 
4932  else
4933  {
4934 #if defined(FEELPP_HAS_MPI)
4935 
4936  if ( functionSpace()->mesh()->comm().size() > 1 )
4937  {
4938  //mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, std::plus<std::vector<int> >() );
4939  mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4940  }
4941 
4942 #endif /* FEELPP_HAS_MPI */
4943  bool found = false;
4944  size_type i = 0;
4945 
4946  for ( ; i < global_found_pt.size(); ++i )
4947  if ( global_found_pt[i] != 0 )
4948  {
4949  DVLOG(2) << "processor " << i << " has the point " << __x << "\n";
4950  found = true;
4951  break;
4952  }
4953 
4954  grad_type g_;
4955 
4956  if ( found )
4957  {
4958  DVLOG(2) << "receiving interpolation context from processor " << i << "\n";
4959 #if defined(FEELPP_HAS_MPI)
4960 
4961  if ( functionSpace()->mesh()->comm().size() > 1 )
4962  mpi::broadcast( functionSpace()->mesh()->comm(), g_, i );
4963 
4964 #endif /* FEELPP_HAS_MPI */
4965 
4966  //DVLOG(2) << "[interpolation] after broadcast id = " << v << "\n";
4967  }
4968 
4969  else
4970  {
4971  Warning() << "no processor seems to have the point " << __x << "\n";
4972  }
4973 
4974  return g_;
4975  }
4976 
4977 } // grad
4978 template<typename A0, typename A1, typename A2, typename A3, typename A4>
4979 template<typename Y, typename Cont>
4980 template<typename ContextType>
4981 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
4982 void
4983 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::grad_( ContextType const & context, grad_array_type& v ) const
4984 {
4985  if ( !this->areGlobalValuesUpdated() )
4986  this->updateGlobalValues();
4987 
4988  size_type elt_id = context.eId();
4989  if ( context.gmContext()->element().mesh()->isSubMeshFrom( this->mesh() ) )
4990  elt_id = context.gmContext()->element().mesh()->subMeshToMesh( context.eId() );
4991  if ( context.gmContext()->element().mesh()->isParentMeshOf( this->mesh() ) )
4992  elt_id = this->mesh()->meshToSubMesh( context.eId() );
4993  if ( elt_id == invalid_size_type_value )
4994  return;
4995 
4996  //std::cout << "coeff=" << coeff << "\n";
4997  //array_type v( boost::extents[nComponents1][nRealDim][context.xRefs().size2()] );
4998  //std::fill( v.data(), v.data()+v.num_elements(), value_type( 0 ) );
4999  for ( int l = 0; l < basis_type::nDof; ++l )
5000  {
5001  const int ncdof = is_product?nComponents1:1;
5002 
5003  for ( int c1 = 0; c1 < ncdof; ++c1 )
5004  {
5005  int ldof = c1*basis_type::nDof+l;
5006  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( elt_id, l, c1 ) );
5007  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5008  gdof < this->lastLocalIndex() )
5009  ( context.eId() )
5010  ( l )( c1 )( ldof )( gdof )
5011  ( this->size() )( this->localSize() )
5012  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5013  .error( "FunctionSpace::Element invalid access index" );
5014 
5015  //value_type v_ = (*this)( gdof );
5016  value_type v_ = this->globalValue( gdof );
5017 
5018  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5019  {
5020  for ( int k = 0; k < nComponents1; ++k )
5021  for ( int j = 0; j < nRealDim; ++j )
5022  {
5023  v[q]( k,j ) += v_*context.grad( ldof, k, j, q );
5024  }
5025  }
5026  }
5027  }
5028 
5029 #if 0
5030  std::cout << "xref = " << context.xRefs() << "\n";
5031  std::cout << "xreal = " << context.xReal() << "\n";
5032  std::cout << "pts = " << context.G() << "\n";
5033 
5034  for ( int l = 0; l < basis_type::nDof; ++l )
5035  {
5036  for ( typename array_type::index i = 0; i < nDim; ++i )
5037  {
5038  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5039  std::cout << "Gref[ " << l << "," << i << "," << q << "] = " << pc.grad( l, 0, i, q ) << "\n";
5040 
5041  }
5042  }
5043 
5044  for ( int l = 0; l < basis_type::nDof; ++l )
5045  {
5046  for ( int j = 0; j < nRealDim; ++j )
5047  for ( typename array_type::index i = 0; i < nDim; ++i )
5048  {
5049  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5050  std::cout << "G[ " << l << "," << j << "," << i << "," << q << "] = " << context.B( q )( j, i )*pc.grad( l, 0, i, q ) << "\n";
5051 
5052  }
5053  }
5054 
5055  for ( int k = 0; k < nComponents1; ++k )
5056  for ( int j = 0; j < nRealDim; ++j )
5057  {
5058  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5059  {
5060  std::cout << "v[ " << k << "," << j << "," << q << "]= " << v[q]( k,j ) << "\n";
5061  std::cout << "B(" << q << ")=" << context.B( q ) << "\n";
5062  std::cout << "J(" << q << ")=" << context.J( q ) << "\n";
5063  std::cout << "K(" << q << ")=" << context.K( q ) << "\n";
5064  }
5065  }
5066 
5067 #endif
5068 }
5069 
5070 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5071 template<typename Y, typename Cont>
5072 void
5073 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::gradInterpolate( matrix_node_type __ptsReal, grad_array_type& v ) const
5074 {
5075  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5076  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5077  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5078 
5079  // create analysys map : id -> List of pt
5080  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5081  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5082  analysis_iterator_type it = __loc->result_analysis_begin();
5083  analysis_iterator_type it_end = __loc->result_analysis_end();
5084  analysis_output_iterator_type itL,itL_end;
5085 
5086  gm_ptrtype __gm = this->functionSpace()->gm();
5087 
5088  //if analysis map is empty : no interpolation
5089  if ( it==it_end ) return;
5090 
5091  size_type nbPtsElt=it->second.size();
5092  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5093  matrix_node_type pts( nbCoord, nbPtsElt );
5094  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5095  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5096 
5097  gmc_ptrtype __c( new gmc_type( __gm,
5098  this->functionSpace()->mesh()->element( it->first ),
5099  __geopc ) );
5100  typedef typename mesh_type::element_type geoelement_type;
5101  typedef typename functionspace_type::fe_type fe_type;
5102  typedef typename fe_type::template Context<vm::JACOBIAN|vm::KB|vm::GRAD|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5103  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5104  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5105  __c,
5106  __pc ) );
5107 
5108  for ( ; it!=it_end; ++it )
5109  {
5110  nbPtsElt = it->second.size();
5111 
5112  //iterate in the list pt for a element
5113  itL=it->second.begin();
5114  itL_end=it->second.end();
5115 
5116  //compute a point matrix with the list of point
5117  pts= matrix_node_type( nbCoord, nbPtsElt );
5118 
5119  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5120  ublas::column( pts, i ) = boost::get<1>( *itL );
5121 
5122  //update geomap context
5123  __geopc->update( pts );
5124 
5125 
5126  //update precompute of basis functions
5127  __pc->update( pts );
5128  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5129  __ctx->update( __c, __pc );
5130 
5131  //evaluate element for these points
5132  grad_type __grad( this->grad( *__ctx ) );
5133 
5134  //update the output data
5135  itL=it->second.begin();
5136 
5137  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5138  {
5139  for ( typename array_type::index i = 0; i < nComponents1; ++i )
5140  for ( uint j = 0; j < nRealDim; ++j )
5141  {
5142  v[boost::get<0>( *itL )]( i,j ) = __grad( i,j,k );
5143  }
5144  }
5145 
5146  }
5147 
5148 }
5149 
5150 
5151 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5152 template<typename Y, typename Cont>
5153 template<typename ContextType>
5154 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
5155 void
5156 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::div_( ContextType const & context, div_array_type& v ) const
5157 {
5158 #if 1
5159 
5160  if ( !this->areGlobalValuesUpdated() )
5161  this->updateGlobalValues();
5162 
5163  size_type elt_id = context.eId();
5164  if ( context.gmContext()->element().mesh()->isSubMeshFrom( this->mesh() ) )
5165  elt_id = context.gmContext()->element().mesh()->subMeshToMesh( context.eId() );
5166  if ( context.gmContext()->element().mesh()->isParentMeshOf( this->mesh() ) )
5167  elt_id = this->mesh()->meshToSubMesh( context.eId() );
5168  if ( elt_id == invalid_size_type_value )
5169  return;
5170 
5171  const size_type Q = context.xRefs().size2();
5172 
5173  for ( int l = 0; l < basis_type::nDof; ++l )
5174  {
5175  const int ncdof = is_product?nComponents1:1;
5176 
5177  for ( int c1 = 0; c1 < ncdof; ++c1 )
5178  {
5179  int ldof = c1*basis_type::nDof+l;
5180  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( elt_id, l, c1 ) );
5181  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5182  gdof < this->lastLocalIndex() )
5183  ( context.eId() )
5184  ( l )( c1 )( ldof )( gdof )
5185  ( this->size() )( this->localSize() )
5186  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5187  .error( "FunctionSpace::Element invalid access index" );
5188  //value_type v_ = (*this)( gdof );
5189  value_type v_ = this->globalValue( gdof );
5190 
5191  for ( size_type q = 0; q < Q; ++q )
5192  {
5193 #if 0
5194  std::cout << "v(" << gdof << ")=" << v_ << "\n";
5195  std::cout << "context.div(" << ldof << "," << q << ")="
5196  << context.div( ldof, 0, 0, q ) << "\n" ;
5197 #endif
5198  v[q]( 0,0 ) += v_*context.div( ldof, 0, 0, q );
5199  }
5200  }
5201  }
5202 
5203 #else
5204 
5205  if ( !this->areGlobalValuesUpdated() )
5206  this->updateGlobalValues();
5207 
5208  for ( int l = 0; l < basis_type::nDof; ++l )
5209  {
5210  const int ncdof = is_product?nComponents1:1;
5211 
5212  for ( int c1 = 0; c1 < ncdof; ++c1 )
5213  {
5214  int ldof = c1*basis_type::nDof+l;
5215  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( context.eId(), l, c1 ) );
5216  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5217  gdof < this->lastLocalIndex() )
5218  ( context.eId() )
5219  ( l )( c1 )( ldof )( gdof )
5220  ( this->size() )( this->localSize() )
5221  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5222  .error( "FunctionSpace::Element invalid access index" );
5223  //value_type v_ = (*this)( gdof );
5224  value_type v_ = this->globalValue( gdof );
5225 
5226  for ( int k = 0; k < nComponents1; ++k )
5227  {
5228  for ( typename array_type::index i = 0; i < nDim; ++i )
5229  {
5230  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5231  {
5232  v[q]( 0,0 ) += v_*context.gmContext()->B( q )( k, i )*context.pc()->grad( ldof, k, i, q );
5233  }
5234  }
5235  }
5236  }
5237  }
5238 
5239 #endif
5240 }
5241 
5242 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5243 template<typename Y, typename Cont>
5244 void
5245 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::divInterpolate( matrix_node_type __ptsReal, div_array_type& v ) const
5246 {
5247 
5248  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5249  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5250  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5251 
5252  // create analysys map : id -> List of pt
5253  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5254  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5255  analysis_iterator_type it = __loc->result_analysis_begin();
5256  analysis_iterator_type it_end = __loc->result_analysis_end();
5257  analysis_output_iterator_type itL,itL_end;
5258 
5259  //geomap
5260  gm_ptrtype __gm = this->functionSpace()->gm();
5261 
5262  //if analysis map is empty : no interpolation
5263  if ( it==it_end ) return;
5264 
5265  //alocate a point matrix
5266  size_type nbPtsElt=it->second.size();
5267  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5268  matrix_node_type pts( nbCoord, nbPtsElt );
5269 
5270  //init the geomap context and precompute basis function
5271  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5272  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5273  gmc_ptrtype __c( new gmc_type( __gm,
5274  this->functionSpace()->mesh()->element( it->first ),
5275  __geopc ) );
5276  typedef typename mesh_type::element_type geoelement_type;
5277  typedef typename functionspace_type::fe_type fe_type;
5278  typedef typename fe_type::template Context<vm::DIV|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5279  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5280  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5281  __c,
5282  __pc ) );
5283 
5284  for ( ; it!=it_end; ++it )
5285  {
5286  nbPtsElt = it->second.size();
5287 
5288  //iterate in the list pt for a element
5289  itL=it->second.begin();
5290  itL_end=it->second.end();
5291 
5292  //compute a point matrix with the list of point
5293  pts= matrix_node_type( nbCoord, nbPtsElt );
5294 
5295  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5296  ublas::column( pts, i ) = boost::get<1>( *itL );
5297 
5298  //update geomap context
5299  __geopc->update( pts );
5300  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5301  //update precompute of basis functions
5302  __pc->update( pts );
5303  __ctx->update( __c, __pc );
5304 
5305  //evaluate element for these points
5306  div_type __div( this->div( *__ctx ) );
5307 
5308  //update the output data
5309  //for( typename array_type::index i = 0; i < nComponents1; ++i )
5310  // {
5311  itL=it->second.begin();
5312 
5313  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5314  v[boost::get<0>( *itL )]( 0,0 ) = __div( 0,0,k );
5315 
5316  // }
5317  }
5318 
5319 }
5320 
5321 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5322 template<typename Y, typename Cont>
5323 template<typename ContextType>
5324 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
5325 void
5326 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curl_( ContextType const & context, curl_array_type& v ) const
5327 {
5328  if ( !this->areGlobalValuesUpdated() )
5329  this->updateGlobalValues();
5330 
5331  size_type elt_id = context.eId();
5332  if ( context.gmContext()->element().mesh()->isSubMeshFrom( this->mesh() ) )
5333  elt_id = context.gmContext()->element().mesh()->subMeshToMesh( context.eId() );
5334  if ( context.gmContext()->element().mesh()->isParentMeshOf( this->mesh() ) )
5335  elt_id = this->mesh()->meshToSubMesh( context.eId() );
5336  if ( elt_id == invalid_size_type_value )
5337  return;
5338 
5339  for ( int l = 0; l < basis_type::nDof; ++l )
5340  {
5341  const int ncdof = is_product?nComponents1:1;
5342 
5343  for ( int c1 = 0; c1 < ncdof; ++c1 )
5344  {
5345  int ldof = c1*basis_type::nDof+l;
5346  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( elt_id, l, c1 ) );
5347  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5348  gdof < this->lastLocalIndex() )
5349  ( context.eId() )
5350  ( l )( c1 )( ldof )( gdof )
5351  ( this->size() )( this->localSize() )
5352  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5353  .error( "FunctionSpace::Element invalid access index" );
5354 
5355  //value_type v_ = (*this)( gdof );
5356  value_type v_ = this->globalValue( gdof );
5357  const uint16_type nq = context.xRefs().size2();
5358 
5359  for ( uint16_type q = 0; q < nq ; ++q )
5360  {
5361  if ( nDim == 3 )
5362  {
5363  for ( typename array_type::index i = 0; i < nDim; ++i )
5364  {
5365  v[q]( i,0 ) += v_*context.curl( ldof, i, 0, q );
5366  }
5367  }
5368 
5369  else if ( nDim == 2 )
5370  {
5371  v[q]( 0,0 ) += v_*context.curl( ldof, 0, 0, q );
5372  }
5373 
5374  }
5375  }
5376  }
5377 
5378 }
5379 
5380 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5381 template<typename Y, typename Cont>
5382 template<typename ContextType>
5383 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
5384 void
5385 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curl_( ContextType const & context, comp_curl_array_type& v, int comp ) const
5386 {
5387  if ( !this->areGlobalValuesUpdated() )
5388  this->updateGlobalValues();
5389 
5390  size_type elt_id = context.eId();
5391  if ( context.gmContext()->element().mesh()->isSubMeshFrom( this->mesh() ) )
5392  elt_id = context.gmContext()->element().mesh()->subMeshToMesh( context.eId() );
5393  if ( context.gmContext()->element().mesh()->isParentMeshOf( this->mesh() ) )
5394  elt_id = this->mesh()->meshToSubMesh( context.eId() );
5395  if ( elt_id == invalid_size_type_value )
5396  return;
5397 
5398 
5399  for ( int l = 0; l < basis_type::nDof; ++l )
5400  {
5401  const int ncdof = is_product?nComponents1:1;
5402 
5403  for ( int c1 = 0; c1 < ncdof; ++c1 )
5404  {
5405  int ldof = c1*basis_type::nDof+l;
5406  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( elt_id, l, c1 ) );
5407  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5408  gdof < this->lastLocalIndex() )
5409  ( context.eId() )
5410  ( l )( c1 )( ldof )( gdof )
5411  ( this->size() )( this->localSize() )
5412  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5413  .error( "FunctionSpace::Element invalid access index" );
5414 
5415  //value_type v_ = (*this)( gdof );
5416  value_type v_ = this->globalValue( gdof );
5417  const uint16_type nq = context.xRefs().size2();
5418 
5419  for ( uint16_type q = 0; q < nq ; ++q )
5420  {
5421  if ( nDim == 3 )
5422  {
5423  v[q]( 0,0 ) += v_*context.curl( ldof, comp, 0, q );
5424  }
5425 
5426  else if ( nDim == 2 )
5427  {
5428  v[q]( 0,0 ) += v_*context.curl( ldof, 2, 0, q );
5429  }
5430 
5431  }
5432  }
5433  }
5434 
5435 }
5436 
5437 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5438 template<typename Y, typename Cont>
5439 void
5440 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlInterpolate( matrix_node_type __ptsReal, curl_array_type& v ) const
5441 {
5442 
5443  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5444  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5445  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5446 
5447  // create analysys map : id -> List of pt
5448  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5449  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5450  analysis_iterator_type it = __loc->result_analysis_begin();
5451  analysis_iterator_type it_end = __loc->result_analysis_end();
5452  analysis_output_iterator_type itL,itL_end;
5453 
5454  //geomap
5455  gm_ptrtype __gm = this->functionSpace()->gm();
5456 
5457  //if analysis map is empty : no interpolation
5458  if ( it==it_end ) return;
5459 
5460  //alocate a point matrix
5461  size_type nbPtsElt=it->second.size();
5462  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5463  matrix_node_type pts( nbCoord, nbPtsElt );
5464 
5465  //init the geomap context and precompute basis function
5466  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5467  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5468  gmc_ptrtype __c( new gmc_type( __gm,
5469  this->functionSpace()->mesh()->element( it->first ),
5470  __geopc ) );
5471  typedef typename mesh_type::element_type geoelement_type;
5472  typedef typename functionspace_type::fe_type fe_type;
5473  typedef typename fe_type::template Context<vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5474  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5475  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5476  __c,
5477  __pc ) );
5478 
5479 
5480  for ( ; it!=it_end; ++it )
5481  {
5482  nbPtsElt = it->second.size();
5483 
5484  //iterate in the list pt for a element
5485  itL=it->second.begin();
5486  itL_end=it->second.end();
5487 
5488  //compute a point matrix with the list of point
5489  pts= matrix_node_type( nbCoord, nbPtsElt );
5490 
5491  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5492  ublas::column( pts, i ) = boost::get<1>( *itL );
5493 
5494  //update geomap context
5495  __geopc->update( pts );
5496  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5497  //update precompute of basis functions
5498  __pc->update( pts );
5499  __ctx->update( __c, __pc );
5500 
5501  //evaluate element for these points
5502  curl_type __curl( this->curl( *__ctx ) );
5503 
5504  //update the output data
5505  itL=it->second.begin();
5506 
5507  if ( nDim == 3 )
5508  {
5509  itL=it->second.begin();
5510 
5511  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5512  {
5513  v[boost::get<0>( *itL )]( 0,0 ) = __curl( 0,0,k );
5514  v[boost::get<0>( *itL )]( 1,0 ) = __curl( 1,0,k );
5515  v[boost::get<0>( *itL )]( 2,0 ) = __curl( 2,0,k );
5516  }
5517  }
5518 
5519  else if ( nDim == 2 )
5520  {
5521  itL=it->second.begin();
5522 
5523  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5524  v[boost::get<0>( *itL )]( 0,0 ) = __curl( 0,0,k );
5525  }
5526  }
5527 
5528 }
5529 
5530 
5531 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5532 template<typename Y, typename Cont>
5533 void
5534 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlxInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const
5535 {
5536 
5537  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5538  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5539  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5540 
5541  // create analysys map : id -> List of pt
5542  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5543  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5544  analysis_iterator_type it = __loc->result_analysis_begin();
5545  analysis_iterator_type it_end = __loc->result_analysis_end();
5546  analysis_output_iterator_type itL,itL_end;
5547 
5548  //geomap
5549  gm_ptrtype __gm = this->functionSpace()->gm();
5550 
5551  //if analysis map is empty : no interpolation
5552  if ( it==it_end ) return;
5553 
5554  //alocate a point matrix
5555  size_type nbPtsElt=it->second.size();
5556  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5557  matrix_node_type pts( nbCoord, nbPtsElt );
5558 
5559  //init the geomap context and precompute basis function
5560  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5561  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5562  gmc_ptrtype __c( new gmc_type( __gm,
5563  this->functionSpace()->mesh()->element( it->first ),
5564  __geopc ) );
5565  typedef typename mesh_type::element_type geoelement_type;
5566  typedef typename functionspace_type::fe_type fe_type;
5567  typedef typename fe_type::template Context<vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5568  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5569  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5570  __c,
5571  __pc ) );
5572 
5573 
5574  for ( ; it!=it_end; ++it )
5575  {
5576  nbPtsElt = it->second.size();
5577 
5578  //iterate in the list pt for a element
5579  itL=it->second.begin();
5580  itL_end=it->second.end();
5581 
5582  //compute a point matrix with the list of point
5583  pts= matrix_node_type( nbCoord, nbPtsElt );
5584 
5585  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5586  ublas::column( pts, i ) = boost::get<1>( *itL );
5587 
5588  //update geomap context
5589  __geopc->update( pts );
5590  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5591  //update precompute of basis functions
5592  __pc->update( pts );
5593  __ctx->update( __c, __pc );
5594 
5595  //evaluate element for these points
5596  curlx_type __curlx( this->curlx( *__ctx ) );
5597 
5598  //update the output data
5599  itL=it->second.begin();
5600 
5601  if ( nDim == 3 )
5602  {
5603  itL=it->second.begin();
5604 
5605  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5606  {
5607  v[boost::get<0>( *itL )]( 0,0 ) = __curlx( 0,0,k );
5608  v[boost::get<0>( *itL )]( 1,0 ) = __curlx( 1,0,k );
5609  v[boost::get<0>( *itL )]( 2,0 ) = __curlx( 2,0,k );
5610  }
5611  }
5612 
5613  else if ( nDim == 2 )
5614  {
5615  itL=it->second.begin();
5616 
5617  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5618  v[boost::get<0>( *itL )]( 0,0 ) = __curlx( 0,0,k );
5619  }
5620  }
5621 
5622 }
5623 
5624 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5625 template<typename Y, typename Cont>
5626 void
5627 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlyInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const
5628 {
5629 
5630  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5631  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5632  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5633 
5634  // create analysys map : id -> List of pt
5635  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5636  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5637  analysis_iterator_type it = __loc->result_analysis_begin();
5638  analysis_iterator_type it_end = __loc->result_analysis_end();
5639  analysis_output_iterator_type itL,itL_end;
5640 
5641  //geomap
5642  gm_ptrtype __gm = this->functionSpace()->gm();
5643 
5644  //if analysis map is empty : no interpolation
5645  if ( it==it_end ) return;
5646 
5647  //alocate a point matrix
5648  size_type nbPtsElt=it->second.size();
5649  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5650  matrix_node_type pts( nbCoord, nbPtsElt );
5651 
5652  //init the geomap context and precompute basis function
5653  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5654  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5655  gmc_ptrtype __c( new gmc_type( __gm,
5656  this->functionSpace()->mesh()->element( it->first ),
5657  __geopc ) );
5658  typedef typename mesh_type::element_type geoelement_type;
5659  typedef typename functionspace_type::fe_type fe_type;
5660  typedef typename fe_type::template Context<vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5661  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5662  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5663  __c,
5664  __pc ) );
5665 
5666  for ( ; it!=it_end; ++it )
5667  {
5668  nbPtsElt = it->second.size();
5669 
5670  //iterate in the list pt for a element
5671  itL=it->second.begin();
5672  itL_end=it->second.end();
5673 
5674  //compute a point matrix with the list of point
5675  pts= matrix_node_type( nbCoord, nbPtsElt );
5676 
5677  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5678  ublas::column( pts, i ) = boost::get<1>( *itL );
5679 
5680  //update geomap context
5681  __geopc->update( pts );
5682  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5683  //update precompute of basis functions
5684  __pc->update( pts );
5685  __ctx->update( __c, __pc );
5686 
5687  //evaluate element for these points
5688  curly_type __curly( this->curly( *__ctx ) );
5689 
5690  //update the output data
5691  itL=it->second.begin();
5692 
5693  if ( nDim == 3 )
5694  {
5695  itL=it->second.begin();
5696 
5697  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5698  {
5699  v[boost::get<0>( *itL )]( 0,0 ) = __curly( 0,0,k );
5700  v[boost::get<0>( *itL )]( 1,0 ) = __curly( 1,0,k );
5701  v[boost::get<0>( *itL )]( 2,0 ) = __curly( 2,0,k );
5702  }
5703  }
5704 
5705  else if ( nDim == 2 )
5706  {
5707  itL=it->second.begin();
5708 
5709  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5710  v[boost::get<0>( *itL )]( 0,0 ) = __curly( 0,0,k );
5711  }
5712  }
5713 
5714 }
5715 
5716 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5717 template<typename Y, typename Cont>
5718 void
5719 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlzInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v ) const
5720 {
5721 
5722  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5723  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5724  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5725 
5726  // create analysys map : id -> List of pt
5727  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5728  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5729  analysis_iterator_type it = __loc->result_analysis_begin();
5730  analysis_iterator_type it_end = __loc->result_analysis_end();
5731  analysis_output_iterator_type itL,itL_end;
5732 
5733  //geomap
5734  gm_ptrtype __gm = this->functionSpace()->gm();
5735 
5736  //if analysis map is empty : no interpolation
5737  if ( it==it_end ) return;
5738 
5739  //alocate a point matrix
5740  size_type nbPtsElt=it->second.size();
5741  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5742  matrix_node_type pts( nbCoord, nbPtsElt );
5743 
5744  //init the geomap context and precompute basis function
5745  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5746  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5747  gmc_ptrtype __c( new gmc_type( __gm,
5748  this->functionSpace()->mesh()->element( it->first ),
5749  __geopc ) );
5750  typedef typename mesh_type::element_type geoelement_type;
5751  typedef typename functionspace_type::fe_type fe_type;
5752  typedef typename fe_type::template Context<vm::CURL|vm::JACOBIAN|vm::KB|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5753  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5754  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5755  __c,
5756  __pc ) );
5757 
5758  for ( ; it!=it_end; ++it )
5759  {
5760  nbPtsElt = it->second.size();
5761 
5762  //iterate in the list pt for a element
5763  itL=it->second.begin();
5764  itL_end=it->second.end();
5765 
5766  //compute a point matrix with the list of point
5767  pts= matrix_node_type( nbCoord, nbPtsElt );
5768 
5769  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5770  ublas::column( pts, i ) = boost::get<1>( *itL );
5771 
5772  //update geomap context
5773  __geopc->update( pts );
5774  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5775  //update precompute of basis functions
5776  __pc->update( pts );
5777  __ctx->update( __c, __pc );
5778 
5779  //evaluate element for these points
5780  curlz_type __curlz( this->curlz( *__ctx ) );
5781 
5782  //update the output data
5783  itL=it->second.begin();
5784 
5785  if ( nDim == 3 )
5786  {
5787  itL=it->second.begin();
5788 
5789  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5790  {
5791  v[boost::get<0>( *itL )]( 0,0 ) = __curlz( 0,0,k );
5792  v[boost::get<0>( *itL )]( 1,0 ) = __curlz( 1,0,k );
5793  v[boost::get<0>( *itL )]( 2,0 ) = __curlz( 2,0,k );
5794  }
5795  }
5796 
5797  else if ( nDim == 2 )
5798  {
5799  itL=it->second.begin();
5800 
5801  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5802  v[boost::get<0>( *itL )]( 0,0 ) = __curlz( 0,0,k );
5803  }
5804  }
5805 
5806 }
5807 
5808 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5809 template<typename Y, typename Cont>
5810 template<typename ContextType>
5811 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
5812 void
5813 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::d_( int N, ContextType const & context, id_array_type& v ) const
5814 {
5815  if ( !this->areGlobalValuesUpdated() )
5816  this->updateGlobalValues();
5817 
5818  for ( int i = 0; i < basis_type::nDof; ++i )
5819  {
5820  const int ncdof = is_product?nComponents1:1;
5821 
5822  for ( int c1 = 0; c1 < ncdof; ++c1 )
5823  {
5824  size_type ldof = basis_type::nDof*c1 + i;
5825  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( context.eId(), i, c1 ) );
5826  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
5827  gdof < this->lastLocalIndex() )
5828  ( context.eId() )
5829  ( i )( c1 )( ldof )( gdof )
5830  ( this->size() )( this->localSize() )
5831  ( this->firstLocalIndex() )( this->lastLocalIndex() )
5832  .error( "FunctionSpace::Element invalid access index" );
5833 
5834  value_type v_ = this->globalValue( gdof );
5835 
5836  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
5837  {
5838  for ( typename array_type::index i = 0; i < nComponents1; ++i )
5839  {
5840  v[q]( i,0 ) += v_*context.d( ldof, i, N, q );
5841  }
5842  }
5843 
5844  }
5845  }
5846 }
5847 
5848 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5849 template<typename Y, typename Cont>
5850 void
5851 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dxInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const
5852 {
5853 
5854  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5855  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5856  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5857 
5858  // create analysys map : id -> List of pt
5859  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5860  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5861  analysis_iterator_type it = __loc->result_analysis_begin();
5862  analysis_iterator_type it_end = __loc->result_analysis_end();
5863  analysis_output_iterator_type itL,itL_end;
5864 
5865  //geomap
5866  gm_ptrtype __gm = this->functionSpace()->gm();
5867 
5868  //if analysis map is empty : no interpolation
5869  if ( it==it_end ) return;
5870 
5871  //alocate a point matrix
5872  size_type nbPtsElt=it->second.size();
5873  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5874  matrix_node_type pts( nbCoord, nbPtsElt );
5875 
5876  //init the geomap context and precompute basis function
5877  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5878  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5879  gmc_ptrtype __c( new gmc_type( __gm,
5880  this->functionSpace()->mesh()->element( it->first ),
5881  __geopc ) );
5882  typedef typename mesh_type::element_type geoelement_type;
5883  typedef typename functionspace_type::fe_type fe_type;
5884  typedef typename fe_type::template Context<vm::JACOBIAN|vm::KB|vm::GRAD|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5885  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5886  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5887  __c,
5888  __pc ) );
5889 
5890  for ( ; it!=it_end; ++it )
5891  {
5892  nbPtsElt = it->second.size();
5893 
5894  //iterate in the list pt for a element
5895  itL=it->second.begin();
5896  itL_end=it->second.end();
5897 
5898  //compute a point matrix with the list of point
5899  pts= matrix_node_type( nbCoord, nbPtsElt );
5900 
5901  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5902  ublas::column( pts, i ) = boost::get<1>( *itL );
5903 
5904  //update geomap context
5905  __geopc->update( pts );
5906  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5907  //update precompute of basis functions
5908  __pc->update( pts );
5909  __ctx->update( __c, __pc );
5910 
5911  //evaluate element for these points
5912  dx_type __dx( this->dx( *__ctx ) );
5913 
5914  //update the output data
5915  itL=it->second.begin();
5916 
5917  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5918  for ( typename array_type::index i = 0; i < nComponents1; ++i )
5919  {
5920  v[boost::get<0>( *itL )]( i,0 ) = __dx( i,0,k );
5921  }
5922  }
5923 
5924 }
5925 
5926 template<typename A0, typename A1, typename A2, typename A3, typename A4>
5927 template<typename Y, typename Cont>
5928 void
5929 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dyInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const
5930 {
5931 
5932  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
5933  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
5934  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
5935 
5936  // create analysys map : id -> List of pt
5937  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
5938  __loc->run_analysis( __ptsReal,invalid_size_type_value );
5939  analysis_iterator_type it = __loc->result_analysis_begin();
5940  analysis_iterator_type it_end = __loc->result_analysis_end();
5941  analysis_output_iterator_type itL,itL_end;
5942 
5943  //geomap
5944  gm_ptrtype __gm = this->functionSpace()->gm();
5945 
5946  //if analysis map is empty : no interpolation
5947  if ( it==it_end ) return;
5948 
5949  //alocate a point matrix
5950  size_type nbPtsElt=it->second.size();
5951  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5952  matrix_node_type pts( nbCoord, nbPtsElt );
5953 
5954  //init the geomap context and precompute basis function
5955  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
5956  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
5957  gmc_ptrtype __c( new gmc_type( __gm,
5958  this->functionSpace()->mesh()->element( it->first ),
5959  __geopc ) );
5960  typedef typename mesh_type::element_type geoelement_type;
5961  typedef typename functionspace_type::fe_type fe_type;
5962  typedef typename fe_type::template Context<vm::JACOBIAN|vm::KB|vm::GRAD|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
5963  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
5964  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
5965  __c,
5966  __pc ) );
5967 
5968  for ( ; it!=it_end; ++it )
5969  {
5970  nbPtsElt = it->second.size();
5971 
5972  //iterate in the list pt for a element
5973  itL=it->second.begin();
5974  itL_end=it->second.end();
5975 
5976  //compute a point matrix with the list of point
5977  pts= matrix_node_type( nbCoord, nbPtsElt );
5978 
5979  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
5980  ublas::column( pts, i ) = boost::get<1>( *itL );
5981 
5982  //update geomap context
5983  __geopc->update( pts );
5984  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5985  //update precompute of basis functions
5986  __pc->update( pts );
5987  __ctx->update( __c, __pc );
5988 
5989  //evaluate element for these points
5990  dy_type __dy( this->dy( *__ctx ) );
5991 
5992  //update the output data
5993  itL=it->second.begin();
5994 
5995  for ( uint k=0; k<nbPtsElt; ++k,++itL )
5996  for ( typename array_type::index i = 0; i < nComponents1; ++i )
5997  {
5998  v[boost::get<0>( *itL )]( i,0 ) = __dy( i,0,k );
5999  }
6000  }
6001 
6002 }
6003 
6004 template<typename A0, typename A1, typename A2, typename A3, typename A4>
6005 template<typename Y, typename Cont>
6006 void
6007 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dzInterpolate( matrix_node_type __ptsReal, id_array_type& v ) const
6008 {
6009 
6010  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
6011  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
6012  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
6013 
6014  // create analysys map : id -> List of pt
6015  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
6016  __loc->run_analysis( __ptsReal,invalid_size_type_value );
6017  analysis_iterator_type it = __loc->result_analysis_begin();
6018  analysis_iterator_type it_end = __loc->result_analysis_end();
6019  analysis_output_iterator_type itL,itL_end;
6020 
6021  //geomap
6022  gm_ptrtype __gm = this->functionSpace()->gm();
6023 
6024  //if analysis map is empty : no interpolation
6025  if ( it==it_end ) return;
6026 
6027  //alocate a point matrix
6028  size_type nbPtsElt=it->second.size();
6029  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
6030  matrix_node_type pts( nbCoord, nbPtsElt );
6031 
6032  //init the geomap context and precompute basis function
6033  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
6034  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
6035  gmc_ptrtype __c( new gmc_type( __gm,
6036  this->functionSpace()->mesh()->element( it->first ),
6037  __geopc ) );
6038  typedef typename mesh_type::element_type geoelement_type;
6039  typedef typename functionspace_type::fe_type fe_type;
6040  typedef typename fe_type::template Context<vm::JACOBIAN|vm::KB|vm::GRAD|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
6041  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
6042  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
6043  __c,
6044  __pc ) );
6045 
6046  for ( ; it!=it_end; ++it )
6047  {
6048  nbPtsElt = it->second.size();
6049 
6050  //iterate in the list pt for a element
6051  itL=it->second.begin();
6052  itL_end=it->second.end();
6053 
6054  //compute a point matrix with the list of point
6055  pts= matrix_node_type( nbCoord, nbPtsElt );
6056 
6057  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
6058  ublas::column( pts, i ) = boost::get<1>( *itL );
6059 
6060  //update geomap context
6061  __geopc->update( pts );
6062  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
6063  //update precompute of basis functions
6064  __pc->update( pts );
6065  __ctx->update( __c, __pc );
6066  //evaluate element for these points
6067  dz_type __dz( this->dz( *__ctx ) );
6068 
6069  //update the output data
6070  itL=it->second.begin();
6071 
6072  for ( uint k=0; k<nbPtsElt; ++k,++itL )
6073  for ( typename array_type::index i = 0; i < nComponents1; ++i )
6074  {
6075  v[boost::get<0>( *itL )]( i,0 ) = __dz( i,0,k );
6076  }
6077  }
6078 
6079 }
6080 
6081 template<typename A0, typename A1, typename A2, typename A3, typename A4>
6082 template<typename Y, typename Cont>
6083 template<typename ContextType>
6084 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
6085 void
6086 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hess_( ContextType const & context, hess_array_type& v ) const
6087 {
6088  hess_( context, v, mpl::int_<rank>() );
6089 }
6090 template<typename A0, typename A1, typename A2, typename A3, typename A4>
6091 template<typename Y, typename Cont>
6092 template<typename ContextType>
6093 //typename FunctionSpace<A0, A1, A2, A3, A4>::template Element<Y,Cont>::array_type
6094 void
6095 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hess_( ContextType const & context, hess_array_type& v, mpl::int_<0> ) const
6096 {
6097  if ( !this->areGlobalValuesUpdated() )
6098  this->updateGlobalValues();
6099 
6100  //FEELPP_ASSERT( comp < nRealDim )( comp )( nRealDim ).error( "[FunctionSpace::Element] grad: invalid component" );
6101  for ( int i = 0; i < basis_type::nDof; ++i )
6102  {
6103  const int ncdof = is_product?nComponents1:1;
6104 
6105  for ( int c1 = 0; c1 < ncdof; ++c1 )
6106  {
6107  size_type ldof = basis_type::nDof*c1 + i;
6108  size_type gdof = boost::get<0>( M_functionspace->dof()->localToGlobal( context.eId(), i, c1 ) );
6109  FEELPP_ASSERT( gdof >= this->firstLocalIndex() &&
6110  gdof < this->lastLocalIndex() )
6111  ( context.eId() )
6112  ( i )( c1 )( ldof )( gdof )
6113  ( this->size() )( this->localSize() )
6114  ( this->firstLocalIndex() )( this->lastLocalIndex() )
6115  .error( "FunctionSpace::Element invalid access index" );
6116 
6117  value_type v_ = this->globalValue( gdof );
6118 
6119  for ( size_type q = 0; q < context.xRefs().size2(); ++q )
6120  {
6121  for ( int i = 0; i < nRealDim; ++i )
6122  for ( int j = 0; j < nRealDim; ++j )
6123  {
6124  v[q]( i,j ) += v_*context.hess( ldof, i, j, q );
6125  } // i,j
6126  } // q
6127 
6128  }
6129  }
6130 
6131 }
6132 
6133 template<typename A0, typename A1, typename A2, typename A3, typename A4>
6134 template<typename Y, typename Cont>
6135 void
6136 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hessInterpolate( matrix_node_type __ptsReal, hess_array_type& v ) const
6137 {
6138 
6139  typedef typename mesh_type::Localization::localization_ptrtype localization_ptrtype;
6140  typedef typename mesh_type::Localization::container_search_iterator_type analysis_iterator_type;
6141  typedef typename mesh_type::Localization::container_output_iterator_type analysis_output_iterator_type;
6142 
6143  // create analysys map : id -> List of pt
6144  localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
6145  __loc->run_analysis( __ptsReal,invalid_size_type_value );
6146  analysis_iterator_type it = __loc->result_analysis_begin();
6147  analysis_iterator_type it_end = __loc->result_analysis_end();
6148  analysis_output_iterator_type itL,itL_end;
6149 
6150  //geomap
6151  gm_ptrtype __gm = this->functionSpace()->gm();
6152 
6153  //if analysis map is empty : no interpolation
6154  if ( it==it_end ) return;
6155 
6156  //alocate a point matrix
6157  size_type nbPtsElt=it->second.size();
6158  uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
6159  matrix_node_type pts( nbCoord, nbPtsElt );
6160 
6161  //init the geomap context and precompute basis function
6162  geopc_ptrtype __geopc( new geopc_type( __gm, pts ) );
6163  pc_ptrtype __pc( new pc_type( this->functionSpace()->fe(), pts ) );
6164  gmc_ptrtype __c( new gmc_type( __gm,
6165  this->functionSpace()->mesh()->element( it->first ),
6166  __geopc ) );
6167  typedef typename mesh_type::element_type geoelement_type;
6168  typedef typename functionspace_type::fe_type fe_type;
6169  typedef typename fe_type::template Context<vm::JACOBIAN|vm::KB|vm::HESSIAN|vm::FIRST_DERIVATIVE|vm::POINT, fe_type, gm_type, geoelement_type,gmc_type::context> fectx_type;
6170  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
6171  fectx_ptrtype __ctx( new fectx_type( this->functionSpace()->fe(),
6172  __c,
6173  __pc ) );
6174 
6175  for ( ; it!=it_end; ++it )
6176  {
6177  nbPtsElt = it->second.size();
6178 
6179  //iterate in the list pt for a element
6180  itL=it->second.begin();
6181  itL_end=it->second.end();
6182 
6183  //compute a point matrix with the list of point
6184  pts= matrix_node_type( nbCoord, nbPtsElt );
6185 
6186  for ( size_type i=0; i<nbPtsElt; ++i,++itL )
6187  ublas::column( pts, i ) = boost::get<1>( *itL );
6188 
6189  //update geomap context
6190  __geopc->update( pts );
6191  __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
6192  //update precompute of basis functions
6193  __pc->update( pts );
6194  __ctx->update( __c, __pc );
6195 
6196  //evaluate element for these points
6197  hess_type __hess( this->hess( *__ctx ) );
6198 
6199  //update the output data
6200  itL=it->second.begin();
6201 
6202  for ( uint k=0; k<nbPtsElt; ++k,++itL )
6203  for ( int i = 0; i < nRealDim; ++i )
6204  for ( int j = 0; j < nRealDim; ++j )
6205  v[boost::get<0>( *itL )]( i,j ) = __hess( i,j,k );
6206  }
6207 
6208 }
6209 
6210 template<typename T,int M,int N>
6211 std::ostream&
6212 operator<<( std::ostream& os, detail::ID<T,M,N> const& id )
6213 {
6214  const size_type* shape = id.M_id.shape();
6215 
6216  for ( size_type i = 0; i < shape[0]; ++i )
6217  {
6218  for ( size_type j = 0; j < shape[1]; ++j )
6219  {
6220  for ( size_type k = 0; k < shape[2]; ++k )
6221  {
6222  os << id( i, j, k ) << ' ';
6223  }
6224 
6225  os << std::endl;
6226  }
6227 
6228  os << std::endl;
6229  }
6230 
6231  return os;
6232 }
6233 
6241 template <typename ElementType>
6242 ElementType
6243 element_product( ElementType const& v1, ElementType const& v2 )
6244 {
6245  FEELPP_ASSERT( v1.functionSpace() == v2.functionSpace() ).error( "incompatible function spaces" );
6246 
6247  typedef typename type_traits<typename ElementType::value_type>::real_type real_type;
6248 
6249  ElementType _t( v1.functionSpace() );
6250  size_type s = v1.localSize();
6251  size_type start = v1.firstLocalIndex();
6252 
6253  for ( size_type i = 0; i < s; ++i )
6254  _t.operator()( start+i ) = v1.operator()( start + i )* v2.operator()( start + i );
6255 
6256  return _t;
6257 }
6258 
6266 template <typename ElementType>
6267 ElementType
6268 element_div( ElementType const& v1, ElementType const& v2 )
6269 {
6270  FEELPP_ASSERT( v1.functionSpace() == v2.functionSpace() ).error( "incompatible function spaces" );
6271 
6272  typedef typename type_traits<typename ElementType::value_type>::real_type real_type;
6273 
6274  ElementType _t( v1.functionSpace() );
6275  size_type s = v1.localSize();
6276  size_type start = v1.firstLocalIndex();
6277 
6278  for ( size_type i = 0; i < s; ++i )
6279  _t.operator()( start+i ) = v1.operator()( start + i )/v2.operator()( start + i );
6280 
6281  return _t;
6282 }
6283 
6284 template<typename MeshType>
6285 auto
6286 measurePointElements( boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<MeshType::nOrder,Scalar> > > >& _space ) -> decltype( _space->element() )
6287 {
6288  auto _fn = _space->element( "measurePointElements" );
6289  _fn.setZero();
6290  std::vector<bool> ptdone( _space->mesh()->numPoints(), false );
6291  auto elit = _space->mesh()->beginElement();
6292  auto elen = _space->mesh()->endElement();
6293 
6294  for ( ; elit != elen; ++ elit )
6295  {
6296  for ( int p = 0; p < elit->numPoints; ++p )
6297  {
6298  if ( ptdone[elit->point( p ).id()] == false )
6299  {
6300  BOOST_FOREACH( auto pt, elit->point( p ).elements() )
6301  {
6302  _fn.plus_assign( elit->id(), p, 0, elit->measure() );
6303  }
6304  ptdone[elit->point( p ).id()] = true;
6305  }
6306  }
6307  }
6308 
6309  return _fn;
6310 }
6311 
6312 template<typename EltType>
6313 typename fusion::result_of::accumulate<typename EltType::functionspace_type::functionspace_vector_type,
6314  fusion::vector<>,
6315  detail::CreateElementVector<EltType> >::type
6316 subelements( EltType const& e, std::vector<std::string> const& n )
6317 {
6318  return fusion::accumulate( e.functionSpaces(), fusion::vector<>(), detail::CreateElementVector<EltType>( e, n ) );
6319 }
6320 
6321 template<typename ElementType, typename CoeffType>
6322 ElementType
6323 expansion( std::vector<ElementType> const& b, CoeffType const& c, int M = -1 )
6324 {
6325  auto res = b[0].functionSpace()->element();
6326  res.zero();
6327  if ( M == -1 ) M = c.size() ;
6328  //LOG_ASSERT( b.size() == c.size() ) << " b.size=" << b.size() << " c.size=" << c.size() << "\n";
6329  for( int i = 0; i < M; ++i )
6330  {
6331  res.add( c[i], b[i] );
6332  }
6333 
6334  return res;
6335 }
6336 
6343 template<int Order,typename MeshType>
6344 inline
6345 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Continuous>>,Periodicity <NoPeriodicity>>>
6346 Pch( boost::shared_ptr<MeshType> mesh )
6347 {
6348  return FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Continuous>>, Periodicity <NoPeriodicity>>::New( mesh );
6349 }
6350 
6357 template<int Order,typename MeshType>
6358 inline
6359 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Discontinuous>>>>
6360 Pdh( boost::shared_ptr<MeshType> mesh )
6361 {
6362  return FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Discontinuous>>>::New( mesh );
6363 }
6364 
6371 template<int Order,typename MeshType>
6372 inline
6373 boost::shared_ptr<FunctionSpace<MeshType,bases<OrthonormalPolynomialSet<Order,Scalar>>>>
6374 Odh( boost::shared_ptr<MeshType> mesh )
6375 {
6376  return FunctionSpace<MeshType,bases<OrthonormalPolynomialSet<Order,Scalar>>>::New( mesh );
6377 }
6378 
6385 template<int Order,typename MeshType>
6386 inline
6387 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Vectorial>>>>
6388 Pchv( boost::shared_ptr<MeshType> mesh )
6389 {
6390  return FunctionSpace<MeshType,bases<Lagrange<Order,Vectorial>>>::New( mesh );
6391 }
6392 
6397 template<int Order,typename MeshType>
6398 inline
6399 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order+1,Vectorial>,Lagrange<Order,Scalar>>>>
6400 THch( boost::shared_ptr<MeshType> mesh )
6401 {
6402  return FunctionSpace<MeshType,bases<Lagrange<Order+1,Vectorial>,Lagrange<Order,Scalar>>>::New( mesh );
6403 }
6404 
6405 
6406 } // Feel
6407 
6408 
6409 
6410 #if 0
6411 template<
6412 typename A0,
6413  typename A1,
6414  typename A2,
6415  typename A3,
6416  typename A4,
6417  typename T,
6418  typename Cont>
6419 struct FSElement: public Feel::FunctionSpace<A0,A1,A2,A3,A4>::template Element<T, Cont>
6420 {
6421 };
6422 
6423 template<
6424 typename A0,
6425  typename A1,
6426  typename A2,
6427  typename A3,
6428  typename A4>
6429 //struct version< typename Feel::FunctionSpace<A0,A1,A2,A3,A4>::template Element<double,Feel::VectorUblas<double> > >
6430 struct version< typename Feel::FunctionSpace<A0,A1,A2,A3,A4>::element_type >
6431 {
6432  //typedef typename version< typename Feel::FunctionSpace<A0,A1,A2,A3,A4>::template Element<double,Feel::VectorUblas<double> > > version_type;
6433  typedef typename version< typename Feel::FunctionSpace<A0,A1,A2,A3,A4>::element_type > version_type;
6434  typedef mpl::int_<2> type;
6435  typedef mpl::integral_c_tag tag;
6436  BOOST_STATIC_CONSTANT( unsigned int, value = version_type::type::value );
6437 };
6438 
6439 #define FEELPP_REGISTER_ELEMENT( element_type ) \
6440  namespace boost { \
6441  namespace serialization { \
6442  template<> \
6443  struct version<element_type> \
6444  { \
6445  typedef mpl::int_<2> type; \
6446  typedef mpl::integral_c_tag tag; \
6447  BOOST_STATIC_CONSTANT(unsigned int, value = version::type::value); \
6448  }; \
6449  } \
6450  }
6451 #
6452 
6453 #endif
6454 
6455 
6456 
6457 
6458 #endif /* __FunctionSpace_H */

Generated on Sun Oct 20 2013 08:24:59 for Feel++ by doxygen 1.8.4