dune-localfunctions  2.3.1
orthonormalcompute.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ORTHONORMALCOMPUTE_HH
4 #define DUNE_ORTHONORMALCOMPUTE_HH
5 
6 #include <cassert>
7 #include <iostream>
8 #include <fstream>
9 #include <iomanip>
10 #include <map>
11 
12 #include <dune/common/fmatrix.hh>
13 
14 #include <dune/geometry/genericgeometry/topologytypes.hh>
15 
20 
21 namespace ONBCompute
22 {
23 
24  template< class scalar_t >
25  scalar_t factorial( int start, int end )
26  {
27  scalar_t ret( 1 );
28  for( int j = start; j <= end; ++j )
29  ret *= scalar_t( j );
30  return ret;
31  }
32 
33 
34 
35  // Integral
36  // --------
37 
38  template< class Topology >
39  struct Integral;
40 
41  template< class Base >
42  struct Integral< Dune::GenericGeometry::Pyramid< Base > >
43  {
44  template< int dim, class scalar_t >
45  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
46  scalar_t &p, scalar_t &q )
47  {
48  const int dimension = Base::dimension+1;
49  int i = alpha.z( Base::dimension );
50  int ord = Integral< Base >::compute( alpha, p, q );
51  p *= factorial< scalar_t >( 1, i );
52  q *= factorial< scalar_t >( dimension + ord, dimension + ord + i );
53  return ord + i;
54  }
55  };
56 
57  template< class Base >
58  struct Integral< Dune::GenericGeometry::Prism< Base > >
59  {
60  template< int dim, class scalar_t >
61  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
62  scalar_t &p, scalar_t &q )
63  {
64  int i = alpha.z( Base::dimension );
65  int ord = Integral< Base >::compute( alpha, p, q );
66  //Integral< Base >::compute( alpha, p, q );
67  //p *= scalar_t( 1 );
68  q *= scalar_t( i+1 );
69  return ord + i;
70  }
71  };
72 
73  template<>
74  struct Integral< Dune::GenericGeometry::Point >
75  {
76  template< int dim, class scalar_t >
77  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
78  scalar_t &p, scalar_t &q )
79  {
80  p = scalar_t( 1 );
81  q = scalar_t( 1 );
82  return 0;
83  }
84  };
85 
86 
87 
88  // ONBMatrix
89  // ---------
90 
91  template< class Topology, class scalar_t >
92  class ONBMatrix
93  : public Dune::LFEMatrix< scalar_t >
94  {
97 
98  public:
99  typedef std::vector< scalar_t > vec_t;
101 
102  explicit ONBMatrix ( unsigned int order )
103  {
104  // get all multiindecies for monomial basis
105  const unsigned int dim = Topology::dimension;
108  const std::size_t size = basis.size();
109  std::vector< Dune::FieldVector< MI, 1 > > y( size );
110  Dune::FieldVector< MI, dim > x;
111  for( unsigned int i = 0; i < dim; ++i )
112  x[ i ].set( i );
113  basis.evaluate( x, y );
114 
115  // set bounds of data
116  Base::resize( size, size );
117  S.resize( size, size );
118  d.resize( size );
119 
120  // setup matrix for bilinear form x^T S y: S_ij = int_A x^(i+j)
121  scalar_t p, q;
122  for( std::size_t i = 0; i < size; ++i )
123  {
124  for( std::size_t j = 0; j < size; ++j )
125  {
126  Integral< Topology >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
127  S( i, j ) = p;
128  S( i, j ) /= q;
129  }
130  }
131 
132  // orthonormalize
133  gramSchmidt();
134  }
135 
136  template< class Vector >
137  void row ( unsigned int row, Vector &vec ) const
138  {
139  // transposed matrix is required
140  assert( row < Base::cols() );
141  for( std::size_t i = 0; i < Base::rows(); ++i )
142  Dune::field_cast( Base::operator()( i, row ), vec[ i ] );
143  }
144 
145  bool test ()
146  {
147  bool ret = true;
148  const std::size_t N = Base::rows();
149 
150  // check that C^T S C = I
151  for( std::size_t i = 0; i < N; ++i )
152  {
153  for( std::size_t j = 0; j < N; ++j )
154  {
155  scalar_t prod( 0 );
156  for( std::size_t k = 0; k < N; ++k )
157  {
158  for( std::size_t l = 0; l < N; ++l )
159  prod += Base::operator()( l, i ) * S( l, k ) * Base::operator()( k, j );
160  }
161  assert( (i == j) || (fabs( Dune::field_cast< double >( prod ) ) < 1e-10) );
162  }
163  }
164  return ret;
165  }
166 
167  private:
168  void sprod ( int col1, int col2, scalar_t &ret )
169  {
170  ret = 0;
171  for( int k = 0; k <= col1; ++k )
172  {
173  for( int l = 0; l <=col2; ++l )
174  ret += Base::operator()( l, col2 ) * S( l, k ) * Base::operator()( k, col1 );
175  }
176  }
177 
178  void vmul ( std::size_t col, std::size_t rowEnd, const scalar_t &s )
179  {
180  for( std::size_t i = 0; i <= rowEnd; ++i )
181  Base::operator()( i, col ) *= s;
182  }
183 
184  void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd, const scalar_t &s )
185  {
186  for( std::size_t i = 0; i <= rowEnd; ++i )
187  Base::operator()( i, coldest ) -= s * Base::operator()( i, colsrc );
188  }
189 
190  void gramSchmidt ()
191  {
192  // setup identity
193  const std::size_t N = Base::rows();
194  for( std::size_t i = 0; i < N; ++i )
195  {
196  for( std::size_t j = 0; j < N; ++j )
197  Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
198  }
199 
200  // perform Gram-Schmidt procedure
201  scalar_t s;
202  sprod( 0, 0, s );
203  vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
204  for( std::size_t i = 1; i < N; ++i )
205  {
206  for( std::size_t k = 0; k < i; ++k )
207  {
208  sprod( i, k, s );
209  vsub( i, k, i, s );
210  }
211  sprod( i, i, s );
212  vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
213  }
214  }
215 
216  vec_t d;
217  mat_t S;
218  };
219 
220 } // namespace ONBCompute
221 
222 #endif // #ifndef DUNE_ORTHONORMALCOMPUTE_HH
unsigned int rows() const
Definition: lfematrix.hh:59
Dune::LFEMatrix< scalar_t > mat_t
Definition: orthonormalcompute.hh:100
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q)
Definition: orthonormalcompute.hh:45
int z(int i) const
Definition: multiindex.hh:89
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:45
unsigned int cols() const
Definition: lfematrix.hh:64
scalar_t factorial(int start, int end)
Definition: orthonormalcompute.hh:25
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:81
bool test()
Definition: orthonormalcompute.hh:145
Definition: multiindex.hh:23
std::vector< scalar_t > vec_t
Definition: orthonormalcompute.hh:99
Definition: orthonormalcompute.hh:39
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q)
Definition: orthonormalcompute.hh:77
Definition: lfematrix.hh:15
void row(unsigned int row, Vector &vec) const
Definition: orthonormalcompute.hh:137
const unsigned int size() const
Definition: monomialbasis.hh:666
Definition: orthonormalcompute.hh:92
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q)
Definition: orthonormalcompute.hh:61
ONBMatrix(unsigned int order)
Definition: orthonormalcompute.hh:102
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Definition: monomialbasis.hh:760
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:689