DUNE-FEM (unstable)

localbasis.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_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH
6#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH
7
8#include <algorithm>
9#include <array>
10#include <bitset>
11#include <numeric>
12#include <vector>
13
16#include <dune/common/math.hh>
19
20#include <dune/localfunctions/common/localbasis.hh>
21
22namespace Dune
23{
35 template<class D, class R, unsigned int dim, unsigned int order>
37 {
38 static_assert( AlwaysFalse<D>::value,
39 "`BDFMCubeLocalBasis` not implemented for chosen `dim` and `order`." );
40 };
41
42
43#ifndef DOXYGEN
44 template<class D, class R, unsigned int dim>
45 class BDFMCubeLocalBasis<D, R, dim, 0>
46 {
47 static_assert( AlwaysFalse<D>::value,
48 "`BDFMCubeLocalBasis` not defined for order 0." );
49 };
50#endif // #ifndef DOXYGEN
51
52
58 template<class D, class R >
59 class BDFMCubeLocalBasis<D, R, 2, 1>
60 {
64
65 public:
67
70 {
71 std::fill(s_.begin(), s_.end(), 1);
72 }
73
79 BDFMCubeLocalBasis (std::bitset<4> s)
80 {
81 for (auto i : range(4))
82 s_[i] = s[i] ? -1 : 1;
83 }
84
86 unsigned int size () const { return 4; }
87
94 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const
95 {
96 out.resize(4);
97
98 out[0] = {s_[0]*(in[0]-1), 0 };
99 out[1] = {s_[1]*(in[0]) , 0 };
100 out[2] = {0, s_[2]*(in[1]-1)};
101 out[3] = {0, s_[3]*(in[1]) };
102 }
103
110 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const
111 {
112 out.resize(4);
113
114 out[0] = {{s_[0], 0}, {0, 0}};
115 out[1] = {{s_[1], 0}, {0, 0}};
116 out[2] = {{0, 0}, {0, s_[2]}};
117 out[3] = {{0, 0}, {0, s_[3]}};
118 }
119
127 void partial (const std::array<unsigned int, 2>& order,
128 const DomainType& in,
129 std::vector<RangeType>& out) const
130 {
131 if (std::accumulate(order.begin(), order.end(), 0) == 0)
132 evaluateFunction(in, out);
133 else
134 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
135 }
136
138 unsigned int order () const { return 1; }
139
140 private:
141 std::array<R, 4> s_;
142 };
143
144
150 template<class D, class R >
151 class BDFMCubeLocalBasis<D, R, 2, 2>
152 {
156
157 public:
159
162 {
163 std::fill(s_.begin(), s_.end(), 1);
164 }
165
171 BDFMCubeLocalBasis (std::bitset<4> s)
172 {
173 for (auto i : range(4))
174 s_[i] = s[i] ? -1 : 1;
175 }
176
178 unsigned int size () const { return 10; }
179
186 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const
187 {
188 out.resize(10);
189
190 const auto& x = in[0];
191 const auto& y = in[1];
192
193 auto pre = 1-x;
194 out[0] = {pre*s_[0]*(3*x-1), 0};
195 out[1] = { pre*3*(2*y-1), 0};
196
197 pre = x;
198 out[2] = {pre*s_[1]*(3*x-2), 0};
199 out[3] = {pre* 3*(2*y-1), 0};
200
201 pre = 1-y;
202 out[4] = {0, pre*s_[2]*(3*y-1)};
203 out[5] = {0, pre*3*(2*x-1)};
204
205 pre = y;
206 out[6] = {0, pre*s_[3]*(3*y-2)};
207 out[7] = {0, pre*3*(2*x-1)};
208
209 out[8] = {6*x*(1-x), 0};
210
211 out[9] = {0, 6*y*(1-y)};
212 }
213
220 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const
221 {
222 out.resize(10);
223
224 const auto& x = in[0];
225 const auto& y = in[1];
226
227 out[0] = {{s_[0]*(4-6*x), 0}, {0, 0}};
228 out[1] = {{3-6*y, 6-6*x}, {0, 0}};
229
230 out[2] = {{s_[1]*(6*x-2), 0}, {0, 0}};
231 out[3] = {{6*y-3, 6*x}, {0, 0}};
232
233 out[4] = {{0, 0}, {0, s_[2]*(4-6*y)}};
234 out[5] = {{0, 0}, {6-6*y, 3-6*x}};
235
236 out[6] = {{0, 0}, {0, s_[3]*(6*y-2)}};
237 out[7] = {{0, 0}, {6*y, 6*x-3}};
238
239 out[8] = {{6-12*x, 0}, {0, 0}};
240
241 out[9] = {{0, 0}, {0, 6-12*y}};
242 }
243
251 void partial (const std::array<unsigned int, 2>& order,
252 const DomainType& in,
253 std::vector<RangeType>& out) const
254 {
255 if (std::accumulate(order.begin(), order.end(), 0) == 0)
256 evaluateFunction(in, out);
257 else
258 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
259 }
260
262 unsigned int order () const { return 2; }
263
264 private:
265 std::array<R, 4> s_;
266 };
267
268
274 template<class D, class R >
275 class BDFMCubeLocalBasis<D, R, 2, 3>
276 {
280
281 public:
283
286 {
287 std::fill(s_.begin(), s_.end(), 1);
288 }
289
295 BDFMCubeLocalBasis (std::bitset<4> s)
296 {
297 for (auto i : range(4))
298 s_[i] = s[i] ? -1 : 1;
299 }
300
302 unsigned int size () const { return 18; }
303
310 inline void evaluateFunction (const DomainType& in, std::vector<RangeType>& out) const
311 {
312 out.resize(18);
313
314 const auto& x = in[0];
315 const auto& y = in[1];
316
317 auto pre = 1-x;
318 out[0] = {pre*s_[0]*-1*(10*x*x-8*x+1), 0};
319 out[1] = {pre* 3*(1-3*x)*(2*y-1), 0};
320 out[2] = {pre* s_[0]*-5*(6*y*y-6*y+1), 0};
321
322 pre = x;
323 out[3] = {pre*s_[1]*(10*x*x-12*x+3), 0};
324 out[4] = {pre* 3*(3*x-2)*(2*y-1), 0};
325 out[5] = {pre*s_[1]*5*(6*y*y-6*y+1), 0};
326
327 pre = 1-y;
328 out[6] = {0, pre*s_[2]*-1*(10*y*y-8*y+1)};
329 out[7] = {0, pre* 3*(1-3*y)*(2*x-1)};
330 out[8] = {0, pre*s_[2]*-5*( 6*x*x-6*x+1)};
331
332 pre = y;
333 out[9] = {0, pre*s_[3]*(10*y*y-12*y+3)};
334 out[10] = {0, pre* 3*(3*y-2)*(2*x-1)};
335 out[11] = {0, pre*s_[3]*5*(6*x*x-6*x+1)};
336
337 pre = x*(1-x);
338 out[12] = {pre* 6 , 0};
339 out[14] = {pre*30*(2*x-1), 0};
340 out[16] = {pre*18*(2*y-1), 0};
341
342 pre = y*(1-y);
343 out[13] = {0, pre* 6 };
344 out[15] = {0, pre*18*(2*x-1)};
345 out[17] = {0, pre*30*(2*y-1)};
346 }
347
354 inline void evaluateJacobian (const DomainType& in, std::vector<JacobianType>& out) const
355 {
356 out.resize(18);
357
358 const auto& x = in[0];
359 const auto& y = in[1];
360
361 out[0] = {{s_[0]*(30*x*x-36*x+9), 0}, {0, 0}};
362 out[1] = {{ 6*(6*x*y-3*x-4*y+2), 6*(3*x*x-4*x+1)}, {0, 0}};
363 out[2] = {{s_[0]*5*(6*y*y-6*y+1), s_[0]*30*(x-1)*(2*y-1)}, {0, 0}};
364
365 out[3] = {{ s_[1]*30*x*x-24*x+3, 0}, {0, 0}};
366 out[4] = {{ 6*(3*x-1)*(2*y-1), 6*x*(3*x-2)}, {0, 0}};
367 out[5] = {{s_[1]*5*(6*y*y-6*y+1), s_[1]*30*x*(2*y-1)}, {0, 0}};
368
369 out[6] = {{0, 0}, { 0, s_[2]*(30*y*y-36*y+9)}};
370 out[7] = {{0, 0}, { 6*(3*y*y-4*y+1), 6*(6*y*x-3*y-4*x+2)}};
371 out[8] = {{0, 0}, {s_[2]*30*(y-1)*(2*x-1), s_[2]*5*(6*x*x-6*x+1)}};
372
373 out[9] = {{0, 0}, { 0, s_[3]*30*y*y-24*y+3}};
374 out[10] = {{0, 0}, { 6*y*(3*y-2), 6*(3*y-1)*(2*x-1)}};
375 out[11] = {{0, 0}, {s_[3]*30*y*(2*x-1), s_[3]*5*(6*x*x-6*x+1)}};
376
377 out[12] = {{ -6*(2*x-1), 0}, {0, 0}};
378 out[14] = {{ -30*(6*x*x-6*x+1), 0}, {0, 0}};
379 out[16] = {{-18*(2*x-1)*(2*y-1), 36*x*(1-x)}, {0, 0}};
380
381 out[13] = {{0, 0}, { 0, -6*(2*y-1)}};
382 out[15] = {{0, 0}, {36*y*(1-y), -18*(2*y-1)*(2*x-1)}};
383 out[17] = {{0, 0}, { 0, -30*(6*y*y-6*y+1)}};
384 }
385
393 void partial (const std::array<unsigned int, 2>& order,
394 const DomainType& in,
395 std::vector<RangeType>& out) const
396 {
397 if (std::accumulate(order.begin(), order.end(), 0) == 0)
398 evaluateFunction(in, out);
399 else
400 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
401 }
402
404 unsigned int order () const { return 3; }
405
406 private:
407 std::array<R, 4> s_;
408 };
409
410} // namespace Dune
411
412#endif // #ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASFORTINMARINI_CUBE_LOCALBASIS_HH
unsigned int size() const
number of shape functions
Definition: localbasis.hh:86
unsigned int order() const
Polynomial order of the shape functions.
Definition: localbasis.hh:138
void evaluateFunction(const DomainType &in, std::vector< RangeType > &out) const
Evaluate all shape functions.
Definition: localbasis.hh:94
BDFMCubeLocalBasis()
Standard constructor.
Definition: localbasis.hh:69
void partial(const std::array< unsigned int, 2 > &order, const DomainType &in, std::vector< RangeType > &out) const
Evaluate all partial derivatives of all shape functions.
Definition: localbasis.hh:127
void evaluateJacobian(const DomainType &in, std::vector< JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: localbasis.hh:110
BDFMCubeLocalBasis(std::bitset< 4 > s)
Make set number s, where 0<= s < 16.
Definition: localbasis.hh:79
void evaluateFunction(const DomainType &in, std::vector< RangeType > &out) const
Evaluate all shape functions.
Definition: localbasis.hh:186
void partial(const std::array< unsigned int, 2 > &order, const DomainType &in, std::vector< RangeType > &out) const
Evaluate all partial derivatives of all shape functions.
Definition: localbasis.hh:251
unsigned int size() const
number of shape functions
Definition: localbasis.hh:178
BDFMCubeLocalBasis()
Standard constructor.
Definition: localbasis.hh:161
void evaluateJacobian(const DomainType &in, std::vector< JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: localbasis.hh:220
unsigned int order() const
Polynomial order of the shape functions.
Definition: localbasis.hh:262
BDFMCubeLocalBasis(std::bitset< 4 > s)
Make set number s, where 0<= s < 16.
Definition: localbasis.hh:171
BDFMCubeLocalBasis()
Standard constructor.
Definition: localbasis.hh:285
void evaluateJacobian(const DomainType &in, std::vector< JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: localbasis.hh:354
BDFMCubeLocalBasis(std::bitset< 4 > s)
Make set number s, where 0<= s < 16.
Definition: localbasis.hh:295
unsigned int size() const
number of shape functions
Definition: localbasis.hh:302
unsigned int order() const
Polynomial order of the shape functions.
Definition: localbasis.hh:404
void evaluateFunction(const DomainType &in, std::vector< RangeType > &out) const
Evaluate all shape functions.
Definition: localbasis.hh:310
void partial(const std::array< unsigned int, 2 > &order, const DomainType &in, std::vector< RangeType > &out) const
Evaluate all partial derivatives of all shape functions.
Definition: localbasis.hh:393
Brezzi-Douglas-Fortin-Marini shape functions on a reference cube.
Definition: localbasis.hh:37
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception for dummy implementations.
Definition: exceptions.hh:263
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.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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
Utilities for reduction like operations on ranges.
template which always yields a false value
Definition: typetraits.hh:124
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)