Dune Core Modules (2.6.0)

brezzidouglasmarini2simplex2dlocalbasis.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_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH
4#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_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{
26 template<class D, class R>
27 class BDM2Simplex2DLocalBasis
28 {
29
30 public:
31 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,
32 R,2,Dune::FieldVector<R,2>,
33 Dune::FieldMatrix<R,2,2> > Traits;
34
36 BDM2Simplex2DLocalBasis()
37 {
38 for (size_t i=0; i<3; i++)
39 sign_[i] = 1.0;
40 }
41
47 BDM2Simplex2DLocalBasis(std::bitset<3> s)
48 {
49 for (size_t i=0; i<3; i++)
50 sign_[i] = s[i] ? -1.0 : 1.0;
51 }
52
54 unsigned int size() const
55 {
56 return 12;
57 }
58
65 inline void evaluateFunction(const typename Traits::DomainType& in,
66 std::vector<typename Traits::RangeType>& out) const
67 {
68 out.resize(size());
69
70 out[0][0] = sign_[0]*(-2*in[0]*in[1] + in[0]*in[0]);
71 out[0][1] = sign_[0]*(-1 + 6*in[1] -2*in[0]*in[1] - 5*in[1]*in[1]);
72
73 out[1][0] = 1.5*in[0] + 3*in[0]*in[1] - 4.5*in[0]*in[0];
74 out[1][1] = -3 + 6*in[0] + 10.5*in[1] - 15*in[0]*in[1] - 7.5*in[1]*in[1];
75
76 out[2][0] = sign_[0]*(-7.5*in[0] + 5*in[0]*in[1] + 12.5*in[0]*in[0]);
77 out[2][1] = sign_[0]*(-5 + 30*in[0] + 7.5*in[1] - 25*in[0]*in[1] - 30*in[0]*in[0] - 2.5*in[1]*in[1]);
78
79
80
81 out[3][0] = sign_[1]*(-1 + 6*in[0] - 2*in[0]*in[1] - 5*in[0]*in[0]);
82 out[3][1] = sign_[1]*(-2*in[0]*in[1] + in[1]*in[1]);
83
84 out[4][0] = 3 - 10.5*in[0] - 6*in[1] + 15*in[0]*in[1] + 7.5*in[0]*in[0];
85 out[4][1] = -1.5*in[1] - 3*in[0]*in[1] + 4.5*in[1]*in[1];
86
87 out[5][0] = sign_[1]*(-5 + 7.5*in[0] + 30*in[1] - 25*in[0]*in[1] - 2.5*in[0]*in[0] - 30*in[1]*in[1]);
88 out[5][1] = sign_[1]*(-7.5*in[1] + 5*in[0]*in[1] + 12.5*in[1]*in[1]);
89
90
91
92 out[6][0] = sign_[2]*(-3*in[0] + 4*in[0]*in[1] + 4*in[0]*in[0]);
93 out[6][1] = sign_[2]*(-3*in[1] + 4*in[0]*in[1] + 4*in[1]*in[1]);
94
95 out[7][0] = -3*in[0] + 6*in[0]*in[0];
96 out[7][1] = 3*in[1] - 6*in[1]*in[1];
97
98 out[8][0] = sign_[2]*(-10*in[0]*in[1] + 5*in[0]*in[0]);
99 out[8][1] = sign_[2]*(-10*in[0]*in[1] + 5*in[1]*in[1]);
100
101
102
103 out[9][0] = 18*in[0] - 12*in[0]*in[1] - 18*in[0]*in[0];
104 out[9][1] = 6*in[1] - 12*in[0]*in[1] - 6*in[1]*in[1];
105
106 out[10][0] = 6*in[0] - 12*in[0]*in[1] - 6*in[0]*in[0];
107 out[10][1] = 18*in[1] - 12*in[0]*in[1] - 18*in[1]*in[1];
108
109 out[11][0] = 90*in[0] - 180*in[0]*in[1] - 90*in[0]*in[0];
110 out[11][1] = -90*in[1] + 180*in[0]*in[1] + 90*in[1]*in[1];
111 }
112
119 inline void evaluateJacobian(const typename Traits::DomainType& in,
120 std::vector<typename Traits::JacobianType>& out) const
121 {
122 out.resize(size());
123
124 out[0][0][0] = sign_[0]*(-2*in[1] + 2*in[0]);
125 out[0][0][1] = sign_[0]*(-2*in[0]);
126
127 out[0][1][0] = sign_[0]*(-2*in[1]);
128 out[0][1][1] = sign_[0]*(6 -2*in[0] - 10*in[1]);
129
130
131 out[1][0][0] = 1.5 + 3*in[1] - 9*in[0];
132 out[1][0][1] = 3*in[0];
133
134 out[1][1][0] = 6 - 15*in[1];
135 out[1][1][1] = 10.5 - 15*in[0] - 15*in[1];
136
137
138 out[2][0][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]);
139 out[2][0][1] = sign_[0]*(5*in[0]);
140
141 out[2][1][0] = sign_[0]*(30 - 25*in[1] - 60*in[0]);
142 out[2][1][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]);
143
144
145
146 out[3][0][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]);
147 out[3][0][1] = sign_[1]*(-2*in[0]);
148
149 out[3][1][0] = sign_[1]*(-2*in[1]);
150 out[3][1][1] = sign_[1]*(-2*in[0] + 2*in[1]);
151
152
153 out[4][0][0] = -10.5 + 15*in[1] + 15*in[0];
154 out[4][0][1] = -6 + 15*in[0];
155
156 out[4][1][0] = -3*in[1];
157 out[4][1][1] = -1.5 - 3*in[0] + 9*in[1];
158
159
160 out[5][0][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]);
161 out[5][0][1] = sign_[1]*(30 - 25*in[0] - 60*in[1]);
162
163 out[5][1][0] = sign_[1]*(5*in[1]);
164 out[5][1][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]);
165
166
167
168 out[6][0][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]);
169 out[6][0][1] = sign_[2]*(4*in[0]);
170
171 out[6][1][0] = sign_[2]*(4*in[1]);
172 out[6][1][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]);
173
174
175 out[7][0][0] = -3 + 12*in[0];
176 out[7][0][1] = 0;
177
178 out[7][1][0] = 0;
179 out[7][1][1] = 3 - 12*in[1];
180
181
182 out[8][0][0] = sign_[2]*(-10*in[1] + 10*in[0]);
183 out[8][0][1] = sign_[2]*(-10*in[0]);
184
185 out[8][1][0] = sign_[2]*(-10*in[1]);
186 out[8][1][1] = sign_[2]*(-10*in[0] + 10*in[1]);
187
188
189 out[9][0][0] = 18 - 12*in[1] - 36*in[0];
190 out[9][0][1] = -12*in[0];
191
192 out[9][1][0] = -12*in[1];
193 out[9][1][1] = 6 - 12*in[0] - 12*in[1];
194
195 out[10][0][0] = 6 - 12*in[1] - 12*in[0];
196 out[10][0][1] = -12*in[0];
197
198 out[10][1][0] = -12*in[1];
199 out[10][1][1] = 18 - 12*in[0] - 36*in[1];
200
201 out[11][0][0] = 90 - 180*in[1] - 180*in[0];
202 out[11][0][1] = -180*in[0];
203
204 out[11][1][0] = 180*in[1];
205 out[11][1][1] = -90 + 180*in[0] + 180*in[1];
206 }
207
209 void partial (const std::array<unsigned int, 2>& order,
210 const typename Traits::DomainType& in, // position
211 std::vector<typename Traits::RangeType>& out) const // return value
212 {
213 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
214 if (totalOrder == 0) {
215 evaluateFunction(in, out);
216 } else if (totalOrder == 1) {
217 out.resize(size());
218 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
219
220 switch (direction) {
221 case 0:
222 out[0][0] = sign_[0]*(-2*in[1] + 2*in[0]);
223 out[0][1] = sign_[0]*(-2*in[1]);
224
225 out[1][0] = 1.5 + 3*in[1] - 9*in[0];
226 out[1][1] = 6 - 15*in[1];
227
228 out[2][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]);
229 out[2][1] = sign_[0]*(30 - 25*in[1] - 60*in[0]);
230
231 out[3][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]);
232 out[3][1] = sign_[1]*(-2*in[1]);
233
234 out[4][0] = -10.5 + 15*in[1] + 15*in[0];
235 out[4][1] = -3*in[1];
236
237 out[5][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]);
238 out[5][1] = sign_[1]*(5*in[1]);
239
240 out[6][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]);
241 out[6][1] = sign_[2]*(4*in[1]);
242
243 out[7][0] = -3 + 12*in[0];
244 out[7][1] = 0;
245
246 out[8][0] = sign_[2]*(-10*in[1] + 10*in[0]);
247 out[8][1] = sign_[2]*(-10*in[1]);
248
249 out[9][0] = 18 - 12*in[1] - 36*in[0];
250 out[9][1] = -12*in[1];
251
252 out[10][0] = 6 - 12*in[1] - 12*in[0];
253 out[10][1] = -12*in[1];
254
255 out[11][0] = 90 - 180*in[1] - 180*in[0];
256 out[11][1] = 180*in[1];
257 break;
258 case 1:
259 out[0][0] = sign_[0]*(-2*in[0]);
260 out[0][1] = sign_[0]*(6 -2*in[0] - 10*in[1]);
261
262 out[1][0] = 3*in[0];
263 out[1][1] = 10.5 - 15*in[0] - 15*in[1];
264
265 out[2][0] = sign_[0]*(5*in[0]);
266 out[2][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]);
267
268 out[3][0] = sign_[1]*(-2*in[0]);
269 out[3][1] = sign_[1]*(-2*in[0] + 2*in[1]);
270
271 out[4][0] = -6 + 15*in[0];
272 out[4][1] = -1.5 - 3*in[0] + 9*in[1];
273
274 out[5][0] = sign_[1]*(30 - 25*in[0] - 60*in[1]);
275 out[5][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]);
276
277 out[6][0] = sign_[2]*(4*in[0]);
278 out[6][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]);
279
280 out[7][0] = 0;
281 out[7][1] = 3 - 12*in[1];
282
283 out[8][0] = sign_[2]*(-10*in[0]);
284 out[8][1] = sign_[2]*(-10*in[0] + 10*in[1]);
285
286 out[9][0] = -12*in[0];
287 out[9][1] = 6 - 12*in[0] - 12*in[1];
288
289 out[10][0] = -12*in[0];
290 out[10][1] = 18 - 12*in[0] - 36*in[1];
291
292 out[11][0] = -180*in[0];
293 out[11][1] = -90 + 180*in[0] + 180*in[1];
294 break;
295 default:
296 DUNE_THROW(RangeError, "Component out of range.");
297 }
298 } else {
299 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
300 }
301 }
302
304 unsigned int order() const
305 {
306 return 2; // TODO: check whether this is not order 3
307 }
308
309 private:
310 std::array<R,3> sign_;
311 };
312} // end namespace Dune
313#endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH
Dune::BDM2Simplex2DLocalBasis
First order Brezzi-Douglas-Marini shape functions on quadrilaterals.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:28
Dune::BDM2Simplex2DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:119
Dune::BDM2Simplex2DLocalBasis::BDM2Simplex2DLocalBasis
BDM2Simplex2DLocalBasis(std::bitset< 3 > s)
Make set number s, where 0 <= s < 8.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:47
Dune::BDM2Simplex2DLocalBasis::partial
void partial(const std::array< unsigned int, 2 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:209
Dune::BDM2Simplex2DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:54
Dune::BDM2Simplex2DLocalBasis::BDM2Simplex2DLocalBasis
BDM2Simplex2DLocalBasis()
Standard constructor.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:36
Dune::BDM2Simplex2DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:65
Dune::BDM2Simplex2DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: brezzidouglasmarini2simplex2dlocalbasis.hh:304
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:68
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
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
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:331
Dune
Dune namespace.
Definition: alignedallocator.hh:10
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 (Nov 24, 23:30, 2024)