24 #ifndef __galerkingraph_H
25 #define __galerkingraph_H 1
27 #define FEELPP_EXPORT_GRAPH 0
28 #if FEELPP_EXPORT_GRAPH
29 #include <feel/feelfilters/exporter.hpp>
32 #include <feel/feelvf/pattern.hpp>
35 #include <feel/feeldiscr/functionspace.hpp>
43 template<
typename BFType,
typename Space1Type>
46 compute_graph3( BFType* bf, boost::shared_ptr<Space1Type>
const& space1,
size_type trial_index,
size_type hints )
51 M_trial_index( trial_index ),
55 template <
typename Space2>
56 void operator()( boost::shared_ptr<Space2>
const& space2 )
const
58 if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
60 if ( M_stencil->isBlockPatternZero( M_test_index,M_trial_index ) )
63 #if !defined(FEELPP_ENABLE_MPI_MODE)
64 const size_type proc_id = M_stencil->testSpace()->mesh()->comm().rank();
65 const size_type n1_dof_on_proc = space2->nLocalDof();
66 const size_type first1_dof_on_proc = space2->dof()->firstDof( proc_id );
67 const size_type last1_dof_on_proc = space2->dof()->lastDof( proc_id );
68 const size_type first2_dof_on_proc = M_space1->dof()->firstDof( proc_id );
69 const size_type last2_dof_on_proc = M_space1->dof()->lastDof( proc_id );
71 const size_type proc_id = M_stencil->testSpace()->worldsComm()[M_test_index].globalRank();
72 const size_type n1_dof_on_proc = space2->nLocalDof();
73 const size_type first1_dof_on_proc = space2->dof()->firstDofGlobalCluster( proc_id );
74 const size_type last1_dof_on_proc = space2->dof()->lastDofGlobalCluster( proc_id );
75 const size_type first2_dof_on_proc = M_space1->dof()->firstDofGlobalCluster( proc_id );
76 const size_type last2_dof_on_proc = M_space1->dof()->lastDofGlobalCluster( proc_id );
79 typename BFType::graph_ptrtype zerograph(
new typename BFType::graph_type( space2->dof(), M_space1->dof() ) );
81 M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ), M_stencil->trialSpace()->nDofStart( M_trial_index ) , zerograph );
86 auto thestencil = stencil( _test=space2, _trial=M_space1,
87 _pattern=M_stencil->blockPattern( M_test_index,M_trial_index ),
88 _pattern_block=M_stencil->blockPattern(),
89 _diag_is_nonzero=
false,
90 _collect_garbage=
false,
92 _range=M_stencil->template subRangeIterator<Space2::basis_type::TAG,Space1Type::basis_type::TAG>
93 (mpl::bool_<BFType::template rangeiteratorType<Space2::basis_type::TAG,Space1Type::basis_type::TAG>::hasnotfindrange_type::value>() )
96 if ( M_stencil->testSpace()->worldComm().globalSize()>1 && M_stencil->testSpace()->hasEntriesForAllSpaces() )
97 M_stencil->mergeGraphMPI( M_test_index, M_trial_index,
98 space2->mapOn(), M_space1->mapOn(),
99 thestencil->graph() );
101 M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
102 M_stencil->trialSpace()->nDofStart( M_trial_index ),
103 thestencil->graph() );
111 mutable BFType* M_stencil;
112 boost::shared_ptr<Space1Type>
const& M_space1;
119 template<
typename BFType,
typename Space1Type>
120 struct compute_graph2
122 compute_graph2( BFType* bf, boost::shared_ptr<Space1Type>
const& space1,
size_type test_index,
size_type hints )
126 M_test_index( test_index ),
131 template <
typename Space2>
132 void operator()( boost::shared_ptr<Space2>
const& space2 )
const
134 if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
136 if ( M_stencil->isBlockPatternZero( M_test_index,M_trial_index ) )
139 #if !defined(FEELPP_ENABLE_MPI_MODE)
140 const size_type proc_id = M_stencil->testSpace()->template mesh<0>()->comm().rank();
141 const size_type n1_dof_on_proc = M_space1->nLocalDof();
142 const size_type first1_dof_on_proc = M_space1->dof()->firstDof( proc_id );
143 const size_type last1_dof_on_proc = M_space1->dof()->lastDof( proc_id );
144 const size_type first2_dof_on_proc = space2->dof()->firstDof( proc_id );
145 const size_type last2_dof_on_proc = space2->dof()->lastDof( proc_id );
147 const size_type proc_id = M_stencil->testSpace()->worldsComm()[M_test_index].globalRank();
148 const size_type n1_dof_on_proc = M_space1->nLocalDof();
149 const size_type first1_dof_on_proc = M_space1->dof()->firstDofGlobalCluster( proc_id );
150 const size_type last1_dof_on_proc = M_space1->dof()->lastDofGlobalCluster( proc_id );
151 const size_type first2_dof_on_proc = space2->dof()->firstDofGlobalCluster( proc_id );
152 const size_type last2_dof_on_proc = space2->dof()->lastDofGlobalCluster( proc_id );
155 typename BFType::graph_ptrtype zerograph(
new typename BFType::graph_type( M_space1->dof(), space2->dof() ) );
157 M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
158 M_stencil->trialSpace()->nDofStart( M_trial_index ),
165 auto thestencil = stencil( _test=M_space1, _trial=space2,
166 _pattern=M_stencil->blockPattern( M_test_index,M_trial_index ),
167 _pattern_block=M_stencil->blockPattern(),
168 _diag_is_nonzero=
false,
169 _collect_garbage=
false,
171 _range=M_stencil->template subRangeIterator<Space1Type::basis_type::TAG,Space2::basis_type::TAG>
172 (mpl::bool_<BFType::template rangeiteratorType<Space1Type::basis_type::TAG,Space2::basis_type::TAG>::hasnotfindrange_type::value>() )
175 if ( M_stencil->testSpace()->worldComm().globalSize()>1 && M_stencil->testSpace()->hasEntriesForAllSpaces() )
176 M_stencil->mergeGraphMPI( M_test_index, M_trial_index,
177 M_space1->mapOn(), space2->mapOn(),
178 thestencil->graph() );
180 M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
181 M_stencil->trialSpace()->nDofStart( M_trial_index ),
182 thestencil->graph() );
190 mutable BFType* M_stencil;
191 boost::shared_ptr<Space1Type>
const& M_space1;
198 template<
typename BFType>
199 struct compute_graph1
201 compute_graph1( BFType* bf,
size_type hints )
208 template <
typename Space1>
209 void operator()( boost::shared_ptr<Space1>
const& space1 )
const
211 fusion::for_each( M_stencil->trialSpace()->functionSpaces(),
212 compute_graph2<BFType,Space1>( M_stencil, space1, M_test_index, M_hints ) );
215 mutable BFType* M_stencil;
225 template <
int I,
int J,
typename IteratorRange>
226 fusion::pair< fusion::pair<mpl::int_<I>,mpl::int_<J> >,IteratorRange >
227 stencilRange( IteratorRange
const& r)
229 return fusion::make_pair< fusion::pair<mpl::int_<I>,mpl::int_<J> > >( r);
233 struct stencilRangeMapTypeBase {};
235 struct stencilRangeMap0Type
237 public stencilRangeMapTypeBase,
240 typedef fusion::map<> super_type;
242 stencilRangeMap0Type()
247 static bool isNullRange() {
return true; }
250 template <
typename ThePair1Type>
251 struct stencilRangeMap1Type
253 public stencilRangeMapTypeBase,
254 public fusion::map< fusion::pair< typename ThePair1Type::first_type, typename ThePair1Type::second_type > >
256 typedef typename ThePair1Type::first_type key1_type;
257 typedef fusion::map< fusion::pair< key1_type, typename ThePair1Type::second_type > > super_type;
259 stencilRangeMap1Type( ThePair1Type
const& p )
264 static bool isNullRange() {
return false; }
267 template <
typename ThePair1Type,
typename ThePair2Type>
268 struct stencilRangeMap2Type
270 public stencilRangeMapTypeBase,
271 public fusion::map< fusion::pair< typename ThePair1Type::first_type, typename ThePair1Type::second_type >,
272 fusion::pair< typename ThePair2Type::first_type, typename ThePair2Type::second_type > >
274 typedef typename ThePair1Type::first_type key1_type;
275 typedef typename ThePair2Type::first_type key2_type;
276 typedef fusion::map< fusion::pair< key1_type, typename ThePair1Type::second_type >,
277 fusion::pair< key2_type, typename ThePair2Type::second_type > > super_type;
279 stencilRangeMap2Type( ThePair1Type
const& p1, ThePair2Type
const& p2 )
284 static bool isNullRange() {
return false; }
287 stencilRangeMap0Type stencilRangeMap();
289 template <
typename ThePair1Type>
290 stencilRangeMap1Type<ThePair1Type>
291 stencilRangeMap( ThePair1Type
const& p)
293 return stencilRangeMap1Type<ThePair1Type>( p.second );
296 template <
typename ThePair1Type,
typename ThePair2Type>
297 stencilRangeMap2Type<ThePair1Type,ThePair2Type>
298 stencilRangeMap( ThePair1Type
const& p1, ThePair2Type
const& p2)
300 return stencilRangeMap2Type<ThePair1Type,ThePair2Type>( p1.second, p2.second );
305 template<
typename X1,
typename X2,
typename RangeIteratorTestType = stencilRangeMap0Type >
309 typedef X1 test_space_ptrtype;
310 typedef X2 trial_space_ptrtype;
311 typedef typename X1::element_type test_space_type;
312 typedef typename X2::element_type trial_space_type;
313 typedef GraphCSR graph_type;
314 typedef boost::shared_ptr<graph_type> graph_ptrtype;
315 typedef Stencil<X1,X2,RangeIteratorTestType> self_type;
317 typedef RangeIteratorTestType rangeiterator_test_type;
319 Stencil( test_space_ptrtype Xh, trial_space_ptrtype Yh,
321 BlocksStencilPattern block_pattern=BlocksStencilPattern(1,1,Pattern::HAS_NO_BLOCK_PATTERN),
322 bool diag_is_nonzero=
false,
324 rangeiterator_test_type r=rangeiterator_test_type())
328 #if !defined(FEELPP_ENABLE_MPI_MODE)
329 M_graph( new graph_type( Xh->nLocalDof(),
330 Xh->nDofStart(), Xh->nDofStart()+ Xh->nLocalDof()-1,
331 Yh->nDofStart(), Yh->nDofStart()+ Yh->nLocalDof()-1 ) ),
333 M_graph( new graph_type( Xh->dof(),Yh->dof() ) ),
335 M_block_pattern( block_pattern ),
336 M_rangeIteratorTest( r )
339 uint16_type nbSubSpace1 = _M_X1->nSubFunctionSpace();
340 uint16_type nbSubSpace2 = _M_X2->nSubFunctionSpace();
342 if ( this->isBlockPatternNoPattern( 0,0 ) )
344 M_block_pattern = BlocksStencilPattern(nbSubSpace1,nbSubSpace2,graph_hints);
347 M_graph = this->computeGraph( graph_hints );
349 if ( diag_is_nonzero && _M_X1->nLocalDofWithoutGhost()>0 && _M_X2->nLocalDofWithoutGhost()>0 ) M_graph->addMissingZeroEntriesDiagonal();
351 if ( close ) M_graph->close();
354 Stencil( test_space_ptrtype Xh, trial_space_ptrtype Yh,
size_type graph_hints, graph_ptrtype g, rangeiterator_test_type r=rangeiterator_test_type() )
359 M_block_pattern(Xh->nSubFunctionSpace(),Yh->nSubFunctionSpace(),
size_type( graph_hints )),
360 M_rangeIteratorTest( r )
364 BlocksStencilPattern blockPattern()
const
366 return M_block_pattern;
369 size_type blockPattern( uint16_type i,uint16_type j )
const
371 return M_block_pattern(i,j);
374 bool isBlockPatternNoPattern( uint16_type i,uint16_type j )
const
377 return ctx.test( HAS_NO_BLOCK_PATTERN );
380 bool isBlockPatternZero( uint16_type i,uint16_type j )
const
383 return ctx.test( ZERO );
386 graph_ptrtype computeGraph(
size_type hints );
388 graph_ptrtype computeGraph(
size_type hints, mpl::bool_<true> );
389 graph_ptrtype computeGraph(
size_type hints, mpl::bool_<false> );
390 graph_ptrtype computeGraph(
size_type hints, mpl::bool_<true>, mpl::bool_<true> );
391 graph_ptrtype computeGraph(
size_type hints, mpl::bool_<false>, mpl::bool_<true> );
392 graph_ptrtype computeGraph(
size_type hints, mpl::bool_<true>, mpl::bool_<false> );
394 graph_ptrtype computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true> );
395 graph_ptrtype computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<false> );
396 graph_ptrtype computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true>, mpl::bool_<true> );
397 graph_ptrtype computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<false>, mpl::bool_<true> );
398 graph_ptrtype computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true>, mpl::bool_<false> );
401 void mergeGraph(
int row,
int col, graph_ptrtype g );
403 DataMap
const& mapOnTest, DataMap
const& mapOnTrial,
406 test_space_ptrtype testSpace()
const
410 trial_space_ptrtype trialSpace()
const
414 graph_ptrtype graph()
const
418 graph_ptrtype graph()
423 std::set<size_type> trialElementId(
size_type test_eid, mpl::int_<0> )
425 std::set<size_type> idsFind;
426 const bool test_related_to_trial = _M_X1->mesh()->isSubMeshFrom( _M_X2->mesh() );
427 const bool trial_related_to_test = _M_X2->mesh()->isSubMeshFrom( _M_X1->mesh() );
428 if ( test_related_to_trial )
430 const size_type domain_eid = _M_X1->mesh()->subMeshToMesh( test_eid );
431 DVLOG(2) <<
"[test_related_to_trial] test element id: " << test_eid <<
" trial element id : " << domain_eid <<
"\n";
434 else if( trial_related_to_test )
436 const size_type domain_eid = _M_X2->mesh()->meshToSubMesh( test_eid );
437 DVLOG(2) <<
"[trial_related_to_test] test element id: " << test_eid <<
" trial element id : " << domain_eid <<
"\n";
442 idsFind.insert(test_eid);
447 std::set<size_type> trialElementId(
size_type test_eid, mpl::int_<1> )
449 std::set<size_type> idsFind;
450 const bool test_related_to_trial = _M_X1->mesh()->isSubMeshFrom( _M_X2->mesh() );
451 const bool trial_related_to_test = _M_X2->mesh()->isSubMeshFrom( _M_X1->mesh() );
452 if ( test_related_to_trial )
454 const size_type domain_eid = _M_X2->mesh()->face(_M_X1->mesh()->subMeshToMesh( test_eid )).element0().id();
455 DVLOG(2) <<
"[test_related_to_trial<1>] test element id: " << test_eid <<
" trial element id : " << domain_eid <<
"\n";
458 else if( trial_related_to_test )
460 auto const& eltTest = _M_X1->mesh()->element(test_eid);
461 for (uint16_type f=0;f< _M_X1->mesh()->numLocalFaces();++f)
463 const size_type idFind = _M_X2->mesh()->meshToSubMesh( eltTest.face(f).id() );
466 DVLOG(2) <<
"[trial_related_to_test<1>] test element id: " << test_eid <<
" idsFind.size() "<< idsFind.size() <<
"\n";
470 CHECK (
false ) <<
"[trial_related_to_test<1>] : test and trial mesh can not be the same here\n";
478 template <
int I,
int J>
479 struct rangeiteratorType
481 typedef typename fusion::result_of::find<rangeiterator_test_type,fusion::pair<mpl::int_<I>,mpl::int_<J> > >::type resultfindrange_it_type;
482 typedef typename fusion::result_of::value_of<resultfindrange_it_type>::type resultfindrange_type;
484 typedef typename boost::is_same<resultfindrange_it_type, typename fusion::result_of::end<rangeiterator_test_type>::type> hasnotfindrange_type;
485 typedef typename boost::tuple<mpl::size_t<MESH_ELEMENTS>,
486 typename MeshTraits<typename test_space_type::mesh_type>::element_const_iterator,
487 typename MeshTraits<typename test_space_type::mesh_type>::element_const_iterator> defaultrange_type;
489 typedef typename mpl::if_< hasnotfindrange_type,
490 mpl::identity< defaultrange_type >,
491 mpl::identity< resultfindrange_type >
495 template <
int I,
int J>
496 typename rangeiteratorType<I,J>::type
497 rangeiterator()
const
499 return rangeiterator<I,J>( mpl::bool_<rangeiteratorType<I,J>::hasnotfindrange_type::value>() );
502 template <
int I,
int J>
503 typename rangeiteratorType<I,J>::defaultrange_type
504 rangeiterator(mpl::bool_<true> )
const
508 template <
int I,
int J>
509 typename rangeiteratorType<I,J>::resultfindrange_type::second_type
510 rangeiterator(mpl::bool_<false> )
const
512 typedef fusion::pair<mpl::int_<I>,mpl::int_<J> > key_type;
513 return fusion::at_key< key_type >(M_rangeIteratorTest);
516 template <
int I,
int J>
518 subRangeIterator( mpl::bool_<true> )
520 return stencilRangeMap0Type();
522 template <
int I,
int J>
523 stencilRangeMap1Type< fusion::pair< fusion::pair<mpl::int_<0>,mpl::int_<0> >,
typename rangeiteratorType<I,J>::resultfindrange_type::second_type > >
524 subRangeIterator( mpl::bool_<false> )
526 return stencilRangeMap( stencilRange<0,0>( rangeiterator<I,J>(mpl::bool_<false>()) ) );
531 test_space_ptrtype _M_X1;
532 trial_space_ptrtype _M_X2;
533 graph_ptrtype M_graph;
534 BlocksStencilPattern M_block_pattern;
535 rangeiterator_test_type M_rangeIteratorTest;
539 template<
typename Args>
540 struct compute_stencil_type
542 typedef typename remove_pointer_const_reference_type<Args,tag::test>::type _test_type;
543 typedef typename remove_pointer_const_reference_type<Args,tag::trial>::type _trial_type;
544 typedef typename remove_pointer_const_reference_default_type<Args,tag::range, stencilRangeMap0Type >::type _range_type;
545 typedef Stencil<_test_type, _trial_type, _range_type> type;
546 typedef boost::shared_ptr<type> ptrtype;
552 class StencilManagerImpl:
553 public std::map<boost::tuple<boost::shared_ptr<FunctionSpaceBase>,
554 boost::shared_ptr<FunctionSpaceBase>,
556 std::vector<size_type>,
557 bool >, boost::shared_ptr<GraphCSR> >,
558 public boost::noncopyable
561 typedef boost::shared_ptr<GraphCSR> graph_ptrtype;
562 typedef boost::tuple<boost::shared_ptr<FunctionSpaceBase>,
563 boost::shared_ptr<FunctionSpaceBase>,
565 std::vector<size_type>,
567 typedef std::map<key_type, graph_ptrtype> graph_manager_type;
578 void stencilManagerAdd(StencilManagerImpl::key_type
const& key,StencilManagerImpl::graph_ptrtype graph);
582 extern BlocksStencilPattern default_block_pattern;
584 BOOST_PARAMETER_FUNCTION(
585 (
typename detail::compute_stencil_type<Args>::ptrtype ),
589 ( test, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
590 ( trial, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
593 ( pattern, *( boost::is_integral<mpl::_> ), Pattern::COUPLED )
594 ( pattern_block, *, default_block_pattern )
595 ( diag_is_nonzero, *( boost::is_integral<mpl::_> ),
false )
596 ( collect_garbage, *( boost::is_integral<mpl::_> ),
true )
597 ( close, *( boost::is_integral<mpl::_> ),
true )
598 ( range, *( boost::is_convertible<mpl::_, stencilRangeMapTypeBase>) , stencilRangeMap0Type() )
602 if ( collect_garbage )
608 Feel::detail::ignore_unused_variable_warning( args );
609 typedef typename detail::compute_stencil_type<Args>::ptrtype stencil_ptrtype;
610 typedef typename detail::compute_stencil_type<Args>::type stencil_type;
613 auto git =
StencilManager::instance().find( boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero ) );
618 auto s = stencil_ptrtype(
new stencil_type( test, trial, pattern, git->second, range ) );
625 auto git_trans =
StencilManager::instance().find( boost::make_tuple( trial, test, pattern, pattern_block.transpose().getSetOfBlocks(), diag_is_nonzero ) );
629 auto g = git_trans->second->transpose(close);
631 StencilManager::instance().operator[]( boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero ) ) = g;
632 auto s = stencil_ptrtype(
new stencil_type( test, trial, pattern, g, range ) );
640 auto s = stencil_ptrtype(
new stencil_type( test, trial, pattern, pattern_block, diag_is_nonzero, close, range ) );
641 if ( range.isNullRange() )
StencilManager::instance().operator[]( boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero ) ) = s->graph();
650 template<
typename B
idirectionalIterator>
653 sortSparsityRow (
const BidirectionalIterator begin,
654 BidirectionalIterator middle,
655 const BidirectionalIterator end )
657 if ( ( begin == middle ) || ( middle == end ) )
return;
659 assert ( std::distance ( begin, middle ) > 0 );
660 assert ( std::distance ( middle, end ) > 0 );
661 FEELPP_ASSERT( std::unique ( begin, middle ) == middle )
662 ( *begin )( *middle ).error(
"duplicate dof(begin,middle)" );
663 FEELPP_ASSERT ( std::unique ( middle, end ) == end )
664 ( *begin )( *middle ).error(
"duplicate dof (middle,end)" );
666 while ( middle != end )
668 BidirectionalIterator
673 while ( !( *a < *b ) )
675 std::swap ( *a, *b );
677 if ( a == begin )
break;
692 BidirectionalIterator
696 for ( ++first; first != end; prev=first, ++first )
697 if ( *first < *prev )
703 assert ( std::unique ( begin, end ) == end );
707 template<
typename X1,
typename X2,
typename RangeItTestType>
709 Stencil<X1,X2,RangeItTestType>::mergeGraph(
int row,
int col, graph_ptrtype g )
712 DVLOG(2) <<
"[merge graph] for composite bilinear form\n";
713 DVLOG(2) <<
"[mergeGraph] row = " << row <<
"\n";
714 DVLOG(2) <<
"[mergeGraph] col = " << col <<
"\n";
720 DVLOG(2) <<
"[merge graph] nothing yet in store, copy graph\n";
724 typename graph_type::const_iterator it = g->begin();
725 typename graph_type::const_iterator en = g->end();
727 for ( ; it != en; ++it )
729 std::vector<size_type>
const& ivec = boost::get<2>( it->second );
731 for (
int i = 0; i < ivec.size(); ++i )
733 DVLOG(2) <<
"[mergeGraph] ivec[" << i <<
"] = " << ivec[i] <<
"\n";
743 DVLOG(2) <<
"[merge graph] already something in store\n";
744 typename graph_type::const_iterator it = g->begin();
745 typename graph_type::const_iterator en = g->end();
747 for ( ; it != en; ++it )
749 int theglobalrow = row+it->first;
752 if ( this->testSpace()->worldComm().globalSize()>1 )
753 thelocalrow = boost::get<1>( it->second );
756 thelocalrow = row + boost::get<1>( it->second );
759 std::set<size_type>& row1_entries = M_graph->row( theglobalrow ).template get<2>();
760 std::set<size_type>
const& row2_entries = boost::get<2>( it->second );
762 DVLOG(2) <<
"[mergeGraph] adding information to global row [" << theglobalrow <<
"], localrow=" << thelocalrow <<
"\n";
763 M_graph->row( theglobalrow ).template get<1>() = thelocalrow;
764 M_graph->row( theglobalrow ).template get<0>() = this->testSpace()->worldComm().mapLocalRankToGlobalRank()[it->second.get<0>()];
766 if ( !row2_entries.empty() )
769 if ( row1_entries.empty() )
774 if ( col==0 ) row1_entries = row2_entries;
778 for (
auto itcol = row2_entries.begin(), encol = row2_entries.end() ; itcol!=encol; ++itcol ) row1_entries.insert( *itcol+col );
786 auto itg = boost::prior( row1_entries.end() );
788 std::for_each( row2_entries.begin(), row2_entries.end(),[&](
size_type o )
790 itg = row1_entries.insert( itg, o+col );
797 DVLOG(2) <<
" -- merge_graph (" << row <<
"," << col <<
") in " << tim.elapsed() <<
"\n";
798 DVLOG(2) <<
"merge graph for composite bilinear form done\n";
801 template<
typename X1,
typename X2,
typename RangeItTestType>
803 Stencil<X1,X2,RangeItTestType>::mergeGraphMPI(
size_type test_index,
size_type trial_index,
804 DataMap
const& mapOnTest, DataMap
const& mapOnTrial,
807 const int myrank = this->testSpace()->worldComm().globalRank();
809 const size_type globalDofRowStart = this->testSpace()->dof()->firstDofGlobalCluster() + this->testSpace()->nLocalDofWithoutGhostOnProcStart( myrank, test_index );
810 const size_type globalDofColStart = this->trialSpace()->dof()->firstDofGlobalCluster() + this->trialSpace()->nLocalDofWithoutGhostOnProcStart( myrank, trial_index );
811 const size_type locdofStart = this->testSpace()->nLocalDofWithGhostOnProcStart( myrank, test_index );
813 typename graph_type::const_iterator it = g->begin();
814 typename graph_type::const_iterator en = g->end();
815 for ( ; it != en; ++it )
817 size_type theglobalrow = globalDofRowStart + ( it->first - mapOnTest.firstDofGlobalCluster() );
818 const size_type thelocalrow = locdofStart + it->second.get<1>();
820 if (it->second.get<0>()!=g->worldComm().globalRank() )
822 const int proc = it->second.get<0>();
823 const size_type realrowStart = this->testSpace()->dof()->firstDofGlobalCluster(proc)
824 + this->testSpace()->nLocalDofWithoutGhostOnProcStart(proc, test_index );
825 theglobalrow = realrowStart+(it->first-mapOnTest.firstDofGlobalCluster(proc));
828 std::set<size_type>& row1_entries = M_graph->row( theglobalrow ).template get<2>();
829 std::set<size_type>
const& row2_entries = boost::get<2>( it->second );
831 DVLOG(2) <<
"[mergeGraph] adding information to global row [" << theglobalrow <<
"], localrow=" << thelocalrow <<
"\n";
832 M_graph->row( theglobalrow ).template get<1>() = thelocalrow;
833 M_graph->row( theglobalrow ).template get<0>() = this->testSpace()->worldComm().mapLocalRankToGlobalRank()[it->second.get<0>()];
835 if ( !row2_entries.empty() )
837 for (
auto itcol = row2_entries.begin(), encol = row2_entries.end() ; itcol!=encol; ++itcol )
839 if (mapOnTrial.dofGlobalClusterIsOnProc(*itcol))
841 const size_type dofcol = globalDofColStart + (*itcol-mapOnTrial.firstDofGlobalCluster());
842 row1_entries.insert( dofcol );
846 const int realproc = mapOnTrial.procOnGlobalCluster(*itcol);
847 const size_type realcolStart = this->trialSpace()->dof()->firstDofGlobalCluster(realproc)
848 + this->trialSpace()->nLocalDofWithoutGhostOnProcStart( realproc, trial_index );
849 const size_type dofcol = realcolStart + (*itcol - mapOnTrial.firstDofGlobalCluster(realproc));
850 row1_entries.insert(dofcol);
859 template<
typename X1,
typename X2,
typename RangeItTestType>
860 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
861 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints )
863 VLOG(2) <<
"computeGraph: deciding whether the mesh are related to optimize the stencil\n";
866 if ( _M_X1->template mesh<0>()->isRelatedTo( _M_X2->template mesh<0>() ) )
868 VLOG(2) <<
"computeGraph: meshes are related\n";
869 return this->computeGraph( hints, mpl::bool_<mpl::and_< mpl::bool_< ( test_space_type::nSpaces == 1 )>,
870 mpl::bool_< ( trial_space_type::nSpaces == 1 )> >::type::value >() );
874 VLOG(2) <<
"computeGraph: meshes are not related\n";
875 return this->computeGraphInCaseOfInterpolate( hints, mpl::bool_<mpl::and_< mpl::bool_< ( test_space_type::nSpaces == 1 )>,
876 mpl::bool_< ( trial_space_type::nSpaces == 1 )> >::type::value >() );
881 template<
typename X1,
typename X2,
typename RangeItTestType>
882 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
883 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<false> )
886 DVLOG(2) <<
"compute graph for composite bilinear form with interpolation\n";
888 auto graph = computeGraph( hints, mpl::bool_< ( test_space_type::nSpaces > 1 )>(), mpl::bool_< ( trial_space_type::nSpaces > 1 )>() );
890 DVLOG(2) <<
"closing graph for composite bilinear form with interpolation done in " << t.elapsed() <<
"s\n";
893 DVLOG(2) <<
"compute graph for composite bilinear form done in " << t.elapsed() <<
"s\n";
898 template<
typename X1,
typename X2,
typename RangeItTestType>
899 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
900 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<true>, mpl::bool_<true> )
902 fusion::for_each( _M_X1->functionSpaces(),
903 detail::compute_graph1<self_type>(
this, hints ) );
907 template<
typename X1,
typename X2,
typename RangeItTestType>
908 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
909 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<true>, mpl::bool_<false> )
911 fusion::for_each( _M_X1->functionSpaces(),
912 detail::compute_graph3<self_type,trial_space_type>(
this, _M_X2, 0, hints ) );
916 template<
typename X1,
typename X2,
typename RangeItTestType>
917 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
918 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<false>, mpl::bool_<true> )
920 fusion::for_each( _M_X2->functionSpaces(),
921 detail::compute_graph2<self_type,test_space_type>(
this, _M_X1, 0, hints ) );
925 template<
typename X1,
typename X2,
typename RangeItTestType>
926 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
927 Stencil<X1,X2,RangeItTestType>::computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<false> )
930 DVLOG(2) <<
"compute graph for composite bilinear form with interpolation\n";
931 auto graph = computeGraphInCaseOfInterpolate( hints,
932 mpl::bool_< ( test_space_type::nSpaces > 1 )>(),
933 mpl::bool_< ( trial_space_type::nSpaces > 1 )>() );
935 DVLOG(2) <<
"closing graph for composite bilinear form with interpolation done in " << t.elapsed() <<
"s\n";
939 template<
typename X1,
typename X2,
typename RangeItTestType>
940 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
941 Stencil<X1,X2,RangeItTestType>::computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true>, mpl::bool_<true> )
943 fusion::for_each( _M_X1->functionSpaces(),
944 detail::compute_graph1<self_type>(
this, hints ) );
948 template<
typename X1,
typename X2,
typename RangeItTestType>
949 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
950 Stencil<X1,X2,RangeItTestType>::computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true>, mpl::bool_<false> )
952 fusion::for_each( _M_X1->functionSpaces(),
953 detail::compute_graph3<self_type,trial_space_type>(
this, _M_X2, 0, hints ) );
957 template<
typename X1,
typename X2,
typename RangeItTestType>
958 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
959 Stencil<X1,X2,RangeItTestType>::computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<false>, mpl::bool_<true> )
961 fusion::for_each( _M_X2->functionSpaces(),
962 detail::compute_graph2<self_type,test_space_type>(
this, _M_X1, 0, hints ) );
970 template<
typename X1,
typename X2,
typename RangeItTestType>
971 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
972 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<true> )
979 const size_type proc_id = _M_X1->mesh()->comm().rank();
980 const size_type n1_dof_on_proc = _M_X1->nLocalDof();
982 const size_type first1_dof_on_proc = _M_X1->dof()->firstDof( proc_id );
983 const size_type last1_dof_on_proc = _M_X1->dof()->lastDof( proc_id );
984 const size_type first2_dof_on_proc = _M_X2->dof()->firstDof( proc_id );
985 const size_type last2_dof_on_proc = _M_X2->dof()->lastDof( proc_id );
986 graph_ptrtype sparsity_graph(
new graph_type( n1_dof_on_proc,
987 first1_dof_on_proc, last1_dof_on_proc,
988 first2_dof_on_proc, last2_dof_on_proc ) );
990 typedef typename mesh_type::element_const_iterator mesh_element_const_iterator;
992 mesh_element_const_iterator elem_it = _M_X1->mesh()->beginElementWithProcessId( proc_id );
993 const mesh_element_const_iterator elem_en = _M_X1->mesh()->endElementWithProcessId( proc_id );
1000 DVLOG(2) <<
"[computeGraph] test : " << ( graph.test ( Pattern::COUPLED ) || graph.test ( Pattern::EXTENDED ) ) <<
"\n";
1001 DVLOG(2) <<
"[computeGraph] : graph.test ( Pattern::COUPLED )=" << graph.test ( Pattern::COUPLED ) <<
"\n";
1002 DVLOG(2) <<
"[computeGraph] : graph.test ( Pattern::EXTENDED)=" << graph.test ( Pattern::EXTENDED ) <<
"\n";
1005 if ( graph.test ( Pattern::COUPLED ) ||
1006 graph.test ( Pattern::EXTENDED ) )
1011 DVLOG(2) <<
"[computeGraph] test (Pattern::COUPLED || Pattern::EXTENDED) ok\n";
1012 std::vector<size_type>
1018 for ( ; elem_it != elem_en; ++elem_it )
1020 #if !defined(NDEBUG)
1021 DVLOG(2) <<
"[Stencil::computePatter] element " << elem_it->id() <<
" on proc " << elem_it->processId() <<
"\n";
1023 mesh_element_type
const& elem = *elem_it;
1026 element_dof1 = _M_X1->dof()->getIndices( elem.id() );
1027 element_dof2 = _M_X2->dof()->getIndices( elem.id() );
1031 std::sort( element_dof1.begin(), element_dof1.end() );
1032 std::sort( element_dof2.begin(), element_dof2.end() );
1034 const uint16_type n1_dof_on_element = element_dof1.size();
1035 const uint16_type n2_dof_on_element = element_dof2.size();
1037 for (
size_type i=0; i<n1_dof_on_element; i++ )
1051 FEELPP_ASSERT ( ig1 >= first1_dof_on_proc )( ig1 )( first1_dof_on_proc ).error (
"invalid dof index" );
1052 FEELPP_ASSERT ( ( ig1 - first1_dof_on_proc ) < sparsity_graph->size() )
1053 ( ig1 )( first1_dof_on_proc )( sparsity_graph->size() ).error(
"invalid dof index" );
1055 graph_type::row_type& row = sparsity_graph->row( ig1 );
1056 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1059 DVLOG(2) <<
"work with row " << ig1 <<
" local index " << ig1 - first1_dof_on_proc <<
"\n";
1063 if ( row.get<2>().empty() )
1065 row.get<2>() = element_dof2;
1076 std::vector<size_type>::iterator
1077 low = std::lower_bound ( row.get<2>().begin(), row.get<2>().end(), element_dof2.front() ),
1078 high = std::upper_bound ( low, row.get<2>().end(), element_dof2.back() );
1080 for (
size_type j=0; j<n2_dof_on_element; j++ )
1086 std::pair<std::vector<size_type>::iterator,
1087 std::vector<size_type>::iterator>
1088 pos = std::equal_range ( low, high, jg );
1091 if ( pos.first == pos.second )
1092 dof_to_add.push_back( jg );
1100 std::pair<std::vector<size_type>::iterator,
1101 std::vector<size_type>::iterator>
1102 pos = std::equal_range ( row.get<2>().begin(), row.get<2>().end(), jg );
1105 if ( pos.first == pos.second )
1106 dof_to_add.push_back( jg );
1112 if ( !dof_to_add.empty() )
1114 const size_type old_size = row.get<2>().size();
1116 row.get<2>().insert ( row.get<2>().end(),
1121 sortSparsityRow ( row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end() );
1128 for ( uint16_type ms=0; ms < elem.nNeighbors(); ms++ )
1130 mesh_element_type
const* neighbor = NULL;
1131 size_type neighbor_id = elem.neighbor( ms ).first;
1132 size_type neighbor_process_id = elem.neighbor( ms ).second;
1139 VLOG(1) <<
"element id " << elem.id()
1140 <<
", element neighbor id " << neighbor_id
1141 <<
" in proc " << neighbor_process_id <<
"\n";
1143 neighbor = boost::addressof( _M_X1->mesh()->element( neighbor_id,
1144 neighbor_process_id ) );
1147 VLOG(1) <<
"found neighbor of element id " << elem.id()
1148 <<
", element neighbor id " << neighbor->id()
1149 <<
" in proc " << neighbor->processId() <<
"\n";
1152 if ( neighbor_id == neighbor->id() )
1154 neighbor_dof = _M_X2->dof()->getIndices( neighbor->id() );
1156 const size_type n_dof_on_neighbor = neighbor_dof.size();
1159 for (
size_type j=0; j<n_dof_on_neighbor; j++ )
1161 DVLOG(2) <<
"neighbor elem id: " << neighbor->id() <<
" dof " << neighbor_dof[j] <<
"\n";
1164 DVLOG(2) <<
"looking for dof " << ig1 <<
"\n";
1167 std::pair<std::vector<size_type>::iterator,
1168 std::vector<size_type>::iterator>
1169 posig = std::equal_range ( neighbor_dof.begin(), neighbor_dof.end(), ig1 );
1171 std::vector<size_type>::iterator posig = std::find( neighbor_dof.begin(), neighbor_dof.end(), ig1 );
1176 if ( posig != neighbor_dof.end() ||
1177 graph.test ( Pattern::EXTENDED ) )
1181 for (
size_type j=0; j<n_dof_on_neighbor; j++ )
1187 std::pair<std::vector<size_type>::iterator,
1188 std::vector<size_type>::iterator>
1189 pos = std::equal_range ( row.get<2>().begin(), row.get<2>().end(), jg );
1191 std::vector<size_type>::iterator pos = std::find( row.get<2>().begin(), row.get<2>().end(), jg );
1196 if ( pos == row.get<2>().end() )
1198 const size_type old_size = row.get<2>().size();
1199 row.get<2>().push_back ( jg );
1201 sortSparsityRow ( row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end() );
1212 for (
int k = 0; k < row.get<2>().size(); ++k )
1213 VLOG(1) <<
"row[ " << ig1 - first1_dof_on_proc <<
","<< k <<
" ]=" << row.get<2>()[k] <<
"\n";
1225 DVLOG(2) <<
"[computeGraph<true>] before calling close in " << t.elapsed() <<
"s\n";
1227 DVLOG(2) <<
"[computeGraph<true>] done in " << t.elapsed() <<
"s\n";
1228 DVLOG(2) <<
"[computeGraph<true>] done in " << t.elapsed() <<
"s\n";
1229 return sparsity_graph;
1232 template<
typename X1,
typename X2,
typename RangeItTestType>
1233 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
1234 Stencil<X1,X2,RangeItTestType>::computeGraph(
size_type hints, mpl::bool_<true> )
1241 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1242 const size_type proc_id = _M_X1->mesh()->comm().rank();
1243 const size_type n1_dof_on_proc = _M_X1->nLocalDof();
1245 const size_type first1_dof_on_proc = _M_X1->dof()->firstDof( proc_id );
1246 const size_type last1_dof_on_proc = _M_X1->dof()->lastDof( proc_id );
1247 const size_type first2_dof_on_proc = _M_X2->dof()->firstDof( proc_id );
1248 const size_type last2_dof_on_proc = _M_X2->dof()->lastDof( proc_id );
1250 const size_type proc_id = _M_X1->worldsComm()[0].localRank();
1251 const size_type n1_dof_on_proc = _M_X1->nLocalDof();
1253 const size_type first1_dof_on_proc = _M_X1->dof()->firstDofGlobalCluster( proc_id );
1254 const size_type last1_dof_on_proc = _M_X1->dof()->lastDofGlobalCluster( proc_id );
1255 const size_type first2_dof_on_proc = _M_X2->dof()->firstDofGlobalCluster( proc_id );
1256 const size_type last2_dof_on_proc = _M_X2->dof()->lastDofGlobalCluster( proc_id );
1259 graph_ptrtype sparsity_graph(
new graph_type( _M_X1->dof(),_M_X2->dof() ) );
1262 if ( graph.test( Pattern::ZERO ) ) { sparsity_graph->zero();
return sparsity_graph; }
1266 auto elem_it = _M_X1->mesh()->beginElementWithProcessId( _M_X1->mesh()->worldComm().localRank() );
1267 auto elem_en = _M_X1->mesh()->endElementWithProcessId( _M_X1->mesh()->worldComm().localRank() );
1273 DVLOG(2) <<
"[computeGraph] : graph.test ( Pattern::DEFAULT )=" << graph.test ( Pattern::DEFAULT ) <<
"\n";
1274 DVLOG(2) <<
"[computeGraph] : graph.test ( Pattern::COUPLED )=" << graph.test ( Pattern::COUPLED ) <<
"\n";
1275 DVLOG(2) <<
"[computeGraph] : graph.test ( Pattern::EXTENDED)=" << graph.test ( Pattern::EXTENDED ) <<
"\n";
1276 bool do_less = ( ( graph.test( Pattern::DEFAULT ) &&
1277 ( _M_X1->dof()->nComponents ==
1278 _M_X2->dof()->nComponents ) ) &&
1279 !graph.test( Pattern::COUPLED ) );
1280 std::vector<size_type>
1281 element_dof2( _M_X2->dof()->getIndicesSize() ),
1284 static const uint16_type nDimTest = test_space_type::mesh_type::nDim;
1285 static const uint16_type nDimTrial = trial_space_type::mesh_type::nDim;
1286 static const uint16_type nDimDiffBetweenTestTrial = ( nDimTest > nDimTrial )? nDimTest-nDimTrial : nDimTrial-nDimTest;
1287 for ( ; elem_it != elem_en; ++elem_it )
1289 DVLOG(4) <<
"[Stencil::computePattern] element " << elem_it->id() <<
" on proc " << elem_it->processId() <<
"\n";
1291 const auto & elem = *elem_it;
1293 auto const domains_eid_set = trialElementId( elem.id(), mpl::int_<nDimDiffBetweenTestTrial>() );
1295 auto it_trial=domains_eid_set.begin();
1296 auto const en_trial=domains_eid_set.end();
1297 for ( ; it_trial!=en_trial ; ++it_trial )
1302 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1303 _M_X2->dof()->getIndicesSet( domain_eid, element_dof2 );
1305 _M_X2->dof()->getIndicesSetOnGlobalCluster( domain_eid, element_dof2 );
1310 std::sort( element_dof2.begin(), element_dof2.end() );
1313 const uint16_type n1_dof_on_element = _M_X1->dof()->getIndicesSize();
1314 const uint16_type n2_dof_on_element = element_dof2.size();
1316 for (
size_type i=0; i<n1_dof_on_element; i++ )
1319 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1320 const size_type ig1 = _M_X1->dof()->localToGlobalId( elem.id(), i );
1322 const size_type ig1 = _M_X1->dof()->mapGlobalProcessToGlobalCluster()[_M_X1->dof()->localToGlobalId( elem.id(), i )];
1323 auto theproc = _M_X1->dof()->procOnGlobalCluster( ig1 );
1325 const size_type il1 = _M_X1->dof()->localToGlobalId( elem.id(), i );
1328 const int ndofpercomponent1 = n1_dof_on_element / _M_X1->dof()->nComponents;
1329 const int ncomp1 = i / ndofpercomponent1;
1330 const int ndofpercomponent2 = n2_dof_on_element / _M_X2->dof()->nComponents;
1338 FEELPP_ASSERT ( ig1 >= first1_dof_on_proc )( ig1 )( first1_dof_on_proc ).error (
"invalid dof index" );
1339 FEELPP_ASSERT ( ( ig1 - first1_dof_on_proc ) < sparsity_graph->size() )
1340 ( ig1 )( first1_dof_on_proc )( sparsity_graph->size() ).error(
"invalid dof index" );
1342 graph_type::row_type& row = sparsity_graph->row( ig1 );
1343 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1344 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1348 row.get<0>() = theproc ;
1351 DVLOG(4) <<
"work with row " << ig1 <<
" local index " << ig1 - first1_dof_on_proc <<
"\n";
1355 if ( ncomp1 == ( _M_X2->dof()->nComponents-1 ) )
1356 row.get<2>().insert( element_dof2.begin()+ncomp1*ndofpercomponent2,
1357 element_dof2.end() );
1360 row.get<2>().insert( element_dof2.begin()+ncomp1*ndofpercomponent2,
1361 element_dof2.begin()+( ncomp1+1 )*ndofpercomponent2 );
1366 row.get<2>().insert( element_dof2.begin(), element_dof2.end() );
1370 if ( graph.test( Pattern::EXTENDED ) )
1372 for ( uint16_type ms=0; ms < elem.nNeighbors(); ms++ )
1374 const auto * neighbor = boost::addressof( elem );
1375 size_type neighbor_id = elem.neighbor( ms ).first;
1376 size_type neighbor_process_id = elem.neighbor( ms ).second;
1380 && neighbor_process_id == proc_id )
1383 neighbor = boost::addressof( _M_X1->mesh()->element( neighbor_id,
1384 neighbor_process_id ) );
1386 if ( neighbor_id == neighbor->id() )
1389 neighbor_dof = _M_X2->dof()->getIndicesOnGlobalCluster( neighbor->id() );
1393 if ( ncomp1 == ( _M_X2->dof()->nComponents-1 ) )
1394 row.get<2>().insert( neighbor_dof.begin()+ncomp1*ndofpercomponent2,
1395 neighbor_dof.end() );
1398 row.get<2>().insert( neighbor_dof.begin()+ncomp1*ndofpercomponent2,
1399 neighbor_dof.begin()+( ncomp1+1 )*ndofpercomponent2 );
1405 row.get<2>().insert( neighbor_dof.begin(), neighbor_dof.end() );
1417 DVLOG(2)<<
"[computeGraph<true>] before calling close in " << t.elapsed() <<
"s\n";
1419 DVLOG(2) <<
"[computeGraph<true>] done in " << t.elapsed() <<
"s\n";
1420 DVLOG(2) <<
"[computeGraph<true>] done in " << t.elapsed() <<
"s\n";
1421 return sparsity_graph;
1427 template<
typename MeshType>
1428 struct gmcDefStencil
1430 typedef typename MeshType::element_type::gm_type::precompute_type pc_type;
1431 typedef typename MeshType::element_type::gm_type::precompute_ptrtype pc_ptrtype;
1432 typedef typename MeshType::element_type::gm_type::template Context<vm::POINT, typename MeshType::element_type> gmc_type;
1433 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1436 template<
typename EltType>
1437 struct gmcDefStencil
1439 typedef typename EltType::gm_type::precompute_type pc_type;
1440 typedef typename EltType::gm_type::precompute_ptrtype pc_ptrtype;
1441 typedef typename EltType::gm_type::template Context<vm::POINT, EltType> gmc_type;
1442 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
1445 template<
typename FaceType>
1446 struct gmcDefFaceStencil
1448 typedef typename boost::remove_reference<FaceType>::type FaceType1;
1449 typedef typename boost::remove_const<FaceType1>::type FaceType2;
1451 typedef typename FaceType2::super2::template Element<FaceType2>::type EltType;
1452 typedef typename gmcDefStencil<EltType>::pc_type pc_type;
1453 typedef typename gmcDefStencil<EltType>::pc_ptrtype pc_ptrtype;
1454 typedef typename gmcDefStencil<EltType>::gmc_type gmc_type;
1455 typedef typename gmcDefStencil<EltType>::gmc_ptrtype gmc_ptrtype;
1458 template<
typename ImType,
typename EltType>
1459 typename gmcDefStencil<EltType>::gmc_ptrtype
1460 gmcStencil( mpl::size_t<MESH_ELEMENTS> , EltType
const& elem )
1462 typedef typename gmcDefStencil<EltType>::pc_type pc_type;
1463 typedef typename gmcDefStencil<EltType>::pc_ptrtype pc_ptrtype;
1464 typedef typename gmcDefStencil<EltType>::gmc_type gmc_type;
1465 typedef typename gmcDefStencil<EltType>::gmc_ptrtype gmc_ptrtype;
1468 pc_ptrtype geopc(
new pc_type( elem.gm(), theim.points() ) );
1469 gmc_ptrtype gmc(
new gmc_type( elem.gm(), elem, geopc ) );
1473 template<
typename ImType,
typename FaceType>
1474 typename gmcDefFaceStencil<FaceType>::gmc_ptrtype
1475 gmcStencil( mpl::size_t<MESH_FACES> , FaceType
const& theface )
1477 typedef typename gmcDefFaceStencil<FaceType>::pc_type pc_type;
1478 typedef typename gmcDefFaceStencil<FaceType>::pc_ptrtype pc_ptrtype;
1479 typedef typename gmcDefFaceStencil<FaceType>::gmc_type gmc_type;
1480 typedef typename gmcDefFaceStencil<FaceType>::gmc_ptrtype gmc_ptrtype;
1482 typedef typename QuadMapped<ImType>::permutation_type permutation_type;
1483 typedef typename QuadMapped<ImType>::permutation_points_type permutation_points_type;
1486 QuadMapped<ImType> qm;
1487 permutation_points_type ppts( qm(theim) );
1489 std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( theim.nFaces() );
1490 for ( uint16_type __f = 0; __f < theim.nFaces(); ++__f )
1492 for ( permutation_type __p( permutation_type::IDENTITY );
1493 __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1495 __geopc[__f][__p] = pc_ptrtype(
new pc_type( theface.element( 0 ).gm(), ppts[__f].find( __p )->second ) );
1499 uint16_type __face_id_in_elt_0 = theface.pos_first();
1500 gmc_ptrtype gmc(
new gmc_type( theface.element( 0 ).gm(),
1501 theface.element( 0 ),
1503 __face_id_in_elt_0 ) );
1506 template<
typename EltType>
1508 gmcUpdateStencil( mpl::size_t<MESH_ELEMENTS> , EltType
const& elem,
typename gmcDefStencil<EltType>::gmc_ptrtype &gmc )
1510 gmc->update( elem );
1512 template<
typename FaceType>
1514 gmcUpdateStencil( mpl::size_t<MESH_FACES> , FaceType
const& theface,
typename gmcDefFaceStencil<FaceType>::gmc_ptrtype &gmc )
1516 gmc->update( theface.element( 0 ), theface.pos_first() );
1518 template<
typename EltType>
1520 idEltStencil( mpl::size_t<MESH_ELEMENTS> , EltType
const& elem )
1524 template<
typename FaceType>
1526 idEltStencil( mpl::size_t<MESH_FACES> , FaceType
const& theface )
1528 return theface.element( 0 ).id();
1533 template<
typename X1,
typename X2,
typename RangeItTestType>
1534 typename Stencil<X1,X2,RangeItTestType>::graph_ptrtype
1535 Stencil<X1,X2,RangeItTestType>::computeGraphInCaseOfInterpolate(
size_type hints, mpl::bool_<true> )
1539 typedef mpl::int_<20> order_1d_type;
1540 typedef mpl::int_<12> order_2d_type;
1541 typedef mpl::int_<8> order_3d_type;
1543 typedef typename test_space_type::mesh_type test_mesh_type;
1544 typedef typename trial_space_type::mesh_type trial_mesh_type;
1545 typedef typename mpl::if_< mpl::equal_to<mpl::int_<test_mesh_type::nDim>,mpl::int_<1> >,
1547 typename mpl::if_<mpl::equal_to<mpl::int_<test_mesh_type::nDim>,mpl::int_<2> >,
1549 order_3d_type >::type>::type order_used_type;
1550 typedef typename mpl::if_<mpl::bool_<test_mesh_type::element_type::is_simplex>,
1551 mpl::identity<typename _Q<order_used_type::value>::template apply<test_mesh_type::element_type::nDim, typename test_mesh_type::value_type, Simplex>::type >,
1552 mpl::identity<typename _Q<order_used_type::value>::template apply<test_mesh_type::element_type::nDim, typename test_mesh_type::value_type, Hypercube>::type >
1553 >::type::type theim_type;
1555 typedef typename test_mesh_type::Localization::matrix_node_type matrix_node_type;
1559 const size_type proc_id = _M_X1->mesh()->comm().rank();
1560 const size_type n1_dof_on_proc = _M_X1->nLocalDof();
1562 const size_type first1_dof_on_proc = _M_X1->dof()->firstDof( proc_id );
1563 const size_type last1_dof_on_proc = _M_X1->dof()->lastDof( proc_id );
1564 const size_type first2_dof_on_proc = _M_X2->dof()->firstDof( proc_id );
1565 const size_type last2_dof_on_proc = _M_X2->dof()->lastDof( proc_id );
1567 graph_ptrtype sparsity_graph(
new graph_type( _M_X1->dof(),_M_X2->dof() ) );
1570 typedef typename test_mesh_type::element_const_iterator mesh_element_const_iterator;
1571 mesh_element_const_iterator elem_it = _M_X1->mesh()->beginElementWithProcessId( proc_id );
1572 const mesh_element_const_iterator elem_en = _M_X1->mesh()->endElementWithProcessId( proc_id );
1573 auto iDimRange = mpl::size_t<MESH_ELEMENTS>();
1576 auto rangeTest = this->rangeiterator<0,0>(mpl::bool_<rangeiteratorType<0,0>::hasnotfindrange_type::value>());
1577 auto iDimRange = rangeTest.template get<0>();
1578 auto elem_it = rangeTest.template get<1>();
1579 auto elem_en = rangeTest.template get<2>();
1581 if ( elem_it==elem_en )
return sparsity_graph;
1585 auto locToolForXh2 = _M_X2->mesh()->tool_localization();
1586 locToolForXh2->updateForUse();
1587 bool doExtrapolationAtStartXh2 = locToolForXh2->doExtrapolation();
1588 if (doExtrapolationAtStartXh2) locToolForXh2->setExtrapolation(
false );
1589 const auto nbNearNeighborAtStartTrial = locToolForXh2->kdtree()->nPtMaxNearNeighbor();
1590 bool notUseOptLocTrial = trial_mesh_type::nDim!=trial_mesh_type::nRealDim;
1591 if (notUseOptLocTrial) locToolForXh2->kdtree()->nbNearNeighbor(trial_mesh_type::element_type::numPoints);
1594 auto locToolForXh1 = _M_X1->mesh()->tool_localization();
1595 locToolForXh1->updateForUse();
1596 bool doExtrapolationAtStartXh1 = locToolForXh1->doExtrapolation();
1597 if (doExtrapolationAtStartXh1) locToolForXh1->setExtrapolation(
false );
1600 matrix_node_type ptsReal( elem_it->vertices().size1(), 1 );
1604 #if FEELPP_EXPORT_GRAPH
1605 std::map<size_type,std::list<size_type> > mapBetweenMeshes;
1608 std::vector<size_type> element_dof1_range, element_dof1, element_dof2;
1609 std::set<size_type> neighLocalizedInXh1;
1611 std::set<size_type > listTup;
1612 auto gmc = Feel::detail::gmcStencil<theim_type>( iDimRange,*elem_it );
1617 if ( _M_X1->nDof()>1 )
1619 for ( ; elem_it != elem_en; ++elem_it )
1621 auto const& elem = *elem_it;
1624 element_dof1_range = _M_X1->dof()->getIndices( elem.id(), iDimRange );
1625 element_dof1 = _M_X1->dof()->getIndices( Feel::detail::idEltStencil( iDimRange,elem ) );
1626 const uint16_type n1_dof_on_element_range = element_dof1_range.size();
1627 const uint16_type n1_dof_on_element = element_dof1.size();
1629 std::vector<boost::tuple<bool,size_type> > hasFinds( n1_dof_on_element_range,boost::make_tuple(
false,
invalid_size_type_value ) );
1631 for (
size_type i=0; i<n1_dof_on_element_range; i++ )
1633 const size_type ig1 = element_dof1_range[i];
1634 auto const ptRealDof = boost::get<0>( _M_X1->dof()->dofPoint( ig1 ) );
1636 ublas::column(ptsReal,0 ) = ptRealDof;
1638 auto resLocalisationInXh2 = locToolForXh2->run_analysis(ptsReal,IdEltInXh2,elem_it->vertices(),mpl::int_<0>());
1639 IdEltInXh2 = resLocalisationInXh2.template get<1>();
1640 bool hasFind = resLocalisationInXh2.template get<0>()[0];
1646 listTup.insert( IdEltInXh2 );
1647 hasFinds[i] = boost::make_tuple(
true,IdEltInXh2 );
1649 auto const& geoelt2 = _M_X2->mesh()->element( IdEltInXh2 );
1651 std::vector<size_type> neighbor_ids;
1652 for ( uint16_type ms=0; ms < geoelt2.nNeighbors(); ms++ )
1654 size_type neighbor_id = geoelt2.neighbor( ms ).first;
1657 neighbor_ids.push_back( neighbor_id );
1660 auto resIsIn = locToolForXh2->isIn( neighbor_ids,ptRealDof );
1663 for (
auto it=resIsIn.template get<1>().begin(),en=resIsIn.template get<1>().end(); it<en; ++it,++cpt )
1667 listTup.insert( neighbor_ids[cpt] );
1671 auto res_it = listTup.begin();
1672 auto res_en = listTup.end();
1673 for ( ; res_it != res_en ; ++res_it )
1675 #if FEELPP_EXPORT_GRAPH
1676 mapBetweenMeshes[*res_it].push_back( elem.id() );
1680 for (
size_type ii=0; ii<n1_dof_on_element; ii++ )
1682 const size_type ig1ongraph = element_dof1[ii];
1684 element_dof2 = _M_X2->dof()->getIndices( *res_it );
1686 graph_type::row_type& row = sparsity_graph->row( ig1ongraph );
1687 bool is_on_proc = ( ig1ongraph >= first1_dof_on_proc ) && ( ig1ongraph <= last1_dof_on_proc );
1691 row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1701 graph_type::row_type& row = sparsity_graph->row( ig1 );
1702 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1705 row.template get<2>().clear();
1710 uint16_type nbFind = hasFinds.size();
1714 for ( uint16_type nF = 0 ; nF<nbFind && !doQ ; ++nF )
1715 if ( hasFinds[nF].
template get<0>() )
1717 id1=hasFinds[nF].template get<1>();
1719 for ( uint16_type nF2 = 0 ; nF2<nbFind && !doQ ; ++nF2 )
1721 if ( hasFinds[nF2].
template get<0>() )
1723 id2=hasFinds[nF2].template get<1>();
1733 Feel::detail::gmcUpdateStencil( iDimRange,elem,gmc );
1735 for (
int q = 0; q < gmc->nPoints(); ++ q )
1737 ublas::column(ptsReal,0 ) = gmc->xReal( q );
1739 auto const resQuad = locToolForXh2->searchElement( gmc->xReal( q ) );
1741 if ( resQuad.template get<0>() )
1743 IdEltInXh2 = resQuad.template get<1>();
1744 #if FEELPP_EXPORT_GRAPH
1745 mapBetweenMeshes[IdEltInXh2].push_back( elem.id() );
1747 element_dof2 = _M_X2->dof()->getIndices( IdEltInXh2 );
1749 for (
size_type i=0; i<n1_dof_on_element; i++ )
1752 graph_type::row_type& row = sparsity_graph->row( ig1 );
1753 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1757 row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1767 for ( ; elem_it != elem_en; ++elem_it )
1769 auto const& elem = *elem_it;
1772 element_dof1 = _M_X1->dof()->getIndices( Feel::detail::idEltStencil( iDimRange,elem ) );
1774 const uint16_type nPtGeo = elem.G().size2();
1775 std::vector<boost::tuple<bool,size_type> > hasFinds( nPtGeo,boost::make_tuple(
false,
invalid_size_type_value ) );
1784 typename matrix_node<typename test_mesh_type::value_type>::type ptReal( test_mesh_type::nRealDim , 1 );
1785 ublas::column( ptReal ,0 ) = ublas::column( elem.G(),i );
1790 auto resTemp = locToolForXh2->searchElement( ublas::column( elem.G(),i ) );
1791 bool hasFind = resTemp.template get<0>();
1792 std::set<size_type > listTup;
1796 listTup.insert( resTemp.template get<1>() );
1797 hasFinds[i] = boost::make_tuple(
true,resTemp.template get<1>() );
1800 size_type idElt2 = resTemp.template get<1>();
1801 auto const& geoelt2 = _M_X2->mesh()->element( idElt2 );
1802 std::vector<size_type> neighbor_ids( geoelt2.nNeighbors() );
1804 for ( uint16_type ms=0; ms < geoelt2.nNeighbors(); ms++ )
1806 size_type neighbor_id = geoelt2.neighbor( ms ).first;
1813 auto resIsIn = locToolForXh2->isIn( neighbor_ids,ublas::column( elem.G(),i ) );
1816 for (
auto it=resIsIn.template get<1>().begin(),en=resIsIn.template get<1>().end(); it<en; ++it,++cpt )
1820 listTup.insert( neighbor_ids[cpt] );
1823 else if ( test_mesh_type::nDim==trial_mesh_type::nDim )
1825 for (
auto it_neigh = neighbor_ids.begin(), en_neigh = neighbor_ids.end() ; it_neigh != en_neigh ; ++it_neigh )
1827 if ( ( sparsity_graph->row( ig1 ) ).
template get<2>().find( *it_neigh )==( sparsity_graph->row( ig1 ) ).
template get<2>().end() )
1829 if ( neighLocalizedInXh1.find( *it_neigh )==neighLocalizedInXh1.end() )
1831 neighLocalizedInXh1.insert( *it_neigh );
1832 bool findNeih=
false;
1833 auto const& geoeltNEW = _M_X2->mesh()->element( *it_neigh );
1834 const uint16_type nPtGeoBis = _M_X2->mesh()->element( *it_neigh ).G().size2();
1836 for (
size_type iii=0; iii<nPtGeoBis && !findNeih ; iii++ )
1838 auto resBis = locToolForXh1->searchElement( ublas::column( geoeltNEW.G(),iii ) );
1840 if ( resBis.template get<0>() )
1842 listTup.insert( *it_neigh );
1852 auto res_it = listTup.begin();
1853 auto res_en = listTup.end();
1854 for ( ; res_it != res_en ; ++res_it )
1856 #if FEELPP_EXPORT_GRAPH
1858 mapBetweenMeshes[*res_it].push_back( elem.id() );
1862 element_dof2 = _M_X2->dof()->getIndices( *res_it );
1864 graph_type::row_type& row = sparsity_graph->row( ig1 );
1865 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1869 row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1874 uint16_type nbFind = hasFinds.size();
1878 for ( uint16_type nF = 0 ; nF<nbFind && !doQ ; ++nF )
1879 if ( hasFinds[nF].
template get<0>() )
1881 id1=hasFinds[nF].template get<1>();
1883 for ( uint16_type nF2 = 0 ; nF2<nbFind && !doQ ; ++nF2 )
1885 if ( hasFinds[nF2].
template get<0>() )
1887 id2=hasFinds[nF2].template get<1>();
1897 Feel::detail::gmcUpdateStencil(iDimRange,elem,gmc);
1899 for (
int q = 0; q < gmc->nPoints(); ++ q )
1901 auto resQuad = locToolForXh2->searchElement( gmc->xReal( q ) );
1903 if ( resQuad.template get<0>() )
1905 #if FEELPP_EXPORT_GRAPH
1906 mapBetweenMeshes[resQuad.template get<1>()].push_back( elem.id() );
1908 element_dof2 = _M_X2->dof()->getIndices( resQuad.template get<1>() );
1912 graph_type::row_type& row = sparsity_graph->row( ig1 );
1913 bool is_on_proc = ( ig1 >= first1_dof_on_proc ) && ( ig1 <= last1_dof_on_proc );
1917 row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1927 if (doExtrapolationAtStartXh1) locToolForXh1->setExtrapolation(
true );
1928 if (doExtrapolationAtStartXh2) locToolForXh2->setExtrapolation(
true );
1929 if (notUseOptLocTrial) locToolForXh2->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTrial);
1935 #if FEELPP_EXPORT_GRAPH
1937 typedef mesh_1_type mesh_export_type;
1938 typedef FunctionSpace<mesh_1_type, bases<Lagrange<0, Scalar,Discontinuous> > > space_disc_type;
1939 auto spaceGraphProj = space_disc_type::New( _M_X1->mesh() );
1940 auto elem_itt = _M_X1->mesh()->beginElementWithProcessId( proc_id );
1941 auto elem_ent = _M_X1->mesh()->endElementWithProcessId( proc_id );
1943 typedef trial_mesh_type mesh_export_type;
1944 typedef FunctionSpace<mesh_export_type, bases<Lagrange<0, Scalar,Discontinuous> > > space_disc_type;
1945 auto spaceGraphProj = space_disc_type::New( _M_X2->mesh() );
1946 auto elem_itt = _M_X2->mesh()->beginElementWithProcessId( proc_id );
1947 auto elem_ent = _M_X2->mesh()->endElementWithProcessId( proc_id );
1949 auto graphProj = spaceGraphProj->element();
1952 for ( ; elem_itt != elem_ent; ++elem_itt )
1954 if ( mapBetweenMeshes.find( elem_itt->id() ) != mapBetweenMeshes.end() )
1956 if ( mapBetweenMeshes.find( elem_itt->id() )->second.size()>0 )
1958 element_dof1 = spaceGraphProj->dof()->getIndices( elem_itt->id() );
1959 const uint16_type n1_dof_on_element = element_dof1.size();
1961 for ( uint16_type i=0; i<n1_dof_on_element; i++ )
1965 graphProj( ig1 ) = 1;
1972 exporter->step( 0 )->setMesh( graphProj.mesh() );
1973 exporter->step( 0 )->add(
"graphProj", graphProj );
1977 return sparsity_graph;