AFEPack
|
00001 00011 #ifndef _FEMSpace_h_ 00012 #define _FEMSpace_h_ 00013 00014 #include <iostream> 00015 #include <string> 00016 #include <vector> 00017 #include <iterator> 00018 #include <algorithm> 00019 00020 #include <base/exceptions.h> 00021 #include <lac/vector.h> 00022 #include <lac/sparse_matrix.h> 00023 00024 #include "Geometry.h" 00025 #include "HGeometry.h" 00026 #include "TemplateElement.h" 00027 #ifdef MULTITHREAD 00028 #include "Thread.h" 00029 #endif 00030 00031 template <class value_type, int DIM, int DOW, int TDIM> class Element; 00032 template <class value_type, int DIM, int DOW, int TDIM> class FEMSpace; 00033 template <class value_type, int DIM, int DOW, int TDIM, typename Number> class FEMFunction; 00034 00036 00043 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM> 00044 class Element 00045 { 00046 public: 00047 enum { dim = DIM, dow = DOW, tdim = TDIM }; 00048 typedef value_type value_t; 00049 typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t; 00050 typedef Element<value_t,DIM,DOW,TDIM> element_t; 00051 typedef TemplateElement<value_t,DOW,TDIM> template_t; 00052 00053 private: 00054 fe_space_t * sp; 00055 int geometry_index; 00056 int template_element_index; 00057 std::vector<int> dof_index; 00058 std::vector<std::vector<int> > geo_img; 00059 public: 00060 explicit Element(fe_space_t& = *((fe_space_t *)NULL)); 00061 Element(const element_t&); 00062 ~Element(); 00063 public: 00064 void reinit(fe_space_t&, int, int); 00065 void reinit(int, int); 00066 void reinit(fe_space_t&, int, int, const std::vector<int>&); 00067 element_t& operator=(const element_t&); 00068 const fe_space_t& femSpace() const; 00069 fe_space_t& femSpace(); 00070 Mesh<DIM,DOW>& mesh() const; 00071 int index() const; 00072 int& index(); 00073 GeometryBM& geometry() const; 00074 void geometry(const GeometryBM&); 00075 void geometry(int); 00076 int templateElementIndex() const; 00077 template_t& templateElement() const; 00078 void templateElement(int); 00079 TemplateGeometry<TDIM>& templateGeometry() const; 00080 const std::vector<int>& dof() const; 00081 std::vector<int>& dof(); 00082 unsigned int n_dof() const; 00083 const std::vector<std::vector<int> >& geometryImage() const; 00084 std::vector<std::vector<int> >& geometryImage(); 00085 void buildGeometryImage(); 00086 void lazyBuildGeometryImage(); 00087 void buildVertexArray(std::vector<afepack::Point<DOW> >&) const; 00088 const double ** buildVertexArray() const; 00090 const typename BasisFunction<value_t,DIM,TDIM>::Identity& basis_function_identity(int i) const; 00091 value_t basis_function_value(int i, const afepack::Point<DOW>&) const; 00092 std::vector<value_t> basis_function_value(const afepack::Point<DOW>&) const; 00093 std::vector<value_t> basis_function_value(int i, const std::vector<afepack::Point<DOW> >&) const; 00094 std::vector<std::vector<value_t> > basis_function_value(const std::vector<afepack::Point<DOW> >&) const; 00095 std::vector<value_t> basis_function_gradient(int i, const afepack::Point<DOW>&) const; 00096 std::vector<std::vector<value_t> > basis_function_gradient(const afepack::Point<DOW>&) const; 00097 std::vector<std::vector<value_t> > basis_function_gradient(int i, const std::vector<afepack::Point<DOW> >&) const; 00098 std::vector<std::vector<std::vector<value_t> > > basis_function_gradient(const std::vector<afepack::Point<DOW> >&) const; 00100 #ifdef QUADRATIC_ELEMENT_SUPPORT 00101 std::vector<std::vector<value_t> > 00102 basis_function_hesse(int i, 00103 const afepack::Point<DOW>&) const; 00104 std::vector<std::vector<std::vector<value_t> > > 00105 basis_function_hesse(const afepack::Point<DOW>&) const; 00106 std::vector<std::vector<std::vector<value_t> > > 00107 basis_function_hesse(int i, 00108 const std::vector<afepack::Point<DOW> >&) const; 00109 std::vector<std::vector<std::vector<std::vector<value_t> > > > 00110 basis_function_gradient(const std::vector<afepack::Point<DOW> >&) const; 00111 #endif // QUADRATIC_ELEMENT_SUPPORT 00112 00113 afepack::Point<DOW> local_to_global(const afepack::Point<TDIM>&) const; 00114 afepack::Point<TDIM> global_to_local(const afepack::Point<DOW>&) const; 00115 double local_to_global_jacobian(const afepack::Point<TDIM>&) const; 00116 double global_to_local_jacobian(const afepack::Point<DOW>&) const; 00117 std::vector<afepack::Point<DOW> > local_to_global(const std::vector<afepack::Point<TDIM> >&) const; 00118 std::vector<afepack::Point<TDIM> > global_to_local(const std::vector<afepack::Point<DOW> >&) const; 00119 std::vector<double> local_to_global_jacobian(const std::vector<afepack::Point<TDIM> >&) const; 00120 std::vector<double> global_to_local_jacobian(const std::vector<afepack::Point<DOW> >&) const; 00122 std::vector<double> unitOutNormal(const afepack::Point<DOW>&, int) const; 00123 std::vector<std::vector<double> > unitOutNormal(const std::vector<afepack::Point<DOW> >&, int) const; 00124 00125 const QuadratureInfo<TDIM>& findQuadratureInfo(int) const; 00127 public: 00128 DeclException1(ExcMeshData, std::string, << "Mesh data uncompatible: " << arg1); 00129 #ifdef __SERIALIZATION__ 00130 template <class T> T * 00131 new_property(const property_id_t<T>& pid) const { 00132 return femSpace().new_property(*this, pid); 00133 } 00134 template <class T> T * 00135 get_property(const property_id_t<T>& pid) const { 00136 return femSpace().get_property(*this, pid); 00137 } 00138 template <class T> void 00139 free_property(const property_id_t<T>& pid) const { 00140 femSpace().free_property(*this, pid); 00141 } 00142 #endif // __SERIALIZATION__ 00143 }; 00144 00153 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM> 00154 class FEMSpace 00155 { 00156 public: 00157 enum { dim = DIM, dow = DOW, tdim = TDIM }; 00158 00159 typedef value_type value_t; 00160 typedef Element<value_t,DIM,DOW,TDIM> element_t; 00161 typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t; 00162 typedef DOFInfo<value_t, DIM, DOW,TDIM> dof_info_t; 00163 typedef typename element_t::template_t template_t; 00164 typedef Mesh<DIM,DOW> mesh_t; 00165 typedef typename mesh_t::bmark_t bmark_t; 00166 typedef typename std::vector<element_t>::iterator ElementIterator; 00167 typedef typename std::vector<element_t>::const_iterator ConstElementIterator; 00168 00169 private: 00170 mesh_t * msh; 00171 std::vector<template_t> * tmp_ele; 00172 std::vector<element_t> ele; 00173 DegreeOfFreedom df; 00174 std::vector<dof_info_t> df_in; 00176 u_int effi_flag; 00178 public: 00179 explicit FEMSpace(mesh_t& = *((mesh_t *)NULL), 00180 std::vector<template_t>& 00181 = *((std::vector<template_t> *)NULL)); 00182 FEMSpace(const fe_space_t&); 00183 virtual ~FEMSpace(); 00184 public: 00185 fe_space_t& operator=(const fe_space_t&); 00186 void reinit(mesh_t& = *((mesh_t *)NULL), 00187 std::vector<template_t>& 00188 = *((std::vector<template_t> *)NULL)); 00190 const mesh_t& mesh() const {return *msh;} 00191 mesh_t& mesh() {return *msh;} 00192 00193 const std::vector<template_t>& templateElement() const {return *tmp_ele;} 00194 std::vector<template_t>& templateElement() {return *tmp_ele;} 00195 const template_t& templateElement(int i) const {return (*tmp_ele)[i];} 00196 template_t& templateElement(int i) {return (*tmp_ele)[i];} 00197 00198 int n_element() {return ele.size();}; 00199 const std::vector<element_t>& element() const {return ele;} 00200 std::vector<element_t>& element() {return ele;} 00201 const element_t& element(int i) const {return ele[i];} 00202 element_t& element(int i) {return ele[i];} 00203 00204 const DegreeOfFreedom& dof() const {return df;} 00205 DegreeOfFreedom& dof() {return df;} 00206 00207 const std::vector<dof_info_t>& dofInfo() const {return df_in;} 00208 std::vector<dof_info_t>& dofInfo() {return df_in;} 00209 const dof_info_t& dofInfo(int i) const {return df_in[i];} 00210 dof_info_t& dofInfo(int i) {return df_in[i];} 00211 00212 unsigned int n_dof() const {return dof().n_dof;} 00213 const std::vector<int>& geometryDof(int i, int j) const {return dof().geometry_dof[i][j];} 00214 const std::vector<int>& geometryDof(int i, const Geometry& g) const {return dof().geometry_dof[i][g.index()];} 00215 00216 const std::vector<DOFIndex>& dofIndex() const {return dof().dof_index;} 00217 const DOFIndex& dofIndex(int i) const {return dof().dof_index[i];} 00218 00219 const bmark_t& dofBoundaryMark(int i) const {return dofInfo(i).boundary_mark;} 00220 bmark_t& dofBoundaryMark(int i) {return dofInfo(i).boundary_mark;} 00221 00222 virtual void buildElement(bool lazy = false); 00224 virtual void buildDofBoundaryMark(); 00226 void buildDof(); 00227 void updateDofInterpPoint(); 00229 ElementIterator beginElement() {return element().begin();}; 00230 ConstElementIterator beginElement() const {return element().begin();}; 00232 ElementIterator endElement() {return element().end();}; 00233 ConstElementIterator endElement() const {return element().end();}; 00234 00235 DeclException1(ExcDOFData, std::string, << "DOF data uncompatible: " << arg1); 00236 00237 const unsigned int& efficiencyFlag() const {return effi_flag;}; 00238 unsigned int& efficiencyFlag() {return effi_flag;} 00239 00240 public: 00241 #ifdef __SERIALIZATION__ 00242 template <class T> T * 00243 new_property(const element_t& e, 00244 const property_id_t<T>& pid) const { 00245 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00246 const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh); 00247 return mesh->template new_property<T,DIM>(e.geometry().index(), pid); 00248 } 00249 template <class T> T * 00250 get_property(const element_t& e, 00251 const property_id_t<T>& pid) const { 00252 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00253 const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh); 00254 return mesh->template get_property<T,DIM>(e.geometry().index(), pid); 00255 } 00256 template <class T> void 00257 free_property(const element_t& e, 00258 const property_id_t<T>& pid) const { 00259 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00260 const reg_mesh_t * mesh = dynamic_cast<const reg_mesh_t *>(msh); 00261 mesh->template free_property<T,DIM>(e.geometry().index(), pid); 00262 } 00263 #endif // __SERIALIZATION__ 00264 00265 friend class Element<value_t,DIM,DOW,TDIM>; 00266 00267 #ifdef MULTITHREAD 00268 void threadBuildElement(bool lazy, int n_thread = 1, int rank = 0); 00269 void threadBuildDof0(std::vector<std::vector<bool> >&, 00270 pthread_mutex_t&, int n_thread = 1, int rank = 0); 00271 void threadBuildDof1(std::vector<std::vector<bool> >&, 00272 pthread_mutex_t&, int n_thread = 1, int rank = 0); 00273 #endif 00274 }; 00275 00281 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double> 00282 class FEMFunction : public Vector<Number> 00283 { 00284 public: 00285 enum { dim = DIM, dow = DOW, tdim = TDIM }; 00286 00287 typedef value_type value_t; 00288 typedef Number number_t; 00289 typedef Element<value_t,DIM,DOW,TDIM> element_t; 00290 typedef FEMSpace<value_t,DIM,DOW,TDIM> fe_space_t; 00291 typedef FEMFunction<value_t,DIM,DOW,TDIM,number_t> fe_func_t; 00292 00293 private: 00294 fe_space_t * sp; 00295 public: 00296 FEMFunction(fe_space_t & = *((fe_space_t *)NULL)); 00297 virtual ~FEMFunction(); 00298 public: 00299 fe_space_t& femSpace(); 00300 const fe_space_t& femSpace() const { return *sp; }; 00301 void reinit(fe_space_t&, bool = false); 00303 value_type value(const afepack::Point<DOW>&, const element_t&) const; 00304 std::vector<value_type> gradient(const afepack::Point<DOW>&, const element_t&) const; 00305 std::vector<value_type> value(const std::vector<afepack::Point<DOW> >&, const element_t&) const; 00306 std::vector<std::vector<value_type> > gradient(const std::vector<afepack::Point<DOW> >&, const element_t&) const; 00313 value_type value(const std::vector<value_type>&, const element_t&) const; 00314 std::vector<value_type> gradient(const std::vector<std::vector<value_type> >&, const element_t&) const; 00315 std::vector<value_type> value(const std::vector<std::vector<value_type> >&, const element_t&) const; 00316 std::vector<std::vector<value_type> > gradient(const std::vector<std::vector<std::vector<value_type> > >&, const element_t&) const; 00318 public: 00319 void loadData(const std::string& filename); 00320 void writeData(const std::string& filename); 00322 void writeEasyMeshData(const std::string& filename); 00332 void writeTecplotData(const std::string& filename); 00344 void writeOpenDXData(const std::string& filename, 00345 int flag = 0) const; 00346 }; 00347 00348 00349 00350 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double> 00351 class LocalFEMFunction : public Vector<Number> 00352 { 00353 private: 00354 Element<value_type,DIM,DOW,TDIM> * ele; 00355 public: 00356 LocalFEMFunction(Element<value_type,DIM,DOW,TDIM> & = *((Element<value_type,DIM,DOW,TDIM> *)NULL)); 00357 virtual ~LocalFEMFunction() {}; 00358 public: 00359 Element<value_type,DIM,DOW,TDIM>& element(); 00360 void reinit(Element<value_type,DIM,DOW,TDIM> &); 00361 value_type value(const afepack::Point<DOW>&) const; 00362 std::vector<value_type> gradient(const afepack::Point<DOW>&) const; 00363 std::vector<value_type> value(const std::vector<afepack::Point<DOW> >&) const; 00364 std::vector<std::vector<value_type> > gradient(const std::vector<afepack::Point<DOW> >&) const; 00365 }; 00366 00367 00368 00376 struct BoundaryConditionInfo 00377 { 00378 public: 00379 typedef int Type; 00380 static Type DIRICHLET; 00381 static Type NEUMANN; 00382 static Type ROBIN; 00383 }; 00384 00385 00386 00387 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double> 00388 class BoundaryCondition : public BoundaryConditionInfo 00389 { 00390 public: 00391 typedef typename FEMSpace<value_type,DIM,DOW,TDIM>::bmark_t bmark_t; 00392 typedef typename BoundaryConditionInfo::Type Type; 00393 private: 00394 Type boundary_type; 00395 bmark_t boundary_mark; 00396 public: 00397 BoundaryCondition() {}; 00398 BoundaryCondition(const Type& t, const bmark_t& m) : 00399 boundary_type(t), boundary_mark(m) {}; 00400 BoundaryCondition(const BoundaryCondition<value_type,DIM>& b) : 00401 boundary_type(b.boundaryType()), boundary_mark(b.boundaryMark()) {}; 00402 virtual ~BoundaryCondition() {}; 00403 public: 00404 const Type& boundaryType() const {return boundary_type;}; 00405 Type& boundaryType() {return boundary_type;}; 00406 const bmark_t& boundaryMark() const {return boundary_mark;}; 00407 bmark_t& boundaryMark() {return boundary_mark;}; 00408 void reinit(const Type& t, const bmark_t& m) { 00409 boundary_type = t; 00410 boundary_mark = m; 00411 } 00412 BoundaryCondition<value_type,DIM,DOW,TDIM>& operator=(const BoundaryCondition<value_type,DIM,DOW,TDIM>& b) { 00413 boundary_type = b.boundaryType(); 00414 boundary_mark = b.boundaryMark(); 00415 return *this; 00416 }; 00417 public: 00418 virtual Number value(const afepack::Point<DOW>&) const {return 0;}; 00419 }; 00420 00421 00422 00427 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double> 00428 class BoundaryFunction : public BoundaryCondition<value_type,DIM,DOW,TDIM,Number> 00429 { 00430 public: 00431 typedef typename BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::bmark_t bmark_t; 00432 typedef typename BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::Type Type; 00433 private: 00434 bool is_newed; 00435 const Function<Number> * f; 00436 public: 00437 BoundaryFunction() : is_newed(false), f(NULL) {}; 00438 BoundaryFunction(const Type& t, 00439 const bmark_t& m, 00440 const Function<Number>& fun) : 00441 BoundaryCondition<value_type,DIM>(t,m), 00442 is_newed(false), 00443 f(&fun) {}; 00444 BoundaryFunction(const Type& t, 00445 const bmark_t& m, 00446 value_type (*fun)(const double *), 00447 std::vector<value_type> (*grad)(const double *) = NULL) : 00448 BoundaryCondition<value_type,DIM,DOW,TDIM,Number>(t,m) 00449 { 00450 is_newed = true; 00451 f = new FunctionFunction<Number>(fun, grad); 00452 }; 00453 BoundaryFunction(const BoundaryFunction<value_type,DIM,DOW,TDIM,Number>& b) : 00454 BoundaryCondition<value_type,DIM,DOW,TDIM,Number>(b), 00455 is_newed(false), 00456 f(b.f) {}; 00457 virtual ~BoundaryFunction() { 00458 if (is_newed && f != NULL) delete f; 00459 }; 00460 public: 00461 void reinit(const Type& t, 00462 const bmark_t& m, 00463 value_type (*fun)(const double *), 00464 std::vector<value_type> (*grad)(const double *) = NULL) { 00465 BoundaryCondition<value_type,DIM,DOW,TDIM,Number>::reinit(t, m); 00466 if (is_newed && f != NULL) delete f; 00467 is_newed = true; 00468 f = new FunctionFunction<Number>(fun, grad); 00469 } 00470 virtual Number value(const afepack::Point<DOW>& p) const {return f->value(p);}; 00471 virtual std::vector<Number> gradient(const afepack::Point<DOW>& p) const {return f->gradient(p);}; 00472 }; 00473 00474 00475 00482 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, typename Number=double> 00483 class BoundaryConditionAdmin : std::vector<const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *> 00484 { 00485 public: 00486 typedef typename FEMSpace<value_type,DIM,DOW,TDIM>::bmark_t bmark_t; 00487 private: 00488 std::vector<int> index_map; 00489 const FEMSpace<value_type,DIM,DOW,TDIM> * fem_space; 00490 public: 00491 BoundaryConditionAdmin(const FEMSpace<value_type,DIM,DOW,TDIM>& sp = 00492 *((const FEMSpace<value_type,DIM,DOW,TDIM> *)(NULL))) {fem_space = &sp;}; 00493 ~BoundaryConditionAdmin() {}; 00494 public: 00495 void reinit(const FEMSpace<value_type,DIM,DOW,TDIM>& sp) {fem_space = &sp; index_map.clear();}; 00496 void setFemSpace(const FEMSpace<value_type,DIM,DOW,TDIM>& sp) {fem_space = &sp;}; 00497 const FEMSpace<value_type,DIM,DOW,TDIM>& femSpace() const {return *fem_space;}; 00498 void apply(SparseMatrix<double>& A, Vector<double>& u, Vector<double>& f, 00499 bool preserve_symmetry = true); 00501 void clearEntry(Vector<double>& f); 00502 void add(const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& b); 00503 00504 bool isValid(const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& bc) const { 00505 return (&bc != NULL); 00506 }; 00507 const BoundaryCondition<value_type,DIM,DOW,TDIM,Number>& find(const bmark_t& bm) const { 00508 if ((unsigned int)bm >= index_map.size() || index_map[bm] == -1) 00509 return *((const BoundaryCondition<value_type,DIM,DOW,TDIM,Number> *)NULL); 00510 else 00511 return *((*this)[index_map[bm]]); 00512 }; 00513 }; 00514 00516 00525 template <typename value_type, int DIM> 00526 struct ElementAdditionalData : public GeometryAdditionalData<DIM> 00527 { 00528 public: 00529 std::vector<std::vector<value_type> > basis_value; 00530 std::vector<std::vector<std::vector<value_type> > > basis_gradient; 00531 }; 00532 00533 #endif //_FEMSpace_h_ 00534