AFEPack
|
00001 00012 #ifndef _HGeometry_templates_h_ 00013 #define _HGeometry_templates_h_ 00014 00015 #include "HGeometry.h" 00016 00017 #define TEMPLATE template <int DIM, int DOW> 00018 #define THIS HGeometry<DIM,DOW> 00019 00020 TEMPLATE 00021 THIS::HGeometry() : HGeometryBase() 00022 { 00023 index = 0; 00024 for (int i = 0;i < THIS::n_vertex;i ++) 00025 vertex[i] = NULL; 00026 for (int i = 0;i < THIS::n_boundary;i ++) 00027 boundary[i] = NULL; 00028 parent = NULL; 00029 for (int i = 0;i < THIS::n_child;i ++) 00030 child[i] = NULL; 00031 bmark = 0; 00032 } 00033 00034 TEMPLATE 00035 bool THIS::isRefined() const 00036 { 00037 return (child[0] != NULL); 00038 } 00039 00040 TEMPLATE 00041 void THIS::refine() 00042 {} 00043 00044 TEMPLATE 00045 void THIS::checkIntegrity() const 00046 { 00051 int k; 00052 for (int i = 0;i < THIS::n_boundary;i ++) { 00053 const bound_t * b = boundary[i]; 00054 for (int j = 0;j < bound_t::n_vertex;j ++) { 00055 for (k = 0;k < THIS::n_vertex;k ++) { 00056 if (b->vertex[j] == vertex[k]) { 00057 break; 00058 } 00059 } 00060 Assert(k < THIS::n_vertex, ExcInternalError()); 00061 } 00062 } 00068 if (isRefined()) { 00069 for (int i = 0;i < THIS::n_child;i ++) { 00070 child[i]->checkIntegrity(); 00071 } 00072 } 00073 } 00074 00075 TEMPLATE 00076 bool THIS::isIncludePoint(const afepack::Point<DOW>& p) const 00077 { 00078 return true; 00079 } 00080 00081 #undef THIS /// HGeometry<DIM,DOW> 00082 00084 00085 #define THIS HGeometryTree<DIM,DOW> 00086 00087 TEMPLATE 00088 void THIS::checkIntegrity() 00089 { 00090 RootIterator the_ele = beginRootElement(); 00091 RootIterator end_ele = endRootElement(); 00092 for (;the_ele != end_ele;++ the_ele) { 00093 the_ele->checkIntegrity(); 00094 } 00095 } 00096 00097 TEMPLATE 00098 void THIS::clear() 00099 { 00100 Tools tools; 00101 RootIterator 00102 the_ele = beginRootElement(), 00103 end_ele = endRootElement(); 00104 for (;the_ele != end_ele;++ the_ele) { 00105 tools.clearIndex(*the_ele); 00106 } 00107 00108 the_ele = beginRootElement(); 00109 for (;the_ele != end_ele;++ the_ele) { 00110 tools.incrIndex(*the_ele); 00111 } 00112 00113 the_ele = beginRootElement(); 00114 for (;the_ele != end_ele;++ the_ele) { 00115 tools.decrIndex(*the_ele); 00116 } 00117 00118 root_element.clear(); 00119 } 00120 00121 TEMPLATE 00122 void THIS::readEasyMesh(const std::string& filename) 00123 { 00124 std::cerr << "THIS::readEasyMesh is only avaiable for 2-dimensional case." 00125 << std::endl; 00126 abort(); 00127 } 00128 00129 TEMPLATE 00130 void THIS::readMesh(const std::string& filename) 00131 { 00132 std::cerr << "Reading in mesh data file " 00133 << filename 00134 << " as geometry tree root ..." 00135 << std::endl; 00136 std::ifstream is(filename.c_str()); 00137 00139 u_int n_point; 00140 is >> n_point; 00141 std::cerr << "\t# points: " << n_point << std::endl; 00142 std::vector<afepack::Point<DOW> > point(n_point); 00143 for (u_int i = 0;i < n_point;i ++) is >> point[i]; 00144 is >> n_point; 00145 std::vector<HGeometry<0,DOW> *> geo_0d(n_point, (HGeometry<0,DOW> *)NULL); 00146 for (u_int i = 0;i < n_point;i ++) { 00147 u_int j, k; 00148 is >> j; geo_0d[j] = new HGeometry<0,DOW>(); 00149 is >> k >> k; *dynamic_cast<afepack::Point<DOW> *>(geo_0d[j]) = point[k]; 00150 is >> k >> k >> geo_0d[j]->bmark; 00151 } 00152 point.clear(); 00153 00154 #define GDIM 1 00155 u_int n_geo_1d; 00156 std::vector<HGeometry<GDIM,DOW>*> geo_1d; 00157 if (DIM >= 1) { 00158 is >> n_geo_1d; 00159 std::cerr << "\t# 1D-geometry: " << n_geo_1d << std::endl; 00160 geo_1d.resize(n_geo_1d, NULL); 00161 for (u_int i = 0;i < n_geo_1d;i ++) { 00162 u_int j, k, l; 00163 is >> j >> k; geo_1d[j] = new HGeometry<GDIM,DOW>(); 00164 for (k = 0;k < GDIM + 1;k ++) { 00165 is >> l; geo_1d[j]->vertex[k] = geo_0d[l]; 00166 } 00167 is >> k; 00168 for (k = 0;k < GDIM + 1;k ++) { 00169 is >> l; 00170 } 00171 is >> geo_1d[j]->bmark; 00172 } 00173 } 00174 #undef GDIM 00175 00176 #define GDIM 2 00177 u_int n_geo_2d; 00178 std::vector<HGeometry<GDIM,DOW>*> geo_2d; 00179 if (DIM >= 2) { 00180 is >> n_geo_2d; 00181 std::cerr << "\t# 2D-geometry: " << n_geo_2d << std::endl; 00182 geo_2d.resize(n_geo_2d, NULL); 00183 for (u_int i = 0;i < n_geo_2d;i ++) { 00184 u_int j, k, l; 00185 is >> j >> k; geo_2d[j] = new HGeometry<GDIM,DOW>(); 00186 for (k = 0;k < GDIM + 1;k ++) { 00187 is >> l; geo_2d[j]->vertex[k] = geo_0d[l]; 00188 } 00189 is >> k; 00190 for (k = 0;k < GDIM + 1;k ++) { 00191 is >> l; geo_2d[j]->boundary[k] = geo_1d[l]; 00192 } 00193 is >> geo_2d[j]->bmark; 00194 } 00195 } 00196 #undef GDIM 00197 00198 #define GDIM 3 00199 u_int n_geo_3d; 00200 std::vector<HGeometry<GDIM,DOW>*> geo_3d; 00201 if (DIM >= 3) { 00202 is >> n_geo_3d; 00203 std::cerr << "\t# 3D-geometry: " << n_geo_3d << std::endl; 00204 geo_3d.resize(n_geo_3d, NULL); 00205 for (u_int i = 0;i < n_geo_3d;i ++) { 00206 u_int j, k, l; 00207 is >> j >> k; geo_3d[j] = new HGeometry<GDIM,DOW>(); 00208 for (k = 0;k < GDIM + 1;k ++) { 00209 is >> l; geo_3d[j]->vertex[k] = geo_0d[l]; 00210 } 00211 is >> k; 00212 for (k = 0;k < GDIM + 1;k ++) { 00213 is >> l; geo_3d[j]->boundary[k] = geo_2d[l]; 00214 } 00215 is >> geo_3d[j]->bmark; 00216 } 00217 #undef GDIM 00218 } 00219 is.close(); 00220 00221 if (DIM == 1) { 00222 for (u_int i = 0;i < n_geo_1d;++ i) { 00223 this->rootElement().push_back((HGeometry<DIM,DOW> *)geo_1d[i]); 00224 } 00225 } else if (DIM == 2) { 00226 for (u_int i = 0;i < n_geo_2d;++ i) { 00227 this->rootElement().push_back((HGeometry<DIM,DOW> *)geo_2d[i]); 00228 } 00229 } else if (DIM == 3) { 00230 for (u_int i = 0;i < n_geo_3d;++ i) { 00231 this->rootElement().push_back((HGeometry<DIM,DOW> *)geo_3d[i]); 00232 } 00233 } 00234 } 00235 00236 #undef THIS /// HGeometryTree<DIM,DOW> 00237 00239 00240 #define THIS HElement<DIM,DOW> 00241 00242 TEMPLATE 00243 THIS::HElement() 00244 : value(-1), parent(NULL), child(THIS::n_child, NULL) {} 00245 00246 TEMPLATE 00247 THIS::HElement(const element_t& ele) : 00248 index(ele.index), value(ele.value), 00249 indicator(ele.indicator), 00250 h_element(ele.h_element) {} 00251 00252 TEMPLATE 00253 THIS::~HElement() 00254 {} 00255 00256 TEMPLATE 00257 typename THIS::element_t& THIS::operator=(const element_t& ele) { 00258 index = ele.index; 00259 value = ele.value; 00260 indicator = ele.indicator; 00261 h_element = ele.h_element; 00262 return *this; 00263 } 00264 00265 TEMPLATE 00266 void THIS::refine() 00267 {} 00268 00269 TEMPLATE 00270 bool THIS::isRefined() const 00271 { 00272 return (child.size() > 0 && child[0] != NULL); 00273 } 00274 00275 TEMPLATE 00276 void THIS::checkIntegrity() const 00277 { 00278 if (!isRefined()) return; 00279 for (int i = 0;i < THIS::n_child;i ++) { 00280 child[i]->checkIntegrity(); 00281 } 00282 } 00283 00284 TEMPLATE 00285 bool THIS::isIncludePoint(const afepack::Point<DOW>& p) const 00286 { 00287 return h_element->isIncludePoint(p); 00288 } 00289 00290 #undef THIS /// HElement<DIM,DOW> 00291 00293 00294 #define THIS IrregularMesh<DIM,DOW> 00295 00296 TEMPLATE 00297 void THIS::semiregularize() 00298 { 00299 if (! geometry_tree->lock()) { // 对几何遗传树进行加锁 00300 std::cerr << "The hierarchy geometry tree is locked, aborting ..."; 00301 abort(); 00302 } 00303 00304 std::cerr << "Semiregularizing the mesh ... " << std::flush; 00305 int round = 0; 00306 const char * timer = "-/|\\"; 00307 bool flag; 00308 int n_element_refined = 0; 00309 prepareSemiregularize(); 00310 do { 00311 std::cerr << "\b" << timer[round] << std::flush; 00312 round = (round + 1)%4; 00313 flag = false; 00314 semiregularizeHelper(flag, n_element_refined); 00315 } while (flag); 00316 std::cerr << "\bOK!\n" 00317 << "\t" << n_element_refined 00318 << " elements refined in semiregularization." 00319 << std::endl; 00320 } 00321 00322 TEMPLATE 00323 void THIS::semiregularizeHelper(bool& flag, 00324 int& n_element_refined) 00325 { 00326 RootIterator the_ele = beginRootElement(); 00327 RootIterator end_ele = endRootElement(); 00328 for (;the_ele != end_ele;++ the_ele) { 00329 semiregularizeHelper(flag, *the_ele, n_element_refined); 00330 } 00331 } 00332 00333 TEMPLATE 00334 void THIS::semiregularizeHelper(bool& flag, 00335 element_t& element, 00336 int& n_element_refined) 00337 { 00338 if (element.value == 0) { // if this element is a leaf element 00339 h_element_t& h_element = *(element.h_element); 00340 flag |= Tools().semiregularizeBoundary(h_element); 00341 if (! Tools().isSemiregular(h_element)) { // if this element is not semiregular 00342 flag = true; 00343 element.refine(); 00344 element.value = 1; 00345 for (int i = 0;i < element_t::n_child;i ++) { 00346 element.child[i]->value = 0; 00347 Tools().setGeometryUsed(*(h_element.child[i])); 00348 } 00349 n_element_refined ++; 00350 } 00351 } else { // if its not a leaf element 00352 assert (element.value == 1); 00353 for (int i = 0;i < element_t::n_child;i ++) { // then check its children 00354 semiregularizeHelper(flag, *element.child[i], n_element_refined); 00355 } 00356 } 00357 } 00358 00359 TEMPLATE 00360 void THIS::prepareSemiregularize() 00361 { 00362 Tools tools; 00363 RootIterator the_root_element = beginRootElement(); 00364 RootIterator end_root_element = endRootElement(); 00365 for (;the_root_element != end_root_element;++ the_root_element) { 00366 tools.setGeometryUnusedRecursively(*(the_root_element->h_element)); 00367 } 00368 RootFirstIterator 00369 the_ele = beginRootFirstElement(), 00370 end_ele = endRootFirstElement(); 00371 for (;the_ele != end_ele;++ the_ele) { 00372 tools.setGeometryUsed(*(the_ele->h_element)); 00373 } 00374 } 00375 00376 TEMPLATE 00377 THIS::IrregularMesh(tree_t& h_geometry_tree) 00378 { 00379 setGeometryTree(&h_geometry_tree); 00380 regular_mesh = NULL; 00381 } 00382 00383 TEMPLATE void 00384 THIS::reinit(tree_t& h_geometry_tree, 00385 bool is_bare) 00386 { 00387 if (is_bare) { 00388 geometry_tree = &h_geometry_tree; 00389 } else { 00390 clear(); 00391 setGeometryTree(&h_geometry_tree); 00392 } 00393 } 00394 00395 TEMPLATE 00396 void THIS::globalRefine(unsigned int i) 00397 { 00398 MeshAdaptor<DIM, DOW> mesh_adaptor(*this); 00399 mesh_adaptor.globalRefine(i); 00400 } 00401 00402 TEMPLATE 00403 void THIS::randomRefine(double percent) 00404 { 00405 MeshAdaptor<DIM, DOW> mesh_adaptor(*this); 00406 mesh_adaptor.randomRefine(percent); 00407 } 00408 00409 TEMPLATE 00410 void THIS::setGeometryTree(tree_t * h_geometry_tree) 00411 { 00412 std::cerr << "Constructing the root mesh from hierarchy geometry tree ..." << std::endl; 00413 int i, j; 00414 geometry_tree = h_geometry_tree; 00415 00416 // construct the root elements at first 00417 std::cerr << "\tconstructing elements ..." << std::flush; 00418 unsigned int n_root_element = geometry_tree->n_rootElement(); 00419 std::vector<element_t *> h_element(n_root_element); 00420 typename tree_t::RootIterator 00421 the_h_element = geometry_tree->beginRootElement(), 00422 end_h_element = geometry_tree->endRootElement(); 00423 for (i = 0;the_h_element != end_h_element;++ the_h_element) 00424 the_h_element->index = i ++; 00425 the_h_element = geometry_tree->beginRootElement(); 00426 for (i = 0;the_h_element != end_h_element;++ the_h_element, ++ i) { 00427 element_t * element = new element_t(); 00428 element->value = 0; 00429 element->h_element = &*the_h_element; 00430 root_element.push_back(element); 00431 h_element[i] = element; 00432 } 00433 std::cerr << " OK!" << std::endl; 00434 } 00435 00436 TEMPLATE 00437 THIS::~IrregularMesh() 00438 { 00439 clear(); 00440 } 00441 00442 TEMPLATE 00443 void THIS::clear() 00444 { 00445 if (geometry_tree != NULL) 00446 geometry_tree = NULL; 00447 00448 RootIterator the_ele = beginRootElement(); 00449 RootIterator end_ele = endRootElement(); 00450 for (;the_ele != end_ele;++ the_ele) { 00451 this->deleteTree(&*the_ele); 00452 } 00453 root_element.clear(); 00454 00455 if (regular_mesh != NULL) { 00456 delete regular_mesh; 00457 regular_mesh = NULL; 00458 } 00459 } 00460 00461 TEMPLATE 00462 void THIS::checkIntegrity() 00463 { 00464 RootIterator the_ele = beginRootElement(); 00465 RootIterator end_ele = endRootElement(); 00466 for (;the_ele != end_ele;++ the_ele) { 00467 the_ele->checkIntegrity(); 00468 } 00469 } 00470 00471 TEMPLATE 00472 void THIS::deleteTree(element_t * element) 00473 { 00474 if (element->isRefined()) { 00475 for (int i = 0;i < element_t::n_child;i ++) { 00476 this->deleteTree(element->child[i]); 00477 } 00478 } 00479 delete element; 00480 } 00481 00482 TEMPLATE 00483 void THIS::copyTree(const element_t * src, 00484 element_t * dst) 00485 { 00486 dst->index = src->index; 00487 dst->value = src->value; 00488 dst->indicator = src->indicator; 00489 dst->h_element = src->h_element; 00490 if (src->isRefined()) { 00491 dst->refine(); 00492 for (int i = 0;i < element_t::n_child;i ++) { 00493 this->copyTree(src->child[i], dst->child[i]); 00494 } 00495 } 00496 } 00497 00498 TEMPLATE 00499 void THIS::copyNonnegtiveSubtree(const element_t * src, 00500 element_t * dst) 00501 { 00502 assert (src->value == 0 || src->value == 1); 00503 dst->value = src->value; 00504 dst->index = src->index; 00505 dst->h_element = src->h_element; 00506 if (src->value == 1) { 00507 dst->refine(); 00508 for (int i = 0;i < element_t::n_child;i ++) { 00509 this->copyNonnegtiveSubtree(src->child[i], dst->child[i]); 00510 } 00511 } 00512 } 00513 00514 TEMPLATE 00515 void THIS::writeFormatted(const std::string& filename) 00516 { 00517 std::ofstream os(filename.c_str()); 00518 os << (*this); 00519 os.close(); 00520 } 00521 00522 00523 00524 TEMPLATE 00525 THIS::IrregularMesh() : 00526 geometry_tree(NULL), regular_mesh(NULL) 00527 {} 00528 00529 00530 TEMPLATE 00531 THIS::IrregularMesh(const ir_mesh_t& mesh) 00532 { 00533 if (mesh.geometry_tree != NULL) { 00534 setGeometryTree(mesh.geometry_tree); 00535 copyNonnegtiveSubtree(mesh); 00536 } 00537 regular_mesh = NULL; 00538 } 00539 00540 00541 00542 TEMPLATE 00543 void THIS::copyTree(const ir_mesh_t& mesh) 00544 { 00545 ConstRootIterator the_ele = mesh.beginRootElement(); 00546 ConstRootIterator end_ele = mesh.endRootElement(); 00547 RootIterator the_ele_1 = beginRootElement(); 00548 for (;the_ele != end_ele;++ the_ele, ++ the_ele_1) { 00549 this->copyTree(&*the_ele, &*the_ele_1); 00550 } 00551 } 00552 00553 00554 00555 TEMPLATE 00556 void THIS::copyNonnegtiveSubtree(const ir_mesh_t& mesh) 00557 { 00558 ConstRootIterator the_ele = mesh.beginRootElement(); 00559 ConstRootIterator end_ele = mesh.endRootElement(); 00560 RootIterator the_ele_1 = beginRootElement(); 00561 for (;the_ele != end_ele;++ the_ele, ++ the_ele_1) { 00562 this->copyNonnegtiveSubtree(&*the_ele, &*the_ele_1); 00563 } 00564 } 00565 00566 00567 00568 TEMPLATE 00569 typename THIS::ir_mesh_t& THIS::operator=(const ir_mesh_t& mesh) 00570 { 00571 clear(); 00572 setGeometryTree(mesh.geometry_tree); 00573 copyNonnegtiveSubtree(mesh); 00574 return *this; 00575 } 00576 00577 00578 00579 TEMPLATE 00580 typename THIS::RootFirstIterator THIS::beginRootFirstElement() 00581 { 00582 typename std::list<element_t *>::iterator element = root_element.begin(); 00583 return RootFirstIterator(this, element, *element); 00584 } 00585 00586 00587 00588 TEMPLATE 00589 typename THIS::RootFirstIterator THIS::endRootFirstElement() 00590 { 00591 typename std::list<element_t *>::iterator element = root_element.end(); 00592 return RootFirstIterator(this, element, NULL); 00593 } 00594 00595 00596 00597 TEMPLATE 00598 typename THIS::ActiveIterator THIS::beginActiveElement() 00599 { 00600 RootFirstIterator it = this->beginRootFirstElement(); 00601 while (it->value > 0) ++ it; 00602 return ActiveIterator(it); 00603 } 00604 00605 00606 00607 TEMPLATE 00608 typename THIS::ActiveIterator THIS::endActiveElement() 00609 { 00610 typename std::list<element_t *>::iterator element = root_element.end(); 00611 return ActiveIterator(this, element, NULL); 00612 } 00613 00614 TEMPLATE 00615 void THIS::refineElement(element_t& ele) { 00616 ele.refine(); 00617 ele.value = 1; 00618 for (int k = 0;k < ele.n_child;k ++) { 00619 ele.child[k]->value = 0; 00620 } 00621 } 00622 00623 TEMPLATE 00624 void THIS::renumerateElement() 00625 { 00626 Assert (regular_mesh != NULL, ExcInternalError()); 00627 int i, j, k, l, m, n; 00628 00629 std::cerr << "Renumerating element of the mesh ..." << std::endl; 00630 int n_ele = regular_mesh->n_geometry(DIM); 00631 std::list<int> element_index; 00632 std::vector<std::list<int>::iterator> element_index_iterator(n_ele); 00633 for (i = 0;i < n_ele;i ++) { 00634 element_index_iterator[i] = 00635 element_index.insert(element_index.end(), 00636 i); 00637 } 00638 00639 std::vector<std::list<std::pair<int,std::list<int>::iterator> > > 00640 element_to_node(regular_mesh->n_point()); 00641 for (i = 0;i < n_ele;i ++) { 00642 GeometryBM& ele = regular_mesh->geometry(DIM,i); 00643 std::list<int>::iterator& the_it = element_index_iterator[i]; 00644 for (j = 0;j < ele.n_vertex();j ++) { 00645 element_to_node[ele.vertex(j)].push_back 00646 (std::pair<int,std::list<int>::iterator>(i, the_it)); 00647 } 00648 } 00649 00650 std::vector<int> value(regular_mesh->n_geometry(DIM), 0); 00651 std::vector<int> new_index(regular_mesh->n_geometry(DIM)); 00652 std::list<std::list<int>::iterator> involved_element_index; 00653 for (i = 0, n = -1;i < n_ele;i ++) { 00654 if (involved_element_index.empty()) { 00655 j = element_index.front(); 00656 element_index.pop_front(); 00657 value[j] ++; 00658 } 00659 else { 00660 std::list<std::list<int>::iterator>::iterator 00661 it = involved_element_index.begin(), 00662 end = involved_element_index.end(), 00663 the_it = it; 00664 int n_the_used_vtx = value[**the_it]; 00665 for (;it != end;++ it) { 00666 int& idx = **it; 00667 int& n_used_vtx = value[idx]; 00668 if (regular_mesh->geometry(DIM,idx).n_vertex() == n_used_vtx) { 00669 the_it = it; break; 00670 } else if (n_the_used_vtx < n_used_vtx) { 00671 the_it = it; 00672 n_the_used_vtx = n_used_vtx; 00673 } 00674 } 00675 j = **the_it; 00676 element_index.erase(*the_it); 00677 involved_element_index.erase(the_it); 00678 } 00679 00680 GeometryBM& ele = regular_mesh->geometry(DIM,j); 00681 for (k = 0;k < ele.n_vertex();k ++) { 00682 m = ele.vertex(k); 00683 std::list<std::pair<int,std::list<int>::iterator> >::iterator 00684 it = element_to_node[m].begin(), 00685 end = element_to_node[m].end(); 00686 for (;it != end;++ it) { 00687 if (value[it->first] == 0) { 00688 involved_element_index.push_back(it->second); 00689 } 00690 value[it->first] ++; 00691 } 00692 } 00693 new_index[i] = j; 00694 if (100*i/n_ele > n) { 00695 n = 100*i/n_ele; 00696 std::cerr << "\r" << n << "% OK!"; 00697 } 00698 } 00699 00700 std::vector<GeometryBM> tmp_ele(regular_mesh->geometry(DIM)); 00701 std::vector<int> old_index(n_ele); 00702 #ifdef __SERIALIZATION__ 00703 std::vector<HGeometryBase*> hgp(regular_mesh->h_geometry_ptr[DIM]); 00704 #endif 00705 for (i = 0;i < n_ele;i ++) { 00706 GeometryBM& ele = regular_mesh->geometry(DIM,i); 00707 ele = tmp_ele[new_index[i]]; 00708 ele.index() = i; 00709 old_index[new_index[i]] = i; 00710 #ifdef __SERIALIZATION__ 00711 regular_mesh->h_geometry_ptr[DIM][i] = hgp[new_index[i]]; 00712 #endif 00713 } 00714 ActiveIterator 00715 the_ele = beginActiveElement(), 00716 end_ele = endActiveElement(); 00717 for (;the_ele != end_ele;++ the_ele) { 00718 the_ele->index = old_index[the_ele->index]; 00719 } 00720 std::cerr << " OK!" << std::endl; 00721 } 00722 00723 #undef THIS /// IrregularMesh<DIM,DOW> 00724 00726 00727 TEMPLATE std::ostream& 00728 operator<<(std::ostream& os, const HGeometry<DIM,DOW>& geometry) 00729 { 00730 for (int i = 0;i < geometry.n_boundary;i ++) 00731 os << *(geometry.boundary[i]); 00732 return os; 00733 } 00734 00735 TEMPLATE 00736 std::ostream& operator<<(std::ostream& os, const HElement<DIM, DOW>& element) 00737 { 00738 if (element.value == 1) { 00739 for (int i = 0;i < element.n_child;i ++) { 00740 os << *(element.child[i]); 00741 } 00742 } 00743 else if (element.value == 0) { 00744 os << *(element.h_element); 00745 } 00746 else { 00747 Assert(false, ExcInternalError()); // what happeden? something must be wrong! 00748 } 00749 return os; 00750 } 00751 00752 TEMPLATE 00753 std::ostream& operator<<(std::ostream& os, IrregularMesh<DIM, DOW>& mesh) 00754 { 00755 // IrregularMesh<DIM, DOW>::RootElementIterator the_ele = mesh.beginRootElement(); 00756 // IrregularMesh<DIM, DOW>::RootElementIterator end_ele = mesh.endRootElement(); 00757 // for (;the_ele != end_ele;++ the_ele) 00758 // os << *the_ele; 00760 // ActiveElementIterator<DIM, DOW> the_ele = mesh.beginActiveElement(); 00761 // ActiveElementIterator<DIM, DOW> end_ele = mesh.endActiveElement(); 00762 // for (;the_ele != end_ele;++ the_ele) 00763 // os << *(the_ele->h_element); 00764 return os; 00765 } 00766 00767 TEMPLATE 00768 ElementIterator<DIM, DOW>::ElementIterator() : 00769 mesh(NULL), element(NULL) {} 00770 00772 00773 #define THIS ElementIterator<DIM,DOW> 00774 00775 TEMPLATE 00776 THIS::ElementIterator(IrregularMesh<DIM, DOW> * m, 00777 root_t& r, 00778 HElement<DIM, DOW> * e) : 00779 mesh(m), root_element(r), element(e) {} 00780 00781 TEMPLATE 00782 THIS::ElementIterator(const ElementIterator<DIM, DOW>& it) : 00783 mesh(it.mesh), root_element(it.root_element), element(it.element) {} 00784 00785 TEMPLATE 00786 THIS::~ElementIterator() 00787 {} 00788 00789 TEMPLATE 00790 ElementIterator<DIM, DOW>& THIS::operator=(const ElementIterator<DIM, DOW>& it) 00791 { 00792 mesh = it.mesh; 00793 root_element = it.root_element; 00794 element = it.element; 00795 return *this; 00796 } 00797 00798 #undef THIS /// ElementIterator<DIM,DOW> 00799 00801 00802 TEMPLATE 00803 bool operator==(const ElementIterator<DIM, DOW>& it0, 00804 ElementIterator<DIM, DOW>& it1) 00805 { 00806 return (it0.mesh == it1.mesh && 00807 it0.root_element == it1.root_element && 00808 it0.element == it1.element); 00809 } 00810 00811 TEMPLATE 00812 bool operator!=(const ElementIterator<DIM, DOW>& it0, 00813 ElementIterator<DIM, DOW>& it1) 00814 { 00815 return (it0.mesh != it1.mesh || 00816 it0.root_element != it1.root_element || 00817 it0.element != it1.element); 00818 } 00819 00821 00825 TEMPLATE 00826 RootFirstElementIterator<DIM, DOW>& 00827 RootFirstElementIterator<DIM, DOW>::operator++() 00828 { 00829 if (element == NULL) return *this; // do nothing 00830 else if (element->value == 1) { 00831 element = element->child[0]; 00832 } else { 00833 assert (element->value == 0); 00834 00835 HElement<DIM, DOW> * child = element; 00836 HElement<DIM, DOW> * parent = child->parent; 00837 while (parent != NULL) { 00838 if (child != parent->child[n_child - 1]) break; 00839 child = parent; 00840 parent = child->parent; 00841 }; 00842 if (parent == NULL) { 00843 ++ root_element; 00844 if (root_element == mesh->rootElement().end()) { 00845 element = NULL; 00846 } else { 00847 element = *(root_element); 00848 } 00849 } else { 00850 int i = 0; 00851 for (;child != parent->child[i];++ i) {} 00852 element = parent->child[++ i]; 00853 } 00854 } 00855 return *this; 00856 } 00857 00858 TEMPLATE 00859 ActiveElementIterator<DIM, DOW>& 00860 ActiveElementIterator<DIM, DOW>::operator++() 00861 { 00862 do { 00863 RootFirstElementIterator<DIM, DOW>::operator++(); 00864 if (this->element == NULL) break; 00865 } while (this->element->value > 0); 00866 return *this; 00867 } 00868 00870 00871 #define THIS IrregularMeshPair<DIM, DOW> 00872 00873 TEMPLATE 00874 THIS::IrregularMeshPair(ir_mesh_t& m0, ir_mesh_t& m1) : 00875 mesh0(&m0), mesh1(&m1) 00876 { 00877 Assert(mesh0->geometry_tree == mesh1->geometry_tree, ExcInternalError()); 00878 } 00879 00880 TEMPLATE 00881 THIS::IrregularMeshPair(ir_mesh_t * m0, ir_mesh_t * m1) : 00882 mesh0(m0), mesh1(m1) 00883 { 00884 Assert(mesh0->geometry_tree == mesh1->geometry_tree, ExcInternalError()); 00885 } 00886 00887 TEMPLATE 00888 THIS::~IrregularMeshPair() 00889 {} 00890 00891 TEMPLATE 00892 typename THIS::iterator THIS::beginActiveElementPair() 00893 { 00894 RootFirstElementIterator<DIM, DOW> 00895 it0 = mesh0->beginRootFirstElement(), 00896 it1 = mesh1->beginRootFirstElement(); 00897 00898 while ((*it0).value == 1 && (*it1).value == 1) { 00899 ++ it0; ++ it1; 00900 } 00901 typename iterator::State state; 00902 if ((*it0).value == 0 && (*it1).value == 0) { 00903 state = iterator::EQUAL; 00904 } 00905 else if ((*it0).value == 0) { // then (*it1).value == 1 00906 while ((*it1).value > 0) ++ it1; 00907 state = iterator::GREAT_THAN; 00908 } 00909 else { // then (*it0).value == 1 && (*it1).value == 0 00910 while ((*it0).value > 0) ++ it0; 00911 state = iterator::LESS_THAN; 00912 } 00913 return iterator(this, state, it0, it1); 00914 } 00915 00916 TEMPLATE 00917 typename THIS::iterator THIS::endActiveElementPair() 00918 { 00919 RootFirstElementIterator<DIM, DOW> 00920 it0 = mesh0->endRootFirstElement(), 00921 it1 = mesh1->endRootFirstElement(); 00922 return iterator(this, iterator::EQUAL, it0, it1); 00923 } 00924 00925 #undef THIS /// IrregularMeshPair<DIM, DOW> 00926 00928 00929 #define THIS ActiveElementPairIterator<DIM, DOW> 00930 00931 TEMPLATE 00932 THIS::ActiveElementPairIterator(const THIS& i) : 00933 mesh_pair(i.mesh_pair), st(i.st), 00934 iterator0(i.iterator0), iterator1(i.iterator1) {} 00935 00936 TEMPLATE 00937 THIS& THIS::operator=(const THIS& i) 00938 { 00939 mesh_pair = i.mesh_pair; 00940 st = i.st; 00941 iterator0 = i.iterator0; 00942 iterator1 = i.iterator1; 00943 00944 return *this; 00945 } 00946 00947 TEMPLATE 00948 THIS& THIS::operator++() 00949 { 00950 if (iterator0.element == NULL && iterator1.element == NULL) { 00951 Assert(st == EQUAL, ExcInternalError()); 00952 } else if (st == EQUAL) { 00953 ++ iterator0; 00954 ++ iterator1; 00955 if (iterator0.element == NULL || iterator1.element == NULL) { 00956 Assert (iterator0.element == NULL && iterator1.element == NULL, ExcInternalError()); 00957 return *this; 00958 } 00959 while (iterator0->value > 0 && iterator1->value > 0) { 00960 ++ iterator0; 00961 ++ iterator1; 00962 if (iterator0.element == NULL || iterator1.element == NULL) { 00963 Assert(iterator0.element == NULL && iterator1.element == NULL, ExcInternalError()); 00964 return *this; 00965 } 00966 }; 00967 if (iterator0->value == 0 && iterator1->value == 0) { 00968 st = EQUAL; 00969 } else if (iterator0->value == 0) { // then iterator1->value > 0 00970 while (iterator1->value > 0) ++ iterator1; 00971 st = GREAT_THAN; 00972 } else { // iterator0->value > 0 && iterator1->value == 0 00973 while (iterator0->value > 0) ++ iterator0; 00974 st = LESS_THAN; 00975 } 00976 } else if (st == GREAT_THAN) { 00977 RootFirstElementIterator<DIM, DOW> next0 = iterator0; 00978 ++ next0; 00979 ++ iterator1; 00980 if (iterator1.element == NULL) { 00981 Assert(next0.element == NULL, ExcInternalError()); 00982 iterator0 = next0; 00983 st = EQUAL; 00984 } else if (next0.element == NULL) { 00985 while (iterator1->value > 0) ++ iterator1; 00986 } else if (next0->h_element == iterator1->h_element) { 00987 iterator0 = next0; 00988 while (iterator0->value > 0 && iterator1->value > 0) { 00989 ++ iterator0; 00990 ++ iterator1; 00991 }; 00992 if (iterator0->value == 0 && iterator1->value == 0) { 00993 st = EQUAL; 00994 } else if (iterator0->value == 0) { // then iterator1->value > 0 00995 while (iterator1->value > 0) ++ iterator1; 00996 st = GREAT_THAN; 00997 } else { // iterator0->value > 0 && iterator1->value == 0 00998 while (iterator0->value > 0) ++ iterator0; 00999 st = LESS_THAN; 01000 } 01001 } else { 01002 while (iterator1->value > 0) ++ iterator1; 01003 } 01004 } else { // st == LESS_THAN 01005 RootFirstElementIterator<DIM, DOW> next1 = iterator1; 01006 ++ next1; 01007 ++ iterator0; 01008 if (iterator0.element == NULL) { 01009 Assert(next1.element == NULL, ExcInternalError()); 01010 iterator1 = next1; 01011 st = EQUAL; 01012 } else if (next1.element == NULL) { 01013 while (iterator0->value > 0) ++ iterator0; 01014 } else if (next1->h_element == iterator0->h_element) { 01015 iterator1 = next1; 01016 while (iterator0->value > 0 && iterator1->value > 0) { 01017 ++ iterator0; 01018 ++ iterator1; 01019 }; 01020 if (iterator0->value == 0 && iterator1->value == 0) { 01021 st = EQUAL; 01022 } else if (iterator0->value == 0) { // then iterator1->value > 0 01023 while (iterator1->value > 0) ++ iterator1; 01024 st = GREAT_THAN; 01025 } else { // iterator0->value > 0 && iterator1->value == 0 01026 while (iterator0->value > 0) ++ iterator0; 01027 st = LESS_THAN; 01028 } 01029 } else { 01030 while (iterator0->value > 0) ++ iterator0; 01031 } 01032 } 01033 return *this; 01034 } 01035 01036 TEMPLATE 01037 bool operator==(const THIS& i0, 01038 THIS& i1) 01039 { 01040 return (i0.mesh_pair == i1.mesh_pair && 01041 i0.st == i1.st && 01042 i0.iterator0 == i1.iterator0 && 01043 i0.iterator1 == i1.iterator1); 01044 } 01045 01046 TEMPLATE 01047 bool operator!=(const THIS& i0, 01048 THIS& i1) 01049 { 01050 return (i0.mesh_pair != i1.mesh_pair || 01051 i0.st != i1.st || 01052 i0.iterator0 != i1.iterator0 || 01053 i0.iterator1 != i1.iterator1); 01054 } 01055 01056 #undef THIS /// ActiveElementPairIterator<DIM, DOW> 01057 01059 01060 #define THIS MeshAdaptor<DIM, DOW> 01061 01062 TEMPLATE 01063 THIS::MeshAdaptor() : 01064 from_mesh(NULL), 01065 to_mesh(NULL), 01066 ind(NULL), 01067 convergence_order(1), 01068 refine_step(1), 01069 refine_threshold(1.33333), 01070 coarse_threshold(0.75), 01071 _is_refine_only(false) 01072 {} 01073 01074 01075 01076 TEMPLATE 01077 THIS::MeshAdaptor(ir_mesh_t& f) : 01078 from_mesh(&f), 01079 to_mesh(&f), 01080 ind(NULL), 01081 convergence_order(1), 01082 refine_step(1), 01083 refine_threshold(1.33333), 01084 coarse_threshold(0.75), 01085 _is_refine_only(false) 01086 {} 01087 01088 01089 01090 TEMPLATE 01091 THIS::MeshAdaptor(ir_mesh_t& f, ir_mesh_t& t) : 01092 from_mesh(&f), 01093 to_mesh(&t), 01094 ind(NULL), 01095 convergence_order(1), 01096 refine_step(1), 01097 refine_threshold(1.33333), 01098 coarse_threshold(0.75), 01099 _is_refine_only(false) 01100 {} 01101 01102 01103 01104 TEMPLATE 01105 THIS::~MeshAdaptor() 01106 {} 01107 01108 01109 01110 TEMPLATE 01111 void THIS::globalRefine(unsigned int s) 01112 { 01113 typedef typename ir_mesh_t::ActiveIterator iterator; 01114 std::cerr << "Global refine the mesh ..." << std::endl; 01115 for (unsigned int i = 0;i < s;i ++) { 01116 std::cerr << "\r\tround " << i+1 << " ..." << std::flush; 01117 iterator 01118 the_active_ele = to_mesh->beginActiveElement(), 01119 end_active_ele = to_mesh->endActiveElement(); 01120 for (;the_active_ele != end_active_ele;) { 01121 iterator the_ele = the_active_ele; 01122 ++ the_active_ele; 01123 the_ele->refine(); 01124 the_ele->value = 1; 01125 for (int k = 0;k < the_ele->n_child;k ++) { 01126 the_ele->child[k]->value = 0; 01127 } 01128 } 01129 } 01130 std::cerr << std::endl; 01131 } 01132 01133 01134 01135 TEMPLATE 01136 void THIS::randomRefine(double percent) 01137 { 01138 typedef typename ir_mesh_t::ActiveIterator iterator; 01139 std::cerr << "Randomly refine the mesh ..." << std::endl; 01140 iterator 01141 the_active_ele = to_mesh->beginActiveElement(), 01142 end_active_ele = to_mesh->endActiveElement(); 01143 for (;the_active_ele != end_active_ele;) { 01144 iterator the_ele = the_active_ele; 01145 ++ the_active_ele; 01146 int r = rand(); 01147 if (100.0*r >= RAND_MAX*percent) continue; 01148 01149 the_ele->refine(); 01150 the_ele->value = 1; 01151 for (int k = 0;k < the_ele->n_child;k ++) { 01152 the_ele->child[k]->value = 0; 01153 } 01154 } 01155 std::cerr << std::endl; 01156 } 01157 01158 01159 01160 TEMPLATE 01161 void THIS::adapt() 01162 { 01163 prepareToMesh(); 01164 collectIndicator(); 01165 implementAdaption(); 01166 } 01167 01168 TEMPLATE 01169 void THIS::collectIndicator(HElement<DIM,DOW>& ele, 01170 double convergenceCoefficient) 01171 { 01172 if (ele.value == 0) { 01173 ele.indicator = indicator(ele.index); 01174 } else { 01175 int n_chd = 0; 01176 ele.indicator = 0.0; 01177 for (int i = 0;i < ele.n_child;++ i) { 01178 HElement<DIM,DOW> * p_chd = ele.child[i]; 01179 collectIndicator(*p_chd, convergenceCoefficient); 01180 ele.indicator += p_chd->indicator; 01181 n_chd += 1; 01182 } 01183 ele.indicator *= convergenceCoefficient*ele.n_child/n_chd; 01184 } 01185 } 01186 01187 TEMPLATE 01188 void THIS::collectIndicator() 01189 { 01190 double c = pow(2.0, convergenceOrder()); 01191 typename ir_mesh_t::RootIterator 01192 the_root = to_mesh->beginRootElement(), 01193 end_root = to_mesh->endRootElement(); 01194 for (;the_root != end_root;++ the_root) { 01195 this->collectIndicator(*the_root, c); 01196 } 01197 } 01198 01199 01200 01201 TEMPLATE 01202 void THIS::prepareToMesh() 01203 { 01204 if (from_mesh == to_mesh) return; 01205 *to_mesh = *from_mesh; 01206 } 01207 01208 TEMPLATE 01209 void THIS::adaptElement(HElement<DIM,DOW>& ele, 01210 double convergenceCoefficient, 01211 int refine_level) { 01212 if ((! is_refine_only()) && is_indicator_underflow(ele.indicator)) { 01213 ele.value = 0; 01214 } else if (ele.value == 0 && 01215 refine_level <= refine_step && 01216 is_indicator_overflow(ele.indicator, convergenceCoefficient)) { 01217 ele.refine(); 01218 ele.value = 1; 01219 for (int k = 0;k < ele.n_child;k ++) { 01220 ele.child[k]->value = 0; 01221 01222 ele.child[k]->indicator = ele.indicator/convergenceCoefficient; 01224 adaptElement(*ele.child[k], convergenceCoefficient, refine_level + 1); 01225 } 01226 } else if (ele.value == 1) { 01227 for (int i = 0;i < ele.n_child;i ++) { 01228 adaptElement(*ele.child[i], convergenceCoefficient, refine_level); 01229 } 01230 } 01231 } 01232 01233 TEMPLATE 01234 void THIS::implementAdaption() 01235 { 01236 std::cerr << "Implementing mesh adaption ..." << std::flush; 01237 double convergenceCoefficient = pow(2.0, DIM + convergenceOrder()); 01238 typename ir_mesh_t::RootIterator 01239 the_root = to_mesh->beginRootElement(), 01240 end_root = to_mesh->endRootElement(); 01241 for (;the_root != end_root;++ the_root) { 01242 this->adaptElement(*the_root, convergenceCoefficient, 0); 01243 } 01244 std::cerr << " OK!" << std::endl; 01245 } 01246 01247 #undef THIS /// MeshAdaptor<DIM,DOW> 01248 01250 01251 #define THIS RegularMesh<DIM,DOW> 01252 01253 TEMPLATE 01254 void THIS::renumerateElementHSFC(void (*f)(double,double,double,double&,double&,double&)) 01255 { 01256 std::cerr << "Renumerating element of the mesh using Hibert space filling curve ..." 01257 << std::flush; 01258 int n_ele = this->n_geometry(DIM); 01259 std::vector<double> x(n_ele, 0.0), y(n_ele, 0.0), z(n_ele, 0.0); 01260 for (int i = 0;i < n_ele;++ i) { 01261 GeometryBM& ele = this->geometry(DIM, i); 01262 int n_vtx = ele.n_vertex(); 01263 for (int j = 0;j < n_vtx;++ j) { 01264 int vtx_idx = this->geometry(0, ele.vertex(j)).vertex(0); 01265 afepack::Point<DOW>& pnt = this->point(vtx_idx); 01266 x[i] += pnt[0]; y[i] += pnt[1]; 01267 if (DOW == 3) z[i] += pnt[2]; 01268 } 01269 x[i] /= n_vtx; y[i] /= n_vtx; 01270 if (DOW == 3) z[i] /= n_vtx; 01271 } 01272 std::vector<int> new_index(n_ele); 01273 if (f == NULL) { 01274 hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0]); 01275 } else { 01276 hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0], f); 01277 } 01278 01279 std::vector<GeometryBM> tmp_ele(this->geometry(DIM)); 01280 std::vector<int> old_index(n_ele); 01281 #ifdef __SERIALIZATION__ 01282 std::vector<HGeometryBase*> hgp(h_geometry_ptr[DIM]); 01283 #endif 01284 for (int i = 0;i < n_ele;i ++) { 01285 GeometryBM& ele = this->geometry(DIM,i); 01286 ele = tmp_ele[new_index[i]]; 01287 ele.index() = i; 01288 old_index[new_index[i]] = i; 01289 #ifdef __SERIALIZATION__ 01290 h_geometry_ptr[DIM][i] = hgp[new_index[i]]; 01291 #endif 01292 } 01293 typename ir_mesh_t::ActiveIterator 01294 the_ele = irregularMesh().beginActiveElement(), 01295 end_ele = irregularMesh().endActiveElement(); 01296 for (;the_ele != end_ele;++ the_ele) { 01297 the_ele->index = old_index[the_ele->index]; 01298 } 01299 std::cerr << " OK!" << std::endl; 01300 } 01301 01302 TEMPLATE 01303 void THIS::writeEasyMesh(const std::string& filename) const 01304 { 01305 Assert (DIM == 2, ExcInternalError()); 01306 } 01307 01308 01309 01310 TEMPLATE 01311 void THIS::writeTecplotData(const std::string& filename) const 01312 {} 01313 01314 #undef TEMPLATE 01315 01316 #endif // _HGeometry_templates_h_ 01317