AFEPack
|
00001 00011 #ifndef __MPI_HGeometry_h__ 00012 #define __MPI_HGeometry_h__ 00013 00014 #include <set> 00015 #include <AFEPack/HGeometry.h> 00016 #include "MPI.h" 00017 00018 namespace MPI { 00019 00020 template <int DIM, int DOW, class MATCHER> class HGeometryForest; 00021 template <class FOREST> class BirdView; 00022 template <class FOREST> class BirdViewSet; 00023 template <class FOREST> class HGeometryMatcher; 00024 template <class FOREST> class HLoadBalance; 00025 00026 template <int DOW> 00027 struct PointDistance { 00028 int value(const afepack::Point<DOW>& p0, 00029 const afepack::Point<DOW>& p1, 00030 double tol) const { 00031 return (distance(p0, p1) < tol)?0:-1; 00032 } 00033 }; 00034 00042 template <int DIM, int DOW=DIM, class MATCHER=PointDistance<DOW> > 00043 class HGeometryForest : public HGeometryTree<DIM,DOW> { 00044 private: 00045 typedef HGeometryForest<DIM,DOW,MATCHER> this_t; 00046 public: 00047 enum {dim = DIM, dow = DOW}; 00048 00049 typedef HGeometryTree<DIM,DOW> tree_t; 00050 typedef MATCHER matcher_t; 00051 00052 template <int D> 00053 const Shared_ptr_list<HGeometry<D,DOW> >& 00054 get_shared_list() const { 00055 const void * p_list; 00056 switch (D) { 00057 case 0: p_list = &_shared_list_0d; break; 00058 case 1: p_list = &_shared_list_1d; break; 00059 case 2: p_list = &_shared_list_2d; break; 00060 case 3: p_list = &_shared_list_3d; break; 00061 } 00062 return *((const Shared_ptr_list<HGeometry<D,DOW> > *)p_list); 00063 } 00064 00070 template <int GDIM> 00071 bool is_on_primary_rank(const HGeometry<GDIM,DOW>& geo) const { 00072 typedef HGeometry<GDIM,DOW> geo_t; 00073 typedef Shared_object<geo_t> obj_t; 00074 obj_t * p_obj = this->get_shared_info(geo); 00075 if (p_obj == NULL) { 00076 return true; 00077 } else { 00078 return p_obj->is_on_primary_rank(_rank); 00079 } 00080 } 00081 00085 template <int GDIM> 00086 int primary_rank(const HGeometry<GDIM,DOW>& geo) const { 00087 typedef HGeometry<GDIM,DOW> geo_t; 00088 typedef Shared_object<geo_t> obj_t; 00089 obj_t * p_obj = this->get_shared_info(geo); 00090 if (p_obj == NULL) { 00091 return _rank; 00092 } else { 00093 return p_obj->primary_rank(_rank); 00094 } 00095 } 00096 00100 template <int GDIM> 00101 bool is_primary_geometry(const HGeometry<GDIM,DOW>& geo) const { 00102 typedef HGeometry<GDIM,DOW> geo_t; 00103 typedef Shared_object<geo_t> obj_t; 00104 obj_t * p_obj = this->get_shared_info(geo); 00105 if (p_obj == NULL) { 00106 return true; 00107 } else { 00108 return p_obj->is_primary_object(_rank); 00109 } 00110 } 00111 00116 bool is_geometry_shared(const RegularMesh<DIM,DOW>& mesh, 00117 int dim, int idx) const { 00118 switch(dim) { 00119 case 0: return (get_shared_info(*mesh.template h_geometry<0>(idx)) != NULL); 00120 case 1: return (get_shared_info(*mesh.template h_geometry<1>(idx)) != NULL); 00121 case 2: return (get_shared_info(*mesh.template h_geometry<2>(idx)) != NULL); 00122 case 3: return (get_shared_info(*mesh.template h_geometry<3>(idx)) != NULL); 00123 } 00124 return false; 00125 } 00126 00127 const matcher_t& matcher() const { return _matcher; } 00128 matcher_t& matcher() { return _matcher; } 00129 00130 private: 00131 MPI_Comm _comm; 00132 int _rank; 00133 int _n_rank; 00134 matcher_t _matcher; 00135 00142 property_id_t<> _pid_shared_info_sent; 00143 template <class GEO> bool 00144 is_shared_info_sent(GEO& geo) const { 00145 return (geo.get_property(_pid_shared_info_sent) != NULL); 00146 } 00147 template <class GEO> void 00148 set_shared_info_sent(GEO& geo) const { 00149 if (! is_shared_info_sent(geo)) { 00150 geo.new_property(_pid_shared_info_sent); 00151 } 00152 } 00153 template <class GEO> void 00154 clear_shared_info_sent(GEO& geo) const { 00155 geo.free_property(_pid_shared_info_sent); 00156 } 00157 00165 property_id_t<> _pid_dummy; 00166 template <class GEO> bool 00167 is_dummy(GEO& geo) const { 00168 return (geo.get_property(_pid_dummy) != NULL); 00169 } 00170 template <class GEO> void 00171 dummy(GEO& geo) const { 00172 geo.new_property(_pid_dummy); 00173 } 00174 template <class GEO> void 00175 undummy(GEO& geo) const { 00176 geo.free_property(_pid_dummy); 00177 } 00178 00187 #define GDIM 0 00188 #define GEO HGeometry<GDIM,DOW> 00189 #define OBJ Shared_object<GEO> 00190 private: 00191 property_id_t<OBJ> _pid_so_0d; 00192 Shared_ptr_list<GEO> _shared_list_0d; 00193 00194 public: 00195 OBJ * get_shared_info(const GEO& geo) const { 00196 return geo.get_property(_pid_so_0d); 00197 } 00198 OBJ * new_shared_info(GEO& geo) { 00199 OBJ * p_obj = geo.new_property(_pid_so_0d); 00200 p_obj->local_pointer() = &geo; 00201 _shared_list_0d.push_back(p_obj); 00202 return p_obj; 00203 } 00204 void erase_shared_info(GEO& geo) { 00205 OBJ * p_obj = geo.get_property(_pid_so_0d); 00206 if (p_obj == NULL) return; 00207 _shared_list_0d.erase(std::find(_shared_list_0d.begin_ptr(), 00208 _shared_list_0d.end_ptr(), 00209 p_obj)); 00210 } 00211 #undef OBJ 00212 #undef GEO 00213 #undef GDIM 00214 00215 #define GDIM 1 00216 #define GEO HGeometry<GDIM,DOW> 00217 #define OBJ Shared_object<GEO> 00218 private: 00219 property_id_t<OBJ> _pid_so_1d; 00220 Shared_ptr_list<GEO> _shared_list_1d; 00221 00222 public: 00223 OBJ * get_shared_info(const GEO& geo) const { 00224 return geo.get_property(_pid_so_1d); 00225 } 00226 OBJ * new_shared_info(GEO& geo) { 00227 OBJ * p_obj = geo.new_property(_pid_so_1d); 00228 p_obj->local_pointer() = &geo; 00229 _shared_list_1d.push_back(p_obj); 00230 return p_obj; 00231 } 00232 void erase_shared_info(GEO& geo) { 00233 OBJ * p_obj = geo.get_property(_pid_so_1d); 00234 if (p_obj == NULL) return; 00235 _shared_list_1d.erase(std::find(_shared_list_1d.begin_ptr(), 00236 _shared_list_1d.end_ptr(), 00237 p_obj)); 00238 } 00239 #undef OBJ 00240 #undef GEO 00241 #undef GDIM 00242 00243 #define GDIM 2 00244 #define GEO HGeometry<GDIM,DOW> 00245 #define OBJ Shared_object<GEO> 00246 private: 00247 property_id_t<OBJ> _pid_so_2d; 00248 Shared_ptr_list<GEO> _shared_list_2d; 00249 00250 public: 00251 OBJ * get_shared_info(const GEO& geo) const { 00252 return geo.get_property(_pid_so_2d); 00253 } 00254 OBJ * new_shared_info(GEO& geo) { 00255 OBJ * p_obj = geo.new_property(_pid_so_2d); 00256 p_obj->local_pointer() = &geo; 00257 _shared_list_2d.push_back(p_obj); 00258 return p_obj; 00259 } 00260 void erase_shared_info(GEO& geo) { 00261 OBJ * p_obj = geo.get_property(_pid_so_2d); 00262 if (p_obj == NULL) return; 00263 _shared_list_2d.erase(std::find(_shared_list_2d.begin_ptr(), 00264 _shared_list_2d.end_ptr(), 00265 p_obj)); 00266 } 00267 #undef OBJ 00268 #undef GEO 00269 #undef GDIM 00270 00271 #define GDIM 3 00272 #define GEO HGeometry<GDIM,DOW> 00273 #define OBJ Shared_object<GEO> 00274 private: 00275 property_id_t<OBJ> _pid_so_3d; 00276 Shared_ptr_list<GEO> _shared_list_3d; 00277 00278 public: 00279 OBJ * get_shared_info(const GEO& geo) const { 00280 return geo.get_property(_pid_so_3d); 00281 } 00282 OBJ * new_shared_info(GEO& geo) { 00283 OBJ * p_obj = geo.new_property(_pid_so_3d); 00284 p_obj->local_pointer() = &geo; 00285 _shared_list_3d.push_back(p_obj); 00286 return p_obj; 00287 } 00288 void erase_shared_info(GEO& geo) { 00289 OBJ * p_obj = geo.get_property(_pid_so_3d); 00290 if (p_obj == NULL) return; 00291 _shared_list_3d.erase(std::find(_shared_list_3d.begin_ptr(), 00292 _shared_list_3d.end_ptr(), 00293 p_obj)); 00294 } 00295 #undef OBJ 00296 #undef GEO 00297 #undef GDIM 00298 00299 bool lock() { return HGeometryTree<DIM,DOW>::lock(); } 00300 void unlock() { HGeometryTree<DIM,DOW>::unlock(); } 00301 00302 public: 00303 HGeometryForest() { new_property(); } 00304 HGeometryForest(MPI_Comm comm) { 00305 set_communicator(comm); 00306 new_property(); 00307 } 00308 virtual ~HGeometryForest() { this->clear(); } 00309 00310 void clear() { 00311 _shared_list_0d.clear(); 00312 _shared_list_1d.clear(); 00313 _shared_list_2d.clear(); 00314 _shared_list_3d.clear(); 00315 00316 free_property(); 00317 00324 typename tree_t::container_t nre; 00325 property_id_t<bool> pid_is_found; 00326 new_property_id(pid_is_found); 00327 typename tree_t::RootIterator 00328 the_ele = this->beginRootElement(), 00329 end_ele = this->endRootElement(); 00330 for (;the_ele != end_ele;++ the_ele) { 00331 this->collect_parent_element(pid_is_found, *the_ele, nre); 00332 } 00333 free_property_id(pid_is_found); 00334 00336 this->rootElement().swap(nre); 00337 nre.clear(); 00338 tree_t::clear(); 00339 00340 new_property(); 00341 } 00342 00343 private: 00344 void collect_parent_element(const property_id_t<bool>& pid, 00345 HGeometry<dim,dow>& geo, 00346 typename tree_t::container_t& nre) const { 00347 if (geo.get_property(pid) != NULL) return; 00348 if (geo.parent != NULL) { 00349 this->collect_parent_element(pid, *geo.parent, nre); 00350 } else { 00351 nre.push_back(&geo); 00352 geo.new_property(pid); 00353 } 00354 } 00355 00356 public: 00357 void new_property() { 00358 new_property_id(_pid_shared_info_sent); 00359 new_property_id(_pid_dummy); 00360 if (DIM >= 0) new_property_id(_pid_so_0d); 00361 if (DIM >= 1) new_property_id(_pid_so_1d); 00362 if (DIM >= 2) new_property_id(_pid_so_2d); 00363 if (DIM >= 3) new_property_id(_pid_so_3d); 00364 } 00365 00366 void free_property() { 00367 free_property_id(_pid_shared_info_sent); 00368 free_property_id(_pid_dummy); 00369 if (DIM >= 0) free_property_id(_pid_so_0d); 00370 if (DIM >= 1) free_property_id(_pid_so_1d); 00371 if (DIM >= 2) free_property_id(_pid_so_2d); 00372 if (DIM >= 3) free_property_id(_pid_so_3d); 00373 } 00374 00375 void set_communicator(MPI_Comm comm) { 00376 _comm = comm; 00377 MPI_Comm_rank(_comm, &_rank); 00378 MPI_Comm_size(_comm, &_n_rank); 00379 } 00380 MPI_Comm communicator() const { return _comm; } 00381 int rank() const { return _rank; } 00382 int n_rank() const { return _n_rank; } 00383 void readMesh(const std::string&); 00384 00385 void eraseRootElement(u_int level = 1); 00386 void renumerateRootElement(void (*f)(double, double, double, 00387 double&,double&,double&) = NULL); 00388 00392 void set_is_used_up(); 00393 template <class HGEO> void 00394 set_is_used_up(HGEO& geo) const; 00395 00401 void set_shared_info_sent() const { 00402 if (dim >= 3) set_dim_shared_info_sent<3>(); 00403 if (dim >= 2) set_dim_shared_info_sent<2>(); 00404 if (dim >= 1) set_dim_shared_info_sent<1>(); 00405 if (dim >= 0) set_dim_shared_info_sent<0>(); 00406 } 00407 template <int D> 00408 void set_dim_shared_info_sent() const { 00409 typedef HGeometry<D,dow> geo_t; 00410 typedef Shared_object<geo_t> obj_t; 00411 typedef Shared_ptr_list<geo_t> list_t; 00412 00413 const list_t& geo_ptr_list = this->template get_shared_list<D>(); 00414 typename list_t::const_iterator 00415 the_geo_ptr = geo_ptr_list.begin(), 00416 end_geo_ptr = geo_ptr_list.end(); 00417 for (;the_geo_ptr != end_geo_ptr;++ the_geo_ptr) { 00418 geo_t * p_geo = the_geo_ptr->local_pointer(); 00419 this->set_shared_info_sent(*p_geo); 00420 } 00421 } 00422 00427 public: 00428 template <int GDIM> void 00429 pack_verify_shared_object(HGeometry<GDIM,dow> * geo, 00430 int remote_rank, 00431 AFEPack::ostream<>& os) { 00432 typedef HGeometry<GDIM,dow> geo_t; 00433 os << geo; 00434 00435 Shared_object<geo_t> * p_info = this->get_shared_info(*geo); 00436 std::size_t n = p_info->size(); 00437 os << n; 00438 typename std::multimap<int,Remote_pointer<geo_t> >::iterator 00439 the_entry = p_info->begin(), end_entry = p_info->end(); 00440 for (;the_entry != end_entry;++ the_entry) { 00441 os << the_entry->first 00442 << the_entry->second.type 00443 << the_entry->second.ptr; 00444 } 00445 } 00446 00447 template <int GDIM> void 00448 unpack_verify_shared_object(HGeometry<GDIM,dow> * geo, 00449 int remote_rank, 00450 AFEPack::istream<>& is) { 00451 typedef HGeometry<GDIM,dow> geo_t; 00452 geo_t * p_remote_geo; 00453 is >> p_remote_geo; 00454 00455 Shared_object<geo_t> * p_info = this->get_shared_info(*geo); 00456 assert (p_info != NULL); 00457 00458 std::size_t n; 00459 is >> n; 00460 Remote_pointer<geo_t> p_remote_ptr; 00461 for (std::size_t i = 0;i < n;++ i) { 00462 int orank; 00463 is >> orank >> p_remote_ptr.type >> p_remote_ptr.ptr; 00464 if ((orank == this->rank()) && 00465 (p_remote_ptr.ptr == geo)) { 00466 orank = remote_rank; 00467 p_remote_ptr.ptr = p_remote_geo; 00468 } 00469 if (! p_info->is_duplicate_entry(orank, p_remote_ptr)) { 00470 MPI_Abort(this->communicator(), 0); 00471 } 00472 } 00473 } 00474 00475 void verify_shared_object() { 00476 if (this->n_rank() <= 1) return; 00477 Shared_type_filter::all type_filter; 00478 #define SYNC_DATA(D) \ 00479 if (dim >= D) { \ 00480 sync_data(this->communicator(), \ 00481 this->template get_shared_list<D>(), *this, \ 00482 &this_t::template pack_verify_shared_object<D>, \ 00483 &this_t::template unpack_verify_shared_object<D>, \ 00484 type_filter); \ 00485 } 00486 SYNC_DATA(0); 00487 SYNC_DATA(1); 00488 SYNC_DATA(2); 00489 SYNC_DATA(3); 00490 #undef SYNC_DATA 00491 } 00492 00493 private: 00494 void eraseRootElementOneLevel(); 00495 template <class GEO> 00496 void nullParent(GEO& geo) { 00497 geo.parent = NULL; 00498 for (u_int i = 0;i < GEO::n_boundary;++ i) { 00499 this->nullParent(*geo.boundary[i]); 00500 } 00501 } 00502 template <class GEO> 00503 void tryDeleteGeometry(GEO * p_geo, 00504 const property_id_t<bool>& pid) { 00505 for (u_int i = 0;i < GEO::n_boundary;++ i) { 00506 if (p_geo->boundary[i]->get_property(pid) != NULL) { 00507 p_geo->boundary[i] = NULL; 00508 } else { 00509 this->tryDeleteGeometry(p_geo->boundary[i], pid); 00510 } 00511 } 00512 bool * p_prp = p_geo->get_property(pid); 00513 assert (p_prp == NULL); 00514 p_geo->new_property(pid); 00515 this->erase_shared_info(*p_geo); 00516 } 00517 template <class GEO> 00518 void deleteGeometry(GEO * p_geo) { 00519 for (u_int i = 0;i < GEO::n_boundary;++ i) { 00520 if (p_geo->boundary[i] != NULL) { 00521 this->deleteGeometry(p_geo->boundary[i]); 00522 } 00523 } 00524 delete p_geo; 00525 } 00526 00527 friend class BirdView<this_t>; 00528 friend class HGeometryMatcher<this_t>; 00529 friend class HLoadBalance<this_t>; 00530 00531 public: 00532 typedef BirdView<this_t> birdview_t; 00533 typedef BirdViewSet<this_t> birdview_set_t; 00534 typedef HLoadBalance<this_t> load_balancer_t; 00535 00536 }; // HGeometryForest 00537 00542 template <class FOREST> 00543 class BirdView : public IrregularMesh<FOREST::dim,FOREST::dow> { 00544 public: 00545 enum {dim = FOREST::dim, dow = FOREST::dow}; 00546 typedef IrregularMesh<dim, dow> base_t; 00547 typedef FOREST forest_t; 00548 public: 00549 BirdView() : base_t() {} 00550 explicit BirdView(forest_t& forest) : base_t(forest) {} 00551 virtual ~BirdView() {} 00552 00553 forest_t& getForest() const { 00554 return dynamic_cast<forest_t&>(this->geometryTree()); 00555 } 00556 virtual void semiregularize(); 00557 void eraseRootElement(u_int level = 1); 00558 private: 00559 void eraseRootElementOneLevel(); 00560 }; 00561 00566 template <class FOREST> 00567 class BirdViewSet : public std::set<BirdView<FOREST>*> { 00568 private: 00569 typedef BirdView<FOREST> birdview_t; 00570 typedef std::set<birdview_t*> base_t; 00571 public: 00572 typedef _Deref_iterator<typename base_t::iterator, birdview_t> iterator; 00573 typedef _Deref_iterator<typename base_t::const_iterator, const birdview_t> const_iterator; 00574 public: 00575 void add(birdview_t& bv) { base_t::insert(&bv); } 00576 void erase(birdview_t& bv) { base_t::erase(&bv); } 00577 00578 iterator begin() { return base_t::begin(); } 00579 iterator end() { return base_t::end(); } 00580 00581 const_iterator begin() const { return base_t::begin(); } 00582 const_iterator end() const { return base_t::end(); } 00583 00584 typename base_t::iterator begin_ptr() { return base_t::begin(); } 00585 typename base_t::iterator end_ptr() { return base_t::end(); } 00586 00587 typename base_t::const_iterator begin_ptr() const { return base_t::begin(); } 00588 typename base_t::const_iterator end_ptr() const { return base_t::end(); } 00589 }; 00590 00591 template <class FOREST> 00592 BirdViewSet<FOREST> 00593 make_set(BirdView<FOREST>& mesh0, 00594 BirdView<FOREST>& mesh1) { 00595 BirdViewSet<FOREST> bvs; 00596 bvs.add(mesh0); 00597 bvs.add(mesh1); 00598 return bvs; 00599 } 00600 template <class FOREST> 00601 BirdViewSet<FOREST> 00602 make_set(BirdView<FOREST>& mesh0, 00603 BirdView<FOREST>& mesh1, 00604 BirdView<FOREST>& mesh2) { 00605 BirdViewSet<FOREST> bvs; 00606 bvs.add(mesh0); 00607 bvs.add(mesh1); 00608 bvs.add(mesh2); 00609 return bvs; 00610 } 00611 template <class FOREST> 00612 BirdViewSet<FOREST> 00613 make_set(BirdView<FOREST>& mesh0, 00614 BirdView<FOREST>& mesh1, 00615 BirdView<FOREST>& mesh2, 00616 BirdView<FOREST>& mesh3) { 00617 BirdViewSet<FOREST> bvs; 00618 bvs.add(mesh0); 00619 bvs.add(mesh1); 00620 bvs.add(mesh2); 00621 bvs.add(mesh3); 00622 return bvs; 00623 } 00624 00648 template <class FOREST> 00649 class HGeometryMatcher { 00650 public: 00651 enum {dim = FOREST::dim, dow = FOREST::dow}; 00652 typedef FOREST forest_t; 00653 private: 00654 typedef HGeometryMatcher<forest_t> this_t; 00655 typedef typename forest_t::matcher_t matcher_t; 00656 00657 HTools _tools; 00658 forest_t * _forest; 00659 bool _is_operated; 00660 00661 const matcher_t& matcher() const { return _forest->matcher(); } 00662 00663 template <class GEO> bool 00664 is_shared_info_sent(GEO& geo) const { 00665 return _forest->is_shared_info_sent(geo); 00666 } 00667 template <class GEO> void 00668 set_shared_info_sent(GEO& geo) const { 00669 _forest->set_shared_info_sent(geo); 00670 } 00671 template <class GEO> Shared_object<GEO> * 00672 new_shared_info(GEO& geo) const { 00673 return _forest->new_shared_info(geo); 00674 } 00675 template <class GEO> Shared_object<GEO> * 00676 get_shared_info(GEO& geo) const { 00677 return _forest->get_shared_info(geo); 00678 } 00679 00680 public: 00681 HGeometryMatcher(forest_t& forest) : _forest(&forest) {} 00682 00689 bool match_geometry(); 00690 template <int GDIM> bool 00691 is_pack_match_geometry(HGeometry<GDIM,dow> * geo) const; 00692 template <int GDIM> void 00693 pack_match_geometry(HGeometry<GDIM,dow> * geo, 00694 int remote_rank, 00695 AFEPack::ostream<>& os); 00696 template <int GDIM> void 00697 unpack_match_geometry(HGeometry<GDIM,dow> * geo, 00698 int remote_rank, 00699 AFEPack::istream<>& is); 00700 00701 public: 00705 void set_matcher(const matcher_t& _matcher) { matcher = _matcher; } 00706 00713 bool sync_is_used(); 00714 template <int GDIM> void 00715 pack_sync_is_used(HGeometry<GDIM,dow> * geo, 00716 int remote_rank, 00717 AFEPack::ostream<>& os); 00718 template <int GDIM> void 00719 unpack_sync_is_used(HGeometry<GDIM,dow> * geo, 00720 int remote_rank, 00721 AFEPack::istream<>& is); 00722 00723 private: 00727 template <class GEO> 00728 void geometry_reference_point(GEO& geo, 00729 afepack::Point<dow>& pnt) const { 00730 int n_vtx = geo.n_vertex; 00731 for (int n = 0;n < dow;++ n) { 00732 pnt[n] = 0.0; 00733 for (int i = 0;i < n_vtx;++ i) { 00734 pnt[n] += (*geo.vertex[i])[n]; 00735 } 00736 pnt[n] /= n_vtx; 00737 } 00738 } 00739 void geometry_reference_point(HGeometry<0,dow>& geo, 00740 afepack::Point<dow>& pnt) const { 00741 for (int n = 0;n < dow;++ n) { 00742 pnt[n] = geo[n]; 00743 } 00744 } 00745 }; 00746 00747 00748 #include "MPI_HGeometry.templates.h" 00749 00750 } // namespace MPI 00751 00752 #endif // __MPI_HGeometry_h__ 00753