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>
25 for (
int j=start;j<=end;j++)
29 template <
class Topology>
32 struct Integral<Dune::GenericGeometry::Pyramid<Base> > {
33 enum {d = Base::dimension+1};
34 template <
int dim,
class scalar_t>
36 scalar_t& p,scalar_t& q) {
38 int ord = Integral<Base>::compute(alpha,p,q);
39 p *= factorial<scalar_t>(1,i);
40 q *= factorial<scalar_t>(d+ord,d+ord+i);
45 struct Integral<Dune::GenericGeometry::Prism<Base> > {
46 enum {d = Base::dimension+1};
47 template <
int dim,
class scalar_t>
49 scalar_t& p,scalar_t& q) {
51 int ord = Integral<Base>::compute(alpha,p,q);
52 Integral<Base>::compute(alpha,p,q);
59 struct Integral<Dune::GenericGeometry::Point> {
60 template <
int dim,
class scalar_t>
62 scalar_t& p,scalar_t& q) {
68 template <
class Topology,
class scalar_t>
73 typedef std::vector< scalar_t >
vec_t;
75 static const unsigned int dim=Topology::dimension;
83 const unsigned int size = basis.size( );
84 std::vector< Dune::FieldVector< MI,1> > y( size );
85 Dune::FieldVector< MI, dim > x;
86 for (
unsigned int i=0;i<
dim;++i) {
89 basis.evaluate( x, y );
96 for(
unsigned int i=0;i<size;++i )
98 for(
unsigned int j=0;j<size;j++ )
100 Integral<Topology>::compute(y[i]*y[j],p,q);
107 template <
class Vector>
108 void row(
const unsigned int row, Vector &vec )
const
120 for (
unsigned int i=0;i<N;++i) {
121 for (
unsigned int j=0;j<N;j++) {
123 for (
unsigned int k=0;k<N;k++) {
124 for (
unsigned int l=0;l<N;l++) {
128 assert((i==j)?1: fabs( Dune::field_cast<double>(prod) )<1e-10);
134 void sprod(
int col1,
int col2, scalar_t& ret)
137 for (
int k=0;k<=col1;k++) {
138 for (
int l=0;l<=col2;l++) {
143 void vmul(
unsigned int col,
unsigned int rowEnd, scalar_t &ret)
145 for (
unsigned int i=0;i<=rowEnd;++i)
146 Base::operator()(i,col) *= ret;
148 void vsub(
unsigned int coldest,
unsigned int colsrc,
149 unsigned int rowEnd, scalar_t &ret)
151 for (
unsigned int i=0;i<=rowEnd;++i)
152 Base::operator()(i,coldest) -= ret*Base::operator()(i,colsrc);
158 for (
unsigned int i=0;i<N;++i)
159 for (
unsigned int j=0;j<N;++j)
160 Base::operator()(i,j) = (i==j)?1:0;
164 for (
unsigned int i=0;i<N-1;i++) {
165 for (
unsigned int k=0;k<=i;++k) {