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