dune-localfunctions
2.2.0
|
00001 #ifndef DUNE_ORTHONORMALCOMPUTE_HH 00002 #define DUNE_ORTHONORMALCOMPUTE_HH 00003 00004 #include <cassert> 00005 #include <iostream> 00006 #include <fstream> 00007 #include <iomanip> 00008 #include <map> 00009 00010 #include <dune/common/fmatrix.hh> 00011 00012 #include <dune/geometry/genericgeometry/topologytypes.hh> 00013 00014 #include <dune/localfunctions/utility/field.hh> 00015 #include <dune/localfunctions/utility/lfematrix.hh> 00016 #include <dune/localfunctions/utility/monomialbasis.hh> 00017 #include <dune/localfunctions/utility/multiindex.hh> 00018 00019 namespace ONBCompute 00020 { 00021 00022 template <class scalar_t> 00023 scalar_t factorial(int start,int end) { 00024 scalar_t ret(1); 00025 for (int j=start;j<=end;j++) 00026 ret*=j; 00027 return ret; 00028 } 00029 template <class Topology> 00030 struct Integral; 00031 template <class Base> 00032 struct Integral<Dune::GenericGeometry::Pyramid<Base> > { 00033 enum {d = Base::dimension+1}; 00034 template <int dim,class scalar_t> 00035 static int compute(const Dune::MultiIndex<dim, scalar_t>& alpha, 00036 scalar_t& p,scalar_t& q) { 00037 int i = alpha.z(d-1); 00038 int ord = Integral<Base>::compute(alpha,p,q); 00039 p *= factorial<scalar_t>(1,i); 00040 q *= factorial<scalar_t>(d+ord,d+ord+i); 00041 return ord+i; 00042 } 00043 }; 00044 template <class Base> 00045 struct Integral<Dune::GenericGeometry::Prism<Base> > { 00046 enum {d = Base::dimension+1}; 00047 template <int dim,class scalar_t> 00048 static int compute(const Dune::MultiIndex<dim,scalar_t>& alpha, 00049 scalar_t& p,scalar_t& q) { 00050 int i = alpha.z(d-1); 00051 int ord = Integral<Base>::compute(alpha,p,q); 00052 Integral<Base>::compute(alpha,p,q); 00053 p *= 1.; 00054 q *= (i+1.); 00055 return ord+i; 00056 } 00057 }; 00058 template <> 00059 struct Integral<Dune::GenericGeometry::Point> { 00060 template <int dim,class scalar_t> 00061 static int compute(const Dune::MultiIndex<dim,scalar_t>& alpha, 00062 scalar_t& p,scalar_t& q) { 00063 p = 1.; 00064 q = 1.; 00065 return 0; 00066 } 00067 }; 00068 template <class Topology, class scalar_t> 00069 struct ONBMatrix : Dune::LFEMatrix<scalar_t> 00070 { 00071 typedef Dune::LFEMatrix<scalar_t> Base; 00072 00073 typedef std::vector< scalar_t > vec_t; 00074 typedef Dune::LFEMatrix< scalar_t > mat_t; 00075 static const unsigned int dim=Topology::dimension; 00076 00077 explicit ONBMatrix( const unsigned int order ) 00078 { 00079 // get all multiindecies for monomial basis 00080 typedef Dune::MultiIndex<dim,scalar_t> MI; 00081 typedef Dune::StandardMonomialBasis< dim, MI > Basis; 00082 Basis basis( order ); 00083 const unsigned int size = basis.size( ); 00084 std::vector< Dune::FieldVector< MI,1> > y( size ); 00085 Dune::FieldVector< MI, dim > x; 00086 for (unsigned int i=0;i<dim;++i) { 00087 x[i].set(i); 00088 } 00089 basis.evaluate( x, y ); 00090 // set bounds of data 00091 Base::resize(size,size); 00092 S.resize(size,size); 00093 d.resize(size); 00094 // Aufstellen der Matrix fuer die Bilinearform xSy: S_ij = int_A x^(i+j) 00095 scalar_t p,q; 00096 for( unsigned int i=0;i<size;++i ) 00097 { 00098 for( unsigned int j=0;j<size;j++ ) 00099 { 00100 Integral<Topology>::compute(y[i]*y[j],p,q); 00101 S(i,j) = p; 00102 S(i,j) /= q; 00103 } 00104 } 00105 gramSchmidt(); 00106 } 00107 template <class Vector> 00108 void row( const unsigned int row, Vector &vec ) const 00109 { 00110 // transposed matrix is required 00111 assert(row<Base::cols()); 00112 for (unsigned int i=0;i<Base::rows();++i) 00113 Dune::field_cast(Base::operator()(i,row), vec[i]); 00114 } 00115 bool test() { 00116 bool ret = true; 00117 const unsigned int N = Base::rows(); 00118 // Nun schauen wir noch, ob wirklich C^T S C = E gilt 00119 scalar_t prod; 00120 for (unsigned int i=0;i<N;++i) { 00121 for (unsigned int j=0;j<N;j++) { 00122 prod = 0; 00123 for (unsigned int k=0;k<N;k++) { 00124 for (unsigned int l=0;l<N;l++) { 00125 prod += Base::operator()(l,i)*S(l,k)*Base::operator()(k,j); 00126 } 00127 } 00128 assert((i==j)?1: fabs( Dune::field_cast<double>(prod) )<1e-10); 00129 } 00130 } 00131 return ret; 00132 } 00133 private: 00134 void sprod(int col1,int col2, scalar_t& ret) 00135 { 00136 ret = 0; 00137 for (int k=0;k<=col1;k++) { 00138 for (int l=0;l<=col2;l++) { 00139 ret += Base::operator()(l,col2)*this->S(l,k)*Base::operator()(k,col1); 00140 } 00141 } 00142 } 00143 void vmul(unsigned int col, unsigned int rowEnd, scalar_t &ret) 00144 { 00145 for (unsigned int i=0;i<=rowEnd;++i) 00146 Base::operator()(i,col) *= ret; 00147 } 00148 void vsub(unsigned int coldest, unsigned int colsrc, 00149 unsigned int rowEnd, scalar_t &ret) 00150 { 00151 for (unsigned int i=0;i<=rowEnd;++i) 00152 Base::operator()(i,coldest) -= ret*Base::operator()(i,colsrc); 00153 } 00154 void gramSchmidt() 00155 { 00156 const unsigned int N = Base::rows(); 00157 scalar_t s; 00158 for (unsigned int i=0;i<N;++i) 00159 for (unsigned int j=0;j<N;++j) 00160 Base::operator()(i,j) = (i==j)?1:0; 00161 sprod(0,0,s); 00162 s = 1./sqrt(s); 00163 vmul(0,0,s); 00164 for (unsigned int i=0;i<N-1;i++) { 00165 for (unsigned int k=0;k<=i;++k) { 00166 sprod(i+1,k,s); 00167 vsub(i+1,k,i+1,s); 00168 } 00169 sprod(i+1,i+1,s); 00170 s = 1./sqrt(s); 00171 vmul(i+1,i+1,s); 00172 } 00173 } 00174 00175 vec_t d; 00176 mat_t S; 00177 }; 00178 00179 } 00180 #endif