DUNE PDELab (2.8)

brezzidouglasmarini1cube3dlocalbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
4#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
5
6#include <array>
7#include <bitset>
8#include <numeric>
9#include <vector>
10
11#include <dune/common/fmatrix.hh>
12
13#include "../../common/localbasis.hh"
14
15namespace Dune
16{
27 template<class D, class R>
28 class BDM1Cube3DLocalBasis
29 {
30
31 public:
32 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,
33 R,3,Dune::FieldVector<R,3>,
34 Dune::FieldMatrix<R,3,3> > Traits;
35
37 BDM1Cube3DLocalBasis()
38 {
39 for (size_t i=0; i<6; i++)
40 sign_[i] = 1.0;
41 }
42
48 BDM1Cube3DLocalBasis(std::bitset<6> s)
49 {
50 for (size_t i=0; i<6; i++)
51 sign_[i] = s[i] ? -1.0 : 1.0;
52 }
53
55 unsigned int size() const
56 {
57 return 18;
58 }
59
66 inline void evaluateFunction(const typename Traits::DomainType& in,
67 std::vector<typename Traits::RangeType>& out) const
68 {
69 out.resize(size());
70
71 out[0][0] = sign_[0] * (in[0] - 1.0);
72 out[0][1] = 0;
73 out[0][2] = 0;
74 out[1][0] = sign_[1] * in[0];
75 out[1][1] = 0;
76 out[1][2] = 0;
77 out[2][0] = 0;
78 out[2][1] = sign_[2] * (in[1] - 1.0);
79 out[2][2] = 0;
80 out[3][0] = 0;
81 out[3][1] = sign_[3] * in[1];
82 out[3][2] = 0;
83 out[4][0] = 0;
84 out[4][1] = 0;
85 out[4][2] = sign_[4] * (in[2] - 1.0);
86 out[5][0] = 0;
87 out[5][1] = 0;
88 out[5][2] = sign_[5] * in[2];
89 out[6][0] = 6.0 * in[0] * in[1] - 3 * in[0]-6 * in[1] + 3.0;
90 out[6][1] = -3.0 * in[1] * in[1] + 3 * in[1];
91 out[6][2] = 0;
92 out[7][0] = -6.0 * in[0] * in[1] + 3 * in[0];
93 out[7][1] = 3.0 * in[1] * in[1] - 3 * in[1];
94 out[7][2] = 0;
95 out[8][0] = 3.0 * in[0] * in[0] - 3 * in[0];
96 out[8][1] = -6.0 * in[0] * in[1] + 3 * in[1]+6 * in[0]-3.0;
97 out[8][2] = 0;
98 out[9][0] = -3.0 * in[0] * in[0] + 3 * in[0];
99 out[9][1] = 6.0 * in[0] * in[1] - 3 * in[1];
100 out[9][2] = 0;
101 out[10][0] = -3.0 * in[0] * in[0] + 3 * in[0];
102 out[10][1] = 0;
103 out[10][2] = 6.0 * in[0] * in[2]-6 * in[0]-3 * in[2] + 3.0;
104 out[11][0] = 3.0 * in[0] * in[0]-3 * in[0];
105 out[11][1] = 0;
106 out[11][2] = -6.0 * in[0] * in[2] + 3 * in[2];
107 out[12][0] = -6.0 * in[0] * in[2]+6 * in[2] + 3 * in[0]-3.0;
108 out[12][1] = 0;
109 out[12][2] = 3.0 * in[2] * in[2]-3 * in[2];
110 out[13][0] = -3 * in[0]+6 * in[0] * in[2];
111 out[13][1] = 0;
112 out[13][2] = -3.0 * in[2] * in[2] + 3 * in[2];
113 out[14][0] = 0;
114 out[14][1] = 6.0 * in[1] * in[2]-3 * in[1]-6 * in[2] + 3.0;
115 out[14][2] = -3 * in[2] * in[2] + 3 * in[2];
116 out[15][0] = 0;
117 out[15][1] = -6.0 * in[1] * in[2] + 3 * in[1];
118 out[15][2] = 3.0 * in[2] * in[2]-3 * in[2];
119 out[16][0] = 0;
120 out[16][1] = 3.0 * in[1] * in[1]-3 * in[1];
121 out[16][2] = -6.0 * in[1] * in[2] + 3 * in[2]+6 * in[1]-3.0;
122 out[17][0] = 0;
123 out[17][1] = -3.0 * in[1] * in[1] + 3 * in[1];
124 out[17][2] = 6.0 * in[1] * in[2] - 3.0 * in[2];
125 }
126
133 inline void evaluateJacobian(const typename Traits::DomainType& in,
134 std::vector<typename Traits::JacobianType>& out) const
135 {
136 out.resize(size());
137
138 out[0][0] = { sign_[0], 0, 0};
139 out[0][1] = { 0, 0, 0};
140 out[0][2] = { 0, 0, 0};
141
142 out[1][0] = { sign_[1], 0, 0};
143 out[1][1] = { 0, 0, 0};
144 out[1][2] = { 0, 0, 0};
145
146 out[2][0] = { 0, 0, 0};
147 out[2][1] = { 0, sign_[2], 0};
148 out[2][2] = { 0, 0, 0};
149
150 out[3][0] = { 0, 0, 0};
151 out[3][1] = { 0, sign_[3], 0};
152 out[3][2] = { 0, 0, 0};
153
154 out[4][0] = { 0, 0, 0};
155 out[4][1] = { 0, 0, 0};
156 out[4][2] = { 0, 0, sign_[4]};
157
158 out[5][0] = { 0, 0, 0};
159 out[5][1] = { 0, 0, 0};
160 out[5][2] = { 0, 0, sign_[5]};
161
162 out[6][0] = { 6*in[1]-3, 6*in[0]-6, 0};
163 out[6][1] = { 0, -6*in[1]+3, 0};
164 out[6][2] = { 0, 0, 0};
165
166 out[7][0] = {-6*in[1]+3, -6*in[0], 0};
167 out[7][1] = { 0, 6*in[1]-3, 0};
168 out[7][2] = { 0, 0, 0};
169
170 out[8][0] = { 6*in[0]-3, 0, 0};
171 out[8][1] = {-6*in[1]+6, -6*in[0]+3, 0};
172 out[8][2] = { 0, 0, 0};
173
174 out[9][0] = {-6*in[0]+3, 0, 0};
175 out[9][1] = { 6*in[1], 6*in[0]-3, 0};
176 out[9][2] = { 0, 0, 0};
177
178 out[10][0] = {-6*in[0]+3, 0, 0};
179 out[10][1] = { 0, 0, 0};
180 out[10][2] = { 6*in[2]-6, 0, 6*in[0]-3};
181
182 out[11][0] = { 6*in[0]-3, 0, 0};
183 out[11][1] = { 0, 0, 0};
184 out[11][2] = { -6*in[2], 0, -6*in[0]+3};
185
186 out[12][0] = {-6*in[2]+3, 0, -6*in[0]+6};
187 out[12][1] = { 0, 0, 0};
188 out[12][2] = { 0, 0, 6*in[2]-3};
189
190 out[13][0] = { 6*in[2]-3, 0, 6*in[0]};
191 out[13][1] = { 0, 0, 0};
192 out[13][2] = { 0, 0, -6*in[2]+3};
193
194 out[14][0] = { 0, 0, 0};
195 out[14][1] = { 0, 6*in[2]-3, 6*in[1]-6};
196 out[14][2] = { 0, 0, -6*in[2]+3};
197
198 out[15][0] = { 0, 0, 0};
199 out[15][1] = { 0, -6*in[2]+3, -6*in[1]};
200 out[15][2] = { 0, 0, 6*in[2]-3};
201
202 out[16][0] = { 0, 0, 0};
203 out[16][1] = { 0, 6*in[1]-3, 0};
204 out[16][2] = { 0, -6*in[2]+6, -6*in[1]+3};
205
206 out[17][0] = { 0, 0, 0};
207 out[17][1] = { 0, -6*in[1]+3, 0};
208 out[17][2] = { 0, 6*in[2], 6*in[1]-3};
209 }
210
212 void partial (const std::array<unsigned int, 3>& order,
213 const typename Traits::DomainType& in, // position
214 std::vector<typename Traits::RangeType>& out) const // return value
215 {
216 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
217 if (totalOrder == 0) {
218 evaluateFunction(in, out);
219 } else if (totalOrder == 1) {
220 out.resize(size());
221 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
222
223 switch (direction) {
224 case 0:
225 out[0] = { sign_[0], 0, 0};
226 out[1] = { sign_[1], 0, 0};
227 out[2] = { 0, 0, 0};
228 out[3] = { 0, 0, 0};
229 out[4] = { 0, 0, 0};
230 out[5] = { 0, 0, 0};
231 out[6] = { 6*in[1]-3, 0, 0};
232 out[7] = {-6*in[1]+3, 0, 0};
233 out[8] = { 6*in[0]-3, -6*in[1]+6, 0};
234 out[9] = {-6*in[0]+3, 6*in[1], 0};
235 out[10] = {-6*in[0]+3, 0, 6*in[2]-6};
236 out[11] = { 6*in[0]-3, 0, -6*in[2]};
237 out[12] = {-6*in[2]+3, 0, 0};
238 out[13] = { 6*in[2]-3, 0, 0};
239 out[14] = { 0, 0, 0};
240 out[15] = { 0, 0, 0};
241 out[16] = { 0, 0, 0};
242 out[17] = { 0, 0, 0};
243 break;
244 case 1:
245 out[0] = { 0, 0, 0};
246 out[1] = { 0, 0, 0};
247 out[2] = { 0, sign_[2], 0};
248 out[3] = { 0, sign_[3], 0};
249 out[4] = { 0, 0, 0};
250 out[5] = { 0, 0, 0};
251 out[6] = { 6*in[0]-6, -6*in[1]+3, 0};
252 out[7] = { -6*in[0], 6*in[1]-3, 0};
253 out[8] = { 0, -6*in[0]+3, 0};
254 out[9] = { 0, 6*in[0]-3, 0};
255 out[10] = { 0, 0, 0};
256 out[11] = { 0, 0, 0};
257 out[12] = { 0, 0, 0};
258 out[13] = { 0, 0, 0};
259 out[14] = { 0, 6*in[2]-3, 0};
260 out[15] = { 0, -6*in[2]+3, 0};
261 out[16] = { 0, 6*in[1]-3, -6*in[2]+6};
262 out[17] = { 0, -6*in[1]+3, 6*in[2]};
263 break;
264 case 2:
265 out[0] = { 0, 0, 0};
266 out[1] = { 0, 0, 0};
267 out[2] = { 0, 0, 0};
268 out[3] = { 0, 0, 0};
269 out[4] = { 0, 0, sign_[4]};
270 out[5] = { 0, 0, sign_[5]};
271 out[6] = { 0, 0, 0};
272 out[7] = { 0, 0, 0};
273 out[8] = { 0, 0, 0};
274 out[9] = { 0, 0, 0};
275 out[10] = { 0, 0, 6*in[0]-3};
276 out[11] = { 0, 0, -6*in[0]+3};
277 out[12] = {-6*in[0]+6, 0, 6*in[2]-3};
278 out[13] = { 6*in[0], 0, -6*in[2]+3};
279 out[14] = { 0, 6*in[1]-6, -6*in[2]+3};
280 out[15] = { 0, -6*in[1], 6*in[2]-3};
281 out[16] = { 0, 0, -6*in[1]+3};
282 out[17] = { 0, 0, 6*in[1]-3};
283 break;
284 default:
285 DUNE_THROW(RangeError, "Component out of range.");
286 }
287 } else {
288 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
289 }
290 }
291
293 unsigned int order() const
294 {
295 return 2;
296 }
297
298 private:
299 std::array<R,6> sign_;
300 };
301} // end namespace Dune
302#endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_CUBE3D_LOCALBASIS_HH
Dune::BDM1Cube3DLocalBasis
First order Brezzi-Douglas-Marini shape functions on the reference hexahedron.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:29
Dune::BDM1Cube3DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:133
Dune::BDM1Cube3DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:66
Dune::BDM1Cube3DLocalBasis::partial
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:212
Dune::BDM1Cube3DLocalBasis::BDM1Cube3DLocalBasis
BDM1Cube3DLocalBasis(std::bitset< 6 > s)
Make set number s, where 0 <= s < 64.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:48
Dune::BDM1Cube3DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:293
Dune::BDM1Cube3DLocalBasis::BDM1Cube3DLocalBasis
BDM1Cube3DLocalBasis()
Standard constructor.
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:37
Dune::BDM1Cube3DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: brezzidouglasmarini1cube3dlocalbasis.hh:55
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:69
Dune::NotImplemented
Default exception for dummy implementations.
Definition: exceptions.hh:261
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:252
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Dune
Dune namespace.
Definition: alignedallocator.hh:11
Dune::LocalBasisTraits
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 21, 23:30, 2024)