AFEPack
|
00001 00011 #ifndef _FEMSpace_templates_h_ 00012 #define _FEMSpace_templates_h_ 00013 00014 #include "FEMSpace.h" 00015 00016 #define TEMPLATE template <class value_type, int DIM, int DOW, int TDIM> 00017 #define THIS Element<value_type, DIM, DOW, TDIM> 00018 00019 TEMPLATE 00020 THIS::Element(fe_space_t& f) : 00021 sp(&f) 00022 {} 00023 00024 TEMPLATE 00025 THIS::Element(const element_t& e) : 00026 sp(e.sp), 00027 geometry_index(e.geometry_index), 00028 template_element_index(e.template_element_index), 00029 dof_index(e.dof_index) 00030 {} 00031 00032 TEMPLATE 00033 THIS::~Element() 00034 {} 00035 00036 TEMPLATE 00037 void THIS::reinit(fe_space_t& f, 00038 int g, 00039 int t) 00040 { 00041 sp = &f; 00042 geometry_index = g; 00043 template_element_index = t; 00044 } 00045 00046 TEMPLATE 00047 void THIS::reinit(int g, int t) 00048 { 00049 Assert (sp != NULL, ExcInternalError()); 00050 geometry_index = g; 00051 template_element_index = t; 00052 } 00053 00054 TEMPLATE 00055 void THIS::reinit(fe_space_t& f, 00056 int g, 00057 int t, 00058 const std::vector<int>& d) 00059 { 00060 sp = &f; 00061 geometry_index = g; 00062 template_element_index = t; 00063 dof_index = d; 00064 } 00065 00066 TEMPLATE 00067 typename THIS::element_t& THIS::operator=(const element_t& e) 00068 { 00069 if (&e != NULL) { 00070 sp = e.sp; 00071 geometry_index = e.geometry_index; 00072 template_element_index = e.template_element_index; 00073 dof_index = e.dof_index; 00074 } 00075 return *this; 00076 } 00077 00078 TEMPLATE 00079 const typename THIS::fe_space_t& THIS::femSpace() const 00080 { 00081 return *sp; 00082 } 00083 00084 TEMPLATE 00085 typename THIS::fe_space_t& THIS::femSpace() 00086 { 00087 return *sp; 00088 } 00089 00090 TEMPLATE 00091 Mesh<DIM,DOW>& THIS::mesh() const 00092 { 00093 return sp->mesh(); 00094 } 00095 00096 TEMPLATE 00097 int THIS::index() const 00098 { 00099 return geometry_index; 00100 } 00101 00102 TEMPLATE 00103 int& THIS::index() 00104 { 00105 return geometry_index; 00106 } 00107 00108 TEMPLATE 00109 GeometryBM& THIS::geometry() const 00110 { 00111 return sp->mesh().geometry(DIM, geometry_index); 00112 } 00113 00114 TEMPLATE 00115 void THIS::geometry(const GeometryBM& g) 00116 { 00117 geometry_index = g.index(); 00118 } 00119 00120 TEMPLATE 00121 void THIS::geometry(int i) 00122 { 00123 geometry_index = i; 00124 } 00125 00126 TEMPLATE 00127 int THIS::templateElementIndex() const 00128 { 00129 return template_element_index; 00130 } 00131 00132 TEMPLATE 00133 typename THIS::template_t& 00134 THIS::templateElement() const 00135 { 00136 return sp->templateElement(template_element_index); 00137 } 00138 00139 TEMPLATE 00140 void THIS::templateElement(int i) 00141 { 00142 template_element_index = i; 00143 } 00144 00145 TEMPLATE 00146 TemplateGeometry<TDIM>& THIS::templateGeometry() const 00147 { 00148 return sp->templateElement(template_element_index).geometry(); 00149 } 00150 00151 TEMPLATE 00152 const std::vector<int>& THIS::dof() const 00153 { 00154 return dof_index; 00155 } 00156 00157 TEMPLATE 00158 std::vector<int>& THIS::dof() 00159 { 00160 return dof_index; 00161 } 00162 00163 TEMPLATE 00164 unsigned int THIS::n_dof() const 00165 { 00166 return dof_index.size(); 00167 } 00168 00169 TEMPLATE 00170 const std::vector<std::vector<int> >& THIS::geometryImage() const 00171 { 00172 return geo_img; 00173 } 00174 00175 TEMPLATE 00176 std::vector<std::vector<int> >& THIS::geometryImage() 00177 { 00178 return geo_img; 00179 } 00180 00181 TEMPLATE 00182 void THIS::buildGeometryImage() 00183 { 00184 Mesh<DIM,DOW>& m = mesh(); 00185 GeometryBM& geo = this->geometry(); 00186 template_t& t_el = templateElement(); 00187 TemplateGeometry<TDIM>& t_geo = t_el.geometry(); 00188 00189 int i, j, k, l, i1, j1, k1, l1; 00190 00191 geo_img.resize(TDIM+1); 00192 for (i = 0;i <= TDIM;i ++) 00193 geo_img[i].resize(t_geo.n_geometry(i), -1); 00194 geo_img[TDIM][0] = geo.index(); 00195 geo_img[0] = geo.vertex(); 00196 i1 = geo_img[0].size(); 00197 for (i = TDIM-1;i > 0;i --) { 00198 // collect the geometry index in dimension \p{i} 00199 j = t_geo.n_geometry(i+1); 00200 std::vector<int> geo_index(t_geo.n_geometry(i), -1); 00201 for (j1 = 0, l = 0;j1 < j;j1 ++) { 00202 const GeometryBM& geo1 = m.geometry(i+1, geo_img[i+1][j1]); 00203 k = geo1.n_boundary(); 00204 for (k1 = 0;k1 < k;k1 ++) { 00205 for (l1 = 0;l1 < l;l1 ++) 00206 if (geo1.boundary(k1) == geo_index[l1]) break; 00207 if (l1 == l) { 00208 assert(l < t_geo.n_geometry(i)); 00209 geo_index[l ++] = geo1.boundary(k1); 00210 } 00211 } 00212 } 00213 assert(l == t_geo.n_geometry(i)); 00214 // reorder those indices according to their position in the template element 00215 j = t_geo.n_geometry(i); 00216 for (j1 = 0;j1 < j;j1 ++) { 00217 const GeometryBM& g = m.geometry(i, geo_index[j1]); 00218 k = g.n_vertex(); 00219 std::vector<int> v(k); 00220 for (k1 = 0;k1 < k;k1 ++) { 00221 for (l1 = 0;l1 < i1;l1 ++) { 00222 if (g.vertex(k1) == geo_img[0][l1]) { 00223 v[k1] = l1; 00224 break; 00225 } 00226 } 00227 assert(l1 < i1); // ExcMeshData("the vertex is not belong to the element.")); 00228 } 00229 sort(v.begin(), v.end()); 00230 l = t_geo.n_geometry(i); 00231 l1 = g.index(); 00232 for (k = 0;k < l;k ++) { 00233 if (v == t_geo.geometry(i, k).vertex()) { 00234 //geo_img[i][j1] = l1; 00238 geo_img[i][k] = l1; 00239 break; 00240 } 00241 } 00242 assert (k < l); //ExcMeshData("no such boundary in the template geometry.")); 00243 } 00244 } 00245 00250 const u_int& ef = sp->efficiencyFlag(); 00251 for (i = 0;i <= TDIM;i ++) { 00252 if (!(ef & (1L<<i))) { 00253 geo_img[i].clear(); 00254 } 00255 } 00256 } 00257 00258 TEMPLATE 00259 void THIS::lazyBuildGeometryImage() 00260 { 00261 Geometry& geo = geometry(); 00262 geo_img.resize(TDIM+1,std::vector<int>(1, 0)); 00263 geo_img[TDIM].resize(1); 00264 geo_img[TDIM][0] = geo.index(); 00265 geo_img[0] = geo.vertex(); 00266 } 00267 00268 TEMPLATE 00269 void THIS::buildVertexArray(std::vector<afepack::Point<DOW> >& arr) const 00270 { 00271 Mesh<DIM,DOW>& m = mesh(); 00272 Geometry& geo = geometry(); 00273 int n_vertex = geo.n_vertex(); 00274 arr.resize(n_vertex); 00275 for (int i = 0;i < n_vertex;i ++) 00276 arr[i] = m.point(m.geometry(0,geo.vertex(i)).vertex(0)); 00277 } 00278 00279 TEMPLATE 00280 const double ** THIS::buildVertexArray() const 00281 { 00282 Mesh<DIM,DOW>& m = mesh(); 00283 Geometry& geo = geometry(); 00284 int n_vertex = geo.n_vertex(); 00285 typedef double * double_pointer; 00286 const double ** arr = (const double **)new double_pointer[n_vertex]; 00287 for (int i = 0;i < n_vertex;i ++) 00288 arr[i] = m.point(m.geometry(0,geo.vertex(i)).vertex(0)); 00289 return arr; 00290 } 00291 00292 TEMPLATE 00293 const typename BasisFunction<value_type,DIM,TDIM>::Identity& THIS::basis_function_identity(int i) const 00294 { 00295 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00296 const BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00297 return basis_function.identity(); 00298 } 00299 00300 TEMPLATE 00301 value_type THIS::basis_function_value(int i, const afepack::Point<DOW>& p) const 00302 { 00303 std::vector<afepack::Point<DOW> > vertex_array; 00304 buildVertexArray(vertex_array); 00305 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00306 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00307 return basis_function.value(p, vertex_array); 00308 } 00309 00310 TEMPLATE 00311 std::vector<value_type> THIS::basis_function_value(const afepack::Point<DOW>& p) const 00312 { 00313 int i, j; 00314 const double ** vertex_array = buildVertexArray(); 00315 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00316 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00317 i = basis_function.size(); 00318 std::vector<value_type> value(i); 00319 for (j = 0;j < i;j ++) 00320 value[j] = basis_function[j].value(p, vertex_array); 00321 delete[] vertex_array; 00322 return value; 00323 } 00324 00325 TEMPLATE 00326 std::vector<value_type> 00327 THIS::basis_function_value(int i, const std::vector<afepack::Point<DOW> >& p) const 00328 { 00329 std::vector<afepack::Point<DOW> > vertex_array; 00330 buildVertexArray(vertex_array); 00331 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00332 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00333 return basis_function.value(p, vertex_array); 00334 } 00335 00336 TEMPLATE 00337 std::vector<std::vector<value_type> > 00338 THIS::basis_function_value(const std::vector<afepack::Point<DOW> >& p) const 00339 { 00340 int i, i1; 00341 const double ** vertex_array = buildVertexArray(); 00342 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00343 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00344 i = basis_function.size(); 00345 std::vector<std::vector<value_type> > value(i); 00346 for (i1 = 0;i1 < i;i1 ++) 00347 value[i1] = basis_function[i1].value(p, vertex_array); 00348 delete[] vertex_array; 00349 return value; 00350 } 00351 00352 TEMPLATE 00353 std::vector<value_type> THIS::basis_function_gradient(int i, const afepack::Point<DOW>& p) const 00354 { 00355 std::vector<afepack::Point<DOW> > vertex_array; 00356 buildVertexArray(vertex_array); 00357 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00358 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00359 return basis_function.gradient(p, vertex_array); 00360 } 00361 00362 TEMPLATE 00363 std::vector<std::vector<value_type> > 00364 THIS::basis_function_gradient(const afepack::Point<DOW>& p) const 00365 { 00366 int i, j; 00367 const double ** vertex_array = buildVertexArray(); 00368 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00369 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00370 i = basis_function.size(); 00371 std::vector<std::vector<value_type> > value(i); 00372 for (j = 0;j < i;j ++) 00373 value[j] = basis_function[j].gradient(p, vertex_array); 00374 delete[] vertex_array; 00375 return value; 00376 } 00377 00378 TEMPLATE 00379 std::vector<std::vector<value_type> > 00380 THIS::basis_function_gradient(int i, const std::vector<afepack::Point<DOW> >& p) const 00381 { 00382 std::vector<afepack::Point<DOW> > vertex_array; 00383 buildVertexArray(vertex_array); 00384 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00385 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00386 return basis_function.gradient(p, vertex_array); 00387 } 00388 00389 TEMPLATE 00390 std::vector<std::vector<std::vector<value_type> > > 00391 THIS::basis_function_gradient(const std::vector<afepack::Point<DOW> >& p) const 00392 { 00393 int i, i1; 00394 const double ** vertex_array = buildVertexArray(); 00395 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00396 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00397 i = basis_function.size(); 00398 std::vector<std::vector<std::vector<value_type> > > value(i); 00399 for (i1 = 0;i1 < i;i1 ++) 00400 value[i1] = basis_function[i1].gradient(p, vertex_array); 00401 delete[] vertex_array; 00402 return value; 00403 } 00404 00405 #ifdef QUADRATIC_ELEMENT_SUPPORT 00406 00407 TEMPLATE 00408 std::vector<std::vector<value_type> > 00409 THIS::basis_function_hesse(int i, 00410 const afepack::Point<DOW>& p) const 00411 { 00412 std::vector<afepack::Point<DOW> > vertex_array; 00413 buildVertexArray(vertex_array); 00414 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00415 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00416 return basis_function.hesse(p, vertex_array); 00417 } 00418 00419 TEMPLATE 00420 std::vector<std::vector<std::vector<value_type> > > 00421 THIS::basis_function_hesse(const afepack::Point<DOW>& p) const 00422 { 00423 int i, j; 00424 const double ** vertex_array = buildVertexArray(); 00425 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00426 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00427 i = basis_function.size(); 00428 std::vector<std::vector<value_type> > value(i); 00429 for (j = 0;j < i;j ++) 00430 value[j] = basis_function[j].hesse(p, vertex_array); 00431 delete[] vertex_array; 00432 return value; 00433 } 00434 00435 TEMPLATE 00436 std::vector<std::vector<std::vector<value_type> > > 00437 THIS::basis_function_hesse(int i, 00438 const std::vector<afepack::Point<DOW> >& p) const 00439 { 00440 std::vector<afepack::Point<DOW> > vertex_array; 00441 buildVertexArray(vertex_array); 00442 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00443 BasisFunction<value_type,DOW,DIM>& basis_function = t_el.basisFunction(i); 00444 return basis_function.hesse(p, vertex_array); 00445 } 00446 00447 TEMPLATE 00448 std::vector<std::vector<std::vector<std::vector<value_type> > > > 00449 THIS::basis_function_hesse(const std::vector<afepack::Point<DOW> >& p) const 00450 { 00451 int i, i1; 00452 const double ** vertex_array = buildVertexArray(); 00453 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00454 std::vector<BasisFunction<value_type,DOW,DIM> >& basis_function = t_el.basisFunction(); 00455 i = basis_function.size(); 00456 std::vector<std::vector<std::vector<value_type> > > value(i); 00457 for (i1 = 0;i1 < i;i1 ++) 00458 value[i1] = basis_function[i1].hesse(p, vertex_array); 00459 delete[] vertex_array; 00460 return value; 00461 } 00462 00463 #endif // QUADRATIC_ELEMENT_SUPPORT 00464 00465 TEMPLATE 00466 afepack::Point<DOW> THIS::local_to_global(const afepack::Point<TDIM>& lp) const 00467 { 00468 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00469 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00470 std::vector<afepack::Point<DOW> > vertex_array; 00471 buildVertexArray(vertex_array); 00472 return coord_transform.local_to_global(lp, t_el.vertexArray(), vertex_array); 00473 } 00474 00475 TEMPLATE 00476 afepack::Point<TDIM> THIS::global_to_local(const afepack::Point<DOW>& gp) const 00477 { 00478 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00479 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00480 std::vector<afepack::Point<DOW> > vertex_array; 00481 buildVertexArray(vertex_array); 00482 return coord_transform.global_to_local(gp, t_el.vertexArray(), vertex_array); 00483 } 00484 00485 TEMPLATE 00486 double THIS::local_to_global_jacobian(const afepack::Point<TDIM>& lp) const 00487 { 00488 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00489 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00490 std::vector<afepack::Point<DOW> > vertex_array; 00491 buildVertexArray(vertex_array); 00492 return coord_transform.local_to_global_jacobian(lp, t_el.vertexArray(), vertex_array); 00493 } 00494 00495 TEMPLATE 00496 double THIS::global_to_local_jacobian(const afepack::Point<DOW>& gp) const 00497 { 00498 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00499 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00500 std::vector<afepack::Point<DOW> > vertex_array; 00501 buildVertexArray(vertex_array); 00502 return coord_transform.global_to_local_jacobian(gp, t_el.vertexArray(), vertex_array); 00503 } 00504 00505 TEMPLATE 00506 std::vector<afepack::Point<DOW> >THIS::local_to_global(const std::vector<afepack::Point<TDIM> >& lp) const 00507 { 00508 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00509 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00510 std::vector<afepack::Point<DOW> > vertex_array; 00511 buildVertexArray(vertex_array); 00512 return coord_transform.local_to_global(lp, t_el.vertexArray(), vertex_array); 00513 } 00514 00515 TEMPLATE 00516 std::vector<afepack::Point<TDIM> > THIS::global_to_local(const std::vector<afepack::Point<DOW> >& gp) const 00517 { 00518 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00519 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00520 std::vector<afepack::Point<DOW> > vertex_array; 00521 buildVertexArray(vertex_array); 00522 return coord_transform.global_to_local(gp, t_el.vertexArray(), vertex_array); 00523 } 00524 00525 TEMPLATE 00526 std::vector<double> THIS::local_to_global_jacobian(const std::vector<afepack::Point<TDIM> >& lp) const 00527 { 00528 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00529 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00530 std::vector<afepack::Point<DOW> > vertex_array; 00531 buildVertexArray(vertex_array); 00532 return coord_transform.local_to_global_jacobian(lp, t_el.vertexArray(), vertex_array); 00533 } 00534 00535 TEMPLATE 00536 std::vector<double> THIS::global_to_local_jacobian(const std::vector<afepack::Point<DOW> >& gp) const 00537 { 00538 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00539 const CoordTransform<TDIM,DOW>& coord_transform = t_el.coordTransform(); 00540 std::vector<afepack::Point<DOW> > vertex_array; 00541 buildVertexArray(vertex_array); 00542 return coord_transform.global_to_local_jacobian(gp, t_el.vertexArray(), vertex_array); 00543 } 00544 00545 TEMPLATE 00546 std::vector<double> THIS::unitOutNormal(const afepack::Point<DOW>& p, int n) const 00547 { 00548 const double ** vertex_array = buildVertexArray(); 00549 std::vector<double> ret(templateElement().unitOutNormal().value(p, vertex_array, n)); 00550 delete[] vertex_array; 00551 return ret; 00552 } 00553 00554 TEMPLATE 00555 std::vector<std::vector<double> > THIS::unitOutNormal(const std::vector<afepack::Point<DOW> >& p, int n) const 00556 { 00557 const double ** vertex_array = buildVertexArray(); 00558 std::vector<std::vector<double> > ret(templateElement().unitOutNormal().value(p, vertex_array, n)); 00559 delete[] vertex_array; 00560 return ret; 00561 } 00562 00563 TEMPLATE 00564 const QuadratureInfo<TDIM>& THIS::findQuadratureInfo(int i) const 00565 { 00566 TemplateElement<value_type,DOW,TDIM>& t_el = templateElement(); 00567 return t_el.findQuadratureInfo(i); 00568 } 00569 00570 #undef THIS 00571 #undef TEMPLATE 00572 00574 00575 #define TEMPLATE template <class value_type, int DIM, int DOW, int TDIM> 00576 #define THIS FEMSpace<value_type,DIM,DOW,TDIM> 00577 00578 TEMPLATE 00579 THIS::FEMSpace(Mesh<DIM,DOW>& m, 00580 std::vector<template_t>& t) : 00581 msh(&m), 00582 tmp_ele(&t), 00583 effi_flag(0xFFFF) 00584 {} 00585 00586 TEMPLATE 00587 THIS::FEMSpace(const fe_space_t& f) : 00588 msh(f.msh), 00589 tmp_ele(f.tmp_ele), 00590 df(f.df), 00591 df_in(f.df_in), 00592 effi_flag(f.effi_flag) 00593 {} 00594 00595 TEMPLATE 00596 THIS::~FEMSpace() 00597 {} 00598 00599 TEMPLATE 00600 typename THIS::fe_space_t& THIS::operator=(const fe_space_t& f) 00601 { 00602 if (&f != NULL) { 00603 msh = f.msh; 00604 tmp_ele = f.tmp_ele; 00605 df = f.df; 00606 df_in = f.df_in; 00607 effi_flag = f.effi_flag; 00608 } 00609 return *this; 00610 } 00611 00612 TEMPLATE 00613 void THIS::reinit(Mesh<DIM,DOW>& m, 00614 std::vector<template_t>& t) 00615 { 00616 msh = &m; 00617 tmp_ele = &t; 00618 } 00619 00620 #ifdef MULTITHREAD 00621 00622 TEMPLATE 00623 void THIS::buildElement(bool lazy) 00624 { 00625 int n_thread = getThread(); 00626 ThreadManager thread; 00627 for (int rank = 1;rank < n_thread;rank ++) { 00628 thread.start(encapsulate(&THIS:: 00629 threadBuildElement).collectArgs(this, lazy, n_thread, rank)); 00630 } 00631 threadBuildElement(lazy, n_thread, 0); 00632 thread.join(encapsulate(&THIS::threadBuildElement)); 00633 } 00634 00635 TEMPLATE 00636 void THIS::threadBuildElement(bool lazy, 00637 int n_thread, 00638 int rank) 00639 { 00640 /* std::cerr << "Entrying threadBuildElement with " */ 00641 /* << "lazy = " << lazy */ 00642 /* << ", n_thread = " << n_thread */ 00643 /* << ", rank = " << rank */ 00644 /* << std::endl; */ 00645 int n = n_element()/n_thread; 00646 ElementIterator it_el = beginElement(); 00647 it_el += n*rank; 00648 ElementIterator end_el = it_el; 00649 if (rank + 1 == n_thread) 00650 end_el = endElement(); 00651 else 00652 end_el += n; 00653 if (lazy) { 00654 for (;it_el < end_el;++ it_el) { 00655 it_el->lazyBuildGeometryImage(); 00656 } 00657 } 00658 else { 00659 for (;it_el < end_el;++ it_el) { 00660 it_el->buildGeometryImage(); 00661 } 00662 } 00663 } 00664 00665 #else 00666 00667 TEMPLATE 00668 void THIS::buildElement(bool lazy) 00669 { 00670 ElementIterator it_el = beginElement(); 00671 ElementIterator end_el = endElement(); 00672 if (lazy) { 00673 for (;it_el != end_el;++ it_el) { 00674 it_el->lazyBuildGeometryImage(); 00675 } 00676 } 00677 else { 00678 for (;it_el != end_el;++ it_el) { 00679 it_el->buildGeometryImage(); 00680 } 00681 } 00682 } 00683 00684 #endif 00685 00686 TEMPLATE 00687 void THIS::buildDofBoundaryMark() 00688 { 00689 int i, j, k; 00690 00691 const DegreeOfFreedom& d = dof(); 00692 int n = n_dof(); 00693 for (i = 0;i < n;i ++) { 00694 j = d.dof_index[i].dimension; 00695 k = d.dof_index[i].geometry_index; 00696 dofBoundaryMark(i) = mesh().boundaryMark(j, k); 00697 } 00698 } 00699 00700 #ifdef MULTITHREAD 00701 00702 TEMPLATE 00703 void THIS::buildDof() 00704 { 00705 int i, j, k, l, m, n, i1, j1; 00706 std::vector<std::vector<bool> > flag; 00707 00708 std::cerr << "Building degree of freedom for the FEM space ..." << std::endl; 00709 dof().n_geometry_dof.resize(TDIM+1); 00710 dof().geometry_dof.resize(TDIM+1); 00711 flag.resize(TDIM+1); 00712 for (i = 0;i <= TDIM;i ++) { 00718 if (!(efficiencyFlag() & (1L<<i))) continue; 00719 00720 j = mesh().n_geometry(i); 00721 dof().n_geometry_dof[i].resize(j, 0); 00722 flag[i].resize(j, false); 00723 dof().geometry_dof[i].resize(j); 00724 } 00725 00726 dof().n_dof = 0; 00727 00728 pthread_mutex_t mutex; 00729 pthread_mutex_init(&mutex, NULL); 00730 int n_thread = getThread(); 00731 ThreadManager thread; 00732 for (int rank = 1;rank < n_thread;rank ++) { 00733 thread.start(encapsulate(&THIS:: 00734 threadBuildDof0).collectArgs(this, flag, mutex, n_thread, rank)); 00735 } 00736 threadBuildDof0(flag, mutex, n_thread, 0); 00737 thread.join(encapsulate(&THIS::threadBuildDof0)); 00738 pthread_mutex_destroy(&mutex); 00739 00740 dof().dof_index.resize(dof().n_dof); 00741 dofInfo().resize(dof().n_dof); 00742 00743 pthread_mutex_init(&mutex, NULL); 00744 for (int rank = 1;rank < n_thread;rank ++) { 00745 thread.start(encapsulate(&THIS:: 00746 threadBuildDof1).collectArgs(this, flag, mutex, n_thread, rank)); 00747 } 00748 threadBuildDof1(flag, mutex, n_thread, 0); 00749 thread.join(encapsulate(&THIS::threadBuildDof1)); 00750 pthread_mutex_destroy(&mutex); 00751 00752 std::cerr << "\ttotal " << dof().n_dof << " degree of freedom found." << std::endl; 00753 } 00754 00755 TEMPLATE 00756 void THIS::threadBuildDof0(std::vector<std::vector<bool> >& flag, 00757 pthread_mutex_t& mutex, 00758 int n_thread, 00759 int rank) 00760 { 00761 int i, j, k, l, m, n; 00762 n = n_element()/n_thread; 00763 ElementIterator it_el = beginElement(); 00764 it_el += n*rank; 00765 ElementIterator end_el = it_el; 00766 if (rank + 1 == n_thread) 00767 end_el = endElement(); 00768 else 00769 end_el += n; 00770 for (;it_el < end_el;++ it_el) { 00771 TemplateElement<value_type,DOW,TDIM>& t_el = it_el->templateElement(); 00772 const TemplateGeometry<TDIM>& t_el_geo = t_el.geometry(); 00773 const TemplateDOF<TDIM>& t_dof = t_el.dof(); 00774 const std::vector<std::vector<int> >& el_geo_img = it_el->geometryImage(); 00775 it_el->dof().resize(t_dof.n_dof, -1); 00776 for (i = 0;i <= TDIM;i ++) { 00782 if (!(efficiencyFlag() & (1L<<i))) continue; 00783 00784 for (j = 0;j < t_el_geo.n_geometry(i);j ++) { 00785 l = el_geo_img[i][j]; 00786 m = t_dof.n_geometry_dof[i][j]; 00787 pthread_mutex_lock(&mutex); 00788 if (flag[i][l]) { 00789 Assert(dof().n_geometry_dof[i][l] == m, ExcDOFData("DOF number mismatch.")); 00790 } 00791 else { 00792 flag[i][l] = true; 00793 dof().n_geometry_dof[i][l] = m; 00794 dof().geometry_dof[i][l].resize(m,0); 00795 for (n = 0;n < m;n ++) { 00796 dof().geometry_dof[i][l][n] = dof().n_dof ++; 00797 } 00798 } 00799 pthread_mutex_unlock(&mutex); 00800 } 00801 } 00802 } 00803 } 00804 00805 00806 TEMPLATE 00807 void THIS::threadBuildDof1(std::vector<std::vector<bool> >& flag, 00808 pthread_mutex_t& mutex, 00809 int n_thread, 00810 int rank) 00811 { 00812 int i, j, k, l, m, n, i1, j1; 00813 n = n_element()/n_thread; 00814 ElementIterator it_el = beginElement(); 00815 it_el += n*rank; 00816 ElementIterator end_el = it_el; 00817 if (rank + 1 == n_thread) 00818 end_el = endElement(); 00819 else 00820 end_el += n; 00821 for (;it_el < end_el;++ it_el) { 00822 TemplateElement<value_type,DOW,TDIM>& t_el = it_el->templateElement(); 00823 std::vector<int>& el_dof = it_el->dof(); 00824 const TemplateGeometry<TDIM>& t_el_geo = t_el.geometry(); 00825 const TemplateDOF<TDIM>& t_dof = t_el.dof(); 00826 const std::vector<BasisFunction<value_type,DOW,TDIM> >& basis_function = t_el.basisFunction(); 00827 const std::vector<std::vector<int> >& el_geo_img = it_el->geometryImage(); 00828 double h = (mesh().point(mesh().geometry(0, it_el->geometry().vertex(0)).vertex(0)) 00829 - mesh().point(mesh().geometry(0, it_el->geometry().vertex(1)).vertex(0))).length(); 00830 for (i = 0;i <= TDIM;i ++) { 00836 if (!(efficiencyFlag() & (1L<<i))) continue; 00837 00838 for (j = 0;j < t_el_geo.n_geometry(i);j ++) { 00839 l = el_geo_img[i][j]; 00840 m = t_dof.n_geometry_dof[i][j]; 00841 pthread_mutex_lock(&mutex); 00842 if (!flag[i][l]) { 00843 pthread_mutex_unlock(&mutex); 00844 for (n = 0;n < m;n ++) { 00845 k = t_dof.geometry_dof[i][j][n]; 00846 afepack::Point<DOW> ip = it_el->local_to_global(basis_function[k].interpPoint()); 00847 const typename BasisFunction<value_type,DIM,TDIM>::Identity& id = it_el->basis_function_identity(k); 00848 for (i1 = 0;i1 < m;i1 ++) { 00849 j1 = dof().geometry_dof[i][l][i1]; 00850 const afepack::Point<DOW>& ip1 = dofInfo(j1).interp_point; 00851 const typename BasisFunction<value_type,DIM,TDIM>::Identity& id1 = dofInfo(j1).identity; 00852 if ((ip - ip1).length() < (1.e-6)*h && id == id1) { 00853 el_dof[k] = j1; 00854 break; 00855 } 00856 } 00857 Assert(i1 < m, ExcDOFData("DOF data mismatch.")); 00858 } 00859 } 00860 else { 00861 flag[i][l] = false; 00862 for (n = 0;n < m;n ++) { 00863 k = t_dof.geometry_dof[i][j][n]; 00864 j1 = dof().geometry_dof[i][l][n]; 00865 el_dof[k] = j1; 00866 dof().dof_index[j1].dimension = i; 00867 dof().dof_index[j1].geometry_index = l; 00868 dof().dof_index[j1].dof_index = n; 00869 dofInfo(j1).interp_point = it_el->local_to_global(basis_function[k].interpPoint()); 00870 dofInfo(j1).identity = it_el->basis_function_identity(k); 00871 } 00872 pthread_mutex_unlock(&mutex); 00873 } 00874 } 00875 } 00876 } 00877 } 00878 00879 #else 00880 00881 TEMPLATE 00882 void THIS::buildDof() 00883 { 00884 int i, j, k, l, m, n, i1, j1; 00885 std::vector<std::vector<bool> > flag; 00886 00887 std::cerr << "Building degree of freedom for the FEM space ..." << std::endl; 00888 dof().n_geometry_dof.resize(TDIM+1); 00889 dof().geometry_dof.resize(TDIM+1); 00890 flag.resize(TDIM+1); 00891 for (i = 0;i <= TDIM;i ++) { 00897 if (!(efficiencyFlag() & (1L<<i))) continue; 00898 00899 j = mesh().n_geometry(i); 00900 dof().n_geometry_dof[i].resize(j, 0); 00901 flag[i].resize(j, false); 00902 dof().geometry_dof[i].resize(j); 00903 } 00904 00905 dof().n_dof = 0; 00906 ElementIterator it_el = beginElement(); 00907 ElementIterator end_el = endElement(); 00908 for (;it_el != end_el;++ it_el) { 00909 TemplateElement<value_type,DOW,TDIM>& t_el = it_el->templateElement(); 00910 const TemplateGeometry<TDIM>& t_el_geo = t_el.geometry(); 00911 const TemplateDOF<TDIM>& t_dof = t_el.dof(); 00912 const std::vector<std::vector<int> >& el_geo_img = it_el->geometryImage(); 00913 it_el->dof().resize(t_dof.n_dof); 00914 for (i = 0;i <= TDIM;i ++) { 00920 if (!(efficiencyFlag() & (1L<<i))) continue; 00921 00922 for (j = 0;j < t_el_geo.n_geometry(i);j ++) { 00923 l = el_geo_img[i][j]; 00924 m = t_dof.n_geometry_dof[i][j]; 00925 if (flag[i][l]) { 00926 Assert(dof().n_geometry_dof[i][l] == m, ExcDOFData("DOF number mismatch.")); 00927 } 00928 else { 00929 flag[i][l] = true; 00930 dof().n_geometry_dof[i][l] = m; 00931 dof().geometry_dof[i][l].resize(m); 00932 for (n = 0;n < m;n ++) { 00933 dof().geometry_dof[i][l][n] = dof().n_dof ++; 00934 } 00935 } 00936 } 00937 } 00938 } 00939 00940 dof().dof_index.resize(dof().n_dof); 00941 dofInfo().resize(dof().n_dof); 00942 00943 it_el = beginElement(); 00944 for (;it_el != end_el;++ it_el) { 00945 TemplateElement<value_type,DOW,TDIM>& t_el = it_el->templateElement(); 00946 std::vector<int>& el_dof = it_el->dof(); 00947 const TemplateGeometry<TDIM>& t_el_geo = t_el.geometry(); 00948 const TemplateDOF<TDIM>& t_dof = t_el.dof(); 00949 const std::vector<BasisFunction<value_type,DOW,TDIM> >& basis_function = t_el.basisFunction(); 00950 const std::vector<std::vector<int> >& el_geo_img = it_el->geometryImage(); 00951 double h = (mesh().point(mesh().geometry(0, it_el->geometry().vertex(0)).vertex(0)) 00952 - mesh().point(mesh().geometry(0, it_el->geometry().vertex(1)).vertex(0))).length(); 00953 for (i = 0;i <= TDIM;i ++) { 00959 if (!(efficiencyFlag() & (1L<<i))) continue; 00960 00961 for (j = 0;j < t_el_geo.n_geometry(i);j ++) { 00962 l = el_geo_img[i][j]; 00963 m = t_dof.n_geometry_dof[i][j]; 00964 if (!flag[i][l]) { 00965 for (n = 0;n < m;n ++) { 00966 k = t_dof.geometry_dof[i][j][n]; 00967 afepack::Point<DOW> ip = it_el->local_to_global(basis_function[k].interpPoint()); 00968 const typename BasisFunction<value_type,DIM,TDIM>::Identity& id = it_el->basis_function_identity(k); 00969 for (i1 = 0;i1 < m;i1 ++) { 00970 j1 = dof().geometry_dof[i][l][i1]; 00971 const afepack::Point<DOW>& ip1 = dofInfo(j1).interp_point; 00972 const typename BasisFunction<value_type,DIM,TDIM>::Identity& id1 = dofInfo(j1).identity; 00973 if ((ip - ip1).length() < (1.e-6)*h && id == id1) { 00974 el_dof[k] = j1; 00975 break; 00976 } 00977 } 00978 Assert(i1 < m, ExcDOFData("DOF data mismatch.")); 00979 } 00980 } 00981 else { 00982 flag[i][l] = false; 00983 for (n = 0;n < m;n ++) { 00984 k = t_dof.geometry_dof[i][j][n]; 00985 j1 = dof().geometry_dof[i][l][n]; 00986 el_dof[k] = j1; 00987 dof().dof_index[j1].dimension = i; 00988 dof().dof_index[j1].geometry_index = l; 00989 dof().dof_index[j1].dof_index = n; 00990 dofInfo(j1).interp_point = it_el->local_to_global(basis_function[k].interpPoint()); 00991 dofInfo(j1).identity = it_el->basis_function_identity(k); 00992 } 00993 } 00994 } 00995 } 00996 } 00997 std::cerr << "\ttotal " << dof().n_dof << " degree of freedom found." << std::endl; 00998 } 00999 01000 #endif 01001 01002 TEMPLATE 01003 void THIS::updateDofInterpPoint() 01004 { 01005 ElementIterator it_el = beginElement(); 01006 ElementIterator end_el = endElement(); 01007 for (;it_el != end_el;++ it_el) { 01008 TemplateElement<value_type,DOW,TDIM>& t_el = it_el->templateElement(); 01009 std::vector<int>& el_dof = it_el->dof(); 01010 int n_el_dof = el_dof.size(); 01011 const std::vector<BasisFunction<value_type,DOW,TDIM> >& basis_function = t_el.basisFunction(); 01012 for (int j = 0;j < n_el_dof;j ++) { 01013 dofInfo(el_dof[j]).interp_point = it_el->local_to_global(basis_function[j].interpPoint()); 01014 } 01015 } 01016 } 01017 01018 #undef THIS 01019 #undef TEMPLATE 01020 01022 01023 #define TEMPLATE template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01024 #define THIS FEMFunction<value_type,DIM,DOW,TDIM,Number> 01025 01026 TEMPLATE 01027 THIS::FEMFunction(fe_space_t & f) : 01028 sp(&f) 01029 { 01030 if (sp != NULL) 01031 Vector<Number>::reinit(sp->n_dof()); 01032 } 01033 01034 TEMPLATE 01035 THIS::~FEMFunction() 01036 {} 01037 01038 TEMPLATE 01039 typename THIS::fe_space_t& THIS::femSpace() 01040 { 01041 return *sp; 01042 } 01043 01044 TEMPLATE 01045 void THIS::reinit(fe_space_t& f, bool is_bare) 01046 { 01047 sp = &f; 01048 if ((sp != NULL) && (! is_bare)) 01049 Vector<Number>::reinit(sp->n_dof()); 01050 } 01051 01052 TEMPLATE 01053 value_type THIS::value(const afepack::Point<DOW>& p, const element_t& e) const 01054 { 01055 int i, j; 01056 value_type val; 01057 const std::vector<int>& dof = e.dof(); 01058 std::vector<value_type> basis_value = e.basis_function_value(p); 01059 j = dof.size(); 01060 for (i = 0, val = 0;i < j;i ++) { 01061 val += (*this)(dof[i])*basis_value[i]; 01062 } 01063 return val; 01064 } 01065 01066 TEMPLATE 01067 std::vector<value_type> THIS::gradient(const afepack::Point<DOW>& p, 01068 const element_t& e) const 01069 { 01070 int i, j, k; 01071 std::vector<value_type> val(DOW, 0); 01072 const std::vector<int>& dof = e.dof(); 01073 std::vector<std::vector<value_type> > basis_gradient = e.basis_function_gradient(p); 01074 j = dof.size(); 01075 for (i = 0;i < j;i ++) { 01076 for (k = 0;k < DOW;k ++) { 01077 val[k] += (*this)(dof[i])*basis_gradient[i][k]; 01078 } 01079 } 01080 return val; 01081 } 01082 01083 TEMPLATE 01084 std::vector<value_type> THIS::value(const std::vector<afepack::Point<DOW> >& p, 01085 const element_t& e) const 01086 { 01087 int i, j, i1, j1; 01088 i = p.size(); 01089 std::vector<value_type> val(i, 0); 01090 const std::vector<int>& dof = e.dof(); 01091 j = dof.size(); 01092 std::vector<std::vector<value_type> > basis_value = e.basis_function_value(p); 01093 for (i1 = 0;i1 < i;i1 ++) 01094 for (j1 = 0;j1 < j;j1 ++) 01095 val[i1] += (*this)(dof[j1])*basis_value[j1][i1]; 01096 return val; 01097 } 01098 01099 TEMPLATE 01100 std::vector<std::vector<value_type> > 01101 THIS::gradient(const std::vector<afepack::Point<DOW> >& p, 01102 const element_t& e) const 01103 { 01104 int i, j, i1, j1, m; 01105 i = p.size(); 01106 std::vector<std::vector<value_type> > val(i, std::vector<value_type>(DOW, 0)); 01107 const std::vector<int>& dof = e.dof(); 01108 std::vector<std::vector<std::vector<value_type> > > basis_gradient = e.basis_function_gradient(p); 01109 j = dof.size(); 01110 for (i1 = 0;i1 < i;i1 ++) { 01111 for (j1 = 0;j1 < j;j1 ++) 01112 for (m = 0;m < DOW;m ++) 01113 val[i1][m] += (*this)(dof[j1])*basis_gradient[j1][i1][m]; 01114 } 01115 return val; 01116 } 01117 01118 01124 TEMPLATE 01125 value_type THIS::value(const std::vector<value_type>& basis_value, const element_t& e) const 01126 { 01127 int i, j; 01128 value_type val; 01129 const std::vector<int>& dof = e.dof(); 01130 j = dof.size(); 01131 for (i = 0, val = 0;i < j;i ++) { 01132 val += (*this)(dof[i])*basis_value[i]; 01133 } 01134 return val; 01135 } 01136 01137 TEMPLATE 01138 std::vector<value_type> THIS::gradient(const std::vector<std::vector<value_type> >& basis_gradient, 01139 const element_t& e) const 01140 { 01141 int i, j, k; 01142 std::vector<value_type> val(DOW, 0); 01143 const std::vector<int>& dof = e.dof(); 01144 j = dof.size(); 01145 for (i = 0;i < j;i ++) { 01146 for (k = 0;k < DOW;k ++) { 01147 val[k] += (*this)(dof[i])*basis_gradient[i][k]; 01148 } 01149 } 01150 return val; 01151 } 01152 01153 TEMPLATE 01154 std::vector<value_type> THIS::value(const std::vector<std::vector<value_type> >& basis_value, 01155 const element_t& e) const 01156 { 01157 int i, j, i1, j1; 01158 i = basis_value[0].size(); 01159 std::vector<value_type> val(i, 0); 01160 const std::vector<int>& dof = e.dof(); 01161 j = dof.size(); 01162 for (i1 = 0;i1 < i;i1 ++) 01163 for (j1 = 0;j1 < j;j1 ++) 01164 val[i1] += (*this)(dof[j1])*basis_value[j1][i1]; 01165 return val; 01166 } 01167 01168 TEMPLATE 01169 std::vector<std::vector<value_type> > 01170 THIS::gradient(const std::vector<std::vector<std::vector<value_type> > >& basis_gradient, 01171 const element_t& e) const 01172 { 01173 int i, j, i1, j1, m; 01174 i = basis_gradient[0].size(); 01175 std::vector<std::vector<value_type> > val(i, std::vector<value_type>(DIM, 0)); 01176 const std::vector<int>& dof = e.dof(); 01177 j = dof.size(); 01178 for (i1 = 0;i1 < i;i1 ++) { 01179 for (j1 = 0;j1 < j;j1 ++) 01180 for (m = 0;m < DOW;m ++) 01181 val[i1][m] += (*this)(dof[j1])*basis_gradient[j1][i1][m]; 01182 } 01183 return val; 01184 } 01185 01186 // end of the optimization code 01187 01188 TEMPLATE 01189 void THIS::loadData(const std::string& filename) 01190 { 01191 std::ifstream is(filename.c_str()); 01192 Vector<Number>::block_read(is); 01193 is.close(); 01194 } 01195 01196 01197 TEMPLATE 01198 void THIS::writeData(const std::string& filename) 01199 { 01200 std::ofstream os(filename.c_str()); 01201 Vector<Number>::block_write(os); 01202 os.close(); 01203 } 01204 01205 #undef THIS 01206 #undef TEMPLATE 01207 01209 01210 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01211 LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::LocalFEMFunction(Element<value_type,DIM,DOW,TDIM> & e) 01212 { 01213 ele = &e; 01214 if (ele != NULL) 01215 Vector<value_type>::reinit(ele->n_dof()); 01216 } 01217 01218 01219 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01220 Element<value_type,DIM,DOW,TDIM>& LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::element() 01221 { 01222 return *ele; 01223 } 01224 01225 01226 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01227 void LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::reinit(Element<value_type,DIM,DOW,TDIM> & e) 01228 { 01229 ele = &e; 01230 if (ele != NULL) 01231 Vector<value_type>::reinit(ele->n_dof()); 01232 } 01233 01234 01235 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01236 value_type LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::value(const afepack::Point<DOW>& p) const 01237 { 01238 int i, j; 01239 value_type val; 01240 j = Vector<Number>::size(); 01241 std::vector<value_type> basis_value = ele->basis_function_value(p); 01242 for (i = 0, val = 0;i < j;i ++) { 01243 val += (*this)(i)*basis_value[i]; 01244 } 01245 return val; 01246 } 01247 01248 01249 01250 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01251 std::vector<value_type> LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::value(const std::vector<afepack::Point<DOW> >& p) const 01252 { 01253 int i, j, i1, j1; 01254 i = p.size(); 01255 std::vector<value_type> val(i, 0); 01256 j = Vector<Number>::size(); 01257 std::vector<std::vector<value_type> > basis_value = ele->basis_function_value(p); 01258 for (i1 = 0;i1 < i;i1 ++) 01259 for (j1 = 0;j1 < j;j1 ++) 01260 val[i1] += (*this)(j1)*basis_value[j1][i1]; 01261 return val; 01262 } 01263 01264 01265 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01266 std::vector<value_type> LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::gradient(const afepack::Point<DOW>& p) const 01267 { 01268 int i, j, k; 01269 std::vector<value_type> val(DOW, 0); 01270 j = Vector<Number>::size(); 01271 std::vector<std::vector<value_type> > basis_gradient = ele->basis_function_gradient(p); 01272 for (i = 0;i < j;i ++) { 01273 for (k = 0;k < DOW;k ++) { 01274 val[k] += (*this)(i)*basis_gradient[i][k]; 01275 } 01276 } 01277 return val; 01278 } 01279 01280 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01281 std::vector<std::vector<value_type> > 01282 LocalFEMFunction<value_type,DIM,DOW,TDIM,Number>::gradient(const std::vector<afepack::Point<DOW> >& p) const 01283 { 01284 int i, j, i1, j1, m; 01285 i = p.size(); 01286 std::vector<std::vector<value_type> > val(i); 01287 j = Vector<Number>::size(); 01288 std::vector<std::vector<std::vector<value_type> > > basis_gradient = ele->basis_function_gradient(p); 01289 for (i1 = 0;i1 < i;i1 ++) { 01290 val[i1].resize(DOW, 0); 01291 for (j1 = 0;j1 < j;j1 ++) 01292 for (m = 0;m < DOW;m ++) 01293 val[i1][m] += (*this)(j1)*basis_gradient[j1][i1][m]; 01294 } 01295 return val; 01296 } 01297 01298 01300 01301 01302 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01303 void BoundaryConditionAdmin<value_type,DIM,DOW,TDIM,Number>::add(const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& b) 01304 { 01305 if (b.boundaryType() != BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::DIRICHLET) { 01306 std::cerr << "Now we can only apply Dirichlet boundary condition." << std::endl; 01307 Assert(false, ExcInternalError()); 01308 } 01309 if (b.boundaryMark() < 0) { 01310 std::cerr << "We now require a boundary mark to be a positive number." 01311 << std::endl; 01312 Assert(false, ExcInternalError()); 01313 } 01314 typename std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>::iterator 01315 the_bc = std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>::begin(); 01316 typename std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>::iterator 01317 end_bc = std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>::end(); 01318 for (;the_bc != end_bc;the_bc ++) { 01319 if ((*the_bc)->boundaryMark() == b.boundaryMark()) { 01320 std::cerr << "There is a boundary condition for the same boundary mark(" 01321 << b.boundaryMark() 01322 << ") already." 01323 << std::endl; 01324 Assert(false, ExcInternalError()); 01325 } 01326 } 01327 push_back(&b); 01328 int i = index_map.size(); 01329 if (i <= b.boundaryMark()) 01330 for (;i <= b.boundaryMark();i ++) 01331 index_map.push_back(-1); 01332 index_map[b.boundaryMark()] = std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *>::size() - 1; 01333 } 01334 01335 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01336 void BoundaryConditionAdmin<value_type,DIM,DOW,TDIM,Number>::apply(SparseMatrix<double>& A, 01337 Vector<double>& u, 01338 Vector<double>& f, 01339 bool preserve_symmetry) 01340 { 01341 unsigned int i, j, k, l, n_dof; 01342 n_dof = fem_space->n_dof(); 01343 Assert (A.n() == A.m(), ExcDimensionMismatch(A.n(), A.m())); 01344 Assert (A.n() == f.size(), ExcDimensionMismatch(A.n(), f.size())); 01345 Assert (A.n() == u.size(), ExcDimensionMismatch(A.n(), u.size())); 01346 Assert (A.n() == n_dof, ExcDimensionMismatch(A.n(), n_dof)); 01347 typename FEMSpace<double,DIM,DOW,TDIM>::bmark_t bm; 01348 const SparsityPattern& spA = A.get_sparsity_pattern(); 01349 const std::size_t * row_start = spA.get_rowstart_indices(); 01350 const unsigned int * column_index = spA.get_column_numbers(); 01351 for (i = 0;i < n_dof;i ++) { 01352 bm = fem_space->dofBoundaryMark(i); 01353 const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& bc = find(bm); 01354 if (!isValid(bc)) continue; 01355 f(i) = A.diag_element(i)*bc.value(fem_space->dofInfo(i).interp_point); 01356 for (j = row_start[i]+1;j < row_start[i+1];j ++) A.global_entry(j) = 0.0; 01357 if (preserve_symmetry) { 01358 for (j = row_start[i]+1;j < row_start[i+1];j ++) { 01359 k = column_index[j]; 01360 const unsigned int * p = std::find(&column_index[row_start[k]+1], 01361 &column_index[row_start[k+1]], i); 01362 if (p != &column_index[row_start[k+1]]) { 01363 l = p - &column_index[row_start[0]]; 01364 f(k) -= A.global_entry(l) * f(i) / A.diag_element(i); 01365 A.global_entry(l) = 0.; 01366 } 01367 } 01368 } 01369 } 01370 } 01371 01372 01373 01374 template <class value_type, int DIM, int DOW, int TDIM, typename Number> 01375 void BoundaryConditionAdmin<value_type,DIM,DOW,TDIM,Number>::clearEntry(Vector<double>& f) 01376 { 01377 unsigned int i, n_dof; 01378 n_dof = fem_space->n_dof(); 01379 typename FEMSpace<double,DIM,DOW,TDIM>::bmark_t bm; 01380 for (i = 0;i < n_dof;i ++) { 01381 bm = fem_space->dofBoundaryMark(i); 01382 if (bm == 0) continue; 01383 f(i) = 0.0; 01384 } 01385 } 01386 01387 01388 01389 #endif //_FEMSpace_templates_h_ 01390