DUNE-FEM (unstable)

orthonormalcompute.hh
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
17
18#include <dune/geometry/type.hh>
19
20#include <dune/localfunctions/utility/field.hh>
21#include <dune/localfunctions/utility/monomialbasis.hh>
22#include <dune/localfunctions/utility/multiindex.hh>
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 >
68 static void computeIntegral ( const Dune::MultiIndex< dim, scalar_t > &alpha,
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::DynamicMatrix< scalar_t >
95 {
96 typedef ONBMatrix< geometryId, scalar_t > This;
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();
108 typedef Dune::MultiIndex< dim, scalar_t > MI;
109 Dune::StandardMonomialBasis< dim, MI > basis( order );
110 const std::size_t size = basis.size();
111 std::vector< Dune::FieldVector< MI, 1 > > y( size );
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
constexpr size_type cols() const
number of columns
Definition: densematrix.hh:715
size_type size() const
size method (number of rows)
Definition: densematrix.hh:200
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:709
constexpr size_type N() const
number of rows
Definition: densematrix.hh:697
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
void resize(size_type r, size_type c, value_type v=value_type())
resize matrix to r × c
Definition: dynmatrix.hh:106
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr unsigned int dim() const
Return dimension of the type.
Definition: type.hh:360
constexpr bool isPrismatic() const
Return true if entity was constructed with a prismatic product in the last step.
Definition: type.hh:342
This file implements a dense matrix with dynamic numbers of rows and columns.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
static constexpr T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition: math.hh:111
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)