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