dune-localfunctions
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil -*- 00002 #ifndef DUNE_MONOMLOCALBASIS_HH 00003 #define DUNE_MONOMLOCALBASIS_HH 00004 00005 #include <cassert> 00006 00007 #include <dune/common/fmatrix.hh> 00008 00009 #include"../common/localbasis.hh" 00010 00011 namespace Dune 00012 { 00013 namespace MonomImp { 00017 template<int d, int k> 00018 struct Size { 00019 enum { val = Size<d,k-1>::val+Size<d-1,k>::val }; 00020 }; 00021 template<int d> 00022 struct Size<d, 0> { 00023 enum { val = 1 }; 00024 }; 00025 template<int k> 00026 struct Size<0, k> { 00027 enum { val = 1 }; 00028 }; 00029 template<> 00030 struct Size<0, 0> { 00031 enum { val = 1 }; 00032 }; 00033 00034 template<class T> 00035 T ipow(T base, int exp) 00036 { 00037 T result(1); 00038 while (exp) 00039 { 00040 if (exp & 1) 00041 result *= base; 00042 exp >>= 1; 00043 base *= base; 00044 } 00045 return result; 00046 } 00047 00049 template <typename Traits> 00050 class EvalAccess { 00051 std::vector<typename Traits::RangeType> &out; 00052 #ifndef NDEBUG 00053 unsigned int first_unused_index; 00054 #endif 00055 00056 public: 00057 EvalAccess(std::vector<typename Traits::RangeType> &out_) 00058 : out(out_) 00059 #ifndef NDEBUG 00060 , first_unused_index(0) 00061 #endif 00062 { } 00063 #ifndef NDEBUG 00064 ~EvalAccess() { 00065 assert(first_unused_index == out.size()); 00066 } 00067 #endif 00068 typename Traits::RangeFieldType &operator[](unsigned int index) 00069 { 00070 assert(index < out.size()); 00071 #ifndef NDEBUG 00072 if(first_unused_index <= index) 00073 first_unused_index = index+1; 00074 #endif 00075 return out[index][0]; 00076 } 00077 }; 00078 00080 template <typename Traits> 00081 class JacobianAccess { 00082 std::vector<typename Traits::JacobianType> &out; 00083 unsigned int row; 00084 #ifndef NDEBUG 00085 unsigned int first_unused_index; 00086 #endif 00087 00088 public: 00089 JacobianAccess(std::vector<typename Traits::JacobianType> &out_, 00090 unsigned int row_) 00091 : out(out_), row(row_) 00092 #ifndef NDEBUG 00093 , first_unused_index(0) 00094 #endif 00095 { } 00096 #ifndef NDEBUG 00097 ~JacobianAccess() { 00098 assert(first_unused_index == out.size()); 00099 } 00100 #endif 00101 typename Traits::RangeFieldType &operator[](unsigned int index) 00102 { 00103 assert(index < out.size()); 00104 #ifndef NDEBUG 00105 if(first_unused_index <= index) 00106 first_unused_index = index+1; 00107 #endif 00108 return out[index][0][row]; 00109 } 00110 }; 00111 00124 template <typename Traits, int c> 00125 struct Evaluate 00126 { 00127 enum { 00129 d = Traits::dimDomain - c 00130 }; 00137 template <typename Access> 00138 static void eval ( 00139 const typename Traits::DomainType &in, 00142 const array<int, Traits::dimDomain> &derivatives, 00145 typename Traits::RangeFieldType prod, 00147 int bound, 00149 int& index, 00151 Access &access) 00152 { 00153 // start with the highest exponent for this dimension, then work down 00154 for (int e = bound; e >= 0; --e) 00155 { 00156 // the rest rest of the available exponents, to be used by the other 00157 // dimensions 00158 int newbound = bound - e; 00159 if(e < derivatives[d]) 00160 Evaluate<Traits,c-1>:: 00161 eval(in, derivatives, 0, newbound, index, access); 00162 else { 00163 int coeff = 1; 00164 for(int i = e - derivatives[d] + 1; i <= e; ++i) 00165 coeff *= i; 00166 // call the evaluator for the next dimension 00167 Evaluate<Traits,c-1>:: 00168 eval(// pass the coordinate and the derivatives unchanged 00169 in, derivatives, 00170 // also pass the product accumulated so far, but also 00171 // include the current dimension 00172 prod * ipow(in[d], e-derivatives[d]) * coeff, 00173 // pass the number of remaining exponents to the next 00174 // dimension 00175 newbound, 00176 // pass the next index to fill and the output access 00177 // wrapper 00178 index, access); 00179 } 00180 } 00181 } 00182 }; 00183 00188 template <typename Traits> 00189 struct Evaluate<Traits, 1> 00190 { 00191 enum { d = Traits::dimDomain-1 }; 00193 template <typename Access> 00194 static void eval (const typename Traits::DomainType &in, 00195 const array<int, Traits::dimDomain> &derivatives, 00196 typename Traits::RangeFieldType prod, 00197 int bound, int& index, Access &access) 00198 { 00199 if(bound < derivatives[d]) 00200 prod = 0; 00201 else { 00202 int coeff = 1; 00203 for(int i = bound - derivatives[d] + 1; i <= bound; ++i) 00204 coeff *= i; 00205 prod *= ipow(in[d], bound-derivatives[d]) * coeff; 00206 } 00207 access[index] = prod; 00208 ++index; 00209 } 00210 }; 00211 00212 } //namespace MonomImp 00213 00228 template<class D, class R, unsigned int d, unsigned int p, unsigned diffOrder = p> 00229 class MonomLocalBasis 00230 { 00231 enum { static_size = MonomImp::Size<d,p>::val }; 00232 00233 public: 00235 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>, 00236 Dune::FieldMatrix<R,1,d>,diffOrder> Traits; 00237 00239 unsigned int size () const 00240 { 00241 return static_size; 00242 } 00243 00245 inline void evaluateFunction (const typename Traits::DomainType& in, 00246 std::vector<typename Traits::RangeType>& out) const 00247 { 00248 evaluate<0>(array<int, 0>(), in, out); 00249 } 00250 00252 template<unsigned int k> 00253 inline void evaluate (const array<int,k>& directions, 00254 const typename Traits::DomainType& in, 00255 std::vector<typename Traits::RangeType>& out) const 00256 { 00257 out.resize(size()); 00258 int index = 0; 00259 array<int, d> derivatives; 00260 for(unsigned int i = 0; i < d; ++i) derivatives[i] = 0; 00261 for(unsigned int i = 0; i < k; ++i) ++derivatives[directions[i]]; 00262 MonomImp::EvalAccess<Traits> access(out); 00263 for(unsigned int lp = 0; lp <= p; ++lp) 00264 MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, 00265 access); 00266 } 00267 00269 inline void 00270 evaluateJacobian (const typename Traits::DomainType& in, // position 00271 std::vector<typename Traits::JacobianType>& out) const // return value 00272 { 00273 out.resize(size()); 00274 array<int, d> derivatives; 00275 for(unsigned int i = 0; i < d; ++i) 00276 derivatives[i] = 0; 00277 for(unsigned int i = 0; i < d; ++i) 00278 { 00279 derivatives[i] = 1; 00280 int index = 0; 00281 MonomImp::JacobianAccess<Traits> access(out, i); 00282 for(unsigned int lp = 0; lp <= p; ++lp) 00283 MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access); 00284 derivatives[i] = 0; 00285 } 00286 } 00287 00289 unsigned int order () const 00290 { 00291 return p; 00292 } 00293 }; 00294 00295 } 00296 00297 #endif // DUNE_MONOMLOCALBASIS_HH