AFEPack
|
00001 00011 #ifndef _Geometry_templates_h_ 00012 #define _Geometry_templates_h_ 00013 00014 #include <set> 00015 #include "Geometry.h" 00016 00017 #define TEMPLATE template <int DIM> 00018 #define THIS afepack::Point<DIM> 00019 00020 TEMPLATE 00021 THIS::Point() 00022 { 00023 for (int i = 0;i < DIM;i ++) 00024 x[i] = 0; 00025 } 00026 00027 TEMPLATE 00028 THIS::Point(const double * data) 00029 { 00030 for (int i = 0;i < DIM;i ++) 00031 x[i] = data[i]; 00032 } 00033 00034 TEMPLATE 00035 THIS::Point(const THIS& p) 00036 { 00037 int i; 00038 for (i = 0;i < DIM;i ++) 00039 x[i] = p.x[i]; 00040 } 00041 00042 TEMPLATE 00043 THIS::Point(double d, ...) 00044 { 00045 x[0] = d; 00046 va_list ap; 00047 va_start(ap, d); 00048 for (int i = 1;i < DIM;i ++) { 00049 x[i] = va_arg(ap, double); 00050 } 00051 va_end(ap); 00052 } 00053 00054 TEMPLATE 00055 THIS::~Point() 00056 {} 00057 00058 TEMPLATE 00059 THIS::operator const double *() const 00060 { 00061 return x; 00062 } 00063 00064 TEMPLATE 00065 THIS::operator double *() 00066 { 00067 return x; 00068 } 00069 00070 TEMPLATE 00071 THIS& THIS::operator=(const THIS& p) 00072 { 00073 int i; 00074 for (i = 0;i < DIM;i ++) 00075 x[i] = p.x[i]; 00076 return *this; 00077 } 00078 00079 TEMPLATE 00080 const double& THIS::operator[](int i) const 00081 { 00082 return x[i]; 00083 } 00084 00085 TEMPLATE 00086 double& THIS::operator[](int i) 00087 { 00088 return x[i]; 00089 } 00090 00091 TEMPLATE 00092 double THIS::length() const 00093 { 00094 double v = 0.; 00095 for (int i = 0;i < DIM;i ++) 00096 v += x[i]*x[i]; 00097 return sqrt(v); 00098 } 00099 00100 TEMPLATE 00101 THIS& THIS::operator+=(const THIS& p) 00102 { 00103 for (int i = 0;i < DIM;i ++) 00104 x[i] += p.x[i]; 00105 return *this; 00106 } 00107 00108 TEMPLATE 00109 THIS& THIS::operator-=(const THIS& p) 00110 { 00111 for (int i = 0;i < DIM;i ++) 00112 x[i] -= p.x[i]; 00113 return *this; 00114 } 00115 00116 TEMPLATE 00117 THIS& THIS::operator*=(const double& s) 00118 { 00119 for (int i = 0;i < DIM;i ++) 00120 x[i] *= s; 00121 return *this; 00122 } 00123 00124 TEMPLATE 00125 THIS& THIS::operator/=(const double& s) 00126 { 00127 for (int i = 0;i < DIM;i ++) 00128 x[i] /= s; 00129 return *this; 00130 } 00131 00132 TEMPLATE 00133 THIS midpoint(const THIS& p1, const THIS& p2) 00134 { 00135 double p[DIM]; 00136 for (int i = 0;i < DIM;i ++) 00137 p[i] = (p1.x[i] + p2.x[i])/2.; 00138 return p; 00139 } 00140 00141 TEMPLATE 00142 double distance(const THIS& p1, const THIS& p2) 00143 { 00144 double d = 0.0; 00145 for (int i = 0;i < DIM;i ++) 00146 d += (p1.x[i] - p2.x[i])*(p1.x[i] - p2.x[i]); 00147 return sqrt(d); 00148 } 00149 00150 00151 TEMPLATE 00152 THIS barycenter(const std::vector<THIS >& p, const double * w) 00153 { 00154 int i, j, k; 00155 double bc[DIM]; 00156 k = p.size(); 00157 if (w == NULL) { 00158 for (i = 0;i < DIM;i ++) { 00159 bc[i] = 0; 00160 for (j = 0;j < k;j ++) 00161 bc[i] += p[j][i]; 00162 bc[i] /= k; 00163 } 00164 } 00165 else { 00166 double sw = 0; 00167 for (j = 0;j < k;j ++) sw += w[j]; 00168 for (i = 0;i < DIM;i ++) { 00169 bc[i] = 0; 00170 for (j = 0;j < k;j ++) 00171 bc[i] += w[j]*p[j][i]; 00172 bc[i] /= sw; 00173 } 00174 } 00175 return bc; 00176 } 00177 00178 TEMPLATE 00179 THIS operator+(const THIS& p1, const THIS& p2) 00180 { 00181 double p[DIM]; 00182 for (int i = 0;i < DIM;i ++) 00183 p[i] = p1.x[i] + p2.x[i]; 00184 return p; 00185 } 00186 00187 TEMPLATE 00188 THIS operator-(const THIS& p1, const THIS& p2) 00189 { 00190 double p[DIM]; 00191 for (int i = 0;i < DIM;i ++) 00192 p[i] = p1.x[i] - p2.x[i]; 00193 return p; 00194 } 00195 00196 TEMPLATE 00197 std::istream& operator>>(std::istream& is, THIS& p) 00198 { 00199 for (int i = 0;i < DIM;i ++) 00200 is >> p.x[i]; 00201 return is; 00202 } 00203 00204 TEMPLATE 00205 std::ostream& operator<<(std::ostream& os, const THIS& p) 00206 { 00207 for (int i = 0;i < DIM;i ++) 00208 os << p.x[i] << "\t"; 00209 return os; 00210 } 00211 00212 #undef THIS 00213 #undef TEMPLATE 00214 00216 00217 #define TEMPLATE template <int DIM, int DOW> 00218 #define THIS Mesh<DIM,DOW> 00219 00220 TEMPLATE 00221 THIS::Mesh() : 00222 geo(DIM+1) 00223 {} 00224 00225 TEMPLATE 00226 THIS::Mesh(const THIS& m) : 00227 pnt(m.pnt), 00228 geo(m.geo) 00229 {} 00230 00231 TEMPLATE 00232 THIS::~Mesh() 00233 {} 00234 00235 TEMPLATE 00236 THIS& THIS::operator=(const THIS& m) 00237 { 00238 if (&m != NULL) { 00239 pnt = m.pnt; 00240 geo = m.geo; 00241 } 00242 return *this; 00243 } 00244 00245 TEMPLATE 00246 unsigned int THIS::n_point() const 00247 { 00248 return pnt.size(); 00249 } 00250 00251 TEMPLATE 00252 unsigned int THIS::n_geometry(int n) const 00253 { 00254 return geo[n].size(); 00255 } 00256 00257 TEMPLATE 00258 const std::vector<afepack::Point<DOW> >& THIS::point() const 00259 { 00260 return pnt; 00261 } 00262 00263 TEMPLATE 00264 std::vector<afepack::Point<DOW> >& THIS::point() 00265 { 00266 return pnt; 00267 } 00268 00269 TEMPLATE 00270 const afepack::Point<DOW>& THIS::point(int i) const 00271 { 00272 return pnt[i]; 00273 } 00274 00275 TEMPLATE 00276 afepack::Point<DOW>& THIS::point(int i) 00277 { 00278 return pnt[i]; 00279 } 00280 00281 TEMPLATE 00282 const std::vector<std::vector<GeometryBM> >& THIS::geometry() const 00283 { 00284 return geo; 00285 } 00286 00287 TEMPLATE 00288 std::vector<std::vector<GeometryBM> >& THIS::geometry() 00289 { 00290 return geo; 00291 } 00292 00293 TEMPLATE 00294 const std::vector<GeometryBM>& THIS::geometry(int n) const 00295 { 00296 return geo[n]; 00297 } 00298 00299 TEMPLATE 00300 std::vector<GeometryBM>& THIS::geometry(int n) 00301 { 00302 return geo[n]; 00303 } 00304 00305 TEMPLATE 00306 const GeometryBM& THIS::geometry(int i, int j) const 00307 { 00308 return geo[i][j]; 00309 } 00310 00311 TEMPLATE 00312 GeometryBM& THIS::geometry(int i, int j) 00313 { 00314 return geo[i][j]; 00315 } 00316 00317 TEMPLATE 00318 typename THIS::bmark_t THIS::boundaryMark(int i, int j) const 00319 { 00320 return geo[i][j].bm; 00321 } 00322 00323 TEMPLATE 00324 typename THIS::bmark_t& THIS::boundaryMark(int i, int j) 00325 { 00326 return geo[i][j].bm; 00327 } 00328 00329 TEMPLATE 00330 void THIS::renumerateElement() 00331 { 00332 int i, j, k, l, m, n; 00333 00334 std::cerr << "Renumerating element of the mesh ..." << std::endl; 00335 int n_ele = n_geometry(DIM); 00336 std::list<int> element_index; 00337 std::vector<std::list<int>::iterator> element_index_iterator(n_ele); 00338 for (i = 0;i < n_ele;i ++) { 00339 element_index_iterator[i] = 00340 element_index.insert(element_index.end(), 00341 i); 00342 } 00343 00344 std::vector<std::list<std::pair<int,std::list<int>::iterator> > > 00345 element_to_node(n_point()); 00346 for (i = 0;i < n_ele;i ++) { 00347 GeometryBM& ele = geometry(DIM,i); 00348 std::list<int>::iterator& the_it = element_index_iterator[i]; 00349 for (j = 0;j < ele.n_vertex();j ++) { 00350 element_to_node[ele.vertex(j)].push_back 00351 (std::pair<int,std::list<int>::iterator>(i, the_it)); 00352 } 00353 } 00354 00355 std::vector<int> value(n_geometry(DIM), 0); 00356 std::vector<int> new_index(n_geometry(DIM)); 00357 std::list<std::list<int>::iterator> involved_element_index; 00358 for (i = 0, n = -1;i < n_ele;i ++) { 00359 if (involved_element_index.empty()) { 00360 j = element_index.front(); 00361 element_index.pop_front(); 00362 value[j] ++; 00363 } 00364 else { 00365 std::list<std::list<int>::iterator>::iterator 00366 it = involved_element_index.begin(), 00367 end = involved_element_index.end(), 00368 the_it = it; 00369 int n_the_used_vtx = value[**the_it]; 00370 for (;it != end;++ it) { 00371 int& idx = **it; 00372 int& n_used_vtx = value[idx]; 00373 if (geometry(DIM,idx).n_vertex() == n_used_vtx) { 00374 the_it = it; break; 00375 } else if (n_the_used_vtx < n_used_vtx) { 00376 the_it = it; 00377 n_the_used_vtx = n_used_vtx; 00378 } 00379 } 00380 j = **the_it; 00381 element_index.erase(*the_it); 00382 involved_element_index.erase(the_it); 00383 } 00384 00385 GeometryBM& ele = geometry(DIM,j); 00386 for (k = 0;k < ele.n_vertex();k ++) { 00387 m = ele.vertex(k); 00388 std::list<std::pair<int,std::list<int>::iterator> >::iterator 00389 it = element_to_node[m].begin(), 00390 end = element_to_node[m].end(); 00391 for (;it != end;++ it) { 00392 if (value[it->first] == 0) { 00393 involved_element_index.push_back(it->second); 00394 } 00395 value[it->first] ++; 00396 } 00397 } 00398 new_index[i] = j; 00399 if (100*i/n_ele > n) { 00400 n = 100*i/n_ele; 00401 std::cerr << "\r" << n << "% OK!"; 00402 } 00403 } 00404 00405 std::vector<GeometryBM> tmp_ele(geometry(DIM)); 00406 for (i = 0;i < n_ele;i ++) { 00407 GeometryBM& ele = geometry(DIM,i); 00408 ele = tmp_ele[new_index[i]]; 00409 ele.index() = i; 00410 } 00411 std::cerr << " OK!" << std::endl; 00412 } 00413 00414 TEMPLATE 00415 void THIS::readData(const std::string& s) 00416 { 00417 std::cerr << "Reading mesh data file " << s << " ..." << std::endl; 00418 std::ifstream is(s.c_str()); 00419 is >> *this; 00420 is.close(); 00421 } 00422 00423 TEMPLATE 00424 void THIS::writeData(const std::string& s) const 00425 { 00426 std::cerr << "Writing mesh data file " << s << " ..." << std::endl; 00427 std::ofstream os(s.c_str()); 00428 os << *this; 00429 os.close(); 00430 } 00431 00432 TEMPLATE 00433 void THIS::readData1d(const std::string& s) 00434 { 00435 Assert(DIM == 1, ExcMeshData("this method can only be used for 1d.")); 00436 Assert(DOW == 1, ExcMeshData("this method can only be used for 1d.")); 00437 int i, j; 00438 std::ifstream is(s.c_str()); 00439 is >> i; 00440 std::vector<double> buffer(i); 00441 for (j = 0;j < i;j ++) is >> buffer[j]; 00442 is.close(); 00443 sort(buffer.begin(), buffer.end()); 00444 point().resize(i); 00445 for (j = 0;j < i;j ++) point(j)[0] = buffer[j]; 00446 geometry().resize(2); 00447 geometry(0).resize(i); 00448 boundaryMark(0, 0) = 1; 00449 boundaryMark(0, i-1) = 1; 00450 for (j = 0;j < i;j ++) { 00451 geometry(0, j).index() = j; 00452 geometry(0, j).vertex().resize(1,j); 00453 geometry(0, j).boundary().resize(1,j); 00454 } 00455 geometry(1).resize(i-1); 00456 for (j = 0;j < i-1;j ++) { 00457 geometry(1, j).index() = j; 00458 geometry(1, j).vertex().resize(2); 00459 geometry(1, j).vertex(0) = j; 00460 geometry(1, j).vertex(1) = j+1; 00461 geometry(1, j).boundary().resize(2); 00462 geometry(1, j).boundary(0) = j; 00463 geometry(1, j).boundary(1) = j+1; 00464 } 00465 } 00466 00467 00468 TEMPLATE 00469 std::istream& operator>>(std::istream& is, THIS& m) 00470 { 00471 int i, j, k, l; 00472 00473 std::cerr << "\tReading points ... " << std::flush; 00474 is >> i; 00475 m.pnt.resize(i); 00476 for (j = 0;j < i;j ++) 00477 is >> m.pnt[j]; 00478 std::cerr << i << " OK!" << std::endl; 00479 00480 for (i = 0;i <= DIM;i ++) { 00481 std:: cerr << "\tReading " << i << "-dim geometries ... " << std::flush; 00482 GeometryBM temp; 00483 std::vector<GeometryBM>& g = m.geo[i]; 00484 is >> j; 00485 g.resize(j); 00486 for (k = 0;k < j;k ++) { 00487 is >> temp; 00488 l = temp.index(); 00489 Assert(l < j, (typename THIS::ExcMeshData("geometry index error."))); 00490 g[l] = temp; 00491 } 00492 std::cerr << j << " OK!" << std::endl; 00493 } 00494 return is; 00495 } 00496 00497 TEMPLATE 00498 std::ostream& operator<<(std::ostream& os, const THIS& m) 00499 { 00500 int i, j, k; 00501 00502 os.width(12); 00503 os.setf(std::ios::scientific); 00504 std::cerr << "\tWriting points ... " << std::flush; 00505 i = m.pnt.size(); 00506 os << i << "\n"; 00507 for (j = 0;j < i;j ++) 00508 os << m.pnt[j] << "\n"; 00509 std::cerr << i << " OK!" << std::endl; 00510 00511 for (i = 0;i <= DIM;i ++) { 00512 std::cerr << "\tWriting " << i << "-dim geometries ... " << std::flush; 00513 const std::vector<GeometryBM>& g = m.geo[i]; 00514 j = g.size(); 00515 os << "\n" << j << "\n"; 00516 for (k = 0;k < j;k ++) { 00517 os << g[k]; 00518 } 00519 std::cerr << j << " OK!" << std::endl; 00520 } 00521 return os; 00522 } 00523 00525 00526 template <int DIM> 00527 QuadratureInfo<DIM>::QuadratureInfo() 00528 {} 00529 00530 template <int DIM> 00531 QuadratureInfo<DIM>::QuadratureInfo(const QuadratureInfo<DIM>& q) : 00532 alg_acc(q.alg_acc), 00533 pnt(q.pnt), 00534 wei(q.wei) 00535 {} 00536 00537 template <int DIM> 00538 QuadratureInfo<DIM>::~QuadratureInfo() 00539 {} 00540 00541 template <int DIM> 00542 QuadratureInfo<DIM>& QuadratureInfo<DIM>::operator=(const QuadratureInfo<DIM>& q) 00543 { 00544 if (&q != NULL) { 00545 alg_acc = q.alg_acc; 00546 pnt = q.pnt; 00547 wei = q.wei; 00548 } 00549 return *this; 00550 } 00551 00552 template <int DIM> 00553 inline int QuadratureInfo<DIM>::n_quadraturePoint() const 00554 { 00555 return pnt.size(); 00556 } 00557 00558 template <int DIM> 00559 inline int QuadratureInfo<DIM>::algebricAccuracy() const 00560 { 00561 return alg_acc; 00562 } 00563 00564 template <int DIM> 00565 inline int& QuadratureInfo<DIM>::algebricAccuracy() 00566 { 00567 return alg_acc; 00568 } 00569 00570 template <int DIM> 00571 inline const std::vector<afepack::Point<DIM> >& QuadratureInfo<DIM>::quadraturePoint() const 00572 { 00573 return pnt; 00574 } 00575 00576 template <int DIM> 00577 inline std::vector<afepack::Point<DIM> >& QuadratureInfo<DIM>::quadraturePoint() 00578 { 00579 return pnt; 00580 } 00581 00582 template <int DIM> 00583 inline const afepack::Point<DIM>& QuadratureInfo<DIM>::quadraturePoint(int i) const 00584 { 00585 return pnt[i]; 00586 } 00587 00588 00589 template <int DIM> 00590 inline afepack::Point<DIM>& QuadratureInfo<DIM>::quadraturePoint(int i) 00591 { 00592 return pnt[i]; 00593 } 00594 00595 template <int DIM> 00596 inline const std::vector<double>& QuadratureInfo<DIM>::weight() const 00597 { 00598 return wei; 00599 } 00600 00601 template <int DIM> 00602 inline std::vector<double>& QuadratureInfo<DIM>::weight() 00603 { 00604 return wei; 00605 } 00606 00607 template <int DIM> 00608 inline const double& QuadratureInfo<DIM>::weight(int i) const 00609 { 00610 return wei[i]; 00611 } 00612 00613 template <int DIM> 00614 inline double& QuadratureInfo<DIM>::weight(int i) 00615 { 00616 return wei[i]; 00617 } 00618 00619 template <int DIM> 00620 filtering_istream& operator>>(filtering_istream& is, QuadratureInfo<DIM>& q) 00621 { 00622 int i, j; 00623 is >> q.alg_acc; 00624 is >> i; 00625 q.pnt.resize(i); 00626 q.wei.resize(i); 00627 for (j = 0;j < i;j ++) { 00628 is >> q.pnt[j]; 00629 is >> q.wei[j]; 00630 } 00631 return is; 00632 } 00633 00634 template <int DIM> 00635 std::ostream& operator<<(std::ostream& os, const QuadratureInfo<DIM>& q) 00636 { 00637 int i, j; 00638 os << q.alg_acc << "\n"; 00639 i = q.n_quadraturePoint(); 00640 os << i << "\n"; 00641 for (j = 0;j < i;j ++) { 00642 os << q.pnt[j]; 00643 os << q.wei[j] << "\n"; 00644 } 00645 return os; 00646 } 00647 00649 00650 template <int DIM> 00651 QuadratureInfoAdmin<DIM>::QuadratureInfoAdmin() 00652 {} 00653 00654 template <int DIM> 00655 QuadratureInfoAdmin<DIM>::QuadratureInfoAdmin(const QuadratureInfoAdmin<DIM>& q) : 00656 acc_tbl(q.acc_tbl) 00657 {} 00658 00659 template <int DIM> 00660 QuadratureInfoAdmin<DIM>::~QuadratureInfoAdmin() 00661 {} 00662 00663 template <int DIM> 00664 QuadratureInfoAdmin<DIM>& QuadratureInfoAdmin<DIM>::operator=(const QuadratureInfoAdmin<DIM>& q) 00665 { 00666 acc_tbl = q.acc_tbl; 00667 return *this; 00668 } 00669 00670 template <int DIM> 00671 const QuadratureInfo<DIM>& QuadratureInfoAdmin<DIM>::find(int i) const 00672 { 00673 int j = i; 00674 if (acc_tbl[j] == -1) { 00675 for (;acc_tbl[j] == -1 && j < acc_tbl.size();j ++); 00676 if ((unsigned int)j > acc_tbl.size()) { 00677 std::cerr << "no such quadrature info, algebric accuracy: " 00678 << i << std::endl; 00679 abort(); 00680 } 00681 /* std::cerr << "no such quadrature info, algebric accuracy: " */ 00682 /* << (unsigned int)i */ 00683 /* << "\na more accurate one with algebric accuracy " */ 00684 /* << j << " found!" << std::endl; */ 00685 } 00686 return std::vector<QuadratureInfo<DIM> >::operator[](acc_tbl[j]); 00687 } 00688 00689 template <int DIM> 00690 QuadratureInfo<DIM>& QuadratureInfoAdmin<DIM>::find(int i) 00691 { 00692 int j = i, n = acc_tbl.size(); 00693 if (j < n || acc_tbl[j] == -1) { 00694 for (;acc_tbl[j] == -1 && j < n;j ++); 00695 } 00696 if (j == n) { 00697 std::cerr << "no such quadrature info, algebric accuracy: " 00698 << (unsigned int)i << std::endl; 00699 abort(); 00700 } 00701 return std::vector<QuadratureInfo<DIM> >::operator[](acc_tbl[j]); 00702 } 00703 00704 template <int DIM> 00705 filtering_istream& operator>>(filtering_istream& is, QuadratureInfoAdmin<DIM>& q) 00706 { 00707 int i, j, k; 00708 00709 is >> i; 00710 q.resize(i); 00711 for (j = 0, k = -1;j < i;j ++) { 00712 is >> q[j]; 00713 k = std::max(k, q[j].algebricAccuracy()); 00714 } 00715 q.acc_tbl.resize(k+1, -1); 00716 for (j = 0;j < i;j ++) { 00717 q.acc_tbl[q[j].algebricAccuracy()] = j; 00718 } 00719 for (j = k;j >= 0;j --) { 00720 if (q.acc_tbl[j] != -1) 00721 i = q.acc_tbl[j]; 00722 else 00723 q.acc_tbl[j] = i; 00724 } 00725 return is; 00726 } 00727 00728 template <int DIM> 00729 std::ostream& operator<<(std::ostream& os, const QuadratureInfoAdmin<DIM>& q) 00730 { 00731 int i, j; 00732 00733 i = q.size(); 00734 os << i << "\n"; 00735 for (j = 0;j < i;j ++) 00736 os << q[j] << "\n"; 00737 return os; 00738 } 00739 00741 00742 template <int DIM> 00743 TemplateGeometry<DIM>::TemplateGeometry() : 00744 handle(NULL) 00745 {} 00746 00747 template <int DIM> 00748 TemplateGeometry<DIM>::TemplateGeometry(const TemplateGeometry<DIM>& t) : 00749 handle(NULL), 00750 library_name(t.library_name), 00751 volume_function_name(t.volume_function_name), 00752 quad_info(t.quad_info) 00753 {} 00754 00755 template <int DIM> 00756 TemplateGeometry<DIM>::~TemplateGeometry() 00757 { 00758 unloadFunction(); 00759 } 00760 00761 template <int DIM> 00762 TemplateGeometry<DIM>& TemplateGeometry<DIM>::operator=(const TemplateGeometry<DIM>& t) 00763 { 00764 if (&t != NULL) { 00765 library_name = t.library_name; 00766 volume_function_name = t.volume_function_name; 00767 quad_info = t.quad_info; 00768 loadFunction(); 00769 } 00770 return *this; 00771 } 00772 00773 template <int DIM> 00774 void TemplateGeometry<DIM>::loadFunction() 00775 { 00776 unloadFunction(); 00777 std::string temp; 00778 if (library_path.length() == 0) 00779 temp = library_name; 00780 else 00781 temp = library_path + "/" + library_name; 00782 handle = AFEPackDLOpen(temp); 00783 00784 void * symbol = dlsym(handle, volume_function_name.c_str()); 00785 Assert(symbol, ExcLoadFunction(volume_function_name.c_str(), library_name.c_str())); 00786 volume_function = (double (*)(const double **))symbol; 00787 } 00788 00789 template <int DIM> 00790 void TemplateGeometry<DIM>::unloadFunction() 00791 { 00792 if (handle != NULL) { 00793 dlclose(handle); 00794 handle = NULL; 00795 }; 00796 } 00797 00798 template <int DIM> 00799 const QuadratureInfoAdmin<DIM>& TemplateGeometry<DIM>::quadratureInfo() const 00800 { 00801 return quad_info; 00802 } 00803 00804 template <int DIM> 00805 QuadratureInfoAdmin<DIM>& TemplateGeometry<DIM>::quadratureInfo() 00806 { 00807 return quad_info; 00808 } 00809 00810 template <int DIM> 00811 inline const QuadratureInfo<DIM>& TemplateGeometry<DIM>::findQuadratureInfo(int i) const 00812 { 00813 return quad_info.find(i); 00814 } 00815 00816 template <int DIM> 00817 inline double TemplateGeometry<DIM>::volume() const 00818 { 00819 int i, k; 00820 const double ** v; 00821 k = Mesh<DIM,DIM>::n_point(); 00822 v = (const double **)new double[k]; 00823 for (i = 0;i < k;i ++) 00824 v[i] = Mesh<DIM,DIM>::point(i); 00825 double val = (*volume_function)(v); 00826 delete[] v; 00827 return val; 00828 } 00829 00830 template <int DIM> 00831 void TemplateGeometry<DIM>::readData(const std::string& s) 00832 { 00833 library_path = FindAFEPackLibraryFilePath(s); 00834 std::string filename(library_path + "/" + s); 00835 ExpandString(filename); 00836 filtering_istream is; 00837 OpenAFEPackLibraryFile(filename, is); 00838 is >> *this; 00839 } 00840 00841 template <int DIM> 00842 void TemplateGeometry<DIM>::writeData(const std::string& s) const 00843 { 00844 std::ofstream os(s.c_str()); 00845 os << *this; 00846 os.close(); 00847 } 00848 00849 template <int DIM> 00850 inline const std::vector<afepack::Point<DIM> >& TemplateGeometry<DIM>::vertexArray() const 00851 { 00852 return Mesh<DIM,DIM>::point(); 00853 } 00854 00855 template <int DIM> 00856 filtering_istream& operator>>(filtering_istream& is, TemplateGeometry<DIM>& t) 00857 { 00858 int i, j, k, l, m; 00859 00860 is >> t.library_name 00861 >> t.volume_function_name; 00862 t.loadFunction(); 00863 00864 is >> i; 00865 t.point().resize(i); 00866 for (j = 0;j < i;j ++) { 00867 is >> t.point(j); 00868 } 00869 00870 for (i = 0;i <= DIM;i ++) { 00871 Geometry temp; 00872 std::vector<GeometryBM>& g = t.geometry(i); 00873 is >> j; 00874 g.resize(j); 00875 for (k = 0;k < j;k ++) { 00876 is >> temp; 00877 // check if the vertices of the geometry are sorted 00878 const std::vector<int>& v = temp.vertex(); 00879 m = v.size() - 1; 00880 for (l = 0;l < m;l ++) 00881 Assert(v[l] < v[l+1], (typename TemplateGeometry<DIM>::ExcTemplateGeometryData("geometry vertices not sorted."))); 00882 l = temp.index(); 00883 Assert(l < j, (typename TemplateGeometry<DIM>::ExcTemplateGeometryData("geometry index error."))); 00884 dynamic_cast<Geometry&>(g[l]) = temp; 00885 } 00886 } 00887 00888 is >> t.quad_info; 00889 return is; 00890 } 00891 00892 template <int DIM> 00893 std::ostream& operator<<(std::ostream& os, const TemplateGeometry<DIM>& t) 00894 { 00895 int i, j, k; 00896 00897 os << t.library_name << "\n\t" 00898 << t.volume_function_name << "\n"; 00899 00900 os.setf(std::ios::scientific); 00901 i = t.point().size(); 00902 os << i << "\n"; 00903 for (j = 0;j < i;j ++) { 00904 os << t.point(j) << "\n"; 00905 } 00906 os << "\n"; 00907 00908 for (i = 0;i <= DIM;i ++) { 00909 const std::vector<GeometryBM>& g = t.geometry(i); 00910 j = g.size(); 00911 os << j << "\n"; 00912 for (k = 0;k < j;k ++) { 00913 os << dynamic_cast<const Geometry&>(g[k]) << "\n"; 00914 } 00915 os << "\n"; 00916 } 00917 os << "\n"; 00918 00919 os << t.quad_info << "\n"; 00920 return os; 00921 } 00922 00923 #undef THIS 00924 #undef TEMPLATE 00925 00927 00928 #if 0 00929 template <int DIM, int DOW> void 00930 SimplestMesh<DIM,DOW>::generateMesh(Mesh<DIM,DOW>& m) 00931 { 00932 int i, j, k, l, n, p; 00933 std::cerr << "Generate mesh structure from the simplest mesh ..." << std::endl; 00934 GeometryBM g; 00935 n = n_element(); 00936 std::vector<std::vector<int> > geo_img(DIM+1); 00937 m.point() = pnt; 00938 for (i = 0;i <= DIM;i ++) m.geometry(i).clear(); 00939 for (i = 0, p = -1;i < n;i ++) { 00940 const TemplateGeometry<DIM>& t_geo = (*tg)[element(i).template_element]; 00941 geo_img[0].resize(t_geo.n_point()); 00942 g.vertex().resize(1); 00943 g.boundary().resize(1); 00944 for (j = 0;j < t_geo.n_point();j ++) { 00945 g.vertex(0) = elementVertex(i, j); 00946 g.boundary(0) = elementVertex(i, j); 00947 for (l = 0;l < m.n_geometry(0);l ++) 00948 if (m.geometry(0,l).vertex(0) == g.vertex(0)) break; 00949 if (l == m.n_geometry(0)) { 00950 g.index() = l; 00951 m.geometry(0).push_back(g); 00952 } 00953 geo_img[0][j] = l; 00954 } 00955 for (j = 1;j <= DIM;j ++) { 00956 geo_img[j].resize(t_geo.n_geometry(j)); 00957 for (k = 0;k < t_geo.n_geometry(j);k ++) { 00958 g.vertex().resize(t_geo.geometry(j,k).n_vertex()); 00959 g.boundary().resize(t_geo.geometry(j,k).n_boundary()); 00960 for (l = 0;l < g.n_vertex();l ++) 00961 g.vertex(l) = geo_img[0][t_geo.geometry(j,k).vertex(l)]; 00962 for (l = 0;l < g.n_boundary();l ++) 00963 g.boundary(l) = geo_img[j-1][t_geo.geometry(j,k).boundary(l)]; 00964 for (l = m.n_geometry(j);l > 0;l --) { 00965 if (isSame(g, m.geometry(j, l - 1))) { 00966 geo_img[j][k] = l - 1; 00967 break; 00968 } 00969 } 00970 if (l == 0) { 00971 l = m.n_geometry(j); 00972 g.index() = l; 00973 m.geometry(j).push_back(g); 00974 geo_img[j][k] = l; 00975 } 00976 } 00977 } 00978 if (100*i/n > p) { 00979 p = 100*i/n; 00980 std::cerr << "\r" << p << "% OK!" << std::flush; 00981 } 00982 } 00983 std::cerr << "\r"; 00984 00986 for (j = 1;j <= DIM;j ++) { 00987 for (k = 0;k < m.n_geometry(j);k ++) { 00988 GeometryBM& geo = m.geometry(j, k); 00989 for (l = 0;l < geo.n_vertex();l ++) { 00990 n = geo.vertex(l); 00991 geo.vertex(l) = m.geometry(0,n).vertex(0); 00992 } 00993 } 00994 } 00996 for (k = 0;k < m.n_geometry(1);k ++) { 00997 GeometryBM& geo = m.geometry(1, k); 00998 for (l = 0;l < geo.n_boundary();l ++) { 00999 n = geo.boundary(l); 01000 geo.boundary(l) = m.geometry(0,n).vertex(0); 01001 } 01002 } 01004 for (k = 0;k < m.n_geometry(0);k ++) { 01005 m.geometry(0, k).vertex(0) = k; 01006 m.geometry(0, k).boundary(0) = k; 01007 }; 01008 } 01009 #endif 01010 01011 template <int DIM, int DOW> void 01012 SimplestMesh<DIM,DOW>::generateMesh(Mesh<DIM,DOW>& m) 01013 { 01014 int i, j, k, l, n, p; 01015 std::cerr << "Generate mesh structure from the simplest mesh ..." << std::endl; 01016 01021 n = n_element(); 01022 std::vector<std::vector<u_int> > pnt_patch(n_point()); 01023 for (i = 0;i < n;++ i) { 01024 for (j = 0;j < element(i).vertex.size();++ j) { 01025 pnt_patch[elementVertex(i, j)].push_back(i); 01026 } 01027 } 01028 std::vector<std::set<u_int, std::less<u_int> > > ele_patch(n); 01029 for (i = 0;i < n;++ i) { 01030 for (j = 0;j < element(i).vertex.size();++ j) { 01031 ele_patch[i].insert(pnt_patch[elementVertex(i, j)].begin(), 01032 pnt_patch[elementVertex(i, j)].end()); 01033 } 01034 } 01035 pnt_patch.clear(); 01036 01037 std::vector<std::vector<std::vector<int> > > 01038 ele_geo(n, std::vector<std::vector<int> >(DIM + 1)); 01039 01040 GeometryBM g; 01041 m.point() = pnt; 01042 for (i = 0;i <= DIM;i ++) m.geometry(i).clear(); 01043 for (i = 0, p = -1;i < n;i ++) { 01044 std::vector<std::vector<int> >& geo_img = ele_geo[i]; 01045 01046 const TemplateGeometry<DIM>& t_geo = (*tg)[element(i).template_element]; 01047 geo_img[0].resize(t_geo.n_point(), -1); 01048 g.vertex().resize(1); 01049 g.boundary().resize(1); 01050 for (j = 0;j < t_geo.n_point();j ++) { 01051 g.vertex(0) = elementVertex(i, j); 01052 g.boundary(0) = elementVertex(i, j); 01053 01054 bool is_found = false; 01055 int geo_idx; 01056 std::set<u_int, std::less<u_int> >::iterator 01057 the_ele = ele_patch[i].begin(), end_ele = ele_patch[i].end(); 01058 for (;the_ele != end_ele;++ the_ele) { 01059 u_int ele_idx = *the_ele; 01060 if (ele_idx >= i) continue; 01061 for (l = 0;l < ele_geo[ele_idx][0].size();++ l) { 01062 geo_idx = ele_geo[ele_idx][0][l]; 01063 if (geo_idx >=0 && m.geometry(0, geo_idx).vertex(0) == g.vertex(0)) { 01064 is_found = true; 01065 break; 01066 } 01067 } 01068 if (is_found) break; 01069 } 01070 if (!is_found) { 01071 geo_idx = m.n_geometry(0); 01072 g.index() = geo_idx; 01073 m.geometry(0).push_back(g); 01074 } 01075 geo_img[0][j] = geo_idx; 01076 } 01077 for (j = 1;j <= DIM;j ++) { 01078 geo_img[j].resize(t_geo.n_geometry(j)); 01079 for (k = 0;k < t_geo.n_geometry(j);k ++) { 01080 g.vertex().resize(t_geo.geometry(j,k).n_vertex()); 01081 g.boundary().resize(t_geo.geometry(j,k).n_boundary()); 01082 for (l = 0;l < g.n_vertex();l ++) 01083 g.vertex(l) = geo_img[0][t_geo.geometry(j,k).vertex(l)]; 01084 for (l = 0;l < g.n_boundary();l ++) 01085 g.boundary(l) = geo_img[j-1][t_geo.geometry(j,k).boundary(l)]; 01086 01087 bool is_found = false; 01088 int geo_idx; 01089 std::set<u_int, std::less<u_int> >::iterator 01090 the_ele = ele_patch[i].begin(), end_ele = ele_patch[i].end(); 01091 for (;the_ele != end_ele;++ the_ele) { 01092 u_int ele_idx = *the_ele; 01093 if (ele_idx >= i) continue; 01094 for (l = 0;l < ele_geo[ele_idx][j].size();++ l) { 01095 geo_idx = ele_geo[ele_idx][j][l]; 01096 if (geo_idx >= 0 && isSame(m.geometry(j, geo_idx), g)) { 01097 is_found = true; 01098 break; 01099 } 01100 } 01101 if (is_found) break; 01102 } 01103 01104 if (!is_found) { 01105 geo_idx = m.n_geometry(j); 01106 g.index() = geo_idx; 01107 m.geometry(j).push_back(g); 01108 } 01109 geo_img[j][k] = geo_idx; 01110 } 01111 } 01112 if (100*i/n > p) { 01113 p = 100*i/n; 01114 std::cerr << "\r" << p << "% OK!" << std::flush; 01115 } 01116 } 01117 std::cerr << "\r"; 01118 01120 for (j = 1;j <= DIM;j ++) { 01121 for (k = 0;k < m.n_geometry(j);k ++) { 01122 GeometryBM& geo = m.geometry(j, k); 01123 for (l = 0;l < geo.n_vertex();l ++) { 01124 n = geo.vertex(l); 01125 geo.vertex(l) = m.geometry(0,n).vertex(0); 01126 } 01127 } 01128 } 01130 for (k = 0;k < m.n_geometry(1);k ++) { 01131 GeometryBM& geo = m.geometry(1, k); 01132 for (l = 0;l < geo.n_boundary();l ++) { 01133 n = geo.boundary(l); 01134 geo.boundary(l) = m.geometry(0,n).vertex(0); 01135 } 01136 } 01138 for (k = 0;k < m.n_geometry(0);k ++) { 01139 m.geometry(0, k).vertex(0) = k; 01140 m.geometry(0, k).boundary(0) = k; 01141 }; 01142 } 01143 01144 template <int DIM, int DOW> 01145 void Mesh<DIM,DOW>::renumerateElementHSFC(void (*f)(double,double,double,double&,double&,double&)) 01146 { 01147 std::cerr << "Renumerating element of the mesh using Hibert space curve filling ..." 01148 << std::flush; 01149 int n_ele = n_geometry(DIM); 01150 std::vector<double> x(n_ele, 0.0), y(n_ele, 0.0), z(n_ele, 0.0); 01151 for (int i = 0;i < n_ele;++ i) { 01152 GeometryBM& ele = geometry(DIM, i); 01153 int n_vtx = ele.n_vertex(); 01154 for (int j = 0;j < n_vtx;++ j) { 01155 int vtx_idx = geometry(0, ele.vertex(j)).vertex(0); 01156 afepack::Point<DOW>& pnt = point(vtx_idx); 01157 x[i] += pnt[0]; y[i] += pnt[1]; 01158 if (DOW == 3) z[i] += pnt[2]; 01159 } 01160 x[i] /= n_vtx; y[i] /= n_vtx; 01161 if (DOW == 3) z[i] /= n_vtx; 01162 } 01163 std::vector<int> new_index(n_ele); 01164 if (f == NULL) { 01165 hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0]); 01166 } else { 01167 hsfc_renumerate(n_ele, &x[0], &y[0], &z[0], &new_index[0], f); 01168 } 01169 01170 std::vector<GeometryBM> tmp_ele(geometry(DIM)); 01171 for (int i = 0;i < n_ele;i ++) { 01172 GeometryBM& ele = geometry(DIM,i); 01173 ele = tmp_ele[new_index[i]]; 01174 ele.index() = i; 01175 } 01176 std::cerr << " OK!" << std::endl; 01177 } 01178 01179 01180 01181 #endif // Geometry_templates_h_ 01182 01183 // 01184 // end of file