AFEPack
HGeometry.templates.h
浏览该文件的文档。
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