AFEPack
|
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