5#ifndef DUNE_ORTHONORMALCOMPUTE_HH
6#define DUNE_ORTHONORMALCOMPUTE_HH
20#include <dune/localfunctions/utility/field.hh>
21#include <dune/localfunctions/utility/monomialbasis.hh>
22#include <dune/localfunctions/utility/multiindex.hh>
27 template<
class scalar_t >
31 for(
int j = start; j <= end; ++j )
41 template< Dune::GeometryType::Id geometryId >
45 static constexpr int dimension = geometry.
dim();
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 )
51 return compute(alpha, p, q, std::make_integer_sequence<int,dimension>{});
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)
62 ((computeIntegral<ints>(alpha,p,q,ord)),...);
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)
71 int i = alpha.z( step );
80 p *= factorial< scalar_t >( 1, i );
81 q *= factorial< scalar_t >( step+1 + ord, step+1 + ord + i );
92 template< Dune::GeometryType::Id geometryId,
class scalar_t >
96 typedef ONBMatrix< geometryId, scalar_t > This;
100 typedef std::vector< scalar_t > vec_t;
103 explicit ONBMatrix (
unsigned int order )
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 )
115 basis.evaluate( x, y );
119 S.resize( size, size );
124 for( std::size_t i = 0; i < size; ++i )
126 for( std::size_t j = 0; j < size; ++j )
128 Integral< geometryId >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
138 template<
class Vector >
139 void row (
unsigned int row, Vector &vec )
const
143 for( std::size_t i = 0; i <
Base::rows(); ++i )
148 void sprod (
int col1,
int col2, scalar_t &ret )
151 for(
int k = 0; k <= col1; ++k )
153 for(
int l = 0; l <=col2; ++l )
154 ret += (*
this)[l][col2] * S[l][k] * (*this)[k][col1];
158 void vmul ( std::size_t col, std::size_t rowEnd,
const scalar_t &s )
160 for( std::size_t i = 0; i <= rowEnd; ++i )
161 (*
this)[i][col] *= s;
164 void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd,
const scalar_t &s )
166 for( std::size_t i = 0; i <= rowEnd; ++i )
167 (*
this)[i][coldest] -= s * (*this)[i][colsrc];
175 for( std::size_t i = 0; i < N; ++i )
177 for( std::size_t j = 0; j < N; ++j )
178 (*
this)[i][j] = scalar_t( i == j ? 1 : 0 );
184 vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
185 for( std::size_t i = 1; i < N; ++i )
187 for( std::size_t k = 0; k < i; ++k )
193 vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
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.