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