DUNE-FEM (unstable)

nedelec1stkindcube.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_NEDELEC1STKINDCUBE_HH
6#define DUNE_LOCALFUNCTIONS_NEDELEC_NEDELEC1STKINDCUBE_HH
7
8#include <numeric>
9
12#include <dune/common/math.hh>
13
14#include <dune/geometry/referenceelements.hh>
15#include <dune/geometry/type.hh>
16
17#include <dune/localfunctions/common/localbasis.hh>
18#include <dune/localfunctions/common/localfiniteelementtraits.hh>
19#include <dune/localfunctions/common/localkey.hh>
20
21namespace Dune
22{
23namespace Impl
24{
37 template<class D, class R, int dim, int k>
38 class Nedelec1stKindCubeLocalBasis
39 {
40 // Number of edges of the reference cube
41 constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim;
42
43 public:
44 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,
45 R,dim,FieldVector<R,dim>,
46 FieldMatrix<R,dim,dim> >;
47
54 Nedelec1stKindCubeLocalBasis()
55 {
56 std::fill(edgeOrientation_.begin(), edgeOrientation_.end(), 1.0);
57 }
58
61 Nedelec1stKindCubeLocalBasis(std::bitset<numberOfEdges> edgeOrientation)
62 : Nedelec1stKindCubeLocalBasis()
63 {
64 for (std::size_t i=0; i<edgeOrientation_.size(); i++)
65 edgeOrientation_[i] *= edgeOrientation[i] ? -1.0 : 1.0;
66 }
67
69 static constexpr unsigned int size()
70 {
71 static_assert(dim==2 || dim==3, "Nedelec shape functions are implemented only for 2d and 3d cubes.");
72 if (dim==2)
73 return 2*k * (k+1);
74 if (dim==3)
75 return 3*k * (k+1) * (k+1);
76 }
77
83 void evaluateFunction(const typename Traits::DomainType& in,
84 std::vector<typename Traits::RangeType>& out) const
85 {
86 static_assert(k==1, "Evaluating Nédélec shape functions is implemented only for first order.");
87 out.resize(size());
88
89 if (dim==2)
90 {
91 // First-order Nédélec shape functions on a square are of the form
92 //
93 // (a, b)^T + (c y, d x)^T, a, b, c, d \in R
94 //
95 // The following coefficients create the four basis vectors
96 // that are dual to the edge degrees of freedom:
97 //
98 // a[0] = 0 b[0] = 1 c[0] = 0 d[0] = -1
99 // a[1] = 0 b[1] = 0 c[1] = 0 d[1] = 1
100 // a[2] = 1 b[2] = 0 c[2] = 0 d[2] = -1
101 // a[3] = 0 b[3] = 0 c[3] = 0 d[3] = 1
102
103 out[0] = { 0, D(1) - in[0]};
104 out[1] = { 0, in[0]};
105 out[2] = { D(1) - in[1], 0};
106 out[3] = { in[1], 0};
107 }
108
109 if constexpr (dim==3)
110 {
111 // First-order Nédélec shape functions on a cube are of the form
112 //
113 // (e1 yz)
114 // a + b x + c y + d z + (e2 xz) , a, b, c, d \in R^3 and b[0]=c[1]=d[2]=0
115 // (e3 xy)
116 //
117 // The following coefficients create the twelve basis vectors
118 // that are dual to the edge degrees of freedom:
119 //
120 // a[0] = { 0, 0, 1} b[0] = { 0, 0, -1} c[0] = { 0, 0, -1} d[0] = { 0, 0, 0} e[0] = { 0, 0, 1}
121 // a[1] = { 0, 0, 0} b[1] = { 0, 0, 1} c[1] = { 0, 0, 0} d[1] = { 0, 0, 0} e[1] = { 0, 0, -1}
122 // a[2] = { 0, 0, 0} b[2] = { 0, 0, 0} c[2] = { 0, 0, 1} d[2] = { 0, 0, 0} e[2] = { 0, 0, -1}
123 // a[3] = { 0, 0, 0} b[3] = { 0, 0, 0} c[3] = { 0, 0, 0} d[3] = { 0, 0, 0} e[3] = { 0, 0, 1}
124 //
125 // The following implementation uses these values, and simply
126 // skips all the zeros.
127
128 for (std::size_t i=0; i<out.size(); i++)
129 out[i] = {0,0,0};
130
131 out[0][2] = { 1 - in[0] - in[1] + in[0]*in[1]};
132 out[1][2] = { in[0] - in[0]*in[1]};
133 out[2][2] = { in[1] - in[0]*in[1]};
134 out[3][2] = { in[0]*in[1]};
135
136 out[4][1] = { 1 - in[0] - in[2] + in[0]*in[2]};
137 out[5][1] = { in[0] - in[0]*in[2]};
138 out[8][1] = { in[2] - in[0]*in[2]};
139 out[9][1] = { in[0]*in[2]};
140
141 out[6][0] = { 1 - in[1] - in[2] + in[1]*in[2]};
142 out[7][0] = { in[1] - in[1]*in[2]};
143 out[10][0] = { in[2] - in[1]*in[2]};
144 out[11][0] = { in[1]*in[2]};
145 }
146
147 for (std::size_t i=0; i<out.size(); i++)
148 out[i] *= edgeOrientation_[i];
149 }
150
156 void evaluateJacobian(const typename Traits::DomainType& in,
157 std::vector<typename Traits::JacobianType>& out) const
158 {
159 out.resize(size());
160 if (dim==2)
161 {
162 for (std::size_t i=0; i<out.size(); i++)
163 for (std::size_t j=0; j<dim; j++)
164 out[i][j] = { 0, 0};
165
166 out[0][1] = { -1, 0};
167 out[1][1] = { 1, 0};
168
169 out[2][0] = { 0, -1};
170 out[3][0] = { 0, 1};
171 }
172 if (dim==3)
173 {
174 for (std::size_t i=0; i<out.size(); i++)
175 for(std::size_t j=0;j<dim; j++)
176 out[i][j] = {0,0,0};
177
178
179 out[0][2] = {-1 +in[1], -1 + in[0], 0};
180 out[1][2] = { 1 -in[1], - in[0], 0};
181 out[2][2] = { -in[1], 1 - in[0], 0};
182 out[3][2] = { in[1], in[0], 0};
183
184 out[4][1] = {-1 +in[2], 0, -1 + in[0]};
185 out[5][1] = { 1 -in[2], 0, - in[0]};
186 out[8][1] = { -in[2], 0, 1 - in[0]};
187 out[9][1] = { in[2], 0, in[0]};
188
189 out[6][0] = { 0, -1 + in[2], -1 + in[1]};
190 out[7][0] = { 0, 1 - in[2], - in[1]};
191 out[10][0] = { 0, - in[2], 1 - in[1]};
192 out[11][0] = { 0, in[2], in[1]};
193
194 }
195
196 for (std::size_t i=0; i<out.size(); i++)
197 out[i] *= edgeOrientation_[i];
198
199 }
200
207 void partial(const std::array<unsigned int, dim>& order,
208 const typename Traits::DomainType& in,
209 std::vector<typename Traits::RangeType>& out) const
210 {
211 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
212 if (totalOrder == 0) {
213 evaluateFunction(in, out);
214 } else if (totalOrder == 1) {
215 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
216 out.resize(size());
217
218 if (dim==2)
219 {
220 if (direction==0)
221 {
222 out[0] = { 0, -1};
223 out[1] = { 0, 1};
224 out[2] = { 0, 0};
225 out[3] = { 0, 0};
226 }
227 else
228 {
229 out[0] = { 0, 0};
230 out[1] = { 0, 0};
231 out[2] = { -1, 0};
232 out[3] = { 1, 0};
233 }
234 }
235
236 if (dim==3)
237 {
238 switch (direction)
239 {
240 case 0:
241 out[0] = { 0, 0, -1 +in[1]};
242 out[1] = { 0, 0, 1 -in[1]};
243 out[2] = { 0, 0, -in[1]};
244 out[3] = { 0, 0, in[1]};
245
246 out[4] = { 0, -1 +in[2], 0};
247 out[5] = { 0, 1 -in[2], 0};
248 out[8] = { 0, -in[2], 0};
249 out[9] = { 0, in[2], 0};
250
251 out[6] = {0,0,0};
252 out[7] = {0,0,0};
253 out[10] = {0,0,0};
254 out[11] = {0,0,0};
255 break;
256
257 case 1:
258 out[0] = { 0, 0, -1 + in[0]};
259 out[1] = { 0, 0, - in[0]};
260 out[2] = { 0, 0, 1 - in[0]};
261 out[3] = { 0, 0, in[0]};
262
263 out[4] = {0,0,0};
264 out[5] = {0,0,0};
265 out[8] = {0,0,0};
266 out[9] = {0,0,0};
267
268 out[6] = { -1 + in[2], 0, 0};
269 out[7] = { 1 - in[2], 0, 0};
270 out[10] = { - in[2], 0, 0};
271 out[11] = { in[2], 0, 0};
272 break;
273
274 case 2:
275 out[0] = {0,0,0};
276 out[1] = {0,0,0};
277 out[2] = {0,0,0};
278 out[3] = {0,0,0};
279
280 out[4] = { 0, -1 + in[0], 0};
281 out[5] = { 0, - in[0], 0};
282 out[8] = { 0, 1 - in[0], 0};
283 out[9] = { 0, in[0], 0};
284
285 out[6] = { -1 + in[1], 0, 0};
286 out[7] = { - in[1], 0, 0};
287 out[10] = { 1 - in[1], 0, 0};
288 out[11] = { in[1], 0, 0};
289 break;
290 }
291 }
292
293 for (std::size_t i=0; i<out.size(); i++)
294 out[i] *= edgeOrientation_[i];
295
296 } else if (totalOrder == 2) {
297 out.resize(size());
298
299 if (dim==2)
300 for (std::size_t i = 0; i < size(); ++i)
301 for (std::size_t j = 0; j < dim; ++j)
302 out[i][j] = 0;
303
304 if (dim==3)
305 {
306 for(size_t i=0; i<out.size(); i++)
307 out[i] = { 0, 0, 0};
308
309 //case (1,1,0):
310 if( order[0] == 1 and order[1]==1)
311 {
312 out[0] = { 0, 0, 1};
313 out[1] = { 0, 0, -1};
314 out[2] = { 0, 0, -1};
315 out[3] = { 0, 0, 1};
316 }
317
318 //case (1,0,1):
319 if( order[0] == 1 and order[2]==1)
320 {
321 out[4] = { 0, 1, 0};
322 out[5] = { 0, -1, 0};
323 out[8] = { 0, -1, 0};
324 out[9] = { 0, 1, 0};
325 }
326
327 //case (0,1,1):
328 if( order[1] == 1 and order[2]==1)
329 {
330 out[6] = { 1, 0, 0};
331 out[7] = { -1, 0, 0};
332 out[10] = { -1, 0, 0};
333 out[11] = { 1, 0, 0};
334 }
335
336 for (std::size_t i=0; i<out.size(); i++)
337 out[i] *= edgeOrientation_[i];
338 }
339
340
341 }else {
342 out.resize(size());
343 for (std::size_t i = 0; i < size(); ++i)
344 for (std::size_t j = 0; j < dim; ++j)
345 out[i][j] = 0;
346 }
347
348 }
349
351 unsigned int order() const
352 {
353 if (dim==2)
354 return 2*k-1;
355 if (dim==3)
356 return 3*k-1;
357 }
358
359 private:
360
361 // Orientations of the cube edges
362 std::array<R,numberOfEdges> edgeOrientation_;
363 };
364
365
370 template <int dim, int k>
371 class Nedelec1stKindCubeLocalCoefficients
372 {
373 public:
375 Nedelec1stKindCubeLocalCoefficients ()
376 : localKey_(size())
377 {
378 static_assert(k==1, "Only first-order Nédélec local coefficients are implemented.");
379 // Assign all degrees of freedom to edges
380 // TODO: This is correct only for first-order Nédélec elements
381 for (std::size_t i=0; i<size(); i++)
382 localKey_[i] = LocalKey(i,dim-1,0);
383 }
384
386 std::size_t size() const
387 {
388 static_assert(dim==2 || dim==3, "Nédélec shape functions are implemented only for 2d and 3d cubes.");
389 return (dim==2) ? 2*k * (k+1)
390 : 3*k * (k+1) * (k+1);
391 }
392
395 const LocalKey& localKey (std::size_t i) const
396 {
397 return localKey_[i];
398 }
399
400 private:
401 std::vector<LocalKey> localKey_;
402 };
403
408 template<class LB>
409 class Nedelec1stKindCubeLocalInterpolation
410 {
411 static constexpr auto dim = LB::Traits::dimDomain;
412 static constexpr auto size = LB::size();
413
414 // Number of edges of the reference cube
415 constexpr static std::size_t numberOfEdges = power(2,dim-1)*dim;
416
417 public:
418
420 Nedelec1stKindCubeLocalInterpolation (std::bitset<numberOfEdges> s = 0)
421 {
422 auto refElement = Dune::referenceElement<double,dim>(GeometryTypes::cube(dim));
423
424 for (std::size_t i=0; i<numberOfEdges; i++)
425 m_[i] = refElement.position(i,dim-1);
426
427 for (std::size_t i=0; i<numberOfEdges; i++)
428 {
429 auto vertexIterator = refElement.subEntities(i,dim-1,dim).begin();
430 auto v0 = *vertexIterator;
431 auto v1 = *(++vertexIterator);
432
433 // By default, edges point from the vertex with the smaller index
434 // to the vertex with the larger index.
435 if (v0>v1)
436 std::swap(v0,v1);
437 edge_[i] = refElement.position(v1,dim) - refElement.position(v0,dim);
438 edge_[i] *= (s[i]) ? -1.0 : 1.0;
439 }
440 }
441
447 template<typename F, typename C>
448 void interpolate (const F& f, std::vector<C>& out) const
449 {
450 out.resize(size);
451
452 for (std::size_t i=0; i<size; i++)
453 {
454 auto y = f(m_[i]);
455 out[i] = 0.0;
456 for (int j=0; j<dim; j++)
457 out[i] += y[j]*edge_[i][j];
458 }
459 }
460
461 private:
462 // Edge midpoints of the reference cube
463 std::array<typename LB::Traits::DomainType, numberOfEdges> m_;
464 // Edges of the reference cube
465 std::array<typename LB::Traits::DomainType, numberOfEdges> edge_;
466 };
467
468}
469
470
494 template<class D, class R, int dim, int k>
496 {
497 public:
499 Impl::Nedelec1stKindCubeLocalCoefficients<dim,k>,
500 Impl::Nedelec1stKindCubeLocalInterpolation<Impl::Nedelec1stKindCubeLocalBasis<D,R,dim,k> > >;
501
502 static_assert(dim==2 || dim==3, "Nedelec elements are only implemented for 2d and 3d elements.");
503 static_assert(k==1, "Nedelec elements of the first kind are currently only implemented for order k==1.");
504
508
514 Nedelec1stKindCubeLocalFiniteElement (std::bitset<power(2,dim-1)*dim> s) :
515 basis_(s),
516 interpolation_(s)
517 {}
518
519 const typename Traits::LocalBasisType& localBasis () const
520 {
521 return basis_;
522 }
523
524 const typename Traits::LocalCoefficientsType& localCoefficients () const
525 {
526 return coefficients_;
527 }
528
529 const typename Traits::LocalInterpolationType& localInterpolation () const
530 {
531 return interpolation_;
532 }
533
534 static constexpr unsigned int size ()
535 {
536 return Traits::LocalBasisType::size();
537 }
538
539 static constexpr GeometryType type ()
540 {
541 return GeometryTypes::cube(dim);
542 }
543
544 private:
545 typename Traits::LocalBasisType basis_;
546 typename Traits::LocalCoefficientsType coefficients_;
547 typename Traits::LocalInterpolationType interpolation_;
548 };
549
550}
551
552#endif
Nédélec elements of the first kind for cube elements.
Definition: nedelec1stkindcube.hh:496
Nedelec1stKindCubeLocalFiniteElement()=default
Default constructor.
Nedelec1stKindCubeLocalFiniteElement(std::bitset< power(2, dim-1) *dim > s)
Constructor with explicitly given edge orientations.
Definition: nedelec1stkindcube.hh:514
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 cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Some useful basic math stuff.
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
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.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 12, 23:30, 2024)