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
operatorinterpolation.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2008-02-01
7 
8  Copyright (C) 2008-2012 Universite 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 3.0 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 */
30 #ifndef __OperatorInterpolation_H
31 #define __OperatorInterpolation_H 1
32 
34 
35 namespace Feel
36 {
37 
38 class InterpolationTypeBase
39 {
40 public :
41  InterpolationTypeBase( bool useComm=true, bool compAreSamePt=true, bool onlyLocalizeOnBoundary=false, int nbNearNeighborInKdTree=15 )
42  :
43  M_searchWithCommunication( useComm ),
44  M_componentsAreSamePoint( compAreSamePt ),
45  M_onlyLocalizeOnBoundary( onlyLocalizeOnBoundary ),
46  M_nbNearNeighborInKdTree( nbNearNeighborInKdTree )
47  {}
48 
49  InterpolationTypeBase( InterpolationTypeBase const& a)
50  :
51  M_searchWithCommunication( a.M_searchWithCommunication ),
52  M_componentsAreSamePoint( a.M_componentsAreSamePoint ),
53  M_onlyLocalizeOnBoundary( a.M_onlyLocalizeOnBoundary ),
54  M_nbNearNeighborInKdTree( a.M_nbNearNeighborInKdTree )
55  {}
56 
57  bool searchWithCommunication() const { return M_searchWithCommunication; }
58  bool componentsAreSamePoint() const { return M_componentsAreSamePoint; }
59  bool onlyLocalizeOnBoundary() const { return M_onlyLocalizeOnBoundary; }
60  int nbNearNeighborInKdTree() const { return M_nbNearNeighborInKdTree; }
61 
62 private :
63 
64  bool M_searchWithCommunication;
65  bool M_componentsAreSamePoint;
66  bool M_onlyLocalizeOnBoundary;
67  int M_nbNearNeighborInKdTree;
68 };
69 
70 struct InterpolationNonConforme : public InterpolationTypeBase
71 {
72  static const uint16_type value=0;
73 
74  InterpolationNonConforme( bool useComm=true, bool compAreSamePt=true, bool onlyLocalizeOnBoundary=false, int nbNearNeighborInKdTree=15 )
75  :
76  InterpolationTypeBase(useComm,compAreSamePt, onlyLocalizeOnBoundary,nbNearNeighborInKdTree)
77  {}
78 };
79 
80 struct InterpolationConforme : public InterpolationTypeBase
81 {
82  static const uint16_type value=1;
83 
84  InterpolationConforme(bool useComm=true, bool compAreSamePt=true, bool onlyLocalizeOnBoundary=false, int nbNearNeighborInKdTree=15 )
85  :
86  InterpolationTypeBase(useComm,compAreSamePt,onlyLocalizeOnBoundary,nbNearNeighborInKdTree)
87  {}
88 };
89 
90 
91 namespace detailsup
92 {
93 
94 template < typename EltType >
96 idElt( EltType & elt,mpl::size_t<MESH_ELEMENTS> )
97 {
98  return elt.id();
99 }
100 
101 template < typename EltType >
102 size_type
103 idElt( EltType & elt,mpl::size_t<MESH_FACES> )
104 {
105  return elt.element0().id();
106 }
107 
108 
109 } //detailsup
110 
111 
120 template<typename DomainSpaceType,
121  typename ImageSpaceType,
122  typename IteratorRange = boost::tuple<mpl::size_t<MESH_ELEMENTS>,
123  typename MeshTraits<typename ImageSpaceType::mesh_type>::element_const_iterator,
124  typename MeshTraits<typename ImageSpaceType::mesh_type>::element_const_iterator>,
125  typename InterpType = InterpolationNonConforme >
126 class OperatorInterpolation : public OperatorLinear<DomainSpaceType, ImageSpaceType >
127 {
129 public:
130 
131 
135 
136  /*
137  * domain
138  */
139  typedef typename super::domain_space_type domain_space_type;
140  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
141  typedef typename domain_space_type::mesh_type domain_mesh_type;
142  typedef typename domain_mesh_type::element_type domain_geoelement_type;
143  typedef typename domain_mesh_type::element_iterator domain_mesh_element_iterator;
144 
145 
146  typedef typename super::backend_ptrtype backend_ptrtype;
147 
148 
149 
150  /*
151  * image
152  */
153  typedef typename super::dual_image_space_type dual_image_space_type;
154  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
155  typedef typename dual_image_space_type::value_type value_type;
156  typedef typename dual_image_space_type::mesh_type image_mesh_type;
157  typedef typename image_mesh_type::element_type image_geoelement_type;
158  typedef typename image_mesh_type::element_iterator image_mesh_element_iterator;
159 
160  // geometric mapping context
161  typedef typename image_mesh_type::gm_type image_gm_type;
162  typedef typename image_mesh_type::gm_ptrtype image_gm_ptrtype;
163  typedef typename image_mesh_type::template gmc<vm::POINT>::type image_gmc_type;
164  typedef typename image_mesh_type::template gmc<vm::POINT>::ptrtype image_gmc_ptrtype;
165 
166  // dof
167  typedef typename dual_image_space_type::dof_type dof_type;
168 
169  // basis
170  typedef typename dual_image_space_type::basis_type image_basis_type;
171  typedef typename domain_space_type::basis_type domain_basis_type;
172 
173  typedef typename boost::tuples::template element<0, IteratorRange>::type idim_type;
174  typedef typename boost::tuples::template element<1, IteratorRange>::type iterator_type;
175  typedef IteratorRange range_iterator;
176 
177 
178  static const uint16_type nLocalDofInDualImageElt = mpl::if_< mpl::equal_to< idim_type ,mpl::size_t<MESH_ELEMENTS> >,
179  mpl::int_< image_basis_type::nLocalDof > ,
180  mpl::int_< image_mesh_type::face_type::numVertices*dual_image_space_type::fe_type::nDofPerVertex +
181  image_mesh_type::face_type::numEdges*dual_image_space_type::fe_type::nDofPerEdge +
182  image_mesh_type::face_type::numFaces*dual_image_space_type::fe_type::nDofPerFace > >::type::value;
183 
184  // type conforme or non conforme
185  typedef InterpType interpolation_type;
186 
187  // matrix graph
188  typedef GraphCSR graph_type;
189  typedef boost::shared_ptr<graph_type> graph_ptrtype;
190 
191  // node type
192  typedef typename matrix_node<typename image_mesh_type::value_type>::type matrix_node_type;
193 
194 
196 
200 
205 
211  OperatorInterpolation( domain_space_ptrtype const& domainspace,
212  dual_image_space_ptrtype const& imagespace,
213  backend_ptrtype const& backend,
214  InterpType const& interptype,
215  bool ddmethod=false);
216 
217  OperatorInterpolation( domain_space_ptrtype const& domainspace,
218  dual_image_space_ptrtype const& imagespace,
219  IteratorRange const& r,
220  backend_ptrtype const& backend,
221  InterpType const& interptype,
222  bool ddmethod=false);
223 
224  OperatorInterpolation( domain_space_ptrtype const& domainspace,
225  dual_image_space_ptrtype const& imagespace,
226  std::list<IteratorRange> const& r,
227  backend_ptrtype const& backend,
228  InterpType const& interptype,
229  bool ddmethod=false);
230 
231 
236  :
237  super( oi ),
238  M_listRange( oi.M_listRange ),
239  M_WorldCommFusion( oi.M_WorldCommFusion ),
240  M_interptype( oi.M_interptype )
241  {}
242 
244 
246 
250 
251 
253 
257 
258  WorldComm const& worldCommFusion() const { return M_WorldCommFusion; }
259 
260  InterpType const& interpolationType() const { return M_interptype; }
261 
262  bool isDomainMeshRelatedToImageMesh() const { return this->domainSpace()->mesh()->isSubMeshFrom( this->dualImageSpace()->mesh() ); }
263 
264  bool isImageMeshRelatedToDomainMesh() const { return this->dualImageSpace()->mesh()->isSubMeshFrom( this->domainSpace()->mesh() ); }
265 
267 
271 
272 
274 
278 
279 
281 
282 
283 
284 protected:
285 
286  virtual void update();
287 
288 private:
289 
290  void updateSameMesh();
291  void updateNoRelationMesh();
292 #if defined(FEELPP_ENABLE_MPI_MODE)
293  void updateNoRelationMeshMPI();
294  void updateNoRelationMeshMPI_run(bool buildNonZeroMatrix=true);
295 
296  typedef std::vector<std::list<boost::tuple<int,
297  typename image_mesh_type::node_type,
298  typename image_mesh_type::node_type,
299  std::vector<std::pair<size_type,size_type> >, // col gdof and gdofInGlobalCluster
300  uint16_type // comp
301  > > > extrapolation_memory_type;
302 
303  // point distrubition
304  boost::tuple<std::vector< std::vector<size_type> >,
305  std::vector< std::vector<uint16_type> >,
306  std::vector<std::vector<typename image_mesh_type::node_type> >,
307  std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > >
308  updateNoRelationMeshMPI_pointDistribution(const std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
309  std::vector<std::set<size_type> > & dof_searchWithProc);
310 
311  // search in my world
312  std::list<boost::tuple<size_type,uint16_type> >
313  updateNoRelationMeshMPI_upWithMyWorld(const std::vector< std::vector<size_type> > & memmapGdof,
314  const std::vector< std::vector<uint16_type> > & memmapComp,
315  const std::vector<std::vector<typename image_mesh_type::node_type> > & pointsSearched,
316  const std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > & memmap_vertices,
317  graph_ptrtype & sparsity_graph,
318  std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
319  std::vector<std::map<size_type,size_type> > & memory_col_globalProcessToGlobalCluster,
320  std::vector<std::set<size_type> > & dof_searchWithProc,
321  bool extrapolation_mode,
322  extrapolation_memory_type & dof_extrapolationData);
323 
324  // search in other world (MPI communication)
325  std::list<boost::tuple<size_type,uint16_type> >
326  updateNoRelationMeshMPI_upWithOtherWorld(const std::vector< std::vector<size_type> > & memmapGdof,
327  const std::vector< std::vector<uint16_type> > & memmapComp,
328  const std::vector<std::vector<typename image_mesh_type::node_type> > & pointsSearched,
329  const std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > & memmap_vertices,
330  graph_ptrtype & sparsity_graph,
331  std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
332  std::vector<std::map<size_type,size_type> > & memory_col_globalProcessToGlobalCluster,
333  std::vector<std::set<size_type> > & dof_searchWithProc,
334  bool extrapolation_mode,
335  extrapolation_memory_type & dof_extrapolationData);
336 #endif // MPI_MODE
337 
338  std::list<range_iterator> M_listRange;
339  WorldComm M_WorldCommFusion;
340  InterpType M_interptype;
341 };
342 
343 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
345  dual_image_space_ptrtype const& imagespace,
346  backend_ptrtype const& backend,
347  InterpType const& interptype,
348  bool ddmethod )
349  :
350  super( domainspace, imagespace, backend, false ),
351  M_listRange(),
352  M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
353  M_interptype(interptype)
354 {
355  M_listRange.push_back( elements( imagespace->mesh() ) );
356  update();
357 }
358 
359 
360 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
362  dual_image_space_ptrtype const& imagespace,
363  IteratorRange const& r,
364  backend_ptrtype const& backend,
365  InterpType const& interptype,
366  bool ddmethod )
367  :
368  super( domainspace, imagespace, backend, false ),
369  M_listRange(),
370  M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
371  M_interptype(interptype)
372 {
373  M_listRange.push_back( r );
374  update();
375 }
376 
377 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
379  dual_image_space_ptrtype const& imagespace,
380  std::list<IteratorRange> const& r,
381  backend_ptrtype const& backend,
382  InterpType const& interptype,
383  bool ddmethod )
384  :
385  super( domainspace, imagespace, backend, false ),
386  M_listRange( r ),
387  M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
388  M_interptype(interptype)
389 {
390  update();
391 }
392 
393 
394 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
395 void
396 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::update()
397 {
398  if ( this->dualImageSpace()->mesh()->numElements() == 0 )
399  {
400  //std::cout << "OperatorInterpolation : update nothing!" << std::endl;
401  this->matPtr() = this->backend()->newZeroMatrix( this->domainSpace()->dofOnOff(),
402  this->dualImageSpace()->dofOn() );
403  return;
404  }
405 
406  // if same mesh but not same function space (e.g. different polynomial
407  // order, different basis) or if the image of domain mesh are related to
408  // each other through an extraction (one of them is the sub mesh of the
409  // other)
410  if ( this->dualImageSpace()->mesh()->isRelatedTo( this->domainSpace()->mesh() ) &&
411  boost::is_same<domain_mesh_type,image_mesh_type>::type::value // warning TODO in this case
412  )
413  {
414  VLOG(2) << "OperatorInterpolation: use same mesh\n";
415  VLOG(2) << "isDomainMeshRelatedToImageMesh: " << isDomainMeshRelatedToImageMesh() << "\n";
416  VLOG(2) << "isImageMeshRelatedToDomainMesh: " << isImageMeshRelatedToDomainMesh() << "\n";
417  this->updateSameMesh();
418  }
419  else // no relation between meshes
420  {
421 #if defined(FEELPP_ENABLE_MPI_MODE)
422  if ( this->dualImageSpace()->worldComm().localSize() > 1 ||
423  this->domainSpace()->worldComm().localSize() > 1 )
424  this->updateNoRelationMeshMPI();
425 
426  else
427  this->updateNoRelationMesh();
428 #else
429  this->updateNoRelationMesh();
430 #endif
431  }
432 
433  // close matrix after build
434  this->mat().close();
435 }
436 
437 //-----------------------------------------------------------------------------------------------------------------//
438 //-----------------------------------------------------------------------------------------------------------------//
439 //-----------------------------------------------------------------------------------------------------------------//
440 
441 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
442 void
443 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateSameMesh()
444 {
445  //std::cout << "OperatorInterpolation::updateSameMesh start " << std::endl;
446 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
447  const size_type proc_id = this->dualImageSpace()->mesh()->worldComm().localRank();
448  const size_type nrow_dof_on_proc = this->dualImageSpace()->nLocalDof();
449  const size_type firstrow_dof_on_proc = this->dualImageSpace()->dof()->firstDof( proc_id );
450  const size_type lastrow_dof_on_proc = this->dualImageSpace()->dof()->lastDof( proc_id );
451  const size_type firstcol_dof_on_proc = this->domainSpace()->dof()->firstDof( proc_id );
452  const size_type lastcol_dof_on_proc = this->domainSpace()->dof()->lastDof( proc_id );
453 #else
454  const size_type proc_id = this->dualImageSpace()->worldsComm()[0].localRank();
455  const size_type nrow_dof_on_proc = this->dualImageSpace()->nLocalDof();
456  const size_type firstrow_dof_on_proc = this->dualImageSpace()->dof()->firstDofGlobalCluster( proc_id );
457  const size_type lastrow_dof_on_proc = this->dualImageSpace()->dof()->lastDofGlobalCluster( proc_id );
458  const size_type firstcol_dof_on_proc = this->domainSpace()->dof()->firstDofGlobalCluster( proc_id );
459  const size_type lastcol_dof_on_proc = this->domainSpace()->dof()->lastDofGlobalCluster( proc_id );
460 #endif
461 #if 0
462  graph_ptrtype sparsity_graph( new graph_type( nrow_dof_on_proc,
463  firstrow_dof_on_proc, lastrow_dof_on_proc,
464  firstcol_dof_on_proc, lastcol_dof_on_proc,
465  this->dualImageSpace()->mesh()->worldComm().subWorldComm() ) );
466 #else
467  graph_ptrtype sparsity_graph( new graph_type( this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
468 #endif
469  auto const* imagedof = this->dualImageSpace()->dof().get();
470  auto const* domaindof = this->domainSpace()->dof().get();
471  auto const* imagebasis = this->dualImageSpace()->basis().get();
472  auto const* domainbasis = this->domainSpace()->basis().get();
473 
474 
475  std::vector<bool> dof_done( nrow_dof_on_proc, false );
476  std::vector< std::list<std::pair<size_type,double> > > memory_valueInMatrix( nrow_dof_on_proc );
477 
478  // Local assembly: compute the Mloc matrix by evaluating
479  // the domain space basis function at the dual image space
480  // dof points (nodal basis) since we have only computation
481  // in the ref elements and the basis and dof points in ref
482  // element are the same, we compute Mloc outside the
483  // element loop.
484 
485  //typename matrix_node<value_type>::type Mloc(domain_basis_type::nLocalDof*domain_basis_type::nComponents1,1);
486  auto const& Mloc = domainbasis->evaluate( imagebasis->dual().points() );
487 
488  DVLOG(2) << "[interpolate] Same mesh but not same space\n";
489 
490  const bool image_related_to_domain = this->dualImageSpace()->mesh()->isSubMeshFrom( this->domainSpace()->mesh() );
491  const bool domain_related_to_image = this->domainSpace()->mesh()->isSubMeshFrom( this->dualImageSpace()->mesh() );
492 
493  auto itListRange = M_listRange.begin();
494  auto const enListRange = M_listRange.end();
495  for ( ; itListRange!=enListRange ; ++itListRange)
496  {
497  iterator_type it, en;
498  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
499  for ( ; it != en; ++ it )
500  {
501  auto idElem = detailsup::idElt( *it,idim_type() );
502  auto domain_eid = idElem;
503  if ( image_related_to_domain )
504  {
505  domain_eid = this->dualImageSpace()->mesh()->subMeshToMesh( idElem );
506  DVLOG(2) << "[image_related_to_domain] image element id: " << idElem << " domain element id : " << domain_eid << "\n";
507  }
508  if( domain_related_to_image )
509  {
510  domain_eid = this->domainSpace()->mesh()->meshToSubMesh( idElem );
511  DVLOG(2) << "[domain_related_to_image] image element id: " << idElem << " domain element id : " << domain_eid << "\n";
512  }
513 
514  if ( domain_eid == invalid_size_type_value )
515  continue;
516  // Global assembly
517  for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
518  {
519  for ( uint16_type comp = 0; comp < image_basis_type::nComponents; ++comp )
520  {
521  size_type i = boost::get<0>( imagedof->localToGlobal( *it, iloc, comp ) );
522 
523  if ( !dof_done[i] )
524  {
525 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
526  const auto ig1 = i;
527  const auto theproc = imagedof->worldComm().localRank();
528 #else // WITH MPI
529  const auto ig1 = imagedof->mapGlobalProcessToGlobalCluster()[i];
530  const auto theproc = imagedof->procOnGlobalCluster( ig1 );
531 #endif
532  auto& row = sparsity_graph->row( ig1 );
533  row.template get<0>() = theproc;
534  const size_type il1 = ig1 - imagedof->firstDofGlobalCluster( theproc );
535  row.template get<1>() = il1;
536  //row.template get<1>() = i;
537 
538  uint16_type ilocprime=imagedof->localDofInElement( *it, iloc, comp ) ;
539 
540  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
541  {
542  // get column
543  const size_type j = boost::get<0>( domaindof->localToGlobal( domain_eid, jloc, comp ) );
544  //up the pattern graph
545 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
546  row.template get<2>().insert( j );
547 #else // WITH MPI
548  row.template get<2>().insert( domaindof->mapGlobalProcessToGlobalCluster()[j] );
549 #endif
550  // get interpolated value
551  const value_type v = Mloc( domain_basis_type::nComponents1*jloc +
552  comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
553  comp,
554  ilocprime );
555  // save in matrux
556  memory_valueInMatrix[i].push_back( std::make_pair( j,v ) );
557  }
558 
559  dof_done[i]=true;
560  }
561  }
562  }
563  } // for ( ; it != en; ++ it )
564  } // for ( ; itListRange!=enListRange ; ++itListRange)
565 
566 
567  //-----------------------------------------
568  // compute graph
569  sparsity_graph->close();
570  //-----------------------------------------
571  // create matrix
572  this->matPtr() = this->backend()->newMatrix( this->domainSpace()->dofOnOff(),
573  this->dualImageSpace()->dofOn(),
574  sparsity_graph );
575  //-----------------------------------------
576 
577  // assemble matrix
578  for ( size_type idx_i=0 ; idx_i<nrow_dof_on_proc; ++idx_i )
579  {
580  for ( auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j )
581  {
582  this->matPtr()->set( idx_i,it_j->first,it_j->second );
583  }
584  }
585 }
586 
587 //-----------------------------------------------------------------------------------------------------------------//
588 //-----------------------------------------------------------------------------------------------------------------//
589 //-----------------------------------------------------------------------------------------------------------------//
590 
591 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
592 void
593 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMesh()
594 {
595  DVLOG(2) << "[interpolate] different meshes\n";
596  //std::cout << "OperatorInterpolation::updateNoRelationMesh start " << std::endl;
597 
598  const size_type proc_id = this->dualImageSpace()->mesh()->worldComm().localRank();
599  const size_type n1_dof_on_proc = this->dualImageSpace()->nLocalDof();
600  const size_type firstrow_dof_on_proc = this->dualImageSpace()->dof()->firstDof( proc_id );
601  const size_type lastrow_dof_on_proc = this->dualImageSpace()->dof()->lastDof( proc_id );
602  const size_type firstcol_dof_on_proc = this->domainSpace()->dof()->firstDof( proc_id );
603  const size_type lastcol_dof_on_proc = this->domainSpace()->dof()->lastDof( proc_id );
604 #if 0
605  graph_ptrtype sparsity_graph( new graph_type( n1_dof_on_proc,
606  firstrow_dof_on_proc, lastrow_dof_on_proc,
607  firstcol_dof_on_proc, lastcol_dof_on_proc,
608  this->dualImageSpace()->mesh()->worldComm().subWorldCommSeq() ) );
609 #else
610  graph_ptrtype sparsity_graph( new graph_type( this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
611 #endif
612 
613  auto const* imagedof = this->dualImageSpace()->dof().get();
614  auto const* domaindof = this->domainSpace()->dof().get();
615  auto const* domainbasis = this->domainSpace()->basis().get();
616 
617  //-----------------------------------------
618  //init the localization tool
619  auto locTool = this->domainSpace()->mesh()->tool_localization();
620  if ( this->interpolationType().onlyLocalizeOnBoundary() ) locTool->updateForUseBoundaryFaces();
621  else locTool->updateForUse();
622  // kdtree parameter
623  locTool->kdtree()->nbNearNeighbor(this->interpolationType().nbNearNeighborInKdTree());
624  bool notUseOptLocTest = domain_mesh_type::nDim!=domain_mesh_type::nRealDim;
625  if (notUseOptLocTest) locTool->kdtree()->nbNearNeighbor(domain_mesh_type::element_type::numPoints);
626 
627  //locTool->kdtree()->nbNearNeighbor(3);
628  //locTool->kdtree()->nbNearNeighbor(this->domainSpace()->mesh()->numElements());
629  //locTool->setExtrapolation(false);
630 
631  //-----------------------------------------
632  // usefull data
633  matrix_node_type ptsReal( image_mesh_type::nRealDim, 1 );
634  matrix_node_type ptsRef( domain_mesh_type::nDim , 1 );
635  typename domain_mesh_type::Localization::container_search_iterator_type itanal,itanal_end;
636  typename domain_mesh_type::Localization::container_output_iterator_type itL,itL_end;
637  matrix_node_type MlocEval( domain_basis_type::nLocalDof*domain_basis_type::nComponents1,1 );
638 
639  std::vector<bool> dof_done( this->dualImageSpace()->nLocalDof(), false );
640  std::vector< std::list<std::pair<size_type,double> > > memory_valueInMatrix( this->dualImageSpace()->nLocalDof() );
641 
642  //-----------------------------------------
643  size_type eltIdLocalised = 0;
644 
645  // for each element in range
646  auto itListRange = M_listRange.begin();
647  auto const enListRange = M_listRange.end();
648  for ( ; itListRange!=enListRange ; ++itListRange)
649  {
650  iterator_type it, en;
651  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
652  for ( ; it != en; ++ it )
653  {
654  for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
655  {
656  for ( uint16_type comp = 0; comp < image_basis_type::nComponents; ++comp )
657  {
658  const auto gdof = boost::get<0>(imagedof->localToGlobal( *it, iloc, comp ));
659  if (!dof_done[gdof])
660  {
661  //------------------------
662  // get the graph row
663 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
664  const auto ig1 = gdof;
665  const auto theproc = imagedof->worldComm().localRank();
666 #else // WITH MPI
667  const auto ig1 = imagedof->mapGlobalProcessToGlobalCluster()[gdof];
668  const auto theproc = imagedof->procOnGlobalCluster( ig1 );
669 #endif
670  auto& row = sparsity_graph->row(ig1);
671  row.template get<0>() = theproc;
672  row.template get<1>() = gdof;
673  //------------------------
674  // the dof point
675  ublas::column(ptsReal,0 ) = boost::get<0>(imagedof->dofPoint(gdof));
676  //------------------------
677  // localisation process
678  if (notUseOptLocTest) eltIdLocalised=invalid_size_type_value;
679  eltIdLocalised = locTool->run_analysis(ptsReal,eltIdLocalised,it->vertices()/*it->G()*/,mpl::int_<interpolation_type::value>()).template get<1>();
680  //------------------------
681  // for each localised points
682  itanal = locTool->result_analysis_begin();
683  itanal_end = locTool->result_analysis_end();
684  for ( ;itanal!=itanal_end;++itanal)
685  {
686  itL=itanal->second.begin();
687  ublas::column( ptsRef, 0 ) = boost::get<1>( *itL );
688 
689  MlocEval = domainbasis->evaluate( ptsRef );
690 
691  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
692  {
693  //get global dof
694  size_type j = boost::get<0>( domaindof->localToGlobal( itanal->first,jloc,comp ) );
695  value_type v = MlocEval( domain_basis_type::nComponents1*jloc
696  + comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof
697  + comp,
698  0 );
699 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
700  row.template get<2>().insert( j );
701 #else // WITH MPI
702  row.template get<2>().insert( domaindof->mapGlobalProcessToGlobalCluster()[j] );
703 #endif
704  memory_valueInMatrix[gdof].push_back( std::make_pair( j,v ) );
705  }
706  }
707  dof_done[gdof]=true;
708  } // if (!dof_done[gdof])
709  } // for ( uint16_type comp = 0; comp < image_basis_type::nComponents; ++comp )
710  } // for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
711  } // for( ; it != en; ++ it )
712  } // for ( ; itListRange!=enListRange ; ++itListRange)
713 
714 
715  //-----------------------------------------
716  // compute graph
717  sparsity_graph->close(); //sparsity_graph->printPython("mygraphpython.py");
718  //-----------------------------------------
719  // create matrix
720  this->matPtr() = this->backend()->newMatrix( this->domainSpace()->dofOnOff(),
721  this->dualImageSpace()->dofOn(),
722  sparsity_graph );
723  //-----------------------------------------
724  // assemble matrix
725  for (size_type idx_i=0 ; idx_i<this->dualImageSpace()->nLocalDof() ;++idx_i)
726  {
727  for (auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j)
728  {
729  this->matPtr()->set(idx_i,it_j->first,it_j->second);
730  }
731  }
732 
733 }
734 
735 //-----------------------------------------------------------------------------------------------------------------//
736 //-----------------------------------------------------------------------------------------------------------------//
737 //-----------------------------------------------------------------------------------------------------------------//
738 
739 #if defined(FEELPP_ENABLE_MPI_MODE) // WITH MPI
740 
741 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
742 void
743 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMeshMPI()
744 {
745  //std::cout << "OperatorInterpolation::updateNoRelationMeshMPI start " << std::endl;
746 
747  auto testCommActivities_image=this->dualImageSpace()->worldComm().hasMultiLocalActivity();
748 
749  if (testCommActivities_image.template get<0>())
750  {
751  //std::cout << "OperatorInterpolation::updateNoRelationMeshMPI hasMultiLocalActivity " << std::endl;
752  // save initial activities
753  std::vector<int> saveActivities_image = this->dualImageSpace()->worldComm().activityOnWorld();
754  // iterate on each local activity
755  const auto colorWhichIsActive = testCommActivities_image.template get<1>();
756  auto it_color=colorWhichIsActive.begin();
757  auto const en_color=colorWhichIsActive.end();
758  for ( ;it_color!=en_color;++it_color )
759  {
760  this->dualImageSpace()->worldComm().applyActivityOnlyOn( *it_color );
761  this->dualImageSpace()->mapOn().worldComm().applyActivityOnlyOn( *it_color );
762  this->dualImageSpace()->mapOnOff().worldComm().applyActivityOnlyOn( *it_color );
763  this->updateNoRelationMeshMPI_run(false);
764  }
765  // revert initial activities
766  this->dualImageSpace()->worldComm().setIsActive(saveActivities_image);
767  this->dualImageSpace()->mapOn().worldComm().setIsActive(saveActivities_image);
768  this->dualImageSpace()->mapOnOff().worldComm().setIsActive(saveActivities_image);
769  }
770  else
771  {
772  //std::cout << "OperatorInterpolation::updateNoRelationMeshMPI has One LocalActivity " << std::endl;
773  if ( !this->dualImageSpace()->worldComm().isActive() && !this->domainSpace()->worldComm().isActive() )
774  {
775  this->matPtr() = this->backend()->newZeroMatrix( this->domainSpace()->dofOnOff(),
776  this->dualImageSpace()->dofOn() );
777  }
778  else
779  {
780  this->updateNoRelationMeshMPI_run(true);
781  }
782  }
783 
784 }
785 
786 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
787 void
788 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMeshMPI_run(bool buildNonZeroMatrix)
789 {
790 
791  //std::cout << "OperatorInterpolation::updateNoRelationMeshMPI_run start " << std::endl;
792 
793  //-----------------------------------------------------------------------------------------
794  // PreProcess : datamap properties and graph
795  //-----------------------------------------------------------------------------------------
796 
797  const size_type proc_id = this->dualImageSpace()->worldComm().localRank();
798  const size_type proc_id_row = this->dualImageSpace()->worldComm().localRank();
799  const size_type proc_id_col = this->domainSpace()->worldComm().localRank();
800  const size_type nProc = this->dualImageSpace()->mesh()->worldComm().size();
801  const size_type nProc_row = this->dualImageSpace()->mesh()->worldComm().localSize();
802  const size_type nProc_col = this->domainSpace()->mesh()->worldComm().localSize();
803  const size_type nrow_dof_on_proc = this->dualImageSpace()->nLocalDof();
804  const size_type firstrow_dof_on_proc = this->dualImageSpace()->dof()->firstDofGlobalCluster( proc_id_row );
805  const size_type lastrow_dof_on_proc = this->dualImageSpace()->dof()->lastDofGlobalCluster( proc_id_row );
806  const size_type firstcol_dof_on_proc = this->domainSpace()->dof()->firstDofGlobalCluster( proc_id_col );
807  const size_type lastcol_dof_on_proc = this->domainSpace()->dof()->lastDofGlobalCluster( proc_id_col );
808 
809 
810 
811  std::vector<size_type> first_col_entry( this->domainSpace()->worldComm().globalSize() );
812  std::vector<size_type> last_col_entry( this->domainSpace()->worldComm().globalSize() );
813  mpi::all_gather( this->domainSpace()->worldComm().globalComm(),
814  firstcol_dof_on_proc,//this->firstColEntryOnProc(),
815  first_col_entry );
816  mpi::all_gather( this->domainSpace()->worldComm().globalComm(),
817  lastcol_dof_on_proc,//this->lastColEntryOnProc(),
818  last_col_entry );
819  size_type thefirstCol = *std::min_element( first_col_entry.begin(),first_col_entry.end() );
820  size_type thelastCol = *std::max_element( last_col_entry.begin(),last_col_entry.end() );
821 
822 #if 0
823  graph_ptrtype sparsity_graph( new graph_type( nrow_dof_on_proc,
824  firstrow_dof_on_proc, lastrow_dof_on_proc,
825  thefirstCol,thelastCol,
826  this->dualImageSpace()->worldComm() ) );
827 #else
828  graph_ptrtype sparsity_graph( new graph_type(this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
829 #endif
830 
831  size_type new_nLocalDofWithoutGhost=this->domainSpace()->nDof()/nProc_row;
832  size_type new_nLocalDofWithoutGhost_tempp=this->domainSpace()->nDof()/nProc_row;
833  size_type new_nLocalDofWithoutGhostMiss=this->domainSpace()->nDof()%nProc_row;
834  if (this->dualImageSpace()->worldComm().globalSize()==this->domainSpace()->worldComm().globalSize() )
835  {
836  new_nLocalDofWithoutGhost = this->domainSpace()->mapOnOff().nLocalDofWithoutGhost();
837  }
838  else
839  {
840  if (this->dualImageSpace()->worldComm().globalRank()==this->dualImageSpace()->worldComm().masterRank()) new_nLocalDofWithoutGhost+=new_nLocalDofWithoutGhostMiss;
841  }
842 
843  size_type new_firstdofcol=0,new_lastdofcol=new_nLocalDofWithoutGhost-1;
844  bool findMyProc=false;
845  int currentProc=0;
846  while(!findMyProc)
847  {
848  if (currentProc==this->dualImageSpace()->worldComm().globalRank())
849  {
850  findMyProc=true;
851  }
852  else if (currentProc==this->dualImageSpace()->worldComm().masterRank())
853  {
854  new_firstdofcol+=new_nLocalDofWithoutGhost_tempp+new_nLocalDofWithoutGhostMiss;
855  new_lastdofcol+=new_nLocalDofWithoutGhost_tempp+new_nLocalDofWithoutGhostMiss;
856  }
857  else
858  {
859  new_firstdofcol+=new_nLocalDofWithoutGhost_tempp;
860  new_lastdofcol+=new_nLocalDofWithoutGhost_tempp;
861  }
862  ++currentProc;
863  }
864 
865 
866  //-----------------------------------------------------------------------------------------
867  // Start localization process
868  //-----------------------------------------------------------------------------------------
869 
870  //memory containers
871  std::vector< std::list<boost::tuple<int,size_type,double> > > memory_valueInMatrix( this->dualImageSpace()->nLocalDof() );
872  for (size_type k=0;k<memory_valueInMatrix.size();++k) memory_valueInMatrix[k].clear();
873 
874  std::vector<std::map<size_type,size_type> > memory_col_globalProcessToGlobalCluster(nProc_col);
875  std::vector<std::set<size_type> > dof_searchWithProc(this->dualImageSpace()->nLocalDof());
876 
877  extrapolation_memory_type dof_extrapolationData(this->dualImageSpace()->nLocalDof());
878 
879  //init the localization tool
880  auto locTool = this->domainSpace()->mesh()->tool_localization();
881  bool doExtrapolationAtStart = locTool->doExtrapolation();
882  // kdtree parameter
883  locTool->kdtree()->nbNearNeighbor(this->interpolationType().nbNearNeighborInKdTree());
884  // points in kdtree
885  if ( this->interpolationType().onlyLocalizeOnBoundary() ) locTool->updateForUseBoundaryFaces();
886  else locTool->updateForUse();
887  // no extrapolation in first
888  if ( doExtrapolationAtStart ) locTool->setExtrapolation(false);
889 
890 
891  uint16_type nMPIsearch=15;//5;
892  if( InterpType::value==1) nMPIsearch=this->domainSpace()->mesh()->worldComm().localSize();
893  else if (this->domainSpace()->mesh()->worldComm().localSize()<nMPIsearch) nMPIsearch=this->domainSpace()->mesh()->worldComm().localSize();
894  // only one int this case
895  if (!this->interpolationType().searchWithCommunication()) nMPIsearch=1;
896  uint16_type counterMPIsearch=1;
897  bool FinishMPIsearch=false;
898  //if (this->interpolationType().searchWithCommunication()) FinishMPIsearch=true;// not run algo1 !!!!
899 
900  boost::mpi::timer mytimer;
901  this->worldCommFusion().globalComm().barrier();
902  if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
903  std::cout << " start while " << std::endl;
904 
905  size_type nbLocalisationFail=1;
906  while(!FinishMPIsearch)
907  {
908  mytimer.restart();
909  auto pointDistribution = this->updateNoRelationMeshMPI_pointDistribution(memory_valueInMatrix,dof_searchWithProc);
910  auto memmapGdof = pointDistribution.template get<0>();
911  auto memmapComp = pointDistribution.template get<1>();
912  auto pointsSearched = pointDistribution.template get<2>();
913  auto memmapVertices = pointDistribution.template get<3>();
914  //std::cout << "proc " << this->worldCommFusion().globalRank() << " pointsSearched.size() " << pointsSearched.size() << std::endl;
915  double t1 = mytimer.elapsed();
916  if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
917  std::cout << "finish-step1 in " << (boost::format("%1%") % t1).str() << std::endl;
918  //this->worldCommFusion().globalComm().barrier();
919  mytimer.restart();
920 
921  auto memory_localisationFail = this->updateNoRelationMeshMPI_upWithMyWorld( memmapGdof, // input
922  memmapComp, // input
923  pointsSearched, // input
924  memmapVertices, // input
925  sparsity_graph, // output
926  memory_valueInMatrix, // output
927  memory_col_globalProcessToGlobalCluster, // output
928  dof_searchWithProc, // output,
929  false,// extrapolation_mode
930  dof_extrapolationData // empty
931  );
932  //std::cout << "proc " << this->worldCommFusion().globalRank() << " memory_localisationFail.size() " << memory_localisationFail.size() << std::endl;
933  double t2 = mytimer.elapsed();
934  if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
935  std::cout << "finish-step2 in " << (boost::format("%1%") % t2).str() << std::endl;
936  //this->worldCommFusion().globalComm().barrier();
937  mytimer.restart();
938 
939 
940  if (this->interpolationType().searchWithCommunication())
941  {
942  auto memory_localisationFail2 = this->updateNoRelationMeshMPI_upWithOtherWorld( memmapGdof, // input
943  memmapComp, // input
944  pointsSearched, // input
945  memmapVertices, // input
946  sparsity_graph, // output
947  memory_valueInMatrix, // output
948  memory_col_globalProcessToGlobalCluster, // output
949  dof_searchWithProc, // output
950  false,// extrapolation_mode
951  dof_extrapolationData // empty
952  );
953 
954  const size_type nbLocalisationFail_loc = memory_localisationFail.size()+memory_localisationFail2.size();
955  mpi::all_reduce( this->worldCommFusion().globalComm(),
956  nbLocalisationFail_loc,
957  nbLocalisationFail,
958  std::plus<double>() );
959  }
960  else nbLocalisationFail = memory_localisationFail.size();
961  //nbLocalisationFail = memory_localisationFail.size()+memory_localisationFail2.size();
962 
963  double t3 = mytimer.elapsed();
964  if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
965  std::cout << "finish-step3 in " << (boost::format("%1%") % t3).str() << std::endl;
966  mytimer.restart();
967 
968  if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
969  std::cout << " it " << counterMPIsearch << " nbLocalisationFail " << nbLocalisationFail << std::endl;
970 
971  //std::cout << "proc " << this->worldCommFusion().globalRank()
972  // << " et " <<nbLocalisationFail << std::endl;
973  if (counterMPIsearch<nMPIsearch && nbLocalisationFail>0) ++counterMPIsearch;
974  else FinishMPIsearch=true;
975  }
976 
977  //std::cout << "\n FINISH SEARCH!!!!!!!!! " << std::endl;
978  if ( doExtrapolationAtStart ) locTool->setExtrapolation(true);
979 
980  if ( doExtrapolationAtStart && nbLocalisationFail>0 )
981  {
982  //std::cout << " Start Extrapolation" << std::endl;
983  std::vector<std::set<size_type> > dof_searchWithProcExtrap(this->dualImageSpace()->nLocalDof());
984  //locTool->setExtrapolation(true);
985  uint16_type nMPIsearchExtrap=5;
986  if (this->domainSpace()->mesh()->worldComm().localSize()<5) nMPIsearchExtrap=this->domainSpace()->mesh()->worldComm().localSize();
987  // only one int this case
988  if (!this->interpolationType().searchWithCommunication()) nMPIsearchExtrap=1;
989  uint16_type counterMPIsearchExtrap=1;
990  bool FinishMPIsearchExtrap=false;
991  // localisation process
992  while(!FinishMPIsearchExtrap)
993  {
994  auto pointDistribution = this->updateNoRelationMeshMPI_pointDistribution(memory_valueInMatrix,dof_searchWithProcExtrap);
995  auto memmapGdof = pointDistribution.template get<0>();
996  auto memmapComp = pointDistribution.template get<1>();
997  auto pointsSearched = pointDistribution.template get<2>();
998  auto memmapVertices = pointDistribution.template get<3>();
999  //std::cout << "proc " << this->worldCommFusion().globalRank() << " pointsSearched.size() " << pointsSearched.size() << std::endl;
1000  auto memory_localisationFail = this->updateNoRelationMeshMPI_upWithMyWorld( memmapGdof, // input
1001  memmapComp, // input
1002  pointsSearched, // input
1003  memmapVertices, // input
1004  sparsity_graph, // output
1005  memory_valueInMatrix, // output
1006  memory_col_globalProcessToGlobalCluster, // output
1007  dof_searchWithProcExtrap, // output,
1008  true,// extrapolation_mode
1009  dof_extrapolationData // output
1010  );
1011  //std::cout << "proc " << this->worldCommFusion().globalRank() << " memory_localisationFail.size() " << memory_localisationFail.size() << std::endl;
1012  if (this->interpolationType().searchWithCommunication())
1013  {
1014  auto memory_localisationFail2 = this->updateNoRelationMeshMPI_upWithOtherWorld( memmapGdof, // input
1015  memmapComp, // input
1016  pointsSearched, // input
1017  memmapVertices, // input
1018  sparsity_graph, // output
1019  memory_valueInMatrix, // output
1020  memory_col_globalProcessToGlobalCluster, // output
1021  dof_searchWithProcExtrap, // output
1022  true, // extrapolation_mode
1023  dof_extrapolationData // output
1024  );
1025  }
1026  //std::cout << "proc " << this->worldCommFusion().globalRank() << " memory_localisationFail2.size() " << memory_localisationFail.size() << std::endl;
1027  if (counterMPIsearchExtrap<nMPIsearchExtrap) ++counterMPIsearchExtrap;
1028  else FinishMPIsearchExtrap=true;
1029  } // while(!FinishMPIsearchExtrap)
1030 
1031  auto const* imagedof = this->dualImageSpace()->dof().get();
1032  auto const* domaindof = this->domainSpace()->dof().get();
1033  auto const* domainbasis = this->domainSpace()->basis().get();
1034  matrix_node_type ptsRef( domain_mesh_type::nRealDim , 1 );
1035  matrix_node_type MlocEval(domain_basis_type::nLocalDof*domain_basis_type::nComponents1,1);
1036 
1037  // analysis result
1038  auto it_extrap = dof_extrapolationData.begin();
1039  auto const en_extrap = dof_extrapolationData.end();
1040  for (size_type cpt_gdof=0 ; it_extrap!=en_extrap ; ++it_extrap,++cpt_gdof )
1041  {
1042  auto const gdof = cpt_gdof;
1043  auto const pointExtrap = imagedof->dofPoint(gdof).template get<0>();
1044 
1045  // search nearer
1046  int procExtrapoled = 0;
1047  double distMin= INT_MAX;
1048  auto it_proc = it_extrap->begin();
1049  auto const en_proc = it_extrap->end();
1050  for ( ; it_proc!=en_proc; ++it_proc)
1051  {
1052  auto const procCurrent = it_proc->template get<0>();
1053  auto const bary = it_proc->template get<1>();
1054  // COMPUTE DISTANCE
1055  double normDist = 0;
1056  for (int q=0;q<image_mesh_type::nRealDim;++q) normDist+=std::pow(bary(q)-pointExtrap(q),2);
1057  //std::cout << std::sqrt(normDist) << " "<< bary.size() << " " << pointExtrap.size() << std::endl;
1058  if (std::sqrt(normDist)<distMin) { distMin=std::sqrt(normDist);procExtrapoled=procCurrent;}
1059  }
1060  it_proc = it_extrap->begin();
1061  for ( ; it_proc!=en_proc; ++it_proc)
1062  {
1063  if ( it_proc->template get<0>()==procExtrapoled)
1064  {
1065  // get the graph row
1066  auto const ig1 = imagedof->mapGlobalProcessToGlobalCluster()[gdof];
1067  auto const theproc = imagedof->procOnGlobalCluster(ig1);
1068  auto& row = sparsity_graph->row(ig1);
1069  row.template get<0>() = theproc;
1070  row.template get<1>() = ig1 - imagedof->firstDofGlobalCluster(theproc);
1071  // get ref point
1072  ublas::column( ptsRef, 0 ) = it_proc->template get<2>();
1073  // evaluate basis functions for this point
1074  MlocEval = domainbasis->evaluate( ptsRef );
1075  //auto it_jdof = it_proc->template get<3>().begin();
1076  //auto const en_jdof = it_proc->template get<3>().end();
1077  auto const comp = it_proc->template get<4>();
1078  auto const vec_jloc = it_proc->template get<3>();
1079  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1080  //for (uint16_type jloc=0 ; it_jdof != en_jdof ; ++it_jdof,++jloc)
1081  {
1082  //const size_type j_gdof=it_jdof->first;
1083  //const size_type j_gdofGlobalCluster=it_jdof->second;
1084  const size_type j_gdof=vec_jloc[jloc].first;
1085  const size_type j_gdofGlobalCluster=vec_jloc[jloc].second;
1086  row.template get<2>().insert(j_gdofGlobalCluster);//domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1087  // get value
1088  value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1089  comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1090  comp, 0 );
1091  // save value
1092  memory_valueInMatrix[gdof].push_back(boost::make_tuple(procExtrapoled,j_gdof,v));
1093  memory_col_globalProcessToGlobalCluster[procExtrapoled][j_gdof]=j_gdofGlobalCluster;//domaindof->mapGlobalProcessToGlobalCluster()[j_gdof];
1094  }
1095  } // if ( it_proc->template get<0>()==procExtrapoled)
1096  } // for ( ; it_proc!=en_proc; ++it_proc)
1097 
1098  } // for (size_type cpt_gdof=0 ; it_extrap!=en_extrap ; ++it_extrap,++cpt_gdof )
1099  } // if ( doExtrapolationAtStart )
1100 
1101 
1102  //-----------------------------------------------------------------------------------------
1103  this->worldCommFusion().barrier();
1104  //std::cout << "Op---1----- " << std::endl;
1105  //-----------------------------------------------------------------------------------------
1106  // compute graph
1107  sparsity_graph->close();//sparsity_graph->printPython("mygraphpythonMPI.py");
1108  //-----------------------------------------------------------------------------------------
1109  //std::cout << "Op---2----- " << std::endl;
1110  this->worldCommFusion().barrier();
1111  //-----------------------------------------------------------------------------------------
1112  size_type mapCol_nLocalDof = 0;
1113  for (int p=0;p<nProc_col;++p)
1114  {
1115  mapCol_nLocalDof += memory_col_globalProcessToGlobalCluster[p].size();
1116  }
1117  //std::cout << "mapCol_nLocalDof " << mapCol_nLocalDof << std::endl;
1118  std::vector<size_type> mapCol_globalProcessToGlobalCluster(mapCol_nLocalDof);
1119  std::vector<size_type> new_mapGlobalClusterToGlobalProcess(new_nLocalDofWithoutGhost);
1120  std::vector<std::map<size_type,size_type> > mapCol_LocalSpaceDofToLocalInterpDof(nProc_col);
1121  size_type currentLocalDof=0;
1122  for (int p = 0 ; p<nProc_col;++p)
1123  {
1124  auto it_map = memory_col_globalProcessToGlobalCluster[p].begin();
1125  auto en_map = memory_col_globalProcessToGlobalCluster[p].end();
1126  for ( ; it_map!=en_map ; ++it_map)
1127  {
1128  mapCol_globalProcessToGlobalCluster[currentLocalDof]=it_map->second;//memory_col_globalProcessToGlobalCluster[proc_id][i];
1129  mapCol_LocalSpaceDofToLocalInterpDof[p][it_map->first]=currentLocalDof;
1130 
1131  if ( new_firstdofcol<=it_map->second && new_lastdofcol>=it_map->second)
1132  new_mapGlobalClusterToGlobalProcess[it_map->second-new_firstdofcol]=currentLocalDof;
1133 
1134  ++currentLocalDof;
1135  }
1136  }
1137 
1138  //-----------------------------------------
1139  //std::cout << "Op---3----- " << this->worldCommFusion().godRank() << std::endl;
1140  this->worldCommFusion().barrier();
1141  //-----------------------------------------
1142  // build data map for the columns
1143  //this->domainSpace()->mapOnOff().showMeMapGlobalProcessToGlobalCluster();
1144  //this->dualImageSpace()->worldComm().showMe();
1145  boost::shared_ptr<DataMap> mapColInterp( new DataMap(this->dualImageSpace()->worldComm()));// this->domainSpace()->mapOnOff().worldComm());
1146  mapColInterp->setNDof(this->domainSpace()->mapOnOff().nDof());
1147 
1148  mapColInterp->setNLocalDofWithoutGhost( proc_id, new_nLocalDofWithoutGhost );// this->domainSpace()->mapOnOff().nLocalDofWithoutGhost() );
1149  mapColInterp->setNLocalDofWithGhost( proc_id, mapCol_nLocalDof/*this->domainSpace()->mapOnOff().nLocalDofWithGhost()*/ );
1150  mapColInterp->setFirstDof( proc_id, this->domainSpace()->mapOnOff().firstDof() );
1151  mapColInterp->setLastDof( proc_id, this->domainSpace()->mapOnOff().lastDof() );
1152  mapColInterp->setFirstDofGlobalCluster( proc_id, new_firstdofcol );
1153  mapColInterp->setLastDofGlobalCluster( proc_id, new_lastdofcol );
1154  mapColInterp->setMapGlobalProcessToGlobalCluster(mapCol_globalProcessToGlobalCluster);
1155  mapColInterp->setMapGlobalClusterToGlobalProcess(new_mapGlobalClusterToGlobalProcess);
1156  //if ( this->dualImageSpace()->worldComm().isActive() ) mapColInterp.showMeMapGlobalProcessToGlobalCluster();
1157 
1158  //-----------------------------------------
1159  //std::cout << "Op---4----- " << this->worldCommFusion().godRank() << " isA " << this->dualImageSpace()->worldComm().isActive() << std::endl;
1160  this->worldCommFusion().barrier();
1161  //-----------------------------------------
1162  // create matrix for active process
1163  if ( this->dualImageSpace()->worldComm().isActive() )
1164  {
1165  this->matPtr() = this->backend()->newMatrix( mapColInterp,//this->domainSpace()->mapOnOff(),
1166  this->dualImageSpace()->dofOn(),
1167  sparsity_graph );
1168  }
1169  //-----------------------------------------
1170  //std::cout << "Op---5----- " << this->worldCommFusion().godRank() << std::endl;
1171  this->worldCommFusion().barrier();
1172  //-----------------------------------------
1173  // create null matrix for inactive process
1174  if ( !this->dualImageSpace()->worldComm().isActive() && buildNonZeroMatrix )
1175  {
1176  this->matPtr() = this->backend()->newZeroMatrix( mapColInterp,//this->domainSpace()->mapOnOff(),
1177  this->dualImageSpace()->dofOn() );
1178  }
1179  //-----------------------------------------
1180  //std::cout << "Op---6----- " << this->worldCommFusion().godRank() << std::endl;
1181  this->worldCommFusion().barrier();
1182  //-----------------------------------------
1183  // assemble matrix
1184  if ( this->dualImageSpace()->worldComm().isActive() )
1185  {
1186  for (size_type idx_i=0 ; idx_i<nrow_dof_on_proc;++idx_i)
1187  {
1188  for (auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j)
1189  {
1190  if (memory_valueInMatrix[idx_i].size()>0)
1191  {
1192  this->matPtr()->set(idx_i,mapCol_LocalSpaceDofToLocalInterpDof[it_j->template get<0>()][it_j->template get<1>()],it_j->template get<2>());
1193  }
1194  }
1195  }
1196  }
1197  //-----------------------------------------
1198  //std::cout << "Op---7----- " << std::endl;
1199  this->worldCommFusion().barrier();
1200  //-----------------------------------------
1201 
1202 
1203 }
1204 
1205 //-----------------------------------------------------------------------------------------------------------------//
1206 
1207 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
1208 std::list<boost::tuple<size_type,uint16_type> >
1209 OperatorInterpolation<DomainSpaceType,
1210  ImageSpaceType,
1211  IteratorRange,
1212  InterpType>::updateNoRelationMeshMPI_upWithMyWorld(const std::vector< std::vector<size_type> > & memmapGdof,
1213  const std::vector< std::vector<uint16_type> > & memmapComp,
1214  const std::vector<std::vector<typename image_mesh_type::node_type> > & pointsSearched,
1215  const std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > & memmap_vertices,
1216  graph_ptrtype & sparsity_graph,
1217  std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
1218  std::vector<std::map<size_type,size_type> > & memory_col_globalProcessToGlobalCluster,
1219  std::vector<std::set<size_type> > & dof_searchWithProc,
1220  bool extrapolation_mode,
1221  extrapolation_memory_type & dof_extrapolationData )
1222 {
1223  std::list<boost::tuple<size_type,uint16_type> > memory_localisationFail;// gdof,comp
1224 
1225  const auto proc_id = this->domainSpace()->mesh()->worldComm().localRank();
1226 
1227  auto const* imagedof = this->dualImageSpace()->dof().get();
1228  auto const* domaindof = this->domainSpace()->dof().get();
1229  auto const* domainbasis = this->domainSpace()->basis().get();
1230 
1231  auto locTool = this->domainSpace()->mesh()->tool_localization();
1232  bool notUseOptLocTest = domain_mesh_type::nDim!=domain_mesh_type::nRealDim;
1233  if (notUseOptLocTest) locTool->kdtree()->nbNearNeighbor(domain_mesh_type::element_type::numPoints);
1234 
1235  matrix_node_type ptsReal( image_mesh_type::nRealDim, 1 );
1236  matrix_node_type ptsRef( domain_mesh_type::nDim , 1 );
1237  matrix_node_type MlocEval(domain_basis_type::nLocalDof*domain_basis_type::nComponents1,1);
1238  matrix_node_type verticesOfEltSearched;
1239 
1240  //size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1241  size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1242  auto const& eltRandom = this->domainSpace()->mesh()->element(eltIdLocalised);
1243 
1244  for ( size_type k=0 ; k<memmapGdof[proc_id].size() ; ++k)
1245  {
1246  //----------------------------------------------------------------
1247  if (this->interpolationType().componentsAreSamePoint())
1248  //------------------------------------------------------------
1249  {
1250  // the searched point
1251  ublas::column(ptsReal,0 ) = pointsSearched[proc_id][k];
1252  // vertice with conforme case
1253  if (InterpType::value==1)
1254  {
1255  const uint16_type sizeVertices = memmap_vertices[proc_id][k].size();
1256  verticesOfEltSearched.resize( image_mesh_type::nRealDim,sizeVertices);
1257  for ( uint16_type v=0;v<sizeVertices;++v)
1258  ublas::column(verticesOfEltSearched,v)=memmap_vertices[proc_id][k][v];
1259  }
1260  else // random
1261  verticesOfEltSearched = eltRandom.vertices();
1262 
1263  // localisation process
1264  if (notUseOptLocTest) eltIdLocalised=invalid_size_type_value;
1265  auto resLocalisation = locTool->run_analysis(ptsReal,eltIdLocalised,verticesOfEltSearched,
1266  mpl::int_<interpolation_type::value>());
1267  if (!resLocalisation.template get<0>()[0]) // not find
1268  {
1269  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1270  {
1271  const auto gdof = memmapGdof[proc_id][k+comp];
1272  memory_localisationFail.push_back(boost::make_tuple(gdof,comp) );
1273  dof_searchWithProc[gdof].insert(proc_id);
1274  }
1275  }
1276  else // point found
1277  {
1278  eltIdLocalised = resLocalisation.template get<1>();
1279 
1280  if (extrapolation_mode)
1281  {
1282  auto const& eltExtrapoled = this->domainSpace()->mesh()->element(eltIdLocalised);
1283  auto const& verticesExtrapoled = eltExtrapoled.vertices();
1284  typename image_mesh_type::node_type bary(verticesExtrapoled.size1());
1285  for (int qi=0;qi<verticesExtrapoled.size1();++qi) bary(qi)=0;//important
1286  for (int qj=0;qj<verticesExtrapoled.size2();++qj)
1287  {
1288  bary(0) += ublas::column(verticesExtrapoled,qj)(0);
1289  if (verticesExtrapoled.size1()>1) bary(1) += ublas::column(verticesExtrapoled,qj)(1);
1290  if (verticesExtrapoled.size1()>2) bary(2) += ublas::column(verticesExtrapoled,qj)(2);
1291  }
1292  bary(0) /= verticesExtrapoled.size2();
1293  if (verticesExtrapoled.size1()>1) bary(1) /= verticesExtrapoled.size2();
1294  if (verticesExtrapoled.size1()>2) bary(2) /= verticesExtrapoled.size2();
1295  typename image_mesh_type::node_type theRefPtExtrap = locTool->result_analysis().begin()->second.begin()->template get<1>();
1296  std::vector<std::pair<size_type,size_type> > j_gdofs(domain_basis_type::nLocalDof);
1297  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1298  {
1299  const auto gdof = memmapGdof[proc_id][k+comp];
1300  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1301  {
1302  const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( eltIdLocalised,jloc,comp ));
1303  j_gdofs[jloc]=std::make_pair(j_gdof,domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1304  }
1305  dof_extrapolationData[gdof].push_back(boost::make_tuple(proc_id,bary,theRefPtExtrap,j_gdofs,comp));
1306  }
1307  }
1308  else // no extrapolation
1309  {
1310  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1311  {
1312  const auto gdof = memmapGdof[proc_id][k+comp];
1313  // get the graph row
1314  auto const ig1 = imagedof->mapGlobalProcessToGlobalCluster()[gdof];
1315  auto const theproc = imagedof->procOnGlobalCluster(ig1);
1316  auto& row = sparsity_graph->row(ig1);
1317  row.template get<0>() = theproc;
1318  row.template get<1>() = ig1 - imagedof->firstDofGlobalCluster(theproc);
1319  // for each localised points
1320  auto itanal = locTool->result_analysis().begin();// result_analysis_begin();
1321  auto const itanal_end = locTool->result_analysis().end();//result_analysis_end();
1322  for ( ;itanal!=itanal_end;++itanal)
1323  {
1324  const auto itL=itanal->second.begin();
1325  ublas::column( ptsRef, 0 ) = boost::get<1>(*itL);
1326  // evaluate basis functions for this point
1327  MlocEval = domainbasis->evaluate( ptsRef );
1328  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1329  {
1330  //get global dof
1331  const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( itanal->first,jloc,comp ));
1332  // up graph
1333  row.template get<2>().insert(domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1334  // get value
1335  const value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1336  comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1337  comp, 0 );
1338  // save value
1339  memory_valueInMatrix[gdof].push_back(boost::make_tuple(proc_id,j_gdof,v));
1340  memory_col_globalProcessToGlobalCluster[proc_id][j_gdof]=domaindof->mapGlobalProcessToGlobalCluster()[j_gdof];
1341  }
1342  }
1343  // dof ok : not anymore localise
1344  dof_searchWithProc[gdof].insert(proc_id);
1345  } // comp
1346  } // no extrapolation
1347  } // else // point found
1348  // change k in for
1349  k=k+image_basis_type::nComponents-1;
1350  } // optimization : this->interpolationType().componentsAreSamePoint()
1351  //----------------------------------------------------
1352  else // component optimization
1353  //------------------------------------------------
1354  {
1355  const auto gdof = memmapGdof[proc_id][k];
1356  const auto comp = memmapComp[proc_id][k];
1357  ublas::column(ptsReal,0 ) = imagedof->dofPoint(gdof).template get<0>();
1358  // vertice with conforme case
1359  if (InterpType::value==1)
1360  {
1361  const uint16_type sizeVertices = memmap_vertices[proc_id][k].size();
1362  verticesOfEltSearched.resize( image_mesh_type::nRealDim,sizeVertices);
1363  for ( uint16_type v=0;v<sizeVertices;++v)
1364  ublas::column(verticesOfEltSearched,v)=memmap_vertices[proc_id][k][v];
1365  }
1366  else // random
1367  verticesOfEltSearched = eltRandom.vertices();
1368 
1369  // localisation process
1370  if (notUseOptLocTest) eltIdLocalised=invalid_size_type_value;
1371  auto resLocalisation = locTool->run_analysis(ptsReal, eltIdLocalised, verticesOfEltSearched, mpl::int_<interpolation_type::value>());
1372  if (!resLocalisation.template get<0>()[0]) // not find
1373  {
1374  memory_localisationFail.push_back(boost::make_tuple(gdof,comp) );
1375  dof_searchWithProc[gdof].insert(proc_id);
1376  }
1377  else // point found
1378  {
1379  eltIdLocalised = resLocalisation.template get<1>();
1380  if (extrapolation_mode)
1381  {
1382  auto const& eltExtrapoled = this->domainSpace()->mesh()->element(eltIdLocalised);
1383  auto const& verticesExtrapoled = eltExtrapoled.vertices();
1384  typename image_mesh_type::node_type bary(verticesExtrapoled.size1());
1385  for (int qi=0;qi<verticesExtrapoled.size1();++qi) bary(qi)=0;//important
1386  for (int qj=0;qj<verticesExtrapoled.size2();++qj)
1387  {
1388  bary(0) += ublas::column(verticesExtrapoled,qj)(0);
1389  if (verticesExtrapoled.size1()>1) bary(1) += ublas::column(verticesExtrapoled,qj)(1);
1390  if (verticesExtrapoled.size1()>2) bary(2) += ublas::column(verticesExtrapoled,qj)(2);
1391  }
1392  bary(0) /= verticesExtrapoled.size2();
1393  if (verticesExtrapoled.size1()>1) bary(1) /= verticesExtrapoled.size2();
1394  if (verticesExtrapoled.size1()>2) bary(2) /= verticesExtrapoled.size2();
1395  typename image_mesh_type::node_type theRefPtExtrap = locTool->result_analysis().begin()->second.begin()->template get<1>();
1396  std::vector<std::pair<size_type,size_type> > j_gdofs(domain_basis_type::nLocalDof);
1397  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1398  {
1399  const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( eltIdLocalised,jloc,comp ));
1400  j_gdofs[jloc]=std::make_pair(j_gdof,domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1401  }
1402  dof_extrapolationData[gdof].push_back(boost::make_tuple(proc_id,bary,theRefPtExtrap,j_gdofs,comp));
1403  }
1404  else
1405  {
1406  // get the graph row
1407  auto const ig1 = imagedof->mapGlobalProcessToGlobalCluster()[gdof];
1408  auto const theproc = imagedof->procOnGlobalCluster(ig1);
1409  auto& row = sparsity_graph->row(ig1);
1410  row.template get<0>() = theproc;
1411  row.template get<1>() = ig1 - imagedof->firstDofGlobalCluster(theproc);
1412  // for each localised points
1413  auto itanal = locTool->result_analysis().begin();// result_analysis_begin();
1414  auto const itanal_end = locTool->result_analysis().end();//result_analysis_end();
1415  for ( ;itanal!=itanal_end;++itanal)
1416  {
1417  const auto itL=itanal->second.begin();
1418  ublas::column( ptsRef, 0 ) = boost::get<1>(*itL);
1419  // evaluate basis functions for this point
1420  MlocEval = domainbasis->evaluate( ptsRef );
1421  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1422  {
1423  //get global dof
1424  const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( itanal->first,jloc,comp ));
1425  // up graph
1426  row.template get<2>().insert(domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1427  // get value
1428  const value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1429  comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1430  comp, 0 );
1431  // save value
1432  memory_valueInMatrix[gdof].push_back(boost::make_tuple(proc_id,j_gdof,v));
1433  memory_col_globalProcessToGlobalCluster[proc_id][j_gdof]=domaindof->mapGlobalProcessToGlobalCluster()[j_gdof];
1434  }
1435  }
1436  // dof ok : not anymore localise
1437  dof_searchWithProc[gdof].insert(proc_id);
1438  }
1439  } // else point found
1440  }// no optimization
1441  } // for ( size_type k=0 ; k<memmapGdof[proc_id].size() ; ++k)
1442 
1443  return memory_localisationFail;
1444 }
1445 
1446 //-----------------------------------------------------------------------------------------------------------------//
1447 
1448 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
1449 std::list<boost::tuple<size_type,uint16_type> >
1450 OperatorInterpolation<DomainSpaceType, ImageSpaceType,
1451  IteratorRange,InterpType>::updateNoRelationMeshMPI_upWithOtherWorld(const std::vector< std::vector<size_type> > & memmapGdof,
1452  const std::vector< std::vector<uint16_type> > & memmapComp,
1453  const std::vector<std::vector<typename image_mesh_type::node_type> > & pointsSearched,
1454  const std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > & memmap_vertices,
1455  graph_ptrtype & sparsity_graph,
1456  std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
1457  std::vector<std::map<size_type,size_type> > & memory_col_globalProcessToGlobalCluster,
1458  std::vector<std::set<size_type> > & dof_searchWithProc,
1459  bool extrapolation_mode,
1460  extrapolation_memory_type & dof_extrapolationData)
1461 {
1462  std::list<boost::tuple<size_type,uint16_type> > memory_localisationFail;// gdof,comp
1463 
1464  auto const* imagedof = this->dualImageSpace()->dof().get();
1465  auto const* domaindof = this->domainSpace()->dof().get();
1466  auto const* domainbasis = this->domainSpace()->basis().get();
1467 
1468  const size_type proc_id = this->dualImageSpace()->worldComm()/*worldsComm()[0]*/.localRank();
1469  const size_type proc_id_image = this->dualImageSpace()->mesh()->worldComm().localRank();
1470  const size_type proc_id_domain = this->domainSpace()->mesh()->worldComm().localRank();
1471  const size_type nProc = this->dualImageSpace()->mesh()->worldComm().size();
1472  const size_type nProc_row = this->dualImageSpace()->mesh()->worldComm().localSize();
1473  const size_type nProc_col = this->domainSpace()->mesh()->worldComm().localSize();
1474  const size_type nProc_image = this->dualImageSpace()->mesh()->worldComm().localSize();
1475  const size_type nProc_domain = this->domainSpace()->mesh()->worldComm().localSize();
1476 
1477  // localisation tool with matrix node
1478  auto locTool = this->domainSpace()->mesh()->tool_localization();
1479  bool notUseOptLocTest = domain_mesh_type::nDim!=domain_mesh_type::nRealDim;
1480  if (notUseOptLocTest) locTool->kdtree()->nbNearNeighbor(domain_mesh_type::element_type::numPoints);
1481 
1482  matrix_node_type ptsReal( image_mesh_type::nRealDim, 1 );
1483  matrix_node_type ptsRef( domain_mesh_type::nDim , 1 );
1484  matrix_node_type MlocEval(domain_basis_type::nLocalDof*domain_basis_type::nComponents1,1);
1485  matrix_node_type verticesOfEltSearched;
1486 
1487  // random (just to start)
1488  //size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1489  size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1490  auto const& eltRandom = this->domainSpace()->mesh()->element(eltIdLocalised);
1491 
1492  std::vector<bool> dof_done( this->dualImageSpace()->nLocalDof(), false);
1493 
1494  // usefull container
1495  std::vector<size_type> pointsSearchedSizeWorld(this->dualImageSpace()->mesh()->worldComm().localComm().size());
1496  std::vector<typename image_mesh_type::node_type> dataToRecv(1);
1497  std::vector<uint16_type> dataToRecv_Comp(1,0);
1498  std::vector< std::vector< typename image_mesh_type::node_type > > dataToRecv_Vertices(1);
1499  std::vector<typename image_mesh_type::node_type> pointsRefFinded(1);
1500  std::vector<typename image_mesh_type::node_type> pointsBaryFinded(1);// extrapolation only
1501  std::vector<bool> pointsRefIsFinded(1,false);
1502  std::vector<int> pointsIdEltFinded(1,0);
1503  std::vector<std::vector<int> > pointsDofsColFinded(1,std::vector<int>(1,0));
1504  std::vector<std::vector<int> > pointsDofsGlobalClusterColFinded(1,std::vector<int>(1,0));
1505  std::vector<uint16_type> pointsComp(1,0);
1506 
1507 
1508  // Attention : marche que si les 2 worldcomms qui s'emboite (mon cas)
1509  std::vector<int> localMeshRankToWorldCommFusion_domain(nProc_col);
1510  mpi::all_gather( this->domainSpace()->mesh()->worldComm().localComm(),
1511  this->worldCommFusion().globalRank(),
1512  localMeshRankToWorldCommFusion_domain );
1513  std::vector<int> localMeshRankToWorldCommFusion_image(nProc_row);
1514  mpi::all_gather( this->dualImageSpace()->mesh()->worldComm().localComm(),
1515  this->worldCommFusion().globalRank(),
1516  localMeshRankToWorldCommFusion_image );
1517 
1518  std::vector<int> domainProcIsActive_fusion(this->worldCommFusion().globalSize());
1519  mpi::all_gather( this->worldCommFusion().globalComm(),
1520  (int)this->domainSpace()->worldComm().isActive(),
1521  domainProcIsActive_fusion );
1522  std::vector<int> imageProcIsActive_fusion(this->worldCommFusion().globalSize());
1523  mpi::all_gather( this->worldCommFusion().globalComm(),
1524  (int)this->dualImageSpace()->worldComm().isActive(),
1525  imageProcIsActive_fusion );
1526 
1527  int firstActiveProc_image=0;
1528  bool findFirstActive_image=false;
1529  while (!findFirstActive_image)
1530  {
1531  if (imageProcIsActive_fusion[firstActiveProc_image])
1532  {
1533  findFirstActive_image=true;
1534  }
1535  else ++firstActiveProc_image;
1536  }
1537  int firstActiveProc_domain=0;
1538  bool findFirstActive_domain=false;
1539  while (!findFirstActive_domain)
1540  {
1541  if (domainProcIsActive_fusion[firstActiveProc_domain])
1542  {
1543  findFirstActive_domain=true;
1544  }
1545  else ++firstActiveProc_domain;
1546  }
1547 
1548  for (int p=0;p<localMeshRankToWorldCommFusion_image.size(); ++p)
1549  {
1550  if (!this->dualImageSpace()->worldComm().isActive()) localMeshRankToWorldCommFusion_image[p]=p%nProc_image+firstActiveProc_image; // FAIRE COMMMUNICATION!!!!!
1551  }
1552  for (int p=0;p<localMeshRankToWorldCommFusion_domain.size(); ++p)
1553  {
1554  if (!this->domainSpace()->worldComm().isActive()) localMeshRankToWorldCommFusion_domain[p]=p%nProc_domain+firstActiveProc_domain; // FAIRE COMMMUNICATION!!!!!
1555  }
1556 
1557 #if 1
1558  std::vector<std::vector<int> > searchDistribution(nProc);
1559  for (int p=0;p<nProc_image;++p)
1560  {
1561  searchDistribution[p].clear();
1562  if ( proc_id_image == p && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1563  {
1564  for (int q=0;q<nProc_domain;++q)
1565  {
1566  if( pointsSearched[q].size()>0 && localMeshRankToWorldCommFusion_image[p]!=localMeshRankToWorldCommFusion_domain[q] )
1567  {
1568  searchDistribution[p].push_back(q);
1569  }
1570  }
1571  }
1572 
1573  mpi::broadcast( this->worldCommFusion().globalComm(), searchDistribution[p], localMeshRankToWorldCommFusion_image[p] );
1574  //mpi::broadcast( this->worldCommFusion().globalComm(), searchDistribution2, localMeshRankToWorldCommFusion_image[proc_id_image] );
1575  }
1576 
1577 #else //OLD
1578  // searchDistribution (no comm with ourself)
1579  std::vector<std::list<int> > searchDistribution(nProc);
1580  for (int p=0;p<nProc_image;++p)
1581  {
1582  searchDistribution[p].clear();
1583  for (int q=0;q<nProc_domain;++q)
1584  {
1585  //if (q!=p)
1586  if( (localMeshRankToWorldCommFusion_image[p])!=localMeshRankToWorldCommFusion_domain[q] )
1587  {
1588  searchDistribution[p].push_back(q);
1589  }
1590  }
1591  }
1592 #endif
1593 
1594 #if 0
1595  this->worldCommFusion().barrier();
1596  for (int p=0;p<this->worldCommFusion().globalSize();++p)
1597  {
1598  if (p==this->worldCommFusion().globalRank())
1599  {
1600  std::cout << "I am proc " << p << " ";// << "\n";
1601  for (int q=0;q<nProc_image;++q)
1602  {
1603  std::cout << "["<< q << "] : ";
1604  auto it_list = searchDistribution[q].begin();
1605  auto en_list = searchDistribution[q].end();
1606  for ( ; it_list!=en_list;++it_list)
1607  {
1608  std::cout << *it_list <<" ";
1609  }
1610  std::cout << std::endl;
1611  }
1612  }
1613  this->worldCommFusion().barrier();
1614  }
1615 #endif
1616 
1617 
1618 
1619 
1620 
1621  // Tag for mpi msg
1622  const int tag_X = 0, tag_Y = 1, tag_Z = 2, tag_IsFind = 3, tag_IdElt = 4, tag_DofsCol = 5, tag_DofsColGC = 6, tag_Comp = 7, tag_Points=8, tag_Bary=9, tag_Vertices=10;
1623  //------------------------------
1624  // proc after proc
1625  for (int proc=0;proc<nProc_row;++proc)
1626  {
1627  if ( proc_id_image == proc && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] ) // send info to rankLocalization
1628  {
1629  for (auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1630  it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
1631  {
1632  const int rankLocalization = *it_rankLocalization;
1633  const int rankToSend = localMeshRankToWorldCommFusion_domain[rankLocalization];
1634  this->worldCommFusion().globalComm().send(rankToSend,tag_Points,pointsSearched[rankLocalization]);
1635  this->worldCommFusion().globalComm().send(rankToSend,tag_Comp,memmapComp[rankLocalization]);
1636  if (InterpType::value==1)
1637  this->worldCommFusion().globalComm().send(rankToSend,tag_Vertices,memmap_vertices[rankLocalization]);
1638  }
1639  }
1640  else
1641  {
1642  for (auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1643  it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
1644  {
1645  const int rankLocalization = *it_rankLocalization;
1646  if ( proc_id_domain == rankLocalization && domainProcIsActive_fusion[this->worldCommFusion().globalRank()] ) // get request of proc
1647  {
1648  const int rankToRecv = localMeshRankToWorldCommFusion_image[proc];
1649  this->worldCommFusion().globalComm().recv(rankToRecv,tag_Points,dataToRecv);
1650  this->worldCommFusion().globalComm().recv(rankToRecv,tag_Comp,dataToRecv_Comp);
1651  if (InterpType::value==1)
1652  this->worldCommFusion().globalComm().recv(rankToRecv,tag_Vertices,dataToRecv_Vertices);
1653 
1654  const size_type nDataRecv = dataToRecv.size();
1655  // init container
1656  pointsRefFinded.resize(nDataRecv);
1657  pointsRefIsFinded.resize(nDataRecv);std::fill(pointsRefIsFinded.begin(),pointsRefIsFinded.end(),false);
1658  pointsIdEltFinded.resize(nDataRecv);
1659  if (extrapolation_mode) pointsBaryFinded.resize(nDataRecv);
1660  pointsDofsColFinded.resize(nDataRecv,std::vector<int>(1,0));
1661  pointsDofsGlobalClusterColFinded.resize(nDataRecv,std::vector<int>(1,0));
1662  // iterate on points
1663  for (size_type k=0;k<nDataRecv;++k)
1664  {
1665  // get real point to search
1666  ublas::column(ptsReal,0) = dataToRecv[k];
1667  // vertice with conforme case
1668  if (InterpType::value==1)
1669  {
1670  const uint16_type sizeVertices = dataToRecv_Vertices[k].size();
1671  verticesOfEltSearched.resize( image_mesh_type::nRealDim,sizeVertices);
1672  for ( uint16_type v=0;v<sizeVertices;++v)
1673  ublas::column(verticesOfEltSearched,v)=dataToRecv_Vertices[k][v];
1674  }
1675  else // random
1676  verticesOfEltSearched = eltRandom.vertices();
1677  // search process
1678  if (notUseOptLocTest) eltIdLocalised=invalid_size_type_value;
1679  auto resLocalisation = locTool->run_analysis(ptsReal,eltIdLocalised,verticesOfEltSearched,mpl::int_<interpolation_type::value>());
1680  if (resLocalisation.template get<0>()[0]) // is find
1681  {
1682  eltIdLocalised = resLocalisation.template get<1>();
1683  //-------------------------------------------------------------------------
1684  if (!this->interpolationType().componentsAreSamePoint() )// all component
1685  //-------------------------------------------------------------------------
1686  {
1687  const uint16_type comp=dataToRecv_Comp[k];
1688  pointsRefIsFinded[k]=true;
1689  pointsIdEltFinded[k]=eltIdLocalised;
1690  // get point in reference element
1691  auto itanal = locTool->result_analysis_begin();
1692  auto const itanal_end = locTool->result_analysis_end();
1693  for ( ;itanal!=itanal_end;++itanal)
1694  {
1695  auto const itL=itanal->second.begin();
1696  //ublas::column( ptsRef, 0 ) = boost::get<1>(*itL);
1697  pointsRefFinded[k] = boost::get<1>(*itL);
1698  }
1699  // get global process dof and global cluster dof : column data map
1700  pointsDofsColFinded[k].resize(domain_basis_type::nLocalDof);
1701  pointsDofsGlobalClusterColFinded[k].resize(domain_basis_type::nLocalDof);
1702  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1703  {
1704  const auto j_gdof = boost::get<0>(domaindof->localToGlobal( eltIdLocalised,jloc,comp ));
1705  pointsDofsColFinded[k][jloc] = j_gdof;
1706  pointsDofsGlobalClusterColFinded[k][jloc] = domaindof->mapGlobalProcessToGlobalCluster()[j_gdof];
1707  }
1708  if (extrapolation_mode)
1709  {
1710  auto const& eltExtrapoled = this->domainSpace()->mesh()->element(eltIdLocalised);
1711  auto const verticesExtrapoled = eltExtrapoled.vertices();
1712  typename image_mesh_type::node_type bary(verticesExtrapoled.size1());
1713  for (int qi=0;qi<verticesExtrapoled.size1();++qi) bary(qi)=0;//important
1714  for (int qj=0;qj<verticesExtrapoled.size2();++qj)
1715  {
1716  bary(0) += ublas::column(verticesExtrapoled,qj)(0);
1717  if (verticesExtrapoled.size1()>1) bary(1) += ublas::column(verticesExtrapoled,qj)(1);
1718  if (verticesExtrapoled.size1()>2) bary(2) += ublas::column(verticesExtrapoled,qj)(2);
1719  }
1720  bary(0) /= verticesExtrapoled.size2();
1721  if (verticesExtrapoled.size1()>1) bary(1) /= verticesExtrapoled.size2();
1722  if (verticesExtrapoled.size1()>2) bary(2) /= verticesExtrapoled.size2();
1723  pointsBaryFinded[k]=bary;
1724  }
1725  }
1726  //-------------------------------------------------------------------------
1727  else // components optimization
1728  //-------------------------------------------------------------------------
1729  {
1730  auto const ptRefFinded = locTool->result_analysis_begin()->second.begin()->template get<1>();
1731  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1732  {
1733  pointsRefIsFinded[k+comp]=true;
1734  pointsIdEltFinded[k+comp]=eltIdLocalised;
1735  // get point in reference element
1736  pointsRefFinded[k+comp] = ptRefFinded;
1737 
1738  pointsDofsColFinded[k+comp].resize(domain_basis_type::nLocalDof);
1739  pointsDofsGlobalClusterColFinded[k+comp].resize(domain_basis_type::nLocalDof);
1740  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1741  {
1742  const auto j_gdof = boost::get<0>(domaindof->localToGlobal( eltIdLocalised,jloc,comp ));
1743  pointsDofsColFinded[k+comp][jloc] = j_gdof;
1744  pointsDofsGlobalClusterColFinded[k+comp][jloc] = domaindof->mapGlobalProcessToGlobalCluster()[j_gdof];
1745  }
1746  }
1747 
1748  if (extrapolation_mode)
1749  {
1750  auto const& eltExtrapoled = this->domainSpace()->mesh()->element(eltIdLocalised);
1751  auto const verticesExtrapoled = eltExtrapoled.vertices();
1752  typename image_mesh_type::node_type bary(verticesExtrapoled.size1());
1753  for (int qi=0;qi<verticesExtrapoled.size1();++qi) bary(qi)=0;//important
1754  for (int qj=0;qj<verticesExtrapoled.size2();++qj)
1755  {
1756  bary(0) += ublas::column(verticesExtrapoled,qj)(0);
1757  if (verticesExtrapoled.size1()>1) bary(1) += ublas::column(verticesExtrapoled,qj)(1);
1758  if (verticesExtrapoled.size1()>2) bary(2) += ublas::column(verticesExtrapoled,qj)(2);
1759  }
1760  bary(0) /= verticesExtrapoled.size2();
1761  if (verticesExtrapoled.size1()>1) bary(1) /= verticesExtrapoled.size2();
1762  if (verticesExtrapoled.size1()>2) bary(2) /= verticesExtrapoled.size2();
1763  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1764  pointsBaryFinded[k+comp]=bary;
1765  }
1766  // change increment in for
1767  k=k+image_basis_type::nComponents-1;
1768  }
1769 
1770  //std::cout << "F";
1771  }
1772  else // Not Find!
1773  {
1774  //memory_localisationFail
1775  //std::cout << "NOT FIND"<<std::endl;
1776  }
1777 
1778  } // for (size_type k=0;k<nDataRecv;++k)
1779 
1780  const int rankToSend = localMeshRankToWorldCommFusion_image[proc];
1781  this->worldCommFusion().globalComm().send(rankToSend,tag_Points,pointsRefFinded);
1782  this->worldCommFusion().globalComm().send(rankToSend,tag_IsFind,pointsRefIsFinded);
1783  this->worldCommFusion().globalComm().send(rankToSend,tag_IdElt,pointsIdEltFinded);
1784  this->worldCommFusion().globalComm().send(rankToSend,tag_DofsCol,pointsDofsColFinded);
1785  this->worldCommFusion().globalComm().send(rankToSend,tag_DofsColGC,pointsDofsGlobalClusterColFinded);
1786  if (extrapolation_mode)
1787  this->worldCommFusion().globalComm().send(rankToSend,tag_Bary,pointsBaryFinded);
1788 
1789  }
1790  } // for (auto it_rankLocalization=...
1791  }
1792 
1793 
1794 
1795 
1796  if ( proc_id_image == proc && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1797  {
1798  for (auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1799  it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
1800  {
1801  const int rankLocalization = *it_rankLocalization;
1802  const int rankToRecv = localMeshRankToWorldCommFusion_domain[rankLocalization];
1803  this->worldCommFusion().globalComm().recv(rankToRecv,tag_Points,pointsRefFinded);
1804  this->worldCommFusion().globalComm().recv(rankToRecv,tag_IsFind,pointsRefIsFinded);
1805  this->worldCommFusion().globalComm().recv(rankToRecv,tag_IdElt,pointsIdEltFinded);
1806  this->worldCommFusion().globalComm().recv(rankToRecv,tag_DofsCol,pointsDofsColFinded);
1807  this->worldCommFusion().globalComm().recv(rankToRecv,tag_DofsColGC,pointsDofsGlobalClusterColFinded);
1808  if (extrapolation_mode)
1809  this->worldCommFusion().globalComm().recv(rankToRecv,tag_Bary,pointsBaryFinded);
1810 
1811  const int rankRecv=rankLocalization;
1812  for ( int k=0;k<pointsRefFinded.size();++k)
1813  {
1814  if (!pointsRefIsFinded[k])
1815  {
1816  memory_localisationFail.push_back(boost::make_tuple(memmapGdof[rankLocalization][k],memmapComp[rankLocalization][k]));
1817  const auto i_gdof = memmapGdof[rankLocalization][k];
1818  dof_searchWithProc[i_gdof].insert(rankRecv);
1819  }
1820  else
1821  {
1822  const auto i_gdof = memmapGdof[rankLocalization][k];
1823  if (!dof_done[i_gdof])
1824  {
1825 
1826  if (extrapolation_mode)
1827  {
1828  auto const bary = pointsBaryFinded[k];
1829  auto const theRefPtExtrap = pointsRefFinded[k];
1830  const auto comp = memmapComp[rankLocalization][k];
1831  std::vector<std::pair<size_type,size_type> > j_gdofs(domain_basis_type::nLocalDof);
1832  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1833  {
1834  const size_type j_gdof = pointsDofsColFinded[k][jloc];
1835  const size_type j_gdof_gc = pointsDofsGlobalClusterColFinded[k][jloc];
1836  j_gdofs[jloc]=std::make_pair(j_gdof,j_gdof_gc);
1837  }
1838  dof_extrapolationData[i_gdof].push_back(boost::make_tuple(rankLocalization,bary,theRefPtExtrap,j_gdofs,comp));
1839  }
1840  else
1841  {
1842  //std::cout << "T";
1843  ublas::column( ptsRef, 0 ) = pointsRefFinded[k];
1844  //evalute point on the reference element
1845  MlocEval = domainbasis->evaluate( ptsRef );
1846 
1847  const auto comp = memmapComp[rankLocalization][k];
1848  const size_type myidElt = pointsIdEltFinded[k];
1849  const auto ig1 = imagedof->mapGlobalProcessToGlobalCluster()[i_gdof];
1850  const auto theproc = imagedof->procOnGlobalCluster(ig1);
1851  auto& row = sparsity_graph->row(ig1);
1852  row.template get<0>() = theproc;
1853  row.template get<1>() = ig1 - imagedof->firstDofGlobalCluster(theproc);
1854 
1855  for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1856  {
1857  //get global process dof
1858  const size_type j_gdof = pointsDofsColFinded[k][jloc];
1859  //get global cluster dof
1860  const size_type j_gdof_gc = pointsDofsGlobalClusterColFinded[k][jloc];
1861  // up graph
1862  row.template get<2>().insert(j_gdof_gc);
1863  // get value
1864  const auto v = MlocEval( domain_basis_type::nComponents1*jloc +
1865  comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1866  comp, 0 );
1867 #if 1
1868  // save value
1869  memory_valueInMatrix[i_gdof].push_back(boost::make_tuple(rankRecv,j_gdof,v));
1870  // usefull to build datamap
1871  memory_col_globalProcessToGlobalCluster[rankRecv][j_gdof]=j_gdof_gc;
1872 #endif
1873  }
1874  // dof ok : not anymore localise
1875  dof_searchWithProc[i_gdof].insert(rankRecv);
1876  }
1877  dof_done[i_gdof]=true;
1878  } // if (!dof_done[i_gdof])
1879  }
1880  }
1881  }
1882  } // if ( proc_id_image == proc && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1883  }
1884 
1885 
1886 
1887 
1888  return memory_localisationFail;
1889 }
1890 
1891 template<typename DomainSpaceType, typename ImageSpaceType,typename IteratorRange,typename InterpType>
1892 boost::tuple<std::vector< std::vector<size_type> >, std::vector< std::vector<uint16_type> >,
1893  std::vector<std::vector<typename OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::image_mesh_type::node_type> >,
1894  std::vector<std::vector< std::vector<typename OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::image_mesh_type::node_type > > > >
1895 OperatorInterpolation<DomainSpaceType, ImageSpaceType,
1896  IteratorRange,InterpType>::updateNoRelationMeshMPI_pointDistribution(const std::vector< std::list<boost::tuple<int,size_type,double> > > & memory_valueInMatrix,
1897  std::vector<std::set<size_type> > & dof_searchWithProc)
1898 {
1899  //std::cout << " pointDistribution--1--- " << this->domainSpace()->mesh()->worldComm().godRank() << std::endl;
1900  //const size_type proc_id = this->dualImageSpace()->worldsComm()[0].localRank();
1901  //const size_type nProc = this->dualImageSpace()->mesh()->worldComm().size();
1902  //const size_type nProc_image = this->dualImageSpace()->mesh()->worldComm().localSize();
1903  const size_type nProc_domain = this->domainSpace()->mesh()->worldComm().localSize();
1904 
1905  auto const* imagedof = this->dualImageSpace()->dof().get();
1906  iterator_type it, en;
1907 
1908  auto locTool = this->domainSpace()->mesh()->tool_localization();
1909 
1910  std::vector<bool> dof_done( this->dualImageSpace()->nLocalDof(), false);
1911  std::vector< std::list<boost::tuple<size_type,uint16_type> > > memSetGdofAndComp( nProc_domain );
1912  std::vector< std::list<matrix_node_type> > memSetVertices_conformeInterp( nProc_domain );
1913 
1914  // Warning communication!!
1915  std::vector<typename image_mesh_type::node_type> vecBarycenter(nProc_domain);
1916  mpi::all_gather( this->domainSpace()->mesh()->worldComm().localComm(),
1917  locTool->barycenter(),
1918  vecBarycenter );
1919  /*std::cout << " proc " << this->domainSpace()->mesh()->worldComm().localRank()
1920  << " procFuion " << this->worldCommFusion().globalRank()
1921  << " bary " << locTool->barycenter()
1922  << std::endl;*/
1923 
1924 
1925  double distanceMin=0,distance=0,distanceSquare=0;
1926  int procForPt=0;
1927 
1928  if ( this->dualImageSpace()->worldComm().isActive() )
1929  {
1930 
1931  auto itListRange = M_listRange.begin();
1932  auto const enListRange = M_listRange.end();
1933  for ( ; itListRange!=enListRange ; ++itListRange)
1934  {
1935  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
1936  for ( ; it!=en;++it )
1937  {
1938  for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
1939  {
1940  for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1941  {
1942  const auto gdof = boost::get<0>(imagedof->localToGlobal( *it, iloc, comp ));
1943  if (!dof_done[gdof] && memory_valueInMatrix[gdof].size()==0)
1944  {
1945  // the dof point
1946  const auto imagePoint = imagedof->dofPoint(gdof).template get<0>();
1947 
1948  if (this->interpolationType().searchWithCommunication()) // mpi communication
1949  {
1950  distanceMin=INT_MAX;
1951  for ( int proc=0 ; proc<nProc_domain; ++proc)
1952  {
1953  const auto bary = vecBarycenter[proc];
1954  distanceSquare = std::pow(imagePoint(0)-bary(0),2);
1955  if (bary.size()>1) distanceSquare += std::pow(imagePoint(1)-bary(1),2);
1956  if (bary.size()>2) distanceSquare += std::pow(imagePoint(2)-bary(2),2);
1957  distance = std::sqrt( distanceSquare );
1958  if (distance<distanceMin && dof_searchWithProc[gdof].find(proc)==dof_searchWithProc[gdof].end() )
1959  {
1960  procForPt = proc;
1961  distanceMin=distance;
1962  }
1963  }
1964  memSetGdofAndComp[procForPt].push_back(boost::make_tuple(gdof,comp));
1965  if (InterpType::value==1)
1966  memSetVertices_conformeInterp[procForPt].push_back(it->vertices());
1967  }
1968  else // only with myself
1969  {
1970  memSetGdofAndComp[this->domainSpace()->worldComm().globalRank()].push_back(boost::make_tuple(gdof,comp));
1971  if (InterpType::value==1) // conforme case
1972  {
1973  memSetVertices_conformeInterp[this->domainSpace()->worldComm().globalRank()].push_back(it->vertices());
1974  }
1975  }
1976 
1977  dof_done[gdof]=true;
1978  }
1979  }
1980  }
1981  }
1982  } //for ( ; itListRange!=enListRange ; ++itListRange)
1983  } // isActive
1984 
1985  // memory map (loc index pt) -> global dofs
1986  std::vector< std::vector<size_type> > memmapGdof( nProc_domain );
1987  // memory map (loc index pt) -> comp
1988  std::vector< std::vector<uint16_type> > memmapComp( nProc_domain );
1989  // points to lacalize
1990  std::vector<std::vector<typename image_mesh_type::node_type> > pointsSearched( nProc_domain );
1991 
1992  std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > memmap_vertices( nProc_domain );
1993 
1994  for (int proc=0; proc<nProc_domain;++proc)
1995  {
1996  const size_type nData = memSetGdofAndComp[proc].size();
1997  memmapGdof[proc].resize(nData);
1998  memmapComp[proc].resize(nData);
1999  pointsSearched[proc].resize(nData);
2000  // conforme case
2001  if(InterpType::value==1) memmap_vertices[proc].resize(nData);
2002 
2003  auto it_GdofAndComp = memSetGdofAndComp[proc].begin();
2004  auto it_vertices = memSetVertices_conformeInterp[proc].begin();
2005  for (int k=0 ; k<nData ; ++k, ++it_GdofAndComp)
2006  {
2007  memmapGdof[proc][k]=it_GdofAndComp->template get<0>();//gdof;
2008  memmapComp[proc][k]=it_GdofAndComp->template get<1>();//comp
2009  pointsSearched[proc][k]=imagedof->dofPoint(it_GdofAndComp->template get<0>()).template get<0>();//node
2010  if(InterpType::value==1) // conforme case
2011  {
2012  memmap_vertices[proc][k].resize(it_vertices->size2());
2013  for (uint16_type v=0;v<it_vertices->size2();++v)
2014  memmap_vertices[proc][k][v]=ublas::column(*it_vertices,v);
2015  ++it_vertices;
2016  }
2017  }
2018  }
2019 
2020  return boost::make_tuple(memmapGdof,memmapComp,pointsSearched,memmap_vertices);
2021 }
2022 
2023 #endif // WITH MPI
2024 
2025 
2026 //-----------------------------------------------------------------------------------------------------------------//
2027 //-----------------------------------------------------------------------------------------------------------------//
2028 //-----------------------------------------------------------------------------------------------------------------//
2029 
2030 namespace detail
2031 {
2032 template <typename RangeType>
2033 struct opinterprangetype
2034 {
2035  typedef typename mpl::if_< boost::is_std_list<RangeType>,
2036  mpl::identity<RangeType>,
2037  mpl::identity<std::list<RangeType> > >::type::type::value_type type;
2038 };
2039 }
2040 
2041 template<typename DomainSpaceType, typename ImageSpaceType, typename IteratorRange, typename InterpType >
2042 boost::shared_ptr<OperatorInterpolation<DomainSpaceType, ImageSpaceType,typename Feel::detail::opinterprangetype<IteratorRange>::type,InterpType> >
2043 opInterp( boost::shared_ptr<DomainSpaceType> const& domainspace,
2044  boost::shared_ptr<ImageSpaceType> const& imagespace,
2045  IteratorRange const& r,
2046  typename OperatorInterpolation<DomainSpaceType, ImageSpaceType,typename Feel::detail::opinterprangetype<IteratorRange>::type,InterpType>::backend_ptrtype const& backend,
2047  InterpType const& interptype,
2048  bool ddmethod
2049  )
2050 {
2051  typedef OperatorInterpolation<DomainSpaceType, ImageSpaceType,typename Feel::detail::opinterprangetype<IteratorRange>::type,InterpType> operatorinterpolation_type;
2052 
2053  boost::shared_ptr<operatorinterpolation_type> opI( new operatorinterpolation_type( domainspace,imagespace,r,backend,interptype,ddmethod ) );
2054 
2055  return opI;
2056 }
2057 
2058 
2059 
2060 template<typename Args>
2061 struct compute_opInterpolation_return
2062 {
2063  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
2064  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
2065 
2066  typedef typename boost::remove_const<
2067  typename boost::remove_reference<
2068  typename parameter::binding<Args,
2069  tag::range,
2070  typename OperatorInterpolation<domain_space_type, image_space_type>::range_iterator
2071  >::type >::type >::type iterator_base_range_type;
2072 
2073  typedef typename Feel::detail::opinterprangetype<iterator_base_range_type>::type iterator_range_type;
2074 
2075  typedef typename boost::remove_const<
2076  typename boost::remove_reference<
2077  typename parameter::binding<Args,
2078  tag::type,
2079  InterpolationNonConforme
2080  >::type >::type >::type interpolation_type;
2081 
2082 
2083  typedef boost::shared_ptr<OperatorInterpolation<domain_space_type, image_space_type,iterator_range_type,interpolation_type> > type;
2084 };
2085 
2086 BOOST_PARAMETER_FUNCTION(
2087  ( typename compute_opInterpolation_return<Args>::type ), // 1. return type
2088  opInterpolation, // 2. name of the function template
2089  tag, // 3. namespace of tag types
2090  ( required
2091  ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
2092  ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
2093  ) // required
2094  ( optional
2095  ( range, *, elements( imageSpace->mesh() ) )
2096  ( backend, *, Backend<typename compute_opInterpolation_return<Args>::domain_space_type::value_type>::build() )
2097  ( type, *, InterpolationNonConforme() )
2098  ( ddmethod, (bool), false )
2099  ) // optionnal
2100 )
2101 {
2102  Feel::detail::ignore_unused_variable_warning( args );
2103 
2104  return opInterp( domainSpace,imageSpace,range,backend,type,ddmethod );
2105 
2106 } // opInterpolation
2107 
2108 
2109 
2110 
2111 } // Feel
2112 #endif /* __OperatorInterpolation_H */
Linear Operator between function spaces, represented by a matrix.
Definition: operatorlinear.hpp:43
elements_type const & elements() const
Definition: elements.hpp:355
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
Global interpolation operator.
Definition: operatorinterpolation.hpp:126
OperatorInterpolation()
Definition: operatorinterpolation.hpp:204
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
OperatorInterpolation(OperatorInterpolation const &oi)
Definition: operatorinterpolation.hpp:235
Graph representation of the Compressed Sparse Row format.
Definition: graphcsr.hpp:57

Generated on Sun Dec 22 2013 13:11:09 for Feel++ by doxygen 1.8.5