5#ifndef DUNE_ORTHONORMALCOMPUTE_HH
6#define DUNE_ORTHONORMALCOMPUTE_HH
15#include <dune/common/dynmatrix.hh>
16#include <dune/common/fmatrix.hh>
18#include <dune/geometry/type.hh>
27 template<
class scalar_t >
31 for(
int j = start; j <= end; ++j )
41 template< Dune::GeometryType::Id geometryId >
44 static constexpr Dune::GeometryType
geometry = geometryId;
47 template<
int dim,
class scalar_t >
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>
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 >
69 scalar_t &p, scalar_t &q,
int& ord)
71 int i = alpha.
z( step );
73 if constexpr (
geometry.isPrismatic(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 >
94 :
public Dune::DynamicMatrix< scalar_t >
97 typedef Dune::DynamicMatrix< scalar_t > Base;
100 typedef std::vector< scalar_t >
vec_t;
101 typedef Dune::DynamicMatrix< scalar_t >
mat_t;
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 )
118 Base::resize( size, size );
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 )
138 template<
class Vector >
139 void row (
unsigned int row, Vector &vec )
const
142 assert(
row < Base::cols() );
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];
174 const std::size_t N = Base::rows();
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 ) );
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