dune-localfunctions 2.10
Loading...
Searching...
No Matches
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 © 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/dynmatrix.hh>
16#include <dune/common/fmatrix.hh>
17
18#include <dune/geometry/type.hh>
19
23
24namespace 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 >
94 : public Dune::DynamicMatrix< scalar_t >
95 {
97 typedef Dune::DynamicMatrix< scalar_t > Base;
98
99 public:
100 typedef std::vector< scalar_t > vec_t;
101 typedef Dune::DynamicMatrix< scalar_t > mat_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( (*this)[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 += (*this)[l][col2] * S[l][k] * (*this)[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 (*this)[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 (*this)[i][coldest] -= s * (*this)[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 (*this)[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
Dune::DynamicMatrix< scalar_t > mat_t
Definition orthonormalcompute.hh:101
void row(unsigned int row, Vector &vec) const
Definition orthonormalcompute.hh:139
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:38
int z(int i) const
Definition multiindex.hh:92