dune-localfunctions  2.2.0
orthonormalcompute.hh
Go to the documentation of this file.
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