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
doftablempi.hpp
Go to the documentation of this file.
1 /*
2  This file is part of the Feel library
3 
4  Copyright (C) 2012 Université de Grenoble 1
5 
6  This library is free software; you can redistribute it and/or
7  modify it under the terms of the GNU Lesser General Public
8  License as published by the Free Software Foundation; either
9  version 3.0 of the License, or (at your option) any later version.
10 
11  This library is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public
17  License along with this library; if not, write to the Free Software
18  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
24 #if !defined(FEELPP_DOFTABLE_MPI_HPP)
25 #define FEELPP_DOFTABLE_MPI_HPP 1
26 
27 namespace Feel
28 {
29 
30 template<typename MeshType, typename FEType, typename PeriodicityType>
31 void
32 DofTable<MeshType, FEType, PeriodicityType>::buildGhostDofMap( mesh_type& mesh )
33 {
34 #if 0 //old
35 
36  // dofs map on interprocess faces : ( dofOnMyRank, tuple< otherProcRank, dofOnOtherRank> )
37  std::map<size_type,boost::tuple<size_type,size_type> > mapInterProcessDof;
38  // dofs no present on interprocess faces : ( dofOnMyRank, tuple< otherProcRank, dofOnOtherRank> )
39  std::map<size_type,boost::tuple<size_type,size_type> > setInterProcessDofNotPresent;
40 
41  if (is_continuous)
42  {
43  DVLOG(2) << "[buildGhostDofMap] call buildGhostInterProcessDofMap() with god rank "<< this->worldComm().godRank() << "\n";
44  buildGhostInterProcessDofMap( mesh,mapInterProcessDof );
45 
46  DVLOG(2) << "[buildGhostDofMap] call buildDofNotPresent() with rank "<< this->worldComm().rank() << "\n";
47  buildDofNotPresent( mapInterProcessDof,setInterProcessDofNotPresent );
48  }
49 
50  DVLOG(2) << "[buildGhostDofMap] call buildGlobalProcessToGlobalClusterDofMap() with rank "<< this->worldComm().rank() << "\n";
51  buildGlobalProcessToGlobalClusterDofMap( mesh,setInterProcessDofNotPresent );
52 
53 #else
54 
55  if (is_continuous)
56  {
57  DVLOG(2) << "[buildGhostDofMap] call buildGlobalProcessToGlobalClusterDofMapContinuous() with god rank "<< this->worldComm().godRank() << "\n";
58  buildGlobalProcessToGlobalClusterDofMapContinuous( mesh );
59  }
60  else
61  {
62  DVLOG(2) << "[buildGhostDofMap] call buildGlobalProcessToGlobalClusterDofMapDiscontinuous() with rank "<< this->worldComm().rank() << "\n";
63  buildGlobalProcessToGlobalClusterDofMapDiscontinuous();
64  }
65 
66 #endif
67 
68 
69  DVLOG(2) << "[buildGhostDofMap] call localtoglobalOnCluster() with rank "<< this->worldComm().rank() << "\n";
70  auto it_elt = mesh.beginElementWithProcessId( this->comm().rank() );
71  auto en_elt = mesh.endElementWithProcessId( this->comm().rank() );
72 
73  for ( ; it_elt != en_elt; ++it_elt )
74  {
75  size_type elid= it_elt->id();
76 
77  for ( int i = 0; i < FEType::nLocalDof; ++i )
78  {
79  int nc1 = ( is_product?nComponents1:1 );
80 
81  for ( int c1 =0; c1 < nc1; ++c1 )
82  {
83  int ind = FEType::nLocalDof*c1+i;
84  auto const& dof = localToGlobalOnCluster( elid, i, c1 );
85 
86  M_locglobOnCluster_indices[elid][ind] = dof.index();
87  M_locglobOnCluster_signs[elid][ind] = dof.sign();
88  }
89  }
90  }
91 
92  DVLOG(2) << "[buildGhostDofMap] finish () with god rank "<< this->worldComm().godRank() << "\n";
93 
94 }
95 
96 
97 
98 //--------------------------------------------------------------------------------------------------------//
99 //--------------------------------------------------------------------------------------------------------//
100 //--------------------------------------------------------------------------------------------------------//
101 //--------------------------------------------------------------------------------------------------------//
102 
103 template<typename MeshType, typename FEType, typename PeriodicityType>
104 void
105 DofTable<MeshType, FEType, PeriodicityType>::buildGlobalProcessToGlobalClusterDofMapContinuous( mesh_type& mesh )
106 {
107  //------------------------------------------------------------------------------//
109  if ( !fe_type::is_modal )
110  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
111  face_type::numEdges * fe_type::nDofPerEdge +
112  face_type::numFaces * fe_type::nDofPerFace );
113  else
114  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
115 
116  if ( nbFaceDof == 0 ) return;
117  //------------------------------------------------------------------------------//
118  // build GlobalProcessToGlobalClusterDofMap for actif dofs and prepare send for ghost dof
119  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > listToSend(this->worldComm().size());
120  this->buildGlobalProcessToGlobalClusterDofMapContinuousActifDof(mesh,listToSend);
121  //------------------------------------------------------------------------------//
122  // update GlobalProcessToGlobalClusterDofMap for ghost dofs
123  this->buildGlobalProcessToGlobalClusterDofMapContinuousGhostDof(mesh,listToSend);
124  //------------------------------------------------------------------------------//
125 }
126 
127 //--------------------------------------------------------------------------------------------------------//
128 //--------------------------------------------------------------------------------------------------------//
129 //--------------------------------------------------------------------------------------------------------//
130 //--------------------------------------------------------------------------------------------------------//
131 
132 namespace detail {
133 
134 template <typename MeshType>
135 boost::tuple<int,size_type>
136 updateDofOnVertices( MeshType const& mesh, typename MeshType::face_type const& theface,
137  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
138  /*typename MeshType::node_type const& thedofnode,*/ const uint16_type locDof )
139 {
140  int procMin = IdProcessOfGhost;
141  size_type idFaceMin = idFaceInPartition;
142 
143  uint16_type iFaEl = ( theface.processId() == theface.proc_first() )? theface.pos_first():theface.pos_second();
144  //local point number (in element)
145  uint16_type iPtEl = MeshType::element_type::fToP( iFaEl, locDof );
146 
147 #if 0
148  //typedef MeshType::face_type
149  bool hasFind=false;
150  for ( uint16_type p=0;p<MeshType::face_type::numVertices && !hasFind; ++p)
151  {
152  auto const& thept = theface.point(p);
153  bool find=true;
154  for (uint16_type d=0;d<MeshType::nRealDim;++d)
155  {
156  find = find && ( std::abs( thept(d)-thedofnode[d] )<1e-9 );
157  }
158  if (find)
159  {
160  hasFind=true;
161  auto const theptId = thept.id();
162  if (eltOnProc.point(iPtEl).id()!=theptId) std::cout << "COOOOYIOKLK"<< std::endl;
163 
164  auto iteltghost=thept.elementsGhost().begin();
165  auto const eneltghost=thept.elementsGhost().end();
166  //size_type idFaceDDD=0;
167  for ( ; iteltghost!=eneltghost ; ++iteltghost)
168  {
169  if (procMin>iteltghost->template get<0>())
170  {
171  procMin=iteltghost->template get<0>();
172  auto const& eltGhost = mesh.element(iteltghost->template get<1>(),procMin);
173  bool findFace=false;
174  for ( uint16_type f = 0; f < MeshType::element_type::numTopologicalFaces && !findFace; ++f )
175  {
176  auto const& faceOnGhost = eltGhost.face(f);
177  //auto const& verticefaceOnGhost =faceOnGhost.vertices();
178  for ( uint16_type vv = 0; vv < MeshType::face_type::numVertices && !findFace ; ++vv )
179  {
180  if (faceOnGhost.point(vv).id()==theptId) {findFace=true;idFaceMin = faceOnGhost.idInOthersPartitions(procMin); }
181  }
182  }
183  if (!findFace) std::cout << "YYYYYYY" << std::endl;
184 
185  }
186  }
187  }
188  }
189  if (!hasFind) std::cout << "JKJKJKJKJKJKJ" << std::endl;
190 #else
191  auto const& thept = eltOnProc.point(iPtEl);
192  auto const theptId = thept.id();
193 
194  auto iteltghost=thept.elementsGhost().begin();
195  auto const eneltghost=thept.elementsGhost().end();
196  for ( ; iteltghost!=eneltghost ; ++iteltghost)
197  {
198  if (procMin>iteltghost->template get<0>())
199  {
200  procMin=iteltghost->template get<0>();
201  auto const& eltGhost = mesh.element(iteltghost->template get<1>(),procMin);
202  bool findFace=false;
203  for ( uint16_type f = 0; f < MeshType::element_type::numTopologicalFaces && !findFace; ++f )
204  {
205  auto const& faceOnGhost = eltGhost.face(f);
206  for ( uint16_type vv = 0; vv < MeshType::face_type::numVertices && !findFace ; ++vv )
207  {
208  if (faceOnGhost.point(vv).id()==theptId) {findFace=true;idFaceMin = faceOnGhost.idInOthersPartitions(procMin); }
209  }
210  }
211  CHECK( findFace ) << "\nPROBLEM with parallel dof table construction \n";
212  }
213  }
214 #endif
215  return boost::make_tuple(procMin,idFaceMin);
216 
217 }
218 
219 //--------------------------------------------------------------------------------------------------------//
220 
221 template <typename MeshType>
222 boost::tuple<int,size_type>
223 updateDofOnEdges( MeshType const& mesh, typename MeshType::face_type const& theface,
224  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
225  const uint16_type idEdgesInFace, mpl::int_<0> )
226 {
227  return boost::make_tuple(0,0);
228 }
229 template <typename MeshType>
230 boost::tuple<int,size_type>
231 updateDofOnEdges( MeshType const& mesh, typename MeshType::face_type const& theface,
232  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
233  const uint16_type idEdgesInFace, mpl::int_<1> )
234 {
235  return boost::make_tuple(0,0);
236 }
237 template <typename MeshType>
238 boost::tuple<int,size_type>
239 updateDofOnEdges( MeshType const& mesh, typename MeshType::face_type const& theface,
240  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
241  const uint16_type idEdgesInFace, mpl::int_<2> )
242 {
243  return boost::make_tuple(0,0);
244 }
245 template <typename MeshType>
246 boost::tuple<int,size_type>
247 updateDofOnEdges( MeshType const& mesh, typename MeshType::face_type const& theface,
248  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
249  const uint16_type idEdgesInFace, mpl::int_<3> )
250 {
251  int procMin = IdProcessOfGhost;
252  size_type idFaceMin = idFaceInPartition;
253 
254  uint16_type iFaEl = ( theface.processId() == theface.proc_first() )? theface.pos_first():theface.pos_second();
255  //local edge number (in element)
256  uint16_type iEdEl = MeshType::element_type::fToE( iFaEl, idEdgesInFace );
257 
258 
259  auto const& theedge = eltOnProc.edge(iEdEl);
260  auto const theedgeId = theedge.id();
261  auto iteltghost=theedge.elementsGhost().begin();
262  auto const eneltghost=theedge.elementsGhost().end();
263 
264  for ( ; iteltghost!=eneltghost ; ++iteltghost)
265  {
266  if (procMin>iteltghost->template get<0>())
267  {
268  procMin=iteltghost->template get<0>();
269  auto const& eltGhost = mesh.element(iteltghost->template get<1>(),procMin);
270  bool findFace=false;
271  for ( uint16_type f = 0; f < MeshType::element_type::numTopologicalFaces && !findFace; ++f )
272  {
273  auto const& faceOnGhost = eltGhost.face(f);
274  for ( uint16_type vv = 0; vv < MeshType::face_type::numEdges && !findFace ; ++vv )
275  {
276  if (faceOnGhost.edge(vv).id()==theedgeId) { findFace=true;idFaceMin = faceOnGhost.idInOthersPartitions(procMin); }
277  }
278  }
279  CHECK( findFace ) << "\nPROBLEM with parallel dof table construction \n";
280  }
281  }
282 
283  return boost::make_tuple(procMin,idFaceMin);
284 }
285 
286 //--------------------------------------------------------------------------------------------------------//
287 
288 template <typename MeshType>
289 boost::tuple<int,size_type>
290 updateDofOnEdges( MeshType const& mesh, typename MeshType::face_type const& theface,
291  const int IdProcessOfGhost, const size_type idFaceInPartition, typename MeshType::element_type const& eltOnProc,
292  const uint16_type idEdgesInFace )
293 {
294  return updateDofOnEdges(mesh,theface,IdProcessOfGhost,idFaceInPartition,eltOnProc,idEdgesInFace,mpl::int_<MeshType::nDim>() );
295 }
296 
297 } // namespace detail
298 
299 //--------------------------------------------------------------------------------------------------------//
300 //--------------------------------------------------------------------------------------------------------//
301 //--------------------------------------------------------------------------------------------------------//
302 //--------------------------------------------------------------------------------------------------------//
303 
304 template<typename MeshType, typename FEType, typename PeriodicityType>
305 void
306 DofTable<MeshType, FEType, PeriodicityType>::buildGlobalProcessToGlobalClusterDofMapContinuousActifDof( mesh_type& mesh,
307  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > & listToSend )
308 {
309  // goal init container listToSend
310  // std::vector< std::map<size_type,std::set<size_type> > > ( proc,( idFace,(globDof,..)), ... )) )
311  const int myRank = this->worldComm().rank();
312  //------------------------------------------------------------------------------//
313  // get nbFaceDof
315  if ( !fe_type::is_modal )
316  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
317  face_type::numEdges * fe_type::nDofPerEdge +
318  face_type::numFaces * fe_type::nDofPerFace );
319  else
320  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
321 
322  DVLOG(2) << "[buildGhostInterProcessDofMap] nbFaceDof " << nbFaceDof << "\n";
323 
324  if ( nbFaceDof == 0 ) return;
325 
326  const uint16_type ncdof = is_product?nComponents:1;
327  DVLOG(2) << "[buildGhostInterProcessDofMap] ncdof " << ncdof << "\n";
328 
329  //------------------------------------------------------------------------------//
330  std::vector<bool> dofdone(this->M_n_localWithGhost_df[myRank],false);
331  size_type nDofNotPresent=0;
332 
333  // iteration on all interprocessfaces in order to send requests to the near proc
334  auto face_it = mesh.interProcessFaces().first;
335  auto const face_en = mesh.interProcessFaces().second;
336  DVLOG(2) << "[buildGhostInterProcessDofMap] nb interprocess faces: " << std::distance( face_it, face_en ) << "\n";
337  for ( ; face_it != face_en ; ++face_it )
338  {
339  DVLOG(2) << "[buildGhostInterProcessDofMap] face id: " << face_it->id() << "\n";
340  auto const& elt0 = face_it->element0();
341  auto const& elt1 = face_it->element1();
342  const bool elt0isGhost = elt0.isGhostCell();
343  auto const& eltOnProc = (elt0isGhost)?elt1:elt0;
344  auto const& eltOffProc = (elt0isGhost)?elt0:elt1;
345  DVLOG(2) << "[buildGhostInterProcessDofMap] (myRank:" << myRank << ") eltOnProc id: " << eltOnProc.id() << "G(): " << eltOnProc.G() << "\n";
346  DVLOG(2) << "[buildGhostInterProcessDofMap] (myRank:" << myRank << ") eltOffProc id: " << eltOffProc.id() << "G(): " << eltOffProc.G() << "\n";
347 
348  //------------------------------------------------------------------------------//
349 
350  const int IdProcessOfGhostIP = eltOffProc.processId();
351  const size_type idFaceInPartitionIP = face_it->idInOthersPartitions( IdProcessOfGhostIP );
352  int IdProcessOfGhost = IdProcessOfGhostIP;
353  size_type idFaceInPartition = idFaceInPartitionIP;
354 
355  //------------------------------------------------------------------------------//
356  // for each dof in face
357  for ( uint16_type locDof = 0; locDof < nbFaceDof; ++locDof )
358  {
359  IdProcessOfGhost = IdProcessOfGhostIP;
360  idFaceInPartition = idFaceInPartitionIP;
361  if ( locDof < face_type::numVertices*fe_type::nDofPerVertex)
362  {
363  const int nDofPerVertexTemp = mpl::if_<boost::is_same<mpl::int_<fe_type::nDofPerVertex>,mpl::int_<0> >,
364  mpl::int_<1>,
365  mpl::int_<fe_type::nDofPerVertex> >::type::value;
366  int pointGetLocDof = locDof / nDofPerVertexTemp;
367 
368  boost::tie( IdProcessOfGhost, idFaceInPartition ) = detail::updateDofOnVertices<mesh_type>(mesh,*face_it, IdProcessOfGhost, idFaceInPartition, eltOnProc, /*thedofnode,*/ pointGetLocDof/*locDof*/ );
369  }
370  else if ( nDim == 3 && locDof < (face_type::numVertices*fe_type::nDofPerVertex + face_type::numEdges*fe_type::nDofPerEdge) )
371  {
372  int locDofInEgde = locDof - face_type::numVertices*fe_type::nDofPerVertex;
373  const int nDofPerEdgeTemp = mpl::if_<boost::is_same<mpl::int_<fe_type::nDofPerEdge>,mpl::int_<0> >,
374  mpl::int_<1>,
375  mpl::int_<fe_type::nDofPerEdge> >::type::value;
376  int edgeGetLocDof = locDofInEgde / nDofPerEdgeTemp;
377 
378  boost::tie( IdProcessOfGhost, idFaceInPartition ) = detail::updateDofOnEdges<mesh_type>(mesh,*face_it, IdProcessOfGhost, idFaceInPartition, eltOnProc, edgeGetLocDof );
379  }
380 
381  if (IdProcessOfGhost<myRank)
382  {
383  for ( uint16_type c = 0; c < ncdof; ++c )
384  {
385  // add dof in subcontainer
386  const size_type theglobdof = faceLocalToGlobal( face_it->id(),locDof,c ).template get<0>();
387  if (!dofdone[theglobdof])
388  {
389  listToSend[IdProcessOfGhost][idFaceInPartition].insert(boost::make_tuple(theglobdof,c));
390  dofdone[theglobdof]=true;
391  ++nDofNotPresent;
392  }
393  }
394  }
395 
396  } // for ( uint16_type locDof = 0; locDof < nbFaceDof; ++locDof )
397 
398 
399  } // for ( ; face_it != face_en ; ++face_it )
400 
401  //------------------------------------------------------------------------------//
402  //------------------------------------------------------------------------------//
403  //------------------------------------------------------------------------------//
404  // update datamap info
405  //const size_type mynDofWithoutGhost = this->M_last_df[myRank] - this->M_first_df[myRank] + 1 - nDofNotPresent;
406  const size_type mynDofWithoutGhost = this->M_n_localWithGhost_df[myRank] - nDofNotPresent;
407  mpi::all_gather( this->worldComm(),
408  mynDofWithoutGhost,
409  this->M_n_localWithoutGhost_df );
410 
411  this->M_n_dofs=0;
412  for ( int proc=0; proc<this->worldComm().size(); ++proc )
413  {
414  this->M_n_dofs+=this->M_n_localWithoutGhost_df[proc];
415  }
416 
417  this->M_first_df_globalcluster=this->M_first_df;
418  this->M_last_df_globalcluster=this->M_last_df;
419  for ( int i=1; i<this->worldComm().size(); ++i )
420  {
421  if ( this->M_n_localWithoutGhost_df[i-1] >0 )
422  this->M_first_df_globalcluster[i]=this->M_last_df_globalcluster[i-1]+1;
423  else
424  this->M_first_df_globalcluster[i]=this->M_last_df_globalcluster[i-1];
425 
426  if ( this->M_n_localWithoutGhost_df[i] >0 )
427  this->M_last_df_globalcluster[i]=this->M_first_df_globalcluster[i]+this->M_n_localWithoutGhost_df[i]-1;
428  else
429  this->M_last_df_globalcluster[i]=this->M_first_df_globalcluster[i];
430  }
431  //------------------------------------------------------------------------------//
432  // init map
433  this->M_mapGlobalProcessToGlobalCluster.resize( this->M_n_localWithGhost_df[myRank],invalid_size_type_value );
434  this->M_mapGlobalClusterToGlobalProcess.resize( this->M_n_localWithoutGhost_df[myRank],invalid_size_type_value );
435  //------------------------------------------------------------------------------//
436  // add in map the dofs presents
437  size_type firstGlobIndex = this->M_first_df_globalcluster[myRank];
438  size_type nextGlobIndex = firstGlobIndex;
439  for ( size_type i=0; i< this->M_n_localWithGhost_df[myRank]; ++i )
440  {
441  //if ( setInterProcessDofNotPresent.find( i )==setInterProcessDofNotPresent.end() )
442  if (!dofdone[i])
443  {
444  this->M_mapGlobalProcessToGlobalCluster[i]=nextGlobIndex;
445  this->M_mapGlobalClusterToGlobalProcess[nextGlobIndex-firstGlobIndex]=i;
446  ++nextGlobIndex;
447  }
448  }
449  //------------------------------------------------------------------------------//
450 } // buildGlobalProcessToGlobalClusterDofMapContinuousActifDof
451 
452 //--------------------------------------------------------------------------------------------------------//
453 //--------------------------------------------------------------------------------------------------------//
454 //--------------------------------------------------------------------------------------------------------//
455 //--------------------------------------------------------------------------------------------------------//
456 
457 template<typename MeshType, typename FEType, typename PeriodicityType>
458 void
459 DofTable<MeshType, FEType, PeriodicityType>::buildGlobalProcessToGlobalClusterDofMapContinuousGhostDof( mesh_type& mesh,
460  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > const& listToSend )
461 {
462  const int myRank = this->worldComm().rank();
463  //--------------------------------------------------------------------------------------------------------//
465  if ( !fe_type::is_modal )
466  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
467  face_type::numEdges * fe_type::nDofPerEdge +
468  face_type::numFaces * fe_type::nDofPerFace );
469 
470  else
471  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
472  //--------------------------------------------------------------------------------------------------------//
473 
474  std::vector<int> nbMsgToSend( this->worldComm().size(), 0 );
475 
476  std::vector< std::vector< std::vector<size_type> > > memoryInitialRequest( this->worldComm().size() );
477 
478  typedef std::vector< boost::tuple<uint16_type, ublas::vector<double> > > dofs_in_face_subcontainer_type;
479  typedef boost::tuple<size_type, dofs_in_face_subcontainer_type > dofs_in_face_container_type;
480 
481  for ( int proc=0; proc<this->worldComm().size(); ++proc )
482  {
483  auto itFaces = listToSend[proc].begin();
484  auto const enFaces = listToSend[proc].end();
485  memoryInitialRequest[proc].resize(std::distance(itFaces,enFaces));
486  for ( int cptFaces=0 ; itFaces!=enFaces ; ++itFaces, ++cptFaces)
487  {
488  auto itDof = itFaces->second.begin();
489  auto const enDof = itFaces->second.end();
490  auto const nDofsInFace = std::distance(itDof,enDof);
491  dofs_in_face_subcontainer_type dofsInFaceContainer(nDofsInFace);
492  memoryInitialRequest[proc][cptFaces].resize(nDofsInFace);
493  for (int cptDof=0 ; itDof!=enDof ; ++itDof,++cptDof)
494  {
495  auto const theglobdof = itDof->get<0>();
496  auto const comp = itDof->get<1>();
497  // save the tag of mpi send
498  memoryInitialRequest[proc][cptFaces][cptDof] = theglobdof;
499  //------------------------------------------------------------------------------//
500  // get info to send
501  ublas::vector<double> nodeDofToSend( nRealDim );
502  nodeDofToSend[0]=dofPoint( theglobdof ).template get<0>()[0];
503  if ( nRealDim>1 )
504  nodeDofToSend[1]=dofPoint( theglobdof ).template get<0>()[1];
505  if ( nRealDim>2 )
506  nodeDofToSend[2]=dofPoint( theglobdof ).template get<0>()[2];
507  // up container
508  dofsInFaceContainer[cptDof] = boost::make_tuple(comp,nodeDofToSend);
509  //------------------------------------------------------------------------------//
510  }
511 
512  this->worldComm().send( proc , nbMsgToSend[proc], boost::make_tuple(itFaces->first,dofsInFaceContainer) );
513  ++nbMsgToSend[proc];
514  } // for ( int cptFaces=0 ; itFaces!=enFaces ; ++itFaces, ++cptFaces)
515  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
516 
517  //--------------------------------------------------------------------------------------------------------//
518  // counter of msg received for each process
519  std::vector<int> nbMsgToRecv;
520  mpi::all_to_all( this->worldComm(),
521  nbMsgToSend,
522  nbMsgToRecv );
523  //--------------------------------------------------------------------------------------------------------//
524  // recv dof asked and re-send
525  for ( int proc=0; proc<this->worldComm().size(); ++proc )
526  {
527  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
528  {
529  dofs_in_face_container_type dataToRecvVec;
530  this->worldComm().recv( proc, cpt, dataToRecvVec );
531 
532  auto const idFaceInMyPartition = dataToRecvVec.template get<0>();
533  DVLOG(2) << "[buildGhostInterProcessDofMap] (myRank:" << myRank << ") "
534  << "idFaceInMyPartition: " << idFaceInMyPartition << "\n";
535 
536  auto itDofInFace = dataToRecvVec.template get<1>().begin();
537  auto const enDofInFace = dataToRecvVec.template get<1>().end();
538  std::vector< size_type > resAskedWithMultiProcess(std::distance(itDofInFace,enDofInFace));
539  for ( int cptDofInFace=0 ; itDofInFace != enDofInFace ; ++itDofInFace,++cptDofInFace )
540  {
541  auto const comp = itDofInFace->template get<0>();
542  auto const nodeDofRecv = itDofInFace->template get<1>();
543 
544  //------------------------------------------------------------------------------//
545 
546  auto const& theface = mesh.face( idFaceInMyPartition );
547  auto const& elt0 = theface.element0();
548  auto const& elt1 = theface.element1();
549 
550  const bool elt0isGhost = elt0.isGhostCell();
551  auto const& eltOnProc = (elt0isGhost)?elt1:elt0;
552  auto const& eltOffProc = (elt0isGhost)?elt0:elt1;
553 
554  //------------------------------------------------------------------------------//
555  // search dof on face recv
556  int locDof = nbFaceDof;
557  bool find=false;
558 
559  if ( nDim==1 )
560  {
561  auto itdofpt = this->dofPointBegin();
562  auto const endofpt = this->dofPointEnd();
563  for ( ; itdofpt!=endofpt && !find ; ++itdofpt )
564  {
565  const auto thedofPt = itdofpt->template get<0>();
566  if ( itdofpt->template get<2>() != comp ) continue;
567 
568  DVLOG(3) << "[buildGhostInterProcessDofMap] (myRank:" << myRank << ") "
569  << "thedofPt: " << thedofPt << "nodeDofRecv: " << nodeDofRecv << "\n";
570 
571  // test equatlity of dofs point
572  bool find2=true;
573  for (uint16_type d=0;d<nRealDim;++d)
574  {
575  find2 = find2 && (std::abs( thedofPt[d]-nodeDofRecv[d] )<1e-9);
576  }
577  // if find else save local dof
578  if (find2) { locDof = itdofpt->template get<1>();find=true; }
579  }
580  // check
581  CHECK( find ) << "\nPROBLEM with parallel dof table construction : Dof point not find on interprocess face " << nodeDofRecv << "\n";
582  //------------------------------------------------------------------------------//
583  // get global dof
584  const auto dofGlobAsked = locDof;
585  // save response
586  resAskedWithMultiProcess[cptDofInFace] = this->M_mapGlobalProcessToGlobalCluster[dofGlobAsked];
587  }
588  else
589  {
590  for ( uint16_type l = 0; ( l < nbFaceDof && !find ) ; ++l )
591  {
592  // dof point in face
593  const auto thedofPtInFace = dofPoint(faceLocalToGlobal( idFaceInMyPartition, l, comp ).template get<0>()).template get<0>();
594  DVLOG(3) << "[buildGhostInterProcessDofMap] (myRank:" << myRank << ") "
595  << "thedofPtInFace: " << thedofPtInFace << "nodeDofRecv: " << nodeDofRecv << "\n";
596 
597  // test equatlity of dofs point
598  bool find2=true;
599  for (uint16_type d=0;d<nRealDim;++d)
600  {
601  find2 = find2 && (std::abs( thedofPtInFace[d]-nodeDofRecv[d] )<1e-9);
602  }
603  // if find else save local dof
604  if (find2)
605  {
606  locDof = l;
607  find=true;
608  }
609  } // for ( uint16_type l = 0; ( l < nbFaceDof && !find ) ; ++l )
610  //------------------------------------------------------------------------------//
611  // check
612  CHECK( find ) << "\nPROBLEM with parallel dof table construction : Dof point not find on interprocess face " << nodeDofRecv << "\n";
613  //------------------------------------------------------------------------------//
614  // get global dof
615  const auto thedof = faceLocalToGlobal( idFaceInMyPartition, locDof, comp );
616  const auto dofGlobAsked = thedof.template get<0>();
617  // save response
618  resAskedWithMultiProcess[cptDofInFace] = this->M_mapGlobalProcessToGlobalCluster[dofGlobAsked];
619  }
620  //------------------------------------------------------------------------------//
621  }
622  this->worldComm().send( proc, cpt, resAskedWithMultiProcess );
623 
624  } // for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
625  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
626 
627  //--------------------------------------------------------------------------------------------------------//
628  // get response to initial request and update Feel::Mesh::Faces data
629  for ( int proc=0; proc<this->worldComm().size(); ++proc )
630  {
631  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
632  {
633  //------------------------------------------------------------------------------//
634  // recv response
635  std::vector< size_type > resultRecvWithMultiProcess;
636  this->worldComm().recv( proc, cpt, resultRecvWithMultiProcess );
637  //------------------------------------------------------------------------------//
638  // iterate on dofs
639  auto itDofRes = resultRecvWithMultiProcess.begin();
640  auto const enDofRes = resultRecvWithMultiProcess.end();
641  for ( int cptDofRes=0 ; itDofRes != enDofRes ; ++itDofRes,++cptDofRes )
642  {
643  auto const myGlobProcessDof = memoryInitialRequest[proc][cpt][cptDofRes];
644  const auto dofGlobRecv = *itDofRes;
645  //update data map
646  this->M_mapGlobalProcessToGlobalCluster[myGlobProcessDof]=dofGlobRecv;
647  }
648  } // for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
649  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
650 
651 
652 
653 } // buildGhostInterProcessDofMapRecursive
654 
655 //--------------------------------------------------------------------------------------------------------//
656 //--------------------------------------------------------------------------------------------------------//
657 //--------------------------------------------------------------------------------------------------------//
658 //--------------------------------------------------------------------------------------------------------//
659 
660 template<typename MeshType, typename FEType, typename PeriodicityType>
661 void
662 DofTable<MeshType, FEType, PeriodicityType>::buildGlobalProcessToGlobalClusterDofMapDiscontinuous()
663 {
664  const int myRank = this->worldComm().rank();
665  //------------------------------------------------------------------------------//
666  // update datamap info
667  const size_type mynDofWithoutGhost = this->M_last_df[myRank] - this->M_first_df[myRank] + 1;// - setInterProcessDofNotPresent.size();
668  mpi::all_gather( this->worldComm(),
669  mynDofWithoutGhost,
670  this->M_n_localWithoutGhost_df );
671 
672  this->M_n_dofs=0;
673  for ( int proc=0; proc<this->worldComm().size(); ++proc )
674  {
675  this->M_n_dofs+=this->M_n_localWithoutGhost_df[proc];
676  }
677 
678  this->M_first_df_globalcluster=this->M_first_df;
679  this->M_last_df_globalcluster=this->M_last_df;
680  for ( int i=1; i<this->worldComm().size(); ++i )
681  {
682  this->M_first_df_globalcluster[i]=this->M_last_df_globalcluster[i-1]+1;
683  this->M_last_df_globalcluster[i]=this->M_first_df_globalcluster[i]+this->M_n_localWithoutGhost_df[i]-1;
684  }
685  //------------------------------------------------------------------------------//
686  // init map
687  this->M_mapGlobalProcessToGlobalCluster.resize( this->M_n_localWithGhost_df[myRank],invalid_size_type_value );
688  this->M_mapGlobalClusterToGlobalProcess.resize( this->M_n_localWithoutGhost_df[myRank],invalid_size_type_value );
689  //------------------------------------------------------------------------------//
690  // add in map the dofs presents
691  size_type firstGlobIndex = this->M_first_df_globalcluster[myRank];
692  size_type nextGlobIndex = firstGlobIndex;
693  for ( size_type i=0; i< this->M_n_localWithGhost_df[myRank]; ++i )
694  {
695  //if ( setInterProcessDofNotPresent.find( i )==setInterProcessDofNotPresent.end() )
696  {
697  this->M_mapGlobalProcessToGlobalCluster[i]=nextGlobIndex;
698  this->M_mapGlobalClusterToGlobalProcess[nextGlobIndex-firstGlobIndex]=i;
699  ++nextGlobIndex;
700  }
701  }
702  //------------------------------------------------------------------------------//
703 
704 } // buildGlobalProcessToGlobalClusterDofMap
705 
706 
707 
708 
709 
710 
711 //--------------------------------------------------------------------------------------------------------//
712 //--------------------------------------------------------------------------------------------------------//
713 //--------------------------------------------------------------------------------------------------------//
714 //--------------------------------------------------------------------------------------------------------//
715 //--------------------------------------------------------------------------------------------------------//
716 //--------------------------------------------------------------------------------------------------------//
717 //--------------------------------------------------------------------------------------------------------//
718 //--------------------------------------------------------------------------------------------------------//
719 
720 template<typename MeshType, typename FEType, typename PeriodicityType>
721 void
722 DofTable<MeshType, FEType, PeriodicityType>::buildGhostInterProcessDofMap( mesh_type& mesh,
723  std::map<size_type,boost::tuple<size_type,size_type> > & mapInterProcessDof )
724 {
725  //------------------------------------------------------------------------------//
726  // if nbFaceDof == 0 else return
728  if ( !fe_type::is_modal )
729  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
730  face_type::numEdges * fe_type::nDofPerEdge +
731  face_type::numFaces * fe_type::nDofPerFace );
732  else
733  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
734 
735  if ( nbFaceDof == 0 ) return;
736 
737  //------------------------------------------------------------------------------//
738  // init dof list to search with interprocess dof
739  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > listToSend(this->worldComm().size());
740  this->buildGhostInterProcessDofMapInit(mesh,listToSend);
741  //------------------------------------------------------------------------------//
742  // search real and ghost dofs and build mapping between them
743  // the search is iterative beacause some dofs are connected over 2 proc
744  // and we need to search all proc around these dofs
745  std::vector< std::set<size_type > > memoryFace(this->worldComm().size());
746  bool allDofAreFind =false;
747  while (!allDofAreFind)
748  {
749  auto const& res = this->buildGhostInterProcessDofMapRecursive(mesh,listToSend,
750  mapInterProcessDof,
751  memoryFace);
752 
753  allDofAreFind = res.template get<0>();
754  if (!allDofAreFind) listToSend = res.template get<1>();
755  }
756  //------------------------------------------------------------------------------//
757 }
758 
759 //--------------------------------------------------------------------------------------------------------//
760 //--------------------------------------------------------------------------------------------------------//
761 //--------------------------------------------------------------------------------------------------------//
762 //--------------------------------------------------------------------------------------------------------//
763 
764 
765 template<typename MeshType, typename FEType, typename PeriodicityType>
766 void
767 DofTable<MeshType, FEType, PeriodicityType>::buildGhostInterProcessDofMapInit( mesh_type& mesh,
768  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > & listToSend )
769 {
770 
771  // goal init container listToSend
772  // std::vector< std::map<size_type,std::set<size_type> > > ( proc,( idFace,(globDof,..)), ... )) )
773  const int myRank = this->worldComm().rank();
774  //------------------------------------------------------------------------------//
775  // get nbFaceDof
777  if ( !fe_type::is_modal )
778  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
779  face_type::numEdges * fe_type::nDofPerEdge +
780  face_type::numFaces * fe_type::nDofPerFace );
781  else
782  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
783 
784  if ( nbFaceDof == 0 ) return;
785  //------------------------------------------------------------------------------//
786  std::vector<bool> dofdone(this->M_n_localWithGhost_df[myRank],false);
787  // iteration on all interprocessfaces in order to send requests to the near proc
788  auto face_it = mesh.interProcessFaces().first;
789  auto const face_en = mesh.interProcessFaces().second;
790  DVLOG(2) << "[buildGhostInterProcessDofMap] nb interprocess faces: " << std::distance( face_it, face_en ) << "\n";
791  for ( ; face_it != face_en ; ++face_it )
792  {
793  DVLOG(2) << "[buildGhostInterProcessDofMap] face id: " << face_it->id() << "\n";
794  auto const& elt0 = face_it->element0();
795  auto const& elt1 = face_it->element1();
796  const bool elt0isGhost = elt0.isGhostCell();
797  auto const& eltOnProc = (elt0isGhost)?elt1:elt0;
798  auto const& eltOffProc = (elt0isGhost)?elt0:elt1;
799 
800 #if 0
801  element_type eltOnProc,eltOffProc;
802  if ( elt0.processId()!=myRank )
803  {
804  eltOffProc = elt0;
805  eltOnProc=elt1;
806  }
807  else if ( elt1.processId()!=myRank )
808  {
809  eltOffProc = elt1;
810  eltOnProc=elt0;
811  }
812  else std::cout << "\n[buildGhostInterProcessDofMapInit] PROBLEME1!!!!" << std::endl;
813 #endif
814  //------------------------------------------------------------------------------//
815 
816  const int IdProcessOfGhost = eltOffProc.processId();
817  const size_type idFaceInPartition = face_it->idInOthersPartitions( IdProcessOfGhost );
818  //------------------------------------------------------------------------------//
819  // for each dof in face
820  const uint16_type ncdof = is_product?nComponents:1;
821  // subcontainer
822  std::set<boost::tuple<size_type,uint16_type> > dofSetAdd;
823 
824  for ( uint16_type l = 0; l < nbFaceDof; ++l )
825  {
826  for ( uint16_type c = 0; c < ncdof; ++c )
827  {
828  // add dof in subcontainer
829  const size_type theglobdof = faceLocalToGlobal( face_it->id(),l,c ).template get<0>();
830  if (!dofdone[theglobdof])
831  {
832  dofSetAdd.insert(boost::make_tuple(theglobdof,c));
833  dofdone[theglobdof]=true;
834  }
835  }
836  // add composante dofs in result container
837  }
838  listToSend[IdProcessOfGhost].insert( std::make_pair( idFaceInPartition/*face_it->id()*/, dofSetAdd ) );
839 
840 
841  } // for ( ; face_it != face_en ; ++face_it )
842 }
843 
844 //--------------------------------------------------------------------------------------------------------//
845 //--------------------------------------------------------------------------------------------------------//
846 //--------------------------------------------------------------------------------------------------------//
847 //--------------------------------------------------------------------------------------------------------//
848 
849 namespace detail {
850 
851 template<typename element_type>
852 void
853 searchPartitionAroundNode( const int myRank, ublas::vector<double> const& nodeSearched,
854  element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
855  std::set<boost::tuple<int,size_type> > & multiProcessRes )
856 {
857  for ( uint16_type f = 0; f < element_type::numTopologicalFaces; ++f )
858  {
859  auto const& theFace = eltOnProc.face(f);
860  size_type theFaceId = theFace.id();
861  if (memoryIdsFaces.find(theFaceId) != memoryIdsFaces.end()) continue;
862 
863  // save face
864  memoryIdsFaces.insert( theFaceId );
865  // face contains node?
866  bool find=false;
867  auto const& facevertices = theFace.vertices();
868  for (uint16_type v=0;v<element_type::face_type::numVertices && !find ;++v)
869  {
870  bool find2=true;
871  for (uint16_type d=0;d<element_type::nDim;++d)
872  {
873  find2 = find2 && ( std::abs( facevertices(d,v)-nodeSearched[d] )<1e-9 );
874  }
875  if (find2) find=true;
876  }
877 
878  if (find)
879  {
880  if ( theFace.isConnectedTo0() )
881  {
882  auto const& eltOnProc0 = theFace.element0();
883  if (eltOnProc0.processId()!=myRank)
884  multiProcessRes.insert( boost::make_tuple(eltOnProc0.processId(),theFace.idInOthersPartitions(eltOnProc0.processId()) ));
885  else if (eltOnProc.id()!=eltOnProc0.id())
886  searchPartitionAroundNode(myRank,nodeSearched,eltOnProc0,memoryIdsFaces,multiProcessRes );
887  }
888  if ( theFace.isConnectedTo1() )
889  {
890  auto const& eltOnProc1 = theFace.element1();
891  if (eltOnProc1.processId()!=myRank)
892  multiProcessRes.insert( boost::make_tuple(eltOnProc1.processId(),theFace.idInOthersPartitions(eltOnProc1.processId()) ));
893  else if (eltOnProc.id()!=eltOnProc1.id())
894  searchPartitionAroundNode(myRank,nodeSearched,eltOnProc1,memoryIdsFaces,multiProcessRes );
895  }
896  } // if (find)
897 
898  } // for ( uint16_type f = 0; f < element_type::numTopologicalFaces; ++f )
899 
900 } // searchPartitionAroundNode
901 
902 
903 template<typename mesh_type>
904 void
905 searchPartitionAroundEdge( const int myRank,
906  typename mesh_type::face_type const& faceContainEdge,const uint16_type idEdgesInFace,
907  typename mesh_type::element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
908  std::set<boost::tuple<int,size_type> > & multiProcessRes, bool firstExp, mpl::int_<0> )
909 {}
910 
911 template<typename mesh_type>
912 void
913 searchPartitionAroundEdge( const int myRank,
914  typename mesh_type::face_type const& faceContainEdge,const uint16_type idEdgesInFace,
915  typename mesh_type::element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
916  std::set<boost::tuple<int,size_type> > & multiProcessRes, bool firstExp, mpl::int_<1> )
917 {}
918 
919 template<typename mesh_type>
920 void
921 searchPartitionAroundEdge( const int myRank,
922  typename mesh_type::face_type const& faceContainEdge,const uint16_type idEdgesInFace,
923  typename mesh_type::element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
924  std::set<boost::tuple<int,size_type> > & multiProcessRes, bool firstExp, mpl::int_<2> )
925 {}
926 
927 template<typename mesh_type>
928 void
929 searchPartitionAroundEdge( const int myRank,
930  typename mesh_type::face_type const& faceContainEdge,const uint16_type idEdgesInFace,
931  typename mesh_type::element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
932  std::set<boost::tuple<int,size_type> > & multiProcessRes, bool firstExp, mpl::int_<3> )
933 {
934 
935  // only usefull whith first pass
936  uint16_type iFaEl=0,iEdEl=0;
937  if (firstExp)
938  {
939  iFaEl = ( faceContainEdge.processId() == faceContainEdge.proc_first() )? faceContainEdge.pos_first():faceContainEdge.pos_second();
940  //local edge number (in element)
941  iEdEl = mesh_type::element_type::fToE( iFaEl, idEdgesInFace );
942  }
943 
944  auto const& theedge = (firstExp)?eltOnProc.edge(iEdEl):faceContainEdge.edge(idEdgesInFace);
945  //auto const& theEdgesVertices = theedge.vertices();
946 
947  for ( uint16_type f = 0; f < mesh_type::element_type::numTopologicalFaces; ++f )
948  {
949  auto const& theFace = eltOnProc.face(f);
950 
951  size_type theFaceId = theFace.id();
952  //auto const& theneighborface = mesh.face( eltOnProc.face(f).id() );
953 
954  if (memoryIdsFaces.find(theFaceId) != memoryIdsFaces.end()) continue;
955 
956  // save face
957  memoryIdsFaces.insert( theFaceId );
958 
959  // face contains edge?
960  bool find=false;uint16_type newEdgeIdInFace=0;
961  for (uint16_type ed=0;ed<mesh_type::element_type::face_type::numEdges && !find ;++ed)
962  {
963  //auto const& otheredges = theFace.edge(ed);
964  if (theFace.edge(ed).id() == theedge.id())
965  {
966  find=true; newEdgeIdInFace=ed;
967 #if 0
968  // check
969  auto const& newEdgesVertices = theFace.edge(ed).vertices();
970  //auto const& theEdgesVertices = theedge.vertices();
971  bool find2=true;
972  for (uint16_type d=0;d<mesh_type::element_type::nDim;++d)
973  {
974  find2 = find2 && ( std::abs( newEdgesVertices(d,0)-theEdgesVertices(d,0) )<1e-9 );
975  }
976  if (!find2) std::cout << "\n IAOIOIOAIOIO"<<std::endl;
977  //if (newEdgesVertices!=theEdgesVertices) std::cout << "\IAOIOIOAIOIO"<<std::endl;
978 #endif
979  }
980  }
981 
982  if (find)
983  {
984  if ( theFace.isConnectedTo0() )
985  {
986  auto const& eltOnProc0 = theFace.element0();
987  if (eltOnProc0.processId()!=myRank)
988  multiProcessRes.insert( boost::make_tuple(eltOnProc0.processId(),theFace.idInOthersPartitions(eltOnProc0.processId()) ));
989  else if (eltOnProc.id()!=eltOnProc0.id())
990  searchPartitionAroundEdge<mesh_type>(myRank,theFace,newEdgeIdInFace,eltOnProc0,memoryIdsFaces,multiProcessRes,false,mpl::int_<3>());
991  }
992  if ( theFace.isConnectedTo1() )
993  {
994  auto const& eltOnProc1 = theFace.element1();
995  if (eltOnProc1.processId()!=myRank)
996  multiProcessRes.insert( boost::make_tuple(eltOnProc1.processId(),theFace.idInOthersPartitions(eltOnProc1.processId()) ));
997  else if (eltOnProc.id()!=eltOnProc1.id())
998  searchPartitionAroundEdge<mesh_type>(myRank,theFace,newEdgeIdInFace,eltOnProc1,memoryIdsFaces,multiProcessRes,false,mpl::int_<3>());
999  }
1000  } // if (find)
1001 
1002  } // for ( uint16_type f = 0; f < element_type::numTopologicalFaces; ++f )
1003 
1004 }
1005 
1006 template<typename mesh_type>
1007 void
1008 searchPartitionAroundEdge( const int myRank,
1009  typename mesh_type::face_type const& faceContainEdge,const uint16_type idEdgesInFace,
1010  typename mesh_type::element_type const& eltOnProc, std::set<size_type> & memoryIdsFaces,
1011  std::set<boost::tuple<int,size_type> > & multiProcessRes )
1012 {
1013  searchPartitionAroundEdge<mesh_type>(myRank,faceContainEdge,idEdgesInFace,eltOnProc,memoryIdsFaces,multiProcessRes,true,mpl::int_<mesh_type::nDim>());
1014 }
1015 
1016 } // namespace detail
1017 
1018 template<typename MeshType, typename FEType, typename PeriodicityType>
1019 boost::tuple<bool, std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > >
1020 DofTable<MeshType, FEType, PeriodicityType>::buildGhostInterProcessDofMapRecursive( mesh_type& mesh,
1021  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > const& listToSend,
1022  std::map<size_type,boost::tuple<size_type,size_type> > & mapInterProcessDof,
1023  std::vector< std::set<size_type > > & memoryFace )
1024 {
1025  const int myRank = this->worldComm().rank();
1026  //--------------------------------------------------------------------------------------------------------//
1027  size_type nbFaceDof = invalid_size_type_value;
1028  if ( !fe_type::is_modal )
1029  nbFaceDof = ( face_type::numVertices * fe_type::nDofPerVertex +
1030  face_type::numEdges * fe_type::nDofPerEdge +
1031  face_type::numFaces * fe_type::nDofPerFace );
1032 
1033  else
1034  nbFaceDof = face_type::numVertices * fe_type::nDofPerVertex;
1035  //--------------------------------------------------------------------------------------------------------//
1036 
1037  std::vector<int> nbMsgToSend( this->worldComm().size(), 0 );
1038 
1039  std::vector< std::map<size_type,std::set<boost::tuple<size_type,uint16_type> > > > listNeedMoreSearch(this->worldComm().size());
1040  int counterListNeedMoreSearch = 0, counterListNeedMoreSearchLoc = 0;
1041 
1042  std::vector< std::vector< std::vector<boost::tuple<uint16_type,int> > > > memoryInitialRequest( this->worldComm().size() );
1043 
1044  typedef std::vector< boost::tuple<uint16_type, ublas::vector<double> > > dofs_in_face_subcontainer_type;
1045  typedef boost::tuple<size_type, dofs_in_face_subcontainer_type > dofs_in_face_container_type;
1046 
1047  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1048  {
1049  auto itFaces = listToSend[proc].begin();
1050  auto const enFaces = listToSend[proc].end();
1051  memoryInitialRequest[proc].resize(std::distance(itFaces,enFaces));
1052  for ( int cptFaces=0 ; itFaces!=enFaces ; ++itFaces, ++cptFaces)
1053  {
1054  auto itDof = itFaces->second.begin();
1055  auto const enDof = itFaces->second.end();
1056  auto const nDofsInFace = std::distance(itDof,enDof);
1057  dofs_in_face_subcontainer_type dofsInFaceContainer(nDofsInFace);
1058  memoryInitialRequest[proc][cptFaces].resize(nDofsInFace);
1059  for (int cptDof=0 ; itDof!=enDof ; ++itDof,++cptDof)
1060  {
1061  auto const theglobdof = itDof->get<0>();
1062  auto const comp = itDof->get<1>();
1063  // save the tag of mpi send
1064  memoryInitialRequest[proc][cptFaces][cptDof] = boost::make_tuple(comp,theglobdof);
1065  //------------------------------------------------------------------------------//
1066  // get info to send
1067  ublas::vector<double> nodeDofToSend( nDim );
1068  nodeDofToSend[0]=dofPoint( theglobdof ).template get<0>()[0];
1069  if ( nDim>1 )
1070  nodeDofToSend[1]=dofPoint( theglobdof ).template get<0>()[1];
1071  if ( nDim>2 )
1072  nodeDofToSend[2]=dofPoint( theglobdof ).template get<0>()[2];
1073  // up container
1074  dofsInFaceContainer[cptDof] = boost::make_tuple(comp,nodeDofToSend);
1075  //------------------------------------------------------------------------------//
1076  }
1077 
1078  this->worldComm().send( proc , nbMsgToSend[proc], boost::make_tuple(itFaces->first,dofsInFaceContainer) );
1079  ++nbMsgToSend[proc];
1080  } // for ( int cptFaces=0 ; itFaces!=enFaces ; ++itFaces, ++cptFaces)
1081  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
1082 
1083  //--------------------------------------------------------------------------------------------------------//
1084  // counter of msg received for each process
1085  std::vector<int> nbMsgToRecv;
1086  mpi::all_to_all( this->worldComm(),
1087  nbMsgToSend,
1088  nbMsgToRecv );
1089  //--------------------------------------------------------------------------------------------------------//
1090  // recv dof asked and re-send
1091  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1092  {
1093  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
1094  {
1095  dofs_in_face_container_type dataToRecvVec;
1096  this->worldComm().recv( proc, cpt, dataToRecvVec );
1097 
1098  auto const idFaceInMyPartition = dataToRecvVec.template get<0>();
1099 
1100  auto itDofInFace = dataToRecvVec.template get<1>().begin();
1101  auto const enDofInFace = dataToRecvVec.template get<1>().end();
1102  std::vector< boost::tuple<size_type,std::vector<boost::tuple<int,size_type> > > > resAskedWithMultiProcess(std::distance(itDofInFace,enDofInFace));
1103  for ( int cptDofInFace=0 ; itDofInFace != enDofInFace ; ++itDofInFace,++cptDofInFace )
1104  {
1105  auto const comp = itDofInFace->template get<0>();
1106  auto const nodeDofRecv = itDofInFace->template get<1>();
1107 
1108  //------------------------------------------------------------------------------//
1109 
1110  auto const& theface = mesh.face( idFaceInMyPartition );
1111  auto const& elt0 = theface.element0();
1112  auto const& elt1 = theface.element1();
1113 
1114  const bool elt0isGhost = elt0.isGhostCell();
1115  auto const& eltOnProc = (elt0isGhost)?elt1:elt0;
1116  auto const& eltOffProc = (elt0isGhost)?elt0:elt1;
1117 
1118 #if 0
1119  element_type eltOnProc, eltOffProc;
1120  uint16_type faceIdInEltOnProc = invalid_uint16_type_value;
1121  if ( elt0.processId()!=myRank )
1122  {
1123  eltOffProc = elt0;
1124  eltOnProc=elt1;
1125  faceIdInEltOnProc = theface.pos_second();
1126  }
1127  else if ( elt1.processId()!=myRank )
1128  {
1129  eltOffProc = elt1;
1130  eltOnProc=elt0;
1131  faceIdInEltOnProc = theface.pos_first();
1132  }
1133  else
1134  {
1135  CHECK( (elt0.processId()==myRank) ||
1136  (elt1.processId()==myRank) )
1137  << "\nPROBLEM with parallel dof table construction\n"
1138  << " elt0.processId() " << elt0.processId()
1139  << " elt1.processId() " << elt1.processId()
1140  << "\n";
1141  }
1142 #endif
1143 
1144  //------------------------------------------------------------------------------//
1145  // search dof on face recv
1146  int locDof = nbFaceDof;
1147  bool find=false;
1148 
1149  for ( uint16_type l = 0; ( l < nbFaceDof && !find ) ; ++l )
1150  {
1151  // dof point in face
1152  const auto thedofPtInFace = dofPoint(faceLocalToGlobal( idFaceInMyPartition, l, comp ).template get<0>()).template get<0>();
1153  // test equatlity of dofs point
1154  bool find2=true;
1155  for (uint16_type d=0;d<nDim;++d)
1156  {
1157  find2 = find2 && (std::abs( thedofPtInFace[d]-nodeDofRecv[d] )<1e-9);
1158  }
1159  // if find else save local dof
1160  if (find2)
1161  {
1162  locDof = l;
1163  find=true;
1164  }
1165  } // for ( uint16_type l = 0; ( l < nbFaceDof && !find ) ; ++l )
1166  //------------------------------------------------------------------------------//
1167  // check
1168  //if ( !find ) std::cout<<"\n [buildGhostInterProcessDofMap] : Dof point not find on interprocess face " << nodeDofRecv << std::endl;
1169  CHECK( find ) << "\nPROBLEM with parallel dof table construction : Dof point not find on interprocess face " << nodeDofRecv << "\n";
1170  //------------------------------------------------------------------------------//
1171  // get global dof
1172  const auto thedof = faceLocalToGlobal( idFaceInMyPartition, locDof, comp );
1173  const auto dofGlobAsked = thedof.template get<0>();
1174  //------------------------------------------------------------------------------//
1175  //------------------------------------------------------------------------------//
1176  //------------------------------------------------------------------------------//
1177  //------------------------------------------------------------------------------//
1178  // search maybe with neighboors elt
1179  std::set<size_type> memoryIdsFaces;
1180  std::set<boost::tuple<int,size_type> > multiProcessRes;
1181  if ( locDof < face_type::numVertices*fe_type::nDofPerVertex)
1182  {
1183  detail::searchPartitionAroundNode( myRank, nodeDofRecv, eltOnProc, memoryIdsFaces, multiProcessRes );
1184  }
1185  else if ( nDim == 3 && locDof < (face_type::numVertices*fe_type::nDofPerVertex + face_type::numEdges*fe_type::nDofPerEdge) )
1186  {
1187  int locDofInEgde = locDof - face_type::numVertices*fe_type::nDofPerVertex;
1188 
1189  const int nDofPerEdgeTemp = mpl::if_<boost::is_same<mpl::int_<fe_type::nDofPerEdge>,mpl::int_<0> >,
1190  mpl::int_<1>,
1191  mpl::int_<fe_type::nDofPerEdge> >::type::value;
1192 
1193  int edgeGetLocDof = locDofInEgde / nDofPerEdgeTemp;
1194 
1195  detail::searchPartitionAroundEdge<mesh_type>( myRank, theface, edgeGetLocDof, eltOnProc, memoryIdsFaces, multiProcessRes );
1196  }
1197  //if (proc==0 ) std::cout << "myRank " << myRank << " multiProcessRes.size() " << multiProcessRes.size() << " " << nodeDofRecv << std::endl;
1198  //------------------------------------------------------------------------------//
1199  // convert in vector for send
1200  auto const nbMultiProcessRes = multiProcessRes.size();
1201  std::vector<boost::tuple<int,size_type> > multiProcessResVec( nbMultiProcessRes );
1202  std::copy( multiProcessRes.begin(),multiProcessRes.end(),
1203  multiProcessResVec.begin() );
1204  //------------------------------------------------------------------------------//
1205  //------------------------------------------------------------------------------//
1206  //------------------------------------------------------------------------------//
1207  // send response
1208  resAskedWithMultiProcess[cptDofInFace] = boost::make_tuple(dofGlobAsked,multiProcessResVec);
1209  //------------------------------------------------------------------------------//
1210  }
1211  this->worldComm().send( proc, cpt, resAskedWithMultiProcess );
1212 
1213  } // for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
1214  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
1215 
1216  //--------------------------------------------------------------------------------------------------------//
1217  // get response to initial request and update Feel::Mesh::Faces data
1218  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1219  {
1220  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1221  {
1222  //------------------------------------------------------------------------------//
1223  // recv response
1224  std::vector< boost::tuple<size_type,std::vector<boost::tuple<int,size_type> > > > resultRecvWithMultiProcess;
1225  this->worldComm().recv( proc, cpt, resultRecvWithMultiProcess );
1226  //------------------------------------------------------------------------------//
1227  // iterate on dofs
1228  auto itDofRes = resultRecvWithMultiProcess.begin();
1229  auto const enDofRes = resultRecvWithMultiProcess.end();
1230  for ( int cptDofRes=0 ; itDofRes != enDofRes ; ++itDofRes,++cptDofRes )
1231  {
1232 
1233  auto const mycomp = memoryInitialRequest[proc][cpt][cptDofRes].template get<0>();
1234  auto const myGlobProcessDof = memoryInitialRequest[proc][cpt][cptDofRes].template get<1>();
1235  const auto dofGlobRecv = itDofRes->template get<0>();
1236  const auto multiProcessRes = itDofRes->template get<1>();
1237 
1238  //------------------------------------------------------------------------------//
1239  // update mapInterProcessDof data
1240  if ( mapInterProcessDof.find( myGlobProcessDof ) == mapInterProcessDof.end() )
1241  {
1242  mapInterProcessDof.insert( std::make_pair( myGlobProcessDof, boost::make_tuple( proc,dofGlobRecv ) ) );
1243  }
1244  else if ( ( int )mapInterProcessDof.find( myGlobProcessDof )->second.template get<0>() > proc )
1245  {
1246  mapInterProcessDof[ myGlobProcessDof ] = boost::make_tuple( proc,dofGlobRecv );
1247  }
1248  //------------------------------------------------------------------------------//
1249  // save global dof with search on proc
1250  memoryFace[proc].insert(myGlobProcessDof);
1251  //------------------------------------------------------------------------------//
1252  // add maybe more search
1253  auto itMP = multiProcessRes.begin();
1254  auto const enMP = multiProcessRes.end();
1255  for ( ; itMP!=enMP ;++itMP )
1256  {
1257  auto const eltMultiProcessProcId = itMP->template get<0>();
1258  auto const eltMultiProcessFaceId = itMP->template get<1>();
1259  if ( memoryFace[eltMultiProcessProcId].find( myGlobProcessDof ) == memoryFace[eltMultiProcessProcId].end() )
1260  {
1261  listNeedMoreSearch[eltMultiProcessProcId][eltMultiProcessFaceId].insert(boost::make_tuple(myGlobProcessDof,mycomp));
1262  ++counterListNeedMoreSearchLoc;
1263  }
1264  }
1265  }
1266  } // for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
1267  } // for ( int proc=0; proc<this->worldComm().size(); ++proc )
1268 
1269  //--------------------------------------------------------------------------------------------------------//
1270  // need to apply an all_reduce
1271  mpi::all_reduce( this->worldComm().globalComm(),
1272  counterListNeedMoreSearchLoc,
1273  counterListNeedMoreSearch,
1274  std::plus<int>() );
1275  //std::cout<< " counterListNeedMoreSearch "<< counterListNeedMoreSearch << " counterListNeedMoreSearchLoc " << counterListNeedMoreSearchLoc << std::endl;
1276  //--------------------------------------------------------------------------------------------------------//
1277  return boost::make_tuple( (bool) (counterListNeedMoreSearch == 0),
1278  listNeedMoreSearch);
1279 
1280 
1281 } // buildGhostInterProcessDofMapRecursive
1282 
1283 //--------------------------------------------------------------------------------------------------------//
1284 //--------------------------------------------------------------------------------------------------------//
1285 //--------------------------------------------------------------------------------------------------------//
1286 //--------------------------------------------------------------------------------------------------------//
1287 
1288 template<typename MeshType, typename FEType, typename PeriodicityType>
1289 void
1290 DofTable<MeshType, FEType, PeriodicityType>::buildDofNotPresent( std::map<size_type,boost::tuple<size_type,size_type> > const & mapInterProcessDof,
1291  std::map<size_type,boost::tuple<size_type,size_type> > & setInterProcessDofNotPresent )
1292 
1293 {
1294  const int myRank = this->worldComm().rank();
1295 
1296  auto it_mapIPDof = mapInterProcessDof.begin();
1297  auto const en_mapIPDof = mapInterProcessDof.end();
1298  for ( ; it_mapIPDof!=en_mapIPDof; ++it_mapIPDof )
1299  {
1300  if ( it_mapIPDof->second.template get<0>() < myRank )
1301  setInterProcessDofNotPresent[it_mapIPDof->first] = it_mapIPDof->second;
1302  }
1303 }
1304 
1305 //--------------------------------------------------------------------------------------------------------//
1306 //--------------------------------------------------------------------------------------------------------//
1307 //--------------------------------------------------------------------------------------------------------//
1308 //--------------------------------------------------------------------------------------------------------//
1309 
1310 template<typename MeshType, typename FEType, typename PeriodicityType>
1311 void
1312 DofTable<MeshType, FEType, PeriodicityType>::buildGlobalProcessToGlobalClusterDofMap( mesh_type& mesh,
1313  std::map<size_type,boost::tuple<size_type,size_type> > const& setInterProcessDofNotPresent )
1314 {
1315  const int myRank = this->worldComm().rank();
1316 
1317  //------------------------------------------------------------------------------//
1318  // update datamap info
1319  //const size_type mynDofWithoutGhost = this->M_last_df[myRank] - this->M_first_df[myRank] + 1 - setInterProcessDofNotPresent.size();
1320  const size_type mynDofWithoutGhost = this->M_n_localWithGhost_df[myRank] - setInterProcessDofNotPresent.size();
1321  mpi::all_gather( this->worldComm(),
1322  mynDofWithoutGhost,
1323  this->M_n_localWithoutGhost_df );
1324 
1325  this->M_n_dofs=0;
1326  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1327  {
1328  this->M_n_dofs+=this->M_n_localWithoutGhost_df[proc];
1329  }
1330 
1331  this->M_first_df_globalcluster=this->M_first_df;
1332  this->M_last_df_globalcluster=this->M_last_df;
1333  for ( int i=1; i<this->worldComm().size(); ++i )
1334  {
1335  if ( this->M_n_localWithoutGhost_df[i-1] >0 )
1336  this->M_first_df_globalcluster[i]=this->M_last_df_globalcluster[i-1]+1;
1337  else
1338  this->M_first_df_globalcluster[i]=this->M_last_df_globalcluster[i-1];
1339 
1340  if ( this->M_n_localWithoutGhost_df[i] >0 )
1341  this->M_last_df_globalcluster[i]=this->M_first_df_globalcluster[i]+this->M_n_localWithoutGhost_df[i]-1;
1342  else
1343  this->M_last_df_globalcluster[i]=this->M_first_df_globalcluster[i];
1344  }
1345 
1346  //------------------------------------------------------------------------------//
1347  // init map
1348  this->M_mapGlobalProcessToGlobalCluster.resize( this->M_n_localWithGhost_df[myRank],invalid_size_type_value );
1349  this->M_mapGlobalClusterToGlobalProcess.resize( this->M_n_localWithoutGhost_df[myRank],invalid_size_type_value );
1350 
1351  //------------------------------------------------------------------------------//
1352  // add in map the dofs presents
1353  size_type firstGlobIndex = this->M_first_df_globalcluster[myRank];
1354  size_type nextGlobIndex = firstGlobIndex;
1355  for ( size_type i=0; i< this->M_n_localWithGhost_df[myRank]; ++i )
1356  {
1357  if ( setInterProcessDofNotPresent.find( i )==setInterProcessDofNotPresent.end() )
1358  {
1359  this->M_mapGlobalProcessToGlobalCluster[i]=nextGlobIndex;
1360  this->M_mapGlobalClusterToGlobalProcess[nextGlobIndex-firstGlobIndex]=i;
1361  ++nextGlobIndex;
1362  }
1363  }
1364 
1365  //------------------------------------------------------------------------------//
1366  // update datamap with dofs non presents (ghosts)
1367  this->updateGhostGlobalDof( setInterProcessDofNotPresent );
1368 
1369 } // buildGlobalProcessToGlobalClusterDofMap
1370 
1371 //--------------------------------------------------------------------------------------------------------//
1372 //--------------------------------------------------------------------------------------------------------//
1373 //--------------------------------------------------------------------------------------------------------//
1374 //--------------------------------------------------------------------------------------------------------//
1375 
1376 template<typename MeshType, typename FEType, typename PeriodicityType>
1377 void
1378 DofTable<MeshType, FEType, PeriodicityType>::updateGhostGlobalDof( std::map<size_type,boost::tuple<size_type,size_type> > const& setInterProcessDofNotPresent )
1379 {
1380  std::vector<int> nbMsgToSend2( this->worldComm().size(), 0 );
1381  std::vector< std::map<int,size_type> > mapMsg2( this->worldComm().size() );
1382 
1383  auto dofNP_it = setInterProcessDofNotPresent.begin();
1384  auto dofNP_en = setInterProcessDofNotPresent.end();
1385  for ( ; dofNP_it!=dofNP_en ; ++dofNP_it )
1386  {
1387  auto const globalProcessDof = dofNP_it->first;
1388  auto const IdProcessOfGhost = dofNP_it->second.template get<0>();
1389  auto const dofFaceInPartition = dofNP_it->second.template get<1>();
1390 
1391  this->worldComm().send( IdProcessOfGhost , nbMsgToSend2[IdProcessOfGhost], dofFaceInPartition );
1392  mapMsg2[IdProcessOfGhost].insert( std::make_pair( nbMsgToSend2[IdProcessOfGhost],globalProcessDof ) );
1393  ++nbMsgToSend2[IdProcessOfGhost];
1394  }
1395 
1396  // counter of msg received for each process
1397  std::vector<int> nbMsgToRecv2;
1398  mpi::all_to_all( this->worldComm(),
1399  nbMsgToSend2,
1400  nbMsgToRecv2 );
1401 
1402  // recv dof asked and re-send dof in this proc
1403  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1404  {
1405  for ( int cpt=0; cpt<nbMsgToRecv2[proc]; ++cpt )
1406  {
1407  //recv
1408  size_type dofRecv;
1409  this->worldComm().recv( proc, cpt, dofRecv );
1410  // send
1411  size_type const dofToSend= this->M_mapGlobalProcessToGlobalCluster[dofRecv];
1412  this->worldComm().send( proc, cpt,dofToSend );
1413  }
1414  }
1415 
1416  // get response to initial request and update Feel::Mesh::Faces data
1417  for ( int proc=0; proc<this->worldComm().size(); ++proc )
1418  {
1419  for ( int cpt=0; cpt<nbMsgToSend2[proc]; ++cpt )
1420  {
1421  //recv
1422  size_type dofGlobClusterRecv;
1423  this->worldComm().recv( proc, cpt, dofGlobClusterRecv );
1424  //update data
1425  this->M_mapGlobalProcessToGlobalCluster[mapMsg2[proc][cpt]]=dofGlobClusterRecv;
1426  }
1427  }
1428 
1429 }
1430 
1431 
1432 
1433 } // namespace Feel
1434 
1435 #endif /* FEELPP_DOFTABLE_MPI_HPP */

Generated on Sun Oct 20 2013 08:24:56 for Feel++ by doxygen 1.8.4