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
createsubmesh.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Vincent Chabannes <vincent.chabannes@imag.fr>
6  Date: 2011-07-21
7 
8  Copyright (C) 2011 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
30 #ifndef __createsubmesh_H
31 #define __createsubmesh_H 1
32 
33 #include <boost/mpl/if.hpp>
34 #include <boost/mpl/identity.hpp>
35 #include <boost/spirit/home/phoenix/stl/algorithm.hpp>
37 
38 
39 namespace Feel
40 {
41 template <typename C, typename V, int T> class Mesh;
42 
43 template <typename MeshType,typename IteratorRange, int TheTag=MeshType::tag>
44 class createSubmeshTool
45 {
46 public :
47 
48  typedef IteratorRange range_type;
49  typedef typename boost::tuples::template element<0, range_type>::type idim_type;
50  typedef typename boost::tuples::template element<1, range_type>::type iterator_type;
51 
52  static const uint16_type tag = TheTag;
53  typedef MeshType mesh_type;
54  typedef typename mesh_type::value_type value_type;
55 
56  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
57 
58  typedef typename mpl::if_<mpl::bool_<mesh_type::shape_type::is_simplex>,
59  mpl::identity< Mesh< Simplex< mesh_type::nDim-1,mesh_type::nOrder,mesh_type::nRealDim>, value_type, tag > >,
60  mpl::identity< Mesh< Hypercube<mesh_type::nDim-1,mesh_type::nOrder,mesh_type::nRealDim>, value_type, tag > > >::type::type mesh_faces_type;
61 
62  typedef boost::shared_ptr<mesh_faces_type> mesh_faces_ptrtype;
63 
64  typedef typename mpl::if_<mpl::bool_<mesh_type::shape_type::is_simplex>,
65  mpl::identity< Mesh< Simplex< (mesh_type::nDim==3)?mesh_type::nDim-2:mesh_type::nDim-1,mesh_type::nOrder,mesh_type::nRealDim>, value_type, tag > >,
66  mpl::identity< Mesh< Hypercube<(mesh_type::nDim==3)?mesh_type::nDim-2:mesh_type::nDim-1,mesh_type::nOrder,mesh_type::nRealDim>, value_type, tag > > >::type::type mesh_edges_type;
67  typedef boost::shared_ptr<mesh_edges_type> mesh_edges_ptrtype;
68 
69  typedef typename mpl::if_< mpl::equal_to< idim_type ,mpl::size_t<MESH_ELEMENTS> >,
70  mpl::identity<mesh_type>,
71  typename mpl::if_< mpl::equal_to< idim_type ,mpl::size_t<MESH_FACES> >,
72  mpl::identity<mesh_faces_type>,
73  mpl::identity<mesh_edges_type> >::type>::type::type mesh_build_type;
74 
75  typedef boost::shared_ptr<mesh_build_type> mesh_build_ptrtype;
76  typedef SubMeshData smd_type;
77  typedef boost::shared_ptr<smd_type> smd_ptrtype;
78 
79  createSubmeshTool( boost::shared_ptr<MeshType> inputMesh,IteratorRange const& range )
80  :
81  M_mesh( inputMesh ),
82  M_listRange(),
83  M_smd( new smd_type( inputMesh ) )
84  {
85  M_listRange.push_back( range );
86  }
87 
88  createSubmeshTool( boost::shared_ptr<MeshType> inputMesh,std::list<IteratorRange> const& range )
89  :
90  M_mesh( inputMesh ),
91  M_listRange( range ),
92  M_smd( new smd_type( inputMesh ) )
93  {}
94 
95 
96  mesh_build_ptrtype
97  build()
98  {
99  DVLOG(2) << "[createSubmeshTool] extracting mesh\n";
100  return build( mpl::int_<idim_type::value>() );
101  }
102 
103  smd_ptrtype subMeshData() { return M_smd; }
104 
105 private:
106 
107  mesh_ptrtype build( mpl::int_<MESH_ELEMENTS> );
108 
109  mesh_faces_ptrtype build( mpl::int_<MESH_FACES> );
110  mesh_edges_ptrtype build( mpl::int_<MESH_EDGES> );
111 
112  mesh_ptrtype M_mesh;
113  std::list<range_type> M_listRange;
114  smd_ptrtype M_smd;
115 
116 };
117 
118 template <typename MeshType,typename IteratorRange,int TheTag>
119 typename createSubmeshTool<MeshType,IteratorRange,TheTag>::mesh_ptrtype
120 createSubmeshTool<MeshType,IteratorRange,TheTag>::build( mpl::int_<MESH_ELEMENTS> )
121 {
122  typedef typename mesh_type::element_type element_type;
123  typedef typename mesh_type::point_type point_type;
124  typedef typename mesh_type::face_type face_type;
125 
126  mesh_ptrtype newMesh( new mesh_type(M_mesh->worldComm()));
127  //mesh_ptrtype newMesh( new mesh_type );
128 
129  //-----------------------------------------------------------//
130 
131  // inherit the table of markersName
132  BOOST_FOREACH( auto itMark, M_mesh->markerNames() )
133  {
134  newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
135  }
136 
137  //-----------------------------------------------------------//
138 
139  // How the nodes on this mesh will be renumbered to nodes
140  // on the new_mesh.
141  std::map<size_type,size_type> new_node_numbers;
142  std::map<size_type, size_type> new_vertex;
143 
144 
145 
146  // the number of nodes on the new mesh, will be incremented
147  unsigned int n_new_nodes = 0;
148  unsigned int n_new_elem = 0;
149  size_type n_new_faces = 0;
150 
151  const int proc_id = Environment::worldComm().localRank();
152  const int nProc = Environment::worldComm().localSize();
153  std::vector< std::list<boost::tuple<size_type,size_type> > > memory_ghostid( nProc );
154  std::vector< std::vector<size_type> > memory_id( nProc );
155  std::vector< std::vector<size_type> > vecToSend( nProc );
156  std::vector< std::vector<size_type> > vecToRecv( nProc );
157 
158  //-----------------------------------------------------------//
159 
160  std::map<size_type,size_type> new_element_id;
161 
162  std::map<int,std::set<boost::tuple<size_type,size_type> > > ghostCellsFind;
163 
164 
165  auto itListRange = M_listRange.begin();
166  auto const enListRange = M_listRange.end();
167  for ( ; itListRange!=enListRange ; ++itListRange)
168  {
169  iterator_type it, en;
170  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
171 
172 
173  for ( ; it != en; ++ it )
174  {
175  element_type const& old_elem = *it;
176  VLOG(2) << "create sub mesh element from " << it->id() << "\n";google::FlushLogFiles(google::GLOG_INFO);
177  // copy element so that we can modify it
178  element_type new_elem = old_elem;
179 
180  /*
181  // get element markers
182  new_elem.setMarker(old_elem.marker().value());
183  new_elem.setMarker2(old_elem.marker2().value());
184  new_elem.setMarker3(old_elem.marker3().value());
185  */
186  // partitioning update
187  /*new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
188  new_elem.setNumberOfPartitions(old_elem.numberOfPartitions());
189  new_elem.setProcessId(old_elem.processId());
190  //new_elem.setIdInPartition( old_elem.pidInPartition(), n_new_elem );
191  new_elem.setNeighborPartitionIds(old_elem.neighborPartitionIds());// TODO
192  */
193 
194  // reset partitioning data
195  new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
196  new_elem.setNumberOfPartitions( 1 );
197  new_elem.setProcessId( old_elem.processId() );
198  new_elem.idInOthersPartitions().clear();
199  new_elem.neighborPartitionIds().clear();
200 
201 
202  // Loop over the nodes on this element.
203  for ( unsigned int n=0; n < old_elem.nPoints(); n++ )
204  {
205  //FEELPP_ASSERT (old_elem.point( n ).id() < new_node_numbers.size()).error( "invalid point id()" );
206  auto const& old_point = old_elem.point( n );
207 
208  if ( new_node_numbers.find( old_point.id() ) == new_node_numbers.end() )
209  {
210  new_node_numbers[old_point.id()] = n_new_nodes;
211 
212  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) << "\n";
213 
214  point_type pt( old_point );
215  pt.setId( n_new_nodes );
216  pt.setProcessId(old_point.processId());
217 
218  // Add this node to the new mesh
219  newMesh->addPoint ( pt );
220 
221  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() << "\n";
222 
223  // Increment the new node counter
224  n_new_nodes++;
225 
226  if ( n < element_type::numVertices )
227  {
228  CHECK( new_vertex.find(old_point.id()) == new_vertex.end() ) << "already seen this point?";
229  new_vertex[old_point.id()]=1;
230  }
231 
232  if (old_point.numberOfElementsGhost()>0)
233  {
234  auto itg = old_point.elementsGhost().begin();
235  auto const eng = old_point.elementsGhost().end();
236  for ( ; itg!=eng ; ++itg )
237  {
238  auto const procIdGhost = itg->template get<0>();
239  auto const eltIdGhost = itg->template get<1>();
240  auto const& ghostElt = M_mesh->element(eltIdGhost,procIdGhost);
241  ghostCellsFind[procIdGhost].insert(boost::make_tuple( ghostElt.id(),
242  ghostElt.idInOthersPartitions(ghostElt.processId())) );
243  }
244  }
245 
246  }
247 
248  // Define this element's connectivity on the new mesh
249  CHECK ( new_node_numbers[old_point.id()] < newMesh->numPoints() ) << "invalid connectivity";
250 
251  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding point old(" << old_point.id()
252  << ") as point new(" << new_node_numbers[old_point.id()]
253  << ") in element " << new_elem.id() << "\n";
254 
255  new_elem.setPoint( n, newMesh->point( new_node_numbers[old_point.id()] ) );
256 
257  } // for (unsigned int n=0 ... )
258 
259  // set id of element
260  new_elem.setId ( n_new_elem );
261 
262  // increment the new element counter
263  n_new_elem++;
264 
265 
266  // Add an equivalent element type to the new_mesh
267  auto const& e = newMesh->addElement( new_elem );
268  new_element_id[it->id()]= e.id();
269  M_smd->bm.insert( typename smd_type::bm_type::value_type( e.id(), it->id() ) );
270 
271  // Maybe add faces for this element
272  for ( unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
273  {
274  if ( !old_elem.facePtr( s ) ) continue;
275  // only add face on the boundary: they have some data
276  // (boundary ids) which cannot be retrieved otherwise
277  //if ( old_elem.neighbor(s) == invalid_size_type_value )
278  size_type global_face_id = old_elem.face( s ).id();
279 
280  if ( M_mesh->hasFace( global_face_id ) )
281  {
282 
283  // get the corresponding face
284  face_type const& old_face = old_elem.face( s );
285  face_type new_face = old_face;
286 
287  // disconnect from elements of old mesh,
288  // the connection will be redone in
289  // \c updateForUse()
290  new_face.disconnect();
291 
292  //this line is very important!!!!!!!!!!
293  //updateForUse put false for internalfaces
294  new_face.setOnBoundary( true );
295 
296  // update points info
297  for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
298  {
299  new_face.setPoint( p, newMesh->point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).id()] ) );
300  }
301 
302  new_face.setId( n_new_faces++ );
303 
304  // reset partitioning data
305  new_face.setProcessIdInPartition( old_face.pidInPartition() );
306  new_face.setNumberOfPartitions( 1 );
307  new_face.setProcessId( old_face.processId() );
308  new_face.idInOthersPartitions().clear();
309  new_face.neighborPartitionIds().clear();
310 
311  // add it to the list of faces
312  auto addFaceRes = newMesh->addFace( new_face );
313 
314  if (addFaceRes.second)
315  {
316  //update the face
317  newMesh->faces().modify( addFaceRes.first, detail::updateIdInOthersPartitions( addFaceRes.first->pidInPartition(), addFaceRes.first->id() ) );
318  }
319 
320 #if 0
321  // maybe have a ghost cell connected to this face
322  if (false && old_face.isInterProcessDomain())
323  {
324  if (old_face.isConnectedTo0())
325  {
326  auto const& eltConnect0 = old_face.element0();
327  if (eltConnect0.isGhostCell())
328  ghostCellsFind[eltConnect0.processId()].insert(boost::make_tuple(eltConnect0.id(),
329  eltConnect0.idInOthersPartitions(eltConnect0.processId())) );
330  }
331  if (old_face.isConnectedTo1())
332  {
333  auto const& eltConnect1 = old_face.element1();
334  if (eltConnect1.isGhostCell())
335  ghostCellsFind[eltConnect1.processId()].insert(boost::make_tuple(eltConnect1.id(),
336  eltConnect1.idInOthersPartitions(eltConnect1.processId())) );
337  }
338  }
339  else
340  {
341  // TODO see points (2D or 3D) or edges (3D)
342  for ( uint16_type p = 0; p < old_face.nPoints(); ++p )
343  {
344  auto const& oldpoint = old_elem.point( old_elem.fToP( s,p ) );
345  //if ( oldpoint.isLinkedToOtherPartitions() ) std::cout << "\n Find point multiprocess"<<std::endl;
346  //if ( oldpoint.idInPartition().size()>0 ) std::cout << "\n Find point multiprocess"<<std::endl;
347  }
348 
349  }
350 #endif
351  }
352 
353  } // for (unsigned int s=0 ... )
354 #if 0
355  if ( it->isGhostCell() )
356  {
357  VLOG(2) << " - is ghostcell \n";google::FlushLogFiles(google::GLOG_INFO);
358  for (auto it_pid=it->idInOthersPartitions().begin(),en_pid=it->idInOthersPartitions().end() ; it_pid!=en_pid ; ++it_pid)
359  {
360  VLOG(2) << " " << it_pid->first << "-" << it_pid->second << "-"<<it->pidInPartition()<<"-"<<Environment::worldComm().localRank() << "\n";
361  const int procToSend=it_pid->first;
362  if (procToSend!=it->pidInPartition())
363  {
364  memory_ghostid[procToSend].push_back(boost::make_tuple(new_elem.id()/*it->id()*/,it_pid->second));
365  }
366  }
367  }
368 #endif
369 
370  } // for( ; it != en; ++ it )
371  } // for ( ; itListRange!=enListRange ; ++itListRange)
372 
373  if ( nProc > 1 )
374  {
375  auto const theWorldCommSize = newMesh->worldComm().localComm().size();
376  std::vector<int> nbMsgToSend( theWorldCommSize , 0 );
377  std::vector< std::map<int,size_type> > mapMsg( theWorldCommSize );
378 
379  auto itGhostFind = ghostCellsFind.begin();
380  auto const enGhostFind = ghostCellsFind.end();
381  for ( ; itGhostFind!=enGhostFind ; ++itGhostFind )
382  {
383  auto const realProcId = itGhostFind->first;
384  auto itIdElt = itGhostFind->second.begin();
385  auto const enIdElt = itGhostFind->second.end();
386  for ( ; itIdElt!=enIdElt ; ++itIdElt)
387  {
388  auto const idEltInMyProc =itIdElt->template get<0>();
389  auto const idEltInRealProc =itIdElt->template get<1>();
390  newMesh->worldComm().localComm().send(realProcId, nbMsgToSend[realProcId], idEltInRealProc);
391  mapMsg[realProcId].insert( std::make_pair( nbMsgToSend[realProcId],idEltInMyProc ) );
392  ++nbMsgToSend[realProcId];
393  }
394  }
395 
396  // counter of msg received for each process
397  std::vector<int> nbMsgToRecv;
398  mpi::all_to_all( newMesh->worldComm().localComm(),
399  nbMsgToSend,
400  nbMsgToRecv );
401 
402  // recv dof asked and re-send dof in this proc
403  for ( int proc=0; proc<theWorldCommSize; ++proc )
404  {
405  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
406  {
407  //recv
408  size_type idEltRecv;
409  newMesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
410  // search id
411  size_type idEltInNewMesh = invalid_size_type_value;
412  auto const itFindId = new_element_id.find(idEltRecv);
413  if ( itFindId != new_element_id.end() )
414  idEltInNewMesh = itFindId->second;
415  // send response
416  newMesh->worldComm().localComm().send( proc, cpt, idEltInNewMesh );
417  }
418  }
419 
420  // get response to initial request and update Feel::Mesh::Faces data
421  for ( int proc=0; proc<theWorldCommSize; ++proc )
422  {
423  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
424  {
425  size_type idEltAsked;
426  newMesh->worldComm().localComm().recv( proc, cpt, idEltAsked );
427 
428  if (idEltAsked != invalid_size_type_value)
429  {
430  element_type const& old_elem = M_mesh->element( mapMsg[proc][cpt], proc );
431 
432  if (new_element_id.find(old_elem.id())!=new_element_id.end() ) continue;
433  // copy element so that we can modify it
434  element_type new_elem = old_elem;
435 
436  // partitioning update
437  new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
438  new_elem.setNumberOfPartitions(2/*old_elem.numberOfPartitions()*/);
439  new_elem.setProcessId(old_elem.processId());
440  new_elem.idInOthersPartitions().clear();
441  std::vector<int> newNeighborPartitionIds(1);
442  newNeighborPartitionIds[0]=old_elem.processId();
443  new_elem.setNeighborPartitionIds( newNeighborPartitionIds );
444 
445  // Loop over the nodes on this element.
446  for ( unsigned int n=0; n < old_elem.nPoints(); n++ )
447  {
448  //FEELPP_ASSERT (old_elem.point( n ).id() < new_node_numbers.size()).error( "invalid point id()" );
449  if ( new_node_numbers.find( old_elem.point( n ).id() ) == new_node_numbers.end() )
450  {
451  new_node_numbers[old_elem.point( n ).id()] = n_new_nodes;
452 
453  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) << "\n";
454 
455  point_type pt( old_elem.point( n ) );
456  pt.setId( n_new_nodes );
457  pt.setProcessId(invalid_uint16_type_value);//old_elem.point( n ).processId());
458 
459  // Add this node to the new mesh
460  newMesh->addPoint ( pt );
461 
462  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() << "\n";
463 
464  // Increment the new node counter
465  n_new_nodes++;
466 
467  if ( n < element_type::numVertices )
468  {
469  CHECK( new_vertex.find(old_elem.point( n ).id()) == new_vertex.end() ) << "already seen this point?";
470  new_vertex[old_elem.point( n ).id()]=1;
471  }
472  }
473 
474  // Define this element's connectivity on the new mesh
475  CHECK ( new_node_numbers[old_elem.point( n ).id()] < newMesh->numPoints() ) << "invalid connectivity";
476 
477  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding point old(" << old_elem.point( n ).id()
478  << ") as point new(" << new_node_numbers[old_elem.point( n ).id()]
479  << ") in element " << new_elem.id() << "\n";
480 
481  new_elem.setPoint( n, newMesh->point( new_node_numbers[old_elem.point( n ).id()] ) );
482 
483  } // for (unsigned int n=0 ... )
484 
485  // set id of element
486  new_elem.setId ( n_new_elem );
487 
488  // increment the new element counter
489  n_new_elem++;
490  // Add an equivalent element type to the new_mesh
491  auto const& e = newMesh->addElement( new_elem );
492  new_element_id[old_elem.id()]= e.id();
493 
494  M_smd->bm.insert( typename smd_type::bm_type::value_type( e.id(), old_elem.id() ) );
495 
496  //idEltAsked;
497  auto elttt = newMesh->elementIterator( e.id(), proc );
498  //newMesh->elements().modify( elttt, detail::updateIdInOthersPartitions( newMesh->worldComm().localRank(), e.id() ) );
499  newMesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, idEltAsked ) );
500 
501  } // if (idEltAsked != invalid_size_type_value)
502 
503  }
504 
505  }
506 
507  } // if ( nProc > 1 )
508 
509 
510  VLOG(2) << "submesh created\n";google::FlushLogFiles(google::GLOG_INFO);
511  newMesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
512  []( int lhs, std::pair<int,int> const& rhs )
513  {
514  return lhs+rhs.second;
515  } ) );
516 
517  VLOG(2) << "[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";google::FlushLogFiles(google::GLOG_INFO);
518  // Prepare the new_mesh for use
519  newMesh->components().set ( MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
520  newMesh->updateForUse();
521  VLOG(2) << "createSubmesh(MESH_ELEMENTS) done\n";google::FlushLogFiles(google::GLOG_INFO);
522 
523  return newMesh;
524 }
525 
526 template <typename MeshType,typename IteratorRange,int TheTag>
527 typename createSubmeshTool<MeshType,IteratorRange,TheTag>::mesh_faces_ptrtype
528 createSubmeshTool<MeshType,IteratorRange,TheTag>::build( mpl::int_<MESH_FACES> )
529 {
530  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] creating new mesh" << "\n";
531  mesh_faces_ptrtype newMesh( new mesh_faces_type( M_mesh->worldComm()) );
532 
533  //-----------------------------------------------------------//
534  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] extraction mesh faces" << "\n";
535  // inherit the table of markersName
536  BOOST_FOREACH( auto itMark, M_mesh->markerNames() )
537  {
538  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding marker " << itMark.first <<"\n";
539  newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
540  }
541 
542  //-----------------------------------------------------------//
543 
544  typedef typename mesh_faces_type::element_type new_element_type;
545  typedef typename mesh_faces_type::face_type new_face_type;
546 
547  std::map<size_type,size_type> new_node_numbers;
548  std::map<size_type,size_type> new_vertex;
549 
550  // the number of nodes on the new mesh, will be incremented
551  unsigned int n_new_nodes = 0;
552  unsigned int n_new_elem = 0;
553  size_type n_new_faces = 0;
554 
555  const int proc_id = newMesh->worldComm().localRank();
556  const int nProc = newMesh->worldComm().localSize();
557  std::map<size_type,size_type> new_element_id;
558  std::map<int,std::set<boost::tuple<size_type,size_type> > > ghostCellsFind;
559 
560  //-----------------------------------------------------------//
561 
562  auto itListRange = M_listRange.begin();
563  auto const enListRange = M_listRange.end();
564  for ( ; itListRange!=enListRange ; ++itListRange)
565  {
566  iterator_type it, en;
567  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
568 
569  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] extracting " << std::distance(it,en) << " faces " << "\n";
570  for ( ; it != en; ++ it )
571  {
572  // create a new element
573  //element_type const& old_elem = *it;
574  auto const& oldElem = *it;
575  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] + face : " << it->id() << "\n";
576 
577  // copy element so that we can modify it
578  new_element_type newElem;// = oldElem;
579 
580  // get element markers
581  newElem.setMarker( oldElem.marker().value() );
582  newElem.setMarker2( oldElem.marker2().value() );
583  newElem.setMarker3( oldElem.marker3().value() );
584  //???newElem.setOnBoundary( true );
585 
586  // reset partitioning data
587  newElem.setProcessIdInPartition( oldElem.pidInPartition() );
588  newElem.setNumberOfPartitions( 1 );
589  newElem.setProcessId( oldElem.processId() );
590  newElem.idInOthersPartitions().clear();
591  newElem.neighborPartitionIds().clear();
592 
593  // Loop over the nodes on this element.
594  for ( unsigned int n=0; n < oldElem.nPoints(); n++ )
595  {
596  auto const& oldPoint = oldElem.point( n );
597 
598  if ( new_node_numbers.find(oldPoint.id()) == new_node_numbers.end() )
599  {
600  new_node_numbers[oldPoint.id()] = n_new_nodes;
601 
602  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << oldPoint << "\n";
603 
604  typename mesh_faces_type::point_type pt( oldPoint );
605  pt.setId( n_new_nodes );
606  pt.elementsGhost().clear();
607  pt.setProcessId( oldPoint.processId() );
608 
609  // Add this node to the new mesh
610  newMesh->addPoint( pt );
611 
612  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() << "\n";
613 
614  // Increment the new node counter
615  n_new_nodes++;
616 
617  if ( n < new_element_type::numVertices )
618  {
619  CHECK( new_vertex.find(oldPoint.id()) == new_vertex.end() ) << "already seen this point?";
620  new_vertex[oldPoint.id()]=1;
621  }
622 
623  if (oldPoint.numberOfElementsGhost()>0)
624  {
625  auto itg = oldPoint.elementsGhost().begin();
626  auto const eng = oldPoint.elementsGhost().end();
627  for ( ; itg!=eng ; ++itg )
628  {
629  auto const procIdGhost = itg->template get<0>();
630  auto const eltIdGhost = itg->template get<1>();
631  auto const& ghostElt = M_mesh->element(eltIdGhost,procIdGhost);
632  //std::cout << " ghostElt.numTopologicalFaces " << ghostElt.numTopologicalFaces << std::endl;
633  for ( unsigned int s=0; s<ghostElt.numTopologicalFaces; s++ )
634  {
635  auto const& ghostFace = ghostElt.face( s );
636  if (ghostFace.isInterProcessDomain()) continue;
637 #if 1
638  ghostCellsFind[procIdGhost].insert(boost::make_tuple( ghostFace.id(),
639  ghostFace.idInOthersPartitions(ghostElt.processId())) );
640 #else
641  auto itop = ghostFace.idInOthersPartitions().begin();
642  auto const enop = ghostFace.idInOthersPartitions().end();
643  for ( ; itop != enop ; ++itop )
644  ghostCellsFind[itop->first].insert(boost::make_tuple( ghostFace.id(),
645  itop->second ) );//ghostFace.idInOthersPartitions(
646 #endif
647  }
648  }
649  }
650 
651 
652  }
653 
654  newElem.setPoint( n, newMesh->point( new_node_numbers[oldPoint.id()] ) );
655 
656  } // end for n
657 
658  // set id of element
659  newElem.setId ( n_new_elem );
660 
661  // increment the new element counter
662  n_new_elem++;
663 
664  // Add an equivalent element type to the new_mesh
665  auto const& e = newMesh->addElement( newElem );
666  new_element_id[it->id()]= e.id();
667  M_smd->bm.insert( typename smd_type::bm_type::value_type( e.id(), it->id() ) );
668 
669 #if 0
670  // Maybe add faces for this element
671  for ( unsigned int s=0; s<oldElem.numTopologicalFaces; s++ )
672  {
673  if ( !oldElem.facePtr( s ) ) continue;
674  // only add face on the boundary: they have some data
675  // (boundary ids) which cannot be retrieved otherwise
676  //if ( old_elem.neighbor(s) == invalid_size_type_value )
677  size_type global_face_id = oldElem.face( s ).id();
678 
679  if ( M_mesh->hasFace( global_face_id ) )
680  {
681  // get the corresponding face
682  auto const& old_face = oldElem.face( s );
683  new_face_type new_face;// = old_face;
684 
685  // disconnect from elements of old mesh,
686  // the connection will be redone in
687  // \c updateForUse()
688  new_face.disconnect();
689 
690  //this line is very important!!!!!!!!!!
691  //updateForUse put false for internalfaces
692  new_face.setOnBoundary( true );
693 
694  // update points info
695  for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
696  {
697  new_face.setPoint( p, newMesh->point( new_node_numbers[oldElem.point( oldElem.fToP( s,p ) ).id()] ) );
698  }
699 
700  new_face.setId( n_new_faces++ );
701 
702  // reset partitioning data
703  new_face.setProcessIdInPartition( old_face.pidInPartition() );
704  new_face.setNumberOfPartitions( 1 );
705  new_face.setProcessId( old_face.processId() );
706  new_face.idInOthersPartitions().clear();
707  new_face.neighborPartitionIds().clear();
708 
709  // add it to the list of faces
710  auto addFaceRes = newMesh->addFace( new_face );
711 
712  if (addFaceRes.second)
713  {
714  //update the face
715  newMesh->faces().modify( addFaceRes.first, detail::updateIdInOthersPartitions( addFaceRes.first->pidInPartition(), addFaceRes.first->id() ) );
716  }
717  }
718  }
719 #endif
720 
721  } // end for it
722  } // for ( ; itListRange!=enListRange ; ++itListRange)
723 
724 
725  if ( nProc > 1 )
726  {
727  auto const theWorldCommSize = newMesh->worldComm().localComm().size();
728  std::vector<int> nbMsgToSend( theWorldCommSize , 0 );
729  std::vector< std::map<int,size_type> > mapMsg( theWorldCommSize );
730 
731  auto itGhostFind = ghostCellsFind.begin();
732  auto const enGhostFind = ghostCellsFind.end();
733  for ( ; itGhostFind!=enGhostFind ; ++itGhostFind )
734  {
735  auto const realProcId = itGhostFind->first;
736  auto itIdElt = itGhostFind->second.begin();
737  auto const enIdElt = itGhostFind->second.end();
738  for ( ; itIdElt!=enIdElt ; ++itIdElt)
739  {
740  auto const idEltInMyProc =itIdElt->template get<0>();
741  auto const idEltInRealProc =itIdElt->template get<1>();
742  newMesh->worldComm().localComm().send(realProcId, nbMsgToSend[realProcId], idEltInRealProc);
743  mapMsg[realProcId].insert( std::make_pair( nbMsgToSend[realProcId],idEltInMyProc ) );
744  ++nbMsgToSend[realProcId];
745  }
746  }
747 
748  // counter of msg received for each process
749  std::vector<int> nbMsgToRecv;
750  mpi::all_to_all( newMesh->worldComm().localComm(),
751  nbMsgToSend,
752  nbMsgToRecv );
753 
754  // recv dof asked and re-send dof in this proc
755  for ( int proc=0; proc<theWorldCommSize; ++proc )
756  {
757  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
758  {
759  //recv
760  size_type idEltRecv;
761  newMesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
762  // search id
763  size_type idEltInNewMesh = invalid_size_type_value;
764  auto const itFindId = new_element_id.find(idEltRecv);
765  if ( itFindId != new_element_id.end() )
766  idEltInNewMesh = itFindId->second;
767  // send response
768  newMesh->worldComm().localComm().send( proc, cpt, idEltInNewMesh );
769  }
770  }
771 
772  // get response to initial request and update Feel::Mesh::Faces data
773  for ( int proc=0; proc<theWorldCommSize; ++proc )
774  {
775  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
776  {
777  size_type idEltAsked;
778  newMesh->worldComm().localComm().recv( proc, cpt, idEltAsked );
779 
780  if (idEltAsked != invalid_size_type_value)
781  {
782  auto const& old_elem = M_mesh->face( mapMsg[proc][cpt] /*, proc*/ );
783 
784  if (new_element_id.find(old_elem.id())!=new_element_id.end() ) continue;
785  // copy element so that we can modify it
786  new_element_type new_elem; // = old_elem;
787 
788  // partitioning update
789  new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
790  new_elem.setNumberOfPartitions(2/*old_elem.numberOfPartitions()*/);
791  new_elem.setProcessId( proc /*old_elem.processId()*/);
792  new_elem.idInOthersPartitions().clear();
793  std::vector<int> newNeighborPartitionIds(1);
794  newNeighborPartitionIds[0]=old_elem.processId();
795  new_elem.setNeighborPartitionIds( newNeighborPartitionIds );
796 
797  // Loop over the nodes on this element.
798  for ( unsigned int n=0; n < old_elem.nPoints(); n++ )
799  {
800  auto const& oldPoint = old_elem.point( n );
801  //FEELPP_ASSERT (old_elem.point( n ).id() < new_node_numbers.size()).error( "invalid point id()" );
802  if ( new_node_numbers.find(oldPoint.id()) == new_node_numbers.end() )
803  {
804  new_node_numbers[oldPoint.id()] = n_new_nodes;
805 
806  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << oldPoint << "\n";
807 
808  typename mesh_faces_type::point_type pt( oldPoint );
809  pt.setId( n_new_nodes );
810  pt.setProcessId( invalid_uint16_type_value );
811  pt.elementsGhost().clear();
812 
813  // Add this node to the new mesh
814  newMesh->addPoint ( pt );
815 
816  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() << "\n";
817 
818  // Increment the new node counter
819  n_new_nodes++;
820 
821  if ( n < new_element_type::numVertices )
822  {
823  CHECK( new_vertex.find(oldPoint.id()) == new_vertex.end() ) << "already seen this point?";
824  new_vertex[oldPoint.id()]=1;
825  }
826  }
827 
828  // Define this element's connectivity on the new mesh
829  CHECK ( new_node_numbers[old_elem.point( n ).id()] < newMesh->numPoints() ) << "invalid connectivity";
830 
831  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding point old(" << old_elem.point( n ).id()
832  << ") as point new(" << new_node_numbers[old_elem.point( n ).id()]
833  << ") in element " << new_elem.id() << "\n";
834 
835  new_elem.setPoint( n, newMesh->point( new_node_numbers[oldPoint.id()] ) );
836 
837  } // for (unsigned int n=0 ... )
838 
839  // set id of element
840  new_elem.setId ( n_new_elem );
841 
842  // increment the new element counter
843  n_new_elem++;
844  // Add an equivalent element type to the new_mesh
845  auto const& e = newMesh->addElement( new_elem );
846  new_element_id[old_elem.id()]= e.id();
847  M_smd->bm.insert( typename smd_type::bm_type::value_type( e.id(), old_elem.id() ) );
848 
849  // save idEltAsked;
850  auto elttt = newMesh->elementIterator( e.id(), proc );
851  newMesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, idEltAsked ) );
852  } // if (idEltAsked != invalid_size_type_value)
853 
854  }
855  }
856 
857  } // if ( nProc > 1 )
858 
859  newMesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
860  []( int lhs, std::pair<int,int> const& rhs )
861  {
862  return lhs+rhs.second;
863  } ) );
864 
865  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";
866  // Prepare the new_mesh for use
867  newMesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
868  newMesh->updateForUse();
869 
870  return newMesh;
871 }
872 
873 template <typename MeshType,typename IteratorRange,int TheTag>
874 typename createSubmeshTool<MeshType,IteratorRange,TheTag>::mesh_edges_ptrtype
875 createSubmeshTool<MeshType,IteratorRange,TheTag>::build( mpl::int_<MESH_EDGES> )
876 {
877  // we don't deal with this situation yet
878  M_smd.reset();
879 
880  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] creating new mesh" << "\n";
881  mesh_edges_ptrtype newMesh( new mesh_edges_type( M_mesh->worldComm()) );
882  //mesh_edges_ptrtype newMesh( new mesh_edges_type );
883 
884  //-----------------------------------------------------------//
885  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] extraction mesh edges" << "\n";
886  // inherit the table of markersName
887  BOOST_FOREACH( auto itMark, M_mesh->markerNames() )
888  {
889  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] adding marker " << itMark.first <<"\n";
890  newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
891  }
892 
893  //-----------------------------------------------------------//
894 
895  typedef typename mesh_edges_type::element_type new_element_type;
896 
897  std::map<size_type,size_type> new_node_numbers;
898  std::map<size_type,size_type> new_vertex;
899 
900  // the number of nodes on the new mesh, will be incremented
901  unsigned int n_new_nodes = 0;
902  unsigned int n_new_elem = 0;
903  size_type n_new_edges = 0;
904 
905  //-----------------------------------------------------------//
906 
907  auto itListRange = M_listRange.begin();
908  auto const enListRange = M_listRange.end();
909  for ( ; itListRange!=enListRange ; ++itListRange)
910  {
911  iterator_type it, en;
912  boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
913 
914  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] extracting " << std::distance(it,en) << " edges " << "\n";
915  for ( ; it != en; ++ it )
916  {
917  // create a new element
918  //element_type const& old_elem = *it;
919  auto oldElem = *it;
920  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] + face : " << it->id() << "\n";
921 
922  // copy element so that we can modify it
923  new_element_type newElem;// = oldElem;
924 
925  // get element markers
926  newElem.setMarker( oldElem.marker().value() );
927  newElem.setMarker2( oldElem.marker2().value() );
928  newElem.setMarker3( oldElem.marker3().value() );
929 
930  // partitioning update
931  newElem.setProcessIdInPartition( oldElem.pidInPartition() );
932  newElem.setNumberOfPartitions(oldElem.numberOfPartitions());
933  newElem.setProcessId(oldElem.processId());
934  //newElem.setIdInPartition( oldElem.pidInPartition(), n_new_elem );
935  newElem.setNeighborPartitionIds(oldElem.neighborPartitionIds());// TODO
936 
937 
938  DVLOG(2) << "\n oldElem.nPoints " << oldElem.nPoints() << "\n";
939  // Loop over the nodes on this element.
940  for ( unsigned int n=0; n < oldElem.nPoints(); n++ )
941  {
942 
943  if ( new_node_numbers.find(oldElem.point( n ).id()) == new_node_numbers.end() )
944  {
945  new_node_numbers[oldElem.point( n ).id()] = n_new_nodes;
946 
947  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] insert point " << oldElem.point(n) << "\n";
948 
949  typename mesh_edges_type::point_type pt( oldElem.point( n ) );
950  pt.setId( n_new_nodes );
951 
952  // Add this node to the new mesh
953  newMesh->addPoint( pt );
954 
955  DVLOG(2) << "[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() << "\n";
956 
957  // Increment the new node counter
958  n_new_nodes++;
959 
960  if ( n < new_element_type::numVertices )
961  {
962  CHECK( new_vertex.find(oldElem.point( n ).id()) == new_vertex.end() ) << "already seen this point?";
963  new_vertex[oldElem.point( n ).id()]=1;
964  }
965 
966  }
967 
968  newElem.setPoint( n, newMesh->point( new_node_numbers[oldElem.point( n ).id()] ) );
969  newElem.setFace( n, newMesh->point( new_node_numbers[oldElem.point( n ).id()] ) );
970  } // end for n
971  CHECK( newElem.pointPtr(0) ) << "invalid point 0 in edge";
972  CHECK( newElem.pointPtr(1) ) << "invalid point 1 in edge";
973  CHECK( newElem.facePtr(0) ) << "invalid face 0 in edge";
974  CHECK( newElem.facePtr(1) ) << "invalid face 1 in edge";
975 
976  // set id of element
977  newElem.setId ( n_new_elem );
978  newElem.setProcessId ( oldElem.processId() );
979 
980  // increment the new element counter
981  n_new_elem++;
982 
983  // Add an equivalent element type to the new_mesh
984  newMesh->addElement( newElem );
985  } // end for it
986  } // for ( ; itListRange!=enListRange ; ++itListRange)
987 
988 
989  newMesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
990  []( int lhs, std::pair<int,int> const& rhs )
991  {
992  return lhs+rhs.second;
993  } ) );
994 
995  DVLOG(2) << "[createSubmesh] update face/edge info if necessary\n";
996  // Prepare the new_mesh for use
997  newMesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
998  newMesh->updateForUse();
999 
1000 
1001 
1002  return newMesh;
1003 }
1004 
1005 namespace detail
1006 {
1007 template <typename RangeType>
1008 struct submeshrangetype
1009 {
1010  typedef typename mpl::if_< boost::is_std_list<RangeType>,
1011  mpl::identity<RangeType>,
1012  mpl::identity<std::list<RangeType> > >::type::type::value_type type;
1013 };
1014 }
1015 template <typename MeshType,typename IteratorRange, int TheTag = MeshType::tag>
1016 typename createSubmeshTool<MeshType,typename detail::submeshrangetype<IteratorRange>::type/*IteratorRange*/,TheTag>::mesh_build_ptrtype
1017 createSubmesh( boost::shared_ptr<MeshType> inputMesh,IteratorRange const& range, size_type ctx = EXTRACTION_KEEP_MESH_RELATION )
1018 {
1019  //DVLOG(2) << "[createSubmesh] extracting " << range.template get<0>() << " nb elements :"
1020  //<< std::distance(range.template get<1>(),range.template get<2>()) << "\n";
1021 
1022  createSubmeshTool<MeshType,typename detail::submeshrangetype<IteratorRange>::type,TheTag> cSmT( inputMesh,range );
1023  auto m = cSmT.build();
1024  Context c( ctx );
1025  if ( c.test( EXTRACTION_KEEP_MESH_RELATION ) )
1026  m->setSubMeshData( cSmT.subMeshData() );
1027  return m;
1028 }
1029 
1030 
1031 } // namespace Feel
1032 
1033 #endif // createsubmesh

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