dune-localfunctions  2.9.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 // SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5 #ifndef DUNE_ORTHONORMALCOMPUTE_HH
6 #define DUNE_ORTHONORMALCOMPUTE_HH
7 
8 #include <cassert>
9 #include <iostream>
10 #include <fstream>
11 #include <iomanip>
12 #include <utility>
13 #include <map>
14 
15 #include <dune/common/fmatrix.hh>
16 
17 #include <dune/geometry/type.hh>
18 
23 
24 namespace ONBCompute
25 {
26 
27  template< class scalar_t >
28  scalar_t factorial( int start, int end )
29  {
30  scalar_t ret( 1 );
31  for( int j = start; j <= end; ++j )
32  ret *= scalar_t( j );
33  return ret;
34  }
35 
36 
37 
38  // Integral
39  // --------
40 
41  template< Dune::GeometryType::Id geometryId >
42  struct Integral
43  {
44  static constexpr Dune::GeometryType geometry = geometryId;
45  static constexpr int dimension = geometry.dim();
46 
47  template< int dim, class scalar_t >
48  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
49  scalar_t &p, scalar_t &q )
50  {
51  return compute(alpha, p, q, std::make_integer_sequence<int,dimension>{});
52  }
53 
54  template< int dim, class scalar_t , int ...ints>
55  static int compute ( const Dune::MultiIndex< dim, scalar_t > &alpha,
56  scalar_t &p, scalar_t &q, std::integer_sequence<int,ints...> intS)
57  {
58  p = scalar_t( 1 );
59  q = scalar_t( 1 );
60 
61  int ord = 0;
62  ((computeIntegral<ints>(alpha,p,q,ord)),...);
63 
64  return ord;
65  }
66 
67  template< int step, int dim, class scalar_t >
69  scalar_t &p, scalar_t &q, int& ord)
70  {
71  int i = alpha.z( step );
72 
73  if constexpr ( geometry.isPrismatic(step))
74  {
75  //p *= scalar_t( 1 );
76  q *= scalar_t( i+1 );
77  }
78  else
79  {
80  p *= factorial< scalar_t >( 1, i );
81  q *= factorial< scalar_t >( step+1 + ord, step+1 + ord + i );
82  }
83  ord +=i;
84  }
85 
86  };
87 
88 
89  // ONBMatrix
90  // ---------
91 
92  template< Dune::GeometryType::Id geometryId, class scalar_t >
93  class ONBMatrix
94  : public Dune::LFEMatrix< scalar_t >
95  {
98 
99  public:
100  typedef std::vector< scalar_t > vec_t;
102 
103  explicit ONBMatrix ( unsigned int order )
104  {
105  // get all multiindecies for monomial basis
106  constexpr Dune::GeometryType geometry = geometryId;
107  constexpr unsigned int dim = geometry.dim();
110  const std::size_t size = basis.size();
111  std::vector< Dune::FieldVector< MI, 1 > > y( size );
112  Dune::FieldVector< MI, dim > x;
113  for( unsigned int i = 0; i < dim; ++i )
114  x[ i ].set( i );
115  basis.evaluate( x, y );
116 
117  // set bounds of data
118  Base::resize( size, size );
119  S.resize( size, size );
120  d.resize( size );
121 
122  // setup matrix for bilinear form x^T S y: S_ij = int_A x^(i+j)
123  scalar_t p, q;
124  for( std::size_t i = 0; i < size; ++i )
125  {
126  for( std::size_t j = 0; j < size; ++j )
127  {
128  Integral< geometryId >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
129  S( i, j ) = p;
130  S( i, j ) /= q;
131  }
132  }
133 
134  // orthonormalize
135  gramSchmidt();
136  }
137 
138  template< class Vector >
139  void row ( unsigned int row, Vector &vec ) const
140  {
141  // transposed matrix is required
142  assert( row < Base::cols() );
143  for( std::size_t i = 0; i < Base::rows(); ++i )
144  Dune::field_cast( Base::operator()( i, row ), vec[ i ] );
145  }
146 
147  private:
148  void sprod ( int col1, int col2, scalar_t &ret )
149  {
150  ret = 0;
151  for( int k = 0; k <= col1; ++k )
152  {
153  for( int l = 0; l <=col2; ++l )
154  ret += Base::operator()( l, col2 ) * S( l, k ) * Base::operator()( k, col1 );
155  }
156  }
157 
158  void vmul ( std::size_t col, std::size_t rowEnd, const scalar_t &s )
159  {
160  for( std::size_t i = 0; i <= rowEnd; ++i )
161  Base::operator()( i, col ) *= s;
162  }
163 
164  void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd, const scalar_t &s )
165  {
166  for( std::size_t i = 0; i <= rowEnd; ++i )
167  Base::operator()( i, coldest ) -= s * Base::operator()( i, colsrc );
168  }
169 
170  void gramSchmidt ()
171  {
172  using std::sqrt;
173  // setup identity
174  const std::size_t N = Base::rows();
175  for( std::size_t i = 0; i < N; ++i )
176  {
177  for( std::size_t j = 0; j < N; ++j )
178  Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
179  }
180 
181  // perform Gram-Schmidt procedure
182  scalar_t s;
183  sprod( 0, 0, s );
184  vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
185  for( std::size_t i = 1; i < N; ++i )
186  {
187  for( std::size_t k = 0; k < i; ++k )
188  {
189  sprod( i, k, s );
190  vsub( i, k, i, s );
191  }
192  sprod( i, i, s );
193  vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
194  }
195  }
196 
197  vec_t d;
198  mat_t S;
199  };
200 
201 } // namespace ONBCompute
202 
203 #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:159
Definition: orthonormalcompute.hh:25
scalar_t factorial(int start, int end)
Definition: orthonormalcompute.hh:28
Definition: orthonormalcompute.hh:43
static int compute(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q)
Definition: orthonormalcompute.hh:48
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:55
static void computeIntegral(const Dune::MultiIndex< dim, scalar_t > &alpha, scalar_t &p, scalar_t &q, int &ord)
Definition: orthonormalcompute.hh:68
static constexpr int dimension
Definition: orthonormalcompute.hh:45
static constexpr Dune::GeometryType geometry
Definition: orthonormalcompute.hh:44
Definition: orthonormalcompute.hh:95
ONBMatrix(unsigned int order)
Definition: orthonormalcompute.hh:103
std::vector< scalar_t > vec_t
Definition: orthonormalcompute.hh:100
void row(unsigned int row, Vector &vec) const
Definition: orthonormalcompute.hh:139
Dune::LFEMatrix< scalar_t > mat_t
Definition: orthonormalcompute.hh:101
Definition: lfematrix.hh:18
const Field & operator()(const unsigned int row, const unsigned int col) const
Definition: lfematrix.hh:44
unsigned int cols() const
Definition: lfematrix.hh:63
void resize(const unsigned int rows, const unsigned int cols)
Definition: lfematrix.hh:80
unsigned int rows() const
Definition: lfematrix.hh:58
unsigned int size() const
Definition: monomialbasis.hh:476
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:498
Definition: monomialbasis.hh:571
Definition: multiindex.hh:37
int z(int i) const
Definition: multiindex.hh:91