DUNE-FEM (unstable)

nedelec1stkindsimplex.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_NEDELEC1STKINDSIMPLEX_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDSIMPLEX_HH
7
8#include <numeric>
9
12
13#include <dune/geometry/referenceelements.hh>
14#include <dune/geometry/type.hh>
15
16#include <dune/localfunctions/common/localbasis.hh>
17#include <dune/localfunctions/common/localfiniteelementtraits.hh>
18#include <dune/localfunctions/common/localkey.hh>
19
20namespace Dune
21{
22namespace Impl
23{
36 template<class D, class R, int dim, int k>
37 class Nedelec1stKindSimplexLocalBasis
38 {
39 // Number of edges of the reference simplex
40 constexpr static std::size_t numberOfEdges = dim*(dim+1)/2;
41
42 public:
43 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,
44 R,dim,FieldVector<R,dim>,
45 FieldMatrix<R,dim,dim> >;
46
53 Nedelec1stKindSimplexLocalBasis()
54 {
55 std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0);
56 }
57
60 Nedelec1stKindSimplexLocalBasis(std::bitset<numberOfEdges> edgeOrientation)
61 : Nedelec1stKindSimplexLocalBasis()
62 {
63 for (std::size_t i=0; i<edgeOrientation_.size(); i++)
64 edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
65 }
66
68 static constexpr unsigned int size()
69 {
70 static_assert(dim==2 || dim==3, "Nedelec shape functions are implemented only for 2d and 3d simplices.");
71 if (dim==2)
72 return k * (k+2);
73 if (dim==3)
74 return k * (k+2) * (k+3) / 2;
75 }
76
82 void evaluateFunction(const typename Traits::DomainType& in,
83 std::vector<typename Traits::RangeType>& out) const
84 {
85 static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order.");
86 out.resize(size());
87
88 if (dim==2)
89 {
90 // First-order Nédélec shape functions on a triangle are of the form
91 //
92 // (a1, a2) + b(-x2, x1)^T, a_1, a_2, b \in R
93 out[0] = {D(1) - in[1], in[0]};
94 out[1] = {in[1], -in[0]+D(1)};
95 out[2] = {-in[1], in[0]};
96 }
97
98 if constexpr (dim==3)
99 {
100 // First-order Nédélec shape functions on a tetrahedron are of the form
101 //
102 // a + b \times x, a, b \in R^3
103 //
104 // The following coefficients create the six basis vectors
105 // that are dual to the edge degrees of freedom:
106 //
107 // a[0] = { 1, 0, 0} b[0] = { 0, -1, 1}
108 // a[1] = { 0, 1, 0} b[1] = { 1, 0, -1}
109 // a[2] = { 0, 0, 0} b[2] = { 0, 0, 1}
110 // a[3] = { 0, 0, 1} b[3] = {-1, 1, 0}
111 // a[4] = { 0, 0, 0} b[4] = { 0, -1, 0}
112 // a[5] = { 0, 0, 0} b[5] = { 1, 0, 0}
113 //
114 // The following implementation uses these values, and simply
115 // skips all the zeros.
116
117 out[0] = { 1 - in[1] - in[2], in[0] , in[0] };
118 out[1] = { in[1] , 1 - in[0] - in[2], in[1]};
119 out[2] = { - in[1] , in[0] , 0 };
120 out[3] = { in[2], in[2], 1 - in[0] - in[1]};
121 out[4] = { -in[2], 0 , in[0] };
122 out[5] = { 0 , -in[2], in[1]};
123 }
124
125 for (std::size_t i=0; i<out.size(); i++)
126 out[i] *= edgeOrientation_[i];
127 }
128
134 void evaluateJacobian(const typename Traits::DomainType& in,
135 std::vector<typename Traits::JacobianType>& out) const
136 {
137 out.resize(size());
138 if (dim==2)
139 {
140 out[0][0] = { 0, -1};
141 out[0][1] = { 1, 0};
142
143 out[1][0] = { 0, 1};
144 out[1][1] = {-1, 0};
145
146 out[2][0] = { 0, -1};
147 out[2][1] = { 1, 0};
148 }
149 if (dim==3)
150 {
151 out[0][0] = { 0,-1,-1};
152 out[0][1] = { 1, 0, 0};
153 out[0][2] = { 1, 0, 0};
154
155 out[1][0] = { 0, 1, 0};
156 out[1][1] = {-1, 0, -1};
157 out[1][2] = { 0, 1, 0};
158
159 out[2][0] = { 0, -1, 0};
160 out[2][1] = { 1, 0, 0};
161 out[2][2] = { 0, 0, 0};
162
163 out[3][0] = { 0, 0, 1};
164 out[3][1] = { 0, 0, 1};
165 out[3][2] = {-1, -1, 0};
166
167 out[4][0] = { 0, 0, -1};
168 out[4][1] = { 0, 0, 0};
169 out[4][2] = { 1, 0, 0};
170
171 out[5][0] = { 0, 0, 0};
172 out[5][1] = { 0, 0, -1};
173 out[5][2] = { 0, 1, 0};
174 }
175
176 for (std::size_t i=0; i<out.size(); i++)
177 out[i] *= edgeOrientation_[i];
178
179 }
180
187 void partial(const std::array<unsigned int, dim>& order,
188 const typename Traits::DomainType& in,
189 std::vector<typename Traits::RangeType>& out) const
190 {
191 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
192 if (totalOrder == 0) {
193 evaluateFunction(in, out);
194 } else if (totalOrder == 1) {
195 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
196 out.resize(size());
197
198 if (dim==2)
199 {
200 if (direction==0)
201 {
202 out[0] = {0, 1};
203 out[1] = {0, -1};
204 out[2] = {0, 1};
205 }
206 else
207 {
208 out[0] = {-1, 0};
209 out[1] = { 1, 0};
210 out[2] = {-1, 0};
211 }
212 }
213
214 if (dim==3)
215 {
216 switch (direction)
217 {
218 case 0:
219 out[0] = { 0, 1, 1};
220 out[1] = { 0,-1, 0};
221 out[2] = { 0, 1, 0};
222 out[3] = { 0, 0,-1};
223 out[4] = { 0, 0, 1};
224 out[5] = { 0, 0, 0};
225 break;
226
227 case 1:
228 out[0] = {-1, 0, 0};
229 out[1] = { 1, 0, 1};
230 out[2] = {-1, 0, 0};
231 out[3] = { 0, 0,-1};
232 out[4] = { 0, 0, 0};
233 out[5] = { 0, 0, 1};
234 break;
235
236 case 2:
237 out[0] = {-1, 0, 0};
238 out[1] = { 0,-1, 0};
239 out[2] = { 0, 0, 0};
240 out[3] = { 1, 1, 0};
241 out[4] = {-1, 0, 0};
242 out[5] = { 0,-1, 0};
243 break;
244 }
245 }
246
247 for (std::size_t i=0; i<out.size(); i++)
248 out[i] *= edgeOrientation_[i];
249
250 } else {
251 out.resize(size());
252 for (std::size_t i = 0; i < size(); ++i)
253 for (std::size_t j = 0; j < dim; ++j)
254 out[i][j] = 0;
255 }
256
257 }
258
260 unsigned int order() const
261 {
262 return k;
263 }
264
265 private:
266
267 // Orientations of the simplex edges
268 std::array<R,numberOfEdges> edgeOrientation_;
269 };
270
271
276 template <int dim, int k>
277 class Nedelec1stKindSimplexLocalCoefficients
278 {
279 public:
281 Nedelec1stKindSimplexLocalCoefficients ()
282 : localKey_(size())
283 {
284 static_assert(k==1, "Only first-order Nédélec local coefficients are implemented.");
285 // Assign all degrees of freedom to edges
286 // TODO: This is correct only for first-order Nédélec elements
287 for (std::size_t i=0; i<size(); i++)
288 localKey_[i] = LocalKey(i,dim-1,0);
289 }
290
292 std::size_t size() const
293 {
294 static_assert(dim==2 || dim==3, "Nédélec shape functions are implemented only for 2d and 3d simplices.");
295 return (dim==2) ? k * (k+2)
296 : k * (k+2) * (k+3) / 2;
297 }
298
301 const LocalKey& localKey (std::size_t i) const
302 {
303 return localKey_[i];
304 }
305
306 private:
307 std::vector<LocalKey> localKey_;
308 };
309
314 template<class LB>
315 class Nedelec1stKindSimplexLocalInterpolation
316 {
317 static constexpr auto dim = LB::Traits::dimDomain;
318 static constexpr auto size = LB::size();
319
320 // Number of edges of the reference simplex
321 constexpr static std::size_t numberOfEdges = dim*(dim+1)/2;
322
323 public:
324
326 Nedelec1stKindSimplexLocalInterpolation (std::bitset<numberOfEdges> s = 0)
327 {
328 auto refElement = Dune::referenceElement<double,dim>(GeometryTypes::simplex(dim));
329
330 for (std::size_t i=0; i<numberOfEdges; i++)
331 m_[i] = refElement.position(i,dim-1);
332
333 for (std::size_t i=0; i<numberOfEdges; i++)
334 {
335 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
336 auto v0 = *vertexIterator;
337 auto v1 = *(++vertexIterator);
338 // By default, edges point from the vertex with the smaller index
339 // to the vertex with the larger index.
340 if (v0>v1)
341 std::swap(v0,v1);
342 edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim);
343 edge_[i] *= (s[i]) ? -1.0 : 1.0;
344 }
345 }
346
352 template<typename F, typename C>
353 void interpolate (const F& f, std::vector<C>& out) const
354 {
355 out.resize(size);
356
357 for (std::size_t i=0; i<size; i++)
358 {
359 auto y = f(m_[i]);
360 out[i] = 0.0;
361 for (int j=0; j<dim; j++)
362 out[i] += y[j]*edge_[i][j];
363 }
364 }
365
366 private:
367 // Edge midpoints of the reference simplex
368 std::array<typename LB::Traits::DomainType, numberOfEdges> m_;
369 // Edges of the reference simplex
370 std::array<typename LB::Traits::DomainType, numberOfEdges> edge_;
371 };
372
373}
374
375
401 template<class D, class R, int dim, int k>
403 {
404 public:
406 Impl::Nedelec1stKindSimplexLocalCoefficients<dim,k>,
407 Impl::Nedelec1stKindSimplexLocalInterpolation<Impl::Nedelec1stKindSimplexLocalBasis<D,R,dim,k> > >;
408
409 static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements.");
410 static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1.");
411
415
421 Nedelec1stKindSimplexLocalFiniteElement (std::bitset<dim*(dim+1)/2> s) :
422 basis_(s),
423 interpolation_(s)
424 {}
425
426 const typename Traits::LocalBasisType& localBasis () const
427 {
428 return basis_;
429 }
430
431 const typename Traits::LocalCoefficientsType& localCoefficients () const
432 {
433 return coefficients_;
434 }
435
436 const typename Traits::LocalInterpolationType& localInterpolation () const
437 {
438 return interpolation_;
439 }
440
441 static constexpr unsigned int size ()
442 {
443 return Traits::LocalBasisType::size();
444 }
445
446 static constexpr GeometryType type ()
447 {
448 return GeometryTypes::simplex(dim);
449 }
450
451 private:
452 typename Traits::LocalBasisType basis_;
453 typename Traits::LocalCoefficientsType coefficients_;
454 typename Traits::LocalInterpolationType interpolation_;
455 };
456
457}
458
459#endif
Nédélec elements of the first kind for simplex elements.
Definition: nedelec1stkindsimplex.hh:403
Nedelec1stKindSimplexLocalFiniteElement()=default
Default constructor.
Nedelec1stKindSimplexLocalFiniteElement(std::bitset< dim *(dim+1)/2 > s)
Constructor with explicitly given edge orientations.
Definition: nedelec1stkindsimplex.hh:421
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
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 (Nov 21, 23:30, 2024)