dune-localfunctions  2.8.0
raviartthomassimplexprebasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_RAVIARTTHOMASPREBASIS_HH
4 #define DUNE_RAVIARTTHOMASPREBASIS_HH
5 
6 #include <fstream>
7 #include <utility>
8 
9 #include <dune/geometry/type.hh>
10 
12 
13 namespace Dune
14 {
15  template < GeometryType::Id geometryId, class Field >
16  struct RTVecMatrix;
17 
18  template <unsigned int dim, class Field>
20  {
22  typedef typename MBasisFactory::Object MBasis;
25 
26  typedef const Basis Object;
27  typedef std::size_t Key;
28 
29  template <unsigned int dd, class FF>
31  {
33  };
34  template< GeometryType::Id geometryId >
35  static Object *create ( const Key &order )
36  {
37  RTVecMatrix<geometryId,Field> vecMatrix(order);
38  MBasis *mbasis = MBasisFactory::template create<geometryId>(order+1);
39  typename std::remove_const<Object>::type *tmBasis = new typename std::remove_const<Object>::type(*mbasis);
40  tmBasis->fill(vecMatrix);
41  return tmBasis;
42  }
43  static void release( Object *object ) { delete object; }
44  };
45 
46  template <GeometryType::Id geometryId, class Field>
47  struct RTVecMatrix
48  {
49  static constexpr GeometryType geometry = geometryId;
50  static const unsigned int dim = geometry.dim();
53  RTVecMatrix(std::size_t order)
54  {
55  /*
56  * Construction of Raviart-Thomas elements in high dimensions see "Mixed Finite Elements in \R^3" by Nedelec, 1980.
57  *
58  * Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$.
59  * The space of Raviart-Thomas functions in $n$ dimensions with index $k$ is defined as
60  *
61  * \begin{equation*}
62  * RT_k := (\P_{k-1})^n \oplus \widetilde \P_k x
63  * \end{equation*}
64  * with $x=(x_1,x_2,\dots, x_n)$ in $n$ dimensions and $\widetilde \P_k$ the homogeneous polynomials of degree $k$.
65  *
66  * For $RT_k$ holds
67  * \begin{equation*}
68  * (\P_{k-1})^n \subset RT_k \subset (\P_k)^n.
69  * \end{equation*}
70  *
71  * We construct $(\P_k)^n$ and and only use the monomials contained in $RT_k$.
72  *
73  */
74 
75  MIBasis basis(order+1);
76  FieldVector< MI, dim > x;
77  /*
78  * Init MultiIndices
79  * x[0]=(1,0,0) x
80  * x[1]=(0,1,0) y
81  * x[2]=(0,0,1) z
82  */
83  for( unsigned int i = 0; i < dim; ++i )
84  x[ i ].set( i, 1 );
85  std::vector< MI > val( basis.size() );
86 
87  // val now contains all monomials in $n$ dimensions with degree $\leq order+1$
88  basis.evaluate( x, val );
89 
90  col_ = basis.size();
91 
92  // get $\dim (\P_{order-1})$
93  unsigned int notHomogen = 0;
94  if (order>0)
95  notHomogen = basis.sizes()[order-1];
96 
97  // get $\dim \widetilde (\P_order)$
98  unsigned int homogen = basis.sizes()[order]-notHomogen;
99 
100  /*
101  *
102  * The set $RT_k$ is defined as
103  *
104  * \begin{equation}
105  * RT_k := (\P_k)^dim + \widetilde \P_k x with x\in \R^n.
106  * \end{equation}
107  *
108  * The space $\P_k$ is split in $\P_k = \P_{k-1} + \widetilde \P_k$.
109  *
110  * \begin{align}
111  * RT_k &= (\P_{k-1} + \widetilde \P_k)^dim + \widetilde \P_k x with x\in \R^n
112  * &= (\P_{k-1})^n + (\widetilde \P_k)^n + \widetilde \P_k x with x\in \R^n
113  * \end{align}
114  *
115  * Thus $\dim RT_k = n * \dim \P_{k-1} + (n+1)*\dim (\widetilde \P_k)$
116  */
117 
118  // row_ = \dim RT_k *dim
119  row_ = (notHomogen*dim+homogen*(dim+1))*dim;
120  mat_ = new Field*[row_];
121  int row = 0;
122 
123  /* Assign the correct values for the nonhomogeneous polymonials $p\in (\P_{oder-1})^dim$
124  * A basis function is represented by $dim$ rows.
125  */
126  for (unsigned int i=0; i<notHomogen+homogen; ++i)
127  {
128  for (unsigned int r=0; r<dim; ++r)
129  {
130  for (unsigned int rr=0; rr<dim; ++rr)
131  {
132  // init row to zero
133  mat_[row] = new Field[col_];
134  for (unsigned int j=0; j<col_; ++j)
135  {
136  mat_[row][j] = 0.;
137  }
138  if (r==rr)
139  mat_[row][i] = 1.;
140  ++row;
141  }
142  }
143  }
144 
145  /* Assign the correct values for the homogeneous polymonials $p\in RT_k \backslash (\P_{oder-1})^dim$
146  * A basis function is represented by $dim$ rows.
147  */
148  for (unsigned int i=0; i<homogen; ++i)
149  {
150  for (unsigned int r=0; r<dim; ++r)
151  {
152  // init rows to zero
153  mat_[row] = new Field[col_];
154  for (unsigned int j=0; j<col_; ++j)
155  {
156  mat_[row][j] = 0.;
157  }
158 
159  unsigned int w;
160  // get a monomial $ p \in \widetilde \P_{order}$
161  MI xval = val[notHomogen+i];
162  xval *= x[r];
163  for (w=homogen+notHomogen; w<val.size(); ++w)
164  {
165  if (val[w] == xval)
166  {
167  mat_[row][w] = 1.;
168  break;
169  }
170  }
171  assert(w<val.size());
172  ++row;
173  }
174  }
175  }
176 
178  {
179  for (unsigned int i=0; i<rows(); ++i) {
180  delete [] mat_[i];
181  }
182  delete [] mat_;
183  }
184 
185  unsigned int cols() const {
186  return col_;
187  }
188 
189  unsigned int rows() const {
190  return row_;
191  }
192 
193  template <class Vector>
194  void row( const unsigned int row, Vector &vec ) const
195  {
196  const unsigned int N = cols();
197  assert( vec.size() == N );
198  for (unsigned int i=0; i<N; ++i)
199  field_cast(mat_[row][i],vec[i]);
200  }
201  unsigned int row_,col_;
202  Field **mat_;
203  };
204 
205 
206 }
207 #endif // DUNE_RAVIARTTHOMASPREBASIS_HH
Definition: bdfmcube.hh:16
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
Definition: raviartthomassimplexprebasis.hh:48
static const unsigned int dim
Definition: raviartthomassimplexprebasis.hh:50
~RTVecMatrix()
Definition: raviartthomassimplexprebasis.hh:177
Field ** mat_
Definition: raviartthomassimplexprebasis.hh:202
RTVecMatrix(std::size_t order)
Definition: raviartthomassimplexprebasis.hh:53
unsigned int cols() const
Definition: raviartthomassimplexprebasis.hh:185
unsigned int row_
Definition: raviartthomassimplexprebasis.hh:201
MultiIndex< dim, Field > MI
Definition: raviartthomassimplexprebasis.hh:51
void row(const unsigned int row, Vector &vec) const
Definition: raviartthomassimplexprebasis.hh:194
MonomialBasis< geometryId, MI > MIBasis
Definition: raviartthomassimplexprebasis.hh:52
static constexpr GeometryType geometry
Definition: raviartthomassimplexprebasis.hh:49
unsigned int rows() const
Definition: raviartthomassimplexprebasis.hh:189
unsigned int col_
Definition: raviartthomassimplexprebasis.hh:201
Definition: raviartthomassimplexprebasis.hh:20
const Basis Object
Definition: raviartthomassimplexprebasis.hh:26
MonomialBasisProvider< dim, Field > MBasisFactory
Definition: raviartthomassimplexprebasis.hh:21
PolynomialBasisWithMatrix< EvalMBasis, SparseCoeffMatrix< Field, dim > > Basis
Definition: raviartthomassimplexprebasis.hh:24
static Object * create(const Key &order)
Definition: raviartthomassimplexprebasis.hh:35
MBasisFactory::Object MBasis
Definition: raviartthomassimplexprebasis.hh:22
std::size_t Key
Definition: raviartthomassimplexprebasis.hh:27
static void release(Object *object)
Definition: raviartthomassimplexprebasis.hh:43
StandardEvaluator< MBasis > EvalMBasis
Definition: raviartthomassimplexprebasis.hh:23
Definition: raviartthomassimplexprebasis.hh:31
MonomialBasisProvider< dd, FF > Type
Definition: raviartthomassimplexprebasis.hh:32
Definition: basisevaluator.hh:129
Definition: monomialbasis.hh:438
unsigned int size() const
Definition: monomialbasis.hh:474
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:496
const unsigned int * sizes(unsigned int order) const
Definition: monomialbasis.hh:463
Definition: monomialbasis.hh:789
Definition: multiindex.hh:35
Definition: polynomialbasis.hh:335