AFEPack
TemplateElement.templates.h
浏览该文件的文档。
00001 
00011 #ifndef _TemplateElement_templates_h_
00012 #define _TemplateElement_templates_h_
00013 
00014 #include "TemplateElement.h"
00015 
00016 template <int DIM>
00017 bool operator==(const BasisFunctionIdentity<DIM>& i0, const BasisFunctionIdentity<DIM>& i1)
00018 {
00019   if (i0.order != i1.order) return false;
00020   for (int i = 0;i < DIM;i ++)
00021     if (i0.alpha[i] != i1.alpha[i])
00022       return false;
00023   if (i0.flag != i1.flag ) return false;
00024   return true;
00025 }
00026 
00028 
00029 template <class value_type, int DIM>
00030   ShapeFunction<value_type,DIM>::ShapeFunction() :
00031   handle(NULL)
00032 {}
00033 
00034 template <class value_type, int DIM>
00035   ShapeFunction<value_type,DIM>::ShapeFunction(const ShapeFunction<value_type,DIM>& s) :
00036   handle(NULL),
00037   library_name(s.library_name),
00038   value_function_name(s.value_function_name),
00039   gradient_function_name(s.gradient_function_name)
00040 {
00041   loadFunction();
00042 }
00043 
00044 template <class value_type, int DIM>
00045   ShapeFunction<value_type,DIM>::~ShapeFunction()
00046 {
00047   unloadFunction();
00048 }
00049 
00050 template <class value_type, int DIM>
00051   ShapeFunction<value_type,DIM>& ShapeFunction<value_type,DIM>::operator=(const ShapeFunction<value_type,DIM>& s)
00052 {
00053   if (&s != NULL) {
00054     library_name = s.library_name;
00055     value_function_name = s.value_function_name;
00056     gradient_function_name = s.gradient_function_name;
00057 #ifdef QUADRATIC_ELEMENT_SUPPORT
00058     hesse_function_name = s.hesse_function_name;
00059 #endif // QUADRATIC_ELEMENT_SUPPORT
00060   }
00061   return *this;
00062 }
00063 
00064   template <class value_type, int DIM>
00065   void ShapeFunction<value_type,DIM>::loadFunction()
00066   {
00067     unloadFunction();
00068         
00069     std::string temp;
00070     if (library_path.length() == 0)
00071       temp = library_name;
00072     else
00073       temp = library_path + "/" + library_name;
00074     handle = AFEPackDLOpen(temp);
00075     if (handle == NULL) return;
00076 
00077     void * symbol = dlsym(handle, value_function_name.c_str());
00078     Assert(symbol, ExcLoadFunction(value_function_name.c_str(), library_name.c_str()));
00079     value_function = (void (*)(const double *, const double **, void *))symbol;
00080 
00081     symbol = dlsym(handle, gradient_function_name.c_str());
00082     Assert(symbol, ExcLoadFunction(gradient_function_name.c_str(), library_name.c_str()));
00083     gradient_function = (void (*)(const double *, const double **, void *))symbol;
00084 
00085 #ifdef QUADRATIC_ELEMENT_SUPPORT
00086     symbol = dlsym(handle, hesse_function_name.c_str());
00087     Assert(symbol, ExcLoadFunction(hesse_function_name.c_str(), library_name.c_str()));
00088     hesse_function = (void (*)(const double *, const double **, void *))symbol;
00089 #endif // QUADRATIC_ELEMENT_SUPPORT
00090   }
00091 
00092 template <class value_type, int DIM>
00093     void ShapeFunction<value_type,DIM>::unloadFunction()
00094 {
00095   if (handle != NULL) {
00096     dlclose(handle);
00097     handle = NULL;
00098   }
00099 }
00100 
00101 template <class value_type, int DIM>
00102   inline value_type ShapeFunction<value_type,DIM>::value(const afepack::Point<DIM>& p, 
00103                                                          const std::vector<afepack::Point<DIM> >& v) const
00104 {
00105   value_type val;
00106   int k = v.size();
00107   const double * v1[k];
00108   for (int i = 0;i < k;i ++) v1[i] = v[i];
00109   (*value_function)(p, v1, (void *)(&val));
00110   return val;
00111 }
00112 
00113 template <class value_type, int DIM>
00114   inline std::vector<value_type> ShapeFunction<value_type,DIM>::gradient(const afepack::Point<DIM>& p, 
00115                                                                          const std::vector<afepack::Point<DIM> >& v) const
00116 {
00117   int k = v.size();
00118   const double * v1[k];
00119   for (int i = 0;i < k;i ++) v1[i] = v[i];
00120   std::vector<value_type> val(DIM);
00121   (*gradient_function)(p, v1, (void *)(&val[0]));
00122   return val;
00123 }
00124 
00125 template <class value_type, int DIM>
00126   inline std::vector<value_type> 
00127   ShapeFunction<value_type,DIM>::value(const std::vector<afepack::Point<DIM> >& p, 
00128                                        const std::vector<afepack::Point<DIM> >& v) const
00129 {
00130   int k = v.size();
00131   int l = p.size();
00132   const double * v1[k];
00133   for (int i = 0;i < k;i ++) v1[i] = v[i];
00134   std::vector<value_type> val(l);
00135   for (int i = 0;i < l;i ++) {
00136     (*value_function)(p[i], v1, (void *)(&val[i]));
00137   }
00138   return val;
00139 }
00140 
00141 template <class value_type, int DIM>
00142   inline std::vector<std::vector<value_type> >
00143   ShapeFunction<value_type,DIM>::gradient(const std::vector<afepack::Point<DIM> >& p, 
00144                                           const std::vector<afepack::Point<DIM> >& v) const
00145 {
00146   int k = v.size();
00147   int l = p.size();
00148   const double * v1[k];
00149   for (int i = 0;i < k;i ++) v1[i] = v[i];
00150   std::vector<std::vector<value_type> > val(l,std::vector<value_type>(DIM));
00151   for (int i = 0;i < l;i ++) {
00152     (*gradient_function)(p[i], v1, (void *)(&val[i][0]));
00153   }
00154   return val;
00155 }
00156 
00157 template <class value_type, int DIM>
00158   inline value_type ShapeFunction<value_type,DIM>::value(const afepack::Point<DIM>& p, const double ** v1) const
00159 {
00160   value_type val;
00161   (*value_function)(p, v1, (void *)(&val));
00162   return val;
00163 }
00164 
00165 template <class value_type, int DIM>
00166   inline std::vector<value_type> 
00167   ShapeFunction<value_type,DIM>::gradient(const afepack::Point<DIM>& p, const double ** v1) const
00168 {
00169   std::vector<value_type> val(DIM);
00170   (*gradient_function)(p, v1, (void *)(&val[0]));
00171   return val;
00172 }
00173 
00174 template <class value_type, int DIM>
00175   inline std::vector<value_type> 
00176   ShapeFunction<value_type,DIM>::value(const std::vector<afepack::Point<DIM> >& p, const double ** v1) const
00177 {
00178   int l = p.size();
00179   std::vector<value_type> val(l);
00180   for (int i = 0;i < l;i ++) {
00181     (*value_function)(p[i], v1, (void *)(&val[i]));
00182   }
00183   return val;
00184 }
00185 
00186 template <class value_type, int DIM>
00187   inline std::vector<std::vector<value_type> >
00188   ShapeFunction<value_type,DIM>::gradient(const std::vector<afepack::Point<DIM> >& p, const double ** v1) const
00189 {
00190   int l = p.size();
00191   std::vector<std::vector<value_type> > val(l, std::vector<value_type>(DIM));
00192   for (int i = 0;i < l;i ++) {
00193     (*gradient_function)(p[i], v1, (void *)(&val[i][0]));
00194   }
00195   return val;
00196 }
00197 
00198 #ifdef QUADRATIC_ELEMENT_SUPPORT
00199 
00200 template <class value_type, int DIM>
00201   inline std::vector<std::vector<value_type> > 
00202   ShapeFunction<value_type,DIM>::hesse(const afepack::Point<DIM>& p, 
00203                                        const std::vector<afepack::Point<DIM> >& v) const
00204 {
00205   value_type val[DIM*DIM];
00206   const double ** v1;
00207         
00208   int i, k;
00209   k = v.size();
00210   typedef double * double_pointer;
00211   v1 = (const double **)new double_pointer[k];
00212   for (i = 0;i < k;i ++) v1[i] = v[i];
00213   (*hesse_function)(p, v1, (void *)(&val));
00214   delete[] v1;
00215         
00216   std::vector<std::vector<value_type> > 
00217     result(DIM, 
00218            std::vector<value_type>(DIM));
00219   for (i = 0;i < DIM;i ++) 
00220     for (k = 0;k < DIM;k ++)
00221       result[i][k] = val[i*DIM+k];
00222   return result;
00223 }
00224 
00225 
00226 template <class value_type, int DIM>
00227   inline std::vector<std::vector<std::vector<value_type> > >
00228   ShapeFunction<value_type,DIM>::hesse(const std::vector<afepack::Point<DIM> >& p, 
00229                                        const std::vector<afepack::Point<DIM> >& v) const
00230 {
00231   value_type val[DIM*DIM];
00232   const double ** v1;
00233         
00234   int i, j, k, l;
00235   k = v.size();
00236   l = p.size();
00237   typedef double * double_pointer;
00238   v1 = (const double **)new double_pointer[k];
00239   for (i = 0;i < k;i ++) v1[i] = v[i];
00240   std::vector<std::vector<std::vector<value_type> > > 
00241     result(l,
00242            std::vector<std::vector<value_type> >(DIM, 
00243                                                  std::vector<value_type>(DIM)));
00244   for (i = 0;i < l;i ++) {
00245     (*hesse_function)(p[i], v1, (void *)(&val));
00246     for (j = 0;j < DIM;j ++)
00247       for (m = 0;m < DIM;m ++)
00248         result[i][j][m] = val[j*DIM+m];
00249   }
00250   delete[] v1;  
00251   return result;
00252 }
00253 
00254 template <class value_type, int DIM>
00255   inline std::vector<std::vector<value_type> >
00256   ShapeFunction<value_type,DIM>::hesse(const afepack::Point<DIM>& p, 
00257                                        const double ** v1) const
00258 {
00259   value_type val[DIM*DIM];
00260   (*hesse_function)(p, v1, (void *)(&val));
00261   std::vector<std::vector<value_type> > result(DIM, std::vector<value_type>(DIM));
00262   for (int i = 0;i < DIM;i ++)
00263     for (int j = 0;j < DIM;j ++)
00264       result[i][j] = val[i*DIM+j];
00265   return result;
00266 }
00267 
00268 template <class value_type, int DIM>
00269   inline std::vector<std::vector<std::vector<value_type> > >
00270   ShapeFunction<value_type,DIM>::hesse(const std::vector<afepack::Point<DIM> >& p, 
00271                                        const double ** v1) const
00272 {
00273   value_type val[DIM*DIM];
00274         
00275   int i, j, l;
00276   l = p.size();
00277   std::vector<std::vector<std::vector<value_type> > > 
00278     result(l,
00279            std::vector<std::vector<value_type> >(DIM, 
00280                                                  std::vector<value_type>(DIM)));
00281   for (i = 0;i < l;i ++) {
00282     (*hesse_function)(p[i], v1, (void *)(&val));
00283     for (j = 0;j < DIM;j ++)
00284       for (m = 0;m < DIM;m ++)
00285         result[i][j][m] = val[j*DIM+m];
00286   }
00287   return result;
00288 }
00289 
00290 #endif // QUADRATIC_ELEMENT_SUPPORT
00291 
00293 
00294 template <class value_type, int DIM, int TDIM>
00295   BasisFunction<value_type,DIM,TDIM>::BasisFunction()
00296 {}
00297 
00298 template <class value_type, int DIM, int TDIM>
00299   BasisFunction<value_type,DIM,TDIM>::BasisFunction(const afepack::Point<TDIM>& p) :
00300   ip(p)
00301 {}
00302 
00303 template <class value_type, int DIM, int TDIM>
00304   BasisFunction<value_type,DIM,TDIM>::BasisFunction(const BasisFunction<value_type,DIM,TDIM>& b) :
00305   ip(b.ip), id(b.id)
00306 {}
00307 
00308 template <class value_type, int DIM, int TDIM>
00309   BasisFunction<value_type,DIM,TDIM>::~BasisFunction()
00310 {}
00311 
00312 template <class value_type, int DIM, int TDIM>
00313   BasisFunction<value_type,DIM,TDIM>& BasisFunction<value_type,DIM,TDIM>::operator=(const BasisFunction<value_type,DIM,TDIM>& b)
00314 {
00315   if (&b != NULL) {
00316     ip = b.ip;
00317     id = b.id;
00318   }
00319   return *this;
00320 }
00321 
00322   template <class value_type, int DIM, int TDIM>
00323   void BasisFunction<value_type,DIM,TDIM>::reinit(const afepack::Point<TDIM>& p)
00324   {
00325     ip = p;
00326   }
00327 
00328 template <class value_type, int DIM, int TDIM>
00329     const afepack::Point<TDIM>& BasisFunction<value_type,DIM,TDIM>::interpPoint() const
00330 {
00331   return ip;
00332 }
00333 
00334 template <class value_type, int DIM, int TDIM>
00335   afepack::Point<TDIM>& BasisFunction<value_type,DIM,TDIM>::interpPoint()
00336 {
00337   return ip;
00338 }
00339 
00340 template <class value_type, int DIM, int TDIM>
00341   const typename BasisFunction<value_type,DIM,TDIM>::Identity& BasisFunction<value_type,DIM,TDIM>::identity() const
00342 {
00343   return id;
00344 }
00345 
00346 template <class value_type, int DIM, int TDIM>
00347   typename BasisFunction<value_type,DIM,TDIM>::Identity& BasisFunction<value_type,DIM,TDIM>::identity()
00348 {
00349   return id;
00350 }
00351 
00353 
00354 template <int DIM>
00355 TemplateDOF<DIM>::TemplateDOF(TemplateGeometry<DIM>& g) : 
00356 geometry(&g)
00357 {
00358   if (geometry == NULL) return;
00359         
00360   int i;
00361   n_geometry_dof.resize(DIM+1);
00362   geometry_dof.resize(DIM+1);
00363   for (i = 0;i <= DIM;i ++) {
00364     n_geometry_dof[i].resize(geometry->n_geometry(i));
00365     geometry_dof[i].resize(geometry->n_geometry(i));
00366   }
00367 }
00368 
00369 template <int DIM>
00370 TemplateDOF<DIM>::TemplateDOF(const TemplateDOF<DIM>& t) :
00371 geometry(t.geometry)
00372 {
00373   if (geometry == NULL) return;
00374         
00375   int i;
00376   n_geometry_dof.resize(DIM+1);
00377   geometry_dof.resize(DIM+1);
00378   for (i = 0;i <= DIM;i ++) {
00379     n_geometry_dof[i].resize(geometry->n_geometry(i));
00380     geometry_dof[i].resize(geometry->n_geometry(i));
00381   }
00382 }
00383 
00384 template <int DIM>
00385 TemplateDOF<DIM>::~TemplateDOF()
00386 {}
00387 
00388 template <int DIM>
00389 TemplateDOF<DIM>& TemplateDOF<DIM>::operator=(const TemplateDOF<DIM>& t)
00390 {
00391   if (&t != NULL) {
00392     n_dof = t.n_dof;
00393     n_geometry_dof = t.n_geometry_dof;
00394     geometry_dof = t.geometry_dof;
00395     dof_index = t.dof_index;
00396     geometry = t.geometry;
00397   }
00398   return *this;
00399 }
00400 
00401   template <int DIM>
00402   void TemplateDOF<DIM>::reinit(TemplateGeometry<DIM>& g)
00403   {
00404     geometry = &g;
00405     if (geometry == NULL) return;
00406         
00407     int i;
00408     n_geometry_dof.resize(DIM+1);
00409     geometry_dof.resize(DIM+1);
00410     for (i = 0;i <= DIM;i ++) {
00411       n_geometry_dof[i].resize(geometry->n_geometry(i));
00412       geometry_dof[i].resize(geometry->n_geometry(i));
00413     }
00414     dof_index.clear();
00415   }
00416 
00417 template <int DIM>
00418 void TemplateDOF<DIM>::readData(const std::string& s)
00419 {
00420   std::string library_path = FindAFEPackLibraryFilePath(s);
00421   std::string filename(library_path + "/" + s);
00422   ExpandString(filename);
00423   filtering_istream is;
00424   OpenAFEPackLibraryFile(filename, is);
00425   is >> *this;
00426 }
00427 
00428 template <int DIM>
00429 void TemplateDOF<DIM>::writeData(const std::string& s) const
00430 {
00431   std::ofstream os(s.c_str());
00432   os << *this;
00433   os.close();
00434 }
00435 
00437 
00438 template <int TDIM, int DIM>
00439   CoordTransform<TDIM,DIM>::CoordTransform() :
00440   handle(NULL)
00441 {}
00442 
00443 template <int TDIM, int DIM>
00444   CoordTransform<TDIM,DIM>::CoordTransform(const CoordTransform<TDIM,DIM>& c) :
00445   handle(NULL),
00446   library_name(c.library_name),
00447   l2g_function_name(c.l2g_function_name),
00448   g2l_function_name(c.g2l_function_name),
00449   l2g_jacobian_function_name(c.l2g_jacobian_function_name),
00450   g2l_jacobian_function_name(c.g2l_jacobian_function_name)
00451 {
00452   loadFunction();
00453 }
00454 
00455 template <int TDIM, int DIM>
00456   CoordTransform<TDIM,DIM>::~CoordTransform()
00457 {
00458   unloadFunction();
00459 }
00460 
00461 template <int TDIM, int DIM>
00462   CoordTransform<TDIM,DIM>& CoordTransform<TDIM,DIM>::operator=(const CoordTransform<TDIM,DIM>& c)
00463 {
00464   if (&c != NULL) {
00465     library_name = c.library_name;
00466     l2g_function_name = c.l2g_function_name;
00467     g2l_function_name = c.g2l_function_name;
00468     l2g_jacobian_function_name = c.l2g_jacobian_function_name;
00469     g2l_jacobian_function_name = c.g2l_jacobian_function_name;
00470   }
00471   loadFunction();
00472   return *this;
00473 }
00474 
00475   template <int TDIM, int DIM>
00476   void CoordTransform<TDIM,DIM>::loadFunction()
00477   {
00478     unloadFunction();
00479         
00480     std::string temp;
00481     if (library_path.length() == 0)
00482       temp = library_name;
00483     else
00484       temp = library_path + "/" + library_name;
00485     handle = AFEPackDLOpen(temp);
00486     if (handle == NULL) return;
00487 
00488     void * symbol = dlsym(handle, l2g_function_name.c_str());
00489     Assert(symbol, ExcLoadFunction(l2g_function_name.c_str(), library_name.c_str()));
00490     l2g_function = (void (*)(const double *, const double **, const double **, double *))symbol;
00491 
00492     symbol = dlsym(handle, g2l_function_name.c_str());
00493     Assert(symbol, ExcLoadFunction(g2l_function_name.c_str(), library_name.c_str()));
00494     g2l_function = (void (*)(const double *, const double **, const double **, double *))symbol;
00495         
00496     symbol = dlsym(handle, l2g_jacobian_function_name.c_str());
00497     Assert(symbol, ExcLoadFunction(l2g_jacobian_function_name.c_str(), library_name.c_str()));
00498     l2g_jacobian_function = (double (*)(const double *, const double **, const double **))symbol;
00499 
00500     symbol = dlsym(handle, g2l_jacobian_function_name.c_str());
00501     Assert(symbol, ExcLoadFunction(g2l_jacobian_function_name.c_str(), library_name.c_str()));
00502     g2l_jacobian_function = (double (*)(const double *, const double **, const double **))symbol;
00503   }
00504 
00505 template <int TDIM, int DIM>
00506     void CoordTransform<TDIM,DIM>::unloadFunction()
00507 {
00508   if (handle != NULL) {
00509     dlclose(handle);
00510     handle = NULL;
00511   }
00512 }
00513 
00514 template <int TDIM, int DIM>
00515   afepack::Point<DIM> CoordTransform<TDIM,DIM>::local_to_global(
00516                                                        const afepack::Point<TDIM>& lp, 
00517                                                        const std::vector<afepack::Point<TDIM> >& lv,
00518                                                        const std::vector<afepack::Point<DIM> >& gv) const
00519 {
00520   double val[DIM];
00521   const double ** v1, ** v2;
00522         
00523   int i, k;
00524   k = lv.size();
00525   typedef double * double_pointer;
00526   v1 = (const double **)new double_pointer[k];
00527   v2 = (const double **)new double_pointer[k];
00528   for (i = 0;i < k;i ++) {
00529     v1[i] = lv[i];
00530     v2[i] = gv[i];
00531   }
00532   (*l2g_function)(lp, v1, v2, val);
00533   delete[] v1;
00534   delete[] v2;
00535   return afepack::Point<DIM>(val);
00536 }
00537 
00538 template <int TDIM, int DIM>
00539   afepack::Point<TDIM> CoordTransform<TDIM,DIM>::global_to_local(
00540                                                         const afepack::Point<DIM>& gp,
00541                                                         const std::vector<afepack::Point<TDIM> >& lv,
00542                                                         const std::vector<afepack::Point<DIM> >& gv) const
00543 {
00544   double val[TDIM];
00545   const double ** v1, ** v2;
00546         
00547   int i, k;
00548   k = lv.size();
00549   typedef double * double_pointer;
00550   v1 = (const double **)new double_pointer[k];
00551   v2 = (const double **)new double_pointer[k];
00552   for (i = 0;i < k;i ++) {
00553     v1[i] = lv[i];
00554     v2[i] = gv[i];
00555   }
00556   (*g2l_function)(gp, v1, v2, val);
00557   delete[] v1;
00558   delete[] v2;
00559   return afepack::Point<TDIM>(val);
00560 }
00561 
00562 template <int TDIM, int DIM>
00563   double CoordTransform<TDIM,DIM>::local_to_global_jacobian(
00564                                                             const afepack::Point<TDIM>& lp, 
00565                                                             const std::vector<afepack::Point<TDIM> >& lv,
00566                                                             const std::vector<afepack::Point<DIM> >& gv) const
00567 {
00568   const double ** v1, ** v2;
00569         
00570   int i, k;
00571   k = lv.size();
00572   typedef double * double_pointer;
00573   v1 = (const double **)new double_pointer[k];
00574   v2 = (const double **)new double_pointer[k];
00575   for (i = 0;i < k;i ++) {
00576     v1[i] = lv[i];
00577     v2[i] = gv[i];
00578   }
00579   double val = (*l2g_jacobian_function)(lp, v1, v2);
00580   delete[] v1;
00581   delete[] v2;
00582   return val;
00583 }
00584 
00585 template <int TDIM, int DIM>
00586   double CoordTransform<TDIM,DIM>::global_to_local_jacobian(
00587                                                             const afepack::Point<DIM>& gp, 
00588                                                             const std::vector<afepack::Point<TDIM> >& lv,
00589                                                             const std::vector<afepack::Point<DIM> >& gv) const
00590 {
00591   const double ** v1, ** v2;
00592         
00593   int i, k;
00594   k = lv.size();
00595   typedef double * double_pointer;
00596   v1 = (const double **)new double_pointer[k];
00597   v2 = (const double **)new double_pointer[k];
00598   for (i = 0;i < k;i ++) {
00599     v1[i] = lv[i];
00600     v2[i] = gv[i];
00601   }
00602   double val = (*g2l_jacobian_function)(gp, v1, v2);
00603   delete[] v1;
00604   return val;
00605 }
00606 
00607 template <int TDIM, int DIM>
00608   std::vector<afepack::Point<DIM> > CoordTransform<TDIM,DIM>::local_to_global(
00609                                                                      const std::vector<afepack::Point<TDIM> >& lp, 
00610                                                                      const std::vector<afepack::Point<TDIM> >& lv,
00611                                                                      const std::vector<afepack::Point<DIM> >& gv) const
00612 {
00613   double val[DIM];
00614   const double ** v1, ** v2;
00615         
00616   int i, k;
00617   k = lv.size();
00618   typedef double * double_pointer;
00619   v1 = (const double **)new double_pointer[k];
00620   v2 = (const double **)new double_pointer[k];
00621   for (i = 0;i < k;i ++) {
00622     v1[i] = lv[i];
00623     v2[i] = gv[i];
00624   }
00625   k = lp.size();
00626   std::vector<afepack::Point<DIM> > ret(k);
00627   for (i = 0;i < k;i ++) {
00628     (*l2g_function)(lp[i], v1, v2, val);
00629     ret[i] = afepack::Point<DIM>(val);
00630   }
00631   delete[] v1;
00632   delete[] v2;
00633   return ret;
00634 }
00635 
00636 template <int TDIM, int DIM>
00637   std::vector<afepack::Point<TDIM> > CoordTransform<TDIM,DIM>::global_to_local(
00638                                                                       const std::vector<afepack::Point<DIM> >& gp,
00639                                                                       const std::vector<afepack::Point<TDIM> >& lv,
00640                                                                       const std::vector<afepack::Point<DIM> >& gv) const
00641 {
00642   double val[TDIM];
00643   const double ** v1, ** v2;
00644         
00645   int i, k;
00646   k = lv.size();
00647   typedef double * double_pointer;
00648   v1 = (const double **)new double_pointer[k];
00649   v2 = (const double **)new double_pointer[k];
00650   for (i = 0;i < k;i ++) {
00651     v1[i] = lv[i];
00652     v2[i] = gv[i];
00653   }
00654   k = gp.size();
00655   std::vector<afepack::Point<TDIM> > ret(k);
00656   for (i = 0;i < k;i ++) {
00657     (*g2l_function)(gp[i], v1, v2, val);
00658     ret[i] = afepack::Point<TDIM>(val);
00659   }
00660   delete[] v1;
00661   delete[] v2;
00662   return ret;
00663 }
00664 
00665 template <int TDIM, int DIM>
00666   std::vector<double> CoordTransform<TDIM,DIM>::local_to_global_jacobian(
00667                                                                          const std::vector<afepack::Point<TDIM> >& lp, 
00668                                                                          const std::vector<afepack::Point<TDIM> >& lv,
00669                                                                          const std::vector<afepack::Point<DIM> >& gv) const
00670 {
00671   const double ** v1, ** v2;
00672         
00673   int i, k;
00674   k = lv.size();
00675   typedef double * double_pointer;
00676   v1 = (const double **)new double_pointer[k];
00677   v2 = (const double **)new double_pointer[k];
00678   for (i = 0;i < k;i ++) {
00679     v1[i] = lv[i];
00680     v2[i] = gv[i];
00681   }
00682   k = lp.size();
00683   std::vector<double> ret(k);
00684   for (i = 0;i < k;i ++)
00685     ret[i] = (*l2g_jacobian_function)(lp[i], v1, v2);
00686   delete[] v1;
00687   delete[] v2;
00688   return ret;
00689 }
00690 
00691 template <int TDIM, int DIM>
00692   std::vector<double> CoordTransform<TDIM,DIM>::global_to_local_jacobian(
00693                                                                          const std::vector<afepack::Point<DIM> >& gp, 
00694                                                                          const std::vector<afepack::Point<TDIM> >& lv,
00695                                                                          const std::vector<afepack::Point<DIM> >& gv) const
00696 {
00697   const double ** v1, ** v2;
00698         
00699   int i, k;
00700   k = lv.size();
00701   typedef double * double_pointer;
00702   v1 = (const double **)new double_pointer[k];
00703   v2 = (const double **)new double_pointer[k];
00704   for (i = 0;i < k;i ++) {
00705     v1[i] = lv[i];
00706     v2[i] = gv[i];
00707   }
00708   k = gp.size();
00709   std::vector<double> ret(k);
00710   for (i = 0;i < k;i ++)
00711     ret[i] = (*g2l_jacobian_function)(gp[i], v1, v2);
00712   delete[] v1;
00713   delete[] v2;
00714   return ret;
00715 }
00716 
00717 template <int TDIM, int DIM>
00718   void CoordTransform<TDIM,DIM>::readData(const std::string& s)
00719 {
00720   library_path = FindAFEPackLibraryFilePath(s);
00721   std::string filename(library_path + "/" + s);
00722   ExpandString(filename);
00723   filtering_istream is;
00724   OpenAFEPackLibraryFile(filename, is);
00725   is >> *this;
00726 }
00727 
00728 template <int TDIM, int DIM>
00729   void CoordTransform<TDIM,DIM>::writeData(const std::string& s) const
00730 {
00731   std::ofstream os(s.c_str());
00732   os << *this;
00733 }
00734 
00736 
00737 template <class value_type, int DIM, int TDIM>
00738   BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin()
00739 {}
00740 
00741 
00742 template <class value_type, int DIM, int TDIM>
00743   BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const int& n) :
00744   std::vector<BasisFunction<value_type,DIM,TDIM> >(n)
00745 {}
00746 
00747 template <class value_type, int DIM, int TDIM>
00748   BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const int& n,
00749                                                               TemplateDOF<TDIM>& t) :
00750   std::vector<BasisFunction<value_type,DIM,TDIM> >(n),
00751   df(&t)
00752 {}
00753 
00754 template <class value_type, int DIM, int TDIM>
00755   BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(TemplateDOF<TDIM>& t) :
00756   df(&t)
00757 {}
00758 
00759 template <class value_type, int DIM, int TDIM>
00760   BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const BasisFunctionAdmin<value_type,DIM,TDIM>& b) :
00761   df(b.df)
00762 {}
00763 
00764 template <class value_type, int DIM, int TDIM>
00765   BasisFunctionAdmin<value_type,DIM,TDIM>::~BasisFunctionAdmin()
00766 {}
00767 
00768 template <class value_type, int DIM, int TDIM>
00769   void BasisFunctionAdmin<value_type,DIM,TDIM>::reinit(TemplateDOF<TDIM>& t)
00770 {
00771   df = &t;
00772 }
00773 
00774 template <class value_type, int DIM, int TDIM>
00775   BasisFunctionAdmin<value_type,DIM,TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::operator=(const BasisFunctionAdmin<value_type,DIM,TDIM>& b)
00776 {
00777   df = b.df;
00778   return *this;
00779 }
00780 
00781   template <class value_type, int DIM, int TDIM>
00782   const TemplateDOF<TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::dof() const
00783   {
00784     return *df;
00785   }
00786 
00787 template <class value_type, int DIM, int TDIM>
00788     TemplateDOF<TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::dof()
00789 {
00790   return *df;
00791 }
00792 
00793 template <class value_type, int DIM, int TDIM>
00794   void BasisFunctionAdmin<value_type,DIM,TDIM>::readData(const std::string& s)
00795 {
00796   library_path = FindAFEPackLibraryFilePath(s);
00797   std::string filename(library_path + "/" + s);
00798   ExpandString(filename);
00799   filtering_istream is;
00800   OpenAFEPackLibraryFile(filename, is);
00801   is >> *this;
00802 }
00803 
00804 template <class value_type, int DIM, int TDIM>
00805   void BasisFunctionAdmin<value_type,DIM,TDIM>::writeData(const std::string& s) const
00806 {
00807   std::ofstream os(s.c_str());
00808   os << *this;
00809   os.close();
00810 }
00811 
00813 
00814 template <int DIM>
00815 UnitOutNormal<DIM>::UnitOutNormal() :
00816 handle(NULL)
00817 {}
00818 
00819 template <int DIM>
00820 UnitOutNormal<DIM>::UnitOutNormal(const UnitOutNormal<DIM>& c) :
00821 handle(NULL),
00822   library_name(c.library_name),
00823   function_name(c.function_name)
00824 {
00825   loadFunction();
00826 }
00827 
00828 template <int DIM>
00829 UnitOutNormal<DIM>::~UnitOutNormal()
00830 {
00831   unloadFunction();
00832 }
00833 
00834 template <int DIM>
00835 UnitOutNormal<DIM>& UnitOutNormal<DIM>::operator=(const UnitOutNormal<DIM>& c)
00836 {
00837   if (&c != NULL) {
00838     library_name = c.library_name;
00839     function_name = c.function_name;
00840   }
00841   return *this;
00842 }
00843 
00844   template <int DIM>
00845   void UnitOutNormal<DIM>::loadFunction()
00846   {
00847     unloadFunction();
00848         
00849     std::string temp;
00850     if (library_path.length() == 0)
00851       temp = library_name;
00852     else
00853       temp = library_path + "/" + library_name;
00854     handle = AFEPackDLOpen(temp);
00855     if (handle == NULL) return;
00856 
00857     void * symbol = dlsym(handle, function_name.c_str());
00858     Assert(symbol, ExcLoadFunction(function_name.c_str(), library_name.c_str()));
00859     function = (void (*)(const double *, const double **, int, double *))symbol;
00860   }
00861 
00862 template <int DIM>
00863 void UnitOutNormal<DIM>::unloadFunction()
00864 {
00865   if (handle != NULL) {
00866     dlclose(handle);
00867     handle = NULL;
00868   }
00869 }
00870 
00871 
00872 template <int DIM>
00873 std::vector<double> UnitOutNormal<DIM>::value(const afepack::Point<DIM>& p,
00874                                               const std::vector<afepack::Point<DIM> >& v, const int& n) const
00875 {
00876   double val[DIM];
00877   const double ** v1;
00878         
00879   int i, k;
00880   k = v.size();
00881   typedef double * double_pointer;
00882   v1 = (const double **)new double_pointer[k];
00883   for (i = 0;i < k;i ++) v1[i] = v[i];
00884   (*function)(p, v1, n, val);
00885   delete[] v1;
00886   return std::vector<double>(&val[0], &val[DIM]);
00887 }
00888 
00889 template <int DIM>
00890 std::vector<std::vector<double> > UnitOutNormal<DIM>::value(const std::vector<afepack::Point<DIM> >& p,
00891                                                             const std::vector<afepack::Point<DIM> >& v, const int& n) const
00892 {
00893   double val[DIM];
00894   const double ** v1;
00895         
00896   int i, k;
00897   k = v.size();
00898   typedef double * double_pointer;
00899   v1 = (const double **)new double_pointer[k];
00900   for (i = 0;i < k;i ++) v1[i] = v[i];
00901   k = p.size();
00902   std::vector<std::vector<double> > ret(k, std::vector<double>(DIM));
00903   for (i = 0;i < k;i ++) {
00904     (*function)(p[i], v1, n, val);
00905     std::copy(&val[0], &val[0]+DIM, ret[i].begin());
00906   }
00907   delete[] v1;
00908   return ret;
00909 }
00910 
00911 template <int DIM>
00912 std::vector<double> UnitOutNormal<DIM>::value(const afepack::Point<DIM>& p,
00913                                               const double ** v1, const int& n) const
00914 {
00915   double val[DIM];
00916   (*function)(p, v1, n, val);
00917   std::vector<double> ret(DIM);
00918   return std::vector<double>(&val[0], &val[DIM]);
00919 }
00920 
00921 template <int DIM>
00922 std::vector<std::vector<double> > UnitOutNormal<DIM>::value(const std::vector<afepack::Point<DIM> >& p,
00923                                                             const double ** v1, const int& n) const
00924 {
00925   double val[DIM];
00926         
00927   int i, k;
00928   k = p.size();
00929   std::vector<std::vector<double> > ret(k, std::vector<double>(DIM, 0.0));
00930   for (i = 0;i < k;i ++) {
00931     (*function)(p[i], v1, n, val);
00932     std::copy(&val[0], &val[0]+DIM, ret[i].begin());
00933   }
00934   return ret;
00935 }
00936 
00937 template <int DIM>
00938 void UnitOutNormal<DIM>::readData(const std::string& s)
00939 {
00940   library_path = FindAFEPackLibraryFilePath(s);
00941   std::string filename(library_path + "/" + s);
00942   ExpandString(filename);
00943   filtering_istream is;
00944   OpenAFEPackLibraryFile(filename, is);
00945   is >> *this;
00946 }
00947 
00948 template <int DIM>
00949 void UnitOutNormal<DIM>::writeData(const std::string& s) const
00950 {
00951   std::ofstream os(s.c_str());
00952   os << *this;
00953 }
00954 
00956 
00957 template <class value_type, int DIM, int TDIM>
00958   TemplateElement<value_type,DIM,TDIM>::TemplateElement(
00959                                                         TemplateGeometry<TDIM>& g,
00960                                                         TemplateDOF<TDIM>& d,
00961                                                         CoordTransform<TDIM,DIM>& c,
00962                                                         BasisFunctionAdmin<value_type,DIM,TDIM>& b,
00963                                                         UnitOutNormal<DIM>& u) :
00964   geo(&g),
00965   df(&d),
00966   ct(&c),
00967   bf(&b),
00968   uon(&u)
00969 {}
00970 
00971 template <class value_type, int DIM, int TDIM>
00972   TemplateElement<value_type,DIM,TDIM>::TemplateElement(const TemplateElement<value_type,DIM,TDIM>& t) :
00973   geo(t.geo),
00974   df(t.df),
00975   ct(t.ct),
00976   bf(t.bf)
00977 {}
00978 
00979 template <class value_type, int DIM, int TDIM>
00980   TemplateElement<value_type,DIM,TDIM>::~TemplateElement()
00981 {}
00982 
00983 template <class value_type, int DIM, int TDIM>
00984   TemplateElement<value_type,DIM,TDIM>& 
00985   TemplateElement<value_type,DIM,TDIM>::operator=(const TemplateElement<value_type,DIM,TDIM>& t)
00986 {
00987   if (&t != NULL) {
00988     geo = t.geo;
00989     df = t.df;
00990     ct = t.ct;
00991     bf = t.bf;
00992   }
00993   return *this;
00994 }
00995 
00996   template <class value_type, int DIM, int TDIM>
00997   void TemplateElement<value_type,DIM,TDIM>::reinit(
00998                                                     TemplateGeometry<TDIM>& g,
00999                                                     TemplateDOF<TDIM>& d,
01000                                                     CoordTransform<TDIM,DIM>& c,
01001                                                     BasisFunctionAdmin<value_type,DIM,TDIM>& b,
01002                                                     UnitOutNormal<DIM>& u)
01003   {
01004     geo = &g;
01005     df = &d;
01006     ct = &c;
01007     bf = &b;
01008     uon = &u;
01009   }
01010 
01011 template <class value_type, int DIM, int TDIM>
01012     const std::vector<afepack::Point<TDIM> >& TemplateElement<value_type,DIM,TDIM>::vertexArray() const
01013 {
01014   return geo->vertexArray();
01015 }
01016 
01017 #endif //_TemplateElement_templates_h_
01018