AFEPack
MPI_HGeometry.templates.h
浏览该文件的文档。
00001 
00011 #define TEMPLATE template <int DIM, int DOW, class MATCHER>
00012 #define THIS HGeometryForest<DIM,DOW,MATCHER>
00013 
00014 TEMPLATE void 
00015 THIS::readMesh(const std::string& filename) {
00016 
00020   const int& rank = this->_rank;
00021   const int& n_rank = this->_n_rank;
00022 
00031   char filename_buf[256];
00032   if (n_rank > 1) {
00033     sprintf(filename_buf, "%s%d.mesh", filename.c_str(), rank);
00034   } else {
00035     sprintf(filename_buf, "%s.mesh", filename.c_str());
00036   }
00037   std::cerr << "Reading in mesh data file " 
00038             << filename_buf 
00039             << " as geometry tree root ..." 
00040             << std::endl;
00041   std::ifstream is(filename_buf);
00042 
00044   u_int n_point;
00045   is >> n_point;
00046   std::cerr << "\t# points: " << n_point << std::endl;
00047   std::vector<afepack::Point<this_t::dow> > point(n_point);
00048   for (u_int i = 0;i < n_point;i ++) is >> point[i];
00049   is >> n_point;
00050   std::vector<HGeometry<0,this_t::dow> *> geo_0d(n_point, (HGeometry<0,this_t::dow> *)NULL);
00051   for (u_int i = 0;i < n_point;i ++) {
00052     u_int j, k;
00053     is >> j; geo_0d[j] = new HGeometry<0,this_t::dow>();
00054     is >> k >> k; *dynamic_cast<afepack::Point<this_t::dow> *>(geo_0d[j]) = point[k];
00055     is >> k >> k >> geo_0d[j]->bmark;
00056   }
00057   point.clear();
00058         
00059 #define GDIM 1
00060   u_int n_geo_1d;
00061   std::vector<HGeometry<GDIM,this_t::dow>*> geo_1d;
00062   if (DIM >= 1) {
00063     is >> n_geo_1d;
00064     std::cerr << "\t# 1D-geometry: " << n_geo_1d << std::endl;
00065     geo_1d.resize(n_geo_1d, NULL);
00066     for (u_int i = 0;i < n_geo_1d;i ++) {
00067       u_int j, k, l;
00068       is >> j >> k; geo_1d[j] = new HGeometry<GDIM,this_t::dow>();
00069       for (k = 0;k < GDIM + 1;k ++) {
00070         is >> l; geo_1d[j]->vertex[k] = geo_0d[l];
00071       }
00072       is >> k;
00073       for (k = 0;k < GDIM + 1;k ++) {
00074         is >> l;
00075       }
00076       is >> geo_1d[j]->bmark;
00077     }
00078   }
00079 #undef GDIM
00080 
00081 #define GDIM 2
00082   u_int n_geo_2d;
00083   std::vector<HGeometry<GDIM,this_t::dow>*> geo_2d;
00084   if (DIM >= 2) {
00085     is >> n_geo_2d;
00086     std::cerr << "\t# 2D-geometry: " << n_geo_2d << std::endl;
00087     geo_2d.resize(n_geo_2d, NULL);
00088     for (u_int i = 0;i < n_geo_2d;i ++) {
00089       u_int j, k, l;
00090       is >> j >> k; geo_2d[j] = new HGeometry<GDIM,this_t::dow>();
00091       for (k = 0;k < GDIM + 1;k ++) {
00092         is >> l; geo_2d[j]->vertex[k] = geo_0d[l];
00093       }
00094       is >> k;
00095       for (k = 0;k < GDIM + 1;k ++) {
00096         is >> l; geo_2d[j]->boundary[k] = geo_1d[l];
00097       }
00098       is >> geo_2d[j]->bmark;
00099     }
00100   }
00101 #undef GDIM
00102 
00103 #define GDIM 3
00104   u_int n_geo_3d;
00105   std::vector<HGeometry<GDIM,this_t::dow>*> geo_3d;
00106   if (DIM >= 3) { 
00107     is >> n_geo_3d;
00108     std::cerr << "\t# 3D-geometry: " << n_geo_3d << std::endl;
00109     geo_3d.resize(n_geo_3d, NULL);
00110     for (u_int i = 0;i < n_geo_3d;i ++) {
00111       u_int j, k, l;
00112       is >> j >> k; geo_3d[j] = new HGeometry<GDIM,this_t::dow>();
00113       for (k = 0;k < GDIM + 1;k ++) {
00114         is >> l; geo_3d[j]->vertex[k] = geo_0d[l];
00115       }
00116       is >> k;
00117       for (k = 0;k < GDIM + 1;k ++) {
00118         is >> l; geo_3d[j]->boundary[k] = geo_2d[l];
00119       }
00120       is >> geo_3d[j]->bmark;
00121     }
00122 #undef GDIM
00123   }
00124   is.close();
00125 
00126   if (DIM == 1) {
00127     for (u_int i = 0;i < n_geo_1d;++ i) {
00128       this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_1d[i]);
00129     }
00130   } else if (DIM == 2) {
00131     for (u_int i = 0;i < n_geo_2d;++ i) {
00132       this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_2d[i]);
00133     }
00134   } else if (DIM == 3) {
00135     for (u_int i = 0;i < n_geo_3d;++ i) {
00136       this->rootElement().push_back((HGeometry<DIM,this_t::dow> *)geo_3d[i]);
00137     }
00138   }
00139 
00145   if (n_rank == 1) return; 
00146 
00147   sprintf(filename_buf, "%s%d.share", filename.c_str(), rank);
00148   is.open(filename_buf); 
00149   std::cerr << "Reading in shared data file " 
00150             << filename_buf
00151             << " as shared information ..." 
00152             << std::endl;
00153 
00155   std::vector<bool> proc_flag(n_rank, false);
00156   std::vector<BinaryBuffer<> > proc_buf_in(n_rank), proc_buf_out(n_rank);
00157   std::vector<AFEPack::ostream<> > proc_os(n_rank);
00158   for (u_int i = 0;i < n_rank;++ i) {
00159     proc_os[i].set_buffer(proc_buf_out[i]);
00160   }
00161 
00163   int n_item;
00164   std::vector<std::pair<std::vector<std::pair<int,int> >, int> > 
00165     shared_geo_0d, shared_geo_1d, shared_geo_2d, shared_geo_3d;        
00166   if (DIM >= 0) {
00167     std::vector<int> n_rank_shared_0d(n_rank, 0);
00168     int n_shared_geo_0d;
00169     is >> n_shared_geo_0d >> n_shared_geo_0d;
00170     std::cerr << "\t# 0d-shared geometry: " << n_shared_geo_0d << std::endl;
00171     shared_geo_0d.resize(n_shared_geo_0d);
00172     for (int i = 0;i < n_shared_geo_0d;++ i) {
00173       is >> n_item; 
00174       shared_geo_0d[i].first.resize(n_item);
00175       for (int j = 0;j < n_item;++ j) {
00176         is >> shared_geo_0d[i].first[j].first   
00177            >> shared_geo_0d[i].first[j].second; 
00178         n_rank_shared_0d[shared_geo_0d[i].first[j].first] += 1;
00179       }
00180       is >> shared_geo_0d[i].second; 
00181     }
00182 
00183     for (int i = 0;i < n_rank;++ i) {
00184       proc_os[i] << n_rank_shared_0d[i];
00185       if (n_rank_shared_0d[i] > 0) {
00186         proc_flag[i] = true;
00187       }
00188     }
00189     for (int i = 0;i < n_shared_geo_0d;++ i) {
00190       n_item = shared_geo_0d[i].first.size();
00191       int local_idx;
00192       for (int j = 0;j < n_item;++ j) {
00193         int rnk = shared_geo_0d[i].first[j].first;
00194         int idx = shared_geo_0d[i].first[j].second;
00195         if (rnk == rank) {
00196           local_idx = idx;
00197           break;
00198         }
00199       }
00200       for (int j = 0;j < n_item;++ j) {
00201         int rnk = shared_geo_0d[i].first[j].first;
00202         int idx = shared_geo_0d[i].first[j].second;
00203         if (rnk != rank) {
00204           proc_os[rnk] << idx << geo_0d[local_idx];
00205         }
00206       }
00207     }
00208   }
00209 
00210   if (DIM >= 1) {
00211     std::vector<int> n_rank_shared_1d(n_rank, 0);
00212     int n_shared_geo_1d;
00213     is >> n_shared_geo_1d >> n_shared_geo_1d;
00214     std::cerr << "\t# 1d-shared geometry: " << n_shared_geo_1d << std::endl;
00215     shared_geo_1d.resize(n_shared_geo_1d);        
00216     for (int i = 0;i < n_shared_geo_1d;++ i) {
00217       is >> n_item; 
00218       shared_geo_1d[i].first.resize(n_item);
00219       for (int j = 0;j < n_item;++ j) {
00220         is >> shared_geo_1d[i].first[j].first   
00221            >> shared_geo_1d[i].first[j].second; 
00222         n_rank_shared_1d[shared_geo_1d[i].first[j].first] += 1;
00223       }
00224       is >> shared_geo_1d[i].second; 
00225     }
00226 
00227     for (int i = 0;i < n_rank;++ i) {
00228       proc_os[i] << n_rank_shared_1d[i];
00229       if (n_rank_shared_1d[i] > 0) {
00230         proc_flag[i] = true;
00231       }
00232     }
00233     for (int i = 0;i < n_shared_geo_1d;++ i) {
00234       n_item = shared_geo_1d[i].first.size();
00235       int local_idx;
00236       for (int j = 0;j < n_item;++ j) {
00237         int rnk = shared_geo_1d[i].first[j].first;
00238         int idx = shared_geo_1d[i].first[j].second;
00239         if (rnk == rank) {
00240           local_idx = idx;
00241           break;
00242         }
00243       }
00244       for (int j = 0;j < n_item;++ j) {
00245         int rnk = shared_geo_1d[i].first[j].first;
00246         int idx = shared_geo_1d[i].first[j].second;
00247         if (rnk != rank) {
00248           proc_os[rnk] << idx << geo_1d[local_idx];
00249         }
00250       }
00251     }
00252   }
00253 
00254   if (DIM >= 2) {
00255     std::vector<int> n_rank_shared_2d(n_rank, 0);
00256     int n_shared_geo_2d;
00257     is >> n_shared_geo_2d >> n_shared_geo_2d;
00258     std::cerr << "\t# 2d-shared geometry: " << n_shared_geo_2d << std::endl;
00259     shared_geo_2d.resize(n_shared_geo_2d);        
00260     for (int i = 0;i < n_shared_geo_2d;++ i) {
00261       is >> n_item; 
00262       shared_geo_2d[i].first.resize(n_item);
00263       for (int j = 0;j < n_item;++ j) {
00264         is >> shared_geo_2d[i].first[j].first   
00265            >> shared_geo_2d[i].first[j].second; 
00266         n_rank_shared_2d[shared_geo_2d[i].first[j].first] += 1;
00267       }
00268       is >> shared_geo_2d[i].second; 
00269     }
00270 
00271     for (int i = 0;i < n_rank;++ i) {
00272       proc_os[i] << n_rank_shared_2d[i];
00273       if (n_rank_shared_2d[i] > 0) {
00274         proc_flag[i] = true;
00275       }
00276     }
00277     for (int i = 0;i < n_shared_geo_2d;++ i) {
00278       n_item = shared_geo_2d[i].first.size();
00279       int local_idx;
00280       for (int j = 0;j < n_item;++ j) {
00281         int rnk = shared_geo_2d[i].first[j].first;
00282         int idx = shared_geo_2d[i].first[j].second;
00283         if (rnk == rank) {
00284           local_idx = idx;
00285           break;
00286         }
00287       }
00288       for (int j = 0;j < n_item;++ j) {
00289         int rnk = shared_geo_2d[i].first[j].first;
00290         int idx = shared_geo_2d[i].first[j].second;
00291         if (rnk != rank) {
00292           proc_os[rnk] << idx << geo_2d[local_idx];
00293         }
00294       }
00295     }
00296   }
00297 
00298   if (DIM >= 3) {
00299     std::vector<int> n_rank_shared_3d(n_rank, 0);
00300     int n_shared_geo_3d;
00301     is >> n_shared_geo_3d >> n_shared_geo_3d;
00302     std::cerr << "\t# 3d-shared geometry: " << n_shared_geo_3d << std::endl;
00303     shared_geo_3d.resize(n_shared_geo_3d);        
00304     for (int i = 0;i < n_shared_geo_3d;++ i) {
00305       is >> n_item; 
00306       shared_geo_3d[i].first.resize(n_item);
00307       for (int j = 0;j < n_item;++ j) {
00308         is >> shared_geo_3d[i].first[j].first   
00309            >> shared_geo_3d[i].first[j].second; 
00310         n_rank_shared_3d[shared_geo_3d[i].first[j].first] += 1;
00311       }
00312       is >> shared_geo_3d[i].second; 
00313     }
00314 
00315     for (int i = 0;i < n_rank;++ i) {
00316       proc_os[i] << n_rank_shared_3d[i];
00317       if (n_rank_shared_3d[i] > 0) {
00318         proc_flag[i] = true;
00319       }
00320     }
00321     for (int i = 0;i < n_shared_geo_3d;++ i) {
00322       n_item = shared_geo_3d[i].first.size();
00323       int local_idx;
00324       for (int j = 0;j < n_item;++ j) {
00325         int rnk = shared_geo_3d[i].first[j].first;
00326         int idx = shared_geo_3d[i].first[j].second;
00327         if (rnk == rank) {
00328           local_idx = idx;
00329           break;
00330         }
00331       }
00332       for (int j = 0;j < n_item;++ j) {
00333         int rnk = shared_geo_3d[i].first[j].first;
00334         int idx = shared_geo_3d[i].first[j].second;
00335         if (rnk != rank) {
00336           proc_os[rnk] << idx << geo_3d[local_idx];
00337         }
00338       }
00339     }
00340   }
00341 
00345   std::cerr << "\ttransfering data ..." << std::endl;
00346   int tag = 97, n_req = 0;
00347   proc_flag[rank] = false; 
00348   MPI_Request request[2*n_rank];
00349   MPI_Status status[2*n_rank];
00350   for (int i = 0;i < n_rank;++ i) {
00351     if (proc_flag[i] == false) continue;
00352     MPI_Isend(proc_buf_out[i].start_address(), proc_buf_out[i].size(), MPI_CHAR,
00353               i, tag, _comm, &request[n_req ++]);
00354   }
00355   for (int i = 0;i < n_rank;++ i) {
00356     MPI_Status mpi_status;
00357     if (proc_flag[i] == false) continue;
00358     MPI_Probe(i, tag, _comm, &mpi_status);
00359     int msg_size;
00360     MPI_Get_count(&mpi_status, MPI_CHAR, &msg_size);
00361     proc_buf_in[i].resize(msg_size);
00362     MPI_Irecv(proc_buf_in[i].start_address(), msg_size, MPI_CHAR,
00363               i, tag, _comm, &request[n_req ++]);
00364   }
00365   MPI_Waitall(n_req, request, status);
00366 
00371   std::cerr << "\tdecoding data received ..." << std::endl;
00372   for (int i = 0;i < n_rank;++ i) { 
00373     if (proc_flag[i] == false) continue;
00374     AFEPack::istream<> proc_is(proc_buf_in[i]);
00375     int n_shared_item, local_idx;
00376 
00377     if (DIM >= 0) { 
00378 #define GDIM 0
00379 #define GEO HGeometry<GDIM,this_t::dow>
00380 #define OBJ Shared_object<GEO>
00381       proc_is >> n_shared_item;
00382       for (int j = 0;j < n_shared_item;++ j) {
00383         GEO * remote_obj, * local_obj;
00384         proc_is >> local_idx >> remote_obj;
00385         local_obj = geo_0d[local_idx];
00386         OBJ * p_info = get_shared_info(*local_obj);
00387         if (p_info == NULL) p_info = new_shared_info(*local_obj);
00388         p_info->add_clone(i, remote_obj);
00389       }
00390 #undef OBJ
00391 #undef GEO
00392 #undef GDIM
00393     }
00394 
00395     if (DIM >= 1) {
00396 #define GDIM 1
00397 #define GEO HGeometry<GDIM,this_t::dow>
00398 #define OBJ Shared_object<GEO>
00399       proc_is >> n_shared_item;
00400       for (int j = 0;j < n_shared_item;++ j) {
00401         GEO * remote_obj, * local_obj;
00402         proc_is >> local_idx >> remote_obj;
00403         local_obj = geo_1d[local_idx];
00404         OBJ * p_info = get_shared_info(*local_obj);
00405         if (p_info == NULL) p_info = new_shared_info(*local_obj);
00406         p_info->add_clone(i, remote_obj);
00407       }
00408 #undef OBJ
00409 #undef GEO
00410 #undef GDIM
00411     }
00412 
00413     if (DIM >= 2) {
00414 #define GDIM 2
00415 #define GEO HGeometry<GDIM,this_t::dow>
00416 #define OBJ Shared_object<GEO>
00417       proc_is >> n_shared_item;
00418       for (int j = 0;j < n_shared_item;++ j) {
00419         GEO * remote_obj, * local_obj;
00420         proc_is >> local_idx >> remote_obj;
00421         local_obj = geo_2d[local_idx];
00422         OBJ * p_info = get_shared_info(*local_obj);
00423         if (p_info == NULL) p_info = new_shared_info(*local_obj);
00424         p_info->add_clone(i, remote_obj);
00425       }
00426 #undef OBJ
00427 #undef GEO
00428 #undef GDIM
00429     }
00430 
00431     if (DIM >= 3) {
00432 #define GDIM 3
00433 #define GEO HGeometry<GDIM,this_t::dow>
00434 #define OBJ Shared_object<GEO>
00435       proc_is >> n_shared_item;
00436       for (int j = 0;j < n_shared_item;++ j) {
00437         GEO * remote_obj, * local_obj;
00438         proc_is >> local_idx >> remote_obj;
00439         local_obj = geo_3d[local_idx];
00440         OBJ * p_info = get_shared_info(*local_obj);
00441         if (p_info == NULL) p_info = new_shared_info(*local_obj);
00442         p_info->add_clone(i, remote_obj);
00443       }
00444 #undef OBJ
00445 #undef GEO
00446 #undef GDIM
00447     }
00448   }
00449 
00450   set_shared_info_sent();
00451 }
00452 
00453 TEMPLATE
00454 void THIS::renumerateRootElement(void (*f)(double, double, double,
00455                                            double&,double&,double&))
00456 {
00457   std::cerr << "Renumerating element of the mesh using Hibert space filling curve ..." 
00458             << std::flush;
00459   int n_ele = this->n_rootElement();
00460   std::vector<HGeometry<DIM,this_t::dow>*> tmp_ele(n_ele);
00461   std::vector<double> x(n_ele, 0.0), y(n_ele, 0.0), z(n_ele, 0.0);
00462   typename std::list<HGeometry<DIM,this_t::dow>*>::iterator
00463     the_ele = this->rootElement().begin();
00464   for (int i = 0;i < n_ele;++ i, ++ the_ele) {
00465     tmp_ele[i] = *the_ele;
00466     HGeometry<DIM,this_t::dow>& ele = **the_ele;
00467     const int& n_vtx = ele.n_vertex;
00468     for (int j = 0;j < n_vtx;++ j) {
00469       afepack::Point<this_t::dow>& pnt = *ele.vertex[j];
00470       x[i] += pnt[0]; y[i] += pnt[1];
00471       if (this_t::dow == 3) z[i] += pnt[2];
00472     }
00473     x[i] /= n_vtx; y[i] /= n_vtx;
00474     if (this_t::dow == 3) z[i] /= n_vtx;
00475   }
00476   std::vector<int> new_index(n_ele);
00477   hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0], f);
00478 
00479   the_ele = this->rootElement().begin();
00480   for (int i = 0;i < n_ele;++ i, ++ the_ele) {
00481     *the_ele = tmp_ele[new_index[i]];
00482   }
00483   std::cerr << " OK!" << std::endl;
00484 }
00485 
00486 
00487 TEMPLATE
00488 void THIS::eraseRootElement(u_int level) {
00489   for (u_int i = 0;i < level;++ i) {
00490     this->eraseRootElementOneLevel();
00491   }
00492 }
00493 
00494 TEMPLATE
00495 void THIS::eraseRootElementOneLevel() {
00497   std::list<HGeometry<DIM,this_t::dow>*> new_root;
00498   typename std::list<HGeometry<DIM,this_t::dow>*>::iterator
00499     the_ele = this->rootElement().begin(),
00500     end_ele = this->rootElement().end();
00501   for (;the_ele != end_ele;++ the_ele) {
00502     for (u_int i = 0;i < HGeometry<DIM,this_t::dow>::n_child;++ i) {
00503       HGeometry<DIM,this_t::dow> * p_geo = (*the_ele)->child[i];
00504       if (p_geo == NULL) {
00505         std::cerr << "Macro element is not refined. The root elements are not erased!" 
00506                   << std::endl;
00507         abort();
00508       }
00509       new_root.push_back(p_geo);
00510       this->nullParent(*p_geo);
00511     }
00512   }
00513   new_root.swap(this->rootElement());
00514 
00521   property_id_t<bool> pid;
00522   new_property_id(pid);
00523   end_ele = new_root.end();
00524   the_ele = new_root.begin();
00525   for (;the_ele != end_ele;++ the_ele) {
00526     this->tryDeleteGeometry(*the_ele, pid);
00527   }
00528   free_property_id(pid);
00529   the_ele = new_root.begin();
00530   for (;the_ele != end_ele;++ the_ele) {
00531     this->deleteGeometry(*the_ele);
00532   }
00533 }
00534 
00535 TEMPLATE
00536 template <class HGEO>
00537 void THIS::set_is_used_up(HGEO& geo) const {
00538   if (geo.parent != NULL) {
00539     set_is_used_up(*geo.parent);
00540   }
00541   u_int n_bound = geo.n_boundary;
00542   for (u_int i = 0;i < n_bound;++ i) {
00543     this->set_is_used_up(*geo.boundary[i]);
00544   }
00545   HTools tools;
00546   tools.setGeometryUsed(geo);
00547 }
00548 
00549 TEMPLATE
00550 void THIS::set_is_used_up() {
00555   typename std::list<HGeometry<this_t::dim,this_t::dow>*>::iterator
00556     the_ele = this->rootElement().begin(),
00557     end_ele = this->rootElement().end();
00558   for (;the_ele != end_ele;++ the_ele) {
00559     if ((*the_ele)->parent != NULL) {
00560       this->set_is_used_up(*(*the_ele)->parent);
00561     }
00562   }
00563 
00564 }
00565 
00566 #undef THIS
00567 #undef TEMPLATE
00568 
00570 
00571 #define TEMPLATE template <class FOREST>
00572 #define THIS HGeometryMatcher<FOREST>
00573 
00574 TEMPLATE
00575 bool THIS::match_geometry() {
00576   typedef THIS this_t;
00577 
00578   std::cerr << "Rank " << _forest->rank()
00579             << ": Matching h-geometry among partitions ..." << std::endl;
00580   bool is_operated = false;
00581   do {
00582     _is_operated = false;
00583 
00584     if (this_t::dim >= 3)
00585       sync_data(_forest->communicator(), _forest->_shared_list_3d, *this, 
00586                 &this_t::template pack_match_geometry<3>,
00587                 &this_t::template unpack_match_geometry<3>,
00588                 &this_t::template is_pack_match_geometry<3>);
00589     if (this_t::dim >= 2)
00590       sync_data(_forest->communicator(), _forest->_shared_list_2d, *this,
00591                 &this_t::template pack_match_geometry<2>,
00592                 &this_t::template unpack_match_geometry<2>,
00593                 &this_t::template is_pack_match_geometry<2>);
00594     if (this_t::dim >= 1)
00595       sync_data(_forest->communicator(), _forest->_shared_list_1d, *this,
00596                 &this_t::template pack_match_geometry<1>,
00597                 &this_t::template unpack_match_geometry<1>,
00598                 &this_t::template is_pack_match_geometry<1>);
00599 
00600     int src = _is_operated?1:0, dst = 0; 
00601     MPI_Allreduce(&src, &dst, 1, MPI_INT, MPI_SUM, _forest->communicator());
00602     if ((_is_operated = (dst > 0))) {
00603       is_operated = true;
00604     }
00605   } while (_is_operated);
00606   _forest->set_shared_info_sent();
00607 
00608   return is_operated;
00609 }
00610 
00611 TEMPLATE
00612 template <int GDIM> bool 
00613 THIS::is_pack_match_geometry(HGeometry<GDIM,dow> * geo) const {
00614   u_int n_unsent_vtx = 0;
00615   u_int n_unsent_bnd = 0;
00616   if (GDIM == 1) { 
00620 
00621     u_int n_vertex = geo->n_vertex;
00622     for (u_int i = 0;i < n_vertex;++ i) {
00623       if (! is_shared_info_sent(*geo->vertex[i])) {
00624         n_unsent_vtx += 1;
00625       }
00626     }
00627   } else { 
00630 
00631     u_int n_boundary = geo->n_boundary;
00632     for (u_int i = 0;i < n_boundary;++ i) {
00633       HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
00634       if (bnd->parent == NULL && 
00635           (! is_shared_info_sent(*bnd))) {
00636         n_unsent_bnd += 1;
00637       }
00638     }
00639   }
00640 
00641   u_int n_child = 0;
00642   if (geo->isRefined() && 
00643       (! is_shared_info_sent(*geo->child[0]))) {
00644     n_child = geo->n_child;
00645   }
00646   return ((n_unsent_vtx > 0) ||
00647           (n_unsent_bnd > 0) ||
00648           (n_child > 0));
00649 }
00650 
00651 TEMPLATE
00652 template <int GDIM> void 
00653 THIS::pack_match_geometry(HGeometry<GDIM,dow> * geo,
00654                           int remote_rank,
00655                           AFEPack::ostream<>& os) {
00656   afepack::Point<dow> grp; 
00657 
00658   if (GDIM == 1) { 
00662 
00663     u_int n_vertex = geo->n_vertex;
00664     u_int n_unsent_vtx = 0;
00665     for (u_int i = 0;i < n_vertex;++ i) {
00666       if (! is_shared_info_sent(*geo->vertex[i])) {
00667         n_unsent_vtx += 1;
00668       }
00669     }
00670 
00671     os << n_unsent_vtx; 
00672     for (u_int i = 0;i < n_vertex;++ i) {
00673       HGeometry<0,dow> * vtx = geo->vertex[i];
00674       if (! is_shared_info_sent(*vtx)) {
00675         grp = (*vtx); 
00676         os << vtx << grp;
00677       }
00678     }
00679   } else { 
00682 
00683     u_int n_boundary = geo->n_boundary;
00684     u_int n_unsent_bnd = 0;
00685     for (u_int i = 0;i < n_boundary;++ i) {
00686       HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
00687       if (bnd->parent == NULL && 
00688           (! is_shared_info_sent(*bnd))) {
00689         n_unsent_bnd += 1;
00690       }
00691     }
00692 
00693     os << n_unsent_bnd; 
00694     for (u_int i = 0;i < n_boundary;++ i) {
00695       HGeometry<GDIM-1,dow> * bnd = geo->boundary[i];
00696       if (bnd->parent == NULL && 
00697           (! is_shared_info_sent(*bnd))) {
00698         geometry_reference_point(*bnd, grp);
00699         os << bnd << grp;
00700       }
00701     }
00702   }
00703 
00707   if (geo->isRefined() && 
00708       (! is_shared_info_sent(*geo->child[0]))) {
00715     u_int n_child = geo->n_child;
00716     os << n_child; 
00717     for (u_int i = 0;i < n_child;++ i) {
00718       HGeometry<GDIM,dow> * chd = geo->child[i]; 
00719       geometry_reference_point(*chd, grp); 
00720       os << chd << grp; 
00721     }
00722   } else { 
00723     u_int n_child = 0;
00724     os << n_child;
00725   }
00726 }
00727 
00728 TEMPLATE
00729 template <int GDIM> void 
00730 THIS::unpack_match_geometry(HGeometry<GDIM,dow> * geo,
00731                             int remote_rank,
00732                             AFEPack::istream<>& is) {
00734   double h = distance(*geo->vertex[0], *geo->vertex[1]);
00735 
00736   if (GDIM == 1) { 
00740 
00741     u_int n_unsent_vtx = 0;
00742     is >> n_unsent_vtx;
00743     std::vector<HGeometry<0,dow>*> remote_vtx(n_unsent_vtx);
00744     std::vector<afepack::Point<dow> > pnt(n_unsent_vtx);
00745     for (u_int i = 0;i < n_unsent_vtx;++ i) {
00746       is >> remote_vtx[i] >> pnt[i];
00747     }
00748 
00750     u_int n_vertex = geo->n_vertex;
00751     for (u_int i = 0;i < n_unsent_vtx;++ i) {
00752       for (u_int j = 0;j < n_vertex;++ j) {
00753         HGeometry<0,dow> * vtx = geo->vertex[j];
00754         int type = matcher().value(pnt[i], *vtx, 1.0e-4*h);
00755         if (type >= 0) {
00756           Shared_object<HGeometry<0,dow> > * p_info = get_shared_info(*vtx);
00757           if (p_info == NULL) p_info = new_shared_info(*vtx);
00758           _is_operated |= p_info->add_clone(remote_rank, type, remote_vtx[i]);
00759           break;
00760         }
00761       }
00762     }
00763   } else { 
00767 
00768     u_int n_unsent_bnd = 0;
00769     is >> n_unsent_bnd;
00770     std::vector<HGeometry<GDIM-1,dow>*> remote_bnd(n_unsent_bnd);
00771     std::vector<afepack::Point<dow> > remote_grp(n_unsent_bnd);
00772     for (u_int i = 0;i < n_unsent_bnd;++ i) {
00773       is >> remote_bnd[i] >> remote_grp[i];
00774     }
00775 
00777     afepack::Point<dow> grp;
00778     u_int n_boundary = geo->n_boundary;
00779     for (u_int i = 0;i < n_unsent_bnd;++ i) {
00780       for (u_int j = 0;j < n_boundary;++ j) {
00781         HGeometry<GDIM-1,dow> * bnd = geo->boundary[j];
00782         if (bnd->parent != NULL) continue;
00783         geometry_reference_point(*bnd, grp);
00784         int type = matcher().value(remote_grp[i], grp, 1.0e-4*h);
00785         if (type >= 0) {
00786           Shared_object<HGeometry<GDIM-1,dow> > * p_info = get_shared_info(*bnd);
00787           if (p_info == NULL) p_info = new_shared_info(*bnd);
00788           _is_operated |= p_info->add_clone(remote_rank, type, remote_bnd[i]);
00789           break;
00790         }
00791       }
00792     }
00793   }
00794       
00795   u_int n_child;
00796   is >> n_child;
00797   if (n_child != 0) { 
00798     assert (n_child == geo->n_child); 
00799         
00801     std::vector<HGeometry<GDIM,dow>*> remote_chd(n_child);
00802     std::vector<afepack::Point<dow> > remote_grp(n_child);
00803     for (u_int i = 0;i < n_child;++ i) {
00804       is >> remote_chd[i] >> remote_grp[i];
00805     }
00806 
00811     if (! (_forest->is_dummy(*geo) || geo->isRefined())) {
00812       geo->refine();
00813     }
00814 
00818     if (geo->isRefined()) { 
00819       afepack::Point<dow> grp;
00820       u_int n_matched_child = 0;
00821       for (u_int i = 0;i < n_child;++ i) {
00822         HGeometry<GDIM,dow> * chd = geo->child[i];
00823         geometry_reference_point(*chd, grp);
00824         for (u_int j = 0;j < n_child;++ j) {
00825           if (remote_chd[j] == NULL) continue; 
00826 
00828           int type = matcher().value(grp, remote_grp[j], 1.0e-4*h);
00829           if (type >= 0) {
00830             Shared_object<HGeometry<GDIM,dow> > * p_info = get_shared_info(*chd); 
00831             if (p_info == NULL) p_info = new_shared_info(*chd);
00832             _is_operated |= p_info->add_clone(remote_rank, type, remote_chd[j]);
00833             remote_chd[j] = NULL;
00834             n_matched_child += 1;
00835             break;
00836           }
00837         }
00838       }
00839       assert (n_matched_child == n_child);
00840     }
00841   }
00842 }
00843 
00844 TEMPLATE
00845 bool THIS::sync_is_used() {
00846   typedef THIS this_t;
00847   std::cerr << "Rank " << _forest->rank()
00848             << ": Sync is_used flag ..." << std::endl;
00849 
00850   _is_operated = false;
00851   if (this_t::dim >= 1)
00852     sync_data(_forest->communicator(), _forest->_shared_list_1d, *this, 
00853               &this_t::template pack_sync_is_used<1>,
00854               &this_t::template unpack_sync_is_used<1>);
00855   if (this_t::dim >= 2)
00856     sync_data(_forest->communicator(), _forest->_shared_list_2d, *this, 
00857               &this_t::template pack_sync_is_used<2>,
00858               &this_t::template unpack_sync_is_used<2>);
00859   if (this_t::dim >= 3)
00860     sync_data(_forest->communicator(), _forest->_shared_list_3d, *this, 
00861               &this_t::template pack_sync_is_used<3>,
00862               &this_t::template unpack_sync_is_used<3>);
00863 
00864   int src = _is_operated?1:0, dst = 0; 
00865   MPI_Allreduce(&src, &dst, 1, MPI_INT, MPI_SUM, _forest->communicator());
00866   return (dst > 0);
00867 }
00868 
00869 TEMPLATE
00870 template <int GDIM> void 
00871 THIS::pack_sync_is_used(HGeometry<GDIM,dow> * geo,
00872                         int remote_rank,
00873                         AFEPack::ostream<>& os) {
00874   bool is_used = _tools.isGeometryUsed(*geo);
00875   os << is_used;
00876 }
00877 
00878 TEMPLATE
00879 template <int GDIM> void 
00880 THIS::unpack_sync_is_used(HGeometry<GDIM,dow> * geo,
00881                           int remote_rank,
00882                           AFEPack::istream<>& is) {
00883   bool is_used;
00884   is >> is_used;
00885 
00886   if (is_used) {
00887     if (! _tools.isGeometryUsed(*geo)) {
00888       _tools.setGeometryUsed(*geo);
00889       _is_operated = true;
00890     }
00891   }
00892 }
00893 
00894 #undef THIS
00895 
00897 
00898 #define THIS BirdView<FOREST>
00899 
00900 TEMPLATE void 
00901 THIS::semiregularize() {
00902   forest_t& forest = this->getForest();
00903 
00904   if (! forest.lock()) { // 对几何遗传树进行加锁
00905     std::cerr << "The hierarchy geometry tree is locked, aborting ...";
00906     abort();
00907   }
00908 
00909   if (forest.rank() == 0)
00910     std::cerr << "Semiregularizing the mesh ...  " << std::flush;
00911   const char * timer = "-/|\\";
00912   int n_element_refined = 0;
00913   this->prepareSemiregularize();
00914 
00915   int round = 0;
00916   bool is_operated;
00917   HGeometryMatcher<forest_t> matcher(forest);
00921   matcher.match_geometry();
00922   matcher.sync_is_used(); 
00923   do {
00924     bool flag;
00925     do {
00926       flag = false;
00927       this->semiregularizeHelper(flag, n_element_refined);
00928     } while (flag);
00929 
00930     if (forest.rank() == 0)
00931       std::cerr << "\b" << timer[round] << std::flush;
00932     round = (round + 1)%4;
00933 
00934     is_operated = false;
00935     if (1) { // forest.n_rank() > 1) {
00936       is_operated |= matcher.match_geometry();
00937       is_operated |= matcher.sync_is_used();
00938     }
00939   } while (is_operated);
00940 
00941 
00942   int n_ele_refined = 0;
00943   MPI_Reduce(&n_element_refined, &n_ele_refined, 1, MPI_INT, MPI_SUM, 0,
00944              forest.communicator());
00945   if (forest.rank() == 0) 
00946     std::cerr << "\bOK!\n"
00947               << "\t" << n_ele_refined 
00948               << " elements refined in semiregularization."
00949               << std::endl;
00950 }
00951 
00952 TEMPLATE
00953 void THIS::eraseRootElement(u_int level) {
00954   for (u_int i = 0;i < level;++ i) {
00955     this->eraseRootElementOneLevel();
00956   }
00957 }
00958 
00959 TEMPLATE
00960 void THIS::eraseRootElementOneLevel() {
00962   typedef typename base_t::container_t container_t;
00963   typedef typename base_t::element_t element_t;
00964   container_t new_root;
00965   typename container_t::iterator
00966     the_ele = this->rootElement().begin(),
00967     end_ele = this->rootElement().end();
00968   for (;the_ele != end_ele;++ the_ele) {
00969     for (u_int i = 0;i < element_t::n_child;++ i) {
00970       element_t * p_geo = (*the_ele)->child[i];
00971       if (p_geo == NULL) {
00972         std::cerr << "Macro element is not refined. The root elements CAN'T be erased!" 
00973                   << std::endl;
00974         abort();
00975       }
00976       new_root.push_back(p_geo);
00977       p_geo->parent = NULL;
00978     }
00979   }
00980   new_root.swap(this->rootElement());
00981 
00982   the_ele = new_root.begin();
00983   end_ele = new_root.end();
00984   for (;the_ele != end_ele;++ the_ele) {
00985     delete *the_ele;
00986   }
00987 }
00988 
00989 #undef THIS
00990 #undef TEMPLATE
00991