dune-localfunctions  2.2.0
raviartthomasprebasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_RAVIARTTHOMASPREBASIS_HH
2 #define DUNE_RAVIARTTHOMASPREBASIS_HH
3 #include <fstream>
4 #include <utility>
5 
6 #include <dune/geometry/genericgeometry/topologytypes.hh>
7 
10 
11 namespace Dune
12 {
13  template <unsigned int dim, class Field>
14  struct RTPreBasisFactory;
15  template <unsigned int dim, class Field>
17  {
18  static const unsigned int dimension = dim;
19 
21  typedef typename MBasisFactory::Object MBasis;
24 
25  typedef const Basis Object;
26  typedef unsigned int Key;
28  };
29 
30  template < class Topology, class Field >
31  struct RTVecMatrix;
32 
33  template <unsigned int dim, class Field>
35  : public TopologyFactory< RTPreBasisFactoryTraits< dim, Field > >
36  {
38  static const unsigned int dimension = dim;
39  typedef typename Traits::Object Object;
40  typedef typename Traits::Key Key;
41  template <unsigned int dd, class FF>
43  {
45  };
46  template< class Topology >
47  static Object *createObject ( const Key &order )
48  {
49  RTVecMatrix<Topology,Field> vecMatrix(order);
50  typename Traits::MBasis *mbasis = Traits::MBasisFactory::template create<Topology>(order+1);
51  typename remove_const<Object>::type *tmBasis =
52  new typename remove_const<Object>::type(*mbasis);
53  tmBasis->fill(vecMatrix);
54  return tmBasis;
55  }
56  };
57  template <class Topology, class Field>
58  struct RTVecMatrix
59  {
60  static const unsigned int dim = Topology::dimension;
63  RTVecMatrix(unsigned int order)
64  {
65  MIBasis basis(order+1);
66  FieldVector< MI, dim > x;
67  for( unsigned int i = 0; i < dim; ++i )
68  x[ i ].set( i, 1 );
69  std::vector< MI > val( basis.size() );
70  basis.evaluate( x, val );
71 
72  col_ = basis.size();
73  unsigned int notHomogen = 0;
74  if (order>0)
75  notHomogen = basis.sizes()[order-1];
76  unsigned int homogen = basis.sizes()[order]-notHomogen;
77  row_ = (notHomogen*dim+homogen*(dim+1))*dim;
78  row1_ = basis.sizes()[order]*dim*dim;
79  mat_ = new Field*[row_];
80  int row = 0;
81  for (unsigned int i=0;i<notHomogen+homogen;++i)
82  {
83  for (unsigned int r=0;r<dim;++r)
84  {
85  for (unsigned int rr=0;rr<dim;++rr)
86  {
87  mat_[row] = new Field[col_];
88  for (unsigned int j=0;j<col_;++j)
89  {
90  mat_[row][j] = 0.;
91  }
92  if (r==rr)
93  mat_[row][i] = 1.;
94  ++row;
95  }
96  }
97  }
98  for (unsigned int i=0;i<homogen;++i)
99  {
100  for (unsigned int r=0;r<dim;++r)
101  {
102  mat_[row] = new Field[col_];
103  for (unsigned int j=0;j<col_;++j)
104  {
105  mat_[row][j] = 0.;
106  }
107  unsigned int w;
108  MI xval = val[notHomogen+i];
109  xval *= x[r];
110  for (w=homogen+notHomogen;w<val.size();++w)
111  {
112  if (val[w] == xval)
113  {
114  mat_[row][w] = 1.;
115  break;
116  }
117  }
118  assert(w<val.size());
119  ++row;
120  }
121  }
122  }
124  {
125  for (unsigned int i=0;i<rows();++i) {
126  delete [] mat_[i];
127  }
128  delete [] mat_;
129  }
130  unsigned int cols() const {
131  return col_;
132  }
133  unsigned int rows() const {
134  return row_;
135  }
136  template <class Vector>
137  void row( const unsigned int row, Vector &vec ) const
138  {
139  const unsigned int N = cols();
140  assert( vec.size() == N );
141  for (unsigned int i=0;i<N;++i)
142  field_cast(mat_[row][i],vec[i]);
143  }
144  unsigned int row_,col_,row1_;
145  Field **mat_;
146  };
147 
148 
149 }
150 #endif // DUNE_RAVIARTTHOMASPREBASIS_HH
151