1 #ifndef DUNE_ORTHONORMALCOMPUTE_HH
2 #define DUNE_ORTHONORMALCOMPUTE_HH
10 #include <dune/common/fmatrix.hh>
12 #include <dune/geometry/genericgeometry/topologytypes.hh>
22 template<
class scalar_t >
26 for(
int j = start; j <= end; ++j )
36 template<
class Topology >
39 template<
class Base >
40 struct Integral< Dune::GenericGeometry::Pyramid< Base > >
42 template<
int dim,
class scalar_t >
44 scalar_t &p, scalar_t &q )
46 const int dimension = Base::dimension+1;
47 int i = alpha.
z( Base::dimension );
49 p *= factorial< scalar_t >( 1, i );
50 q *= factorial< scalar_t >( dimension + ord, dimension + ord + i );
55 template<
class Base >
56 struct Integral< Dune::GenericGeometry::Prism< Base > >
58 template<
int dim,
class scalar_t >
60 scalar_t &p, scalar_t &q )
62 int i = alpha.
z( Base::dimension );
72 struct Integral< Dune::GenericGeometry::Point >
74 template<
int dim,
class scalar_t >
76 scalar_t &p, scalar_t &q )
89 template<
class Topology,
class scalar_t >
97 typedef std::vector< scalar_t >
vec_t;
103 const unsigned int dim = Topology::dimension;
106 const std::size_t size = basis.
size();
107 std::vector< Dune::FieldVector< MI, 1 > > y( size );
108 Dune::FieldVector< MI, dim > x;
109 for(
unsigned int i = 0; i < dim; ++i )
120 for( std::size_t i = 0; i < size; ++i )
122 for( std::size_t j = 0; j < size; ++j )
134 template<
class Vector >
135 void row (
unsigned int row, Vector &vec )
const
139 for( std::size_t i = 0; i <
Base::rows(); ++i )
149 for( std::size_t i = 0; i < N; ++i )
151 for( std::size_t j = 0; j < N; ++j )
154 for( std::size_t k = 0; k < N; ++k )
156 for( std::size_t l = 0; l < N; ++l )
159 assert( (i == j) || (fabs( Dune::field_cast< double >( prod ) ) < 1e-10) );
166 void sprod (
int col1,
int col2, scalar_t &ret )
169 for(
int k = 0; k <= col1; ++k )
171 for(
int l = 0; l <=col2; ++l )
172 ret += Base::operator()( l, col2 ) * S( l, k ) *
Base::operator()( k, col1 );
176 void vmul ( std::size_t col, std::size_t rowEnd,
const scalar_t &s )
178 for( std::size_t i = 0; i <= rowEnd; ++i )
179 Base::operator()( i, col ) *= s;
182 void vsub ( std::size_t coldest, std::size_t colsrc, std::size_t rowEnd,
const scalar_t &s )
184 for( std::size_t i = 0; i <= rowEnd; ++i )
185 Base::operator()( i, coldest ) -= s * Base::operator()( i, colsrc );
192 for( std::size_t i = 0; i < N; ++i )
194 for( std::size_t j = 0; j < N; ++j )
195 Base::operator()( i, j ) = scalar_t( i == j ? 1 : 0 );
201 vmul( 0, 0, scalar_t( 1 ) / sqrt( s ) );
202 for( std::size_t i = 1; i < N; ++i )
204 for( std::size_t k = 0; k < i; ++k )
210 vmul( i, i, scalar_t( 1 ) / sqrt( s ) );
220 #endif // #ifndef DUNE_ORTHONORMALCOMPUTE_HH