AFEPack
DGFEMSpace.templates.h
浏览该文件的文档。
00001 
00002 // DGFEMSpace.templates.h : by R.Lie
00003 //
00004 
00005 #ifndef _DGFEMSpace_templates_h_
00006 #define _DGFEMSpace_templates_h_
00007 
00008 #include "DGFEMSpace.h"
00009 
00010 #define TEMPLATE template <class value_type, int DIM, int DOW, int TDIM, int TDIM1>
00011 #define THIS DGElement<value_type,DIM,DOW,TDIM,TDIM1>
00012 
00013 TEMPLATE
00014 THIS::DGElement(fe_space_t& f) : sp(&f) 
00015 {
00016   neigh[0] = NULL;
00017   neigh[1] = NULL;
00018 };
00019 
00020 TEMPLATE
00021 THIS::DGElement(const dg_element_t& e) :
00022 sp(e.sp),
00023   geometry_index(e.geometry_index),
00024   template_element_index(e.template_element_index),
00025   geo_img(e.geo_img)
00026 {
00027   neigh[0] = e.neigh[0];
00028   neigh[1] = e.neigh[1];
00029 };
00030 
00031 TEMPLATE
00032 THIS::~DGElement()
00033 {};
00034 
00035 TEMPLATE
00036 void THIS::reinit(fe_space_t& f, 
00037                   const int& g, 
00038                   const int& t)
00039 {
00040   sp = &f;
00041   geometry_index = g;
00042   template_element_index = t;
00043 };
00044 
00045 TEMPLATE
00046 void THIS::reinit(const int& g, const int& t)
00047 {
00048   Assert (sp != NULL, ExcInternalError());
00049   geometry_index = g;
00050   template_element_index = t;
00051 };
00052 
00053 TEMPLATE
00054 typename THIS::dg_element_t& THIS::operator=(const dg_element_t& e)
00055 {
00056   if (&e != NULL) {
00057     sp = e.sp;
00058     geometry_index = e.geometry_index;
00059     template_element_index = e.template_element_index;
00060     geo_img = e.geo_img;
00061     neigh[0] = e.neigh[0];
00062     neigh[1] = e.neigh[1];
00063   }
00064   return *this;
00065 };
00066 
00067 TEMPLATE
00068 const typename THIS::fe_space_t& THIS::femSpace() const
00069 {
00070   return *sp;
00071 };
00072 
00073 TEMPLATE
00074 typename THIS::fe_space_t& THIS::femSpace()
00075 {
00076   return *sp;
00077 };
00078 
00079 TEMPLATE
00080 const int& THIS::index() const
00081 {
00082   return geometry_index;
00083 };
00084 
00085 TEMPLATE
00086 int& THIS::index()
00087 {
00088   return geometry_index;
00089 };
00090 
00091 TEMPLATE
00092 const int& THIS::boundaryIndex(const u_int& i) const
00093 {
00094   return bnd_idx[i];
00095 };
00096 
00097 TEMPLATE
00098 int& THIS::boundaryIndex(const u_int& i)
00099 {
00100   return bnd_idx[i];
00101 };
00102 
00103 TEMPLATE
00104 GeometryBM& THIS::geometry() const
00105 {
00106   return sp->mesh().geometry(TDIM1, geometry_index);
00107 };
00108 
00109 TEMPLATE
00110 void THIS::geometry(const GeometryBM& g)
00111 {
00112   geometry_index = g.index();
00113 };
00114 
00115 TEMPLATE
00116 void THIS::geometry(const int& i)
00117 {
00118   geometry_index = i;
00119 };
00120 
00121 TEMPLATE
00122 typename THIS::bmark_t THIS::boundaryMark() const
00123 {
00124   return sp->mesh().geometry(TDIM1, geometry_index).boundaryMark();
00125 };
00126 
00127 TEMPLATE
00128 TemplateDGElement<TDIM1,DOW>& THIS::templateElement() const
00129 {
00130   return sp->templateDGElement(template_element_index);
00131 };
00132 
00133 TEMPLATE
00134 void THIS::templateElement(const int& i)
00135 {
00136   template_element_index = i;
00137 };
00138 
00139 TEMPLATE
00140 TemplateGeometry<TDIM1>& THIS::templateGeometry() const
00141 {
00142   return sp->templateDGElement(template_element_index).geometry();
00143 };
00144 
00145 TEMPLATE
00146 const std::vector<std::vector<int> >& THIS::geometryImage() const
00147 {
00148   return geo_img;
00149 };
00150 
00151 TEMPLATE
00152 std::vector<std::vector<int> >& THIS::geometryImage()
00153 {
00154   return geo_img;
00155 };
00156 
00157 TEMPLATE
00158 void THIS::buildGeometryImage()
00159 {
00160   Mesh<DIM,DOW>& m = sp->mesh();
00161   Geometry& geo = this->geometry();
00162   TemplateDGElement<TDIM1,DOW>& t_el = templateElement();
00163   TemplateGeometry<TDIM1>& t_geo = t_el.geometry();
00164         
00165   int i, j, k, l, i1, j1, k1, l1;
00166         
00167   geo_img.resize(TDIM1+1);
00168   for (i = 0;i <= TDIM1;i ++)
00169     geo_img[i].resize(t_geo.n_geometry(i), -1);
00170   geo_img[TDIM1][0] = geo.index();
00171   geo_img[0] = geo.vertex();
00172   i1 = geo_img[0].size();
00173   for (i = TDIM1-1;i > 0;i --) {
00174     // collect the geometry index in dimension \p{i}
00175     j = t_geo.n_geometry(i+1);
00176     std::vector<int> geo_index(t_geo.n_geometry(i), -1);
00177     for (j1 = 0, l = 0;j1 < j;j1 ++) {
00178       Geometry& geo1 = m.geometry(i+1, geo_img[i+1][j1]);
00179       k = geo1.n_boundary();
00180       for (k1 = 0;k1 < k;k1 ++) {
00181         for (l1 = 0;l1 < l;l1 ++)
00182           if (geo1.boundary(k1) == geo_index[l1]) break;
00183         if (l1 == l) {
00184           Assert(l < t_geo.n_geometry(i), ExcInternalError());
00185           geo_index[l ++] = geo1.boundary(k1);
00186         }
00187       }
00188     }
00189     Assert(l == t_geo.n_geometry(i), ExcInternalError());
00190     // reorder those indices according their position in the template element
00191     j = t_geo.n_geometry(i);
00192     for (j1 = 0;j1 < j;j1 ++) {
00193       const Geometry& g = m.geometry(i, geo_index[j1]);
00194       k = g.n_vertex();
00195       std::vector<int> v(k);
00196       for (k1 = 0;k1 < k;k1 ++) {
00197         for (l1 = 0;l1 < i1;l1 ++) {
00198           if (g.vertex(k1) == geo_img[0][l1]) {
00199             v[k1] = l1;
00200             break;
00201           }
00202         }
00203         Assert(l1 < i1, ExcMeshData("the vertex is not belong to the element."));
00204       }
00205       sort(v.begin(), v.end());
00206       l = t_geo.n_geometry(i);
00207       l1 = g.index();
00208       for (k = 0;k < l;k ++) {
00209         if (v == t_geo.geometry(i, k).vertex()) {
00210           geo_img[i][j1] = l1;
00211           break;
00212         }
00213       }
00214       Assert(k < l,ExcMeshData("no such boundary in the template geometry."));
00215     }
00216   }
00217 };
00218 
00219 TEMPLATE
00220 void THIS::buildVertexArray(std::vector<afepack::Point<DOW> >& arr) const
00221 {
00222   Mesh<DIM,DOW>& m = sp->mesh();
00223   Geometry& geo = geometry();
00224   int n_vertex = geo.n_vertex();
00225   arr.resize(n_vertex);
00226   for (int i = 0;i < n_vertex;i ++)
00227     arr[i] = m.point(m.geometry(0,geo.vertex(i)).vertex(0));
00228 };
00229 
00230 TEMPLATE
00231 const double ** THIS::buildVertexArray() const
00232 {
00233   Mesh<DIM,DOW>& m = sp->mesh();
00234   Geometry& geo = geometry();
00235   int n_vertex = geo.n_vertex();
00236   typedef double * double_pointer;
00237   const double ** arr = (const double **)new double_pointer[n_vertex];
00238   for (int i = 0;i < n_vertex;i ++)
00239     arr[i] = m.point(m.geometry(0,geo.vertex(i)).vertex(0));
00240   return arr;
00241 };
00242 
00243 TEMPLATE
00244 afepack::Point<DOW> THIS::local_to_global(const afepack::Point<TDIM1>& lp) const
00245 {
00246   dg_template_t& t_el = templateElement();
00247   const coord_trans_t& ct = t_el.coordTransform();
00248   std::vector<afepack::Point<DOW> > vertex_array;
00249   buildVertexArray(vertex_array);
00250   return ct.local_to_global(lp, t_el.vertexArray(), vertex_array);
00251 };
00252 
00253 TEMPLATE
00254 afepack::Point<TDIM1> THIS::global_to_local(const afepack::Point<DOW>& gp) const
00255 {
00256   dg_template_t& t_el = templateElement();
00257   const coord_trans_t& ct = t_el.coordTransform();
00258   std::vector<afepack::Point<DOW> > vertex_array;
00259   buildVertexArray(vertex_array);
00260   return ct.global_to_local(gp, t_el.vertexArray(), vertex_array);
00261 };
00262 
00263 TEMPLATE
00264 double THIS::local_to_global_jacobian(const afepack::Point<TDIM1>& lp) const
00265 {
00266   dg_template_t& t_el = templateElement();
00267   const coord_trans_t& ct = t_el.coordTransform();
00268   std::vector<afepack::Point<DOW> > vertex_array;
00269   buildVertexArray(vertex_array);
00270   return ct.local_to_global_jacobian(lp, t_el.vertexArray(), vertex_array);
00271 };
00272 
00273 TEMPLATE
00274 double THIS::global_to_local_jacobian(const afepack::Point<DOW>& gp) const
00275 {
00276   dg_template_t& t_el = templateElement();
00277   const coord_trans_t& ct = t_el.coordTransform();
00278   std::vector<afepack::Point<DOW> > vertex_array;
00279   buildVertexArray(vertex_array);
00280   return ct.global_to_local_jacobian(gp, t_el.vertexArray(), vertex_array);
00281 };
00282 
00283 TEMPLATE
00284 std::vector<afepack::Point<DOW> >THIS::local_to_global(const std::vector<afepack::Point<TDIM1> >& lp) const
00285 {
00286   dg_template_t& t_el = templateElement();
00287   const coord_trans_t& ct = t_el.coordTransform();
00288   std::vector<afepack::Point<DOW> > vertex_array;
00289   buildVertexArray(vertex_array);
00290   return ct.local_to_global(lp, t_el.vertexArray(), vertex_array);
00291 };
00292 
00293 TEMPLATE
00294 std::vector<afepack::Point<TDIM1> > THIS::global_to_local(const std::vector<afepack::Point<DOW> >& gp) const
00295 {
00296   dg_template_t& t_el = templateElement();
00297   const coord_trans_t& ct = t_el.coordTransform();
00298   std::vector<afepack::Point<DOW> > vertex_array;
00299   buildVertexArray(vertex_array);
00300   return ct.global_to_local(gp, t_el.vertexArray(), vertex_array);
00301 };
00302 
00303 TEMPLATE
00304 std::vector<double> THIS::local_to_global_jacobian(const std::vector<afepack::Point<TDIM1> >& lp) const
00305 {
00306   dg_template_t& t_el = templateElement();
00307   const coord_trans_t& ct = t_el.coordTransform();
00308   std::vector<afepack::Point<DOW> > vertex_array;
00309   buildVertexArray(vertex_array);
00310   return ct.local_to_global_jacobian(lp, t_el.vertexArray(), vertex_array);
00311 };
00312 
00313 TEMPLATE
00314 std::vector<double> THIS::global_to_local_jacobian(const std::vector<afepack::Point<DOW> >& gp) const
00315 {
00316   dg_template_t& t_el = templateElement();
00317   const coord_trans_t& ct = t_el.coordTransform();
00318   std::vector<afepack::Point<DOW> > vertex_array;
00319   buildVertexArray(vertex_array);
00320   return ct.global_to_local_jacobian(gp, t_el.vertexArray(), vertex_array);
00321 };
00322 
00323 TEMPLATE
00324 const typename THIS::element_t& THIS::neighbourElement(const int& i) const
00325 {
00326   Assert(i >= 0 && i < 2, ExcInternalError());
00327   return *neigh[i];
00328 };
00329 
00330 TEMPLATE
00331 typename THIS::element_t& THIS::neighbourElement(const int& i)
00332 {
00333   Assert(i >= 0 && i < 2, ExcInternalError());
00334   return *neigh[i];
00335 };
00336 
00337 TEMPLATE
00338 const typename THIS::element_t * THIS::p_neighbourElement(const int& i) const
00339 {
00340   return neigh[i];
00341 };
00342 
00343 TEMPLATE
00344 typename THIS::element_t * THIS::p_neighbourElement(const int& i)
00345 {
00346   return neigh[i];
00347 };
00348 
00349 TEMPLATE
00350 const QuadratureInfo<TDIM1>& THIS::findQuadratureInfo(const int& i) const
00351 {
00352   dg_template_t& t_el = templateElement();
00353   return t_el.findQuadratureInfo(i);
00354 };
00355 
00356 TEMPLATE 
00357 std::vector<double> 
00358 unitOutNormal(const afepack::Point<DIM>& p, 
00359               const Element<value_type,DIM,DOW,TDIM>& ele, 
00360               const DGElement<value_type,DIM,DOW,TDIM,TDIM1>& dgele)
00361 {
00362   if (dgele.neigh[0] == &ele)
00363     return ele.unitOutNormal(p, dgele.bnd_idx[0]);
00364   else if (dgele.neigh[1] == &ele)
00365     return ele.unitOutNormal(p, dgele.bnd_idx[1]);
00366 };
00367 
00368 TEMPLATE 
00369 std::vector<std::vector<double> > 
00370 unitOutNormal(const std::vector<afepack::Point<DIM> >& p, 
00371               const Element<value_type,DIM,DOW,TDIM>& ele, 
00372               const DGElement<value_type,DIM,DOW,TDIM,TDIM1>& dgele)
00373 {
00374   if (dgele.neigh[0] == &ele)
00375     return ele.unitOutNormal(p, dgele.bnd_idx[0]);
00376   else if (dgele.neigh[1] == &ele)
00377     return ele.unitOutNormal(p, dgele.bnd_idx[1]);
00378 };
00379 
00380 #undef THIS
00381 #undef TEMPLATE
00382 
00384 
00385 #define TEMPLATE template <class value_type, int DIM, int DOW, int TDIM, int TDIM1>
00386 #define THIS DGFEMSpace<value_type,DIM,DOW,TDIM,TDIM1>
00387 
00388 TEMPLATE
00389 THIS::DGFEMSpace(Mesh<DIM,DOW>& m,
00390                  std::vector<template_t>& t,
00391                  std::vector<dg_template_t>& tdg) :
00392 base_t(m, t), tmp_dgele(&tdg)
00393 {};
00394 
00395 TEMPLATE
00396 THIS::DGFEMSpace(const fe_space_t& f) :
00397 base_t(f), tmp_dgele(f.tmp_dgele)
00398 {};
00399 
00400 TEMPLATE
00401 THIS::~DGFEMSpace()
00402 {};
00403 
00404 TEMPLATE
00405 THIS& THIS::operator=(const THIS& f)
00406 {
00407   base_t::operator=(f);
00408   if (&f != NULL) 
00409     tmp_dgele = f.tmp_dgele;
00410   return *this;
00411 };
00412 
00413 TEMPLATE
00414 void THIS::reinit(Mesh<DIM,DOW>& m,
00415                   std::vector<template_t>& t,
00416                   std::vector<dg_template_t>& tdg)
00417 {
00418   base_t::reinit(m,t);
00419   tmp_dgele = &tdg;
00420 };
00421 
00422 TEMPLATE
00423 void THIS::buildDGElement()
00424 {
00425   int i, j;
00426   std::vector<int> index(base_t::mesh().n_geometry(DIM-1), -1);
00427   DGElementIterator it_el = beginDGElement();
00428   DGElementIterator end_el = endDGElement();
00429   for (i = 0;it_el != end_el;++ it_el) {
00430     it_el->buildGeometryImage();
00431     it_el->neigh[0] = NULL;
00432     it_el->neigh[1] = NULL;
00433     index[it_el->index()] = i ++;
00434   }
00435   typename base_t::ElementIterator 
00436     the_ele = base_t::beginElement(),
00437     end_ele = base_t::endElement();
00438   for (;the_ele != end_ele;++ the_ele) {
00439     GeometryBM& g = the_ele->geometry();
00440     for (i = 0;i < g.n_boundary();i ++) {
00441       j = index[g.boundary(i)];
00442       if (j != -1) {
00443         if (dgele[j].neigh[0] == NULL) {
00444           dgele[j].neigh[0] = &(*the_ele);
00445           dgele[j].bnd_idx[0] = i;
00446         }
00447         else if (dgele[j].neigh[1] == NULL) {
00448           dgele[j].neigh[1] = &(*the_ele);
00449           dgele[j].bnd_idx[1] = i;
00450         }
00451         else {
00452           Assert(false, ExcInternalError());
00453         }
00454       }
00455     }
00456   }
00457 };
00458 
00459 #undef THIS
00460 #undef TEMPLATE
00461 
00462 #endif //_DGFEMSpace_templates_h_
00463 
00464 //
00465 // end of file