Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
stencil.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2011-12-23
7 
8  Copyright (C) 2011 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
24 #ifndef __galerkingraph_H
25 #define __galerkingraph_H 1
26 
27 #define FEELPP_EXPORT_GRAPH 0
28 #if FEELPP_EXPORT_GRAPH
29 #include <feel/feelfilters/exporter.hpp>
30 #endif
31 
32 #include <feel/feelvf/pattern.hpp>
33 #include <feel/feelvf/block.hpp>
35 #include <feel/feeldiscr/functionspace.hpp>
36 
37 
38 #if 1
39 namespace Feel
40 {
41 namespace detail
42 {
43 template<typename BFType, typename Space1Type>
44 struct compute_graph3
45 {
46  compute_graph3( BFType* bf, boost::shared_ptr<Space1Type> const& space1, size_type trial_index, size_type hints )
47  :
48  M_stencil( bf ),
49  M_space1( space1 ),
50  M_test_index( 0 ),
51  M_trial_index( trial_index ),
52  M_hints( hints )
53  {}
54 
55  template <typename Space2>
56  void operator()( boost::shared_ptr<Space2> const& space2 ) const
57  {
58  if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
59  {
60  if ( M_stencil->isBlockPatternZero( M_test_index,M_trial_index ) )
61  {
62 #if 0
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 );
70 #else
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 );
77 #endif
78 #endif
79  typename BFType::graph_ptrtype zerograph( new typename BFType::graph_type( space2->dof(), M_space1->dof() ) );
80  zerograph->zero();
81  M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ), M_stencil->trialSpace()->nDofStart( M_trial_index ) , zerograph );
82  }
83 
84  else
85  {
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,
91  _close=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>() )
94  );
95 
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() );
100  else
101  M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
102  M_stencil->trialSpace()->nDofStart( M_trial_index ),
103  thestencil->graph() );
104  }
105  } // if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
106 
107  ++M_test_index;
108  }
109 
110 
111  mutable BFType* M_stencil;
112  boost::shared_ptr<Space1Type> const& M_space1;
113  mutable size_type M_test_index;
114  size_type M_trial_index;
115  size_type M_hints;
116 };
117 
118 
119 template<typename BFType, typename Space1Type>
120 struct compute_graph2
121 {
122  compute_graph2( BFType* bf, boost::shared_ptr<Space1Type> const& space1, size_type test_index, size_type hints )
123  :
124  M_stencil( bf ),
125  M_space1( space1 ),
126  M_test_index( test_index ),
127  M_trial_index( 0 ),
128  M_hints( hints )
129  {}
130 
131  template <typename Space2>
132  void operator()( boost::shared_ptr<Space2> const& space2 ) const
133  {
134  if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
135  {
136  if ( M_stencil->isBlockPatternZero( M_test_index,M_trial_index ) )
137  {
138 #if 0
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 );
146 #else
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 );
153 #endif
154 #endif
155  typename BFType::graph_ptrtype zerograph( new typename BFType::graph_type( M_space1->dof(), space2->dof() ) );
156  zerograph->zero();
157  M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
158  M_stencil->trialSpace()->nDofStart( M_trial_index ),
159  zerograph );
160  }
161 
162  else
163  {
164 
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,
170  _close=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>() )
173  );
174 
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() );
179  else
180  M_stencil->mergeGraph( M_stencil->testSpace()->nDofStart( M_test_index ),
181  M_stencil->trialSpace()->nDofStart( M_trial_index ),
182  thestencil->graph() );
183  }
184  } // if ( M_stencil->testSpace()->worldsComm()[M_test_index].isActive() )
185 
186  ++M_trial_index;
187  }
188 
189 
190  mutable BFType* M_stencil;
191  boost::shared_ptr<Space1Type> const& M_space1;
192  size_type M_test_index;
193  mutable size_type M_trial_index;
194  size_type M_hints;
195 };
196 
197 
198 template<typename BFType>
199 struct compute_graph1
200 {
201  compute_graph1( BFType* bf, size_type hints )
202  :
203  M_stencil( bf ),
204  M_test_index( 0 ),
205  M_hints( hints )
206  {}
207 
208  template <typename Space1>
209  void operator()( boost::shared_ptr<Space1> const& space1 ) const
210  {
211  fusion::for_each( M_stencil->trialSpace()->functionSpaces(),
212  compute_graph2<BFType,Space1>( M_stencil, space1, M_test_index, M_hints ) );
213  ++M_test_index;
214  }
215  mutable BFType* M_stencil;
216  mutable size_type M_test_index;
217  size_type M_hints;
218 };
219 
220 }
221 
222 
223 
224 
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)
228 {
229  return fusion::make_pair< fusion::pair<mpl::int_<I>,mpl::int_<J> > >( r);
230 }
231 
232 
233 struct stencilRangeMapTypeBase {};
234 
235 struct stencilRangeMap0Type
236  :
237  public stencilRangeMapTypeBase,
238  public fusion::map<>
239 {
240  typedef fusion::map<> super_type;
241 
242  stencilRangeMap0Type()
243  :
244  super_type()
245  {}
246 
247  static bool isNullRange() { return true; }
248 };
249 
250 template <typename ThePair1Type>
251 struct stencilRangeMap1Type
252  :
253  public stencilRangeMapTypeBase,
254  public fusion::map< fusion::pair< typename ThePair1Type::first_type, typename ThePair1Type::second_type > >
255 {
256  typedef typename ThePair1Type::first_type key1_type;
257  typedef fusion::map< fusion::pair< key1_type, typename ThePair1Type::second_type > > super_type;
258 
259  stencilRangeMap1Type( ThePair1Type const& p )
260  :
261  super_type( p )
262  {}
263 
264  static bool isNullRange() { return false; }
265 };
266 
267 template <typename ThePair1Type,typename ThePair2Type>
268 struct stencilRangeMap2Type
269  :
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 > >
273 {
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;
278 
279  stencilRangeMap2Type( ThePair1Type const& p1, ThePair2Type const& p2 )
280  :
281  super_type( p1,p2 )
282  {}
283 
284  static bool isNullRange() { return false; }
285 };
286 
287 stencilRangeMap0Type stencilRangeMap();
288 
289 template <typename ThePair1Type>
290 stencilRangeMap1Type<ThePair1Type>
291 stencilRangeMap( ThePair1Type const& p)
292 {
293  return stencilRangeMap1Type<ThePair1Type>( p.second );
294 }
295 
296 template <typename ThePair1Type,typename ThePair2Type>
297 stencilRangeMap2Type<ThePair1Type,ThePair2Type>
298 stencilRangeMap( ThePair1Type const& p1, ThePair2Type const& p2)
299 {
300  return stencilRangeMap2Type<ThePair1Type,ThePair2Type>( p1.second, p2.second );
301 }
302 
303 
304 
305 template<typename X1, typename X2,typename RangeIteratorTestType = stencilRangeMap0Type >
306 class Stencil
307 {
308 public:
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;
316 
317  typedef RangeIteratorTestType rangeiterator_test_type;
318 
319  Stencil( test_space_ptrtype Xh, trial_space_ptrtype Yh,
320  size_type graph_hints,
321  BlocksStencilPattern block_pattern=BlocksStencilPattern(1,1,Pattern::HAS_NO_BLOCK_PATTERN),
322  bool diag_is_nonzero=false,
323  bool close=true,
324  rangeiterator_test_type r=rangeiterator_test_type())
325  :
326  _M_X1( Xh ),
327  _M_X2( Yh ),
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 ) ),
332 #else
333  M_graph( new graph_type( Xh->dof(),Yh->dof() ) ),
334 #endif
335  M_block_pattern( block_pattern ),
336  M_rangeIteratorTest( r )
337  {
338  // init block_pattern if empty
339  uint16_type nbSubSpace1 = _M_X1->nSubFunctionSpace();
340  uint16_type nbSubSpace2 = _M_X2->nSubFunctionSpace();
341 
342  if ( this->isBlockPatternNoPattern( 0,0 ) )
343  {
344  M_block_pattern = BlocksStencilPattern(nbSubSpace1,nbSubSpace2,graph_hints);
345  }
346 
347  M_graph = this->computeGraph( graph_hints );
348 
349  if ( diag_is_nonzero && _M_X1->nLocalDofWithoutGhost()>0 && _M_X2->nLocalDofWithoutGhost()>0 ) M_graph->addMissingZeroEntriesDiagonal();
350 
351  if ( close ) M_graph->close();
352  }
353 
354  Stencil( test_space_ptrtype Xh, trial_space_ptrtype Yh, size_type graph_hints, graph_ptrtype g, rangeiterator_test_type r=rangeiterator_test_type() )
355  :
356  _M_X1( Xh ),
357  _M_X2( Yh ),
358  M_graph( g ),
359  M_block_pattern(Xh->nSubFunctionSpace(),Yh->nSubFunctionSpace(),size_type( graph_hints/*Pattern::HAS_NO_BLOCK_PATTERN*/ )),
360  M_rangeIteratorTest( r )
361  {}
362 
363 
364  BlocksStencilPattern blockPattern() const
365  {
366  return M_block_pattern;
367  }
368 
369  size_type blockPattern( uint16_type i,uint16_type j ) const
370  {
371  return M_block_pattern(i,j);
372  }
373 
374  bool isBlockPatternNoPattern( uint16_type i,uint16_type j ) const
375  {
376  Feel::Context ctx( M_block_pattern(i,j) );
377  return ctx.test( HAS_NO_BLOCK_PATTERN );
378 
379  }
380  bool isBlockPatternZero( uint16_type i,uint16_type j ) const
381  {
382  Feel::Context ctx( M_block_pattern(i,j) );
383  return ctx.test( ZERO );
384  }
385 
386  graph_ptrtype computeGraph( size_type hints );
387 
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> );
393 
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> );
399 
400 
401  void mergeGraph( int row, int col, graph_ptrtype g );
402  void mergeGraphMPI( size_type test_index, size_type trial_index,
403  DataMap const& mapOnTest, DataMap const& mapOnTrial,
404  graph_ptrtype g);
405 
406  test_space_ptrtype testSpace() const
407  {
408  return _M_X1;
409  }
410  trial_space_ptrtype trialSpace() const
411  {
412  return _M_X2;
413  }
414  graph_ptrtype graph() const
415  {
416  return M_graph;
417  }
418  graph_ptrtype graph()
419  {
420  return M_graph;
421  }
422 private:
423  std::set<size_type> trialElementId( size_type test_eid, mpl::int_<0> )
424  {
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 )
429  {
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";
432  if ( domain_eid != invalid_size_type_value ) idsFind.insert( domain_eid );
433  }
434  else if( trial_related_to_test )
435  {
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";
438  if ( domain_eid != invalid_size_type_value ) idsFind.insert( domain_eid );
439  }
440  else // same mesh
441  {
442  idsFind.insert(test_eid);
443  }
444 
445  return idsFind;
446  }
447  std::set<size_type> trialElementId( size_type test_eid, mpl::int_<1> )
448  {
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 )
453  {
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";
456  if ( domain_eid != invalid_size_type_value ) idsFind.insert( domain_eid );
457  }
458  else if( trial_related_to_test )
459  {
460  auto const& eltTest = _M_X1->mesh()->element(test_eid);
461  for (uint16_type f=0;f< _M_X1->mesh()->numLocalFaces();++f)
462  {
463  const size_type idFind = _M_X2->mesh()->meshToSubMesh( eltTest.face(f).id() );
464  if ( idFind != invalid_size_type_value ) idsFind.insert( idFind );
465  }
466  DVLOG(2) << "[trial_related_to_test<1>] test element id: " << test_eid << " idsFind.size() "<< idsFind.size() << "\n";
467  }
468  else
469  {
470  CHECK ( false ) << "[trial_related_to_test<1>] : test and trial mesh can not be the same here\n";
471  }
472  return idsFind;
473  }
474 
475 
476 
477 public :
478  template <int I,int J>
479  struct rangeiteratorType
480  {
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;
483 
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;
488 
489  typedef typename mpl::if_< hasnotfindrange_type,
490  mpl::identity< defaultrange_type >,
491  mpl::identity< resultfindrange_type >
492  >::type::type type;
493  };
494 #if 0
495  template <int I,int J>
496  typename rangeiteratorType<I,J>::type
497  rangeiterator() const
498  {
499  return rangeiterator<I,J>( mpl::bool_<rangeiteratorType<I,J>::hasnotfindrange_type::value>() );
500  }
501 #endif
502  template <int I,int J>
503  typename rangeiteratorType<I,J>::defaultrange_type
504  rangeiterator(mpl::bool_<true> ) const
505  {
506  return elements(_M_X1->mesh());
507  }
508  template <int I,int J>
509  typename rangeiteratorType<I,J>::resultfindrange_type::second_type
510  rangeiterator(mpl::bool_<false> ) const
511  {
512  typedef fusion::pair<mpl::int_<I>,mpl::int_<J> > key_type;
513  return fusion::at_key< key_type >(M_rangeIteratorTest);//.second;
514  //return *fusion::find< key_type >(M_rangeIteratorTest);//->second;
515  }
516  template <int I,int J>
517  stencilRangeMap0Type
518  subRangeIterator( mpl::bool_<true> )
519  {
520  return stencilRangeMap0Type();
521  }
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> )
525  {
526  return stencilRangeMap( stencilRange<0,0>( rangeiterator<I,J>(mpl::bool_<false>()) ) );
527  }
528 
529 private:
530 
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;
536 };
537 namespace detail
538 {
539 template<typename Args>
540 struct compute_stencil_type
541 {
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;
547 
548 };
549 
550 }
551 
552 class StencilManagerImpl:
553  public std::map<boost::tuple<boost::shared_ptr<FunctionSpaceBase>,
554  boost::shared_ptr<FunctionSpaceBase>,
555  size_type,
556  std::vector<size_type>,
557  bool >, boost::shared_ptr<GraphCSR> >,
558 public boost::noncopyable
559 {
560 public:
561  typedef boost::shared_ptr<GraphCSR> graph_ptrtype;
562  typedef boost::tuple<boost::shared_ptr<FunctionSpaceBase>,
563  boost::shared_ptr<FunctionSpaceBase>,
564  size_type,
565  std::vector<size_type>,
566  bool > key_type;
567  typedef std::map<key_type, graph_ptrtype> graph_manager_type;
568 
569 };
570 
571 typedef Feel::Singleton<StencilManagerImpl> StencilManager;
572 
576 void stencilManagerGarbage(StencilManagerImpl::key_type const& key);
578 void stencilManagerAdd(StencilManagerImpl::key_type const& key,StencilManagerImpl::graph_ptrtype graph);
580 void stencilManagerPrint();
581 
582 extern BlocksStencilPattern default_block_pattern;
583 
584 BOOST_PARAMETER_FUNCTION(
585  ( typename detail::compute_stencil_type<Args>::ptrtype ), // 1. return type
586  stencil, // 2. name of the function template
587  tag, // 3. namespace of tag types
588  ( required // 4. one required parameter, and
589  ( test, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
590  ( trial, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
591  )
592  ( optional // four optional parameters, with defaults
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() )
599  )
600 )
601 {
602  if ( collect_garbage )
603  {
604  // cleanup memory before doing anything
606  }
607 
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;
611 
612  // we look into the spaces dictionary for existing graph
613  auto git = StencilManager::instance().find( boost::make_tuple( test, trial, pattern, pattern_block.getSetOfBlocks(), diag_is_nonzero ) );
614 
615  if ( git != StencilManager::instance().end() && range.isNullRange() )
616  {
617  //std::cout << "Found a stencil in manager (" << test.get() << "," << trial.get() << "," << pattern << ")\n";
618  auto s = stencil_ptrtype( new stencil_type( test, trial, pattern, git->second, range ) );
619  return s;
620  }
621 
622  else
623  {
624  // look for transposed stencil if it exist and transpose it to get the stencil
625  auto git_trans = StencilManager::instance().find( boost::make_tuple( trial, test, pattern, pattern_block.transpose().getSetOfBlocks(), diag_is_nonzero ) );
626 
627  if ( git_trans != StencilManager::instance().end() && range.isNullRange() )
628  {
629  auto g = git_trans->second->transpose(close);
630  //auto g = git_trans->second->transpose();
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 ) );
633  //std::cout << "Found a transposed stencil in manager (" << test.get() << "," << trial.get() << "," << pattern << ")\n";
634  return s;
635  }
636 
637  else
638  {
639  //std::cout << "Creating a new stencil in manager (" << test.get() << "," << trial.get() << "," << pattern << ")\n";
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();
642  return s;
643  }
644  }
645 }
646 
647 
648 namespace detail
649 {
650 template<typename BidirectionalIterator>
651 inline
652 void
653 sortSparsityRow ( const BidirectionalIterator begin,
654  BidirectionalIterator middle,
655  const BidirectionalIterator end )
656 {
657  if ( ( begin == middle ) || ( middle == end ) ) return;
658 
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)" );
665 
666  while ( middle != end )
667  {
668  BidirectionalIterator
669  b = middle,
670  a = b-1;
671 
672  // Bubble-sort the middle value downward
673  while ( !( *a < *b ) ) // *a & *b are less-than comparable, so use <
674  {
675  std::swap ( *a, *b );
676 
677  if ( a == begin ) break;
678 
679  b=a;
680  --a;
681  }
682 
683  ++middle;
684  }
685 
686  // Assure the algorithm worked if we are in DEBUG mode
687 #ifdef DEBUG
688  {
689  // SGI STL extension!
690  // assert (std::is_sorted(begin,end));
691 
692  BidirectionalIterator
693  prev = begin,
694  first = begin;
695 
696  for ( ++first; first != end; prev=first, ++first )
697  if ( *first < *prev )
698  assert( false );
699  }
700 #endif
701 
702  // Make sure the two ranges did not contain any common elements
703  assert ( std::unique ( begin, end ) == end );
704 } //
705 
706 }
707 template<typename X1, typename X2, typename RangeItTestType>
708 void
709 Stencil<X1,X2,RangeItTestType>::mergeGraph( int row, int col, graph_ptrtype g )
710 {
711  boost::timer tim;
712  DVLOG(2) << "[merge graph] for composite bilinear form\n";
713  DVLOG(2) << "[mergeGraph] row = " << row << "\n";
714  DVLOG(2) << "[mergeGraph] col = " << col << "\n";
715 
716  // nothing yet in store
717  //if ( !M_graph || M_graph->empty() )
718  if ( 0 )
719  {
720  DVLOG(2) << "[merge graph] nothing yet in store, copy graph\n";
721  M_graph = g;
722 
723 #if 0
724  typename graph_type::const_iterator it = g->begin();
725  typename graph_type::const_iterator en = g->end();
726 
727  for ( ; it != en; ++it )
728  {
729  std::vector<size_type> const& ivec = boost::get<2>( it->second );
730 
731  for ( int i = 0; i < ivec.size(); ++i )
732  {
733  DVLOG(2) << "[mergeGraph] ivec[" << i << "] = " << ivec[i] << "\n";
734  }
735  }
736 
737 #endif // 9
738  }
739 
740  else
741  {
742  //std::cout << "\n row " << row << " col " << col << " with god rank" << this->testSpace()->worldComm().godRank() << std::endl;
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();
746 
747  for ( ; it != en; ++it )
748  {
749  int theglobalrow = row+it->first;
750  int thelocalrow;// warning : not the same in parallel
751 
752  if ( this->testSpace()->worldComm().globalSize()>1 )
753  thelocalrow = /*row +*/ boost::get<1>( it->second );
754 
755  else
756  thelocalrow = row + boost::get<1>( it->second );
757 
758  //auto row1_entries = boost::unwrap_ref( boost::ref( M_graph->row(theglobalrow).template get<2>() ) );
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 );
761 
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>()];
765 
766  if ( !row2_entries.empty() )
767  {
768 
769  if ( row1_entries.empty() )
770  {
771  // if row is empty then no need to shift the dof in
772  // composite case since the merge in done block-row-wise
773  //row1_entries = row2_entries;
774  if ( col==0 ) row1_entries = row2_entries;
775 
776  else
777  {
778  for ( auto itcol = row2_entries.begin(), encol = row2_entries.end() ; itcol!=encol; ++itcol ) row1_entries.insert( *itcol+col );
779  }
780 
781  }
782 
783  else
784  {
785  // ensure unique sorted ids
786  auto itg = boost::prior( row1_entries.end() );
787  // shift dofs in case of composite spaces
788  std::for_each( row2_entries.begin(), row2_entries.end(),[&]( size_type o )
789  {
790  itg = row1_entries.insert( itg, o+col );
791  } );
792  }
793  }
794  } // for( ; it != en; ++it )
795  }
796 
797  DVLOG(2) << " -- merge_graph (" << row << "," << col << ") in " << tim.elapsed() << "\n";
798  DVLOG(2) << "merge graph for composite bilinear form done\n";
799 }
800 
801 template<typename X1, typename X2, typename RangeItTestType>
802 void
803 Stencil<X1,X2,RangeItTestType>::mergeGraphMPI( size_type test_index, size_type trial_index,
804  DataMap const& mapOnTest, DataMap const& mapOnTrial,
805  graph_ptrtype g )
806 {
807  const int myrank = this->testSpace()->worldComm().globalRank();
808 
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 );
812 
813  typename graph_type::const_iterator it = g->begin();
814  typename graph_type::const_iterator en = g->end();
815  for ( ; it != en; ++it )
816  {
817  size_type theglobalrow = globalDofRowStart + ( it->first - mapOnTest.firstDofGlobalCluster() );
818  const size_type thelocalrow = locdofStart + it->second.get<1>();
819 
820  if (it->second.get<0>()!=g->worldComm().globalRank() )
821  {
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));
826  }
827 
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 );
830 
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>()];
834 
835  if ( !row2_entries.empty() )
836  {
837  for ( auto itcol = row2_entries.begin(), encol = row2_entries.end() ; itcol!=encol; ++itcol )
838  {
839  if (mapOnTrial.dofGlobalClusterIsOnProc(*itcol))
840  {
841  const size_type dofcol = globalDofColStart + (*itcol-mapOnTrial.firstDofGlobalCluster());
842  row1_entries.insert( dofcol );
843  }
844  else
845  {
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);
851  }
852  }
853  }
854 
855  } // for( ; it != en; ++it )
856 
857 }
858 
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 )
862 {
863  VLOG(2) << "computeGraph: deciding whether the mesh are related to optimize the stencil\n";
864  //if ( (is_shared_ptr<typename test_space_type::mesh_ptrtype>::value && is_shared_ptr<typename trial_space_type::mesh_ptrtype>::value ) &&
865  //dynamic_cast<void*>( _M_X1->template mesh<0>().get() ) == dynamic_cast<void*>( _M_X2->template mesh<0>().get() ) )
866  if ( _M_X1->template mesh<0>()->isRelatedTo( _M_X2->template mesh<0>() ) )
867  {
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 >() );
871  }
872  else
873  {
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 >() );
877  }
878 }
879 
880 
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> )
884 {
885  boost::timer t;
886  DVLOG(2) << "compute graph for composite bilinear form with interpolation\n";
887 
888  auto graph = computeGraph( hints, mpl::bool_< ( test_space_type::nSpaces > 1 )>(), mpl::bool_< ( trial_space_type::nSpaces > 1 )>() );
889 
890  DVLOG(2) << "closing graph for composite bilinear form with interpolation done in " << t.elapsed() << "s\n";
891  t.restart();
892  //graph->close();
893  DVLOG(2) << "compute graph for composite bilinear form done in " << t.elapsed() << "s\n";
894 
895  return graph;
896 }
897 
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> )
901 {
902  fusion::for_each( _M_X1->functionSpaces(),
903  detail::compute_graph1<self_type>( this, hints ) );
904  return M_graph;
905 }
906 
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> )
910 {
911  fusion::for_each( _M_X1->functionSpaces(),
912  detail::compute_graph3<self_type,trial_space_type>( this, _M_X2, 0, hints ) );
913  return M_graph;
914 }
915 
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> )
919 {
920  fusion::for_each( _M_X2->functionSpaces(),
921  detail::compute_graph2<self_type,test_space_type>( this, _M_X1, 0, hints ) );
922  return M_graph;
923 }
924 
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> )
928 {
929  boost::timer t;
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 )>() );
934 
935  DVLOG(2) << "closing graph for composite bilinear form with interpolation done in " << t.elapsed() << "s\n";
936  return graph;
937 }
938 
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> )
942 {
943  fusion::for_each( _M_X1->functionSpaces(),
944  detail::compute_graph1<self_type>( this, hints ) );
945  return M_graph;
946 }
947 
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> )
951 {
952  fusion::for_each( _M_X1->functionSpaces(),
953  detail::compute_graph3<self_type,trial_space_type>( this, _M_X2, 0, hints ) );
954  return M_graph;
955 }
956 
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> )
960 {
961  fusion::for_each( _M_X2->functionSpaces(),
962  detail::compute_graph2<self_type,test_space_type>( this, _M_X1, 0, hints ) );
963  return M_graph;
964 }
965 
966 
967 
968 
969 #if 0
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> )
973 {
974  boost::timer t;
975  // Compute the sparsity structure of the global matrix. This can be
976  // fed into a PetscMatrix to allocate exacly the number of nonzeros
977  // necessary to store the matrix. This algorithm should be linear
978  // in the (# of elements)*(# nodes per element)
979  const size_type proc_id = _M_X1->mesh()->comm().rank();
980  const size_type n1_dof_on_proc = _M_X1->nLocalDof();
981  //const size_type n2_dof_on_proc = _M_X2->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 ) );
989 
990  typedef typename mesh_type::element_const_iterator mesh_element_const_iterator;
991 
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 );
994 
995  Feel::Context graph( hints );
996  // If the user did not explicitly specify the DOF coupling
997  // then all the DOFS are coupled to each other. Furthermore,
998  // we can take a shortcut and do this more quickly here. So
999  // we use an if-test.
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";
1003 #if 0
1004 
1005  if ( graph.test ( Pattern::COUPLED ) ||
1006  graph.test ( Pattern::EXTENDED ) )
1007 #else
1008  if ( 1 )
1009 #endif
1010  {
1011  DVLOG(2) << "[computeGraph] test (Pattern::COUPLED || Pattern::EXTENDED) ok\n";
1012  std::vector<size_type>
1013  element_dof1,
1014  element_dof2,
1015  neighbor_dof,
1016  dof_to_add;
1017 
1018  for ( ; elem_it != elem_en; ++elem_it )
1019  {
1020 #if !defined(NDEBUG)
1021  DVLOG(2) << "[Stencil::computePatter] element " << elem_it->id() << " on proc " << elem_it->processId() << "\n";
1022 #endif /* NDEBUG */
1023  mesh_element_type const& elem = *elem_it;
1024 
1025  // Get the global indices of the DOFs with support on this element
1026  element_dof1 = _M_X1->dof()->getIndices( elem.id() );
1027  element_dof2 = _M_X2->dof()->getIndices( elem.id() );
1028 
1029  // We can be more efficient if we sort the element DOFs
1030  // into increasing order
1031  std::sort( element_dof1.begin(), element_dof1.end() );
1032  std::sort( element_dof2.begin(), element_dof2.end() );
1033 
1034  const uint16_type n1_dof_on_element = element_dof1.size();
1035  const uint16_type n2_dof_on_element = element_dof2.size();
1036 
1037  for ( size_type i=0; i<n1_dof_on_element; i++ )
1038  {
1039  const size_type ig1 = element_dof1[i];
1040 
1041  // Only bother if this matrix row will be stored
1042  // on this processor.
1043  //if ((ig1 >= first1_dof_on_proc) &&
1044  //(ig1 <= last1_dof_on_proc))
1045  {
1046  // This is what I mean
1047  // assert ((ig - first_dof_on_proc) >= 0);
1048  // but do the test like this because ig and
1049  // first_dof_on_proc are size_types
1050 #if 0
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" );
1054 #endif
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 );
1057  row.get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1058  row.get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1059  DVLOG(2) << "work with row " << ig1 << " local index " << ig1 - first1_dof_on_proc << "\n";
1060 
1061  // If the row is empty we will add *all* the element DOFs,
1062  // so just do that.
1063  if ( row.get<2>().empty() )
1064  {
1065  row.get<2>() = element_dof2;
1066  }
1067 
1068  else
1069  {
1070  // Build a list of the DOF indices not found in the
1071  // sparsity graph
1072  dof_to_add.clear();
1073 
1074  // Cache iterators. Low will move forward, subsequent
1075  // searches will be on smaller ranges
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() );
1079 
1080  for ( size_type j=0; j<n2_dof_on_element; j++ )
1081  {
1082  const size_type jg = element_dof2[j];
1083  //VLOG(1) << "computeGraph : ig:" << ig1 << ", lig: " << ig1-first1_dof_on_proc << ", jg = " << jg << "\n";
1084 #if 0
1085  // See if jg is in the sorted range
1086  std::pair<std::vector<size_type>::iterator,
1087  std::vector<size_type>::iterator>
1088  pos = std::equal_range ( low, high, jg );
1089 
1090  // Must add jg if it wasn't found
1091  if ( pos.first == pos.second )
1092  dof_to_add.push_back( jg );
1093 
1094  // pos.first is now a valid lower bound for any
1095  // remaining element Dof. (That's why we sorted them.)
1096  // Use it for the next search
1097  low = pos.first;
1098 #else
1099  // See if jg is in the sorted range
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 );
1103 
1104  // Insert jg if it wasn't found
1105  if ( pos.first == pos.second )
1106  dof_to_add.push_back( jg );
1107 
1108 #endif
1109  }
1110 
1111  // Add to the sparsity graph
1112  if ( !dof_to_add.empty() )
1113  {
1114  const size_type old_size = row.get<2>().size();
1115 
1116  row.get<2>().insert ( row.get<2>().end(),
1117  dof_to_add.begin(),
1118  dof_to_add.end() );
1119 
1120  //std::inplace_merge (row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end());
1121  sortSparsityRow ( row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end() );
1122  }
1123 
1124  }
1125 
1126  // Now (possibly) add dof from neighboring elements
1127  //if ( graph.test( Pattern::EXTENDED ) )
1128  for ( uint16_type ms=0; ms < elem.nNeighbors(); ms++ )
1129  {
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;
1133 
1134  if ( neighbor_id != invalid_size_type_value )
1135  //&& neighbor_process_id != proc_id )
1136  {
1137 
1138 #if 0
1139  VLOG(1) << "element id " << elem.id()
1140  << ", element neighbor id " << neighbor_id
1141  << " in proc " << neighbor_process_id << "\n";
1142 #endif
1143  neighbor = boost::addressof( _M_X1->mesh()->element( neighbor_id,
1144  neighbor_process_id ) );
1145 
1146 #if 0
1147  VLOG(1) << "found neighbor of element id " << elem.id()
1148  << ", element neighbor id " << neighbor->id()
1149  << " in proc " << neighbor->processId() << "\n";
1150 #endif
1151 
1152  if ( neighbor_id == neighbor->id() )
1153  {
1154  neighbor_dof = _M_X2->dof()->getIndices( neighbor->id() );
1155 
1156  const size_type n_dof_on_neighbor = neighbor_dof.size();
1157 #if 0
1158 
1159  for ( size_type j=0; j<n_dof_on_neighbor; j++ )
1160  {
1161  DVLOG(2) << "neighbor elem id: " << neighbor->id() << " dof " << neighbor_dof[j] << "\n";
1162  }
1163 
1164  DVLOG(2) << "looking for dof " << ig1 << "\n";
1165 #endif
1166 #if 0
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 );
1170 #else
1171  std::vector<size_type>::iterator posig = std::find( neighbor_dof.begin(), neighbor_dof.end(), ig1 );
1172 #endif
1173 
1174  // Insert jg if it wasn't found
1175  //if (posig.first != posig.second)
1176  if ( posig != neighbor_dof.end() ||
1177  graph.test ( Pattern::EXTENDED ) )
1178 
1179  {
1180  //VLOG(1) << "found element in proc " << neighbor_process_id << " that shares dof\n";
1181  for ( size_type j=0; j<n_dof_on_neighbor; j++ )
1182  {
1183  const size_type jg = neighbor_dof[j];
1184 
1185 #if 0
1186  // See if jg is in the sorted range
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 );
1190 #else
1191  std::vector<size_type>::iterator pos = std::find( row.get<2>().begin(), row.get<2>().end(), jg );
1192 
1193 #endif
1194 
1195  // Insert jg if it wasn't found
1196  if ( pos == row.get<2>().end() )
1197  {
1198  const size_type old_size = row.get<2>().size();
1199  row.get<2>().push_back ( jg );
1200  //std::inplace_merge (row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end());
1201  sortSparsityRow ( row.get<2>().begin(), row.get<2>().begin()+old_size, row.get<2>().end() );
1202  }
1203  }
1204  }
1205  }
1206  }
1207 
1208  } // neighbor graph
1209 
1210 #if 0
1211 
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";
1214 
1215 #endif
1216  } // only dof on proc
1217 
1218  }// dof loop
1219  } // element iterator loop
1220  }
1221 
1222  else
1223  {}
1224 
1225  DVLOG(2) << "[computeGraph<true>] before calling close in " << t.elapsed() << "s\n";
1226  //sparsity_graph->close();
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;
1230 }
1231 #else
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> )
1235 {
1236  boost::timer t;
1237  // Compute the sparsity structure of the global matrix. This can be
1238  // fed into a PetscMatrix to allocate exacly the number of nonzeros
1239  // necessary to store the matrix. This algorithm should be linear
1240  // in the (# of elements)*(# nodes per element)
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();
1244  //const size_type n2_dof_on_proc = _M_X2->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 );
1249 #else // MPI
1250  const size_type proc_id = _M_X1->worldsComm()[0].localRank();
1251  const size_type n1_dof_on_proc = _M_X1->nLocalDof();
1252  //const size_type n2_dof_on_proc = _M_X2->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 );
1257 #endif
1258 
1259  graph_ptrtype sparsity_graph( new graph_type( _M_X1->dof(),_M_X2->dof() ) );
1260 
1261  Feel::Context graph( hints );
1262  if ( graph.test( Pattern::ZERO ) ) { sparsity_graph->zero(); return sparsity_graph; }
1263 
1264  //if (_M_X1->nLocalDofWithoutGhost()==0 && _M_X2->nLocalDofWithoutGhost()==0 ) return sparsity_graph;
1265 
1266  auto elem_it = _M_X1->mesh()->beginElementWithProcessId( _M_X1->mesh()->worldComm().localRank() /*proc_id*/ );
1267  auto elem_en = _M_X1->mesh()->endElementWithProcessId( _M_X1->mesh()->worldComm().localRank() /*proc_id*/ );
1268 
1269  // If the user did not explicitly specify the DOF coupling
1270  // then all the DOFS are coupled to each other. Furthermore,
1271  // we can take a shortcut and do this more quickly here. So
1272  // we use an if-test.
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() ),
1282  neighbor_dof;
1283 
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 )
1288  {
1289  DVLOG(4) << "[Stencil::computePattern] element " << elem_it->id() << " on proc " << elem_it->processId() << "\n";
1290 
1291  const auto & elem = *elem_it;
1292 
1293  auto const domains_eid_set = trialElementId( elem.id(), mpl::int_<nDimDiffBetweenTestTrial>() );
1294 
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 )
1298  {
1299  const size_type domain_eid = *it_trial;
1300 
1301  // Get the global indices of the DOFs with support on this element
1302 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1303  _M_X2->dof()->getIndicesSet( domain_eid, element_dof2 );
1304 #else // MPI
1305  _M_X2->dof()->getIndicesSetOnGlobalCluster( domain_eid, element_dof2 );
1306 #endif
1307  // We can be more efficient if we sort the element DOFs
1308  // into increasing order
1309  //std::sort(element_dof1.begin(), element_dof1.end());
1310  std::sort( element_dof2.begin(), element_dof2.end() );
1311 
1312  //const uint16_type n1_dof_on_element = element_dof1.size();
1313  const uint16_type n1_dof_on_element = _M_X1->dof()->getIndicesSize();
1314  const uint16_type n2_dof_on_element = element_dof2.size();
1315 
1316  for ( size_type i=0; i<n1_dof_on_element; i++ )
1317  //BOOST_FOREACH( auto ig1, _M_X1->dof()->getIndices( elem.id() ) )
1318  {
1319 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
1320  const size_type ig1 = _M_X1->dof()->localToGlobalId( elem.id(), i );
1321 #else // MPI
1322  const size_type ig1 = _M_X1->dof()->mapGlobalProcessToGlobalCluster()[_M_X1->dof()->localToGlobalId( elem.id(), i )];
1323  auto theproc = _M_X1->dof()->procOnGlobalCluster( ig1 );
1324  // numLocal without ghosts ! very important for the graph with petsc
1325  const size_type il1 = _M_X1->dof()->localToGlobalId( elem.id(), i );// ig1 - _M_X1->dof()->firstDofGlobalCluster( theproc );
1326 #endif
1327  //const size_type ig1 = element_dof1[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;
1331 
1332  {
1333  // This is what I mean
1334  // assert ((ig - first_dof_on_proc) >= 0);
1335  // but do the test like this because ig and
1336  // first_dof_on_proc are size_types
1337 #if 0
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" );
1341 #endif
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 );
1345  row.get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1346  row.get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1347 #else // MPI
1348  row.get<0>() = theproc ;
1349  row.get<1>() = il1;
1350 #endif
1351  DVLOG(4) << "work with row " << ig1 << " local index " << ig1 - first1_dof_on_proc << "\n";
1352 
1353  if ( do_less )
1354  {
1355  if ( ncomp1 == ( _M_X2->dof()->nComponents-1 ) )
1356  row.get<2>().insert( element_dof2.begin()+ncomp1*ndofpercomponent2,
1357  element_dof2.end() );
1358 
1359  else
1360  row.get<2>().insert( element_dof2.begin()+ncomp1*ndofpercomponent2,
1361  element_dof2.begin()+( ncomp1+1 )*ndofpercomponent2 );
1362  }
1363 
1364  else
1365  {
1366  row.get<2>().insert( element_dof2.begin(), element_dof2.end() );
1367  }
1368 
1369  // Now (possibly) add dof from neighboring elements
1370  if ( graph.test( Pattern::EXTENDED ) )
1371  {
1372  for ( uint16_type ms=0; ms < elem.nNeighbors(); ms++ )
1373  {
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;
1377 
1378  // warning ! the last condition is a temporary solution
1379  if ( neighbor_id != invalid_size_type_value
1380  && neighbor_process_id == proc_id )
1381  {
1382 
1383  neighbor = boost::addressof( _M_X1->mesh()->element( neighbor_id,
1384  neighbor_process_id ) );
1385 
1386  if ( neighbor_id == neighbor->id() )
1387  {
1388  //neighbor_dof = _M_X2->dof()->getIndices( neighbor->id() );
1389  neighbor_dof = _M_X2->dof()->getIndicesOnGlobalCluster( neighbor->id() );
1390 
1391  if ( do_less )
1392  {
1393  if ( ncomp1 == ( _M_X2->dof()->nComponents-1 ) )
1394  row.get<2>().insert( neighbor_dof.begin()+ncomp1*ndofpercomponent2,
1395  neighbor_dof.end() );
1396 
1397  else
1398  row.get<2>().insert( neighbor_dof.begin()+ncomp1*ndofpercomponent2,
1399  neighbor_dof.begin()+( ncomp1+1 )*ndofpercomponent2 );
1400 
1401  }
1402 
1403  else
1404  {
1405  row.get<2>().insert( neighbor_dof.begin(), neighbor_dof.end() );
1406  }
1407 
1408  } // neighbor_id
1409  }
1410  } // neighbor graph
1411  } // if ( graph.test( Pattern::EXTENDED ) )
1412  } // only dof on proc
1413  } // dof loop
1414  } // trial id loop
1415  } // element iterator loop
1416 
1417  DVLOG(2)<< "[computeGraph<true>] before calling close in " << t.elapsed() << "s\n";
1418  //sparsity_graph->close();
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;
1422 }
1423 #endif
1424 namespace detail
1425 {
1426 #if 0
1427 template<typename MeshType>
1428 struct gmcDefStencil
1429 {
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;
1434 };
1435 #endif
1436 template<typename EltType>
1437 struct gmcDefStencil
1438 {
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;
1443 };
1444 
1445 template<typename FaceType>
1446 struct gmcDefFaceStencil
1447 {
1448  typedef typename boost::remove_reference<FaceType>::type FaceType1;
1449  typedef typename boost::remove_const<FaceType1>::type FaceType2;
1450 
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;
1456 };
1457 
1458 template<typename ImType,typename EltType>
1459 typename gmcDefStencil<EltType>::gmc_ptrtype
1460 gmcStencil( mpl::size_t<MESH_ELEMENTS> , EltType const& elem )
1461 {
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;
1466 
1467  ImType theim;
1468  pc_ptrtype geopc( new pc_type( elem.gm(), theim.points() ) );
1469  gmc_ptrtype gmc( new gmc_type( elem.gm(), elem, geopc ) );
1470  return gmc;
1471 }
1472 
1473 template<typename ImType,typename FaceType>
1474 typename gmcDefFaceStencil<FaceType>::gmc_ptrtype
1475 gmcStencil( mpl::size_t<MESH_FACES> , FaceType const& theface )
1476 {
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;
1481 
1482  typedef typename QuadMapped<ImType>::permutation_type permutation_type;
1483  typedef typename QuadMapped<ImType>::permutation_points_type permutation_points_type;
1484 
1485  ImType theim;
1486  QuadMapped<ImType> qm;
1487  permutation_points_type ppts( qm(theim) );
1488 
1489  std::vector<std::map<permutation_type, pc_ptrtype> > __geopc( theim.nFaces() );
1490  for ( uint16_type __f = 0; __f < theim.nFaces(); ++__f )
1491  {
1492  for ( permutation_type __p( permutation_type::IDENTITY );
1493  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1494  {
1495  __geopc[__f][__p] = pc_ptrtype( new pc_type( theface.element( 0 ).gm(), ppts[__f].find( __p )->second ) );
1496  }
1497  }
1498 
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 ),
1502  __geopc,
1503  __face_id_in_elt_0 ) );
1504  return gmc;
1505 }
1506 template<typename EltType>
1507 void
1508 gmcUpdateStencil( mpl::size_t<MESH_ELEMENTS> , EltType const& elem, typename gmcDefStencil<EltType>::gmc_ptrtype &gmc )
1509 {
1510  gmc->update( elem );
1511 }
1512 template<typename FaceType>
1513 void
1514 gmcUpdateStencil( mpl::size_t<MESH_FACES> , FaceType const& theface, typename gmcDefFaceStencil<FaceType>::gmc_ptrtype &gmc )
1515 {
1516  gmc->update( theface.element( 0 ), theface.pos_first() );
1517 }
1518 template<typename EltType>
1519 size_type
1520 idEltStencil( mpl::size_t<MESH_ELEMENTS> , EltType const& elem )
1521 {
1522  return elem.id();
1523 }
1524 template<typename FaceType>
1525 size_type
1526 idEltStencil( mpl::size_t<MESH_FACES> , FaceType const& theface )
1527 {
1528  return theface.element( 0 ).id();
1529 }
1530 
1531 } // namespace detail
1532 
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> )
1536 {
1537  //std::cout << "\n start graphInterp "<< std::endl;
1538 
1539  typedef mpl::int_<20/*50*/> order_1d_type;
1540  typedef mpl::int_<12/*20*/> order_2d_type;
1541  typedef mpl::int_<8/*10*/> order_3d_type;
1542 
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> >,
1546  order_1d_type,
1547  typename mpl::if_<mpl::equal_to<mpl::int_<test_mesh_type::nDim>,mpl::int_<2> >,
1548  order_2d_type,
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;
1554 
1555  typedef typename test_mesh_type::Localization::matrix_node_type matrix_node_type;
1556 
1557  //-----------------------------------------------------------------------//
1558 
1559  const size_type proc_id = _M_X1->mesh()->comm().rank();
1560  const size_type n1_dof_on_proc = _M_X1->nLocalDof();
1561  //const size_type n2_dof_on_proc = _M_X2->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 );
1566 
1567  graph_ptrtype sparsity_graph( new graph_type( _M_X1->dof(),_M_X2->dof() ) );
1568 
1569 #if 0
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>();
1574 #endif
1575 
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>();
1580 
1581  if ( elem_it==elem_en ) return sparsity_graph;
1582 
1583  //-----------------------------------------------------------------------//
1584  // init localisation tools
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);
1592 
1593  //locTool->kdtree()->nbNearNeighbor(_M_X2->mesh()->numElements() );
1594  auto locToolForXh1 = _M_X1->mesh()->tool_localization();
1595  locToolForXh1->updateForUse();
1596  bool doExtrapolationAtStartXh1 = locToolForXh1->doExtrapolation();
1597  if (doExtrapolationAtStartXh1) locToolForXh1->setExtrapolation( false );
1598 
1599 
1600  matrix_node_type ptsReal( elem_it->vertices().size1(), 1 );
1601  size_type IdEltInXh2 = invalid_size_type_value;
1602  //node_type trialNodeRef,testNodeRef;
1603 
1604 #if FEELPP_EXPORT_GRAPH
1605  std::map<size_type,std::list<size_type> > mapBetweenMeshes;
1606 #endif
1607 
1608  std::vector<size_type> element_dof1_range, element_dof1, element_dof2;
1609  std::set<size_type> neighLocalizedInXh1;
1610 
1611  std::set<size_type > listTup;
1612  auto gmc = Feel::detail::gmcStencil<theim_type>( iDimRange,*elem_it );
1613 
1614 
1615  //-----------------------------------------------------------------------//
1616 
1617  if ( _M_X1->nDof()>1 )
1618  {
1619  for ( ; elem_it != elem_en; ++elem_it )
1620  {
1621  auto const& elem = *elem_it;
1622 
1623  // Get the global indices of the DOFs with support on this element
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();
1628 
1629  std::vector<boost::tuple<bool,size_type> > hasFinds( n1_dof_on_element_range,boost::make_tuple( false,invalid_size_type_value ) );
1630 
1631  for ( size_type i=0; i<n1_dof_on_element_range; i++ )
1632  {
1633  const size_type ig1 = element_dof1_range[i];
1634  auto const ptRealDof = boost::get<0>( _M_X1->dof()->dofPoint( ig1 ) );
1635 
1636  ublas::column(ptsReal,0 ) = ptRealDof;
1637  if (notUseOptLocTrial) IdEltInXh2=invalid_size_type_value;
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];
1641 
1642  listTup.clear();
1643 
1644  if ( hasFind )
1645  {
1646  listTup.insert( IdEltInXh2 );
1647  hasFinds[i] = boost::make_tuple( true,IdEltInXh2 );
1648  // maybe is on boundary->more elts
1649  auto const& geoelt2 = _M_X2->mesh()->element( IdEltInXh2 );
1650 
1651  std::vector<size_type> neighbor_ids;
1652  for ( uint16_type ms=0; ms < geoelt2.nNeighbors(); ms++ )
1653  {
1654  size_type neighbor_id = geoelt2.neighbor( ms ).first;
1655 
1656  if ( neighbor_id!=invalid_size_type_value )
1657  neighbor_ids.push_back( neighbor_id );
1658  }
1659 
1660  auto resIsIn = locToolForXh2->isIn( neighbor_ids,ptRealDof );
1661 
1662  uint16_type cpt=0;
1663  for ( auto it=resIsIn.template get<1>().begin(),en=resIsIn.template get<1>().end(); it<en; ++it,++cpt )
1664  {
1665  if ( *it )
1666  {
1667  listTup.insert( neighbor_ids[cpt] );
1668  }
1669  }
1670 
1671  auto res_it = listTup.begin();
1672  auto res_en = listTup.end();
1673  for ( ; res_it != res_en ; ++res_it )
1674  {
1675 #if FEELPP_EXPORT_GRAPH
1676  mapBetweenMeshes[*res_it].push_back( elem.id() );
1677  //std::cout << " test id " << *res_it << " trial id " << elem.id() << std::endl;
1678 #endif
1679  // not efficient but sometimes necessary
1680  for ( size_type ii=0; ii<n1_dof_on_element; ii++ )
1681  {
1682  const size_type ig1ongraph = element_dof1[ii];
1683 
1684  element_dof2 = _M_X2->dof()->getIndices( *res_it );
1685 
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 );
1688  row.template get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1689  row.template get<1>() = is_on_proc?ig1ongraph - first1_dof_on_proc:invalid_size_type_value;
1690  //if ( do_less ) {}
1691  row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1692  }
1693  }//res
1694  } // if (hasFind)
1695 
1696  else
1697  {
1698 #if 0
1699  //std::cout << "\n not find"<<std::endl;
1700  // row empty
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 );
1703  row.template get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1704  row.template get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1705  row.template get<2>().clear();
1706 #endif
1707  }
1708  } // for (size_type i=0; i<n1_dof_on_element_range; i++)
1709 
1710  uint16_type nbFind = hasFinds.size();
1711  bool doQ=false;
1712  size_type id1,id2;
1713 
1714  for ( uint16_type nF = 0 ; nF<nbFind && !doQ ; ++nF )
1715  if ( hasFinds[nF].template get<0>() )
1716  {
1717  id1=hasFinds[nF].template get<1>();
1718 
1719  for ( uint16_type nF2 = 0 ; nF2<nbFind && !doQ ; ++nF2 )
1720  {
1721  if ( hasFinds[nF2].template get<0>() )
1722  {
1723  id2=hasFinds[nF2].template get<1>();
1724 
1725  if ( id1!=id2 )
1726  doQ=true; //if( !_M_X2->mesh()->element(id1).isNeighbor(_M_X2->mesh()->element(id2) ) ) doQ=true;
1727  }
1728  }
1729  }
1730 
1731  if ( doQ )
1732  {
1733  Feel::detail::gmcUpdateStencil( iDimRange,elem,gmc );
1734 
1735  for ( int q = 0; q < gmc->nPoints(); ++ q )
1736  {
1737  ublas::column(ptsReal,0 ) = gmc->xReal( q );
1738  //auto const resQuad = locToolForXh2->run_analysis(ptsReal,IdEltInXh2,elem_it->vertices(),mpl::int_<0>());
1739  auto const resQuad = locToolForXh2->searchElement( gmc->xReal( q ) );
1740 
1741  if ( resQuad.template get<0>() )
1742  {
1743  IdEltInXh2 = resQuad.template get<1>();
1744 #if FEELPP_EXPORT_GRAPH
1745  mapBetweenMeshes[IdEltInXh2].push_back( elem.id() );
1746 #endif
1747  element_dof2 = _M_X2->dof()->getIndices( IdEltInXh2 );
1748 
1749  for ( size_type i=0; i<n1_dof_on_element; i++ )
1750  {
1751  const size_type ig1 = element_dof1[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 );
1754  row.template get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1755  row.template get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1756  //if ( do_less ) {}
1757  row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1758  }
1759  }
1760  } // for ( int q = 0; q < gmc->nPoints(); ++ q )
1761  } // if (doQ)
1762  } // for ( ; elem_it ... )
1763  } //Xh1->nof >1
1764 
1765  else //case Xh1->nof == 1
1766  {
1767  for ( ; elem_it != elem_en; ++elem_it )
1768  {
1769  auto const& elem = *elem_it;
1770 
1771  // Get the global indices of the DOFs with support on this element
1772  element_dof1 = _M_X1->dof()->getIndices( Feel::detail::idEltStencil( iDimRange,elem ) );
1773 
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 ) );
1776 
1777  for ( size_type i=0; i<nPtGeo; i++ )
1778  {
1779  const size_type ig1 = element_dof1[0];
1780  //const int ndofpercomponent1 = n1_dof_on_element / _M_X1->dof()->nComponents;
1781  //const int ncomp1 = i / ndofpercomponent1;
1782  //const int ndofpercomponent2 = n2_dof_on_element / _M_X2->dof()->nComponents;
1783 
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 );
1786 
1787  //auto res = locTool->searchElements(ublas::column( elem.G(),i )/*ptReal*/);
1788  //auto hasFind = res.template get<0>();
1789 
1790  auto resTemp = locToolForXh2->searchElement( ublas::column( elem.G(),i ) );
1791  bool hasFind = resTemp.template get<0>();
1792  std::set<size_type > listTup;
1793 
1794  if ( hasFind )
1795  {
1796  listTup.insert( resTemp.template get<1>() );
1797  hasFinds[i] = boost::make_tuple( true,resTemp.template get<1>() );
1798  // maybe is on boundary->more elts
1799  //size_type idElt1 = elem.id();
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() ); //neighbor_ids.clear();//(geoelt2.nNeighbors());
1803 
1804  for ( uint16_type ms=0; ms < geoelt2.nNeighbors(); ms++ )
1805  {
1806  size_type neighbor_id = geoelt2.neighbor( ms ).first;
1807 
1808  if ( neighbor_id!=invalid_size_type_value ) neighbor_ids.push_back( neighbor_id );
1809 
1810  //neighbor_ids[ms]=neighbor_id;
1811  }
1812 
1813  auto resIsIn = locToolForXh2->isIn( neighbor_ids,ublas::column( elem.G(),i ) );
1814  uint16_type cpt=0;
1815 
1816  for ( auto it=resIsIn.template get<1>().begin(),en=resIsIn.template get<1>().end(); it<en; ++it,++cpt )
1817  {
1818  if ( *it )
1819  {
1820  listTup.insert( neighbor_ids[cpt] );
1821  }
1822 
1823  else if ( test_mesh_type::nDim==trial_mesh_type::nDim ) // pt n'est pas chez le voison, peut-etre les pts G()(si dim=dim!!!)
1824  {
1825  for ( auto it_neigh = neighbor_ids.begin(), en_neigh = neighbor_ids.end() ; it_neigh != en_neigh ; ++it_neigh )
1826  {
1827  if ( ( sparsity_graph->row( ig1 ) ).template get<2>().find( *it_neigh )==( sparsity_graph->row( ig1 ) ).template get<2>().end() )
1828  {
1829  if ( neighLocalizedInXh1.find( *it_neigh )==neighLocalizedInXh1.end() )
1830  {
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();
1835 
1836  for ( size_type iii=0; iii<nPtGeoBis && !findNeih ; iii++ )
1837  {
1838  auto resBis = locToolForXh1->searchElement( ublas::column( geoeltNEW.G(),iii ) );
1839 
1840  if ( resBis.template get<0>() )
1841  {
1842  listTup.insert( *it_neigh );
1843  findNeih=true;
1844  }
1845  }
1846  }
1847  }
1848  }
1849  } //else
1850  } // for (auto it=resIsIn...
1851 
1852  auto res_it = listTup.begin();
1853  auto res_en = listTup.end();
1854  for ( ; res_it != res_en ; ++res_it )
1855  {
1856 #if FEELPP_EXPORT_GRAPH
1857  //mapBetweenMeshes[elem.id()].push_back(*res_it);
1858  mapBetweenMeshes[*res_it/*->get<0>()*/].push_back( elem.id() );
1859 #endif
1860 
1861 
1862  element_dof2 = _M_X2->dof()->getIndices( *res_it/*->get<0>()*/ );
1863 
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 );
1866  row.template get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1867  row.template get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1868  //if ( do_less ) {}
1869  row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1870  }
1871  } // hasFind
1872  } // for (size_type i=0; i<nPtGeo; i++)
1873 
1874  uint16_type nbFind = hasFinds.size();
1875  bool doQ=false;
1876  size_type id1,id2;
1877 
1878  for ( uint16_type nF = 0 ; nF<nbFind && !doQ ; ++nF )
1879  if ( hasFinds[nF].template get<0>() )
1880  {
1881  id1=hasFinds[nF].template get<1>();
1882 
1883  for ( uint16_type nF2 = 0 ; nF2<nbFind && !doQ ; ++nF2 )
1884  {
1885  if ( hasFinds[nF2].template get<0>() )
1886  {
1887  id2=hasFinds[nF2].template get<1>();
1888 
1889  if ( id1!=id2 )
1890  doQ=true; //if( !_M_X2->mesh()->element(id1).isNeighbor(_M_X2->mesh()->element(id2) ) ) doQ=true;
1891  }
1892  }
1893  }
1894 
1895  if ( doQ )
1896  {
1897  Feel::detail::gmcUpdateStencil(iDimRange,elem,gmc);
1898 
1899  for ( int q = 0; q < gmc->nPoints(); ++ q )
1900  {
1901  auto resQuad = locToolForXh2->searchElement( gmc->xReal( q ) );
1902 
1903  if ( resQuad.template get<0>() )
1904  {
1905 #if FEELPP_EXPORT_GRAPH
1906  mapBetweenMeshes[resQuad.template get<1>()].push_back( elem.id() );
1907 #endif
1908  element_dof2 = _M_X2->dof()->getIndices( resQuad.template get<1>() );
1909  //for (size_type i=0; i<n1_dof_on_element; i++)
1910  // {
1911  const size_type ig1 = element_dof1[0];
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 );
1914  row.template get<0>() = is_on_proc?proc_id:invalid_size_type_value;
1915  row.template get<1>() = is_on_proc?ig1 - first1_dof_on_proc:invalid_size_type_value;
1916  //if ( do_less ) {}
1917  row.template get<2>().insert( element_dof2.begin(), element_dof2.end() );
1918  // }
1919  }
1920  }
1921  } // if (doQ)
1922 
1923 
1924  } // for ( ; elem_it != elem_en; ++elem_it)
1925  } // M_X1->nDof==1
1926 
1927  if (doExtrapolationAtStartXh1) locToolForXh1->setExtrapolation( true );
1928  if (doExtrapolationAtStartXh2) locToolForXh2->setExtrapolation( true );
1929  if (notUseOptLocTrial) locToolForXh2->kdtree()->nbNearNeighbor(nbNearNeighborAtStartTrial);
1930 
1931  //locTool->kdtree()->nbNearNeighbor( 15 );
1932 
1933  //sparsity_graph->close();
1934 
1935 #if FEELPP_EXPORT_GRAPH
1936 #if 0
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 );
1942 #else
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 );
1948 #endif
1949  auto graphProj = spaceGraphProj->element();
1950  graphProj.zero();
1951 
1952  for ( ; elem_itt != elem_ent; ++elem_itt )
1953  {
1954  if ( mapBetweenMeshes.find( elem_itt->id() ) != mapBetweenMeshes.end() )
1955  {
1956  if ( mapBetweenMeshes.find( elem_itt->id() )->second.size()>0 )
1957  {
1958  element_dof1 = spaceGraphProj->dof()->getIndices( elem_itt->id() );
1959  const uint16_type n1_dof_on_element = element_dof1.size();
1960  //std::cout << "\nn1_dof_on_element " << n1_dof_on_element << std::endl;
1961  for ( uint16_type i=0; i<n1_dof_on_element; i++ )
1962  {
1963  const size_type ig1 = element_dof1[i];
1964  //std::cout << "ig1 "<<ig1 << " graphProj.nDof() "<< graphProj.nDof() <<std::endl;
1965  graphProj( ig1 ) = 1; //graphProj.set(0, 1. );
1966  }
1967  }
1968  }
1969  }
1970 
1971  auto exporter = boost::shared_ptr<Feel::Exporter<mesh_export_type> >( Feel::Exporter<mesh_export_type>::New( "ensight", "ExportGraph" ) );
1972  exporter->step( 0 )->setMesh( graphProj.mesh() );
1973  exporter->step( 0 )->add( "graphProj", graphProj );
1974  exporter->save();
1975 #endif
1976  //std::cout << "\n finish graphInterp "<< std::endl;
1977  return sparsity_graph;
1978 }
1979 
1980 
1981 }
1982 #endif
1983 #endif /* __galerkingraph_H */

Generated on Sun Oct 20 2013 08:25:05 for Feel++ by doxygen 1.8.4