dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_RAVIARTTHOMASPREBASIS_HH 00002 #define DUNE_RAVIARTTHOMASPREBASIS_HH 00003 #include <fstream> 00004 #include <utility> 00005 00006 #include <dune/geometry/genericgeometry/topologytypes.hh> 00007 00008 #include <dune/localfunctions/utility/lfematrix.hh> 00009 #include <dune/localfunctions/utility/polynomialbasis.hh> 00010 00011 namespace Dune 00012 { 00013 template <unsigned int dim, class Field> 00014 struct RTPreBasisFactory; 00015 template <unsigned int dim, class Field> 00016 struct RTPreBasisFactoryTraits 00017 { 00018 static const unsigned int dimension = dim; 00019 00020 typedef MonomialBasisProvider<dim,Field> MBasisFactory; 00021 typedef typename MBasisFactory::Object MBasis; 00022 typedef StandardEvaluator<MBasis> EvalMBasis; 00023 typedef PolynomialBasisWithMatrix<EvalMBasis,SparseCoeffMatrix<Field,dim> > Basis; 00024 00025 typedef const Basis Object; 00026 typedef unsigned int Key; 00027 typedef RTPreBasisFactory<dim,Field> Factory; 00028 }; 00029 00030 template < class Topology, class Field > 00031 struct RTVecMatrix; 00032 00033 template <unsigned int dim, class Field> 00034 struct RTPreBasisFactory 00035 : public TopologyFactory< RTPreBasisFactoryTraits< dim, Field > > 00036 { 00037 typedef RTPreBasisFactoryTraits< dim, Field > Traits; 00038 static const unsigned int dimension = dim; 00039 typedef typename Traits::Object Object; 00040 typedef typename Traits::Key Key; 00041 template <unsigned int dd, class FF> 00042 struct EvaluationBasisFactory 00043 { 00044 typedef MonomialBasisProvider<dd,FF> Type; 00045 }; 00046 template< class Topology > 00047 static Object *createObject ( const Key &order ) 00048 { 00049 RTVecMatrix<Topology,Field> vecMatrix(order); 00050 typename Traits::MBasis *mbasis = Traits::MBasisFactory::template create<Topology>(order+1); 00051 typename remove_const<Object>::type *tmBasis = 00052 new typename remove_const<Object>::type(*mbasis); 00053 tmBasis->fill(vecMatrix); 00054 return tmBasis; 00055 } 00056 }; 00057 template <class Topology, class Field> 00058 struct RTVecMatrix 00059 { 00060 static const unsigned int dim = Topology::dimension; 00061 typedef MultiIndex<dim,Field> MI; 00062 typedef MonomialBasis<Topology,MI> MIBasis; 00063 RTVecMatrix(unsigned int order) 00064 { 00065 MIBasis basis(order+1); 00066 FieldVector< MI, dim > x; 00067 for( unsigned int i = 0; i < dim; ++i ) 00068 x[ i ].set( i, 1 ); 00069 std::vector< MI > val( basis.size() ); 00070 basis.evaluate( x, val ); 00071 00072 col_ = basis.size(); 00073 unsigned int notHomogen = 0; 00074 if (order>0) 00075 notHomogen = basis.sizes()[order-1]; 00076 unsigned int homogen = basis.sizes()[order]-notHomogen; 00077 row_ = (notHomogen*dim+homogen*(dim+1))*dim; 00078 row1_ = basis.sizes()[order]*dim*dim; 00079 mat_ = new Field*[row_]; 00080 int row = 0; 00081 for (unsigned int i=0;i<notHomogen+homogen;++i) 00082 { 00083 for (unsigned int r=0;r<dim;++r) 00084 { 00085 for (unsigned int rr=0;rr<dim;++rr) 00086 { 00087 mat_[row] = new Field[col_]; 00088 for (unsigned int j=0;j<col_;++j) 00089 { 00090 mat_[row][j] = 0.; 00091 } 00092 if (r==rr) 00093 mat_[row][i] = 1.; 00094 ++row; 00095 } 00096 } 00097 } 00098 for (unsigned int i=0;i<homogen;++i) 00099 { 00100 for (unsigned int r=0;r<dim;++r) 00101 { 00102 mat_[row] = new Field[col_]; 00103 for (unsigned int j=0;j<col_;++j) 00104 { 00105 mat_[row][j] = 0.; 00106 } 00107 unsigned int w; 00108 MI xval = val[notHomogen+i]; 00109 xval *= x[r]; 00110 for (w=homogen+notHomogen;w<val.size();++w) 00111 { 00112 if (val[w] == xval) 00113 { 00114 mat_[row][w] = 1.; 00115 break; 00116 } 00117 } 00118 assert(w<val.size()); 00119 ++row; 00120 } 00121 } 00122 } 00123 ~RTVecMatrix() 00124 { 00125 for (unsigned int i=0;i<rows();++i) { 00126 delete [] mat_[i]; 00127 } 00128 delete [] mat_; 00129 } 00130 unsigned int cols() const { 00131 return col_; 00132 } 00133 unsigned int rows() const { 00134 return row_; 00135 } 00136 template <class Vector> 00137 void row( const unsigned int row, Vector &vec ) const 00138 { 00139 const unsigned int N = cols(); 00140 assert( vec.size() == N ); 00141 for (unsigned int i=0;i<N;++i) 00142 field_cast(mat_[row][i],vec[i]); 00143 } 00144 unsigned int row_,col_,row1_; 00145 Field **mat_; 00146 }; 00147 00148 00149 } 00150 #endif // DUNE_RAVIARTTHOMASPREBASIS_HH 00151