AFEPack
|
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