33 #include <boost/unordered_map.hpp>
35 #include <boost/foreach.hpp>
36 #include <boost/signals2/signal.hpp>
37 #include <boost/serialization/vector.hpp>
38 #include <boost/serialization/array.hpp>
39 #include <boost/serialization/base_object.hpp>
40 #include <boost/archive/text_oarchive.hpp>
41 #include <boost/archive/text_iarchive.hpp>
42 #include <boost/archive/binary_oarchive.hpp>
43 #include <boost/archive/binary_iarchive.hpp>
45 #include <feel/feelcore/context.hpp>
52 #include <feel/feelmesh/meshutil.hpp>
54 #include <feel/feelmesh/enums.hpp>
56 #include <feel/feelpoly/geomap.hpp>
63 #include <boost/preprocessor/comparison/less.hpp>
64 #include <boost/preprocessor/logical/and.hpp>
65 #include <boost/preprocessor/control/if.hpp>
66 #include <boost/preprocessor/list/at.hpp>
67 #include <boost/preprocessor/list/cat.hpp>
68 #include <boost/preprocessor/list/for_each_product.hpp>
69 #include <boost/preprocessor/logical/or.hpp>
70 #include <boost/preprocessor/tuple/to_list.hpp>
71 #include <boost/preprocessor/tuple/eat.hpp>
72 #include <boost/preprocessor/facilities/empty.hpp>
73 #include <boost/preprocessor/punctuation/comma.hpp>
74 #include <boost/preprocessor/facilities/identity.hpp>
76 #include <boost/enable_shared_from_this.hpp>
80 const size_type EXTRACTION_KEEP_POINTS_IDS = ( 1<<0 );
81 const size_type EXTRACTION_KEEP_EDGES_IDS = ( 1<<1 );
82 const size_type EXTRACTION_KEEP_FACES_IDS = ( 1<<2 );
83 const size_type EXTRACTION_KEEP_VOLUMES_IDS = ( 1<<3 );
84 const size_type EXTRACTION_KEEP_ALL_IDS = ( EXTRACTION_KEEP_POINTS_IDS |
85 EXTRACTION_KEEP_EDGES_IDS |
86 EXTRACTION_KEEP_FACES_IDS |
87 EXTRACTION_KEEP_VOLUMES_IDS );
88 const size_type EXTRACTION_KEEP_MESH_RELATION = ( 1<<4 );
102 std::vector<MeshMarkerName> markerMap(
int Dim );
119 template <
typename GeoShape,
typename T =
double,
int Tag = 0>
122 public mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<0> >,
123 mpl::identity<Mesh0D<GeoShape > >,
124 typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<1> >,
125 mpl::identity<Mesh1D<GeoShape > >,
126 typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<2> >,
127 mpl::identity<Mesh2D<GeoShape> >,
128 mpl::identity<Mesh3D<GeoShape> > >::type>::type>::type::type,
129 public boost::addable<Mesh<GeoShape,T,Tag> >,
130 public boost::enable_shared_from_this< Mesh<GeoShape,T,Tag> >
132 typedef typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<0> >,
133 mpl::identity<Mesh0D<GeoShape> >,
134 typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<1> >,
135 mpl::identity<Mesh1D<GeoShape> >,
136 typename mpl::if_<mpl::equal_to<mpl::int_<GeoShape::nDim>,mpl::int_<2> >,
137 mpl::identity<Mesh2D<GeoShape> >,
138 mpl::identity<Mesh3D<GeoShape> > >::type>::type>::type::type super;
146 static const uint16_type nDim = GeoShape::nDim;
147 static const uint16_type nRealDim = GeoShape::nRealDim;
148 static const uint16_type Shape = GeoShape::Shape;
149 static const uint16_type nOrder = GeoShape::nOrder;
150 static const uint16_type tag = Tag;
156 typedef T value_type;
157 typedef GeoShape shape_type;
158 typedef typename super::return_type return_type;
161 typedef typename super::super_elements super_elements;
162 typedef typename super::elements_type elements_type;
163 typedef typename super::element_type element_type;
164 typedef typename super::element_iterator element_iterator;
165 typedef typename super::element_const_iterator element_const_iterator;
167 typedef typename super::super_faces super_faces;
168 typedef typename super::faces_type faces_type;
169 typedef typename super::face_type face_type;
170 typedef typename super::face_iterator face_iterator;
171 typedef typename super::face_const_iterator face_const_iterator;
173 typedef typename super::edge_type edge_type;
175 typedef typename super::points_type points_type;
176 typedef typename super::point_type point_type;
177 typedef typename super::point_iterator point_iterator;
178 typedef typename super::point_const_iterator point_const_iterator;
180 typedef typename element_type::gm_type gm_type;
181 typedef boost::shared_ptr<gm_type> gm_ptrtype;
183 typedef typename element_type::gm1_type gm1_type;
184 typedef boost::shared_ptr<gm1_type> gm1_ptrtype;
186 template<
size_type ContextID>
190 typedef boost::shared_ptr<type> ptrtype;
194 typedef self_type mesh_type;
195 typedef boost::shared_ptr<self_type> self_ptrtype;
196 typedef self_ptrtype mesh_ptrtype;
198 typedef typename element_type::template reference_convex<T>::type reference_convex_type;
203 typedef typename super::face_processor_type face_processor_type;
204 typedef typename super::face_processor_type element_edge_type;
206 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
207 mpl::identity< Mesh< Simplex< GeoShape::nDim,1,GeoShape::nRealDim>, value_type, Tag > >,
208 mpl::identity< Mesh< Hypercube<GeoShape::nDim,1,GeoShape::nRealDim>,value_type, Tag > > >::type::type P1_mesh_type;
210 typedef boost::shared_ptr<P1_mesh_type> P1_mesh_ptrtype;
215 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
216 mpl::identity<
Mesh<
Simplex< GeoShape::nDim-1,nOrder,GeoShape::nRealDim>, value_type, TheTag > >,
217 mpl::identity<
Mesh< Hypercube<GeoShape::nDim-1,nOrder,GeoShape::nRealDim>,value_type, TheTag > > >::type::type type;
218 typedef boost::shared_ptr<type> ptrtype;
219 typedef boost::shared_ptr<const type> const_ptrtype;
221 typedef typename trace_mesh<Tag>::type trace_mesh_type;
222 typedef typename trace_mesh<Tag>::ptrtype trace_mesh_ptrtype;
226 struct trace_trace_mesh
228 static const uint16_type nDim = (GeoShape::nDim==1)?GeoShape::nDim-1:GeoShape::nDim-2;
229 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
230 mpl::identity< Mesh< Simplex<nDim,nOrder,GeoShape::nRealDim>, value_type, TheTag > >,
231 mpl::identity< Mesh< Hypercube<nDim,nOrder,GeoShape::nRealDim>,value_type, TheTag > > >::type::type type;
233 typedef boost::shared_ptr<type> ptrtype;
234 typedef boost::shared_ptr<const type> const_ptrtype;
236 typedef typename trace_trace_mesh<Tag>::type trace_trace_mesh_type;
237 typedef typename trace_trace_mesh<Tag>::ptrtype trace_trace_mesh_ptrtype;
245 Mesh( WorldComm
const& worldComm = Environment::worldComm() );
251 M_tool_localization.reset();
262 self_type& operator+=( self_type
const& m );
275 return std::distance( this->beginElementWithProcessId( this->worldComm().rank() ),
276 this->endElementWithProcessId( this->worldComm().rank() ) );
284 int ne = std::distance( this->beginElementWithProcessId( this->worldComm().rank() ),
285 this->endElementWithProcessId( this->worldComm().rank() ) );
287 mpi::all_reduce( this->worldComm(), ne, gne, [] (
int x,
int y )
304 gm_ptrtype
const&
gm()
const
312 gm1_ptrtype
const&
gm1()
const
339 return reference_convex_type();
348 return M_e2f.find(std::make_pair(e.id(),n))->second;
357 return M_e2f.find(std::make_pair(e,n))->second;
363 WorldComm
const& worldComm()
const
368 void setWorldComm( WorldComm
const& _worldComm )
370 M_worldComm = _worldComm;
373 mpi::communicator
const& comm()
const
375 return M_worldComm.localComm();
398 return M_e2f[std::make_pair(e.id(),n)];
407 return M_e2f[std::make_pair(e,n)];
415 auto mit = M_markername.find( marker );
416 if ( mit != M_markername.end() )
417 return mit->second[0];
425 auto mit = M_markername.find( marker );
426 if ( mit != M_markername.end() )
427 return mit->second[1];
434 std::map<std::string, std::vector<size_type> >
markerNames()
const
443 boost::shared_ptr<typename self_type::Localization> tool_localization()
445 return M_tool_localization;
473 void addMarkerName( std::pair<std::string, std::vector<size_type> >
const& marker )
475 M_markername.insert( marker );
483 std::vector<size_type> data(2);
486 M_markername[__name]=data;
489 flag_type markerId( boost::any
const& marker );
501 element_iterator
eraseElement( element_iterator position,
bool modify=
true );
504 template<
typename RangeT,
int TheTag>
505 typename trace_mesh<TheTag>::ptrtype
506 trace( RangeT
const& range, mpl::int_<TheTag> )
const;
520 template<
int TheTag=Tag>
521 typename trace_mesh<TheTag>::ptrtype
527 template<
typename RangeT>
528 typename trace_mesh<Tag>::ptrtype
529 trace( RangeT
const& range )
const;
532 template<
int TheTag=Tag>
533 typename trace_mesh<TheTag>::ptrtype
536 return wireBasket(
boundaryfaces(this->shared_from_this()),mpl::int_<TheTag>());
539 template<
typename RangeT,
int TheTag>
540 typename trace_trace_mesh<TheTag>::ptrtype
541 wireBasket( RangeT
const& range, mpl::int_<TheTag> )
const;
543 template<
typename RangeT>
544 typename trace_trace_mesh<Tag>::ptrtype
545 wireBasket( RangeT
const& range )
const;
548 template<
typename Iterator>
549 void createSubmesh( self_type& mesh,
550 Iterator
const& begin_elt,
551 Iterator
const& end_elt,
552 size_type extraction_policies = EXTRACTION_KEEP_ALL_IDS )
const;
563 this->createSubmesh( mesh,
564 this->beginElementWithProcessId( pid ),
565 this->endElementWithProcessId( pid ) );
578 template<
typename IteratorRange>
581 const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
588 template<
typename IteratorRange>
591 const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
592 this->updateMarker2WithRangeElements( range,flag );
598 template<
typename IteratorRange>
601 this->updateMarker2WithRangeFaces( range,flag );
609 template<
typename IteratorRange>
612 const size_type iDim = boost::tuples::template element<0, IteratorRange>::type::value;
619 template<
typename IteratorRange>
622 this->updateMarker3WithRangeElements( range,flag );
628 template<
typename IteratorRange>
631 this->updateMarker3WithRangeFaces( range,flag );
637 void partition (
const uint16_type n_parts = 1 );
649 renumber( mpl::bool_<( nDim > 1 )>() );
651 void renumber( std::vector<size_type>
const& node_map, mpl::int_<1> );
652 void renumber( std::vector<size_type>
const& node_map, mpl::int_<2> );
653 void renumber( std::vector<size_type>
const& node_map, mpl::int_<3> );
672 void send(
int p,
int tag );
677 void recv(
int p,
int tag );
693 BOOST_PARAMETER_MEMBER_FUNCTION( (
void ),
697 ( name,(std::string) )
700 ( type,( std::string ),std::string(
"binary" ) )
701 ( suffix,( std::string ),std::string(
"" ) )
702 ( sep,( std::string ),std::string(
"" ) )
705 Feel::detail::ignore_unused_variable_warning( args );
707 if ( !fs::exists( fs::path( path ) ) )
709 fs::create_directories( fs::path( path ) );
712 std::ostringstream os1;
713 os1 << name << sep << suffix <<
"-" << this->worldComm().globalSize() <<
"." << this->worldComm().globalRank() <<
".fdb";
714 fs::path p = fs::path( path ) / os1.str();
715 fs::ofstream ofs( p );
717 if ( type ==
"binary" )
719 boost::archive::binary_oarchive oa( ofs );
723 else if ( type ==
"text" )
725 boost::archive::text_oarchive oa( ofs );
729 else if ( type ==
"xml" )
735 BOOST_PARAMETER_MEMBER_FUNCTION(
740 ( name,(std::string) )
743 ( update,(
size_type), MESH_CHECK|MESH_UPDATE_EDGES|MESH_UPDATE_FACES )
744 ( type,( std::string ),std::string(
"binary" ) )
745 ( suffix,( std::string ),std::string(
"" ) )
746 ( sep,( std::string ),std::string(
"" ) )
750 Feel::detail::ignore_unused_variable_warning( args );
751 std::ostringstream os1;
752 os1 << name << sep << suffix <<
"-" << this->worldComm().globalSize() <<
"." << this->worldComm().globalRank() <<
".fdb";
753 fs::path p = fs::path( path ) / os1.str();
754 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
755 std::cout <<
"try loading " << p.native() <<
"\n";
756 if ( !fs::exists( p ) )
758 LOG(INFO) <<
"[mesh::load] failed loading " << p.native() <<
"\n";
759 std::ostringstream os2;
760 os2 << name << sep << suffix <<
"-" << this->worldComm().globalSize() <<
"." << this->worldComm().globalRank();
761 p = fs::path( path ) / os2.str();
762 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
763 std::cout <<
" now try loading " << p.native() <<
"\n";
765 if ( !fs::exists( p ) )
767 LOG(INFO) <<
"[mesh::load] failed loading " << p.native() <<
"\n";
772 if ( !fs::is_regular_file( p ) )
774 LOG(INFO) <<
"[mesh::load] failed loading " << p.native() <<
"\n";
778 fs::ifstream ifs( p );
780 if ( type ==
"binary" )
782 boost::archive::binary_iarchive ia( ifs );
786 else if ( type ==
"text" )
788 boost::archive::text_iarchive ia( ifs );
792 else if ( type ==
"xml" )
799 this->components().reset();
800 this->components().set( update );
806 this->components().reset();
847 for(
auto fit = this->beginFace(), fen = this->endFace(); fit != fen; ++fit )
849 if ( fit->isOnBoundary() && !fit->isConnectedTo0() )
851 std::cout <<
"erase boundary face...\n";
852 this->eraseFace( fit );
854 if ( !fit->isOnBoundary() && fit->isConnectedTo0() && !fit->isConnectedTo1() )
856 std::cout <<
"found boundary face...\n";
857 this->
faces().modify( fit, []( face_type& f ){ f.setOnBoundary(
true ); } );
860 this->setUpdatedForUse(
false );
865 void propagateMarkers( mpl::int_<1> ) {}
866 void propagateMarkers( mpl::int_<2> );
867 void propagateMarkers( mpl::int_<3> );
869 friend class boost::serialization::access;
870 template<
class Archive>
871 void serialize( Archive & ar,
const unsigned int version )
873 if ( Archive::is_saving::value )
875 DVLOG(2) <<
"Serializing mesh(saving) ...\n";
876 DVLOG(2) <<
"encoding...\n";
878 DVLOG(2) <<
"loading markers...\n";
880 DVLOG(2) <<
"loading pts...\n";
882 DVLOG(2) <<
"loading faces...\n";
884 DVLOG(2) <<
"loading elts...\n";
887 if ( Archive::is_loading::value )
889 DVLOG(2) <<
"Serializing mesh(loading) ...\n";
890 DVLOG(2) <<
"loading markers...\n";
892 DVLOG(2) <<
"loading pts...\n";
894 DVLOG(2) <<
"loading faces...\n";
896 DVLOG(2) <<
"loading elts...\n";
906 public mpl::if_<mpl::bool_<GeoShape::is_simplex>,
907 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
908 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type
910 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
911 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
912 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type super;
913 typedef typename super::gic_type gic_type;
915 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
916 mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Simplex> >,
917 mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Hypercube> > >::type::type super1;
918 typedef typename super1::gic_type gic1_type;
920 Inverse( boost::shared_ptr<self_type>
const& m )
931 return M_pts_cvx[i].size();
933 void pointsInConvex(
size_type i, std::vector<boost::tuple<size_type, uint16_type > > &itab )
const
935 itab.resize( M_pts_cvx[i].size() );
938 for ( map_iterator it = M_pts_cvx[i].begin(); it != M_pts_cvx[i].end(); ++it )
939 itab[j++] = boost::make_tuple( it->first, it->second );
942 const boost::unordered_map<size_type,node_type> &referenceCoords(
void )
955 void distribute(
bool extrapolation =
false );
960 boost::shared_ptr<self_type> M_mesh;
961 std::vector<std::map<size_type,uint16_type > > M_pts_cvx;
962 typedef typename std::map<size_type, uint16_type >::const_iterator map_iterator;
964 boost::unordered_map<size_type,node_type> M_ref_coords;
965 boost::unordered_map<size_type,double> M_dist;
966 boost::unordered_map<size_type,size_type> M_cvx_pts;
973 typedef Localization localization_type;
974 typedef boost::shared_ptr<localization_type> localization_ptrtype;
976 typedef typename matrix_node<typename node_type::value_type>::type matrix_node_type;
978 typedef KDTree kdtree_type;
979 typedef typename boost::shared_ptr<KDTree> kdtree_ptrtype;
982 typedef boost::tuple<node_type, std::list<size_type> > node_elem_type;
984 typedef typename std::list<boost::tuple<size_type,node_type> > container_output_type;
985 typedef typename std::list<boost::tuple<size_type,node_type> >::iterator container_output_iterator_type;
988 typedef std::map<size_type, container_output_type > container_search_type;
989 typedef typename container_search_type::const_iterator container_search_const_iterator_type;
990 typedef typename container_search_type::iterator container_search_iterator_type;
993 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
994 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Simplex> >,
995 mpl::identity<GeoMapInverse<nDim,nOrder,nRealDim,T,Hypercube> > >::type::type gm_inverse_type;
996 typedef typename gm_inverse_type::gic_type gmc_inverse_type;
998 typedef typename mpl::if_<mpl::bool_<GeoShape::is_simplex>,
999 mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Simplex> >,
1000 mpl::identity<GeoMapInverse<nDim,1,nRealDim,T,Hypercube> > >::type::type gm1_inverse_type;
1001 typedef typename gm1_inverse_type::gic_type gmc1_inverse_type;
1004 typedef typename self_type::gm_type::reference_convex_type ref_convex_type;
1005 typedef typename self_type::gm1_type::reference_convex_type ref_convex1_type;
1012 M_kd_tree( new kdtree_type() ),
1014 IsInitBoundaryFaces( false ),
1015 M_doExtrapolation( true )
1017 DVLOG(2) <<
"[Mesh::Localization] create Localization tool\n";
1018 M_kd_tree->nbNearNeighbor( 15 );
1019 M_resultAnalysis.clear();
1020 DVLOG(2) <<
"[Mesh::Localization] create Localization tool done\n";
1023 Localization( boost::shared_ptr<self_type> m,
bool init_b =
true ) :
1026 IsInitBoundaryFaces( false ),
1027 M_doExtrapolation( true )
1032 M_kd_tree->nbNearNeighbor( 15 );
1033 M_resultAnalysis.clear();
1036 Localization( Localization
const & L ) :
1038 M_kd_tree( new kdtree_type( *( L.M_kd_tree ) ) ),
1039 M_geoGlob_Elts( L.M_geoGlob_Elts ),
1041 IsInitBoundaryFaces( L.IsInitBoundaryFaces ),
1042 M_resultAnalysis( L.M_resultAnalysis ),
1043 M_doExtrapolation( L.M_doExtrapolation )
1050 setMesh( boost::shared_ptr<self_type> m,
bool b=
true )
1059 M_resultAnalysis.clear();
1066 setExtrapolation(
bool b )
1068 M_doExtrapolation = b;
1076 if ( IsInit==
false )
1083 void updateForUseBoundaryFaces()
1085 if ( IsInitBoundaryFaces==
false )
1086 initBoundaryFaces();
1097 bool isInitBoundaryFaces()
const
1099 return IsInitBoundaryFaces;
1102 bool doExtrapolation()
const
1104 return M_doExtrapolation;
1107 boost::shared_ptr<self_type> mesh()
1111 boost::shared_ptr<self_type>
const& mesh()
const
1116 kdtree_ptrtype kdtree()
1121 kdtree_ptrtype
const& kdtree()
const
1126 node_type barycenter()
const;
1127 node_type barycenter(mpl::int_<1> )
const;
1128 node_type barycenter(mpl::int_<2> )
const;
1129 node_type barycenter(mpl::int_<3> )
const;
1131 container_search_type
const & result_analysis()
const {
return M_resultAnalysis;}
1133 container_search_iterator_type result_analysis_begin()
1135 return M_resultAnalysis.begin();
1137 container_search_iterator_type result_analysis_end()
1139 return M_resultAnalysis.end();
1145 boost::tuple<bool,node_type,double> isIn(
size_type _id,
const node_type & _pt )
const;
1146 boost::tuple<bool,node_type,double> isIn(
size_type _id,
const node_type & _pt,
const matrix_node_type & setPoints, mpl::int_<1> )
const;
1147 boost::tuple<uint16_type,std::vector<bool> > isIn( std::vector<size_type> _ids,
const node_type & _pt );
1152 boost::tuple<bool, size_type,node_type> searchElement(
const node_type & p);
1156 boost::tuple<bool, size_type,node_type> searchElement(
const node_type & p,
1157 const matrix_node_type & setPoints,
1160 return searchElement( p );
1166 boost::tuple<bool, size_type,node_type> searchElement(
const node_type & p,
1167 const matrix_node_type & setPoints,
1173 boost::tuple<bool, std::list<boost::tuple<size_type,node_type> > > searchElements(
const node_type & p );
1179 boost::tuple<std::vector<bool>,
size_type> run_analysis(
const matrix_node_type & m,
1186 boost::tuple<std::vector<bool>,
size_type> run_analysis(
const matrix_node_type & m,
1188 const matrix_node_type & setPoints,
1191 return run_analysis( m,eltHypothetical );
1198 boost::tuple<std::vector<bool>,
size_type> run_analysis(
const matrix_node_type & m,
1200 const matrix_node_type & setPoints,
1212 void resetBoundaryFaces()
1214 IsInitBoundaryFaces=
false;
1215 initBoundaryFaces();
1228 void initBoundaryFaces();
1233 void searchInKdTree(
const node_type & p,
1234 std::list< std::pair<size_type, uint> > & listTri );
1238 boost::shared_ptr<self_type> M_mesh;
1240 kdtree_ptrtype M_kd_tree;
1242 std::map<size_type, node_elem_type > M_geoGlob_Elts;
1243 bool IsInit,IsInitBoundaryFaces;
1244 container_search_type M_resultAnalysis;
1245 bool M_doExtrapolation;
1247 ref_convex_type M_refelem;
1248 ref_convex1_type M_refelem1;
1262 template<
typename Observer>
1263 void addObserver( Observer& obs )
1268 void removeFacesFromBoundary( std::initializer_list<uint16_type> markers );
1298 void renumber( mpl::bool_<false> ) {}
1308 void modifyEdgesOnBoundary( face_iterator& face, mpl::bool_<true> );
1313 void modifyEdgesOnBoundary( face_iterator& face, mpl::bool_<false> );
1318 bool modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<false> );
1323 bool modifyElementOnBoundaryFromEdge( element_iterator& e, mpl::bool_<true> );
1328 void updateOnBoundary();
1342 value_type M_measbdy;
1348 std::vector<uint16_type> M_neighboring_processors;
1355 boost::unordered_map<std::pair<int,int>,face_processor_type> M_e2f;
1360 boost::multi_array<element_edge_type,2> M_e2e;
1367 std::map<std::string, std::vector<size_type> > M_markername;
1372 std::map<uint64_type, boost::tuple<bool, std::vector<int>, std::vector<double> > > M_enc_pts;
1377 std::map<uint64_type, std::vector<int> > M_enc_faces;
1382 std::map<uint64_type, std::vector<int> > M_enc_elts;
1387 boost::shared_ptr<Localization> M_tool_localization;
1391 template<
typename Shape,
typename T,
int Tag>
1392 const uint16_type Mesh<Shape, T, Tag>::nDim;
1393 template<
typename Shape,
typename T,
int Tag>
1394 const uint16_type Mesh<Shape, T, Tag>::nOrder;
1396 template<
typename Shape,
typename T,
int Tag>
1397 template<
typename RangeT>
1398 typename Mesh<Shape, T, Tag>::template trace_mesh<Tag>::ptrtype
1401 DVLOG(2) <<
"[trace] extracting " << range.template get<0>() <<
" nb elements :"
1402 << std::distance(range.template get<1>(),range.template get<2>()) <<
"\n";
1403 return Feel::createSubmesh<const mesh_type,RangeT,Tag>( this->shared_from_this(), range );
1407 template<
typename Shape,
typename T,
int Tag>
1408 template<
typename RangeT>
1409 typename Mesh<Shape, T, Tag>::template trace_trace_mesh<Tag>::ptrtype
1410 Mesh<Shape, T, Tag>::wireBasket( RangeT
const& range )
const
1412 DVLOG(2) <<
"[trace] extracting " << range.template get<0>() <<
" nb elements :"
1413 << std::distance(range.template get<1>(),range.template get<2>()) <<
"\n";
1414 return Feel::createSubmesh<const mesh_type,RangeT,Tag>( this->shared_from_this(), range );
1419 template<
int TheTag>
struct Tag :
public mpl::int_<TheTag> {};
1421 template<
typename Shape,
typename T,
int Tag>
1422 template<
typename RangeT,
int TheTag>
1423 typename Mesh<Shape, T, Tag>::template trace_mesh<TheTag>::ptrtype
1426 DVLOG(2) <<
"[trace] extracting " << range.template get<0>() <<
" nb elements :"
1427 << std::distance(range.template get<1>(),range.template get<2>()) <<
"\n";
1428 return Feel::createSubmesh<const mesh_type,RangeT,TheTag>( this->shared_from_this(), range );
1432 template<
typename Shape,
typename T,
int Tag>
1433 template<
typename RangeT,
int TheTag>
1434 typename Mesh<Shape, T, Tag>::template trace_trace_mesh<TheTag>::ptrtype
1435 Mesh<Shape, T, Tag>::wireBasket( RangeT
const& range, mpl::int_<TheTag> )
const
1437 DVLOG(2) <<
"[trace] extracting " << range.template get<0>() <<
" nb elements :"
1438 << std::distance(range.template get<1>(),range.template get<2>()) <<
"\n";
1439 return Feel::createSubmesh<const mesh_type,RangeT,TheTag>( this->shared_from_this(), range );
1443 template<
typename Shape,
typename T,
int Tag>
1444 template<
typename Iterator>
1446 Mesh<Shape, T, Tag>::createSubmesh( self_type& new_mesh,
1447 Iterator
const& begin_elt,
1448 Iterator
const& end_elt,
1452 new_mesh = Feel::createSubmesh( this->shared_from_this(), boost::make_tuple(mpl::int_<MESH_ELEMENTS>(),begin_elt, end_elt), extraction_policies );
1454 Context policies( extraction_policies );
1456 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] start\n";
1460 new_mesh.M_markername = this->markerNames();
1461 BOOST_FOREACH(
auto marker, new_mesh.M_markername )
1463 LOG(INFO) <<
"marker name " << marker.first
1464 <<
" id: " << marker.second[0]
1465 <<
" geoe: " << marker.second[1] <<
"\n";
1470 std::vector<size_type> new_node_numbers ( this->numPoints() );
1471 std::vector<size_type> new_vertex ( this->numPoints() );
1473 std::fill ( new_node_numbers.begin(),
1474 new_node_numbers.end(),
1477 std::fill ( new_vertex.begin(),
1484 unsigned int n_new_nodes = 0;
1485 unsigned int n_new_elem = 0;
1488 for ( Iterator it = begin_elt; it != end_elt; ++it )
1491 element_type
const& old_elem = *it;
1494 element_type new_elem = old_elem;
1497 for (
unsigned int n=0; n < old_elem.nPoints(); n++ )
1499 FEELPP_ASSERT ( old_elem.point( n ).id() < new_node_numbers.size() ).error(
"invalid point id()" );
1503 new_node_numbers[old_elem.point( n ).id()] = n_new_nodes;
1505 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] insert point " << old_elem.point( n ) <<
"\n";
1507 point_type pt( old_elem.point( n ) );
1508 pt.setId( n_new_nodes );
1511 new_mesh.addPoint ( pt );
1513 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << new_mesh.numPoints() <<
"\n";
1518 if ( n < element_type::numVertices )
1520 FEELPP_ASSERT( new_vertex[old_elem.point( n ).id()] == 0 ).error(
"already seen this point?" );
1521 new_vertex[old_elem.point( n ).id()]=1;
1526 FEELPP_ASSERT ( new_node_numbers[old_elem.point( n ).id()] < new_mesh.numPoints() ).error(
"invalid connectivity" );
1528 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] adding point old(" << old_elem.point( n ).id()
1529 <<
") as point new(" << new_node_numbers[old_elem.point( n ).id()]
1530 <<
") in element " << new_elem.id() <<
"\n";
1532 new_elem.setPoint( n, new_mesh.point( new_node_numbers[old_elem.point( n ).id()] ) );
1537 new_elem.setId ( n_new_elem );
1544 new_mesh.addElement( new_elem );
1548 for (
unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
1550 if ( !old_elem.facePtr( s ) )
continue;
1553 std::cout <<
"local face id: " << s
1554 <<
" global face id: " << old_elem.face( s ).id() <<
"\n";
1559 size_type global_face_id = old_elem.face( s ).id();
1561 if ( this->hasFace( global_face_id ) )
1566 face_type
const& old_face = old_elem.face( s );
1567 face_type new_face = old_face;
1572 new_face.disconnect();
1576 for ( uint16_type p = 0; p < new_face.nPoints(); ++p )
1579 std::cout <<
"add new point " << new_face.point( p ).id() <<
" to face \n";
1580 std::cout <<
"add old point " << old_face.point( p ).id() <<
" to face \n";
1581 std::cout <<
"new point id " << new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).
id()] <<
"\n";
1583 new_face.setPoint( p, new_mesh.point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).
id()] ) );
1587 new_face.setId( n_new_faces++ );
1589 std::cout <<
"face id" << new_face.id()
1590 <<
" marker1 : " << new_face.marker()
1591 <<
" old marker1 : " << old_face.marker()
1595 new_mesh.addFace( new_face );
1600 new_mesh.setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0 ) );
1602 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] update face/edge info if necessary\n";
1604 new_mesh.components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
1605 new_mesh.updateForUse();
1607 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] stop\n";
1612 template<
typename Shape,
typename T,
int Tag>
1613 typename Mesh<Shape, T, Tag>::P1_mesh_ptrtype
1617 P1_mesh_ptrtype new_mesh(
new P1_mesh_type( this->worldComm() ) );
1621 boost::unordered_map<size_type,size_type> new_node_numbers;
1622 boost::unordered_map<size_type,int> new_vertex;
1623 std::vector<size_type> new_element_numbers ( this->numElements() );
1626 unsigned int n_new_nodes = 0;
1627 unsigned int n_new_elem = 0;
1631 BOOST_FOREACH(
auto itMark, this->markerNames() )
1633 new_mesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
1638 const int proc_id = this->worldComm().localRank();
1639 const int nProc = this->worldComm().localSize();
1640 std::vector< std::list<boost::tuple<size_type,size_type> > > memory_ghostid( nProc );
1641 std::vector< std::vector<size_type> > memory_id( nProc );
1642 std::vector< std::vector<size_type> > vecToSend( nProc );
1643 std::vector< std::vector<size_type> > vecToRecv( nProc );
1648 auto it = this->beginElement();
1649 auto en = this->endElement();
1651 for ( ; it != en; ++it )
1653 element_type
const& old_elem = *it;
1656 typename P1_mesh_type::element_type new_elem;
1658 new_elem.setId ( n_new_elem );
1659 new_element_numbers[old_elem.id()]= n_new_elem;
1661 new_elem.setMarker( old_elem.marker().value() );
1662 new_elem.setMarker2( old_elem.marker2().value() );
1663 new_elem.setMarker3( old_elem.marker3().value() );
1665 new_elem.setProcessIdInPartition( old_elem.pidInPartition() );
1666 new_elem.setNumberOfPartitions(old_elem.numberOfPartitions());
1667 new_elem.setProcessId(old_elem.processId());
1669 new_elem.setNeighborPartitionIds(old_elem.neighborPartitionIds());
1672 for ( uint16_type n=0; n < element_type::numVertices; n++ )
1674 auto const& old_point = old_elem.point( n );
1677 if ( new_node_numbers.find( old_point.id() ) == new_node_numbers.end() )
1679 new_node_numbers[old_point.id()] = n_new_nodes;
1680 DVLOG(2) <<
"[Mesh<Shape,T>::createP1mesh] insert point " << old_point <<
"\n";
1681 typename P1_mesh_type::point_type pt( old_point );
1682 pt.setId( n_new_nodes );
1683 pt.setProcessId(old_point.processId());
1685 new_mesh->addPoint( pt );
1686 DVLOG(2) <<
"[Mesh<Shape,T>::createSubmesh] number of points " << new_mesh->numPoints() <<
"\n";
1689 FEELPP_ASSERT( !new_vertex[old_point.id()] ).error(
"already seen this point?" );
1690 new_vertex[old_point.id()]=1;
1694 DVLOG(2) <<
"[Mesh<Shape,T>::createP1mesh] adding point old(" << old_point.id()
1695 <<
") as point new(" << new_node_numbers[old_point.id()]
1696 <<
") in element " << new_elem.id() <<
"\n";
1698 new_elem.setPoint( n, new_mesh->point( new_node_numbers[old_point.id()] ) );
1702 new_mesh->addElement( new_elem );
1707 for (
unsigned int s=0; s<old_elem.numTopologicalFaces; s++ )
1709 if ( !old_elem.facePtr( s ) )
continue;
1712 const size_type global_face_id = old_elem.face( s ).id();
1713 if ( this->hasFace( global_face_id ) )
1716 face_type
const& old_face = old_elem.face( s );
1717 typename P1_mesh_type::face_type new_face;
1720 new_face.disconnect();
1722 new_face.setOnBoundary( old_face.isOnBoundary() );
1724 new_face.setId( n_new_faces );
1726 new_face.setMarker( old_face.marker().value() );
1727 new_face.setMarker2( old_face.marker2().value() );
1728 new_face.setMarker2( old_face.marker3().value() );
1730 new_face.setProcessIdInPartition( old_face.pidInPartition() );
1731 new_face.setNumberOfPartitions(old_face.numberOfPartitions());
1732 new_face.setProcessId(old_face.processId());
1734 new_face.setNeighborPartitionIds(old_face.neighborPartitionIds());
1736 for ( uint16_type p = 0; p < face_type::numVertices; ++p )
1738 new_face.setPoint( p, new_mesh->point( new_node_numbers[old_elem.point( old_elem.fToP( s,p ) ).
id()] ) );
1741 new_mesh->addFace( new_face );
1748 if ( it->isGhostCell() )
1750 DVLOG(2) <<
"element " << it->id() <<
" is a ghost cell\n";
1751 for (
auto it_pid=it->idInOthersPartitions().begin(),en_pid=it->idInOthersPartitions().end() ; it_pid!=en_pid ; ++it_pid)
1753 DVLOG(2) <<
" " << it_pid->first <<
"-" << it_pid->second <<
"-"<<it->pidInPartition()<<
"-"<<new_mesh->worldComm().localRank();
1754 const int procToSend=it_pid->first;
1755 if (procToSend!=it->pidInPartition())
1757 memory_ghostid[procToSend].push_back(boost::make_tuple(new_elem.id(),it_pid->second));
1767 for (
int proc=0; proc<nProc; ++proc )
1769 vecToSend[proc].resize(memory_ghostid[proc].size());
1770 memory_id[proc].resize(memory_ghostid[proc].size());
1771 auto it_mem = memory_ghostid[proc].begin();
1772 auto en_mem = memory_ghostid[proc].end();
1773 for (
int k = 0 ; it_mem!=en_mem ; ++k,++it_mem )
1775 vecToSend[proc][k] = it_mem->get<1>();
1776 memory_id[proc][k] = it_mem->get<0>();
1781 this->worldComm().localComm().barrier();
1784 std::vector<size_type> nDataSize_vec(nProc);
1785 for (
int proc=0; proc<nProc; ++proc )
1787 const auto nDataSize = vecToSend[proc].size();
1788 mpi::all_gather( this->worldComm().localComm(),
1792 for (
int proc2=0; proc2<nProc; ++proc2 )
1794 if ( nDataSize_vec[proc2] > 0 && proc!=proc2 )
1798 this->worldComm().localComm().send( proc, 0, vecToSend[proc] );
1800 else if (proc_id==proc)
1802 this->worldComm().localComm().recv( proc2, 0, vecToRecv[proc2] );
1809 this->worldComm().localComm().barrier();
1812 for (
int proc=0 ; proc<nProc; ++proc )
1814 vecToSend[proc].resize(vecToRecv[proc].size());
1815 for (
int k=0;k<vecToRecv[proc].size();++k)
1817 const size_type idRequest = vecToRecv[proc][k];
1818 vecToSend[proc][k]= new_element_numbers[idRequest];
1823 this->worldComm().localComm().barrier();
1826 for (
int proc=0; proc<nProc; ++proc )
1828 const auto nDataSize = vecToSend[proc].size();
1829 mpi::all_gather( this->worldComm().localComm(),
1833 for (
int proc2=0; proc2<nProc; ++proc2 )
1835 if ( nDataSize_vec[proc2] > 0 && proc!=proc2 )
1839 this->worldComm().localComm().send( proc, 0, vecToSend[proc] );
1841 else if (proc_id==proc)
1843 this->worldComm().localComm().recv( proc2, 0, vecToRecv[proc2] );
1851 this->worldComm().localComm().barrier();
1854 for (
int proc=0 ; proc<nProc; ++proc )
1857 for (
int k=0;k<vecToRecv[proc].size();++k)
1859 auto elttt = new_mesh->elementIterator( memory_id[proc][k], proc );
1860 new_mesh->elements().modify( elttt, detail::updateIdInOthersPartitions( proc, vecToRecv[proc][k] ) );
1869 new_mesh->setNumVertices( std::accumulate( new_vertex.begin(), new_vertex.end(), 0,
1871 [](
size_type lhs, std::pair<size_type,int>
const& rhs )
1873 return lhs+rhs.second;
1877 DVLOG(2) <<
"[Mesh<Shape,T>::createP1mesh] update face/edge info if necessary\n";
1879 new_mesh->components().set ( MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
1881 new_mesh->updateForUse();
1887 template<
typename T>
1890 template<
typename MeshType,
typename IteratorType>
1891 MeshPoints( MeshType*, IteratorType it, IteratorType en,
const bool outer =
false,
const bool renumber =
false );
1893 std::vector<int> ids;
1894 std::map<int,int> new2old;
1895 std::map<int,int> old2new;
1896 std::map<int,int> nodemap;
1897 std::vector<T> coords;
1899 template<
typename T>
1900 template<
typename MeshType,
typename IteratorType>
1901 MeshPoints<T>::MeshPoints( MeshType* mesh, IteratorType it, IteratorType en,
const bool outer,
const bool renumber )
1903 std::set<int> nodeset;
1905 for( ; it != en; ++it )
1907 for (
size_type j = 0; j < it->numLocalVertices; j++ )
1909 int pid = it->point( j ).id();
1910 auto ins = nodeset.insert( pid );
1914 ids.push_back( p+1 );
1916 ids.push_back( pid );
1917 old2new[pid]=ids[p];
1918 new2old[ids[p]]=pid;
1924 CHECK( p == ids.size() ) <<
"Invalid number of points " << ids.size() <<
"!=" << p;
1925 int nv = ids.size();
1927 coords.resize( 3*nv, 0 );
1929 auto pit = ids.begin();
1930 auto pen = ids.end();
1932 for(
int i = 0; pit != pen; ++pit, ++i )
1934 CHECK( *pit > 0 ) <<
"invalid id " << *pit;
1939 auto const& p = mesh->point( new2old[*pit] );
1941 coords[i] = ( T ) p.node()[0];
1943 coords[3*i] = ( T ) p.node()[0];
1945 if ( MeshType::nRealDim >= 2 )
1948 coords[nv+i] = ( T ) p.node()[1];
1950 coords[3*i+1] = ( T ) p.node()[1];
1952 if ( MeshType::nRealDim >= 3 )
1955 coords[2*nv+i] = T( p.node()[2] );
1957 coords[3*i+2] = T( p.node()[2] );