AFEPack
MPI_ULoadBalance.h
浏览该文件的文档。
00001 
00053 #ifndef __MPI_ULoadBalance_h__
00054 #define __MPI_ULoadBalance_h__
00055 
00056 #ifdef __MPI_LoadBalance_h__
00057 #error "MPI_LoadBalance.h 和 MPI_ULoadBalance.h 是冲突的,不能一起用。"
00058 #endif
00059 
00060 #include <sys/types.h>
00061 #include <sys/stat.h>
00062 #include <unistd.h>
00063 #include <cstdio>
00064 
00065 #include <boost/program_options.hpp>
00066 #include <boost/archive/binary_iarchive.hpp>
00067 #include <boost/archive/binary_oarchive.hpp>
00068 #include <boost/archive/text_iarchive.hpp>
00069 #include <boost/archive/text_oarchive.hpp>
00070 
00071 #include <AFEPack/Serialization.h>
00072 #include "MPI.h"
00073 #include "MPI_UGeometry_archive.h"
00074 #include "MPI_Migration.h"
00075 
00076 namespace MPI {
00077 
00078   template <class FOREST>
00079     class HLoadBalance {
00080   public:
00081     typedef FOREST forest_t;
00082     enum {dim = FOREST::dim, dow = FOREST::dow};
00083     typedef typename forest_t::matcher_t matcher_t;
00084     
00085     struct rank_map_val {
00086       int old_rank; 
00087       int priority; 
00092     };
00093     typedef std::map<int,rank_map_val> rank_map_t;
00094 
00095     template <class GEO> bool 
00096       is_refined_on_new_rank(const GEO& geo, 
00097                              int new_rank) const {
00098       bool result = false;
00099       if (geo.isRefined() &&
00100           is_on_this_new_rank(*geo.child[0], new_rank)) {
00101         result = true;
00102       }
00103       return result;
00104     }
00105 
00106   private:
00107   // 在计算每个几何体上的负载是使用的性质
00108   property_id_t<u_int> _pid_loading;
00109 
00110   // 标记一个共享几何体是否已经存储了的性质
00111   property_id_t<> _pid_is_saved;
00112   public:
00113   template <class GEO> void
00114     set_is_saved(const GEO& geo) const {
00115     if (geo.get_property(_pid_is_saved) == NULL) {
00116       geo.new_property(_pid_is_saved);
00117     }
00118   }
00119 
00120   private:
00125   property_id_t<rank_map_t> _pid_rank_map;
00126 
00127   // 部分需要拆分和合并的几何体的全局标号
00128   property_id_t<unsigned long> _pid_global_idx;
00129 
00130   // 元素为(全局标号,对象指针)的对
00131   std::map<unsigned long, void *> _global_pointer_to_merge;
00132   std::map<unsigned long, std::list<void *> > _global_pointer_to_merge_later;
00133 
00134   // (进程,列表(全局标号,维数,对象指针)组)的映射
00135   typedef std::map<unsigned long, std::pair<int, void *> > tuple_map_t;
00136   std::map<int, tuple_map_t> _global_pointer_to_share;
00137 
00138   typedef BirdView<forest_t> birdview_t;
00139   typedef BirdViewSet<forest_t> birdview_set_t;
00140   typedef HLoadBalance<forest_t> this_t;
00141 
00142   forest_t * _forest; 
00143   std::vector<u_int> _cut_point; 
00144   std::vector<int> _new_rank; 
00145 
00146   typename forest_t::container_t _nre; 
00147 
00148   public:
00149   HLoadBalance(forest_t& forest) : _forest(&forest) {}
00150   ~HLoadBalance() {}
00151 
00152   forest_t& forest() const { return *_forest; }
00153   int rank() const { return _forest->rank(); }
00154 
00155   private:
00156   template <class GEO> rank_map_t *
00157   new_rank_map(const GEO& geo) const {
00158     return geo.new_property(_pid_rank_map);
00159   }
00160   template <class GEO> unsigned long *
00161   new_global_idx(const GEO& geo) const {
00162     return geo.new_property(_pid_global_idx);
00163   }
00164 
00165   public:
00166   void clear() {
00167     free_property_id(_pid_loading);
00168     free_property_id(_pid_rank_map);
00169     free_property_id(_pid_global_idx);
00170     _global_pointer_to_merge.clear();
00171     _global_pointer_to_merge_later.clear();
00172     _global_pointer_to_share.clear();
00173     _cut_point.clear();
00174     _new_rank.clear();
00175   }
00176 
00177   template <class GEO> rank_map_t *
00178   get_rank_map(const GEO& geo) const {
00179     return geo.get_property(_pid_rank_map);
00180   }
00181   template <class GEO> unsigned long *
00182   get_global_idx(const GEO& geo) const {
00183     return geo.get_property(_pid_global_idx);
00184   }
00185 
00187 
00190   template <class GEO> void
00191     print_geo_rank_map(const GEO& geo) const {
00192     std::cerr << "[" << _forest->rank() << ":" << geo.dimension;
00193     rank_map_t * p_map = get_rank_map(geo);
00194     if (p_map != NULL) print_rank_map(p_map);
00195     std::cerr << "]";
00196   }
00197 
00198   void print_rank_map(const rank_map_t * p_map) const {
00199     typename rank_map_t::const_iterator
00200       the_pair = p_map->begin(),
00201       end_pair = p_map->end();
00202     for (;the_pair != end_pair;++ the_pair) {
00203       std::cerr << "(" << the_pair->first << ":"
00204                 << the_pair->second.old_rank << ":"
00205                 << the_pair->second.priority << ")";
00206     }
00207   }
00209 
00213   template <class GEO> bool 
00214     is_on_this_new_rank(const GEO& geo,
00215                         int new_rank) const {
00216     rank_map_t * p_map = this->get_rank_map(geo);
00217     if (p_map == NULL) return false;
00218     typename rank_map_t::iterator
00219       the_pair = p_map->find(new_rank);
00220     return (the_pair != p_map->end());
00221   }
00222 
00231   template <class GEO> bool
00232     is_save_on_this_rank(const GEO& geo,
00233                          int new_rank) const {
00234     rank_map_t * p_map = this->get_rank_map(geo);
00235     assert (p_map != NULL);
00236     typename rank_map_t::iterator
00237       the_pair = p_map->find(new_rank);
00239     if (the_pair == p_map->end()) return false;
00240 
00242     return (the_pair->second.old_rank == _forest->rank()); 
00243   }
00244 
00256   template <class GEO>
00257   void merge_global_pointer(int type,
00258                             unsigned long global_idx,
00259                             GEO *& p_geo) {
00260     std::map<unsigned long, void *>::iterator
00261     the_pair = _global_pointer_to_merge.find(global_idx);
00262     if (type == 2) {
00263       if (the_pair == _global_pointer_to_merge.end()) {
00264         _global_pointer_to_merge.insert(std::pair<unsigned long,void*>(global_idx, p_geo));
00265 
00269         std::map<unsigned long, std::list<void *> >::iterator
00270           the_entry = _global_pointer_to_merge_later.find(global_idx);
00271         if (the_entry != _global_pointer_to_merge_later.end()) {
00272           std::list<void *>::iterator
00273             the_ptr = the_entry->second.begin(),
00274             end_ptr = the_entry->second.end();
00275           for (;the_ptr != end_ptr;++ the_ptr) {
00276             GEO ** geo_ptr_ptr = (GEO **)(*the_ptr);
00277             (*geo_ptr_ptr) = p_geo;
00278           }
00279           _global_pointer_to_merge_later.erase(the_entry);
00280         }
00281       } else {
00282         assert (the_pair->second == p_geo);
00283       }
00284     } else {
00285       assert (type == 4);
00286       if (the_pair != _global_pointer_to_merge.end()) {
00287         p_geo = (GEO *)(the_pair->second);
00288       } else {
00289         _global_pointer_to_merge_later[global_idx].push_back((void *)(&p_geo));
00290       }
00291     }
00292   }
00293 
00302   template <class GEO>
00303   void share_global_pointer(int rank, 
00304                             unsigned long global_idx,
00305                             GEO *& p_geo) {
00306     _global_pointer_to_share[rank][global_idx] = 
00307       std::pair<int,void *>(p_geo->dimension, (void *)(p_geo));
00308   }
00309 
00310   private:
00314   void share_global_pointer() {
00315     std::cerr << "Exchanging shared global pointers ..." << std::endl;
00316 
00317     typedef std::map<int, tuple_map_t> map_t;
00318 
00319     int n = 0;
00320     std::list<int> target;
00321     std::list<BinaryBuffer<> > out_buf, in_buf;
00322     typename map_t::iterator
00323     the_pair = _global_pointer_to_share.begin(),
00324     end_pair = _global_pointer_to_share.end();
00325     for (;the_pair != end_pair;++ the_pair, ++ n) {
00326       target.push_back(the_pair->first);
00327 
00328       out_buf.push_back(BinaryBuffer<>());
00329       AFEPack::ostream<> os(out_buf.back());
00330       
00331       in_buf.push_back(BinaryBuffer<>());
00332 
00333       u_int n_item = the_pair->second.size();
00334       os << n_item;
00335 #if 0
00336       std::cerr << "sending " << the_pair->second.size()
00337                 << " items from " << _forest->rank()
00338                 << " to " << the_pair->first << std::endl;
00339 #endif
00340       typename tuple_map_t::iterator
00341         the_tuple = the_pair->second.begin(),
00342         end_tuple = the_pair->second.end();
00343       for (;the_tuple != end_tuple;++ the_tuple) {
00344         const unsigned long& global_idx = the_tuple->first;
00345         int& dimension = the_tuple->second.first;
00346         void *& p_obj = the_tuple->second.second;
00347         os << global_idx << dimension << p_obj;
00348       }
00349     }
00350 
00351     MPI_Barrier(_forest->communicator());
00352     sendrecv_data(_forest->communicator(), n, out_buf.begin(), 
00353                   in_buf.begin(), target.begin());
00354 
00355     typename std::list<int>::iterator 
00356     the_rank = target.begin(),
00357     end_rank = target.end();
00358     typename std::list<BinaryBuffer<> >::iterator 
00359     the_buf = in_buf.begin();
00360     for (;the_rank != end_rank;++ the_rank, ++ the_buf) {
00361       AFEPack::istream<> is(*the_buf);
00362       u_int n_item;
00363       is >> n_item;
00364 #if 0
00365       std::cerr << "received " << n_item
00366                 << " items from " << *the_rank
00367                 << " to " << _forest->rank() << std::endl;
00368 #endif
00369       for (u_int i = 0;i < n_item;++ i) {
00370         unsigned long global_idx;
00371         int dimension;
00372         void * remote_obj;
00373         is >> global_idx >> dimension >> remote_obj;
00374         std::pair<int,void *>& the_pair = 
00375           _global_pointer_to_share[*the_rank][global_idx];
00376         assert ((dimension == the_pair.first) && 
00377                 (dimension >= 0) && 
00378                 (dimension <= 3));
00379 
00380         if (dimension == 0) {
00381 #define GDIM 0
00382 #define GEO HGeometry<GDIM,dow>
00383 #define OBJ Shared_object<GEO>
00384           GEO * p_geo = (GEO *)(the_pair.second);
00385           OBJ * p_info = _forest->get_shared_info(*p_geo);
00386           if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
00387           p_info->add_clone(*the_rank, (GEO *)remote_obj);
00388 #undef OBJ
00389 #undef GEO
00390 #undef GDIM
00391         } else if (dimension == 1) {
00392 #define GDIM 1
00393 #define GEO HGeometry<GDIM,dow>
00394 #define OBJ Shared_object<GEO>
00395           GEO * p_geo = (GEO *)(the_pair.second);
00396           OBJ * p_info = _forest->get_shared_info(*p_geo);
00397           if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
00398           p_info->add_clone(*the_rank, (GEO *)remote_obj);
00399 #undef OBJ
00400 #undef GEO
00401 #undef GDIM
00402         } else if (dimension == 2) {
00403 #define GDIM 2
00404 #define GEO HGeometry<GDIM,dow>
00405 #define OBJ Shared_object<GEO>
00406           GEO * p_geo = (GEO *)(the_pair.second);
00407           OBJ * p_info = _forest->get_shared_info(*p_geo);
00408           if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
00409           p_info->add_clone(*the_rank, (GEO *)remote_obj);
00410 #undef OBJ
00411 #undef GEO
00412 #undef GDIM
00413         } else if (dimension == 3) {
00414 #define GDIM 3
00415 #define GEO HGeometry<GDIM,dow>
00416 #define OBJ Shared_object<GEO>
00417           GEO * p_geo = (GEO *)(the_pair.second);
00418           OBJ * p_info = _forest->get_shared_info(*p_geo);
00419           if (p_info == NULL) p_info = _forest->new_shared_info(*p_geo);
00420           p_info->add_clone(*the_rank, (GEO *)remote_obj);
00421 #undef OBJ
00422 #undef GEO
00423 #undef GDIM
00424         }
00425       }
00426     }
00427 
00428 #if 1 
00429 
00432     the_pair = _global_pointer_to_share.begin();
00433     for (;the_pair != end_pair;++ the_pair) {
00434       typename tuple_map_t::iterator
00435         the_tuple = the_pair->second.begin(),
00436         end_tuple = the_pair->second.end();
00437       for (;the_tuple != end_tuple;++ the_tuple) {
00438         int& dimension = the_tuple->second.first;
00439         void *& p_obj = the_tuple->second.second;
00440         if (dimension == 0) {
00441 #define GDIM 0
00442 #define GEO HGeometry<GDIM,dow>
00443           GEO * p_geo = (GEO *)p_obj;
00444           assert (_forest->get_shared_info(*p_geo) != NULL);
00445 #undef GEO
00446 #undef GDIM
00447         } else if (dimension == 1) {
00448 #define GDIM 1
00449 #define GEO HGeometry<GDIM,dow>
00450           GEO * p_geo = (GEO *)p_obj;
00451           assert (_forest->get_shared_info(*p_geo) != NULL);
00452 #undef GEO
00453 #undef GDIM
00454         } else if (dimension == 2) {
00455 #define GDIM 2
00456 #define GEO HGeometry<GDIM,dow>
00457           GEO * p_geo = (GEO *)p_obj;
00458           assert (_forest->get_shared_info(*p_geo) != NULL);
00459 #undef GEO
00460 #undef GDIM
00461         } else if (dimension == 3) {
00462 #define GDIM 3
00463 #define GEO HGeometry<GDIM,dow>
00464           GEO * p_geo = (GEO *)p_obj;
00465           assert (_forest->get_shared_info(*p_geo) != NULL);
00466 #undef GEO
00467 #undef GDIM
00468         }
00469       }
00470     }
00471 #endif
00472   }
00473 
00475   public:
00481   template <class LOADER>
00482   void config(birdview_t& ir_mesh,
00483               LOADER& loader,
00484               u_int (LOADER::*value)(const GeometryBM&)) {
00485     if (&(ir_mesh.getForest()) != _forest) {
00486       std::cerr << "The mesh should be on the same forest of the HLoadBalance." 
00487                 << std::endl;
00488       abort();
00489     }
00490 
00491     RegularMesh<dim,dow>& mesh = ir_mesh.regularMesh();
00492     new_property_id(_pid_loading); 
00493     u_int n_ele = mesh.n_geometry(dim);
00494     for (u_int i = 0;i < n_ele;++ i) {
00495       HGeometry<dim,dow> * p_geo = mesh.template h_geometry<dim>(i);
00496       u_int * p_loading = p_geo->get_property(_pid_loading);
00497       if (p_loading == NULL) {
00498         p_loading = p_geo->new_property(_pid_loading);
00499       }
00500       (*p_loading) = (loader.*value)(mesh.geometry(dim,i));
00501     }
00502   }
00503 
00504   void config(birdview_t& ir_mesh,
00505               u_int (*value)(const GeometryBM&) = &_default_element_loading) {
00506     if (&(ir_mesh.getForest()) != _forest) {
00507       std::cerr << "The mesh should be on the same forest of the HLoadBalance." 
00508                 << std::endl;
00509       abort();
00510     }
00511 
00512     RegularMesh<dim,dow>& mesh = ir_mesh.regularMesh();
00513     new_property_id(_pid_loading); 
00514     u_int n_ele = mesh.n_geometry(dim);
00515     for (u_int i = 0;i < n_ele;++ i) {
00516       const HGeometry<dim,dow> * p_geo = mesh.template h_geometry<dim>(i);
00517       u_int * p_loading = p_geo->get_property(_pid_loading);
00518       if (p_loading == NULL) {
00519         p_loading = p_geo->new_property(_pid_loading);
00520       }
00521       (*p_loading) = (*value)(mesh.geometry(dim,i));
00522     }
00523   }
00524 
00525   static u_int _default_element_loading(const GeometryBM&) { return 1; }
00527 
00529   public:
00537   void partition(u_int n_new_rank = 0,
00538                  double tol_percent = 0.05,
00539                  bool is_renumerate_root = false,
00540                  void (*f)(double, double, double,
00541                            double&,double&,double&) = NULL) {
00542     int rank =  _forest->rank();
00543     int n_rank = _forest->n_rank();
00544 
00545     u_int loading = this->lump_loading();
00546     u_int total_loading, partial_loading;
00547     MPI_Comm comm = _forest->communicator();
00548     MPI_Scan(&loading, &partial_loading, 1, MPI_UNSIGNED, MPI_SUM, comm);
00549 
00550     if (rank == n_rank - 1) {
00551       total_loading = partial_loading;
00552     }
00553     if (n_new_rank == 0) n_new_rank = n_rank; 
00554 
00555     MPI_Bcast(&total_loading, 1, MPI_UNSIGNED, n_rank - 1, comm);
00556     double mean_loading = n_new_rank;
00557     mean_loading = total_loading/mean_loading;
00558 
00562     this->refine_root_element(mean_loading*tol_percent, is_renumerate_root, f);
00563 
00564     _cut_point.clear();
00565     _new_rank.clear();
00566 
00567     u_int idx = 0;
00568     _cut_point.push_back(idx); 
00569 
00570     loading  = partial_loading - loading;
00571     u_int current_rank = (u_int)(loading/mean_loading);
00572     _new_rank.push_back(current_rank);
00573 
00574     typename forest_t::RootIterator
00575     the_ele = _forest->beginRootElement(),
00576     end_ele = _forest->endRootElement();
00577     for (;the_ele != end_ele;++ the_ele) {
00578       idx += 1;
00579       loading += this->get_loading(*the_ele);
00580       u_int the_rank = (u_int)(loading/mean_loading);
00581       the_rank = std::min(the_rank, n_new_rank - 1); 
00582       if (the_rank > current_rank) {
00583         _cut_point.push_back(idx);
00584         _new_rank.push_back(++ current_rank);
00585       }
00586     }
00587     if (_cut_point.back() != idx) {
00588       _cut_point.push_back(idx);
00589     } else {
00590       _new_rank.pop_back();
00591     }
00592 
00593     free_property_id(_pid_loading); 
00594   }
00595 
00596   private:
00601   u_int lump_loading() {
00602     u_int loading = 0;
00603     typename forest_t::RootIterator
00604     the_ele = _forest->beginRootElement(),
00605     end_ele = _forest->endRootElement();
00606     for (;the_ele != end_ele;++ the_ele) {
00607       loading += this->get_loading(*the_ele);
00608     }
00609     return loading;
00610   }
00611 
00612   private:
00613 
00618   template <class GEO> u_int
00619   get_loading(const GEO& geo) const {
00620     if (&geo == NULL) return 0; 
00621     u_int * p_loading = geo.get_property(_pid_loading);
00622     if (p_loading == NULL) {
00623       p_loading = geo.new_property(_pid_loading);
00624       (*p_loading) = 0;
00625       for (u_int i = 0;i < geo.n_child;++ i) {
00626         (*p_loading) += get_loading(*geo.child[i]);
00627       }
00628     }
00629     return (*p_loading);
00630   }
00632 
00634 
00644   void refine_root_element(double tol, bool is_renum_root, 
00645                            void (*f)(double,double,double,
00646                                      double&,double&,double&)) {
00647     _nre.clear();
00648 
00649     typename forest_t::RootIterator
00650     the_ele = _forest->beginRootElement(),
00651     end_ele = _forest->endRootElement();
00652     for (;the_ele != end_ele;++ the_ele) {
00653       this->refine_root_element(tol, *the_ele, _nre);
00654     }
00655     _forest->rootElement().swap(_nre);
00656     if (is_renum_root) {
00657       _forest->renumerateRootElement(f);
00658     }
00659   }
00660 
00661   void refine_root_element(double tol,
00662                            HGeometry<dim,dow>& geo,
00663                            typename forest_t::container_t& nre) const {
00664     u_int * p_loading = geo.get_property(_pid_loading);
00665     u_int loading = 0;
00666     if (p_loading != NULL) loading = *p_loading;
00667     if ((loading > tol) &&  
00668         (geo.isRefined()) && 
00669         (geo.child[0]->get_property(_pid_loading) != NULL)) { 
00670       for (u_int i = 0;i < geo.n_child;++ i) {
00671         this->refine_root_element(tol, *geo.child[i], nre);
00672       }
00673     } else {
00674       nre.push_back(&geo);
00675     }
00676   }
00677             
00682   void coarse_root_element() {
00683     property_id_t<bool> pid_is_found;
00684     new_property_id(pid_is_found);
00685 
00686     _nre.clear();
00687     typename forest_t::RootIterator
00688     the_ele = _forest->beginRootElement(),
00689     end_ele = _forest->endRootElement();
00690     for (;the_ele != end_ele;++ the_ele) {
00691       this->coarse_root_element(pid_is_found, *the_ele, _nre);
00692     }
00693     _forest->rootElement().swap(_nre);
00694   }
00695 
00696   void coarse_root_element(const property_id_t<bool>& pid,
00697                            HGeometry<dim,dow>& geo,
00698                            typename forest_t::container_t& nre) const {
00699     if (geo.get_property(pid) != NULL) return;
00700     if ((geo.parent != NULL) &&
00701         (_forest->get_shared_info(*geo.parent) == NULL)) {
00702       this->coarse_root_element(pid, *geo.parent, nre);
00703     } else {
00704       nre.push_back(&geo);
00705       geo.new_property(pid);
00706     }
00707   }
00709 
00710 
00712   private:
00718   void set_new_rank() {
00719     new_property_id(_pid_rank_map);
00720 
00721     set_new_rank_down();
00722     set_new_rank_up();
00723 
00724     if (_forest->n_rank() <= 1) return;
00725 
00729     Shared_type_filter::only<0> type_filter;
00730     MPI_Comm comm = _forest->communicator();
00731 #define SYNC_DATA(D) \
00732     if (dim >= D) { \
00733       sync_data(comm, _forest->template get_shared_list<D>(), *this, \
00734                 &this_t::template pack_set_new_rank<D>, \
00735                 &this_t::template unpack_set_new_rank<D>, \
00736                 type_filter); \
00737     }
00738     SYNC_DATA(0);
00739     SYNC_DATA(1);
00740     SYNC_DATA(2);
00741     SYNC_DATA(3);
00742 #undef SYNC_DATA
00743 
00744     //check_new_rank(); /// 调试时用
00745   }
00746 
00747   public:
00748   template <int GDIM> void
00749   pack_set_new_rank(HGeometry<GDIM,dow> * geo,
00750                     int remote_rank,
00751                     AFEPack::ostream<>& os) {
00752     assert ((_forest->get_shared_info(*geo) != NULL));
00753 
00754     rank_map_t * p_map = get_rank_map(*geo);
00755     if (p_map == NULL) {
00756       p_map = new_rank_map(*geo);
00757     }
00758 
00759     std::size_t n = p_map->size();
00760     os << n;
00761     typename rank_map_t::iterator
00762     the_pair = p_map->begin(),
00763     end_pair = p_map->end();
00764     for (;the_pair != end_pair;++ the_pair) {
00765       const int& new_rank = the_pair->first;
00766       const int& old_rank = the_pair->second.old_rank;
00767       int& priority = the_pair->second.priority;
00768       if (is_refined_on_new_rank(*geo, new_rank)) {
00769         priority = 2;
00770       } else if (_forest->is_dummy(*geo)) {
00771         priority = 0;
00772       } else {
00773         priority = 1;
00774       }
00775 
00776       os << new_rank << old_rank << priority;
00777     }
00778   }
00779 
00780   template <int GDIM> void
00781   unpack_set_new_rank(HGeometry<GDIM,dow> * geo,
00782                       int remote_rank,
00783                       AFEPack::istream<>& is) {
00784     assert ((_forest->get_shared_info(*geo) != NULL));
00785 
00794     rank_map_t * p_map = get_rank_map(*geo);
00795     assert (p_map != NULL);
00796 
00797     int new_rank;
00798     std::size_t n;
00799     is >> n;
00800     for (u_int i = 0;i < n;++ i) {
00801       rank_map_val remote_rank_map;
00802       is >> new_rank 
00803          >> remote_rank_map.old_rank
00804          >> remote_rank_map.priority;
00805       typename rank_map_t::iterator 
00806         the_pair = p_map->find(new_rank);
00807       if (the_pair == p_map->end()) {
00808         (*p_map)[new_rank] = remote_rank_map;
00809       } else {
00810         if (the_pair->second.priority == remote_rank_map.priority) {
00812           if (the_pair->second.old_rank > remote_rank_map.old_rank) {
00813             the_pair->second.old_rank = remote_rank_map.old_rank;
00814           }
00815         } else if (remote_rank_map.priority > the_pair->second.priority) {
00816           the_pair->second = remote_rank_map;
00817         }
00818       }
00819     }
00820   }
00821 
00822   private:
00827   void check_new_rank() {
00828     Shared_type_filter::only<0> type_filter;
00829     MPI_Comm comm = _forest->communicator();
00830 
00831 #define SYNC_DATA(D) \
00832     if (dim >= D) { \
00833       sync_data(comm, _forest->template get_shared_list<D>(), *this,    \
00834                 &this_t::template pack_check_new_rank<D>, \
00835                 &this_t::template unpack_check_new_rank<D>, \
00836                 type_filter); \
00837     }
00838     SYNC_DATA(0);
00839     SYNC_DATA(1);
00840     SYNC_DATA(2);
00841     SYNC_DATA(3);
00842 #undef SYNC_DATA
00843   }
00844 
00845   public:
00846   template <int GDIM> void
00847   pack_check_new_rank(HGeometry<GDIM,dow> * geo,
00848                       int remote_rank,
00849                       AFEPack::ostream<>& os) {
00850     rank_map_t * p_map = get_rank_map(*geo);
00851     assert (p_map != NULL);
00852 
00853     os << p_map->size();
00854     typename rank_map_t::iterator
00855       the_pair = p_map->begin(),
00856       end_pair = p_map->end();
00857     for (;the_pair != end_pair;++ the_pair) {
00858       const int& new_rank = the_pair->first;
00859       const int& old_rank = the_pair->second.old_rank;
00860       int& priority = the_pair->second.priority;
00861       os << new_rank << old_rank << priority;
00862     }
00863   }
00864 
00865   template <int GDIM> void
00866   unpack_check_new_rank(HGeometry<GDIM,dow> * geo,
00867                         int remote_rank,
00868                         AFEPack::istream<>& is) {
00869     assert ((_forest->get_shared_info(*geo) != NULL));
00870 
00871     rank_map_t * p_map = get_rank_map(*geo);
00872     assert (p_map != NULL);
00873 
00874     std::size_t n;
00875     is >> n;
00876 
00877     rank_map_t remote_rank_map;
00878     for (u_int i = 0;i < n;++ i) {
00879       int new_rank;
00880       rank_map_val remote_rank_map_val;
00881       is >> new_rank 
00882          >> remote_rank_map_val.old_rank
00883          >> remote_rank_map_val.priority;
00884       remote_rank_map[new_rank] = remote_rank_map_val;
00885     }
00886 
00887     typename rank_map_t::iterator
00888       the_remote_rank_map = remote_rank_map.begin(),
00889       end_remote_rank_map = remote_rank_map.end();
00890     for (;the_remote_rank_map != end_remote_rank_map;++ the_remote_rank_map) {
00891       const int& remote_new_rank = the_remote_rank_map->first;
00892       const int& remote_old_rank = the_remote_rank_map->second.old_rank;
00893       const int& remote_priority = the_remote_rank_map->second.priority;
00894 
00895       typename rank_map_t::iterator 
00896         the_pair = p_map->find(remote_new_rank);
00897       int errcode = 1;
00898       if (the_pair == p_map->end()) {
00899         errcode *= 2;
00900       } else {
00901         if (the_pair->second.old_rank != remote_old_rank) {
00902           errcode *= 3;
00903         }
00904         if (the_pair->second.priority != remote_priority) {
00905           errcode *= 5;
00906         }
00907       }
00908       if (errcode != 1) {
00909         std::cerr << " [errcode:" << errcode << "]";
00910         print_geo_rank_map(*geo);
00911         std::cerr << "-[" << remote_rank 
00912                   << ":" << GDIM;
00913         print_rank_map(&remote_rank_map);
00914         std::cerr << "] ";
00915         MPI_Abort(_forest->communicator(), errcode);
00916       }
00917     }
00918   }
00919 
00920   private:
00925   template <class GEO> void
00926     geometry_set_new_rank_down(const GEO& geo, int new_rank) const {
00927     rank_map_t * p_map = get_rank_map(geo);
00928     if (p_map == NULL) {
00929       p_map = new_rank_map(geo);
00930     } else if (p_map->find(new_rank) != p_map->end()) {
00931       return; 
00932     }
00933 
00934     (*p_map)[new_rank].old_rank = _forest->rank(); 
00935     
00937     for (u_int i = 0;i < geo.n_vertex;++ i) {
00938       geometry_set_new_rank_down(*geo.vertex[i], new_rank);
00939     }
00940     for (u_int i = 0;i < geo.n_boundary;++ i) {
00941       geometry_set_new_rank_down(*geo.boundary[i], new_rank);
00942     }
00943     if (geo.isRefined()) {
00944       for (u_int i = 0;i < geo.n_child;++ i) {
00945         geometry_set_new_rank_down(*geo.child[i], new_rank);
00946       }
00947     }
00948   }
00949 
00953   void set_new_rank_down() {
00954     typename forest_t::RootIterator
00955     the_ele = _forest->beginRootElement();
00956     u_int n_part = _new_rank.size();
00957     for (u_int i = 0;i < n_part;++ i) {
00958       for (u_int j = _cut_point[i];j < _cut_point[i + 1];++ j, ++ the_ele) {
00959         geometry_set_new_rank_down(*the_ele, _new_rank[i]);
00960       }
00961     }
00962     assert (the_ele == _forest->endRootElement());
00963   }
00964 
00971   template <class GEO> void
00972     geometry_set_new_rank_up(const GEO& geo, int new_rank) const {
00973     rank_map_t * p_map = get_rank_map(geo);
00974     if (p_map == NULL) p_map = new_rank_map(geo);
00975 
00976     (*p_map)[new_rank].old_rank = _forest->rank(); 
00977 
00979     for (u_int i = 0;i < geo.n_vertex;++ i) {
00980       geometry_set_new_rank_up(*geo.vertex[i], new_rank);
00981     }
00982     for (u_int i = 0;i < geo.n_boundary;++ i) {
00983       geometry_set_new_rank_up(*geo.boundary[i], new_rank);
00984     }
00985     GEO* const& p_parent = geo.parent;
00986     if (p_parent != NULL) {
00987       if (! is_on_this_new_rank(*p_parent, new_rank)) {
00988         geometry_set_new_rank_up(*p_parent, new_rank);
00989       }
00990 
00991       for (u_int i = 0;i < p_parent->n_child;++ i) {
00992         GEO* const& p_sibling = p_parent->child[i];
00993         if (p_sibling == &geo) continue; 
00994         if (! is_on_this_new_rank(*p_sibling, new_rank)) {
00995           geometry_set_new_rank_up(*p_sibling, new_rank);
00996         }
00997       }
00998     }
00999   }
01000 
01004   void set_new_rank_up() {
01005     typename forest_t::RootIterator
01006     the_ele = _forest->beginRootElement();
01007     u_int n_part = _new_rank.size();
01008     for (u_int i = 0;i < n_part;++ i) {
01009       for (u_int j = _cut_point[i];j < _cut_point[i + 1];++ j, ++ the_ele) {
01010         geometry_set_new_rank_up(*the_ele, _new_rank[i]);
01011       }
01012     }
01013     assert (the_ele == _forest->endRootElement());
01014   }
01016 
01018   private:
01034   void global_index() {
01038     new_property_id(_pid_global_idx);
01039     unsigned long idx = 0;
01040     typename forest_t::RootIterator
01041     the_ele = _forest->beginRootElement(),
01042     end_ele = _forest->endRootElement();
01043     for (;the_ele != end_ele;++ the_ele) {
01044       HGeometry<dim,dow> * p_geo = &(*the_ele);
01045       do {
01046         if (p_geo->parent != NULL) {
01047           p_geo = p_geo->parent;
01048         } else break;
01049       } while (true);
01050       geometry_global_index(*p_geo, idx);
01051     }
01052     free_property_id(_pid_global_idx); 
01053 
01057     MPI_Comm comm = _forest->communicator();
01058     unsigned long global_idx = idx;
01059     MPI_Scan(&idx, &global_idx, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
01060     idx = global_idx - idx;
01061 
01065     new_property_id(_pid_global_idx);
01066     the_ele = _forest->beginRootElement();
01067     for (;the_ele != end_ele;++ the_ele) {
01068       HGeometry<dim,dow> * p_geo = &(*the_ele);
01069       do {
01070         if (p_geo->parent != NULL) {
01071           p_geo = p_geo->parent;
01072         } else break;
01073       } while (true);
01074       geometry_global_index(*p_geo, idx);
01075     }
01076     
01077     if (_forest->n_rank() <= 1) return;
01082     Shared_type_filter::only<0> type_filter;
01083 #define SYNC_DATA(D) \
01084     if (dim >= D) { \
01085       sync_data(comm, _forest->template get_shared_list<D>(), *this, \
01086                 &this_t::template pack_global_index<D>, \
01087                 &this_t::template unpack_global_index<D>, \
01088                 type_filter); \
01089     }
01090     SYNC_DATA(0);
01091     SYNC_DATA(1);
01092     SYNC_DATA(2);
01093     SYNC_DATA(3);
01094 #undef SYNC_DATA
01095   }
01096 
01097   public:
01098   template <int GDIM> void
01099   pack_global_index(HGeometry<GDIM,dow> * geo,
01100                     int remote_rank,
01101                     AFEPack::ostream<>& os) {
01102     unsigned long * p_idx = get_global_idx(*geo);
01103     assert (p_idx != NULL);
01104 
01105     os << *p_idx; 
01106   }
01107 
01113   template <int GDIM> void
01114   unpack_global_index(HGeometry<GDIM,dow> * geo,
01115                       int remote_rank,
01116                       AFEPack::istream<>& is) {
01117     assert ((_forest->get_shared_info(*geo) != NULL));
01118 
01119     unsigned long * p_idx = get_global_idx(*geo);
01120     assert (p_idx != NULL);
01121 
01122     unsigned long remote_idx;
01123     is >> remote_idx;
01124     if (*p_idx > remote_idx) {
01125       *p_idx = remote_idx; 
01126     }
01127   }
01128 
01129   private:
01133   template <class GEO> void
01134   geometry_global_index(GEO& geo, 
01135                         unsigned long& idx) const {
01136     unsigned long * p_idx = get_global_idx(geo);
01137     if (p_idx != NULL) return; 
01138 
01139     rank_map_t * p_map = get_rank_map(geo);
01140     if (_forest->get_shared_info(geo) != NULL || 
01141         (p_map != NULL && p_map->size() > 1)) { 
01142       p_idx = new_global_idx(geo);
01143       *p_idx = idx ++;
01144     }
01145 
01147     for (u_int i = 0;i < geo.n_vertex;++ i) {
01148       geometry_global_index(*geo.vertex[i], idx);
01149     }
01150     for (u_int i = 0;i < geo.n_boundary;++ i) {
01151       geometry_global_index(*geo.boundary[i], idx);
01152     }
01153     if (geo.isRefined()) {
01154       for (u_int i = 0;i < geo.n_child;++ i) {
01155         geometry_global_index(*geo.child[i], idx);
01156       }
01157     }
01158   }
01160 
01162   public:
01163   void write_config_file(const std::string& dirname) {
01164     char filename[1024];
01165     sprintf(filename, "%s/.config", dirname.c_str());
01166     std::ofstream os(filename);
01167     os << _forest->n_rank() << "\t# number of old rank" << std::endl;
01168     os.close();
01169     Migration::save_config(dirname);
01170   }
01171 
01172   private:
01173   void export_birdview(birdview_t& ir_mesh, 
01174                        Migration::data_id_t data_id, 
01175                        u_int idx) {
01176     typename birdview_t::ActiveIterator
01177     the_ele = ir_mesh.beginActiveElement(),
01178     end_ele = ir_mesh.endActiveElement();
01179     for (;the_ele != end_ele;++ the_ele) {
01180       HGeometry<dim,dow> * p_geo = the_ele->h_element;
01181       AFEPack::ostream<> os(p_geo->buffer[data_id]);
01182       os << idx;
01183     }
01184   }
01185 
01186   std::vector<u_int> _n_orphans;
01187   std::vector<std::list<void *> > _orphans_ptr;
01188 
01189   template <class GEO> void
01190     geometry_collect_orphans(GEO& geo, 
01191                              const property_id_t<>& pid,
01192                              int new_rank) {
01193     if (geo.get_property(pid) != NULL) return;
01194     geo.new_property(pid);
01195 
01196     if ((geo.get_property(_pid_is_saved) == NULL) &&
01197         is_save_on_this_rank(geo, new_rank)) {
01198       _n_orphans[GEO::dim] += 1;
01199       _orphans_ptr[GEO::dim].push_back((void *)(&geo));
01200     }
01201 
01202     for (u_int i = 0;i < geo.n_vertex;++ i) {
01203       geometry_collect_orphans(*geo.vertex[i], pid, new_rank);
01204     }
01205     for (u_int i = 0;i < geo.n_boundary;++ i) {
01206       geometry_collect_orphans(*geo.boundary[i], pid, new_rank);
01207     }
01208     if (geo.isRefined()) {
01209       for (u_int i = 0;i < geo.n_child;++ i) {
01210         geometry_collect_orphans(*geo.child[i], pid, new_rank);
01211       }
01212     }
01213   }    
01214 
01215   void collect_orphans(int new_rank) {
01216     _n_orphans.clear();
01217     _n_orphans.resize(dim + 1, 0);
01218 
01219     _orphans_ptr.clear();
01220     _orphans_ptr.resize(dim + 1);
01221 
01222     property_id_t<> pid_is_collected;
01223     new_property_id(pid_is_collected);
01224 
01225     typename forest_t::RootIterator
01226     the_ele = _forest->beginRootElement(),
01227     end_ele = _forest->endRootElement();
01228     for (;the_ele != end_ele;++ the_ele) {
01229       HGeometry<dim,dow> * p_geo = &(*the_ele);
01230       do {
01231         if (p_geo->parent != NULL) {
01232           p_geo = p_geo->parent;
01233         } else break;
01234       } while (true);
01235       geometry_collect_orphans(*p_geo, pid_is_collected, new_rank);
01236     }
01237   }
01238 
01239   template <int D> void 
01240     save_dim_orphans(boost::archive::HGeometry_oarchive<this_t>& oa) {
01241     typedef HGeometry<D,dow> GEO;
01242     u_int n_orphans = _n_orphans[D];
01243 #if 1
01244     std::cerr << D << "D("  << n_orphans << ") ";
01245 #endif
01246     assert (n_orphans == _orphans_ptr[D].size());
01247     oa & boost::serialization::make_binary_object(&n_orphans, sizeof(u_int));
01248     typename std::list<void *>::iterator
01249       the_ptr = _orphans_ptr[D].begin(), 
01250       end_ptr = _orphans_ptr[D].end();
01251     for (;the_ptr != end_ptr;++ the_ptr) {
01252       GEO * p_geo = (GEO *)(*the_ptr);
01253       oa & p_geo;
01254     }
01255   }
01256 
01257   void save_orphans(boost::archive::HGeometry_oarchive<this_t>& oa) {
01258 #if 1
01259     std::cerr << "saving orphans: ";
01260 #endif
01261     collect_orphans(oa.new_rank());
01262     if (dim >= 0) this->template save_dim_orphans<0>(oa);
01263     if (dim >= 1) this->template save_dim_orphans<1>(oa);
01264     if (dim >= 2) this->template save_dim_orphans<2>(oa);
01265     if (dim >= 3) this->template save_dim_orphans<3>(oa);
01266 #if 1
01267     std::cerr << std::endl;
01268 #endif
01269   }
01270 
01271   public:
01272   void save_data(const std::string& dirname, 
01273                  birdview_set_t& bvs) {
01274     set_new_rank();
01275     global_index();
01276 
01277     int dummy;
01278     char command[1024], *filename = command;
01279 
01280     Migration::data_id_t id = Migration::name_to_id("__internal_birdviewflag__"); 
01281     if (! Migration::is_valid(id)) {
01282       id = Migration::register_data_name("__internal_birdviewflag__");
01283     }
01284     typename birdview_set_t::iterator
01285     the_bv = bvs.begin(), end_bv = bvs.end();
01286     for (u_int idx = 0;the_bv != end_bv;++ the_bv) {
01287       export_birdview(*the_bv, id, idx ++);
01288     }
01289     Migration::sync_data_buffer(*_forest);
01290 
01291     typedef boost::archive::HGeometry_oarchive<this_t> archive_t;
01292     typename forest_t::RootIterator
01293     the_ele = _forest->beginRootElement();
01294     u_int n_part = _new_rank.size();
01295     for (u_int i = 0;i < n_part;++ i) {
01296       sprintf(command, "mkdir -p %s && mkdir -p %s/%d", 
01297               dirname.c_str(), dirname.c_str(), _new_rank[i]);
01298       dummy = system(command);
01299       sprintf(filename, "%s/%d/%d.dat", dirname.c_str(), 
01300               _new_rank[i], _forest->rank());
01301       std::ofstream os(filename, std::ios::binary);
01302 
01303       {
01304         archive_t oa(*this, _new_rank[i], os);
01305         new_property_id(_pid_is_saved);
01306 
01307         u_int n_ele = _cut_point[i + 1] - _cut_point[i];
01308         oa & boost::serialization::make_binary_object(&n_ele, sizeof(u_int));
01309         std::cerr << "saving data in " << filename 
01310                   << ", # macro element = " << n_ele << std::endl;
01311 
01312         for (u_int j = _cut_point[i];j < _cut_point[i + 1];++ j) {
01313           HGeometry<dim,dow> * p_ele = *(the_ele ++);
01314           oa & p_ele;
01315         }
01316 
01317         save_orphans(oa);
01318         free_property_id(_pid_is_saved);
01319       }
01320 
01321       os.close();
01322     }
01323 
01324     if (_forest->rank() == 0) write_config_file(dirname);
01325     _forest->rootElement().swap(_nre); 
01326     _nre.clear();
01327     MPI_Barrier(_forest->communicator()); 
01328   }
01330 
01332   public:
01333   int load_config_file(const std::string& dirname) {
01334     Migration::load_config(dirname);
01335 
01336     char filename[1024];
01337     sprintf(filename, "%s/.config", dirname.c_str());
01338     filtering_istream is;
01339     Migration::ensured_open_filtering_stream(filename, is);
01340 
01341     int old_rank;
01342     is >> old_rank;
01343     return old_rank;
01344   }
01345 
01346   private:
01347   void reconstruct_birdview(HElement<dim,dow>& ele,
01348                             Migration::data_id_t data_id,
01349                             u_int idx) {
01350     typedef Migration::buffer_t buffer_t;
01351     buffer_t& hbuf = ele.h_element->buffer;
01352     typename buffer_t::iterator it = hbuf.find(data_id);
01353     bool is_active = true;
01354     ele.value = 0;
01355     if (it == hbuf.end()) {
01356       is_active = false;
01357     } else {
01358       BinaryBuffer<>& buf = it->second;
01359       int n = buf.size()/sizeof(u_int), i;
01360       u_int idx1;
01361       AFEPack::istream<> is(buf);
01362       for (i = 0;i < n;++ i) {
01363         is >> idx1;
01364         if (idx1 == idx) break;
01365       }
01366       if (i == n) is_active = false;
01367     }
01368     if (! is_active) {
01369       assert (ele.h_element->isRefined());
01370       ele.value = 1;
01371       ele.refine();
01372       for (int i = 0;i < ele.n_child;++ i) {
01373         reconstruct_birdview(*ele.child[i], data_id, idx);
01374       }
01375     }
01376   }
01377 
01378   void reconstruct_birdview(birdview_t& ir_mesh,
01379                             Migration::data_id_t data_id,
01380                             u_int idx) {
01381     typename birdview_t::RootIterator
01382     the_ele = ir_mesh.beginRootElement(),
01383     end_ele = ir_mesh.endRootElement();
01384     for (;the_ele != end_ele;++ the_ele) {
01385       reconstruct_birdview(*the_ele, data_id, idx);
01386     }
01387   }
01388 
01393   template <class GEO>
01394     void set_dummy_flag(GEO& geo) const {
01395     if ((_forest->get_shared_info(geo) != NULL) &&
01396         (! _forest->is_dummy(geo))) {
01397       _forest->dummy(geo);
01398     }
01399 
01400     for (u_int i = 0;i < geo.n_boundary;++ i) {
01401       this->set_dummy_flag(*geo.boundary[i]);
01402     }
01403     if (geo.isRefined()) {
01404       for (u_int i = 0;i < geo.n_child;++ i) {
01405         this->set_dummy_flag(*geo.child[i]);
01406       }
01407     }
01408   }
01409 
01413   template <class GEO>
01414     void clear_dummy_flag(GEO& geo) const {
01415     if (_forest->is_dummy(geo)) {
01416       _forest->undummy(geo);
01417     }
01418 
01419     for (u_int i = 0;i < geo.n_boundary;++ i) {
01420       this->clear_dummy_flag(*geo.boundary[i]);
01421     }
01422     if (geo.isRefined()) {
01423       for (u_int i = 0;i < geo.n_child;++ i) {
01424         this->clear_dummy_flag(*geo.child[i]);
01425       }
01426     }
01427   }
01428 
01439   void set_dummy_flag() {
01440     typename forest_t::RootIterator
01441       the_ele = _forest->beginRootElement(),
01442       end_ele = _forest->endRootElement();
01443     for (;the_ele != end_ele;++ the_ele) {
01444       const HGeometry<dim,dow> * p_geo = &(*the_ele);
01445       do {
01446         if (p_geo->parent != NULL) {
01447           p_geo = p_geo->parent;
01448         } else break;
01449       } while (true);
01450       set_dummy_flag(*p_geo);
01451     }
01452 
01453     the_ele = _forest->beginRootElement();
01454     for (;the_ele != end_ele;++ the_ele) {
01455       clear_dummy_flag(*the_ele);
01456     }
01457   }
01458 
01459   private:
01460   template <int D>
01461     void load_dim_orphans(boost::archive::HGeometry_iarchive<this_t>& ia) const {
01462     typedef HGeometry<D,dow> GEO;
01463     u_int n_orphans;
01464     ia & boost::serialization::make_binary_object(&n_orphans, sizeof(u_int));
01465 #if 1
01466     std::cerr << D << "D("  << n_orphans << "), ";
01467 #endif
01468     for (u_int i = 0;i < n_orphans;++ i) {
01469       GEO * p_geo;
01470       ia & p_geo;
01471     }
01472   }
01473 
01474   void load_orphans(boost::archive::HGeometry_iarchive<this_t>& ia) const {
01475 #if 1
01476     std::cerr << "loading orphans: ";
01477 #endif
01478     if (dim >= 0) this->template load_dim_orphans<0>(ia);
01479     if (dim >= 1) this->template load_dim_orphans<1>(ia);
01480     if (dim >= 2) this->template load_dim_orphans<2>(ia);
01481     if (dim >= 3) this->template load_dim_orphans<3>(ia);
01482 #if 1
01483     std::cerr << std::endl;
01484 #endif
01485   }
01486 
01488   template <class GEO> void
01489     geometry_check_integrity(const GEO& geo) const {
01493     if (geo.parent != NULL) {
01494       u_int i = 0;
01495       for (;i < geo.parent->n_child;++ i) {
01496         if (&geo == geo.parent->child[i]) break;
01497       }
01498       assert ((i < geo.parent->n_child));
01499     }
01500 
01504     for (u_int i = 0;i < geo.n_boundary;++ i) {
01505       assert ((geo.boundary[i] != NULL));
01506       geometry_check_integrity(*geo.boundary[i]);
01507     }
01508     if (geo.isRefined()) {
01513       for (u_int i = 0;i < geo.n_child;++ i) {
01514         assert ((geo.child[i] != NULL));
01515         assert ((&geo == geo.child[i]->parent));
01516         geometry_check_integrity(*geo.child[i]);
01517       }
01518     } else {
01522       for (u_int i = 0;i < geo.n_child;++ i) {
01523         assert (geo.child[i] == NULL);
01524       }
01525     }
01526   }
01527 
01528   public:
01529   void check_integrity() const {
01530     typename forest_t::RootIterator
01531     the_ele = _forest->beginRootElement(),
01532     end_ele = _forest->endRootElement();
01533     for (;the_ele != end_ele;++ the_ele) {
01534       HGeometry<dim,dow> * p_geo = &(*the_ele);
01535       do {
01536         if (p_geo->parent != NULL) {
01537           p_geo = p_geo->parent;
01538         } else break;
01539       } while (true);
01540       geometry_check_integrity(*p_geo);
01541     }
01542   }
01543 
01544   public:
01545   void load_data(const std::string& dirname,
01546                  birdview_set_t& bvs,
01547                  bool is_has_orphans = false) {
01548     clear();
01549 
01550     if (rank() == 6)
01551       std::cerr << "is_has_orphans=" << is_has_orphans << "\n";
01552 
01553     int n_old_rank = load_config_file(dirname);
01554 
01556     MPI_Barrier(_forest->communicator()); 
01557 
01559     char filename[1024];
01560 
01561     typedef boost::archive::HGeometry_iarchive<this_t> archive_t;
01562 
01563     int rank = 0;
01564     do {
01565       std::ifstream is;
01566       do {
01567         sprintf(filename, "%s/%d/%d.dat", dirname.c_str(), 
01568                 _forest->rank(), rank ++);
01569         if (rank > n_old_rank) {
01570           if (! _global_pointer_to_merge_later.empty()) {
01571             std::cerr << "Rank " << _forest->rank() << ": ";
01572             std::map<unsigned long, std::list<void *> >::iterator
01573               the_entry = _global_pointer_to_merge_later.begin(),
01574               end_entry = _global_pointer_to_merge_later.end();
01575             for (;the_entry != end_entry;++ the_entry) {
01576               std::cerr << "(" << the_entry->first 
01577                         << "->" << the_entry->second.size()
01578                 << ")";
01579             }
01580             std::cerr << std::endl;
01581             assert (false);
01582           }
01583 
01584           this->share_global_pointer();
01585           this->coarse_root_element();
01586           this->set_dummy_flag();
01587           _forest->set_is_used_up();
01588 
01589           //_forest->verify_shared_object();
01590           //this->check_integrity(); /// 调试时候使用
01591 
01593           Migration::data_id_t id = Migration::name_to_id("__internal_birdviewflag__"); 
01594           typename birdview_set_t::iterator
01595             the_bv = bvs.begin(), end_bv = bvs.end();
01596           for (u_int idx = 0;the_bv != end_bv;++ the_bv) {
01597             the_bv->reinit(*_forest);
01598             reconstruct_birdview(*the_bv, id, idx ++);
01599           }
01600           return;
01601         }
01602         is.open(filename, std::ios::binary); 
01603         if (is.good()) break; 
01604       } while (1);
01605 
01606       {
01607         u_int n_ele;
01608         std::cerr << "loading data from " << filename;
01609         archive_t ia(*this, _forest->rank(), is); 
01610         ia & boost::serialization::make_binary_object(&n_ele, sizeof(u_int));
01611         std::cerr << ", # macro element = " << n_ele << std::endl;
01612 
01613         for (u_int i = 0;i < n_ele;++ i) {
01614           HGeometry<dim,dow> * p_ele;
01615           ia & p_ele;
01616           _forest->rootElement().push_back(p_ele);
01617         }
01618 
01619         if (is_has_orphans) load_orphans(ia);
01620       }
01621 
01622       is.close();
01623     } while (1);
01624 
01625   }
01626 
01627   };
01628 
01629   template <class FOREST>
01630     void load_forest(const std::string& dirname,
01631                      FOREST& forest,
01632                      bool is_has_orphans = false) {
01633     HLoadBalance<FOREST> hlb(forest);
01634     BirdViewSet<FOREST> bvs;
01635     hlb.load_data(dirname, bvs, is_has_orphans);
01636   }
01637 
01638   template <class FOREST>
01639     void load_mesh(const std::string& dirname,
01640                    FOREST& forest,
01641                    BirdView<FOREST>& mesh,
01642                    bool is_has_orphans = false) {
01643     HLoadBalance<FOREST> hlb(forest);
01644     BirdViewSet<FOREST> bvs;
01645     bvs.add(mesh);
01646     hlb.load_data(dirname, bvs, is_has_orphans);
01647   }
01648 
01649   template <class FOREST>
01650     void load_mesh(const std::string& dirname,
01651                    bool is_has_orphans,
01652                    FOREST& forest,
01653                    u_int n_mesh, ...) {
01654     HLoadBalance<FOREST> hlb(forest);
01655     typedef BirdView<FOREST> * mesh_ptr_t;
01656 
01657     BirdViewSet<FOREST> bvs;
01658     va_list ap;
01659     va_start(ap, n_mesh);
01660     for (u_int i = 0;i < n_mesh;++ i) {
01661       mesh_ptr_t p_mesh = va_arg(ap, mesh_ptr_t);
01662       bvs.add(*p_mesh);
01663     }
01664     va_end(ap);
01665 
01666     hlb.load_data(dirname, bvs, is_has_orphans);
01667   }
01668 
01669   template <class FOREST>
01670     void load_mesh_set(const std::string& dirname,
01671                        FOREST& forest,
01672                        BirdViewSet<FOREST>& bvs, 
01673                        bool is_has_orphans = false) {
01674     HLoadBalance<FOREST> hlb(forest);
01675     hlb.load_data(dirname, bvs, is_has_orphans);
01676   }
01677 
01685   void lb_collect_local_data_dir(MPI_Comm comm,
01686                                  const std::string& src_dir,
01687                                  const std::string& dst_dir);
01693   void lb_sync_local_data_dir(MPI_Comm comm,
01694                               const std::string& src_dir,
01695                               const std::string& dst_dir);
01696 }
01697 
01698 #endif // __MPI_ULoadBalance_h__
01699