30 #ifndef __createsubmesh_H
31 #define __createsubmesh_H 1
33 #include <boost/mpl/if.hpp>
34 #include <boost/mpl/identity.hpp>
35 #include <boost/spirit/home/phoenix/stl/algorithm.hpp>
41 template <
typename C,
typename V,
int T>
class Mesh;
43 template <
typename MeshType,
typename IteratorRange,
int TheTag=MeshType::tag>
44 class createSubmeshTool
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;
52 static const uint16_type tag = TheTag;
54 typedef typename mesh_type::value_type value_type;
56 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
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;
62 typedef boost::shared_ptr<mesh_faces_type> mesh_faces_ptrtype;
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;
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;
75 typedef boost::shared_ptr<mesh_build_type> mesh_build_ptrtype;
77 typedef boost::shared_ptr<smd_type> smd_ptrtype;
79 createSubmeshTool( boost::shared_ptr<MeshType> inputMesh,IteratorRange
const& range )
83 M_smd( new smd_type( inputMesh ) )
85 M_listRange.push_back( range );
88 createSubmeshTool( boost::shared_ptr<MeshType> inputMesh,std::list<IteratorRange>
const& range )
92 M_smd( new smd_type( inputMesh ) )
99 DVLOG(2) <<
"[createSubmeshTool] extracting mesh\n";
100 return build( mpl::int_<idim_type::value>() );
103 smd_ptrtype subMeshData() {
return M_smd; }
107 mesh_ptrtype build( mpl::int_<MESH_ELEMENTS> );
109 mesh_faces_ptrtype build( mpl::int_<MESH_FACES> );
110 mesh_edges_ptrtype build( mpl::int_<MESH_EDGES> );
113 std::list<range_type> M_listRange;
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> )
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;
126 mesh_ptrtype newMesh(
new mesh_type(M_mesh->worldComm()));
132 BOOST_FOREACH(
auto itMark, M_mesh->markerNames() )
134 newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
141 std::map<size_type,size_type> new_node_numbers;
142 std::map<size_type, size_type> new_vertex;
147 unsigned int n_new_nodes = 0;
148 unsigned int n_new_elem = 0;
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 );
160 std::map<size_type,size_type> new_element_id;
162 std::map<int,std::set<boost::tuple<size_type,size_type> > > ghostCellsFind;
165 auto itListRange = M_listRange.begin();
166 auto const enListRange = M_listRange.end();
167 for ( ; itListRange!=enListRange ; ++itListRange)
169 iterator_type it, en;
170 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
173 for ( ; it != en; ++ it )
175 element_type
const& old_elem = *it;
176 VLOG(2) <<
"create sub mesh element from " << it->id() <<
"\n";google::FlushLogFiles(google::GLOG_INFO);
178 element_type new_elem = old_elem;
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();
203 for (
unsigned int n=0; n < old_elem.nPoints(); n++ )
206 auto const& old_point = old_elem.point( n );
208 if ( new_node_numbers.find( old_point.id() ) == new_node_numbers.end() )
210 new_node_numbers[old_point.id()] = n_new_nodes;
212 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) <<
"\n";
214 point_type pt( old_point );
215 pt.setId( n_new_nodes );
216 pt.setProcessId(old_point.processId());
219 newMesh->addPoint ( pt );
221 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() <<
"\n";
226 if ( n < element_type::numVertices )
228 CHECK( new_vertex.find(old_point.id()) == new_vertex.end() ) <<
"already seen this point?";
229 new_vertex[old_point.id()]=1;
232 if (old_point.numberOfElementsGhost()>0)
234 auto itg = old_point.elementsGhost().begin();
235 auto const eng = old_point.elementsGhost().end();
236 for ( ; itg!=eng ; ++itg )
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())) );
249 CHECK ( new_node_numbers[old_point.id()] < newMesh->numPoints() ) <<
"invalid connectivity";
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";
255 new_elem.setPoint( n, newMesh->point( new_node_numbers[old_point.id()] ) );
260 new_elem.setId ( n_new_elem );
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() ) );
272 for (
unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
274 if ( !old_elem.facePtr( s ) )
continue;
278 size_type global_face_id = old_elem.face( s ).id();
280 if ( M_mesh->hasFace( global_face_id ) )
284 face_type
const& old_face = old_elem.face( s );
285 face_type new_face = old_face;
290 new_face.disconnect();
294 new_face.setOnBoundary(
true );
297 for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
299 new_face.setPoint( p, newMesh->point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).
id()] ) );
302 new_face.setId( n_new_faces++ );
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();
312 auto addFaceRes = newMesh->addFace( new_face );
314 if (addFaceRes.second)
317 newMesh->faces().modify( addFaceRes.first, detail::updateIdInOthersPartitions( addFaceRes.first->pidInPartition(), addFaceRes.first->id() ) );
322 if (
false && old_face.isInterProcessDomain())
324 if (old_face.isConnectedTo0())
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())) );
331 if (old_face.isConnectedTo1())
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())) );
342 for ( uint16_type p = 0; p < old_face.nPoints(); ++p )
344 auto const& oldpoint = old_elem.point( old_elem.fToP( s,p ) );
355 if ( it->isGhostCell() )
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)
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())
364 memory_ghostid[procToSend].push_back(boost::make_tuple(new_elem.id(),it_pid->second));
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 );
379 auto itGhostFind = ghostCellsFind.begin();
380 auto const enGhostFind = ghostCellsFind.end();
381 for ( ; itGhostFind!=enGhostFind ; ++itGhostFind )
383 auto const realProcId = itGhostFind->first;
384 auto itIdElt = itGhostFind->second.begin();
385 auto const enIdElt = itGhostFind->second.end();
386 for ( ; itIdElt!=enIdElt ; ++itIdElt)
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];
397 std::vector<int> nbMsgToRecv;
398 mpi::all_to_all( newMesh->worldComm().localComm(),
403 for (
int proc=0; proc<theWorldCommSize; ++proc )
405 for (
int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
409 newMesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
412 auto const itFindId = new_element_id.find(idEltRecv);
413 if ( itFindId != new_element_id.end() )
414 idEltInNewMesh = itFindId->second;
416 newMesh->worldComm().localComm().send( proc, cpt, idEltInNewMesh );
421 for (
int proc=0; proc<theWorldCommSize; ++proc )
423 for (
int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
426 newMesh->worldComm().localComm().recv( proc, cpt, idEltAsked );
430 element_type
const& old_elem = M_mesh->element( mapMsg[proc][cpt], proc );
432 if (new_element_id.find(old_elem.id())!=new_element_id.end() )
continue;
434 element_type new_elem = old_elem;
437 new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
438 new_elem.setNumberOfPartitions(2);
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 );
446 for (
unsigned int n=0; n < old_elem.nPoints(); n++ )
449 if ( new_node_numbers.find( old_elem.point( n ).id() ) == new_node_numbers.end() )
451 new_node_numbers[old_elem.point( n ).id()] = n_new_nodes;
453 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) <<
"\n";
455 point_type pt( old_elem.point( n ) );
456 pt.setId( n_new_nodes );
460 newMesh->addPoint ( pt );
462 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() <<
"\n";
467 if ( n < element_type::numVertices )
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;
475 CHECK ( new_node_numbers[old_elem.point( n ).id()] < newMesh->numPoints() ) <<
"invalid connectivity";
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";
481 new_elem.setPoint( n, newMesh->point( new_node_numbers[old_elem.point( n ).id()] ) );
486 new_elem.setId ( n_new_elem );
491 auto const& e = newMesh->addElement( new_elem );
492 new_element_id[old_elem.id()]= e.id();
494 M_smd->bm.insert(
typename smd_type::bm_type::value_type( e.id(), old_elem.id() ) );
497 auto elttt = newMesh->elementIterator( e.id(), proc );
499 newMesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, idEltAsked ) );
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 )
514 return lhs+rhs.second;
517 VLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";google::FlushLogFiles(google::GLOG_INFO);
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);
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> )
530 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] creating new mesh" <<
"\n";
531 mesh_faces_ptrtype newMesh(
new mesh_faces_type( M_mesh->worldComm()) );
534 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] extraction mesh faces" <<
"\n";
536 BOOST_FOREACH(
auto itMark, M_mesh->markerNames() )
538 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] adding marker " << itMark.first <<
"\n";
539 newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
544 typedef typename mesh_faces_type::element_type new_element_type;
545 typedef typename mesh_faces_type::face_type new_face_type;
547 std::map<size_type,size_type> new_node_numbers;
548 std::map<size_type,size_type> new_vertex;
551 unsigned int n_new_nodes = 0;
552 unsigned int n_new_elem = 0;
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;
562 auto itListRange = M_listRange.begin();
563 auto const enListRange = M_listRange.end();
564 for ( ; itListRange!=enListRange ; ++itListRange)
566 iterator_type it, en;
567 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
569 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] extracting " << std::distance(it,en) <<
" faces " <<
"\n";
570 for ( ; it != en; ++ it )
574 auto const& oldElem = *it;
575 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] + face : " << it->id() <<
"\n";
578 new_element_type newElem;
581 newElem.setMarker( oldElem.marker().value() );
582 newElem.setMarker2( oldElem.marker2().value() );
583 newElem.setMarker3( oldElem.marker3().value() );
587 newElem.setProcessIdInPartition( oldElem.pidInPartition() );
588 newElem.setNumberOfPartitions( 1 );
589 newElem.setProcessId( oldElem.processId() );
590 newElem.idInOthersPartitions().clear();
591 newElem.neighborPartitionIds().clear();
594 for (
unsigned int n=0; n < oldElem.nPoints(); n++ )
596 auto const& oldPoint = oldElem.point( n );
598 if ( new_node_numbers.find(oldPoint.id()) == new_node_numbers.end() )
600 new_node_numbers[oldPoint.id()] = n_new_nodes;
602 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << oldPoint <<
"\n";
604 typename mesh_faces_type::point_type pt( oldPoint );
605 pt.setId( n_new_nodes );
606 pt.elementsGhost().clear();
607 pt.setProcessId( oldPoint.processId() );
610 newMesh->addPoint( pt );
612 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() <<
"\n";
617 if ( n < new_element_type::numVertices )
619 CHECK( new_vertex.find(oldPoint.id()) == new_vertex.end() ) <<
"already seen this point?";
620 new_vertex[oldPoint.id()]=1;
623 if (oldPoint.numberOfElementsGhost()>0)
625 auto itg = oldPoint.elementsGhost().begin();
626 auto const eng = oldPoint.elementsGhost().end();
627 for ( ; itg!=eng ; ++itg )
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);
633 for (
unsigned int s=0; s<ghostElt.numTopologicalFaces; s++ )
635 auto const& ghostFace = ghostElt.face( s );
636 if (ghostFace.isInterProcessDomain())
continue;
638 ghostCellsFind[procIdGhost].insert(boost::make_tuple( ghostFace.id(),
639 ghostFace.idInOthersPartitions(ghostElt.processId())) );
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(),
654 newElem.setPoint( n, newMesh->point( new_node_numbers[oldPoint.id()] ) );
659 newElem.setId ( n_new_elem );
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() ) );
671 for (
unsigned int s=0; s<oldElem.numTopologicalFaces; s++ )
673 if ( !oldElem.facePtr( s ) )
continue;
677 size_type global_face_id = oldElem.face( s ).id();
679 if ( M_mesh->hasFace( global_face_id ) )
682 auto const& old_face = oldElem.face( s );
683 new_face_type new_face;
688 new_face.disconnect();
692 new_face.setOnBoundary(
true );
695 for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
697 new_face.setPoint( p, newMesh->point( new_node_numbers[oldElem.point( oldElem.fToP( s,p ) ).
id()] ) );
700 new_face.setId( n_new_faces++ );
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();
710 auto addFaceRes = newMesh->addFace( new_face );
712 if (addFaceRes.second)
715 newMesh->faces().modify( addFaceRes.first, detail::updateIdInOthersPartitions( addFaceRes.first->pidInPartition(), addFaceRes.first->id() ) );
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 );
731 auto itGhostFind = ghostCellsFind.begin();
732 auto const enGhostFind = ghostCellsFind.end();
733 for ( ; itGhostFind!=enGhostFind ; ++itGhostFind )
735 auto const realProcId = itGhostFind->first;
736 auto itIdElt = itGhostFind->second.begin();
737 auto const enIdElt = itGhostFind->second.end();
738 for ( ; itIdElt!=enIdElt ; ++itIdElt)
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];
749 std::vector<int> nbMsgToRecv;
750 mpi::all_to_all( newMesh->worldComm().localComm(),
755 for (
int proc=0; proc<theWorldCommSize; ++proc )
757 for (
int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
761 newMesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
764 auto const itFindId = new_element_id.find(idEltRecv);
765 if ( itFindId != new_element_id.end() )
766 idEltInNewMesh = itFindId->second;
768 newMesh->worldComm().localComm().send( proc, cpt, idEltInNewMesh );
773 for (
int proc=0; proc<theWorldCommSize; ++proc )
775 for (
int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
778 newMesh->worldComm().localComm().recv( proc, cpt, idEltAsked );
782 auto const& old_elem = M_mesh->face( mapMsg[proc][cpt] );
784 if (new_element_id.find(old_elem.id())!=new_element_id.end() )
continue;
786 new_element_type new_elem;
789 new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
790 new_elem.setNumberOfPartitions(2);
791 new_elem.setProcessId( proc );
792 new_elem.idInOthersPartitions().clear();
793 std::vector<int> newNeighborPartitionIds(1);
794 newNeighborPartitionIds[0]=old_elem.processId();
795 new_elem.setNeighborPartitionIds( newNeighborPartitionIds );
798 for (
unsigned int n=0; n < old_elem.nPoints(); n++ )
800 auto const& oldPoint = old_elem.point( n );
802 if ( new_node_numbers.find(oldPoint.id()) == new_node_numbers.end() )
804 new_node_numbers[oldPoint.id()] = n_new_nodes;
806 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << oldPoint <<
"\n";
808 typename mesh_faces_type::point_type pt( oldPoint );
809 pt.setId( n_new_nodes );
811 pt.elementsGhost().clear();
814 newMesh->addPoint ( pt );
816 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() <<
"\n";
821 if ( n < new_element_type::numVertices )
823 CHECK( new_vertex.find(oldPoint.id()) == new_vertex.end() ) <<
"already seen this point?";
824 new_vertex[oldPoint.id()]=1;
829 CHECK ( new_node_numbers[old_elem.point( n ).id()] < newMesh->numPoints() ) <<
"invalid connectivity";
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";
835 new_elem.setPoint( n, newMesh->point( new_node_numbers[oldPoint.id()] ) );
840 new_elem.setId ( n_new_elem );
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() ) );
850 auto elttt = newMesh->elementIterator( e.id(), proc );
851 newMesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, idEltAsked ) );
859 newMesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
860 [](
int lhs, std::pair<int,int>
const& rhs )
862 return lhs+rhs.second;
865 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";
867 newMesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
868 newMesh->updateForUse();
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> )
880 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] creating new mesh" <<
"\n";
881 mesh_edges_ptrtype newMesh(
new mesh_edges_type( M_mesh->worldComm()) );
885 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] extraction mesh edges" <<
"\n";
887 BOOST_FOREACH(
auto itMark, M_mesh->markerNames() )
889 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] adding marker " << itMark.first <<
"\n";
890 newMesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
895 typedef typename mesh_edges_type::element_type new_element_type;
897 std::map<size_type,size_type> new_node_numbers;
898 std::map<size_type,size_type> new_vertex;
901 unsigned int n_new_nodes = 0;
902 unsigned int n_new_elem = 0;
907 auto itListRange = M_listRange.begin();
908 auto const enListRange = M_listRange.end();
909 for ( ; itListRange!=enListRange ; ++itListRange)
911 iterator_type it, en;
912 boost::tie( boost::tuples::ignore, it, en ) = *itListRange;
914 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] extracting " << std::distance(it,en) <<
" edges " <<
"\n";
915 for ( ; it != en; ++ it )
920 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] + face : " << it->id() <<
"\n";
923 new_element_type newElem;
926 newElem.setMarker( oldElem.marker().value() );
927 newElem.setMarker2( oldElem.marker2().value() );
928 newElem.setMarker3( oldElem.marker3().value() );
931 newElem.setProcessIdInPartition( oldElem.pidInPartition() );
932 newElem.setNumberOfPartitions(oldElem.numberOfPartitions());
933 newElem.setProcessId(oldElem.processId());
935 newElem.setNeighborPartitionIds(oldElem.neighborPartitionIds());
938 DVLOG(2) <<
"\n oldElem.nPoints " << oldElem.nPoints() <<
"\n";
940 for (
unsigned int n=0; n < oldElem.nPoints(); n++ )
943 if ( new_node_numbers.find(oldElem.point( n ).id()) == new_node_numbers.end() )
945 new_node_numbers[oldElem.point( n ).id()] = n_new_nodes;
947 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << oldElem.point(n) <<
"\n";
949 typename mesh_edges_type::point_type pt( oldElem.point( n ) );
950 pt.setId( n_new_nodes );
953 newMesh->addPoint( pt );
955 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << newMesh->numPoints() <<
"\n";
960 if ( n < new_element_type::numVertices )
962 CHECK( new_vertex.find(oldElem.point( n ).id()) == new_vertex.end() ) <<
"already seen this point?";
963 new_vertex[oldElem.point( n ).id()]=1;
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()] ) );
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";
977 newElem.setId ( n_new_elem );
978 newElem.setProcessId ( oldElem.processId() );
984 newMesh->addElement( newElem );
989 newMesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
990 [](
int lhs, std::pair<int,int>
const& rhs )
992 return lhs+rhs.second;
995 DVLOG(2) <<
"[createSubmesh] update face/edge info if necessary\n";
997 newMesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
998 newMesh->updateForUse();
1007 template <
typename RangeType>
1008 struct submeshrangetype
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;
1015 template <
typename MeshType,
typename IteratorRange,
int TheTag = MeshType::tag>
1016 typename createSubmeshTool<MeshType,typename detail::submeshrangetype<IteratorRange>::type,TheTag>::mesh_build_ptrtype
1017 createSubmesh( boost::shared_ptr<MeshType> inputMesh,IteratorRange
const& range,
size_type ctx = EXTRACTION_KEEP_MESH_RELATION )
1022 createSubmeshTool<MeshType,typename detail::submeshrangetype<IteratorRange>::type,TheTag> cSmT( inputMesh,range );
1023 auto m = cSmT.build();
1025 if ( c.test( EXTRACTION_KEEP_MESH_RELATION ) )
1026 m->setSubMeshData( cSmT.subMeshData() );
1033 #endif // createsubmesh
simplex of dimension Dim
Definition: simplex.hpp:73
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
unifying mesh class
Definition: createsubmesh.hpp:41
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
data structure storing sub mesh data
Definition: submeshdata.hpp:45