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