AFEPack
|
00001 00002 // DGFEMSpace.h : by R.Lie, July 7, 2003 00003 // 00004 00005 #ifndef _DGFEMSpace_h_ 00006 #define _DGFEMSpace_h_ 00007 00008 #include <lac/full_matrix.h> 00009 00010 #include "Miscellaneous.h" 00011 #include "Geometry.h" 00012 #include "TemplateElement.h" 00013 #include "FEMSpace.h" 00014 00015 template <int TDIM, int DIM> class TemplateDGElement; 00016 template <class value_type, int DIM, int DOW, int TDIM, int TDIM1> class DGElement; 00017 template <class value_type, int DIM, int DOW, int TDIM, int TDIM1> class DGFEMSpace; 00018 00019 template <class value_type, int DIM, int DOW, int TDIM, int TDIM1> std::vector<double> unitOutNormal(const afepack::Point<DIM>&, const Element<value_type,DIM,DOW,TDIM>&, const DGElement<value_type,DIM,DOW,TDIM,TDIM1>&); 00020 template <class value_type, int DIM, int DOW, int TDIM, int TDIM1> std::vector<std::vector<double> > unitOutNormal(const std::vector<afepack::Point<DIM> >&, const Element<value_type,DIM,DOW,TDIM>&, const DGElement<value_type,DIM,DOW,TDIM,TDIM1>&); 00021 00023 00024 template <int TDIM, int DIM=TDIM+1> 00025 class TemplateDGElement 00026 { 00027 private: 00028 TemplateGeometry<TDIM> * geo; 00029 CoordTransform<TDIM,DIM> * ct; 00030 public: 00031 TemplateDGElement(TemplateGeometry<TDIM>& _geo = *((TemplateGeometry<TDIM> *)NULL), 00032 CoordTransform<TDIM,DIM>& _ct = *((CoordTransform<TDIM,DIM> *)NULL)) : 00033 geo(&_geo), ct(&_ct) {}; 00034 TemplateDGElement(const TemplateDGElement<TDIM,DIM>& other) : 00035 geo(other.geo), ct(other.ct) {}; 00036 ~TemplateDGElement() {}; 00037 public: 00038 void reinit(TemplateGeometry<TDIM>& _geo = *((TemplateGeometry<TDIM> *)NULL), 00039 CoordTransform<TDIM,DIM>& _ct = *((CoordTransform<TDIM,DIM> *)NULL)) {geo = &_geo; ct = &_ct;}; 00040 const TemplateGeometry<TDIM>& geometry() const {return *geo;}; 00041 TemplateGeometry<TDIM>& geometry() {return *geo;}; 00042 const CoordTransform<TDIM,DIM>& coordTransform() const {return *ct;}; 00043 CoordTransform<TDIM,DIM>& coordTransform() {return *ct;}; 00044 const std::vector<afepack::Point<TDIM> >& vertexArray() const {return geo->vertexArray();}; 00045 const QuadratureInfoAdmin<TDIM>& quadratureInfo() const {return geo->quadratureInfo();}; 00046 QuadratureInfoAdmin<TDIM>& quadratureInfo() {return geo->quadratureInfo();}; 00047 const QuadratureInfo<TDIM>& findQuadratureInfo(const int& i) const {return geo->findQuadratureInfo(i);}; 00048 double volume() const {return geo->volume();}; 00049 }; 00050 00051 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, int TDIM1=DIM-1> 00052 class DGElement 00053 { 00054 public: 00055 enum { dim = DIM, dow = DOW, tdim = TDIM, tdim1 = TDIM1 }; 00056 typedef value_type value_t; 00057 typedef DGFEMSpace<value_t,DIM,DOW,TDIM,TDIM1> fe_space_t; 00058 typedef Element<value_t,DIM,DOW,TDIM> element_t; 00059 typedef TemplateDGElement<TDIM1,DOW> dg_template_t; 00060 typedef DGElement<value_t,DIM,DOW,TDIM,TDIM1> dg_element_t; 00061 typedef typename Mesh<DIM,DOW>::bmark_t bmark_t; 00062 typedef CoordTransform<TDIM1,DOW> coord_trans_t; 00063 00064 private: 00065 fe_space_t * sp; 00066 int geometry_index; 00067 int template_element_index; 00068 std::vector<std::vector<int> > geo_img; 00069 element_t * neigh[2]; 00070 int bnd_idx[2]; 00072 public: 00073 DGElement(fe_space_t& = *((fe_space_t *)NULL)); 00074 DGElement(const dg_element_t&); 00075 ~DGElement(); 00076 public: 00077 void reinit(fe_space_t&, const int&, const int&); 00078 void reinit(const int&, const int&); 00079 dg_element_t& operator=(const dg_element_t&); 00080 const fe_space_t& femSpace() const; 00081 fe_space_t& femSpace(); 00082 const int& index() const; 00083 int& index(); 00084 const int& boundaryIndex(const u_int&) const; 00085 int& boundaryIndex(const u_int&); 00086 GeometryBM& geometry() const; 00087 bmark_t boundaryMark() const; 00088 void geometry(const GeometryBM&); 00089 void geometry(const int&); 00090 TemplateDGElement<TDIM1,DOW>& templateElement() const; 00091 void templateElement(const int&); 00092 TemplateGeometry<TDIM1>& templateGeometry() const; 00093 const std::vector<std::vector<int> >& geometryImage() const; 00094 std::vector<std::vector<int> >& geometryImage(); 00095 void buildGeometryImage(); 00096 void buildVertexArray(std::vector<afepack::Point<DOW> >&) const; 00097 const double ** buildVertexArray() const; 00099 afepack::Point<DOW> local_to_global(const afepack::Point<TDIM1>&) const; 00100 afepack::Point<TDIM1> global_to_local(const afepack::Point<DOW>&) const; 00101 double local_to_global_jacobian(const afepack::Point<TDIM1>&) const; 00102 double global_to_local_jacobian(const afepack::Point<DOW>&) const; 00103 std::vector<afepack::Point<DOW> > local_to_global(const std::vector<afepack::Point<TDIM1> >&) const; 00104 std::vector<afepack::Point<TDIM1> > global_to_local(const std::vector<afepack::Point<DOW> >&) const; 00105 std::vector<double> local_to_global_jacobian(const std::vector<afepack::Point<TDIM1> >&) const; 00106 std::vector<double> global_to_local_jacobian(const std::vector<afepack::Point<DOW> >&) const; 00108 const element_t * p_neighbourElement(const int& i) const; 00109 element_t * p_neighbourElement(const int& i); 00110 const element_t& neighbourElement(const int& i) const; 00111 element_t& neighbourElement(const int& i); 00112 00113 const QuadratureInfo<TDIM1>& findQuadratureInfo(const int&) const; 00114 std::vector<double> unitNormal() const; 00115 00116 public: 00117 DeclException1(ExcMeshData, std::string, << "Mesh data uncompatible: " << arg1); 00118 public: 00119 friend std::vector<double> unitOutNormal <>(const afepack::Point<DIM>&, 00120 const element_t&, 00121 const dg_element_t&); 00122 friend std::vector<std::vector<double> > unitOutNormal <>(const std::vector<afepack::Point<DIM> >&, 00123 const element_t&, 00124 const dg_element_t&); 00125 friend class DGFEMSpace<value_type, DIM, DOW, TDIM, TDIM1>; 00126 #ifdef __SERIALIZATION__ 00127 template <class T> T * 00128 new_property(const property_id_t<T>& pid) const { 00129 return femSpace().new_property(*this, pid); 00130 } 00131 template <class T> T * 00132 get_property(const property_id_t<T>& pid) const { 00133 return femSpace().get_property(*this, pid); 00134 } 00135 template <class T> void 00136 free_property(const property_id_t<T>& pid) const { 00137 femSpace().free_property(*this, pid); 00138 } 00139 #endif // __SERIALIZATION__ 00140 }; 00141 00142 template <class value_type, int DIM, int DOW=DIM, int TDIM=DIM, int TDIM1=DIM-1> 00143 class DGFEMSpace : public FEMSpace<value_type, DIM, DOW, TDIM> 00144 { 00145 public: 00146 enum { dim = DIM, dow = DOW, tdim = TDIM, tdim1 = TDIM1 }; 00147 typedef value_type value_t; 00148 typedef DGFEMSpace<value_t,DIM,DOW,TDIM,TDIM1> fe_space_t; 00149 typedef DGElement<value_t,DIM,DOW,TDIM,TDIM1> dg_element_t; 00150 typedef FEMSpace<value_t, DIM, DOW, TDIM> base_t; 00151 typedef TemplateDGElement<TDIM1,DOW> dg_template_t; 00152 typedef typename base_t::template_t template_t; 00153 typedef typename base_t::element_t element_t; 00154 00155 typedef typename std::vector<dg_element_t>::iterator DGElementIterator; 00156 typedef typename std::vector<dg_element_t>::const_iterator ConstDGElementIterator; 00157 00158 private: 00159 std::vector<dg_template_t> * tmp_dgele; 00160 std::vector<dg_element_t> dgele; 00161 public: 00162 explicit DGFEMSpace(Mesh<DIM,DOW>& = *((Mesh<DIM,DOW> *)NULL), 00163 std::vector<template_t>& = *((std::vector<template_t> *)NULL), 00164 std::vector<dg_template_t>& = *((std::vector<dg_template_t> *)NULL)); 00165 DGFEMSpace(const fe_space_t&); 00166 virtual ~DGFEMSpace(); 00167 public: 00168 void reinit(Mesh<DIM,DOW>&, 00169 std::vector<template_t>& = *((std::vector<template_t> *)NULL), 00170 std::vector<dg_template_t>& = *((std::vector<dg_template_t> *)NULL)); 00171 fe_space_t& operator=(const fe_space_t&); 00172 const std::vector<dg_template_t>& templateDGElement() const {return *tmp_dgele;}; 00173 std::vector<dg_template_t>& templateDGElement() {return *tmp_dgele;} 00174 const dg_template_t& templateDGElement(const int& i) const {return (*tmp_dgele)[i];}; 00175 dg_template_t& templateDGElement(const int& i) {return (*tmp_dgele)[i];}; 00176 int n_DGElement() const {return dgele.size();}; 00177 const std::vector<dg_element_t>& dgElement() const {return dgele;}; 00178 std::vector<dg_element_t>& dgElement() {return dgele;}; 00179 const dg_element_t& dgElement(const int& i) const {return dgele[i];}; 00180 dg_element_t& dgElement(const int& i) {return dgele[i];}; 00181 virtual void buildDGElement(); 00182 DGElementIterator beginDGElement() {return dgElement().begin();}; 00183 DGElementIterator endDGElement() {return dgElement().end();}; 00184 ConstDGElementIterator beginDGElement() const {return dgElement().begin();}; 00185 ConstDGElementIterator endDGElement() const {return dgElement().end();}; 00186 00187 #ifdef __SERIALIZATION__ 00188 template <class T> T * 00189 new_property(const dg_element_t& ele, 00190 const property_id_t<T>& pid) const { 00191 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00192 const reg_mesh_t& mesh = dynamic_cast<const reg_mesh_t&>(this->mesh()); 00193 return mesh.template new_property<T,DIM-1>(ele.geometry().index(), pid); 00194 } 00195 template <class T> T * 00196 get_property(const dg_element_t& ele, 00197 const property_id_t<T>& pid) const { 00198 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00199 const reg_mesh_t& mesh = dynamic_cast<const reg_mesh_t&>(this->mesh()); 00200 return mesh.template get_property<T,DIM-1>(ele.geometry().index(), pid); 00201 } 00202 template <class T> void 00203 free_property(const dg_element_t& ele, 00204 const property_id_t<T>& pid) const { 00205 typedef RegularMesh<DIM,DOW> reg_mesh_t; 00206 const reg_mesh_t& mesh = dynamic_cast<const reg_mesh_t&>(this->mesh()); 00207 mesh.template free_property<T,DIM-1>(ele.geometry().index(), pid); 00208 } 00209 00210 using base_t::new_property; 00211 using base_t::get_property; 00212 using base_t::free_property; 00213 #endif // __SERIALIZATION__ 00214 public: 00215 friend class DGElement<value_t,DIM,DOW,TDIM,TDIM1>; 00216 }; 00217 00218 00220 00229 template <typename value_type, int DIM> 00230 struct DGElementAdditionalData : public GeometryAdditionalData<DIM> 00231 { 00232 public: 00233 std::vector<std::vector<value_type> > basis_value0; 00234 std::vector<std::vector<value_type> > basis_value1; 00235 std::vector<std::vector<double> > unit_normal; 00236 }; 00237 00238 #endif // _DGFEMSpace_h_ 00239 00240 // 00241 // end of file