dune-localfunctions  2.8.0
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 <utility>
11 #include <map>
12 
13 #include <dune/common/fmatrix.hh>
14 
15 #include <dune/geometry/type.hh>
16 
21 
22 namespace ONBCompute
23 {
24 
25  template< class scalar_t >
26  scalar_t factorial( int start, int end )
27  {
28  scalar_t ret( 1 );
29  for( int j = start; j <= end; ++j )
30  ret *= scalar_t( j );
31  return ret;
32  }
33 
34 
35 
36  // Integral
37  // --------
38 
39  template< Dune::GeometryType::Id geometryId >
40  struct Integral
41  {
42  static constexpr Dune::GeometryType geometry = geometryId;
43  static constexpr int dimension = geometry.dim();
44 
45  template< int dim, class scalar_t >
46  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
47  scalar_t &p, scalar_t &q )
48  {
49  return compute(alpha, p, q, std::make_integer_sequence<int,dimension>{});
50  }
51 
52  template< int dim, class scalar_t , int ...ints>
53  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
54  scalar_t &p, scalar_t &q, std::integer_sequence<int,ints...> intS)
55  {
56  p = scalar_t( 1 );
57  q = scalar_t( 1 );
58 
59  int ord = 0;
60  ((computeIntegral<ints>(alpha,p,q,ord)),...);
61 
62  return ord;
63  }
64 
65  template< int step, int dim, class scalar_t >
67  scalar_t &p, scalar_t &q, int& ord)
68  {
69  int i = alpha.z( step );
70 
71  if constexpr ( geometry.isPrismatic(step))
72  {
73  //p *= scalar_t( 1 );
74  q *= scalar_t( i+1 );
75  }
76  else
77  {
78  p *= factorial< scalar_t >( 1, i );
79  q *= factorial< scalar_t >( step+1 + ord, step+1 + ord + i );
80  }
81  ord +=i;
82  }
83 
84  };
85 
86 
87  // ONBMatrix
88  // ---------
89 
90  template< Dune::GeometryType::Id geometryId, class scalar_t >
91  class ONBMatrix
92  : public Dune::LFEMatrix< scalar_t >
93  {
96 
97  public:
98  typedef std::vector< scalar_t > vec_t;
100 
101  explicit ONBMatrix ( unsigned int order )
102  {
103  // get all multiindecies for monomial basis
104  constexpr Dune::GeometryType geometry = geometryId;
105  constexpr unsigned int dim = geometry.dim();
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< geometryId >::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  using std::sqrt;
171  // setup identity
172  const std::size_t N = Base::rows();
173  for( std::size_t i = 0; i < N; ++i )
174  {
175  for( std::size_t j = 0; j < N; ++j )
176  Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
177  }
178 
179  // perform Gram-Schmidt procedure
180  scalar_t s;
181  sprod( 0, 0, s );
182  vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
183  for( std::size_t i = 1; i < N; ++i )
184  {
185  for( std::size_t k = 0; k < i; ++k )
186  {
187  sprod( i, k, s );
188  vsub( i, k, i, s );
189  }
190  sprod( i, i, s );
191  vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
192  }
193  }
194 
195  vec_t d;
196  mat_t S;
197  };
198 
199 } // namespace ONBCompute
200 
201 #endif // #ifndef DUNE_ORTHONORMALCOMPUTE_HH
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Definition: orthonormalcompute.hh:23
scalar_t factorial(int start, int end)
Definition: orthonormalcompute.hh:26
Definition: orthonormalcompute.hh:41
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q)
Definition: orthonormalcompute.hh:46
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q, std::integer_sequence< int, ints... > intS)
Definition: orthonormalcompute.hh:53
static void computeIntegral(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q, int &ord)
Definition: orthonormalcompute.hh:66
static constexpr int dimension
Definition: orthonormalcompute.hh:43
static constexpr Dune::GeometryType geometry
Definition: orthonormalcompute.hh:42
Definition: orthonormalcompute.hh:93
ONBMatrix(unsigned int order)
Definition: orthonormalcompute.hh:101
std::vector< scalar_t > vec_t
Definition: orthonormalcompute.hh:98
void row(unsigned int row, Vector &vec) const
Definition: orthonormalcompute.hh:137
Dune::LFEMatrix< scalar_t > mat_t
Definition: orthonormalcompute.hh:99
Definition: lfematrix.hh:16
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:42
unsigned int cols() const
Definition: lfematrix.hh:61
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:78
unsigned int rows() const
Definition: lfematrix.hh:56
unsigned int size() const
Definition: monomialbasis.hh:474
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:496
Definition: monomialbasis.hh:569
Definition: multiindex.hh:35
int z(int i) const
Definition: multiindex.hh:89