30 #ifndef __FunctionSpace_H
31 #define __FunctionSpace_H 1
33 #include <boost/static_assert.hpp>
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>
45 #include <boost/serialization/vector.hpp>
46 #include <boost/serialization/array.hpp>
47 #include <boost/serialization/version.hpp>
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>
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>
60 #include <boost/smart_ptr/enable_shared_from_this.hpp>
70 #include <feel/feelcore/parameter.hpp>
77 #include <feel/feelmesh/regiontree.hpp>
78 #include <feel/feelpoly/geomap.hpp>
86 #include <feel/feeldiscr/parameter.hpp>
95 namespace fusion = boost::fusion;
96 namespace parameter = boost::parameter;
104 std::vector<T> operator()( std::vector<T>
const& v1, std::vector<T>
const& v2 )
const
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() );
109 for (
size_type i = 0; i < v1.size(); ++i )
115 template<
typename T,
int M,
int N>
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;
127 typedef typename array_type::template array_view<1>::type type;
139 ID( array_type
const&
id )
145 template<
typename Elem,
typename ContextType>
146 ID( Elem
const& elem, ContextType
const & context )
148 M_id( elem.idExtents( context ) )
150 elem.id_( context, M_id );
165 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q )
const
167 return M_id[q]( c1,c2 );
170 template <
typename ExtentList>
171 void resize(
const ExtentList& extents )
173 M_id.resize( extents );
178 template<
class Archive>
179 void save( Archive & ar,
const unsigned int )
const
182 DVLOG(2) <<
"saving in archive e1= " << e1 <<
"\n";
185 DVLOG(2) <<
"saving in archive e2= " << e2 <<
"\n";
188 DVLOG(2) <<
"saving in archive e3= " << e3 <<
"\n";
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";
194 template<
class Archive>
195 void load( Archive & ar,
const unsigned int )
199 DVLOG(2) <<
"loading from archive e1= " << e1 <<
"\n";
201 DVLOG(2) <<
"loading from archive e2= " << e2 <<
"\n";
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";
210 BOOST_SERIALIZATION_SPLIT_MEMBER()
212 template<
typename T,
int M,
int N>
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;
223 typedef array_type type;
231 template<
typename Elem,
typename ContextType>
232 DD( Elem
const& elem, ContextType
const & context )
234 M_grad( elem.gradExtents( context ) )
236 elem.grad_( context, M_grad );
239 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q )
const
241 return M_grad[q]( c1,c2 );
245 template<
class Archive>
246 void save( Archive & ar,
const unsigned int )
const
249 DVLOG(2) <<
"saving in archive e1= " << e1 <<
"\n";
252 DVLOG(2) <<
"saving in archive e2= " << e2 <<
"\n";
255 DVLOG(2) <<
"saving in archive e3= " << e3 <<
"\n";
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";
261 template<
class Archive>
262 void load( Archive & ar,
const unsigned int )
266 DVLOG(2) <<
"loading from archive e1= " << e1 <<
"\n";
268 DVLOG(2) <<
"loading from archive e2= " << e2 <<
"\n";
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";
277 BOOST_SERIALIZATION_SPLIT_MEMBER()
280 template<
typename T,
int N,
int M,
int P>
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;
291 typedef typename array_type::template array_view<1>::type type;
299 template<
typename Elem,
typename ContextType>
300 D( Elem
const& elem, ContextType
const & context )
302 M_grad( elem.dExtents( context ) )
304 elem.d_( N, context, M_grad );
307 value_type operator()( uint16_type c1, uint16_type , uint16_type q )
const
309 return M_grad[q]( c1,0 );
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;
324 typedef typename array_type::template array_view<1>::type type;
332 template<
typename Elem,
typename ContextType>
333 Div( Elem
const& elem, ContextType
const & context )
335 M_div( elem.divExtents( context ) )
337 elem.div_( context, M_div );
339 uint16_type nComponents1 = elem.nComponents1;
340 std::fill( M_div.data(), M_div.data()+M_div.num_elements(), value_type( 0 ) );
342 M_grad( elem.div_( context, pc, M_div ) ),
345 const uint16_type nq = context.xRefs().size2();
347 for (
int c1 = 0; c1 < nComponents1; ++c1 )
348 for ( uint16_type q = 0; q < nq ; ++q )
350 M_div[q]( 0,0 ) += M_grad[q]( c1,c1 );
356 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q )
const
358 return M_div[q]( c1,c2 );
362 template<
typename T,
int N,
int D>
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;
372 typedef typename array_type::template array_view<1>::type type;
380 template<
typename Elem,
typename ContextType>
381 Curl( Elem
const& elem, ContextType
const & context )
383 M_curl( elem.curlExtents( context ) )
385 init( elem, context, boost::is_same<mpl::int_<N>, mpl::int_<-1> >() );
387 template<
typename Elem,
typename ContextType>
389 init( Elem
const& elem, ContextType
const& context, mpl::bool_<true> )
391 elem.curl_( context, M_curl );
393 template<
typename Elem,
typename ContextType>
395 init( Elem
const& elem, ContextType
const& context, mpl::bool_<false> )
397 elem.curl_( context, M_curl, N );
400 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q )
const
402 return this->operator()( c1, c2, q, mpl::int_<N>() );
404 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<-1> )
const
406 return M_curl[q]( c1,c2 );
408 value_type operator()( uint16_type , uint16_type , uint16_type q, mpl::int_<0> )
const
410 return M_curl[q]( 0,0 );
412 value_type operator()( uint16_type , uint16_type , uint16_type q, mpl::int_<1> )
const
414 return M_curl[q]( 0,0 );
416 value_type operator()( uint16_type , uint16_type , uint16_type q, mpl::int_<2> )
const
418 return M_curl[q]( 0,0 );
423 template<
typename T,
int M,
int N>
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;
435 typedef typename array_type::template array_view<1>::type type;
443 template<
typename Elem,
typename ContextType>
444 H( Elem
const& elem, ContextType
const & context )
446 M_hess( elem.hessExtents( context ) )
448 elem.hess_( context, M_hess );
451 value_type operator()( uint16_type c1, uint16_type c2, uint16_type q )
const
453 return M_hess[q]( c1,c2 );
458 template<
typename MeshPtrType,
typename PeriodicityType = Periodicity<NoPeriodicity> >
459 struct InitializeSpace
461 InitializeSpace( MeshPtrType
const& mesh,
462 PeriodicityType
const& periodicity,
463 std::vector<Dof>
const& dofindices,
464 std::vector<WorldComm>
const & worldsComm )
467 M_worldsComm( worldsComm ),
469 M_dofindices( dofindices ),
470 M_periodicity( periodicity )
473 template <
typename T>
474 void operator()( boost::shared_ptr<T> & x )
const
476 operator()( x, is_shared_ptr<MeshPtrType>() );
478 template <
typename T>
479 void operator()( boost::shared_ptr<T> & x, mpl::bool_<true> )
const
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" );
488 template <
typename T>
489 void operator()( boost::shared_ptr<T> & x, mpl::bool_<false> )
const
491 auto p = *fusion::find<typename T::periodicity_0_type>(M_periodicity);
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" );
501 mutable uint16_type M_cursor;
502 std::vector<WorldComm> M_worldsComm;
504 std::vector<Dof>
const& M_dofindices;
505 PeriodicityType M_periodicity;
507 template<
typename DofType>
508 struct updateDataMapProcess
510 updateDataMapProcess( std::vector<WorldComm>
const & worldsComm,
511 WorldComm
const& worldCommFusion,
512 uint16_type lastCursor )
516 M_lastCursor( lastCursor ),
517 M_worldsComm( worldsComm ),
518 M_dm( new DofType( worldCommFusion ) ),
519 M_dmOnOff( new DofType( worldCommFusion ) )
522 template <
typename T>
523 void operator()( boost::shared_ptr<T> & x )
const
526 if ( M_worldsComm[M_cursor].isActive() )
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() );
537 M_dm->resizeMapGlobalProcessToGlobalCluster( nLocWithGhost );
538 M_dm->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhost );
540 for (
size_type i=0; i<nLocWithGhost; ++i )
541 M_dm->setMapGlobalProcessToGlobalCluster( i, M_start_index + x->dof()->mapGlobalProcessToGlobalCluster( i ) );
543 for (
size_type i=0; i<nLocWithoutGhost; ++i )
544 M_dm->setMapGlobalClusterToGlobalProcess( i, x->dof()->mapGlobalClusterToGlobalProcess( i ) );
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(),
555 M_dmOnOff->setNLocalDofWithGhost( M_dmOnOff->worldComm().globalRank(),
559 if ( M_cursor==M_lastCursor )
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() );
568 size_type nLocWithoutGhostOnOff= M_dmOnOff->nLocalDofWithoutGhost() + x->dof()->nLocalDofWithoutGhost();
569 size_type nLocWithGhostOnOff= M_dmOnOff->nLocalDofWithGhost() + x->dof()->nLocalDofWithGhost();
571 M_dmOnOff->setNLocalDofWithoutGhost( M_dmOnOff->worldComm().globalRank(),
572 nLocWithoutGhostOnOff );
573 M_dmOnOff->setNLocalDofWithGhost( M_dmOnOff->worldComm().globalRank(),
574 nLocWithGhostOnOff );
577 M_dmOnOff->resizeMapGlobalProcessToGlobalCluster( nLocWithGhostOnOff );
578 M_dmOnOff->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhostOnOff );
580 size_type startGlobClusterDof = M_dmOnOff->nLocalDofWithoutGhost() - x->dof()->nLocalDofWithoutGhost();
581 size_type startGlobProcessDof = M_dmOnOff->nLocalDofWithGhost() - x->dof()->nLocalDofWithGhost();
583 for (
size_type i=0; i<x->dof()->nLocalDofWithGhost(); ++i )
585 M_dmOnOff->setMapGlobalProcessToGlobalCluster( startGlobProcessDof + i, M_start_index + x->dof()->mapGlobalProcessToGlobalCluster( i ) );
588 for (
size_type i=0; i<x->dof()->nLocalDofWithoutGhost(); ++i )
590 M_dmOnOff->setMapGlobalClusterToGlobalProcess( startGlobClusterDof + i, x->dof()->mapGlobalClusterToGlobalProcess( i ) );
594 M_start_index+=x->nDof();
599 boost::shared_ptr<DofType> dataMap()
const
603 boost::shared_ptr<DofType> dataMapOnOff()
const
608 mutable uint16_type M_cursor;
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;
618 template<
typename DofType>
619 struct updateDataMapProcessStandard
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)
628 M_start_index( worldCommFusion.globalSize(), 0 ),
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)
638 const int myrank = M_dm->worldComm().globalRank();
639 const int worldsize = M_dm->worldComm().globalSize();
640 for (
int proc = 0 ; proc < worldsize ; ++proc)
642 M_dm->setNLocalDofWithoutGhost( proc, nLocWithoutGhost[proc] );
643 M_dm->setNLocalDofWithGhost( proc, nLocWithGhost[proc] );
644 M_start_index[proc] = M_startDofGlobalCluster[proc];
647 M_dm->resizeMapGlobalProcessToGlobalCluster( nLocWithGhost[myrank] );
648 M_dm->resizeMapGlobalClusterToGlobalProcess( nLocWithoutGhost[myrank] );
651 template <
typename T>
652 void operator()( boost::shared_ptr<T> & x )
const
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);
659 for (
int proc = 0 ; proc < worldsize ; ++proc)
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] );
669 M_dm->setLastDofGlobalCluster( proc, M_startDofGlobalCluster[proc] + x->dof()->nLocalDofWithoutGhost(proc) - 1 );
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) );
678 const size_type firstDofGC = M_dm->firstDofGlobalCluster();
679 const size_type firstDofGCSubSpace = x->dof()->firstDofGlobalCluster();
681 for (
size_type gdof = x->dof()->firstDof(myrank) ; gdof < x->dof()->nLocalDofWithGhost(myrank) ; ++gdof )
683 const size_type localDof = M_startIndexWithGhost+gdof;
684 const size_type gdofGC = x->dof()->mapGlobalProcessToGlobalCluster(gdof);
686 if ( x->dof()->dofGlobalClusterIsOnProc( gdofGC ) )
688 const size_type globalDof = M_start_index[myrank] + ( gdofGC-firstDofGCSubSpace );
689 M_dm->setMapGlobalProcessToGlobalCluster( localDof, globalDof );
690 M_dm->setMapGlobalClusterToGlobalProcess( globalDof-firstDofGC ,localDof );
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 );
700 for (
int proc = 0 ; proc < worldsize ; ++proc)
701 M_start_index[proc] += x->dof()->nLocalDofWithoutGhost(proc);
702 M_startIndexWithGhost+=nLocWithGhost;
707 boost::shared_ptr<DofType> dataMap()
const
711 boost::shared_ptr<DofType> dataMapOnOff()
const
716 mutable uint16_type M_cursor;
717 mutable std::vector<size_type> M_start_index;
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;
738 template<
typename Sig>
741 template<
typename T,
typename S>
742 #if BOOST_VERSION < 104200
743 struct result<NbDof( T,S )>
745 struct result<NbDof( S,T )>
748 boost::remove_reference<S>
750 template <
typename T>
752 operator()( T
const& x,
size_type s )
const
759 if ( M_cursor < M_finish )
766 template <
typename T>
768 operator()(
size_type s, T
const& x )
const
770 return this->operator()( x, s );
785 template<
typename Sig>
788 template<
typename T,
typename S>
789 #if BOOST_VERSION < 104200
790 struct result<NLocalDof( T,S )>
792 struct result<NLocalDof( S,T )>
795 boost::remove_reference<S>
797 template <
typename T>
799 operator()( T
const& x,
size_type s )
const
803 if ( M_cursor < M_finish )
804 ret += x->nLocalDof();
809 template <
typename T>
811 operator()(
size_type s, T
const& x )
const
813 return this->operator()( x, s );
820 template<
typename IsWithGhostType>
824 NLocalDof( std::vector<WorldComm>
const & worldsComm = std::vector<WorldComm>( 1,
Environment::worldComm() ),
825 bool useOffSubSpace =
false,
830 M_worldsComm( worldsComm ),
831 M_useOffSubSpace( useOffSubSpace )
833 template<
typename Sig>
836 template<
typename T,
typename S>
837 #if BOOST_VERSION < 104200
838 struct result<NLocalDof( T,S )>
840 struct result<NLocalDof( S,T )>
843 boost::remove_reference<S>
846 template <
typename T>
848 nLocalDof( T
const& x, mpl::bool_<true> )
const
850 return x->nLocalDofWithGhost();
853 template <
typename T>
855 nLocalDof( T
const& x, mpl::bool_<false> )
const
857 return x->nLocalDofWithoutGhost();
860 template <
typename T>
862 operator()( T
const& x,
size_type s )
const
866 if ( M_cursor < M_finish )
868 if ( M_useOffSubSpace )
870 ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
875 if ( M_worldsComm[M_cursor].isActive() )
876 ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
884 template <
typename T>
886 operator()(
size_type s, T
const& x )
const
888 return this->operator()( x, s );
893 std::vector<WorldComm>
const& M_worldsComm;
894 bool M_useOffSubSpace;
899 template<
typename IsWithGhostType>
900 struct NLocalDofOnProc
903 NLocalDofOnProc(
const int proc,
905 bool useOffSubSpace =
false,
911 M_worldsComm( worldsComm ),
912 M_useOffSubSpace( useOffSubSpace )
915 template<
typename Sig>
918 template<
typename T,
typename S>
919 #if BOOST_VERSION < 104200
920 struct result<NLocalDofOnProc( T,S )>
922 struct result<NLocalDofOnProc( S,T )>
925 boost::remove_reference<S>
928 template <
typename T>
930 nLocalDof( T
const& x, mpl::bool_<true> )
const
932 return x->nLocalDofWithGhostOnProc(M_proc);
935 template <
typename T>
937 nLocalDof( T
const& x, mpl::bool_<false> )
const
939 return x->nLocalDofWithoutGhostOnProc(M_proc);
942 template <
typename T>
944 operator()( T
const& x,
size_type s )
const
948 if ( M_cursor < M_finish )
950 if ( M_useOffSubSpace )
952 ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
957 if ( M_worldsComm[M_cursor].isActive() )
958 ret += nLocalDof( x, mpl::bool_<IsWithGhostType::value>() );
966 template <
typename T>
968 operator()(
size_type s, T
const& x )
const
970 return this->operator()( x, s );
976 std::vector<WorldComm>
const& M_worldsComm;
977 bool M_useOffSubSpace;
981 template<
int i,
typename SpaceCompositeType>
982 struct InitializeContainersOff
984 InitializeContainersOff( boost::shared_ptr<SpaceCompositeType>
const& _space )
989 template <
typename T>
990 void operator()( boost::shared_ptr<T> & x )
const
992 if ( M_cursor==i && !x )
993 x = boost::shared_ptr<T>(
new T( M_space->template functionSpace<i>()->dof() ) );
997 mutable uint16_type M_cursor;
998 boost::shared_ptr<SpaceCompositeType> M_space;
1002 template<
int i,
typename SpaceCompositeType>
1003 struct SendContainersOn
1005 SendContainersOn( boost::shared_ptr<SpaceCompositeType>
const& _space,
1006 std::vector<double>
const& _dataToSend )
1010 M_dataToSend( _dataToSend )
1012 template <
typename T>
1013 void operator()( boost::shared_ptr<T> & x )
const
1017 int locRank=M_space->worldComm().localRank();
1018 int globRank=M_space->worldComm().localColorToGlobalRank( M_cursor,locRank );
1022 M_space->worldComm().globalComm().send( globRank,tag,M_dataToSend );
1027 mutable uint16_type M_cursor;
1028 boost::shared_ptr<SpaceCompositeType> M_space;
1029 std::vector<double> M_dataToSend;
1033 template<
int i,
typename SpaceCompositeType>
1034 struct RecvContainersOff
1036 RecvContainersOff( boost::shared_ptr<SpaceCompositeType>
const& _space )
1041 template <
typename T>
1042 void operator()( boost::shared_ptr<T> & x )
const
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 );
1052 M_space->worldComm().globalComm().recv( globRank,tag,dataToRecv );
1053 std::copy( dataToRecv.begin(), dataToRecv.end(), x->begin() );
1058 mutable uint16_type M_cursor;
1059 boost::shared_ptr<SpaceCompositeType> M_space;
1064 template<
typename map_type >
1065 struct searchIndicesBySpace
1067 searchIndicesBySpace()
1070 searchIndicesBySpace( map_type& )
1073 template<
typename T>
1074 searchIndicesBySpace( T
const& fspace, map_type& u )
1076 u = getIndicesFromSpace( fspace,u );
1079 template<
typename Sig>
1082 template<
typename T,
typename M>
1083 #if BOOST_VERSION < 104200
1084 struct result<searchIndicesBySpace( T,M )>
1086 struct result<searchIndicesBySpace( M,T )>
1089 boost::remove_reference<M>
1092 template <
typename T >
1093 map_type getIndicesFromSpace( T
const& fspace, map_type t )
const
1095 if ( fspace->mesh()->numElements() == 0 )
1098 size_type nProc = fspace->dof()->nProcessors();
1101 std::vector< size_type > max_per_space;
1108 max_per_space.push_back( t[j][_end-1] );
1115 max_index = *max_element( max_per_space.begin(), max_per_space.end() ) + 1;
1128 size_type _first = fspace->dof()->firstDof( i );
1129 size_type _last = fspace->dof()->lastDof( i );
1132 for (
size_type j=_first; j<=_last; j++ )
1133 t[i].push_back( max_index + j );
1138 template <
typename T>
1140 #if BOOST_VERSION < 104200
1141 operator()( T
const& fspace, map_type t )
const
1143 operator()( map_type t, T
const& fspace )
const
1146 return getIndicesFromSpace( fspace,t );
1151 struct computeStartOfFieldSplit
1153 typedef boost::tuple< uint16_type , size_type > result_type;
1155 template<
typename T>
1156 result_type operator()( result_type
const & previousRes, T
const& t )
1158 auto cptSpaces = previousRes.get<0>();
1159 auto start = previousRes.get<1>();
1161 for (
int proc=0;proc<t->dof()->worldComm().globalSize();++proc)
1163 if (proc < t->dof()->worldComm().globalRank())
1164 start+=t->dof()->nLocalDofWithoutGhost(proc);
1166 return boost::make_tuple( ++cptSpaces, start );
1172 struct computeNDofForEachSpace
1175 computeNDofForEachSpace(
size_type startSplit)
1177 M_startSplit(startSplit)
1180 typedef boost::tuple< uint16_type, size_type, std::vector<std::vector<size_type> > > result_type;
1182 template<
typename T>
1183 result_type operator()( result_type
const & previousRes, T
const& t )
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>();
1195 is[cptSpaces].resize(nDofWithoutGhost);
1197 for (
size_type i=0; i<nDofWithGhost; ++i )
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);
1204 return boost::make_tuple( ++cptSpaces, ( start+nDofWithoutGhost ), is );
1210 struct rebuildDofPointsTool
1213 template <
typename T>
1214 void operator()( boost::shared_ptr<T> & x )
const
1216 x->dof()->rebuildDofPoints( *x->mesh() );
1222 typedef std::string result_type;
1224 template<
typename T>
1225 result_type operator()( result_type
const & previousRes, T
const& t )
1227 std::ostringstream os;
1229 if ( previousRes.size() )
1230 os << previousRes <<
"_" << t->basis()->familyName();
1233 os << t->basis()->familyName();
1241 typedef std::vector<int> result_type;
1243 template<
typename T>
1244 result_type operator()( result_type
const & previousRes, T
const& t )
1246 std::vector<int> res( previousRes );
1247 res.push_back( t->nSubFunctionSpace() );
1252 template<
typename ElementType>
1253 struct CreateElementVector
1256 template<
typename Sig>
1259 template<
typename Lhs,
typename Rhs>
1260 struct result<CreateElementVector( Lhs,Rhs )>
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;
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;
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;
1278 template<
typename Lhs,
typename Rhs>
1279 typename result<CreateElementVector( Lhs,Rhs )>::type
1280 operator()( Lhs
const& lhs, Rhs
const& rhs )
const
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 )
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] );
1296 else if ( ( M_names.size() == 1 ) && s > 1 )
1298 elt.setName( (boost::format(
"%1%-%2%" ) % M_names[0] % index::value ).str() );
1302 elt.setName( (boost::format(
"%1%-%2%" ) % M_e.name() % index::value ).str() );
1304 return boost::fusion::as_vector( boost::fusion::push_back( lhs, elt ) );
1323 enum FunctionSpaceType
1329 template<uint16_type PN,
1333 static const uint16_type PolynomialOrder = PN;
1334 static const uint16_type GeometricOrder = GN;
1336 static const bool is_isoparametric = ( PN == GN );
1337 static const bool is_subparametric = ( PN > GN );
1338 static const bool is_surparametric = ( PN < GN );
1341 typedef parameter::parameters<
1343 parameter::required<tag::mesh_type, boost::is_base_and_derived<MeshBase,_> >
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<_> > > >
1349 , parameter::optional<parameter::deduced<tag::bases_list>, boost::is_base_and_derived<detail::bases_base,_> >
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;
1370 typename A1 = parameter::void_,
1371 typename A2 = parameter::void_,
1372 typename A3 = parameter::void_,
1373 typename A4 = parameter::void_>
1377 public boost::enable_shared_from_this<FunctionSpace<A0,A1,A2,A3,A4> >
1380 typedef typename functionspace_signature::bind<A0,A1,A2,A3,A4>::type args;
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;
1387 BOOST_MPL_ASSERT_NOT( ( boost::is_same<mpl::at<bases_list,mpl::int_<0> >, mpl::void_> ) );
1391 template<
typename ThePeriodicityType,
int pos>
1392 struct GetPeriodicity
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;
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;
1406 template<
typename BasisType>
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;
1416 detail::bases<BasisType>,value_type,
1417 Periodicity<typename GetPeriodicity<periodicity_type,pos::value>::type > > _type;
1418 typedef boost::shared_ptr<_type> type;
1422 template<
typename BasisType>
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;
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;
1436 typedef typename mpl::transform<bases_list, ChangeBasis<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type functionspace_vector_type;
1438 template<
typename BasisType>
1439 struct ChangeBasisToComponentBasis
1441 typedef typename BasisType::component_basis_type component_basis_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;
1448 typedef typename mpl::transform<bases_list,
1449 ChangeBasisToComponentBasis<mpl::_1>,
1450 mpl::back_inserter<fusion::vector<> > >::type component_functionspace_vector_type;
1452 template<
typename BasisType>
1453 struct GetComponentBasis
1455 typedef typename BasisType::component_basis_type type;
1458 typedef detail::bases2<
typename mpl::transform<bases_list,
1459 GetComponentBasis<mpl::_1>,
1460 mpl::back_inserter<fusion::vector<> > >::type> component_basis_vector_type;
1469 static const bool is_composite = ( mpl::size<bases_list>::type::value > 1 );
1471 template<
typename MeshListType,
int N>
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;
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;
1486 template<
typename MeshType>
1487 struct ChangeToMeshPtr
1489 typedef boost::shared_ptr<MeshType> type;
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;
1499 template<
typename BasisType>
1500 struct GetNComponents
1503 typedef mpl::int_<BasisType::template apply<mesh_type::nDim,
1504 mesh_type::nRealDim,
1506 typename mesh_type::element_type>::type::nComponents> type;
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;
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,
1522 typename GetMesh<meshes_list,0>::type::element_type>::type basis_0_type;
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 );
1533 typedef typename basis_0_type::continuity_type continuity_type;
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;
1541 typedef typename GetPeriodicity<periodicity_type,0>::type periodicity_0_type;
1542 static const bool is_periodic = periodicity_0_type::is_periodic;
1550 typedef typename ublas::type_traits<value_type>::real_type real_type;
1555 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
1556 typedef boost::shared_ptr<functionspace_type> pointer_type;
1559 typedef boost::shared_ptr<component_functionspace_type> component_functionspace_ptrtype;
1563 typedef bases_list BasisType;
1567 typedef typename mpl::at_c<bases_list,N>::type type;
1569 typedef typename mpl::if_<mpl::bool_<is_composite>,
1570 mpl::identity<bases_list>,
1571 mpl::identity<basis_0_type> >::type::type basis_type;
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;
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;
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,
1595 convex_type> >::type::type component_basis_type;
1597 typedef boost::shared_ptr<component_basis_type> component_basis_ptrtype;
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;
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;
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;
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;
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;
1646 typedef boost::shared_ptr<dof_type> dof_ptrtype;
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;
1655 typedef boost::function<return_type ( node_type const& )> function_type;
1659 struct sub_functionspace
1661 typedef typename mpl::at_c<functionspace_vector_type,i>::type type;
1673 template<
typename T =
double,
typename Cont = VectorUblas<T> >
1676 public Cont,boost::addable<Element<T,Cont> >, boost::subtractable<Element<T,Cont> >
1679 typedef T value_type;
1681 template<
typename BasisType>
1682 struct ChangeElement
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;
1690 typedef typename mpl::transform<bases_list, ChangeElement<mpl::_1>, mpl::back_inserter<mpl::vector<> > >::type element_vector_type;
1699 template<
typename BasisType>
1702 BOOST_MPL_ASSERT_NOT( ( boost::is_same<BasisType,mpl::void_> ) );
1703 typedef boost::shared_ptr<Cont> type;
1705 typedef typename mpl::transform<bases_list, AddOffContainer<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type container_vector_type;
1710 typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
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;
1729 typedef typename ublas::type_traits<value_type>::real_type real_type;
1730 typedef T element_type;
1733 typedef Cont container_type;
1734 typedef container_type vector_temporary_type;
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;
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;
1745 typedef typename functionspace_type::return_type return_type;
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;
1753 typedef Element<T,Cont> this_type;
1757 typedef typename mpl::at_c<element_vector_type,i>::type type;
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;
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;
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;
1782 Element( Element
const& __e );
1786 Element( functionspace_ptrtype
const& __functionspace,
1787 std::string
const& __name =
"unknown",
1789 ComponentType __ct = NO_COMPONENT );
1791 Element( functionspace_ptrtype
const& __functionspace,
1792 container_type
const& __c,
1793 std::string
const& __name =
"unknown",
1795 ComponentType __ct = NO_COMPONENT );
1799 void initFromSpace( functionspace_ptrtype
const& __functionspace,
1800 container_type
const& __c );
1808 Element&
operator=( Element
const& __e );
1810 template<
typename ExprT>
1811 Element&
operator=( vf::Expr<ExprT>
const& __expr )
1813 *
this =
project( this->functionspace(),
elements( this->functionspace()->
mesh() ), __expr );
1817 template<
typename VectorExpr>
1818 Element&
operator=( VectorExpr
const& __v );
1824 super
const& container()
const
1843 template<ComponentType THECOMP>
1848 return this->
template comp<THECOMP>(
typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1850 template<ComponentType THECOMP>
1852 comp( mpl::bool_<true> )
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 );
1858 typename component_type::container_type( this->vec().data().expression(), s ),
1863 template<ComponentType THECOMP>
1865 comp( mpl::bool_<false> )
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 );
1871 typename component_type::container_type( ( VectorUblas<value_type>& )*
this, s ),
1885 comp( ComponentType i )
const
1888 return comp( i,
typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1891 comp( ComponentType i, mpl::bool_<true> )
const
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 );
1898 typename component_type::container_type( this->vec().data().expression(), s ),
1905 comp( ComponentType i, mpl::bool_<false> )
const
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 );
1912 typename component_type::container_type( ( VectorUblas<value_type>& )*
this, s ),
1927 comp( ComponentType i )
1930 return comp( i,
typename mpl::not_<boost::is_same<container_type,VectorUblas<value_type> > >::type() );
1933 comp( ComponentType i, mpl::bool_<true> )
1935 FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1936 auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1938 std::string __name = this->name() +
"_" + componentToString( i );
1940 typename component_type::container_type( this->vec().data().expression(), s ),
1947 comp( ComponentType i, mpl::bool_<false> )
1949 FEELPP_ASSERT( i >= X && i < N_COMPONENTS );
1950 auto s = ublas::slice( i, N_COMPONENTS, M_functionspace->nDofPerComponent() );
1952 std::string __name = this->name() +
"_" + componentToString( i );
1954 typename component_type::container_type( ( VectorUblas<value_type>& )*
this, s ),
1963 size_type index=start()+boost::get<0>( M_functionspace->dof()->localToGlobal( e, l, c ) );
1964 return super::operator()( index );
1969 size_type index=start()+boost::get<0>( M_functionspace->dof()->localToGlobal( e, l, c ) );
1970 return super::operator()( index );
1975 return super::operator()( i );
1979 return super::operator()( i );
1981 Element& operator+=( Element
const& _e )
1983 for (
int i=0; i < _e.nLocalDof(); ++i )
1984 this->
operator()( i ) += _e( i );
1988 Element& operator-=( Element
const& _e )
1990 for (
size_type i=0; i < _e.nLocalDof(); ++i )
1991 this->
operator()( i ) -= _e( i );
2000 DVLOG(2) <<
"Update element after a change in the mesh\n";
2003 template<
typename AE>
2004 container_type& assign(
const ublas::vector_expression<AE> &ae )
2006 return super::assign( ae );
2008 void assign(
size_type ie, uint16_type il, uint16_type c, value_type
const& __v )
2010 size_type index=start()+ M_functionspace->dof()->localToGlobal( ie, il, c ).index();
2011 this->operator[]( index ) = __v;
2013 void plus_assign(
size_type ie, uint16_type il, uint16_type c, value_type
const& __v )
2015 size_type index=start()+ M_functionspace->dof()->localToGlobal( ie, il, c ).index();
2016 this->operator[]( index ) += __v;
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;
2041 DataMap
const&
map()
const
2043 return M_functionspace->map();
2051 return M_functionspace->mesh();
2057 mesh_ptrtype
mesh()
const
2059 return M_functionspace->mesh();
2065 Element<T,Cont> sqrt()
const
2067 Element<T,Cont> _e( M_functionspace );
2069 for (
int i=0; i < _e.nLocalDof(); ++i ) _e( i ) = math::sqrt( this->
operator()( i ) );
2077 Element<T,Cont> pow(
int n )
const
2079 Element<T,Cont> _e( M_functionspace );
2081 for (
int i=0; i < _e.nLocalDof(); ++i ) _e( i ) = math::pow( this->
operator()( i ),n );
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 )
2091 FEELPP_ASSERT( P0h->mesh() == this->
mesh() ).error(
"mesh is not the same" );
2092 FEELPP_ASSERT( is_scalar ).error(
"only works for scalar fields" );
2094 typename p0_space_type::element_type p0Element( P0h );
2096 for (
auto elt_it = P0h->mesh()->beginElement(); elt_it != P0h->mesh()->endElement(); ++elt_it )
2100 size_type dofp0 = boost::get<0>( P0h->dof()->localToGlobal( eid, 0, 0 ) );
2101 std::vector<value_type> values ( functionspace_type::fe_type::nLocalDof );
2105 for ( uint16_type local_id=0; local_id < functionspace_type::fe_type::nLocalDof; ++local_id )
2107 dofpn = boost::get<0>( this->
functionSpace()->dof()->localToGlobal( eid, local_id, 0 ) );
2108 values[local_id] = this->
operator()( dofpn );
2111 if ( extreme ==
"max" )
2112 p0Element.assign( eid, 0, 0, *( std::max_element( values.begin(), values.end() ) ) );
2115 p0Element.assign( eid, 0, 0, *( std::min_element( values.begin(), values.end() ) ) );
2121 value_type max()
const
2123 return super::max();
2126 template <
typename p0_space_type >
2127 typename p0_space_type::element_type max( boost::shared_ptr<p0_space_type>
const& P0h )
2129 return this->extremeValue( P0h,
"max" );
2132 value_type min()
const
2134 return super::min();
2137 template <
typename p0_space_type >
2138 typename p0_space_type::element_type min( boost::shared_ptr<p0_space_type>
const& P0h )
2140 return this->extremeValue( P0h,
"min" );
2149 typedef detail::ID<value_type,nComponents1,nComponents2> id_type;
2156 template<
typename ContextType>
2157 boost::array<typename array_type::index, 1>
2158 idExtents( ContextType
const & context )
const
2160 boost::array<typename array_type::index, 1> shape;
2161 shape[0] = context.xRefs().size2();
2167 template<
typename Context_t>
2169 id( Context_t
const & context )
const
2171 return id_type( *
this, context );
2185 template<
typename Context_t>
2187 id_( Context_t
const & context, id_array_type& v )
const;
2189 template<
typename Context_t>
2191 id( Context_t
const & context, id_array_type& v )
const
2198 idInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const;
2208 template<
typename Context_t>
2210 ptsInContext( Context_t
const & context, mpl::int_<1> )
const
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;
2216 gmc_interp_ptrtype __c_interp(
new gmc_interp_type( context.geometricMapping(), context.element_c(), context.pc() ) );
2218 return __c_interp->xReal();
2226 template<
typename Context_t>
2228 ptsInContext( Context_t
const & context, mpl::int_<2> )
const
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;
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;
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() ) );
2243 return __c_interp->xReal();
2250 gmc_ptrtype geomapPtr( geoelement_type
const& elt )
const
2253 geopc_ptrtype __geopc(
new geopc_type( __gm, elt.G() ) );
2254 return gmc_ptrtype(
new gmc_type( __gm, elt, __geopc ) );
2260 pc_ptrtype pcPtr( geoelement_type
const& elt )
const
2265 id_type
operator()( Eigen::Matrix<value_type,nDim,1>
const& __x,
bool extrapolate =
false )
const
2267 node_type n( nDim );
2268 for(
int i = 0; i < nDim; ++i ) n[i]=__x[i];
2277 id_type
operator()( node_type
const& __x,
bool extrapolate =
false )
const;
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;
2289 template<
typename ContextType>
2290 boost::array<typename array_type::index, 1>
2291 gradExtents( ContextType
const & context )
const
2293 boost::array<typename array_type::index, 1> shape;
2294 shape[0] = context.xRefs().size2();
2300 template<
typename ContextType>
2302 grad( ContextType
const & context )
const
2304 return grad_type( *
this, context );
2306 template<
typename ElementType>
2308 grad( ElementType
const& elt, grad_array_type& v )
2310 gmc_ptrtype gmc( geomapPtr( elt ) );
2311 v.resize( gradExtents( *gmc ) );
2312 grad_( *gmc, *pcPtr( elt ), v );
2329 template<
typename ContextType>
2330 void grad_( ContextType
const & context, grad_array_type& v )
const;
2332 template<
typename ContextType>
2333 void grad( ContextType
const & context, grad_array_type& v )
const
2335 grad_( context, v );
2339 gradInterpolate( matrix_node_type __ptsReal, grad_array_type& v )
const;
2346 grad_type grad( node_type
const& __x )
const;
2348 template<
typename ContextType>
2350 dx( ContextType
const & context )
const
2352 return dx_type( *
this, context );
2355 template<
typename ContextType>
2357 dy( ContextType
const & context )
const
2359 return dy_type( *
this, context );
2361 template<
typename ContextType>
2363 dz( ContextType
const & context )
const
2365 return dz_type( *
this, context );
2368 template<
typename ContextType>
2369 boost::array<typename array_type::index, 1>
2370 dxExtents( ContextType
const & context )
const
2372 boost::array<typename array_type::index, 1> shape;
2373 shape[0] = context.xRefs().size2();
2378 template<
typename ContextType>
2379 boost::array<typename array_type::index, 1>
2380 dyExtents( ContextType
const & context )
const
2382 return dxExtents( context );
2384 template<
typename ContextType>
2385 boost::array<typename array_type::index, 1>
2386 dzExtents( ContextType
const & context )
const
2388 return dxExtents( context );
2390 template<
typename ContextType>
2391 boost::array<typename array_type::index, 1>
2392 dExtents( ContextType
const & context )
const
2394 return dxExtents( context );
2397 template<
typename ContextType>
2398 void d_(
int N, ContextType
const & context, id_array_type& v )
const;
2400 template<
typename ContextType>
2401 void dx( ContextType
const & context, id_array_type& v )
const
2403 d_( 0, context, v );
2405 template<
typename ContextType>
2406 void dy( ContextType
const & context, id_array_type& v )
const
2408 d_( 1, context, v );
2410 template<
typename ContextType>
2411 void dz( ContextType
const & context, id_array_type& v )
const
2413 d_( 2, context, v );
2417 dxInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const;
2420 dyInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const;
2423 dzInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const;
2428 typedef detail::Div<value_type> div_type;
2430 template<
typename ContextType>
2431 boost::array<typename array_type::index, 1>
2432 divExtents( ContextType
const & context )
const
2434 boost::array<typename array_type::index, 1> shape;
2435 shape[0] = context.xRefs().size2();
2440 template<
typename ContextType>
2442 div( ContextType
const & context )
const
2445 return div_type( *
this, context );
2449 template<
typename ContextType>
2451 div( ContextType
const & context, div_array_type& v )
const
2456 template<
typename ContextType>
2457 void div_( ContextType
const & context, div_array_type& v )
const;
2460 divInterpolate( matrix_node_type __ptsReal, div_array_type& v )
const;
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;
2467 template<
typename ContextType>
2469 curl( ContextType
const & context )
const
2471 return curl_type( *
this, context );
2474 template<
typename ContextType>
2476 curlx( ContextType
const & context )
const
2478 return curlx_type( *
this, context );
2481 template<
typename ContextType>
2483 curly( ContextType
const & context )
const
2485 return curly_type( *
this, context );
2488 template<
typename ContextType>
2490 curlz( ContextType
const & context )
const
2492 return curlz_type( *
this, context );
2495 template<
typename ContextType>
2496 boost::array<typename array_type::index, 1>
2497 curlExtents( ContextType
const & context )
const
2501 boost::array<typename array_type::index, 1> shape;
2502 shape[0] = context.xRefs().size2();
2507 template<
typename ContextType>
2509 curl( ContextType
const & context, curl_array_type& v )
const
2513 curl_( context, v );
2517 curlInterpolate( matrix_node_type __ptsReal, curl_array_type& v )
const;
2519 template<
typename ContextType>
2520 void curl_( ContextType
const & context, curl_array_type& v )
const;
2522 template<
typename ContextType>
2523 void curl_( ContextType
const & context, comp_curl_array_type& v,
int comp )
const;
2526 template<
typename ContextType>
2527 boost::array<typename array_type::index, 1>
2528 curlxExtents( ContextType
const & context )
const
2530 boost::array<typename array_type::index, 1> shape;
2531 shape[0] = context.xRefs().size2();
2536 template<
typename ContextType>
2537 boost::array<typename array_type::index, 1>
2538 curlyExtents( ContextType
const & context )
const
2540 return curlxExtents( context );
2542 template<
typename ContextType>
2543 boost::array<typename array_type::index, 1>
2544 curlzExtents( ContextType
const & context )
const
2546 return curlxExtents( context );
2549 template<
typename ContextType>
2551 curlx( ContextType
const & context, comp_curl_array_type& v )
const
2554 curl_( context, v, 0 );
2556 template<
typename ContextType>
2558 curly( ContextType
const & context, comp_curl_array_type& v )
const
2561 curl_( context, v, 1 );
2563 template<
typename ContextType>
2565 curlz( ContextType
const & context, comp_curl_array_type& v )
const
2568 curl_( context, v, 2 );
2572 curlxInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const;
2575 curlyInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const;
2578 curlzInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const;
2580 template<
typename ContextType>
2581 boost::array<typename array_type::index, 1>
2582 hessExtents( ContextType
const & context )
const
2584 BOOST_STATIC_ASSERT( rank == 0 );
2586 boost::array<typename array_type::index, 1> shape;
2587 shape[0] = context.xRefs().size2();
2605 template<
typename ContextType>
2607 hess_( ContextType
const & context, hess_array_type& v )
const;
2609 template<
typename ContextType>
2611 hess_( ContextType
const & context, hess_array_type& v, mpl::int_<0> )
const;
2614 hessInterpolate( matrix_node_type __ptsReal, hess_array_type& v )
const;
2616 typedef detail::H<value_type,nRealDim,nRealDim> hess_type;
2618 template<
typename ContextType>
2620 hess( ContextType
const & context )
const
2622 return hess_type( *
this, context );
2625 template<
typename ContextType>
2627 hess( ContextType
const & context, hess_array_type& v )
const
2629 hess_( context, v );
2631 template<
typename ElementType>
2633 hess( ElementType
const& elt, hess_array_type& v )
2635 gmc_ptrtype gmc( geomapPtr( elt ) );
2636 v.resize( hessExtents( *gmc ) );
2637 hess_( *gmc, *pcPtr( elt ), v );
2645 return M_functionspace;
2654 functionspace_vector_type
const&
2657 return M_functionspace->functionSpaces();
2664 typename mpl::at_c<functionspace_vector_type,i>::type
2667 return M_functionspace->template functionSpace<i>();
2670 Eigen::Matrix<value_type, Eigen::Dynamic, 1>
2671 extractValuesWithMarker( std::string
const& m )
2674 Eigen::Matrix<value_type, Eigen::Dynamic, 1> res( std::distance( r.first, r.second ) );
2676 for(
auto it = r.first, en = r.second; it != en; ++it, ++i )
2677 res( i ) = this->operator[]( it->second );
2680 Eigen::Matrix<value_type, Eigen::Dynamic, 1>
2681 extractValuesWithoutMarker( std::string
const& 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 );
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 );
2696 typename mpl::at_c<element_vector_type,i>::type
2697 element( std::string
const& name =
"u",
bool updateOffViews=
true )
2699 if ( this->worldComm().globalSize()>1 ) this->worldComm().globalComm().barrier();
2703 detail::NLocalDof<mpl::bool_<true> >( this->worldsComm(),
false, 0, i ) );
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";
2711 if ( this->
functionSpace()->
template functionSpace<i>()->worldComm().isActive() )
2713 ct_type ct( *
this, ublas::range( nbdof_start, nbdof_start+space->nLocalDof() ),
2714 M_functionspace->template functionSpace<i>()->dof() );
2716 #if defined(FEELPP_ENABLE_MPI_MODE)
2719 if ( this->worldComm().globalSize()>1 && updateOffViews && !this->
functionSpace()->hasEntriesForAllSpaces() )
2721 std::vector<double> dataToSend( space->nLocalDof() );
2722 std::copy( ct.begin(), ct.end(), dataToSend.begin() );
2724 if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2726 fusion::for_each( *M_containersOffProcess, detail::SendContainersOn<i,functionspace_type>( this->
functionSpace(), dataToSend ) );
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 );
2738 if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2740 fusion::for_each( *M_containersOffProcess, detail::InitializeContainersOff<i,functionspace_type>( this->
functionSpace() ) );
2742 #if defined(FEELPP_ENABLE_MPI_MODE)
2745 if ( this->worldComm().globalSize()>1 && updateOffViews && !this->
functionSpace()->hasEntriesForAllSpaces() )
2747 fusion::for_each( *M_containersOffProcess, detail::RecvContainersOff<i,functionspace_type>( this->
functionSpace() ) );
2752 ct_type ct( *fusion::at_c<i>( *M_containersOffProcess ),
2753 ublas::range( 0, space->nLocalDof() ),
2754 M_functionspace->template functionSpace<i>()->dof() );
2756 DVLOG(2) <<
"Element <" << i <<
">::range.size : "<< ct.size()<<
"\n";
2757 DVLOG(2) <<
"Element <" << i <<
">::range.start : "<< ct.start()<<
"\n";
2759 return typename mpl::at_c<element_vector_type,i>::type( space, ct, name );
2764 typename mpl::at_c<element_vector_type,i>::type
2765 element( std::string
const& name =
"u",
bool updateOffViews=
true )
const
2767 size_type nbdof_start = fusion::accumulate( M_functionspace->functionSpaces(),
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>() );
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";
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() );
2783 #if defined(FEELPP_ENABLE_MPI_MODE)
2786 if ( this->worldComm().globalSize()>1 && updateOffViews && !this->
functionSpace()->hasEntriesForAllSpaces() )
2788 std::vector<double> dataToSend( space->nLocalDof() );
2789 std::copy( ct.begin(), ct.end(), dataToSend.begin() );
2791 if ( !M_containersOffProcess ) M_containersOffProcess = boost::in_place();
2793 fusion::for_each( *M_containersOffProcess, detail::SendContainersOn<i,functionspace_type>( this->
functionSpace(), dataToSend ) );
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 );
2805 #if defined(FEELPP_ENABLE_MPI_MODE)
2808 if ( this->worldComm().globalSize()>1 && updateOffViews && !this->
functionSpace()->hasEntriesForAllSpaces() )
2810 fusion::for_each( *M_containersOffProcess, detail::RecvContainersOff<i,functionspace_type>( this->
functionSpace() ) );
2816 ct_type ct( *fusion::at_c<i>( *M_containersOffProcess ),
2817 ublas::range( 0, space->nLocalDof() ),
2818 M_functionspace->template functionSpace<i>()->dof() );
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 );
2830 component_functionspace_ptrtype
const&
compSpace()
const
2832 return M_functionspace->compSpace();
2838 std::vector<WorldComm>
const& worldsComm()
const
2840 return M_functionspace->worldsComm();
2846 WorldComm
const& worldComm()
const
2848 return M_functionspace->worldComm();
2856 return M_functionspace->nDof();
2864 return M_functionspace->nLocalDof();
2872 return M_functionspace->nDofPerComponent();
2878 std::string
const& name()
const
2888 bool isAComponent()
const
2890 return M_ct >= X && M_ct <= Z;
2893 ComponentType component()
const
2895 if ( M_ct < X || M_ct > Z )
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() );
2905 std::string componentToString( ComponentType ct )
const
2907 if ( ct < X || ct > Z )
2909 std::ostringstream __str;
2910 __str <<
"invalid component extraction (should be 0 or 1 or 2): " << ct;
2911 throw std::logic_error( __str.str() );
2926 std::ostringstream __str;
2927 __str <<
"invalid component extraction (should be 0 or 1 or 2): " << ct;
2928 throw std::logic_error( __str.str() );
2942 void setName( std::string
const & __name )
2947 void setFunctionSpace( functionspace_ptrtype space )
2949 M_functionspace = space;
2950 super::init( M_functionspace->dof() );
2954 BOOST_PARAMETER_MEMBER_FUNCTION( (
void ),
2960 ( type,( std::string ),std::string(
"binary" ) )
2961 ( suffix,( std::string ),std::string(
"" ) )
2962 ( sep,( std::string ),std::string(
"" ) )
2965 Feel::detail::ignore_unused_variable_warning( args );
2967 if ( !fs::exists( fs::path( path ) ) )
2969 fs::create_directories( fs::path( path ) );
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 );
2978 if ( type ==
"binary" )
2980 boost::archive::binary_oarchive oa( ofs );
2984 else if ( type ==
"text" )
2986 boost::archive::text_oarchive oa( ofs );
2990 else if ( type ==
"xml" )
2996 BOOST_PARAMETER_MEMBER_FUNCTION(
3003 ( type,( std::string ),std::string(
"binary" ) )
3004 ( suffix,( std::string ),std::string(
"" ) )
3005 ( sep,( std::string ),std::string(
"" ) )
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);
3015 fs::path full_path_dir_sol(fs::current_path());
3016 full_path_dir_sol = full_path_dir_sol/partial_path;
3018 if ( !fs::exists( p ) )
3020 std::ostringstream os2;
3021 os2 << M_name << sep << suffix<<
"-" << this->worldComm().globalSize() <<
"." << this->worldComm().globalRank();
3022 p = fs::path( path ) / os2.str();
3024 if ( !fs::exists( p ) )
3026 LOG(WARNING) <<
"[load] :" << full_path_dir_sol <<
" FILE : " << os1.str() <<
" OR " << os2.str() <<
" DO NOT EXIST" << std::endl ;
3031 LOG(INFO) << p <<
" exists, is is a regular file : " << fs::is_regular_file( p ) <<
"\n";
3032 if ( !fs::is_regular_file( p ) )
3035 LOG(WARNING) <<
"[load] : " << full_path_dir_sol << p <<
" is not a regular_file !" << std::endl;
3039 fs::ifstream ifs( p );
3041 if ( type ==
"binary" )
3043 boost::archive::binary_iarchive ia( ifs );
3047 else if ( type ==
"text" )
3049 boost::archive::text_iarchive ia( ifs );
3053 else if ( type ==
"xml" )
3062 void printMatlab( std::string fname,
bool gmsh =
false )
const
3064 container_type m( *
this );
3067 auto relation = this->
functionSpace()->dof()->pointIdToDofRelation();
3068 for(
size_type i = 0; i < this->localSize(); ++i )
3070 m[relation.first[i]-1] = this->operator[]( i );
3073 m.printMatlab( fname );
3079 friend class boost::serialization::access;
3081 template<
class Archive>
3082 void serialize( Archive & ar,
const unsigned int version )
3085 ar & boost::serialization::make_nvp(
"name", M_name );
3086 DVLOG(2) <<
"got name " << M_name <<
"\n";
3088 if ( Archive::is_saving::value )
3092 ar & boost::serialization::make_nvp(
"size", s );
3094 std::vector<int> no = M_functionspace->basisOrder();
3096 std::string family = M_functionspace->basisName();
3099 ar & boost::serialization::make_nvp(
"order", no );
3101 ar & boost::serialization::make_nvp(
"family", family );
3104 typename container_type::const_iterator it = this->begin();
3105 typename container_type::const_iterator en = this->end();
3107 for (
size_type i = 0; it != en; ++it, ++i )
3109 T value = this->operator[]( i );
3110 std::ostringstream v_str;
3111 v_str <<
"value_" << i;
3115 ar & boost::serialization::make_nvp( v_str.str().c_str(), value );
3119 if ( Archive::is_loading::value )
3124 ar & boost::serialization::make_nvp(
"size", s );
3127 DVLOG(2) <<
"loading ublas::vector of size " << s <<
"\n";
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() );
3133 std::vector<int> order;
3135 ar & boost::serialization::make_nvp(
"order", order );
3138 ar & boost::serialization::make_nvp(
"family", family );
3141 auto orders = M_functionspace->basisOrder();
3143 if ( order != orders )
3144 throw std::logic_error( ( boost::format(
"load function: invalid polynomial order, read %1% but has %2%" ) % order % orders ).str() );
3146 std::string bname = M_functionspace->basisName();
3148 if ( family != bname )
3149 throw std::logic_error( ( boost::format(
"load function: invalid polynomial family, read %1% but has %2%" ) % family % bname ).str() );
3155 value_type value( 0 );
3156 std::ostringstream v_str;
3157 v_str <<
"value_" << i;
3159 ar & boost::serialization::make_nvp( v_str.str().c_str(), value );
3161 this->operator[]( i ) = value;
3173 functionspace_ptrtype M_functionspace;
3183 mutable boost::optional<container_vector_type> M_containersOffProcess;
3191 typedef Element<value_type> element_type;
3192 typedef Element<value_type> real_element_type;
3193 typedef boost::shared_ptr<element_type> element_ptrtype;
3195 typedef std::map< size_type, std::vector< size_type > > proc_dist_map_type;
3209 size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3210 periodicity_type periodicity = periodicity_type(),
3211 std::vector<WorldComm>
const& _worldsComm = Environment::worldsComm(nSpaces) )
3213 M_worldsComm( _worldsComm ),
3214 M_worldComm( new WorldComm( _worldsComm[0] ) )
3216 this->
init( mesh, mesh_components, periodicity );
3220 std::vector<Dof >
const& dofindices,
3221 periodicity_type periodicity = periodicity_type(),
3222 std::vector<WorldComm>
const& _worldsComm = Environment::worldsComm(nSpaces) )
3224 M_worldsComm( _worldsComm ),
3225 M_worldComm( new WorldComm( _worldsComm[0] ) )
3227 this->
init( mesh, 0, dofindices, periodicity );
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 )
3237 return pointer_type(
new functionspace_type( __m, mesh_components ) );
3241 static pointer_type
New( mesh_ptrtype
const& __m, std::vector<Dof >
const& dofindices )
3245 #if !defined( FEELPP_ENABLE_MPI_MODE)
3246 BOOST_PARAMETER_MEMBER_FUNCTION( ( pointer_type ),
3253 ( components, (
size_type ), MESH_RENUMBER | MESH_CHECK )
3254 ( periodicity,*,periodicity_type() )
3258 return NewImpl( mesh, components, periodicity );
3260 static pointer_type NewImpl( mesh_ptrtype
const& __m,
3261 size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3262 periodicity_type periodicity = periodicity_type() )
3264 return pointer_type(
new functionspace_type( __m, mesh_components, periodicity ) );
3267 BOOST_PARAMETER_MEMBER_FUNCTION( ( pointer_type ),
3274 ( worldscomm, *, Environment::worldsComm(nSpaces) )
3275 ( components, (
size_type ), MESH_RENUMBER | MESH_CHECK )
3276 ( periodicity,*,periodicity_type() )
3280 return NewImpl( mesh, worldscomm, components, periodicity );
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() )
3289 return pointer_type(
new functionspace_type( __m, mesh_components, periodicity, worldsComm ) );
3296 void init( mesh_ptrtype
const& mesh,
3297 size_type mesh_components = MESH_RENUMBER | MESH_CHECK,
3298 periodicity_type periodicity = periodicity_type() )
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";
3306 this->
init( mesh, mesh_components, periodicity, std::vector<Dof >(), mpl::bool_<is_composite>() );
3310 void init( mesh_ptrtype
const& mesh,
3312 std::vector<Dof >
const& dofindices,
3313 periodicity_type periodicity = periodicity_type() )
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";
3322 this->
init( mesh, mesh_components, periodicity, dofindices, mpl::bool_<is_composite>() );
3335 void setWorldsComm( std::vector<WorldComm>
const& _worldsComm )
3337 M_worldsComm=_worldsComm;
3339 void setWorldComm( WorldComm
const& _worldComm )
3341 M_worldComm.reset(
new WorldComm( _worldComm ) );
3344 std::vector<WorldComm>
const& worldsComm()
const
3346 return M_worldsComm;
3348 WorldComm
const& worldComm()
const
3350 return *M_worldComm;
3353 bool hasEntriesForAllSpaces()
3355 return (this->
template mesh<0>()->worldComm().localSize() == this->
template mesh<0>()->worldComm().globalSize() );
3368 DVLOG(2) <<
"Update function space after a change in the mesh\n";
3385 template<
typename Context_t>
3387 ptsInContext( Context_t
const & context, mpl::int_<2> )
const
3391 typedef boost::shared_ptr<gmc_interp_type> gmc_interp_ptrtype;
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() ) );
3402 return __c_interp->xReal();
3411 return N_COMPONENTS;
3419 return this->
nDof( mpl::bool_<is_composite>() );
3427 return this->
nLocalDof( mpl::bool_<is_composite>() );
3432 return this->nLocalDofWithGhost( mpl::bool_<is_composite>() );
3437 return this->nLocalDofWithoutGhost( mpl::bool_<is_composite>() );
3440 size_type nLocalDofWithGhostOnProc(
const int proc )
const
3442 return this->nLocalDofWithGhostOnProc( proc, mpl::bool_<is_composite>() );
3445 size_type nLocalDofWithoutGhostOnProc(
const int proc)
const
3447 return this->nLocalDofWithoutGhostOnProc( proc, mpl::bool_<is_composite>() );
3487 size_type start = fusion::accumulate( this->
functionSpaces(),
size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<true> >( proc, this->worldsComm(),
true,0,i ) );
3493 size_type start = fusion::accumulate( this->
functionSpaces(),
size_type( 0 ), detail::NLocalDofOnProc<mpl::bool_<false> >( proc, this->worldsComm(),
true,0,i ) );
3497 uint16_type nSubFunctionSpace()
const
3528 typename GetMesh<mesh_ptrtype,i>::type
3531 return mesh<i>( is_shared_ptr<mesh_ptrtype>() );
3534 typename GetMesh<mesh_ptrtype,i>::type
3535 mesh( mpl::bool_<true> )
const
3540 typename GetMesh<mesh_ptrtype,i>::type
3541 mesh( mpl::bool_<false> )
const
3543 return fusion::at_c<i>(M_mesh);
3562 return basisName( mpl::bool_<( nSpaces>1 )>() );
3564 std::string
basisName( mpl::bool_<true> )
const
3566 return fusion::accumulate( this->
functionSpaces(), std::string(), detail::BasisName() );
3568 std::string
basisName( mpl::bool_<false> )
const
3570 return this->
basis()->familyName();
3581 return basisOrder( mpl::bool_<( nSpaces>1 )>() );
3583 std::vector<int>
basisOrder( mpl::bool_<true> )
const
3585 return fusion::accumulate( this->
functionSpaces(), std::vector<int>(), detail::BasisOrder() );
3587 std::vector<int>
basisOrder( mpl::bool_<false> )
const
3589 std::vector<int> o( 1 );
3590 o[0]=basis_type::nOrder;
3597 reference_element_ptrtype
const&
fe()
const
3605 gm_ptrtype
const&
gm()
const
3607 return M_mesh->gm();
3613 gm1_ptrtype
const&
gm1()
const
3615 return M_mesh->gm1();
3669 dof_ptrtype
const&
dof()
const
3696 functionspace_vector_type
const&
3699 return M_functionspaces;
3702 std::vector<std::vector<size_type> > dofIndexSplit()
3706 uint16_type cptSpaces=0;
3708 std::vector<std::vector<size_type> > is(nSpaces);
3709 auto initial = boost::make_tuple( cptSpaces,start,is );
3712 auto startSplit = boost::fusion::fold(
functionSpaces(), boost::make_tuple(0,0), detail::computeStartOfFieldSplit() ).template get<1>();
3714 auto result = boost::fusion::fold(
functionSpaces(), initial, detail::computeNDofForEachSpace(startSplit) );
3715 is = result.template get<2>();
3719 for (
int proc = 0; proc<this->worldComm().globalSize(); ++proc )
3721 this->worldComm().globalComm().barrier();
3722 if ( proc==this->worldComm().globalRank() )
3724 std::cout <<
"proc " << proc <<
"\n";
3725 std::cout <<
"split size=" << result.template get<2>().size() <<
" nspace=" << nSpaces <<
"\n";
3726 std::cout <<
"split:\n";
3731 for (
int s= 0; s < nSpaces; ++s )
3733 std::cout <<
"space: " << is[s].size() <<
"\n";
3735 for (
int i = 0; i < is[s].size(); ++i )
3737 std::cout << is[s][i] <<
" ";
3740 std::cout <<
"\n\n";
3749 std::vector<std::vector<size_type> > is;
3750 is.push_back( std::vector<size_type>(
nLocalDof() ) );
3753 BOOST_FOREACH(
auto& i, is[0] )
3766 element_type u( this->shared_from_this(), name );
3777 element_ptrtype u(
new element_type( this->shared_from_this(), name ) );
3786 typename mpl::at_c<functionspace_vector_type,i>::type
3789 return fusion::at_c<i>( M_functionspaces );
3804 trace_functionspace_ptrtype
3810 template<
typename RangeT>
3811 trace_functionspace_ptrtype
3812 trace( RangeT range )
const
3814 return trace_functionspace_type::New(
mesh()->
trace( range ) );
3821 trace_trace_functionspace_ptrtype
3827 template<
typename RangeT>
3828 trace_trace_functionspace_ptrtype
3831 return trace_trace_functionspace_type::New(
mesh()->
wireBasket( range ) );
3883 void printInfo()
const
3885 LOG(INFO) <<
" number of components : " <<
qDim() <<
"\n";
3886 LOG(INFO) <<
" n Global Dof : " <<
nDof() <<
"\n";
3887 LOG(INFO) <<
" n Local Dof : " <<
nLocalDof() <<
"\n";
3893 FunctionSpace( FunctionSpace
const& __fe )
3895 M_worldsComm( __fe.M_worldsComm ),
3896 M_worldComm( __fe.M_worldComm ),
3897 M_mesh( __fe.M_mesh ),
3904 DVLOG(2) <<
"copying FunctionSpace\n";
3911 void init( mesh_ptrtype
const& mesh,
3913 periodicity_type
const& periodicity,
3914 std::vector<Dof >
const& dofindices,
3915 mpl::bool_<false> );
3916 void init( mesh_ptrtype
const& mesh,
3918 periodicity_type
const& periodicity,
3919 std::vector<Dof >
const& dofindices,
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;
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;
3941 friend class ComponentSpace;
3942 class ComponentSpace
3945 typedef FunctionSpace<A0,A1,A2,A3,A4> functionspace_type;
3946 typedef functionspace_type* functionspace_ptrtype;
3947 typedef functionspace_type
const* functionspace_cptrtype;
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;
3953 ComponentSpace( FunctionSpace<A0,A1,A2,A3,A4> * __functionspace,
3956 M_functionspace( __functionspace ),
3961 return operator()( mpl::bool_<functionspace_type::is_scalar>() );
3963 component_functionspace_ptrtype
operator()( mpl::bool_<true> )
3965 FEELPP_ASSERT( 0 ).error(
"invalid call for component space extraction" );
3967 return component_functionspace_ptrtype();
3969 component_functionspace_ptrtype
operator()( mpl::bool_<false> )
3976 FunctionSpace<A0,A1,A2,A3,A4> * M_functionspace;
3977 mesh_ptrtype M_mesh;
3981 template<
typename ElementRange,
typename OnExpr>
3983 on( ElementRange
const& _range,
3984 OnExpr
const& _expr )
3986 M_constraints.push_back(
vf::project( this->shared_from_this(), _range, _expr ) );
3995 std::vector<WorldComm> M_worldsComm;
3996 boost::shared_ptr<WorldComm> M_worldComm;
3999 mesh_ptrtype M_mesh;
4014 mutable boost::optional<region_tree_ptrtype>
M_rt;
4022 functionspace_vector_type M_functionspaces;
4024 proc_dist_map_type procDistMap;
4026 std::list<Element> M_constraints;
4030 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4031 template<
typename T,
typename Cont>
4034 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4035 template<
typename T,
typename Cont>
4038 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4039 template<
typename T,
typename Cont>
4042 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4043 template<
typename T,
typename Cont>
4046 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4050 periodicity_type
const& periodicity,
4051 std::vector<Dof>
const& dofindices,
4054 DVLOG(2) <<
"calling init(<space>) begin\n";
4055 DVLOG(2) <<
"calling init(<space>) is_periodic: " << is_periodic <<
"\n";
4059 if ( basis_type::nDofPerEdge || nDim >= 3 )
4060 mesh_components |= MESH_UPDATE_EDGES;
4067 if ( basis_type::nDofPerFace || is_continuous || nDim >= 3 )
4068 mesh_components |= MESH_UPDATE_FACES;
4070 M_mesh->components().set( mesh_components );
4072 M_mesh->updateForUse();
4076 M_mesh->removeFacesFromBoundary( { periodicity.tag1(), periodicity.tag2() } );
4079 M_ref_fe = basis_ptrtype(
new basis_type );
4081 M_dof = dof_ptrtype(
new dof_type( M_ref_fe, fusion::at_c<0>(periodicity), this->worldsComm()[0] ) );
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";
4087 M_dof->build( M_mesh );
4096 M_comp_space = component_functionspace_ptrtype(
new component_functionspace_type( M_mesh,
4097 MESH_COMPONENTS_DEFAULTS,
4099 std::vector<WorldComm>( 1,this->worldsComm()[0] ) ) );
4102 DVLOG(2) <<
"nb dim : " << qDim() <<
"\n";
4103 DVLOG(2) <<
"nb dof : " << nDof() <<
"\n";
4104 DVLOG(2) <<
"nb dof per component: " << nDofPerComponent() <<
"\n";
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";
4115 DVLOG(2) <<
"calling init(<space>) end\n";
4119 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4123 periodicity_type
const& periodicity,
4124 std::vector<Dof>
const& dofindices,
4127 DVLOG(2) <<
"calling init(<composite>) begin\n";
4131 fusion::for_each( M_functionspaces, detail::InitializeSpace<mesh_ptrtype,periodicity_type>( __m, periodicity,
4132 dofindices, this->worldsComm() ) );
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";
4139 proc_dist_map_type emptyMap;
4140 procDistMap = fusion::accumulate( M_functionspaces,
4142 detail::searchIndicesBySpace<proc_dist_map_type>() );
4144 #else // new version with MPI
4146 if ( this->worldComm().globalSize()>1 )
4148 if ( this->hasEntriesForAllSpaces() )
4152 DVLOG(2) <<
"init(<composite>) type hasEntriesForAllSpaces\n";
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 )
4159 nLocalDofWithoutGhostWorld[proc] = this->nLocalDofWithoutGhostOnProc(proc);
4160 nLocalDofWithGhostWorld[proc] = this->nLocalDofWithGhostOnProc(proc);
4163 startDofGlobalCluster[0]=0;
4164 for (
int p=1;p<worldsize;++p)
4166 startDofGlobalCluster[p] = startDofGlobalCluster[p-1] + nLocalDofWithoutGhostWorld[p-1];
4170 auto dofInitTool=detail::updateDataMapProcessStandard<dof_type>( this->worldsComm(), this->worldComm(),
4171 this->nSubFunctionSpace()-1,
4172 startDofGlobalCluster,
4173 nLocalDofWithoutGhostWorld,
4174 nLocalDofWithGhostWorld
4176 fusion::for_each( M_functionspaces, dofInitTool );
4178 M_dof = dofInitTool.dataMap();
4179 M_dof->setNDof( this->nDof() );
4186 DVLOG(2) <<
"init(<composite>) type Not hasEntriesForAllSpaces\n";
4189 WorldComm mixSpaceWorldComm = this->worldsComm()[0];
4191 if ( this->worldsComm().size()>1 )
4192 for (
int i=1; i<( int )this->worldsComm().size(); ++i )
4194 mixSpaceWorldComm = mixSpaceWorldComm + this->worldsComm()[i];
4197 this->setWorldComm( mixSpaceWorldComm );
4201 auto dofInitTool=detail::updateDataMapProcess<dof_type>( this->worldsComm(), mixSpaceWorldComm, this->nSubFunctionSpace()-1 );
4202 fusion::for_each( M_functionspaces, dofInitTool );
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();
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() );
4227 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
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";
4237 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4241 return M_dof->nDof();
4244 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
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";
4254 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4259 return M_dof->nLocalDofWithGhost();
4262 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4264 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhost( mpl::bool_<true> )
const
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";
4272 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4274 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhost( mpl::bool_<false> )
const
4276 return M_dof->nLocalDofWithGhost();
4279 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4281 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhostOnProc(
const int proc, mpl::bool_<true> )
const
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";
4289 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4291 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithGhostOnProc(
const int proc, mpl::bool_<false> )
const
4293 return M_dof->nLocalDofWithGhost(proc);
4296 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4298 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhost( mpl::bool_<true> )
const
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";
4306 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4308 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhost( mpl::bool_<false> )
const
4310 return M_dof->nLocalDofWithoutGhost();
4313 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4315 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhostOnProc(
const int proc, mpl::bool_<true> )
const
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";
4323 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4325 FunctionSpace<A0, A1, A2, A3, A4>::nLocalDofWithoutGhostOnProc(
const int proc, mpl::bool_<false> )
const
4327 return M_dof->nLocalDofWithoutGhost(proc);
4330 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4334 M_dof->rebuildDofPoints( *M_mesh );
4337 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4341 fusion::for_each( M_functionspaces, detail::rebuildDofPointsTool() );
4344 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4348 scalar_type EPS=1E-13;
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() );
4359 for (
size_type __i = 0; it != en; ++__i, ++it )
4361 __bb.make( it->G() );
4363 for (
unsigned k=0; k < __bb.min.size(); ++k )
4369 __rt->addBox( __bb.min, __bb.max, it->id() );
4376 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4382 scalar_type EPS=1E-13;
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() );
4393 for (
size_type __i = 0; it != en; ++__i, ++it )
4395 __bb.make( it->G() );
4397 for (
unsigned k=0; k < __bb.min.size(); ++k )
4403 __rt->addBox( __bb.min, __bb.max, it->id() );
4413 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4417 if ( !hasRegionTree() )
4424 region_tree_type::pbox_set_type boxlst;
4425 __rt->findBoxesAtPoint( pt, boxlst );
4427 typedef typename gm_type::Inverse inv_trans_type;
4428 typename gm_type::reference_convex_type refelem;
4432 region_tree_type::pbox_set_type::const_iterator it = boxlst.begin();
4433 region_tree_type::pbox_set_type::const_iterator ite = boxlst.end();
4435 for ( ; it != ite; ++it )
4437 inv_trans_type __git( M_mesh->gm(), M_mesh->element( ( *it )->id ), this->worldComm().subWorldCommSeq() );
4442 DVLOG(2) <<
"[FunctionSpace::findPoint] id : " << cv_stored <<
"\n";
4444 __git.setXReal( pt );
4450 boost::tie( isin, dmin ) = refelem.isIn( ptr );
4451 DVLOG(2) <<
"[FunctionSpace::findPoint] isin: " << isin <<
" dmin: " << dmin <<
"\n";
4453 closest = ( dmin > closest.second )?std::make_pair( cv_stored, dmin ):closest;
4457 DVLOG(2) <<
"[FunctionSpace::findPoint] id of the convex where " << pt <<
" belongs : " << cv_stored <<
"\n";
4458 DVLOG(2) <<
"[FunctionSpace::findPoint] ref coordinate: " << ptr <<
"\n";
4473 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4474 template<
typename Y,
typename Cont>
4479 M_ct( NO_COMPONENT ),
4480 M_containersOffProcess( boost::none )
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 )
4488 M_functionspace( __e.M_functionspace ),
4489 M_name( __e.M_name ),
4490 M_start( __e.M_start ),
4492 M_containersOffProcess( __e.M_containersOffProcess )
4494 DVLOG(2) <<
"Element<copy>::range::start = " << this->start() <<
"\n";
4495 DVLOG(2) <<
"Element<copy>::range::size = " << this->size() <<
"\n";
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,
4504 ComponentType __ct )
4506 super( __functionspace->dof() ),
4507 M_functionspace( __functionspace ),
4511 M_containersOffProcess( boost::none )
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";
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,
4525 ComponentType __ct )
4528 M_functionspace( __functionspace ),
4532 M_containersOffProcess( boost::none )
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();
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()
4548 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4549 template<
typename Y,
typename Cont>
4551 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::initFromSpace( functionspace_ptrtype
const& __functionspace,
4552 container_type
const& __c )
4554 M_functionspace = __functionspace;
4555 ( container_type )*
this = __c;
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>&
4565 M_functionspace = __e.M_functionspace;
4567 if ( __e.M_name !=
"unknown" )
4568 M_name = __e.M_name;
4570 M_start = __e.M_start;
4572 M_containersOffProcess = __e.M_containersOffProcess;
4573 this->resize( M_functionspace->nLocalDof() );
4575 this->outdateGlobalValues();
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>&
4587 if ( __v.size() != this->size() )
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() );
4595 this->outdateGlobalValues();
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
4609 this->updateGlobalValues();
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 );
4618 if ( functionSpace()->findPoint( __x, __cv_id, __x_ref ) || extrapolate )
4620 #if !defined( NDEBUG )
4621 DVLOG(2) <<
"Point " << __x <<
" is in element " << __cv_id <<
" pt_ref=" << __x_ref <<
"\n";
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 ) );
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 ),
4636 pc_ptrtype pc(
new pc_type( this->functionSpace()->fe(), pts ) );
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(),
4645 found_pt[ rank ] = 1;
4647 #if defined(FEELPP_HAS_MPI)
4652 mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4656 global_found_pt[ 0 ] = found_pt[ 0 ];
4659 id_type __id( this->
id( *fectx ) );
4662 #if defined(FEELPP_HAS_MPI)
4663 DVLOG(2) <<
"sending interpolation context to all processors from " << functionSpace()->mesh()->comm().rank() <<
"\n";
4665 if ( functionSpace()->mesh()->comm().size() > 1 )
4667 mpi::broadcast( functionSpace()->mesh()->comm(), __id, functionSpace()->mesh()->comm().rank() );
4677 #if defined(FEELPP_HAS_MPI)
4679 if ( functionSpace()->mesh()->comm().size() > 1 )
4682 mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4689 for ( ; i < global_found_pt.size(); ++i )
4690 if ( global_found_pt[i] != 0 )
4692 DVLOG(2) <<
"processor " << i <<
" has the point " << __x <<
"\n";
4701 DVLOG(2) <<
"receiving interpolation context from processor " << i <<
"\n";
4702 #if defined(FEELPP_HAS_MPI)
4704 if ( functionSpace()->mesh()->comm().size() > 1 )
4705 mpi::broadcast( functionSpace()->mesh()->comm(), __id, i );
4712 LOG(WARNING) <<
"no processor seems to have the point " << __x <<
"\n";
4719 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4720 template<
typename Y,
typename Cont>
4721 template<
typename Context_t>
4724 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::id_( Context_t
const & context, id_array_type& v )
const
4726 if ( !this->areGlobalValuesUpdated() )
4727 this->updateGlobalValues();
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() );
4737 const uint16_type nq = context.xRefs().size2();
4740 for (
int l = 0; l < basis_type::nDof; ++l )
4742 const int ncdof = is_product?nComponents1:1;
4744 for (
typename array_type::index c1 = 0; c1 < ncdof; ++c1 )
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 ) );
4751 FEELPP_ASSERT( gdof < this->size() )
4753 ( l )( c1 )( ldof )( gdof )
4754 ( this->size() )( this->localSize() )
4755 ( this->firstLocalIndex() )( this->lastLocalIndex() )
4756 .warn(
"FunctionSpace::Element invalid access index" );
4758 value_type v_ = this->globalValue( gdof );
4762 for ( uint16_type q = 0; q < nq; ++q )
4764 for (
typename array_type::index i = 0; i < nComponents1; ++i )
4767 v[q]( i,0 ) += v_*context.id( ldof, i, 0, q );
4777 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4778 template<
typename Y,
typename Cont>
4780 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::idInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const
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;
4788 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
4795 gm_ptrtype __gm = this->functionSpace()->gm();
4798 if ( it==it_end )
return;
4802 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
4803 matrix_node_type pts( nbCoord, nbPtsElt );
4806 geopc_ptrtype __geopc(
new geopc_type( __gm, pts ) );
4807 pc_ptrtype __pc(
new pc_type( this->functionSpace()->fe(), pts ) );
4809 gmc_ptrtype __c(
new gmc_type( __gm,
4810 this->functionSpace()->mesh()->element( it->first ),
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(),
4821 for ( ; it!=it_end; ++it )
4823 nbPtsElt = it->second.size();
4826 itL=it->second.begin();
4827 itL_end=it->second.end();
4830 pts= matrix_node_type( nbCoord, nbPtsElt );
4832 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
4833 ublas::column( pts, i ) = boost::get<1>( *itL );
4836 __geopc->update( pts );
4839 __pc->update( pts );
4841 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
4842 __ctx->update( __c, __pc );
4845 id_type __id( this->
id( *__ctx ) );
4848 itL=it->second.begin();
4850 for ( uint k=0; k<nbPtsElt; ++k,++itL )
4852 for (
typename array_type::index i = 0; i < nComponents1; ++i )
4854 v[boost::get<0>( *itL )]( i,0 ) = __id( i,0,k );
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
4869 this->updateGlobalValues();
4873 std::vector<int> found_pt( functionSpace()->mesh()->comm().size(), 0 );
4874 std::vector<int> global_found_pt( functionSpace()->mesh()->comm().size(), 0 );
4876 if ( functionSpace()->findPoint( __x, __cv_id, __x_ref ) )
4878 #if !defined( NDEBUG )
4879 DVLOG(2) <<
"Point " << __x <<
" is in element " << __cv_id <<
" pt_ref=" << __x_ref <<
"\n";
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 ) );
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 ),
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(),
4903 found_pt[ functionSpace()->mesh()->comm().rank() ] = 1;
4905 #if defined(FEELPP_HAS_MPI)
4907 if ( functionSpace()->mesh()->comm().size() > 1 )
4910 mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4914 global_found_pt[ 0 ] = found_pt[ 0 ];
4917 grad_type g_( this->grad( *fectx ) );
4919 #if defined(FEELPP_HAS_MPI)
4920 DVLOG(2) <<
"sending interpolation context to all processors from " << functionSpace()->mesh()->comm().rank() <<
"\n";
4922 if ( functionSpace()->mesh()->comm().size() > 1 )
4924 mpi::broadcast( functionSpace()->mesh()->comm(), g_, functionSpace()->mesh()->comm().rank() );
4934 #if defined(FEELPP_HAS_MPI)
4936 if ( functionSpace()->mesh()->comm().size() > 1 )
4939 mpi::all_reduce( functionSpace()->mesh()->comm(), found_pt, global_found_pt, detail::vector_plus<int>() );
4946 for ( ; i < global_found_pt.size(); ++i )
4947 if ( global_found_pt[i] != 0 )
4949 DVLOG(2) <<
"processor " << i <<
" has the point " << __x <<
"\n";
4958 DVLOG(2) <<
"receiving interpolation context from processor " << i <<
"\n";
4959 #if defined(FEELPP_HAS_MPI)
4961 if ( functionSpace()->mesh()->comm().size() > 1 )
4962 mpi::broadcast( functionSpace()->mesh()->comm(), g_, i );
4971 Warning() <<
"no processor seems to have the point " << __x <<
"\n";
4978 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
4979 template<
typename Y,
typename Cont>
4980 template<
typename ContextType>
4983 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::grad_( ContextType
const & context, grad_array_type& v )
const
4985 if ( !this->areGlobalValuesUpdated() )
4986 this->updateGlobalValues();
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() );
4999 for (
int l = 0; l < basis_type::nDof; ++l )
5001 const int ncdof = is_product?nComponents1:1;
5003 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5010 ( l )( c1 )( ldof )( gdof )
5011 ( this->size() )( this->localSize() )
5012 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5013 .error(
"FunctionSpace::Element invalid access index" );
5016 value_type v_ = this->globalValue( gdof );
5018 for (
size_type q = 0; q < context.xRefs().size2(); ++q )
5020 for (
int k = 0; k < nComponents1; ++k )
5021 for (
int j = 0; j < nRealDim; ++j )
5023 v[q]( k,j ) += v_*context.grad( ldof, k, j, q );
5030 std::cout <<
"xref = " << context.xRefs() <<
"\n";
5031 std::cout <<
"xreal = " << context.xReal() <<
"\n";
5032 std::cout <<
"pts = " << context.G() <<
"\n";
5034 for (
int l = 0; l < basis_type::nDof; ++l )
5036 for (
typename array_type::index i = 0; i < nDim; ++i )
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";
5044 for (
int l = 0; l < basis_type::nDof; ++l )
5046 for (
int j = 0; j < nRealDim; ++j )
5047 for (
typename array_type::index i = 0; i < nDim; ++i )
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";
5055 for (
int k = 0; k < nComponents1; ++k )
5056 for (
int j = 0; j < nRealDim; ++j )
5058 for (
size_type q = 0; q < context.xRefs().size2(); ++q )
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";
5070 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5071 template<
typename Y,
typename Cont>
5073 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::gradInterpolate( matrix_node_type __ptsReal, grad_array_type& v )
const
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;
5080 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5086 gm_ptrtype __gm = this->functionSpace()->gm();
5089 if ( it==it_end )
return;
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 ) );
5097 gmc_ptrtype __c(
new gmc_type( __gm,
5098 this->functionSpace()->mesh()->element( it->first ),
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(),
5108 for ( ; it!=it_end; ++it )
5110 nbPtsElt = it->second.size();
5113 itL=it->second.begin();
5114 itL_end=it->second.end();
5117 pts= matrix_node_type( nbCoord, nbPtsElt );
5119 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5120 ublas::column( pts, i ) = boost::get<1>( *itL );
5123 __geopc->update( pts );
5127 __pc->update( pts );
5128 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5129 __ctx->update( __c, __pc );
5132 grad_type __grad( this->grad( *__ctx ) );
5135 itL=it->second.begin();
5137 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5139 for (
typename array_type::index i = 0; i < nComponents1; ++i )
5140 for ( uint j = 0; j < nRealDim; ++j )
5142 v[boost::get<0>( *itL )]( i,j ) = __grad( i,j,k );
5151 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5152 template<
typename Y,
typename Cont>
5153 template<
typename ContextType>
5156 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::div_( ContextType
const & context, div_array_type& v )
const
5160 if ( !this->areGlobalValuesUpdated() )
5161 this->updateGlobalValues();
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() );
5171 const size_type Q = context.xRefs().size2();
5173 for (
int l = 0; l < basis_type::nDof; ++l )
5175 const int ncdof = is_product?nComponents1:1;
5177 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5184 ( l )( c1 )( ldof )( gdof )
5185 ( this->size() )( this->localSize() )
5186 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5187 .error(
"FunctionSpace::Element invalid access index" );
5189 value_type v_ = this->globalValue( gdof );
5194 std::cout <<
"v(" << gdof <<
")=" << v_ <<
"\n";
5195 std::cout <<
"context.div(" << ldof <<
"," << q <<
")="
5196 << context.div( ldof, 0, 0, q ) <<
"\n" ;
5198 v[q]( 0,0 ) += v_*context.div( ldof, 0, 0, q );
5205 if ( !this->areGlobalValuesUpdated() )
5206 this->updateGlobalValues();
5208 for (
int l = 0; l < basis_type::nDof; ++l )
5210 const int ncdof = is_product?nComponents1:1;
5212 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5219 ( l )( c1 )( ldof )( gdof )
5220 ( this->size() )( this->localSize() )
5221 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5222 .error(
"FunctionSpace::Element invalid access index" );
5224 value_type v_ = this->globalValue( gdof );
5226 for (
int k = 0; k < nComponents1; ++k )
5228 for (
typename array_type::index i = 0; i < nDim; ++i )
5230 for (
size_type q = 0; q < context.xRefs().size2(); ++q )
5232 v[q]( 0,0 ) += v_*context.gmContext()->B( q )( k, i )*context.pc()->grad( ldof, k, i, q );
5242 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5243 template<
typename Y,
typename Cont>
5245 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::divInterpolate( matrix_node_type __ptsReal, div_array_type& v )
const
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;
5253 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5260 gm_ptrtype __gm = this->functionSpace()->gm();
5263 if ( it==it_end )
return;
5267 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5268 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5284 for ( ; it!=it_end; ++it )
5286 nbPtsElt = it->second.size();
5289 itL=it->second.begin();
5290 itL_end=it->second.end();
5293 pts= matrix_node_type( nbCoord, nbPtsElt );
5295 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5296 ublas::column( pts, i ) = boost::get<1>( *itL );
5299 __geopc->update( pts );
5300 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5302 __pc->update( pts );
5303 __ctx->update( __c, __pc );
5306 div_type __div( this->div( *__ctx ) );
5311 itL=it->second.begin();
5313 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5314 v[boost::get<0>( *itL )]( 0,0 ) = __div( 0,0,k );
5321 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5322 template<
typename Y,
typename Cont>
5323 template<
typename ContextType>
5326 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curl_( ContextType
const & context, curl_array_type& v )
const
5328 if ( !this->areGlobalValuesUpdated() )
5329 this->updateGlobalValues();
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() );
5339 for (
int l = 0; l < basis_type::nDof; ++l )
5341 const int ncdof = is_product?nComponents1:1;
5343 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5350 ( l )( c1 )( ldof )( gdof )
5351 ( this->size() )( this->localSize() )
5352 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5353 .error(
"FunctionSpace::Element invalid access index" );
5356 value_type v_ = this->globalValue( gdof );
5357 const uint16_type nq = context.xRefs().size2();
5359 for ( uint16_type q = 0; q < nq ; ++q )
5363 for (
typename array_type::index i = 0; i < nDim; ++i )
5365 v[q]( i,0 ) += v_*context.curl( ldof, i, 0, q );
5369 else if ( nDim == 2 )
5371 v[q]( 0,0 ) += v_*context.curl( ldof, 0, 0, q );
5380 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5381 template<
typename Y,
typename Cont>
5382 template<
typename ContextType>
5385 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curl_( ContextType
const & context, comp_curl_array_type& v,
int comp )
const
5387 if ( !this->areGlobalValuesUpdated() )
5388 this->updateGlobalValues();
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() );
5399 for (
int l = 0; l < basis_type::nDof; ++l )
5401 const int ncdof = is_product?nComponents1:1;
5403 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5410 ( l )( c1 )( ldof )( gdof )
5411 ( this->size() )( this->localSize() )
5412 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5413 .error(
"FunctionSpace::Element invalid access index" );
5416 value_type v_ = this->globalValue( gdof );
5417 const uint16_type nq = context.xRefs().size2();
5419 for ( uint16_type q = 0; q < nq ; ++q )
5423 v[q]( 0,0 ) += v_*context.curl( ldof, comp, 0, q );
5426 else if ( nDim == 2 )
5428 v[q]( 0,0 ) += v_*context.curl( ldof, 2, 0, q );
5437 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5438 template<
typename Y,
typename Cont>
5440 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlInterpolate( matrix_node_type __ptsReal, curl_array_type& v )
const
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;
5448 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5455 gm_ptrtype __gm = this->functionSpace()->gm();
5458 if ( it==it_end )
return;
5462 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5463 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5480 for ( ; it!=it_end; ++it )
5482 nbPtsElt = it->second.size();
5485 itL=it->second.begin();
5486 itL_end=it->second.end();
5489 pts= matrix_node_type( nbCoord, nbPtsElt );
5491 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5492 ublas::column( pts, i ) = boost::get<1>( *itL );
5495 __geopc->update( pts );
5496 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5498 __pc->update( pts );
5499 __ctx->update( __c, __pc );
5502 curl_type __curl( this->
curl( *__ctx ) );
5505 itL=it->second.begin();
5509 itL=it->second.begin();
5511 for ( uint k=0; k<nbPtsElt; ++k,++itL )
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 );
5519 else if ( nDim == 2 )
5521 itL=it->second.begin();
5523 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5524 v[boost::get<0>( *itL )]( 0,0 ) = __curl( 0,0,k );
5531 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5532 template<
typename Y,
typename Cont>
5534 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlxInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const
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;
5542 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5549 gm_ptrtype __gm = this->functionSpace()->gm();
5552 if ( it==it_end )
return;
5556 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5557 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5574 for ( ; it!=it_end; ++it )
5576 nbPtsElt = it->second.size();
5579 itL=it->second.begin();
5580 itL_end=it->second.end();
5583 pts= matrix_node_type( nbCoord, nbPtsElt );
5585 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5586 ublas::column( pts, i ) = boost::get<1>( *itL );
5589 __geopc->update( pts );
5590 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5592 __pc->update( pts );
5593 __ctx->update( __c, __pc );
5596 curlx_type __curlx( this->curlx( *__ctx ) );
5599 itL=it->second.begin();
5603 itL=it->second.begin();
5605 for ( uint k=0; k<nbPtsElt; ++k,++itL )
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 );
5613 else if ( nDim == 2 )
5615 itL=it->second.begin();
5617 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5618 v[boost::get<0>( *itL )]( 0,0 ) = __curlx( 0,0,k );
5624 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5625 template<
typename Y,
typename Cont>
5627 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlyInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const
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;
5635 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5642 gm_ptrtype __gm = this->functionSpace()->gm();
5645 if ( it==it_end )
return;
5649 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5650 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5666 for ( ; it!=it_end; ++it )
5668 nbPtsElt = it->second.size();
5671 itL=it->second.begin();
5672 itL_end=it->second.end();
5675 pts= matrix_node_type( nbCoord, nbPtsElt );
5677 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5678 ublas::column( pts, i ) = boost::get<1>( *itL );
5681 __geopc->update( pts );
5682 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5684 __pc->update( pts );
5685 __ctx->update( __c, __pc );
5688 curly_type __curly( this->curly( *__ctx ) );
5691 itL=it->second.begin();
5695 itL=it->second.begin();
5697 for ( uint k=0; k<nbPtsElt; ++k,++itL )
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 );
5705 else if ( nDim == 2 )
5707 itL=it->second.begin();
5709 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5710 v[boost::get<0>( *itL )]( 0,0 ) = __curly( 0,0,k );
5716 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5717 template<
typename Y,
typename Cont>
5719 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::curlzInterpolate( matrix_node_type __ptsReal, comp_curl_array_type& v )
const
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;
5727 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5734 gm_ptrtype __gm = this->functionSpace()->gm();
5737 if ( it==it_end )
return;
5741 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5742 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5758 for ( ; it!=it_end; ++it )
5760 nbPtsElt = it->second.size();
5763 itL=it->second.begin();
5764 itL_end=it->second.end();
5767 pts= matrix_node_type( nbCoord, nbPtsElt );
5769 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5770 ublas::column( pts, i ) = boost::get<1>( *itL );
5773 __geopc->update( pts );
5774 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5776 __pc->update( pts );
5777 __ctx->update( __c, __pc );
5780 curlz_type __curlz( this->curlz( *__ctx ) );
5783 itL=it->second.begin();
5787 itL=it->second.begin();
5789 for ( uint k=0; k<nbPtsElt; ++k,++itL )
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 );
5797 else if ( nDim == 2 )
5799 itL=it->second.begin();
5801 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5802 v[boost::get<0>( *itL )]( 0,0 ) = __curlz( 0,0,k );
5808 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5809 template<
typename Y,
typename Cont>
5810 template<
typename ContextType>
5813 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::d_(
int N, ContextType
const & context, id_array_type& v )
const
5815 if ( !this->areGlobalValuesUpdated() )
5816 this->updateGlobalValues();
5818 for (
int i = 0; i < basis_type::nDof; ++i )
5820 const int ncdof = is_product?nComponents1:1;
5822 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
5829 ( i )( c1 )( ldof )( gdof )
5830 ( this->size() )( this->localSize() )
5831 ( this->firstLocalIndex() )( this->lastLocalIndex() )
5832 .error(
"FunctionSpace::Element invalid access index" );
5834 value_type v_ = this->globalValue( gdof );
5836 for (
size_type q = 0; q < context.xRefs().size2(); ++q )
5838 for (
typename array_type::index i = 0; i < nComponents1; ++i )
5840 v[q]( i,0 ) += v_*context.d( ldof, i, N, q );
5848 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5849 template<
typename Y,
typename Cont>
5851 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dxInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const
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;
5859 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5866 gm_ptrtype __gm = this->functionSpace()->gm();
5869 if ( it==it_end )
return;
5873 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5874 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5890 for ( ; it!=it_end; ++it )
5892 nbPtsElt = it->second.size();
5895 itL=it->second.begin();
5896 itL_end=it->second.end();
5899 pts= matrix_node_type( nbCoord, nbPtsElt );
5901 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5902 ublas::column( pts, i ) = boost::get<1>( *itL );
5905 __geopc->update( pts );
5906 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5908 __pc->update( pts );
5909 __ctx->update( __c, __pc );
5912 dx_type __dx( this->
dx( *__ctx ) );
5915 itL=it->second.begin();
5917 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5918 for (
typename array_type::index i = 0; i < nComponents1; ++i )
5920 v[boost::get<0>( *itL )]( i,0 ) = __dx( i,0,k );
5926 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
5927 template<
typename Y,
typename Cont>
5929 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dyInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const
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;
5937 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
5944 gm_ptrtype __gm = this->functionSpace()->gm();
5947 if ( it==it_end )
return;
5951 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
5952 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
5968 for ( ; it!=it_end; ++it )
5970 nbPtsElt = it->second.size();
5973 itL=it->second.begin();
5974 itL_end=it->second.end();
5977 pts= matrix_node_type( nbCoord, nbPtsElt );
5979 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
5980 ublas::column( pts, i ) = boost::get<1>( *itL );
5983 __geopc->update( pts );
5984 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
5986 __pc->update( pts );
5987 __ctx->update( __c, __pc );
5990 dy_type __dy( this->
dy( *__ctx ) );
5993 itL=it->second.begin();
5995 for ( uint k=0; k<nbPtsElt; ++k,++itL )
5996 for (
typename array_type::index i = 0; i < nComponents1; ++i )
5998 v[boost::get<0>( *itL )]( i,0 ) = __dy( i,0,k );
6004 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
6005 template<
typename Y,
typename Cont>
6007 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::dzInterpolate( matrix_node_type __ptsReal, id_array_type& v )
const
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;
6015 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
6022 gm_ptrtype __gm = this->functionSpace()->gm();
6025 if ( it==it_end )
return;
6029 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
6030 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
6046 for ( ; it!=it_end; ++it )
6048 nbPtsElt = it->second.size();
6051 itL=it->second.begin();
6052 itL_end=it->second.end();
6055 pts= matrix_node_type( nbCoord, nbPtsElt );
6057 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
6058 ublas::column( pts, i ) = boost::get<1>( *itL );
6061 __geopc->update( pts );
6062 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
6064 __pc->update( pts );
6065 __ctx->update( __c, __pc );
6067 dz_type __dz( this->
dz( *__ctx ) );
6070 itL=it->second.begin();
6072 for ( uint k=0; k<nbPtsElt; ++k,++itL )
6073 for (
typename array_type::index i = 0; i < nComponents1; ++i )
6075 v[boost::get<0>( *itL )]( i,0 ) = __dz( i,0,k );
6081 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
6082 template<
typename Y,
typename Cont>
6083 template<
typename ContextType>
6086 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hess_( ContextType
const & context, hess_array_type& v )
const
6088 hess_( context, v, mpl::int_<rank>() );
6090 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
6091 template<
typename Y,
typename Cont>
6092 template<
typename ContextType>
6095 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hess_( ContextType
const & context, hess_array_type& v, mpl::int_<0> )
const
6097 if ( !this->areGlobalValuesUpdated() )
6098 this->updateGlobalValues();
6101 for (
int i = 0; i < basis_type::nDof; ++i )
6103 const int ncdof = is_product?nComponents1:1;
6105 for (
int c1 = 0; c1 < ncdof; ++c1 )
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() )
6112 ( i )( c1 )( ldof )( gdof )
6113 ( this->size() )( this->localSize() )
6114 ( this->firstLocalIndex() )( this->lastLocalIndex() )
6115 .error(
"FunctionSpace::Element invalid access index" );
6117 value_type v_ = this->globalValue( gdof );
6119 for (
size_type q = 0; q < context.xRefs().size2(); ++q )
6121 for (
int i = 0; i < nRealDim; ++i )
6122 for (
int j = 0; j < nRealDim; ++j )
6124 v[q]( i,j ) += v_*context.hess( ldof, i, j, q );
6133 template<
typename A0,
typename A1,
typename A2,
typename A3,
typename A4>
6134 template<
typename Y,
typename Cont>
6136 FunctionSpace<A0, A1, A2, A3, A4>::Element<Y,Cont>::hessInterpolate( matrix_node_type __ptsReal, hess_array_type& v )
const
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;
6144 localization_ptrtype __loc = this->functionSpace()->mesh()->tool_localization();
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;
6151 gm_ptrtype __gm = this->functionSpace()->gm();
6154 if ( it==it_end )
return;
6158 uint nbCoord=boost::get<1>( *( it->second.begin() ) ).size();
6159 matrix_node_type pts( nbCoord, nbPtsElt );
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 ),
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(),
6175 for ( ; it!=it_end; ++it )
6177 nbPtsElt = it->second.size();
6180 itL=it->second.begin();
6181 itL_end=it->second.end();
6184 pts= matrix_node_type( nbCoord, nbPtsElt );
6186 for (
size_type i=0; i<nbPtsElt; ++i,++itL )
6187 ublas::column( pts, i ) = boost::get<1>( *itL );
6190 __geopc->update( pts );
6191 __c->update( this->functionSpace()->mesh()->element( it->first ), __geopc );
6193 __pc->update( pts );
6194 __ctx->update( __c, __pc );
6197 hess_type __hess( this->hess( *__ctx ) );
6200 itL=it->second.begin();
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 );
6210 template<
typename T,
int M,
int N>
6212 operator<<( std::ostream& os, detail::ID<T,M,N>
const& id )
6214 const size_type* shape =
id.M_id.shape();
6216 for (
size_type i = 0; i < shape[0]; ++i )
6218 for (
size_type j = 0; j < shape[1]; ++j )
6220 for (
size_type k = 0; k < shape[2]; ++k )
6222 os << id( i, j, k ) <<
' ';
6241 template <
typename ElementType>
6245 FEELPP_ASSERT( v1.functionSpace() == v2.functionSpace() ).error(
"incompatible function spaces" );
6247 typedef typename type_traits<typename ElementType::value_type>::real_type real_type;
6249 ElementType _t( v1.functionSpace() );
6254 _t.operator()( start+i ) = v1.operator()( start + i )* v2.operator()( start + i );
6266 template <
typename ElementType>
6270 FEELPP_ASSERT( v1.functionSpace() == v2.functionSpace() ).error(
"incompatible function spaces" );
6272 typedef typename type_traits<typename ElementType::value_type>::real_type real_type;
6274 ElementType _t( v1.functionSpace() );
6279 _t.operator()( start+i ) = v1.operator()( start + i )/v2.operator()( start + i );
6284 template<
typename MeshType>
6286 measurePointElements( boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<MeshType::nOrder,Scalar> > > >& _space ) -> decltype( _space->element() )
6288 auto _fn = _space->element(
"measurePointElements" );
6290 std::vector<bool> ptdone( _space->mesh()->numPoints(), false );
6291 auto elit = _space->mesh()->beginElement();
6292 auto elen = _space->mesh()->endElement();
6294 for ( ; elit != elen; ++ elit )
6296 for (
int p = 0; p < elit->numPoints; ++p )
6298 if ( ptdone[elit->point( p ).id()] == false )
6300 BOOST_FOREACH(
auto pt, elit->point( p ).elements() )
6302 _fn.plus_assign( elit->id(), p, 0, elit->measure() );
6304 ptdone[elit->point( p ).id()] =
true;
6312 template<
typename EltType>
6313 typename fusion::result_of::accumulate<
typename EltType::functionspace_type::functionspace_vector_type,
6315 detail::CreateElementVector<EltType> >::type
6316 subelements( EltType
const& e, std::vector<std::string>
const& n )
6318 return fusion::accumulate( e.functionSpaces(), fusion::vector<>(), detail::CreateElementVector<EltType>( e, n ) );
6321 template<
typename ElementType,
typename CoeffType>
6323 expansion( std::vector<ElementType>
const& b, CoeffType
const& c,
int M = -1 )
6325 auto res = b[0].functionSpace()->element();
6327 if ( M == -1 ) M = c.size() ;
6329 for(
int i = 0; i < M; ++i )
6331 res.add( c[i], b[i] );
6343 template<
int Order,
typename MeshType>
6345 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Continuous>>,Periodicity <NoPeriodicity>>>
6346 Pch( boost::shared_ptr<MeshType> mesh )
6348 return FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Continuous>>, Periodicity <NoPeriodicity>>::New( mesh );
6357 template<
int Order,
typename MeshType>
6359 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Discontinuous>>>>
6360 Pdh( boost::shared_ptr<MeshType> mesh )
6362 return FunctionSpace<MeshType,bases<Lagrange<Order,Scalar,Discontinuous>>>::New( mesh );
6371 template<
int Order,
typename MeshType>
6373 boost::shared_ptr<FunctionSpace<MeshType,bases<OrthonormalPolynomialSet<Order,Scalar>>>>
6374 Odh( boost::shared_ptr<MeshType> mesh )
6376 return FunctionSpace<MeshType,bases<OrthonormalPolynomialSet<Order,Scalar>>>::New( mesh );
6385 template<
int Order,
typename MeshType>
6387 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order,Vectorial>>>>
6388 Pchv( boost::shared_ptr<MeshType> mesh )
6390 return FunctionSpace<MeshType,bases<Lagrange<Order,Vectorial>>>::New( mesh );
6397 template<
int Order,
typename MeshType>
6399 boost::shared_ptr<FunctionSpace<MeshType,bases<Lagrange<Order+1,Vectorial>,Lagrange<Order,Scalar>>>>
6400 THch( boost::shared_ptr<MeshType> mesh )
6402 return FunctionSpace<MeshType,bases<Lagrange<Order+1,Vectorial>,Lagrange<Order,Scalar>>>::New( mesh );
6430 struct version< typename Feel::FunctionSpace<A0,A1,A2,A3,A4>::element_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 );
6439 #define FEELPP_REGISTER_ELEMENT( element_type ) \
6441 namespace serialization { \
6443 struct version<element_type> \
6445 typedef mpl::int_<2> type; \
6446 typedef mpl::integral_c_tag tag; \
6447 BOOST_STATIC_CONSTANT(unsigned int, value = version::type::value); \