dune-localfunctions  2.6-git
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/type.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::Impl::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::Impl::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::Impl::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  private:
146  void sprod ( int col1, int col2, scalar_t &ret )
147  {
148  ret = 0;
149  for( int k = 0; k <= col1; ++k )
150  {
151  for( int l = 0; l <=col2; ++l )
152  ret += Base::operator()( l, col2 ) * S( l, k ) * Base::operator()( k, col1 );
153  }
154  }
155 
156  void vmul ( std::size_t col, std::size_t rowEnd, const scalar_t &s )
157  {
158  for( std::size_t i = 0; i <= rowEnd; ++i )
159  Base::operator()( i, col ) *= s;
160  }
161 
162  void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd, const scalar_t &s )
163  {
164  for( std::size_t i = 0; i <= rowEnd; ++i )
165  Base::operator()( i, coldest ) -= s * Base::operator()( i, colsrc );
166  }
167 
168  void gramSchmidt ()
169  {
170  // setup identity
171  const std::size_t N = Base::rows();
172  for( std::size_t i = 0; i < N; ++i )
173  {
174  for( std::size_t j = 0; j < N; ++j )
175  Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
176  }
177 
178  // perform Gram-Schmidt procedure
179  scalar_t s;
180  sprod( 0, 0, s );
181  vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
182  for( std::size_t i = 1; i < N; ++i )
183  {
184  for( std::size_t k = 0; k < i; ++k )
185  {
186  sprod( i, k, s );
187  vsub( i, k, i, s );
188  }
189  sprod( i, i, s );
190  vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
191  }
192  }
193 
194  vec_t d;
195  mat_t S;
196  };
197 
198 } // namespace ONBCompute
199 
200 #endif // #ifndef DUNE_ORTHONORMALCOMPUTE_HH
Dune::LFEMatrix< scalar_t >::cols
unsigned int cols() const
Definition: lfematrix.hh:61
lfematrix.hh
ONBCompute::Integral
Definition: orthonormalcompute.hh:39
Dune::LFEMatrix< scalar_t >::operator()
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:42
ONBCompute
Definition: orthonormalcompute.hh:21
ONBCompute::ONBMatrix::row
void row(unsigned int row, Vector &vec) const
Definition: orthonormalcompute.hh:137
field.hh
ONBCompute::ONBMatrix
Definition: orthonormalcompute.hh:92
ONBCompute::ONBMatrix::mat_t
Dune::LFEMatrix< scalar_t > mat_t
Definition: orthonormalcompute.hh:100
Dune::LFEMatrix< scalar_t >::resize
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:78
ONBCompute::ONBMatrix::ONBMatrix
ONBMatrix(unsigned int order)
Definition: orthonormalcompute.hh:102
ONBCompute::factorial
scalar_t factorial(int start, int end)
Definition: orthonormalcompute.hh:25
Dune::MultiIndex::z
int z(int i) const
Definition: multiindex.hh:89
Dune::LFEMatrix< scalar_t >::rows
unsigned int rows() const
Definition: lfematrix.hh:56
Dune::field_cast
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Dune::MultiIndex
Definition: multiindex.hh:23
Dune::LFEMatrix
Definition: lfematrix.hh:15
Dune::StandardMonomialBasis
Definition: monomialbasis.hh:766
ONBCompute::ONBMatrix::vec_t
std::vector< scalar_t > vec_t
Definition: orthonormalcompute.hh:99
monomialbasis.hh
Dune
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
multiindex.hh