|
AFEPack
|
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 Point<DIM>& p, 00103 const std::vector<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 Point<DIM>& p, 00115 const std::vector<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<Point<DIM> >& p, 00128 const std::vector<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<Point<DIM> >& p, 00144 const std::vector<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 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 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<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<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 Point<DIM>& p, 00203 const std::vector<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<Point<DIM> >& p, 00229 const std::vector<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 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<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 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 Point<TDIM>& p) 00324 { 00325 ip = p; 00326 } 00327 00328 template <class value_type, int DIM, int TDIM> 00329 const 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 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 Point<DIM> CoordTransform<TDIM,DIM>::local_to_global( 00516 const Point<TDIM>& lp, 00517 const std::vector<Point<TDIM> >& lv, 00518 const std::vector<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 Point<DIM>(val); 00536 } 00537 00538 template <int TDIM, int DIM> 00539 Point<TDIM> CoordTransform<TDIM,DIM>::global_to_local( 00540 const Point<DIM>& gp, 00541 const std::vector<Point<TDIM> >& lv, 00542 const std::vector<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 Point<TDIM>(val); 00560 } 00561 00562 template <int TDIM, int DIM> 00563 double CoordTransform<TDIM,DIM>::local_to_global_jacobian( 00564 const Point<TDIM>& lp, 00565 const std::vector<Point<TDIM> >& lv, 00566 const std::vector<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 Point<DIM>& gp, 00588 const std::vector<Point<TDIM> >& lv, 00589 const std::vector<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<Point<DIM> > CoordTransform<TDIM,DIM>::local_to_global( 00609 const std::vector<Point<TDIM> >& lp, 00610 const std::vector<Point<TDIM> >& lv, 00611 const std::vector<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<Point<DIM> > ret(k); 00627 for (i = 0;i < k;i ++) { 00628 (*l2g_function)(lp[i], v1, v2, val); 00629 ret[i] = Point<DIM>(val); 00630 } 00631 delete[] v1; 00632 delete[] v2; 00633 return ret; 00634 } 00635 00636 template <int TDIM, int DIM> 00637 std::vector<Point<TDIM> > CoordTransform<TDIM,DIM>::global_to_local( 00638 const std::vector<Point<DIM> >& gp, 00639 const std::vector<Point<TDIM> >& lv, 00640 const std::vector<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<Point<TDIM> > ret(k); 00656 for (i = 0;i < k;i ++) { 00657 (*g2l_function)(gp[i], v1, v2, val); 00658 ret[i] = 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<Point<TDIM> >& lp, 00668 const std::vector<Point<TDIM> >& lv, 00669 const std::vector<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<Point<DIM> >& gp, 00694 const std::vector<Point<TDIM> >& lv, 00695 const std::vector<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 return ret; 00714 } 00715 00716 template <int TDIM, int DIM> 00717 void CoordTransform<TDIM,DIM>::readData(const std::string& s) 00718 { 00719 library_path = FindAFEPackLibraryFilePath(s); 00720 std::string filename(library_path + "/" + s); 00721 ExpandString(filename); 00722 filtering_istream is; 00723 OpenAFEPackLibraryFile(filename, is); 00724 is >> *this; 00725 } 00726 00727 template <int TDIM, int DIM> 00728 void CoordTransform<TDIM,DIM>::writeData(const std::string& s) const 00729 { 00730 std::ofstream os(s.c_str()); 00731 os << *this; 00732 } 00733 00735 00736 template <class value_type, int DIM, int TDIM> 00737 BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin() 00738 {} 00739 00740 00741 template <class value_type, int DIM, int TDIM> 00742 BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const int& n) : 00743 std::vector<BasisFunction<value_type,DIM,TDIM> >(n) 00744 {} 00745 00746 template <class value_type, int DIM, int TDIM> 00747 BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const int& n, 00748 TemplateDOF<TDIM>& t) : 00749 std::vector<BasisFunction<value_type,DIM,TDIM> >(n), 00750 df(&t) 00751 {} 00752 00753 template <class value_type, int DIM, int TDIM> 00754 BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(TemplateDOF<TDIM>& t) : 00755 df(&t) 00756 {} 00757 00758 template <class value_type, int DIM, int TDIM> 00759 BasisFunctionAdmin<value_type,DIM,TDIM>::BasisFunctionAdmin(const BasisFunctionAdmin<value_type,DIM,TDIM>& b) : 00760 df(b.df) 00761 {} 00762 00763 template <class value_type, int DIM, int TDIM> 00764 BasisFunctionAdmin<value_type,DIM,TDIM>::~BasisFunctionAdmin() 00765 {} 00766 00767 template <class value_type, int DIM, int TDIM> 00768 void BasisFunctionAdmin<value_type,DIM,TDIM>::reinit(TemplateDOF<TDIM>& t) 00769 { 00770 df = &t; 00771 } 00772 00773 template <class value_type, int DIM, int TDIM> 00774 BasisFunctionAdmin<value_type,DIM,TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::operator=(const BasisFunctionAdmin<value_type,DIM,TDIM>& b) 00775 { 00776 df = b.df; 00777 return *this; 00778 } 00779 00780 template <class value_type, int DIM, int TDIM> 00781 const TemplateDOF<TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::dof() const 00782 { 00783 return *df; 00784 } 00785 00786 template <class value_type, int DIM, int TDIM> 00787 TemplateDOF<TDIM>& BasisFunctionAdmin<value_type,DIM,TDIM>::dof() 00788 { 00789 return *df; 00790 } 00791 00792 template <class value_type, int DIM, int TDIM> 00793 void BasisFunctionAdmin<value_type,DIM,TDIM>::readData(const std::string& s) 00794 { 00795 library_path = FindAFEPackLibraryFilePath(s); 00796 std::string filename(library_path + "/" + s); 00797 ExpandString(filename); 00798 filtering_istream is; 00799 OpenAFEPackLibraryFile(filename, is); 00800 is >> *this; 00801 } 00802 00803 template <class value_type, int DIM, int TDIM> 00804 void BasisFunctionAdmin<value_type,DIM,TDIM>::writeData(const std::string& s) const 00805 { 00806 std::ofstream os(s.c_str()); 00807 os << *this; 00808 os.close(); 00809 } 00810 00812 00813 template <int DIM> 00814 UnitOutNormal<DIM>::UnitOutNormal() : 00815 handle(NULL) 00816 {} 00817 00818 template <int DIM> 00819 UnitOutNormal<DIM>::UnitOutNormal(const UnitOutNormal<DIM>& c) : 00820 handle(NULL), 00821 library_name(c.library_name), 00822 function_name(c.function_name) 00823 { 00824 loadFunction(); 00825 } 00826 00827 template <int DIM> 00828 UnitOutNormal<DIM>::~UnitOutNormal() 00829 { 00830 unloadFunction(); 00831 } 00832 00833 template <int DIM> 00834 UnitOutNormal<DIM>& UnitOutNormal<DIM>::operator=(const UnitOutNormal<DIM>& c) 00835 { 00836 if (&c != NULL) { 00837 library_name = c.library_name; 00838 function_name = c.function_name; 00839 } 00840 return *this; 00841 } 00842 00843 template <int DIM> 00844 void UnitOutNormal<DIM>::loadFunction() 00845 { 00846 unloadFunction(); 00847 00848 std::string temp; 00849 if (library_path.length() == 0) 00850 temp = library_name; 00851 else 00852 temp = library_path + "/" + library_name; 00853 handle = AFEPackDLOpen(temp); 00854 if (handle == NULL) return; 00855 00856 void * symbol = dlsym(handle, function_name.c_str()); 00857 Assert(symbol, ExcLoadFunction(function_name.c_str(), library_name.c_str())); 00858 function = (void (*)(const double *, const double **, int, double *))symbol; 00859 } 00860 00861 template <int DIM> 00862 void UnitOutNormal<DIM>::unloadFunction() 00863 { 00864 if (handle != NULL) { 00865 dlclose(handle); 00866 handle = NULL; 00867 } 00868 } 00869 00870 00871 template <int DIM> 00872 std::vector<double> UnitOutNormal<DIM>::value(const Point<DIM>& p, 00873 const std::vector<Point<DIM> >& v, const int& n) const 00874 { 00875 double val[DIM]; 00876 const double ** v1; 00877 00878 int i, k; 00879 k = v.size(); 00880 typedef double * double_pointer; 00881 v1 = (const double **)new double_pointer[k]; 00882 for (i = 0;i < k;i ++) v1[i] = v[i]; 00883 (*function)(p, v1, n, val); 00884 delete[] v1; 00885 return std::vector<double>(&val[0], &val[DIM]); 00886 } 00887 00888 template <int DIM> 00889 std::vector<std::vector<double> > UnitOutNormal<DIM>::value(const std::vector<Point<DIM> >& p, 00890 const std::vector<Point<DIM> >& v, const int& n) const 00891 { 00892 double val[DIM]; 00893 const double ** v1; 00894 00895 int i, k; 00896 k = v.size(); 00897 typedef double * double_pointer; 00898 v1 = (const double **)new double_pointer[k]; 00899 for (i = 0;i < k;i ++) v1[i] = v[i]; 00900 k = p.size(); 00901 std::vector<std::vector<double> > ret(k, std::vector<double>(DIM)); 00902 for (i = 0;i < k;i ++) { 00903 (*function)(p[i], v1, n, val); 00904 std::copy(&val[0], &val[0]+DIM, ret[i].begin()); 00905 } 00906 delete[] v1; 00907 return ret; 00908 } 00909 00910 template <int DIM> 00911 std::vector<double> UnitOutNormal<DIM>::value(const Point<DIM>& p, 00912 const double ** v1, const int& n) const 00913 { 00914 double val[DIM]; 00915 (*function)(p, v1, n, val); 00916 std::vector<double> ret(DIM); 00917 return std::vector<double>(&val[0], &val[DIM]); 00918 } 00919 00920 template <int DIM> 00921 std::vector<std::vector<double> > UnitOutNormal<DIM>::value(const std::vector<Point<DIM> >& p, 00922 const double ** v1, const int& n) const 00923 { 00924 double val[DIM]; 00925 00926 int i, k; 00927 k = p.size(); 00928 std::vector<std::vector<double> > ret(k, std::vector<double>(DIM, 0.0)); 00929 for (i = 0;i < k;i ++) { 00930 (*function)(p[i], v1, n, val); 00931 std::copy(&val[0], &val[0]+DIM, ret[i].begin()); 00932 } 00933 return ret; 00934 } 00935 00936 template <int DIM> 00937 void UnitOutNormal<DIM>::readData(const std::string& s) 00938 { 00939 library_path = FindAFEPackLibraryFilePath(s); 00940 std::string filename(library_path + "/" + s); 00941 ExpandString(filename); 00942 filtering_istream is; 00943 OpenAFEPackLibraryFile(filename, is); 00944 is >> *this; 00945 } 00946 00947 template <int DIM> 00948 void UnitOutNormal<DIM>::writeData(const std::string& s) const 00949 { 00950 std::ofstream os(s.c_str()); 00951 os << *this; 00952 } 00953 00955 00956 template <class value_type, int DIM, int TDIM> 00957 TemplateElement<value_type,DIM,TDIM>::TemplateElement( 00958 TemplateGeometry<TDIM>& g, 00959 TemplateDOF<TDIM>& d, 00960 CoordTransform<TDIM,DIM>& c, 00961 BasisFunctionAdmin<value_type,DIM,TDIM>& b, 00962 UnitOutNormal<DIM>& u) : 00963 geo(&g), 00964 df(&d), 00965 ct(&c), 00966 bf(&b), 00967 uon(&u) 00968 {} 00969 00970 template <class value_type, int DIM, int TDIM> 00971 TemplateElement<value_type,DIM,TDIM>::TemplateElement(const TemplateElement<value_type,DIM,TDIM>& t) : 00972 geo(t.geo), 00973 df(t.df), 00974 ct(t.ct), 00975 bf(t.bf) 00976 {} 00977 00978 template <class value_type, int DIM, int TDIM> 00979 TemplateElement<value_type,DIM,TDIM>::~TemplateElement() 00980 {} 00981 00982 template <class value_type, int DIM, int TDIM> 00983 TemplateElement<value_type,DIM,TDIM>& 00984 TemplateElement<value_type,DIM,TDIM>::operator=(const TemplateElement<value_type,DIM,TDIM>& t) 00985 { 00986 if (&t != NULL) { 00987 geo = t.geo; 00988 df = t.df; 00989 ct = t.ct; 00990 bf = t.bf; 00991 } 00992 return *this; 00993 } 00994 00995 template <class value_type, int DIM, int TDIM> 00996 void TemplateElement<value_type,DIM,TDIM>::reinit( 00997 TemplateGeometry<TDIM>& g, 00998 TemplateDOF<TDIM>& d, 00999 CoordTransform<TDIM,DIM>& c, 01000 BasisFunctionAdmin<value_type,DIM,TDIM>& b, 01001 UnitOutNormal<DIM>& u) 01002 { 01003 geo = &g; 01004 df = &d; 01005 ct = &c; 01006 bf = &b; 01007 uon = &u; 01008 } 01009 01010 template <class value_type, int DIM, int TDIM> 01011 const std::vector<Point<TDIM> >& TemplateElement<value_type,DIM,TDIM>::vertexArray() const 01012 { 01013 return geo->vertexArray(); 01014 } 01015 01016 #endif //_TemplateElement_templates_h_ 01017
1.7.4