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