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