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