dune-localfunctions  2.2.0
orthonormalcompute.hh
Go to the documentation of this file.
1 #ifndef DUNE_ORTHONORMALCOMPUTE_HH
2 #define DUNE_ORTHONORMALCOMPUTE_HH
3 
4 #include <cassert>
5 #include <iostream>
6 #include <fstream>
7 #include <iomanip>
8 #include <map>
9 
10 #include <dune/common/fmatrix.hh>
11 
12 #include <dune/geometry/genericgeometry/topologytypes.hh>
13 
18 
19 namespace ONBCompute
20 {
21 
22  template <class scalar_t>
23  scalar_t factorial(int start,int end) {
24  scalar_t ret(1);
25  for (int j=start;j<=end;j++)
26  ret*=j;
27  return ret;
28  }
29  template <class Topology>
30  struct Integral;
31  template <class Base>
32  struct Integral<Dune::GenericGeometry::Pyramid<Base> > {
33  enum {d = Base::dimension+1};
34  template <int dim,class scalar_t>
35  static int compute(const Dune::MultiIndex<dim, scalar_t>& alpha,
36  scalar_t& p,scalar_t& q) {
37  int i = alpha.z(d-1);
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);
41  return ord+i;
42  }
43  };
44  template <class Base>
45  struct Integral<Dune::GenericGeometry::Prism<Base> > {
46  enum {d = Base::dimension+1};
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) {
50  int i = alpha.z(d-1);
51  int ord = Integral<Base>::compute(alpha,p,q);
52  Integral<Base>::compute(alpha,p,q);
53  p *= 1.;
54  q *= (i+1.);
55  return ord+i;
56  }
57  };
58  template <>
59  struct Integral<Dune::GenericGeometry::Point> {
60  template <int dim,class scalar_t>
61  static int compute(const Dune::MultiIndex<dim,scalar_t>& alpha,
62  scalar_t& p,scalar_t& q) {
63  p = 1.;
64  q = 1.;
65  return 0;
66  }
67  };
68  template <class Topology, class scalar_t>
69  struct ONBMatrix : Dune::LFEMatrix<scalar_t>
70  {
72 
73  typedef std::vector< scalar_t > vec_t;
75  static const unsigned int dim=Topology::dimension;
76 
77  explicit ONBMatrix( const unsigned int order )
78  {
79  // get all multiindecies for monomial basis
82  Basis basis( order );
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) {
87  x[i].set(i);
88  }
89  basis.evaluate( x, y );
90  // set bounds of data
91  Base::resize(size,size);
92  S.resize(size,size);
93  d.resize(size);
94  // Aufstellen der Matrix fuer die Bilinearform xSy: S_ij = int_A x^(i+j)
95  scalar_t p,q;
96  for( unsigned int i=0;i<size;++i )
97  {
98  for( unsigned int j=0;j<size;j++ )
99  {
100  Integral<Topology>::compute(y[i]*y[j],p,q);
101  S(i,j) = p;
102  S(i,j) /= q;
103  }
104  }
105  gramSchmidt();
106  }
107  template <class Vector>
108  void row( const unsigned int row, Vector &vec ) const
109  {
110  // transposed matrix is required
111  assert(row<Base::cols());
112  for (unsigned int i=0;i<Base::rows();++i)
113  Dune::field_cast(Base::operator()(i,row), vec[i]);
114  }
115  bool test() {
116  bool ret = true;
117  const unsigned int N = Base::rows();
118  // Nun schauen wir noch, ob wirklich C^T S C = E gilt
119  scalar_t prod;
120  for (unsigned int i=0;i<N;++i) {
121  for (unsigned int j=0;j<N;j++) {
122  prod = 0;
123  for (unsigned int k=0;k<N;k++) {
124  for (unsigned int l=0;l<N;l++) {
125  prod += Base::operator()(l,i)*S(l,k)*Base::operator()(k,j);
126  }
127  }
128  assert((i==j)?1: fabs( Dune::field_cast<double>(prod) )<1e-10);
129  }
130  }
131  return ret;
132  }
133  private:
134  void sprod(int col1,int col2, scalar_t& ret)
135  {
136  ret = 0;
137  for (int k=0;k<=col1;k++) {
138  for (int l=0;l<=col2;l++) {
139  ret += Base::operator()(l,col2)*this->S(l,k)*Base::operator()(k,col1);
140  }
141  }
142  }
143  void vmul(unsigned int col, unsigned int rowEnd, scalar_t &ret)
144  {
145  for (unsigned int i=0;i<=rowEnd;++i)
146  Base::operator()(i,col) *= ret;
147  }
148  void vsub(unsigned int coldest, unsigned int colsrc,
149  unsigned int rowEnd, scalar_t &ret)
150  {
151  for (unsigned int i=0;i<=rowEnd;++i)
152  Base::operator()(i,coldest) -= ret*Base::operator()(i,colsrc);
153  }
154  void gramSchmidt()
155  {
156  const unsigned int N = Base::rows();
157  scalar_t s;
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;
161  sprod(0,0,s);
162  s = 1./sqrt(s);
163  vmul(0,0,s);
164  for (unsigned int i=0;i<N-1;i++) {
165  for (unsigned int k=0;k<=i;++k) {
166  sprod(i+1,k,s);
167  vsub(i+1,k,i+1,s);
168  }
169  sprod(i+1,i+1,s);
170  s = 1./sqrt(s);
171  vmul(i+1,i+1,s);
172  }
173  }
174 
175  vec_t d;
176  mat_t S;
177  };
178 
179 }
180 #endif