30 #ifndef __OperatorInterpolation_H
31 #define __OperatorInterpolation_H 1
38 class InterpolationTypeBase
41 InterpolationTypeBase(
bool useComm=
true,
bool compAreSamePt=
true,
bool onlyLocalizeOnBoundary=
false,
int nbNearNeighborInKdTree=15 )
43 M_searchWithCommunication( useComm ),
44 M_componentsAreSamePoint( compAreSamePt ),
45 M_onlyLocalizeOnBoundary( onlyLocalizeOnBoundary ),
46 M_nbNearNeighborInKdTree( nbNearNeighborInKdTree )
49 InterpolationTypeBase( InterpolationTypeBase
const& a)
51 M_searchWithCommunication( a.M_searchWithCommunication ),
52 M_componentsAreSamePoint( a.M_componentsAreSamePoint ),
53 M_onlyLocalizeOnBoundary( a.M_onlyLocalizeOnBoundary ),
54 M_nbNearNeighborInKdTree( a.M_nbNearNeighborInKdTree )
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; }
64 bool M_searchWithCommunication;
65 bool M_componentsAreSamePoint;
66 bool M_onlyLocalizeOnBoundary;
67 int M_nbNearNeighborInKdTree;
70 struct InterpolationNonConforme :
public InterpolationTypeBase
72 static const uint16_type value=0;
74 InterpolationNonConforme(
bool useComm=
true,
bool compAreSamePt=
true,
bool onlyLocalizeOnBoundary=
false,
int nbNearNeighborInKdTree=15 )
76 InterpolationTypeBase(useComm,compAreSamePt, onlyLocalizeOnBoundary,nbNearNeighborInKdTree)
80 struct InterpolationConforme :
public InterpolationTypeBase
82 static const uint16_type value=1;
84 InterpolationConforme(
bool useComm=
true,
bool compAreSamePt=
true,
bool onlyLocalizeOnBoundary=
false,
int nbNearNeighborInKdTree=15 )
86 InterpolationTypeBase(useComm,compAreSamePt,onlyLocalizeOnBoundary,nbNearNeighborInKdTree)
94 template <
typename EltType >
96 idElt( EltType & elt,mpl::size_t<MESH_ELEMENTS> )
101 template <
typename EltType >
103 idElt( EltType & elt,mpl::size_t<MESH_FACES> )
105 return elt.element0().id();
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 >
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;
146 typedef typename super::backend_ptrtype backend_ptrtype;
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;
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;
167 typedef typename dual_image_space_type::dof_type dof_type;
170 typedef typename dual_image_space_type::basis_type image_basis_type;
171 typedef typename domain_space_type::basis_type domain_basis_type;
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;
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;
185 typedef InterpType interpolation_type;
189 typedef boost::shared_ptr<graph_type> graph_ptrtype;
192 typedef typename matrix_node<typename image_mesh_type::value_type>::type matrix_node_type;
212 dual_image_space_ptrtype
const& imagespace,
213 backend_ptrtype
const& backend,
214 InterpType
const& interptype,
215 bool ddmethod=
false);
218 dual_image_space_ptrtype
const& imagespace,
219 IteratorRange
const& r,
220 backend_ptrtype
const& backend,
221 InterpType
const& interptype,
222 bool ddmethod=
false);
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);
238 M_listRange( oi.M_listRange ),
239 M_WorldCommFusion( oi.M_WorldCommFusion ),
240 M_interptype( oi.M_interptype )
258 WorldComm
const& worldCommFusion()
const {
return M_WorldCommFusion; }
260 InterpType
const& interpolationType()
const {
return M_interptype; }
262 bool isDomainMeshRelatedToImageMesh()
const {
return this->domainSpace()->mesh()->isSubMeshFrom( this->dualImageSpace()->mesh() ); }
264 bool isImageMeshRelatedToDomainMesh()
const {
return this->dualImageSpace()->mesh()->isSubMeshFrom( this->domainSpace()->mesh() ); }
286 virtual void update();
290 void updateSameMesh();
291 void updateNoRelationMesh();
292 #if defined(FEELPP_ENABLE_MPI_MODE)
293 void updateNoRelationMeshMPI();
294 void updateNoRelationMeshMPI_run(
bool buildNonZeroMatrix=
true);
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> >,
301 > > > extrapolation_memory_type;
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);
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);
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);
338 std::list<range_iterator> M_listRange;
339 WorldComm M_WorldCommFusion;
340 InterpType M_interptype;
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,
350 super( domainspace, imagespace, backend, false ),
352 M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
353 M_interptype(interptype)
355 M_listRange.push_back(
elements( imagespace->mesh() ) );
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,
368 super( domainspace, imagespace, backend, false ),
370 M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
371 M_interptype(interptype)
373 M_listRange.push_back( r );
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,
385 super( domainspace, imagespace, backend, false ),
387 M_WorldCommFusion( (ddmethod) ? this->domainSpace()->worldComm() : this->domainSpace()->worldComm()+this->dualImageSpace()->worldComm() ),
388 M_interptype(interptype)
394 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
396 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::update()
398 if ( this->dualImageSpace()->mesh()->numElements() == 0 )
401 this->matPtr() = this->backend()->newZeroMatrix( this->domainSpace()->dofOnOff(),
402 this->dualImageSpace()->dofOn() );
410 if ( this->dualImageSpace()->mesh()->isRelatedTo( this->domainSpace()->mesh() ) &&
411 boost::is_same<domain_mesh_type,image_mesh_type>::type::value
414 VLOG(2) <<
"OperatorInterpolation: use same mesh\n";
415 VLOG(2) <<
"isDomainMeshRelatedToImageMesh: " << isDomainMeshRelatedToImageMesh() <<
"\n";
416 VLOG(2) <<
"isImageMeshRelatedToDomainMesh: " << isImageMeshRelatedToDomainMesh() <<
"\n";
417 this->updateSameMesh();
421 #if defined(FEELPP_ENABLE_MPI_MODE)
422 if ( this->dualImageSpace()->worldComm().localSize() > 1 ||
423 this->domainSpace()->worldComm().localSize() > 1 )
424 this->updateNoRelationMeshMPI();
427 this->updateNoRelationMesh();
429 this->updateNoRelationMesh();
441 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
443 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateSameMesh()
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 );
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 );
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() ) );
467 graph_ptrtype sparsity_graph(
new graph_type( this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
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();
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 );
486 auto const& Mloc = domainbasis->evaluate( imagebasis->dual().points() );
488 DVLOG(2) <<
"[interpolate] Same mesh but not same space\n";
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() );
493 auto itListRange = M_listRange.begin();
494 auto const enListRange = M_listRange.end();
495 for ( ; itListRange!=enListRange ; ++itListRange)
497 iterator_type it, en;
498 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
499 for ( ; it != en; ++ it )
501 auto idElem = detailsup::idElt( *it,idim_type() );
502 auto domain_eid = idElem;
503 if ( image_related_to_domain )
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";
508 if( domain_related_to_image )
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";
517 for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
519 for ( uint16_type comp = 0; comp < image_basis_type::nComponents; ++comp )
521 size_type i = boost::get<0>( imagedof->localToGlobal( *it, iloc, comp ) );
525 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
527 const auto theproc = imagedof->worldComm().localRank();
529 const auto ig1 = imagedof->mapGlobalProcessToGlobalCluster()[i];
530 const auto theproc = imagedof->procOnGlobalCluster( ig1 );
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;
538 uint16_type ilocprime=imagedof->localDofInElement( *it, iloc, comp ) ;
540 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
543 const size_type j = boost::get<0>( domaindof->localToGlobal( domain_eid, jloc, comp ) );
545 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
546 row.template get<2>().insert( j );
548 row.template get<2>().insert( domaindof->mapGlobalProcessToGlobalCluster()[j] );
551 const value_type v = Mloc( domain_basis_type::nComponents1*jloc +
552 comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
556 memory_valueInMatrix[i].push_back( std::make_pair( j,v ) );
569 sparsity_graph->close();
572 this->matPtr() = this->backend()->newMatrix( this->domainSpace()->dofOnOff(),
573 this->dualImageSpace()->dofOn(),
578 for (
size_type idx_i=0 ; idx_i<nrow_dof_on_proc; ++idx_i )
580 for (
auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j )
582 this->matPtr()->set( idx_i,it_j->first,it_j->second );
591 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
593 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMesh()
595 DVLOG(2) <<
"[interpolate] different meshes\n";
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 );
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() ) );
610 graph_ptrtype sparsity_graph(
new graph_type( this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
613 auto const* imagedof = this->dualImageSpace()->dof().get();
614 auto const* domaindof = this->domainSpace()->dof().get();
615 auto const* domainbasis = this->domainSpace()->basis().get();
619 auto locTool = this->domainSpace()->mesh()->tool_localization();
620 if ( this->interpolationType().onlyLocalizeOnBoundary() ) locTool->updateForUseBoundaryFaces();
621 else locTool->updateForUse();
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);
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 );
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() );
646 auto itListRange = M_listRange.begin();
647 auto const enListRange = M_listRange.end();
648 for ( ; itListRange!=enListRange ; ++itListRange)
650 iterator_type it, en;
651 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
652 for ( ; it != en; ++ it )
654 for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
656 for ( uint16_type comp = 0; comp < image_basis_type::nComponents; ++comp )
658 const auto gdof = boost::get<0>(imagedof->localToGlobal( *it, iloc, comp ));
663 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
664 const auto ig1 = gdof;
665 const auto theproc = imagedof->worldComm().localRank();
667 const auto ig1 = imagedof->mapGlobalProcessToGlobalCluster()[gdof];
668 const auto theproc = imagedof->procOnGlobalCluster( ig1 );
670 auto& row = sparsity_graph->row(ig1);
671 row.template get<0>() = theproc;
672 row.template get<1>() = gdof;
675 ublas::column(ptsReal,0 ) = boost::get<0>(imagedof->dofPoint(gdof));
679 eltIdLocalised = locTool->run_analysis(ptsReal,eltIdLocalised,it->vertices(),mpl::int_<interpolation_type::value>()).template get<1>();
682 itanal = locTool->result_analysis_begin();
683 itanal_end = locTool->result_analysis_end();
684 for ( ;itanal!=itanal_end;++itanal)
686 itL=itanal->second.begin();
687 ublas::column( ptsRef, 0 ) = boost::get<1>( *itL );
689 MlocEval = domainbasis->evaluate( ptsRef );
691 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
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
699 #if !defined(FEELPP_ENABLE_MPI_MODE) // NOT MPI
700 row.template get<2>().insert( j );
702 row.template get<2>().insert( domaindof->mapGlobalProcessToGlobalCluster()[j] );
704 memory_valueInMatrix[gdof].push_back( std::make_pair( j,v ) );
717 sparsity_graph->close();
720 this->matPtr() = this->backend()->newMatrix( this->domainSpace()->dofOnOff(),
721 this->dualImageSpace()->dofOn(),
725 for (
size_type idx_i=0 ; idx_i<this->dualImageSpace()->nLocalDof() ;++idx_i)
727 for (
auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j)
729 this->matPtr()->set(idx_i,it_j->first,it_j->second);
739 #if defined(FEELPP_ENABLE_MPI_MODE) // WITH MPI
741 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
743 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMeshMPI()
747 auto testCommActivities_image=this->dualImageSpace()->worldComm().hasMultiLocalActivity();
749 if (testCommActivities_image.template get<0>())
753 std::vector<int> saveActivities_image = this->dualImageSpace()->worldComm().activityOnWorld();
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 )
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);
766 this->dualImageSpace()->worldComm().setIsActive(saveActivities_image);
767 this->dualImageSpace()->mapOn().worldComm().setIsActive(saveActivities_image);
768 this->dualImageSpace()->mapOnOff().worldComm().setIsActive(saveActivities_image);
773 if ( !this->dualImageSpace()->worldComm().isActive() && !this->domainSpace()->worldComm().isActive() )
775 this->matPtr() = this->backend()->newZeroMatrix( this->domainSpace()->dofOnOff(),
776 this->dualImageSpace()->dofOn() );
780 this->updateNoRelationMeshMPI_run(
true);
786 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
788 OperatorInterpolation<DomainSpaceType, ImageSpaceType,IteratorRange,InterpType>::updateNoRelationMeshMPI_run(
bool buildNonZeroMatrix)
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 );
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,
816 mpi::all_gather( this->domainSpace()->worldComm().globalComm(),
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() );
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() ) );
828 graph_ptrtype sparsity_graph(
new graph_type(this->dualImageSpace()->dof(), this->domainSpace()->dof() ) );
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() )
836 new_nLocalDofWithoutGhost = this->domainSpace()->mapOnOff().nLocalDofWithoutGhost();
840 if (this->dualImageSpace()->worldComm().globalRank()==this->dualImageSpace()->worldComm().masterRank()) new_nLocalDofWithoutGhost+=new_nLocalDofWithoutGhostMiss;
843 size_type new_firstdofcol=0,new_lastdofcol=new_nLocalDofWithoutGhost-1;
844 bool findMyProc=
false;
848 if (currentProc==this->dualImageSpace()->worldComm().globalRank())
852 else if (currentProc==this->dualImageSpace()->worldComm().masterRank())
854 new_firstdofcol+=new_nLocalDofWithoutGhost_tempp+new_nLocalDofWithoutGhostMiss;
855 new_lastdofcol+=new_nLocalDofWithoutGhost_tempp+new_nLocalDofWithoutGhostMiss;
859 new_firstdofcol+=new_nLocalDofWithoutGhost_tempp;
860 new_lastdofcol+=new_nLocalDofWithoutGhost_tempp;
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();
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());
877 extrapolation_memory_type dof_extrapolationData(this->dualImageSpace()->nLocalDof());
880 auto locTool = this->domainSpace()->mesh()->tool_localization();
881 bool doExtrapolationAtStart = locTool->doExtrapolation();
883 locTool->kdtree()->nbNearNeighbor(this->interpolationType().nbNearNeighborInKdTree());
885 if ( this->interpolationType().onlyLocalizeOnBoundary() ) locTool->updateForUseBoundaryFaces();
886 else locTool->updateForUse();
888 if ( doExtrapolationAtStart ) locTool->setExtrapolation(
false);
891 uint16_type nMPIsearch=15;
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();
895 if (!this->interpolationType().searchWithCommunication()) nMPIsearch=1;
896 uint16_type counterMPIsearch=1;
897 bool FinishMPIsearch=
false;
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;
906 while(!FinishMPIsearch)
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>();
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;
921 auto memory_localisationFail = this->updateNoRelationMeshMPI_upWithMyWorld( memmapGdof,
926 memory_valueInMatrix,
927 memory_col_globalProcessToGlobalCluster,
930 dof_extrapolationData
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;
940 if (this->interpolationType().searchWithCommunication())
942 auto memory_localisationFail2 = this->updateNoRelationMeshMPI_upWithOtherWorld( memmapGdof,
947 memory_valueInMatrix,
948 memory_col_globalProcessToGlobalCluster,
951 dof_extrapolationData
954 const size_type nbLocalisationFail_loc = memory_localisationFail.size()+memory_localisationFail2.size();
955 mpi::all_reduce( this->worldCommFusion().globalComm(),
956 nbLocalisationFail_loc,
958 std::plus<double>() );
960 else nbLocalisationFail = memory_localisationFail.size();
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;
968 if ( this->worldCommFusion().globalRank() == this->dualImageSpace()->worldComm().masterRank() )
969 std::cout <<
" it " << counterMPIsearch <<
" nbLocalisationFail " << nbLocalisationFail << std::endl;
973 if (counterMPIsearch<nMPIsearch && nbLocalisationFail>0) ++counterMPIsearch;
974 else FinishMPIsearch=
true;
978 if ( doExtrapolationAtStart ) locTool->setExtrapolation(
true);
980 if ( doExtrapolationAtStart && nbLocalisationFail>0 )
983 std::vector<std::set<size_type> > dof_searchWithProcExtrap(this->dualImageSpace()->nLocalDof());
985 uint16_type nMPIsearchExtrap=5;
986 if (this->domainSpace()->mesh()->worldComm().localSize()<5) nMPIsearchExtrap=this->domainSpace()->mesh()->worldComm().localSize();
988 if (!this->interpolationType().searchWithCommunication()) nMPIsearchExtrap=1;
989 uint16_type counterMPIsearchExtrap=1;
990 bool FinishMPIsearchExtrap=
false;
992 while(!FinishMPIsearchExtrap)
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>();
1000 auto memory_localisationFail = this->updateNoRelationMeshMPI_upWithMyWorld( memmapGdof,
1005 memory_valueInMatrix,
1006 memory_col_globalProcessToGlobalCluster,
1007 dof_searchWithProcExtrap,
1009 dof_extrapolationData
1012 if (this->interpolationType().searchWithCommunication())
1014 auto memory_localisationFail2 = this->updateNoRelationMeshMPI_upWithOtherWorld( memmapGdof,
1019 memory_valueInMatrix,
1020 memory_col_globalProcessToGlobalCluster,
1021 dof_searchWithProcExtrap,
1023 dof_extrapolationData
1027 if (counterMPIsearchExtrap<nMPIsearchExtrap) ++counterMPIsearchExtrap;
1028 else FinishMPIsearchExtrap=
true;
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);
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 )
1042 auto const gdof = cpt_gdof;
1043 auto const pointExtrap = imagedof->dofPoint(gdof).template get<0>();
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)
1052 auto const procCurrent = it_proc->template get<0>();
1053 auto const bary = it_proc->template get<1>();
1055 double normDist = 0;
1056 for (
int q=0;q<image_mesh_type::nRealDim;++q) normDist+=std::pow(bary(q)-pointExtrap(q),2);
1058 if (std::sqrt(normDist)<distMin) { distMin=std::sqrt(normDist);procExtrapoled=procCurrent;}
1060 it_proc = it_extrap->begin();
1061 for ( ; it_proc!=en_proc; ++it_proc)
1063 if ( it_proc->template get<0>()==procExtrapoled)
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);
1072 ublas::column( ptsRef, 0 ) = it_proc->template get<2>();
1074 MlocEval = domainbasis->evaluate( ptsRef );
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 )
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);
1088 value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1089 comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1092 memory_valueInMatrix[gdof].push_back(boost::make_tuple(procExtrapoled,j_gdof,v));
1093 memory_col_globalProcessToGlobalCluster[procExtrapoled][j_gdof]=j_gdofGlobalCluster;
1103 this->worldCommFusion().barrier();
1107 sparsity_graph->close();
1110 this->worldCommFusion().barrier();
1113 for (
int p=0;p<nProc_col;++p)
1115 mapCol_nLocalDof += memory_col_globalProcessToGlobalCluster[p].size();
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);
1122 for (
int p = 0 ; p<nProc_col;++p)
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)
1128 mapCol_globalProcessToGlobalCluster[currentLocalDof]=it_map->second;
1129 mapCol_LocalSpaceDofToLocalInterpDof[p][it_map->first]=currentLocalDof;
1131 if ( new_firstdofcol<=it_map->second && new_lastdofcol>=it_map->second)
1132 new_mapGlobalClusterToGlobalProcess[it_map->second-new_firstdofcol]=currentLocalDof;
1140 this->worldCommFusion().barrier();
1145 boost::shared_ptr<DataMap> mapColInterp(
new DataMap(this->dualImageSpace()->worldComm()));
1146 mapColInterp->setNDof(this->domainSpace()->mapOnOff().nDof());
1148 mapColInterp->setNLocalDofWithoutGhost( proc_id, new_nLocalDofWithoutGhost );
1149 mapColInterp->setNLocalDofWithGhost( proc_id, mapCol_nLocalDof );
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);
1160 this->worldCommFusion().barrier();
1163 if ( this->dualImageSpace()->worldComm().isActive() )
1165 this->matPtr() = this->backend()->newMatrix( mapColInterp,
1166 this->dualImageSpace()->dofOn(),
1171 this->worldCommFusion().barrier();
1174 if ( !this->dualImageSpace()->worldComm().isActive() && buildNonZeroMatrix )
1176 this->matPtr() = this->backend()->newZeroMatrix( mapColInterp,
1177 this->dualImageSpace()->dofOn() );
1181 this->worldCommFusion().barrier();
1184 if ( this->dualImageSpace()->worldComm().isActive() )
1186 for (
size_type idx_i=0 ; idx_i<nrow_dof_on_proc;++idx_i)
1188 for (
auto it_j=memory_valueInMatrix[idx_i].begin(),en_j=memory_valueInMatrix[idx_i].end() ; it_j!=en_j ; ++it_j)
1190 if (memory_valueInMatrix[idx_i].size()>0)
1192 this->matPtr()->set(idx_i,mapCol_LocalSpaceDofToLocalInterpDof[it_j->template get<0>()][it_j->template get<1>()],it_j->template get<2>());
1199 this->worldCommFusion().barrier();
1207 template<
typename DomainSpaceType,
typename ImageSpaceType,
typename IteratorRange,
typename InterpType>
1208 std::list<boost::tuple<size_type,uint16_type> >
1209 OperatorInterpolation<DomainSpaceType,
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 )
1223 std::list<boost::tuple<size_type,uint16_type> > memory_localisationFail;
1225 const auto proc_id = this->domainSpace()->mesh()->worldComm().localRank();
1227 auto const* imagedof = this->dualImageSpace()->dof().get();
1228 auto const* domaindof = this->domainSpace()->dof().get();
1229 auto const* domainbasis = this->domainSpace()->basis().get();
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);
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;
1241 size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1242 auto const& eltRandom = this->domainSpace()->mesh()->element(eltIdLocalised);
1244 for (
size_type k=0 ; k<memmapGdof[proc_id].size() ; ++k)
1247 if (this->interpolationType().componentsAreSamePoint())
1251 ublas::column(ptsReal,0 ) = pointsSearched[proc_id][k];
1253 if (InterpType::value==1)
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];
1261 verticesOfEltSearched = eltRandom.vertices();
1265 auto resLocalisation = locTool->run_analysis(ptsReal,eltIdLocalised,verticesOfEltSearched,
1266 mpl::int_<interpolation_type::value>());
1267 if (!resLocalisation.template get<0>()[0])
1269 for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
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);
1278 eltIdLocalised = resLocalisation.template get<1>();
1280 if (extrapolation_mode)
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;
1286 for (
int qj=0;qj<verticesExtrapoled.size2();++qj)
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);
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 )
1299 const auto gdof = memmapGdof[proc_id][k+comp];
1300 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
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]);
1305 dof_extrapolationData[gdof].push_back(boost::make_tuple(proc_id,bary,theRefPtExtrap,j_gdofs,comp));
1310 for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1312 const auto gdof = memmapGdof[proc_id][k+comp];
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);
1320 auto itanal = locTool->result_analysis().begin();
1321 auto const itanal_end = locTool->result_analysis().end();
1322 for ( ;itanal!=itanal_end;++itanal)
1324 const auto itL=itanal->second.begin();
1325 ublas::column( ptsRef, 0 ) = boost::get<1>(*itL);
1327 MlocEval = domainbasis->evaluate( ptsRef );
1328 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1331 const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( itanal->first,jloc,comp ));
1333 row.template get<2>().insert(domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1335 const value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1336 comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
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];
1344 dof_searchWithProc[gdof].insert(proc_id);
1349 k=k+image_basis_type::nComponents-1;
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>();
1359 if (InterpType::value==1)
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];
1367 verticesOfEltSearched = eltRandom.vertices();
1371 auto resLocalisation = locTool->run_analysis(ptsReal, eltIdLocalised, verticesOfEltSearched, mpl::int_<interpolation_type::value>());
1372 if (!resLocalisation.template get<0>()[0])
1374 memory_localisationFail.push_back(boost::make_tuple(gdof,comp) );
1375 dof_searchWithProc[gdof].insert(proc_id);
1379 eltIdLocalised = resLocalisation.template get<1>();
1380 if (extrapolation_mode)
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;
1386 for (
int qj=0;qj<verticesExtrapoled.size2();++qj)
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);
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 )
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]);
1402 dof_extrapolationData[gdof].push_back(boost::make_tuple(proc_id,bary,theRefPtExtrap,j_gdofs,comp));
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);
1413 auto itanal = locTool->result_analysis().begin();
1414 auto const itanal_end = locTool->result_analysis().end();
1415 for ( ;itanal!=itanal_end;++itanal)
1417 const auto itL=itanal->second.begin();
1418 ublas::column( ptsRef, 0 ) = boost::get<1>(*itL);
1420 MlocEval = domainbasis->evaluate( ptsRef );
1421 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1424 const size_type j_gdof = boost::get<0>(domaindof->localToGlobal( itanal->first,jloc,comp ));
1426 row.template get<2>().insert(domaindof->mapGlobalProcessToGlobalCluster()[j_gdof]);
1428 const value_type v = MlocEval( domain_basis_type::nComponents1*jloc +
1429 comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
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];
1437 dof_searchWithProc[gdof].insert(proc_id);
1443 return memory_localisationFail;
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)
1462 std::list<boost::tuple<size_type,uint16_type> > memory_localisationFail;
1464 auto const* imagedof = this->dualImageSpace()->dof().get();
1465 auto const* domaindof = this->domainSpace()->dof().get();
1466 auto const* domainbasis = this->domainSpace()->basis().get();
1468 const size_type proc_id = this->dualImageSpace()->worldComm().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();
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);
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;
1489 size_type eltIdLocalised = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank())->id();
1490 auto const& eltRandom = this->domainSpace()->mesh()->element(eltIdLocalised);
1492 std::vector<bool> dof_done( this->dualImageSpace()->nLocalDof(),
false);
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);
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);
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 );
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 );
1527 int firstActiveProc_image=0;
1528 bool findFirstActive_image=
false;
1529 while (!findFirstActive_image)
1531 if (imageProcIsActive_fusion[firstActiveProc_image])
1533 findFirstActive_image=
true;
1535 else ++firstActiveProc_image;
1537 int firstActiveProc_domain=0;
1538 bool findFirstActive_domain=
false;
1539 while (!findFirstActive_domain)
1541 if (domainProcIsActive_fusion[firstActiveProc_domain])
1543 findFirstActive_domain=
true;
1545 else ++firstActiveProc_domain;
1548 for (
int p=0;p<localMeshRankToWorldCommFusion_image.size(); ++p)
1550 if (!this->dualImageSpace()->worldComm().isActive()) localMeshRankToWorldCommFusion_image[p]=p%nProc_image+firstActiveProc_image;
1552 for (
int p=0;p<localMeshRankToWorldCommFusion_domain.size(); ++p)
1554 if (!this->domainSpace()->worldComm().isActive()) localMeshRankToWorldCommFusion_domain[p]=p%nProc_domain+firstActiveProc_domain;
1558 std::vector<std::vector<int> > searchDistribution(nProc);
1559 for (
int p=0;p<nProc_image;++p)
1561 searchDistribution[p].clear();
1562 if ( proc_id_image == p && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1564 for (
int q=0;q<nProc_domain;++q)
1566 if( pointsSearched[q].size()>0 && localMeshRankToWorldCommFusion_image[p]!=localMeshRankToWorldCommFusion_domain[q] )
1568 searchDistribution[p].push_back(q);
1573 mpi::broadcast( this->worldCommFusion().globalComm(), searchDistribution[p], localMeshRankToWorldCommFusion_image[p] );
1579 std::vector<std::list<int> > searchDistribution(nProc);
1580 for (
int p=0;p<nProc_image;++p)
1582 searchDistribution[p].clear();
1583 for (
int q=0;q<nProc_domain;++q)
1586 if( (localMeshRankToWorldCommFusion_image[p])!=localMeshRankToWorldCommFusion_domain[q] )
1588 searchDistribution[p].push_back(q);
1595 this->worldCommFusion().barrier();
1596 for (
int p=0;p<this->worldCommFusion().globalSize();++p)
1598 if (p==this->worldCommFusion().globalRank())
1600 std::cout <<
"I am proc " << p <<
" ";
1601 for (
int q=0;q<nProc_image;++q)
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)
1608 std::cout << *it_list <<
" ";
1610 std::cout << std::endl;
1613 this->worldCommFusion().barrier();
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;
1625 for (
int proc=0;proc<nProc_row;++proc)
1627 if ( proc_id_image == proc && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1629 for (
auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1630 it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
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]);
1642 for (
auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1643 it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
1645 const int rankLocalization = *it_rankLocalization;
1646 if ( proc_id_domain == rankLocalization && domainProcIsActive_fusion[this->worldCommFusion().globalRank()] )
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);
1654 const size_type nDataRecv = dataToRecv.size();
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));
1666 ublas::column(ptsReal,0) = dataToRecv[k];
1668 if (InterpType::value==1)
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];
1676 verticesOfEltSearched = eltRandom.vertices();
1679 auto resLocalisation = locTool->run_analysis(ptsReal,eltIdLocalised,verticesOfEltSearched,mpl::int_<interpolation_type::value>());
1680 if (resLocalisation.template get<0>()[0])
1682 eltIdLocalised = resLocalisation.template get<1>();
1684 if (!this->interpolationType().componentsAreSamePoint() )
1687 const uint16_type comp=dataToRecv_Comp[k];
1688 pointsRefIsFinded[k]=
true;
1689 pointsIdEltFinded[k]=eltIdLocalised;
1691 auto itanal = locTool->result_analysis_begin();
1692 auto const itanal_end = locTool->result_analysis_end();
1693 for ( ;itanal!=itanal_end;++itanal)
1695 auto const itL=itanal->second.begin();
1697 pointsRefFinded[k] = boost::get<1>(*itL);
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 )
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];
1708 if (extrapolation_mode)
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;
1714 for (
int qj=0;qj<verticesExtrapoled.size2();++qj)
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);
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;
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 )
1733 pointsRefIsFinded[k+comp]=
true;
1734 pointsIdEltFinded[k+comp]=eltIdLocalised;
1736 pointsRefFinded[k+comp] = ptRefFinded;
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 )
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];
1748 if (extrapolation_mode)
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;
1754 for (
int qj=0;qj<verticesExtrapoled.size2();++qj)
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);
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;
1767 k=k+image_basis_type::nComponents-1;
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);
1796 if ( proc_id_image == proc && imageProcIsActive_fusion[this->worldCommFusion().globalRank()] )
1798 for (
auto it_rankLocalization=searchDistribution[proc].begin(),en_rankLocalization=searchDistribution[proc].end();
1799 it_rankLocalization!=en_rankLocalization;++it_rankLocalization)
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);
1811 const int rankRecv=rankLocalization;
1812 for (
int k=0;k<pointsRefFinded.size();++k)
1814 if (!pointsRefIsFinded[k])
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);
1822 const auto i_gdof = memmapGdof[rankLocalization][k];
1823 if (!dof_done[i_gdof])
1826 if (extrapolation_mode)
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 )
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);
1838 dof_extrapolationData[i_gdof].push_back(boost::make_tuple(rankLocalization,bary,theRefPtExtrap,j_gdofs,comp));
1843 ublas::column( ptsRef, 0 ) = pointsRefFinded[k];
1845 MlocEval = domainbasis->evaluate( ptsRef );
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);
1855 for ( uint16_type jloc = 0; jloc < domain_basis_type::nLocalDof; ++jloc )
1858 const size_type j_gdof = pointsDofsColFinded[k][jloc];
1860 const size_type j_gdof_gc = pointsDofsGlobalClusterColFinded[k][jloc];
1862 row.template get<2>().insert(j_gdof_gc);
1864 const auto v = MlocEval( domain_basis_type::nComponents1*jloc +
1865 comp*domain_basis_type::nComponents1*domain_basis_type::nLocalDof +
1869 memory_valueInMatrix[i_gdof].push_back(boost::make_tuple(rankRecv,j_gdof,v));
1871 memory_col_globalProcessToGlobalCluster[rankRecv][j_gdof]=j_gdof_gc;
1875 dof_searchWithProc[i_gdof].insert(rankRecv);
1877 dof_done[i_gdof]=
true;
1888 return memory_localisationFail;
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)
1903 const size_type nProc_domain = this->domainSpace()->mesh()->worldComm().localSize();
1905 auto const* imagedof = this->dualImageSpace()->dof().get();
1906 iterator_type it, en;
1908 auto locTool = this->domainSpace()->mesh()->tool_localization();
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 );
1915 std::vector<typename image_mesh_type::node_type> vecBarycenter(nProc_domain);
1916 mpi::all_gather( this->domainSpace()->mesh()->worldComm().localComm(),
1917 locTool->barycenter(),
1925 double distanceMin=0,distance=0,distanceSquare=0;
1928 if ( this->dualImageSpace()->worldComm().isActive() )
1931 auto itListRange = M_listRange.begin();
1932 auto const enListRange = M_listRange.end();
1933 for ( ; itListRange!=enListRange ; ++itListRange)
1935 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
1936 for ( ; it!=en;++it )
1938 for ( uint16_type iloc = 0; iloc < nLocalDofInDualImageElt; ++iloc )
1940 for ( uint16_type comp = 0;comp < image_basis_type::nComponents;++comp )
1942 const auto gdof = boost::get<0>(imagedof->localToGlobal( *it, iloc, comp ));
1943 if (!dof_done[gdof] && memory_valueInMatrix[gdof].size()==0)
1946 const auto imagePoint = imagedof->dofPoint(gdof).template get<0>();
1948 if (this->interpolationType().searchWithCommunication())
1950 distanceMin=INT_MAX;
1951 for (
int proc=0 ; proc<nProc_domain; ++proc)
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() )
1961 distanceMin=distance;
1964 memSetGdofAndComp[procForPt].push_back(boost::make_tuple(gdof,comp));
1965 if (InterpType::value==1)
1966 memSetVertices_conformeInterp[procForPt].push_back(it->vertices());
1970 memSetGdofAndComp[this->domainSpace()->worldComm().globalRank()].push_back(boost::make_tuple(gdof,comp));
1971 if (InterpType::value==1)
1973 memSetVertices_conformeInterp[this->domainSpace()->worldComm().globalRank()].push_back(it->vertices());
1977 dof_done[gdof]=
true;
1986 std::vector< std::vector<size_type> > memmapGdof( nProc_domain );
1988 std::vector< std::vector<uint16_type> > memmapComp( nProc_domain );
1990 std::vector<std::vector<typename image_mesh_type::node_type> > pointsSearched( nProc_domain );
1992 std::vector<std::vector< std::vector<typename image_mesh_type::node_type > > > memmap_vertices( nProc_domain );
1994 for (
int proc=0; proc<nProc_domain;++proc)
1996 const size_type nData = memSetGdofAndComp[proc].size();
1997 memmapGdof[proc].resize(nData);
1998 memmapComp[proc].resize(nData);
1999 pointsSearched[proc].resize(nData);
2001 if(InterpType::value==1) memmap_vertices[proc].resize(nData);
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)
2007 memmapGdof[proc][k]=it_GdofAndComp->template get<0>();
2008 memmapComp[proc][k]=it_GdofAndComp->template get<1>();
2009 pointsSearched[proc][k]=imagedof->dofPoint(it_GdofAndComp->template get<0>()).template get<0>();
2010 if(InterpType::value==1)
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);
2020 return boost::make_tuple(memmapGdof,memmapComp,pointsSearched,memmap_vertices);
2032 template <
typename RangeType>
2033 struct opinterprangetype
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;
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,
2051 typedef OperatorInterpolation<DomainSpaceType, ImageSpaceType,typename Feel::detail::opinterprangetype<IteratorRange>::type,InterpType> operatorinterpolation_type;
2053 boost::shared_ptr<operatorinterpolation_type> opI(
new operatorinterpolation_type( domainspace,imagespace,r,backend,interptype,ddmethod ) );
2060 template<
typename Args>
2061 struct compute_opInterpolation_return
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;
2066 typedef typename boost::remove_const<
2067 typename boost::remove_reference<
2068 typename parameter::binding<Args,
2070 typename OperatorInterpolation<domain_space_type, image_space_type>::range_iterator
2071 >::type >::type >::type iterator_base_range_type;
2073 typedef typename Feel::detail::opinterprangetype<iterator_base_range_type>::type iterator_range_type;
2075 typedef typename boost::remove_const<
2076 typename boost::remove_reference<
2077 typename parameter::binding<Args,
2079 InterpolationNonConforme
2080 >::type >::type >::type interpolation_type;
2083 typedef boost::shared_ptr<OperatorInterpolation<domain_space_type, image_space_type,iterator_range_type,interpolation_type> > type;
2086 BOOST_PARAMETER_FUNCTION(
2087 (
typename compute_opInterpolation_return<Args>::type ),
2091 ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
2092 ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
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 )
2102 Feel::detail::ignore_unused_variable_warning( args );
2104 return opInterp( domainSpace,imageSpace,range,backend,type,ddmethod );