Dune Core Modules (unstable)

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