dune-localfunctions  2.8.0
nedelecsimplexprebasis.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_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_HH
4 #define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_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 NedelecVecMatrix;
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 
35  template< GeometryType::Id geometryId >
36  static Object *create ( Key order )
37  {
38  /*
39  * The nedelec parameter begins at 1.
40  * This is the numbering used by J.C. Nedelec himself.
41  * See "Mixed Finite Elements in \R^3" published in 1980.
42  *
43  * This construction is based on the construction of Raviart-Thomas elements.
44  * There the numbering starts at 0.
45  * Because of this we reduce the order internally by 1.
46  */
47  order--;
48  NedelecVecMatrix<geometryId,Field> vecMatrix(order);
49  MBasis *mbasis = MBasisFactory::template create<geometryId>(order+1);
50  std::remove_const_t<Object>* tmBasis = new std::remove_const_t<Object>(*mbasis);
51  tmBasis->fill(vecMatrix);
52  return tmBasis;
53  }
54  static void release( Object *object ) { delete object; }
55  };
56 
57  template <GeometryType::Id geometryId, class Field>
59  {
60  static constexpr GeometryType geometry = geometryId;
61  static const unsigned int dim = geometry.dim();
64  NedelecVecMatrix(std::size_t order)
65  {
66  /*
67  * Construction of Nedelec elements see "Mixed Finite Elements in \R^3" by Nedelec, 1980.
68  *
69  * Let $\P_{n,k}$ be the space of polynomials in $n$ variables with degree $\leq k$.
70  * The space of Nedelec functions in $n$ dimensions with index $k$ is defined as
71  *
72  * \begin{equation*}
73  * Ned_k := (\P_{n,k-1})^n \oplus \{p \in (\P_{n,k})^n: <p,x>=0 \}
74  * \end{equation*}
75  * with $x=(x,y)$ in two dimensions and $x=(x,y,z)$ in three dimensions.
76  *
77  * For $Ned_k$ holds
78  * \begin{equation*}
79  * (\P_{n,k-1})^n \subset Ned_k \subset (\P_{n,k})^n.
80  * \end{equation*}
81  *
82  * We construct $(\P_{n,k})^n$ and and only use the monomials contained in $Ned$.
83  *
84  */
85  if( dim!=2 && dim!=3 || !geometry.isSimplex())
86  DUNE_THROW(Dune::NotImplemented,"High order nedelec elements are only supported by simplices in 2d and 3d");
87 
88  MIBasis basis(order+1);
89  FieldVector< MI, dim > x;
90  /*
91  * Init MultiIndices
92  * x[0]=(1,0,0) x
93  * x[1]=(0,1,0) y
94  * x[2]=(0,0,1) z
95  */
96  for( unsigned int i = 0; i < dim; ++i )
97  x[i].set(i,1);
98  std::vector< MI > val( basis.size() );
99 
100  // val now contains all monomials in $n$ dimensions with degree $\leq order+1$
101  basis.evaluate( x, val );
102 
103  col_ = basis.size();
104 
105  // get $\dim (\P_{n,order-1})$
106  unsigned int notHomogen = 0;
107  if (order>0)
108  notHomogen = basis.sizes()[order-1];
109 
110  // the number of basis functions for the set of homogeneous polynomials with degree $order$
111  unsigned int homogen = basis.sizes()[order]-notHomogen;
112 
113  /*
114  * 2D:
115  * \begin{equation*}
116  * Ned_{order} = (\P_{order-1})^2 \oplus (-y,x)^T \widetilde \P_{order-1}
117  * \end{equation*}
118  *
119  * It gets more complicated in higher dimensions.
120  *
121  * 3D:
122  * \begin{equation*}
123  * Ned_{order} = (\P_{n,order-1})^3 \oplus (z,0,-x)^T \widetilde \P_{n,order-1} \oplus (-y,x,0)^T \widetilde \P_{n,order-1} \oplus (0,-z,y)^T \widetilde \P_{n-1,order-1}
124  * \end{equation*}
125  *
126  * Note the last term. The index $n-1$ is on purpose.
127  * Else i.e. k=2
128  *
129  * (0,z,-y)^T x = (z,0,-x)^T y - (y,-x,0)^T z.
130  *
131  */
132 
133  /*
134  * compute the number of rows for the coefficient matrix
135  *
136  * row_ = dim* \dim Ned_{order}
137  */
138  if (dim == 2)
139  row_ = (notHomogen*dim+homogen*(dim+1))*dim;
140  else if (dim==3)
141  {
142  // get dim \P_{n-1,order-1}
143  int homogenTwoVariables = 0;
144  for( int w = notHomogen; w<notHomogen + homogen; w++)
145  if (val[w].z(0)==0)
146  homogenTwoVariables++;
147  row_ = (notHomogen*dim+homogen*(dim+2) + homogenTwoVariables)*dim;
148  }
149 
150  mat_ = new Field*[row_];
151  int row = 0;
152 
153  /* Assign the correct values for the nonhomogeneous polymonials $p\in (\P_{n,order-1})^dim$
154  * A basis function is represented by $dim$ rows.
155  */
156  for (unsigned int i=0; i<notHomogen+homogen; ++i)
157  {
158  for (unsigned int r=0; r<dim; ++r)
159  {
160  for (unsigned int rr=0; rr<dim; ++rr)
161  {
162  // init row to zero
163  mat_[row] = new Field[col_];
164  for (unsigned int j=0; j<col_; ++j)
165  mat_[row][j] = 0.;
166 
167  if (r==rr)
168  mat_[row][i] = 1.;
169  ++row;
170  }
171  }
172  }
173 
174  /* Assign the correct values for the homogeneous polymonials $p\in Ned_{order} \backslash (\P_{n,order-1})^dim$
175  * A basis function is represented by $dim$ rows.
176  */
177  for (unsigned int i=0; i<homogen; ++i)
178  {
179  // get a monomial $ p \in \P_{n,order}\backslash \P_{n,order-1}$
180  MI xval = val[notHomogen+i];
181  if(dim==2)
182  {
183  for (unsigned int r=0; r<dim; ++r)
184  {
185  // init rows to zero
186  mat_[row+r] = new Field[col_];
187  for (unsigned int j=0; j<col_; ++j)
188  mat_[row+r][j] = 0.;
189  }
190 
191  /* set $(-y,x)^T p$ with a homogeneous monomial $p$
192  *
193  * The loop over the monomials is needed to obtain the corresponding column index.
194  */
195  for (int w=homogen+notHomogen; w<val.size(); ++w)
196  {
197  if (val[w] == xval*x[0])
198  mat_[row+1][w] = 1.;
199  if (val[w] == xval*x[1])
200  mat_[row][w] = -1.;
201  }
202  row +=dim;
203  }
204  else if(dim==3)
205  {
206  int skipLastDim = xval.z(0)>0;
207  /*
208  * Init 9 rows to zero.
209  * If the homogeneous monomial has a positive x-exponent (0,-z,y) gets skipped ( see example for the Nedelec space in 3D )
210  * In this case only 6 rows get initialised.
211  */
212  for (unsigned int r=0; r<dim*(dim-skipLastDim); ++r)
213  {
214  // init rows to zero
215  mat_[row+r] = new Field[col_];
216  for (unsigned int j=0; j<col_; ++j)
217  mat_[row+r][j] = 0.;
218  }
219 
220  /*
221  * first $dim$ rows are for (z,0,-x)
222  *
223  * second $dim$ rows are for (-y,x,0)
224  *
225  * third $dim$ rows are for (0,-z,y)
226  *
227  */
228  for (unsigned int r=0; r<dim - skipLastDim; ++r)
229  {
230  int index = (r+dim-1)%dim;
231  for (int w=homogen+notHomogen; w<val.size(); ++w)
232  {
233  if (val[w] == xval*x[index])
234  mat_[row+r][w] = 1.;
235  if (val[w] == xval*x[r])
236  mat_[row+index][w] = -1.;
237  }
238  row +=dim;
239  }
240 
241  }
242  }
243  }
244 
246  {
247  for (unsigned int i=0; i<rows(); ++i) {
248  delete [] mat_[i];
249  }
250  delete [] mat_;
251  }
252 
253  unsigned int cols() const {
254  return col_;
255  }
256 
257  unsigned int rows() const {
258  return row_;
259  }
260 
261  template <class Vector>
262  void row( const unsigned int row, Vector &vec ) const
263  {
264  const unsigned int N = cols();
265  assert( vec.size() == N );
266  for (unsigned int i=0; i<N; ++i)
267  field_cast(mat_[row][i],vec[i]);
268  }
269 
270  unsigned int row_,col_;
271  Field **mat_;
272  };
273 
274 
275 }
276 #endif // #ifndef DUNE_LOCALFUNCTIONS_NEDELEC_NEDELECSIMPLEX_NEDELECSIMPLEXPREBASIS_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: nedelecsimplexprebasis.hh:59
NedelecVecMatrix(std::size_t order)
Definition: nedelecsimplexprebasis.hh:64
MultiIndex< dim, Field > MI
Definition: nedelecsimplexprebasis.hh:62
unsigned int row_
Definition: nedelecsimplexprebasis.hh:270
unsigned int cols() const
Definition: nedelecsimplexprebasis.hh:253
~NedelecVecMatrix()
Definition: nedelecsimplexprebasis.hh:245
MonomialBasis< geometryId, MI > MIBasis
Definition: nedelecsimplexprebasis.hh:63
unsigned int col_
Definition: nedelecsimplexprebasis.hh:270
static const unsigned int dim
Definition: nedelecsimplexprebasis.hh:61
void row(const unsigned int row, Vector &vec) const
Definition: nedelecsimplexprebasis.hh:262
static constexpr GeometryType geometry
Definition: nedelecsimplexprebasis.hh:60
unsigned int rows() const
Definition: nedelecsimplexprebasis.hh:257
Field ** mat_
Definition: nedelecsimplexprebasis.hh:271
Definition: nedelecsimplexprebasis.hh:20
static void release(Object *object)
Definition: nedelecsimplexprebasis.hh:54
MBasisFactory::Object MBasis
Definition: nedelecsimplexprebasis.hh:22
static Object * create(Key order)
Definition: nedelecsimplexprebasis.hh:36
PolynomialBasisWithMatrix< EvalMBasis, SparseCoeffMatrix< Field, dim > > Basis
Definition: nedelecsimplexprebasis.hh:24
const Basis Object
Definition: nedelecsimplexprebasis.hh:26
StandardEvaluator< MBasis > EvalMBasis
Definition: nedelecsimplexprebasis.hh:23
MonomialBasisProvider< dim, Field > MBasisFactory
Definition: nedelecsimplexprebasis.hh:21
std::size_t Key
Definition: nedelecsimplexprebasis.hh:27
Definition: nedelecsimplexprebasis.hh:31
MonomialBasisProvider< dd, FF > Type
Definition: nedelecsimplexprebasis.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
int z(int i) const
Definition: multiindex.hh:89
Definition: polynomialbasis.hh:335