3 #ifndef DUNE_ORTHONORMALCOMPUTE_HH
4 #define DUNE_ORTHONORMALCOMPUTE_HH
12 #include <dune/common/fmatrix.hh>
14 #include <dune/geometry/genericgeometry/topologytypes.hh>
24 template<
class scalar_t >
28 for(
int j = start; j <= end; ++j )
38 template<
class Topology >
41 template<
class Base >
42 struct Integral< Dune::GenericGeometry::Pyramid< Base > >
44 template<
int dim,
class scalar_t >
46 scalar_t &p, scalar_t &q )
48 const int dimension = Base::dimension+1;
49 int i = alpha.
z( Base::dimension );
50 int ord = Integral< Base >::compute( alpha, p, q );
51 p *= factorial< scalar_t >( 1, i );
52 q *= factorial< scalar_t >( dimension + ord, dimension + ord + i );
57 template<
class Base >
58 struct Integral< Dune::GenericGeometry::Prism< Base > >
60 template<
int dim,
class scalar_t >
62 scalar_t &p, scalar_t &q )
64 int i = alpha.
z( Base::dimension );
65 int ord = Integral< Base >::compute( alpha, p, q );
74 struct Integral< Dune::GenericGeometry::Point >
76 template<
int dim,
class scalar_t >
78 scalar_t &p, scalar_t &q )
91 template<
class Topology,
class scalar_t >
99 typedef std::vector< scalar_t >
vec_t;
105 const unsigned int dim = Topology::dimension;
108 const std::size_t size = basis.
size();
109 std::vector< Dune::FieldVector< MI, 1 > > y( size );
110 Dune::FieldVector< MI, dim > x;
111 for(
unsigned int i = 0; i < dim; ++i )
122 for( std::size_t i = 0; i < size; ++i )
124 for( std::size_t j = 0; j < size; ++j )
126 Integral< Topology >::compute( y[ i ][ 0 ] * y[ j ][ 0 ], p, q );
136 template<
class Vector >
137 void row (
unsigned int row, Vector &vec )
const
141 for( std::size_t i = 0; i <
Base::rows(); ++i )
151 for( std::size_t i = 0; i < N; ++i )
153 for( std::size_t j = 0; j < N; ++j )
156 for( std::size_t k = 0; k < N; ++k )
158 for( std::size_t l = 0; l < N; ++l )
161 assert( (i == j) || (fabs( Dune::field_cast< double >( prod ) ) < 1e-10) );
168 void sprod (
int col1,
int col2, scalar_t &ret )
171 for(
int k = 0; k <= col1; ++k )
173 for(
int l = 0; l <=col2; ++l )
174 ret += Base::operator()( l, col2 ) * S( l, k ) *
Base::operator()( k, col1 );
178 void vmul ( std::size_t col, std::size_t rowEnd,
const scalar_t &s )
180 for( std::size_t i = 0; i <= rowEnd; ++i )
181 Base::operator()( i, col ) *= s;
184 void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd,
const scalar_t &s )
186 for( std::size_t i = 0; i <= rowEnd; ++i )
187 Base::operator()( i, coldest ) -= s * Base::operator()( i, colsrc );
194 for( std::size_t i = 0; i < N; ++i )
196 for( std::size_t j = 0; j < N; ++j )
197 Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
203 vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
204 for( std::size_t i = 1; i < N; ++i )
206 for( std::size_t k = 0; k < i; ++k )
212 vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
222 #endif // #ifndef DUNE_ORTHONORMALCOMPUTE_HH