DUNE PDELab (2.8)

brezzidouglasmarini2cube2dlocalbasis.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_CUBE2D_LOCALBASIS_HH
4#define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_CUBE2D_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 BDM2Cube2DLocalBasis
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 BDM2Cube2DLocalBasis()
37 {
38 for (size_t i=0; i<4; i++)
39 sign_[i] = 1.0;
40 }
41
47 BDM2Cube2DLocalBasis(std::bitset<4> s)
48 {
49 for (size_t i=0; i<4; i++)
50 sign_[i] = s[i] ? -1.0 : 1.0;
51 }
52
54 unsigned int size() const
55 {
56 return 14;
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.25 + 5.25*in[0] + 7.5*in[1] - 7.5*in[0]*in[1] - 3.0*in[0]*in[0] - 7.5*in[1]*in[1] + 7.5*in[0]*in[1]*in[1]);
71 out[0][1] = sign_[0]*(-1.25*in[1] + 3.75*in[1]*in[1] - 2.5*in[1]*in[1]*in[1]);
72 out[1][0] = 3.0 - 3.0*in[0]-6.0*in[1] + 6.0*in[0]*in[1];
73 out[1][1] = 0.0;
74 out[2][0] = sign_[0]*(-3.75 + 3.75*in[0] + 22.5*in[1] - 22.5*in[0]*in[1] - 22.5*in[1]*in[1] + 22.5*in[0]*in[1]*in[1]);
75 out[2][1] = sign_[0]*(-3.75*in[1] + 11.25*in[1]*in[1] - 7.5*in[1]*in[1]*in[1]);
76 out[3][0] = sign_[1]*(-0.75*in[0] - 7.5*in[0]*in[1] + 3.0*in[0]*in[0] + 7.5*in[0]*in[1]*in[1]);
77 out[3][1] = sign_[1]*(-1.25*in[1] + 3.75*in[1]*in[1] - 2.5*in[1]*in[1]*in[1]);
78 out[4][0] = 3.0*in[0] - 6.0*in[0]*in[1];
79 out[4][1] = 0.0;
80 out[5][0] = sign_[1]*(+3.75*in[0] - 22.5*in[0]*in[1] + 22.5*in[0]*in[1]*in[1]);
81 out[5][1] = sign_[1]*(-3.75*in[1] + 11.25*in[1]*in[1] - 7.5*in[1]*in[1]*in[1]);
82 out[6][0] = sign_[2]*(-1.25*in[0] + 3.75*in[0]*in[0] - 2.5*in[0]*in[0]*in[0]);
83 out[6][1] = sign_[2]*(-2.25 + 7.5*in[0] + 5.25*in[1] - 7.5*in[0]*in[1] - 7.5*in[0]*in[0] - 3.0*in[1]*in[1] + 7.5*in[0]*in[0]*in[1]);
84 out[7][0] = 0.0;
85 out[7][1] = -3.0 + 6.0*in[0] + 3.0*in[1] - 6.0*in[0]*in[1];
86 out[8][0] = sign_[2]*(-3.75*in[0] + 11.25*in[0]*in[0] - 7.5*in[0]*in[0]*in[0]);
87 out[8][1] = sign_[2]*(-3.75 + 22.5*in[0] + 3.75*in[1] - 22.5*in[0]*in[1] - 22.5*in[0]*in[0] + 22.5*in[0]*in[0]*in[1]);
88 out[9][0] = sign_[3]*(-1.25*in[0] + 3.75*in[0]*in[0] - 2.5*in[0]*in[0]*in[0]);
89 out[9][1] = sign_[3]*(-0.75*in[1] - 7.5*in[0]*in[1] + 3.0*in[1]*in[1] + 7.5*in[0]*in[0]*in[1]);
90 out[10][0] = 0.0;
91 out[10][1] = -3.0*in[1] + 6.0*in[0]*in[1];
92 out[11][0] = sign_[3]*(-3.75*in[0] + 11.25*in[0]*in[0] - 7.5*in[0]*in[0]*in[0]);
93 out[11][1] = sign_[3]*(3.75*in[1] - 22.5*in[0]*in[1] + 22.5*in[0]*in[0]*in[1]);
94 out[12][0] = 6.0*in[0] - 6.0*in[0]*in[0];
95 out[12][1] = 0.0;
96 out[13][0] = 0.0;
97 out[13][1] = 6.0*in[1] - 6.0*in[1]*in[1];
98 }
99
106 inline void evaluateJacobian(const typename Traits::DomainType& in,
107 std::vector<typename Traits::JacobianType>& out) const
108 {
109 out.resize(size());
110
111 out[0][0][0] = sign_[0]*(5.25 - 7.5*in[1] - 6.0*in[0] + 7.5*in[1]*in[1]);
112 out[0][0][1] = sign_[0]*(7.5 - 7.5*in[0] - 15.0*in[1] + 15.0*in[0]*in[1]);
113 out[0][1][0] = 0.0;
114 out[0][1][1] = sign_[0]*(-1.25 + 7.5*in[1] - 7.5*in[1]*in[1]);
115
116 out[1][0][0] = -3.0 + 6.0*in[1];
117 out[1][0][1] = -6.0 + 6.0*in[0];
118 out[1][1][0] = 0.0;
119 out[1][1][1] = 0.0;
120
121 out[2][0][0] = sign_[0]*(3.75 - 22.5*in[1] + 22.5*in[1]*in[1]);
122 out[2][0][1] = sign_[0]*(22.5 - 22.5*in[0] - 45.0*in[1] + 45.0*in[0]*in[1]);
123 out[2][1][0] = 0.0;
124 out[2][1][1] = sign_[0]*(-3.75 + 22.5*in[1] - 22.5*in[1]*in[1]);
125
126 out[3][0][0] = sign_[1]*(-0.75 - 7.5*in[1] + 6.0*in[0] + 7.5*in[1]*in[1]);
127 out[3][0][1] = sign_[1]*(-7.5*in[0] + 15.0*in[0]*in[1]);
128 out[3][1][0] = 0.0;
129 out[3][1][1] = sign_[1]*(-1.25 + 7.5*in[1] - 7.5*in[1]*in[1]);
130
131 out[4][0][0] = 3.0 - 6.0*in[1];
132 out[4][0][1] = -6.0*in[0];
133 out[4][1][0] = 0.0;
134 out[4][1][1] = 0.0;
135
136 out[5][0][0] = sign_[1]*(3.75 - 22.5*in[1] + 22.5*in[1]*in[1]);
137 out[5][0][1] = sign_[1]*(-22.5*in[0] + 45.0*in[0]*in[1]);
138 out[5][1][0] = 0.0;
139 out[5][1][1] = sign_[1]*(-3.75 + 22.5*in[1] - 22.5*in[1]*in[1]);
140
141 out[6][0][0] = sign_[2]*(-1.25 + 7.5*in[0] - 7.5*in[0]*in[0]);
142 out[6][0][1] = 0.0;
143 out[6][1][0] = sign_[2]*(7.5 - 7.5*in[1] - 15.0*in[0] + 15.0*in[0]*in[1]);
144 out[6][1][1] = sign_[2]*(5.25 - 7.5*in[0]- 6.0*in[1] + 7.5*in[0]*in[0]);
145
146 out[7][0][0] = 0.0;
147 out[7][0][1] = 0.0;
148 out[7][1][0] = 6.0 - 6.0*in[1];
149 out[7][1][1] = 3.0 - 6.0*in[0];
150
151 out[8][0][0] = sign_[2]*(-3.75 + 22.5*in[0] - 22.5*in[0]*in[0]);
152 out[8][0][1] = 0.0;
153 out[8][1][0] = sign_[2]*(22.5 - 22.5*in[1] - 45.0*in[0] + 45.0*in[0]*in[1]);
154 out[8][1][1] = sign_[2]*(3.75 - 22.5*in[0] + 22.5*in[0]*in[0]);
155
156 out[9][0][0] = sign_[3]*(-1.25 + 7.5*in[0] - 7.5*in[0]*in[0]);
157 out[9][0][1] = 0.0;
158 out[9][1][0] = sign_[3]*(-7.5*in[1] + 15.0*in[0]*in[1]);
159 out[9][1][1] = sign_[3]*(-0.75 - 7.5*in[0] + 6.0*in[1] + 7.5*in[0]*in[0]);
160
161 out[10][0][0] = 0.0;
162 out[10][0][1] = 0.0;
163 out[10][1][0] = 6.0*in[1];
164 out[10][1][1] = -3.0 + 6.0*in[0];
165
166 out[11][0][0] = sign_[3]*(-3.75 + 22.5*in[0] - 22.5*in[0]*in[0]);
167 out[11][0][1] = 0.0;
168 out[11][1][0] = sign_[3]*(-22.5*in[1] + 45*in[0]*in[1]);
169 out[11][1][1] = sign_[3]*(3.75 - 22.5*in[0] + 22.5*in[0]*in[0]);
170
171 out[12][0][0] = 6.0 - 12.0*in[0];
172 out[12][0][1] = 0.0;
173 out[12][1][0] = 0.0;
174 out[12][1][1] = 0.0;
175
176 out[13][0][0] = 0.0;
177 out[13][0][1] = 0.0;
178 out[13][1][0] = 0.0;
179 out[13][1][1] = 6.0 - 12.0*in[1];
180 }
181
183 void partial (const std::array<unsigned int, 2>& order,
184 const typename Traits::DomainType& in, // position
185 std::vector<typename Traits::RangeType>& out) const // return value
186 {
187 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
188 if (totalOrder == 0) {
189 evaluateFunction(in, out);
190 } else if (totalOrder == 1) {
191 out.resize(size());
192 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
193
194 switch (direction) {
195 case 0:
196 out[0][0] = sign_[0]*(5.25 - 7.5*in[1] - 6.0*in[0] + 7.5*in[1]*in[1]);
197 out[0][1] = 0.0;
198
199 out[1][0] = -3.0 + 6.0*in[1];
200 out[1][1] = 0.0;
201
202 out[2][0] = sign_[0]*(3.75 - 22.5*in[1] + 22.5*in[1]*in[1]);
203 out[2][1] = 0.0;
204
205 out[3][0] = sign_[1]*(-0.75 - 7.5*in[1] + 6.0*in[0] + 7.5*in[1]*in[1]);
206 out[3][1] = 0.0;
207
208 out[4][0] = 3.0 - 6.0*in[1];
209 out[4][1] = 0.0;
210
211 out[5][0] = sign_[1]*(3.75 - 22.5*in[1] + 22.5*in[1]*in[1]);
212 out[5][1] = 0.0;
213
214 out[6][0] = sign_[2]*(-1.25 + 7.5*in[0] - 7.5*in[0]*in[0]);
215 out[6][1] = sign_[2]*(7.5 - 7.5*in[1] - 15.0*in[0] + 15.0*in[0]*in[1]);
216
217 out[7][0] = 0.0;
218 out[7][1] = 6.0 - 6.0*in[1];
219
220 out[8][0] = sign_[2]*(-3.75 + 22.5*in[0] - 22.5*in[0]*in[0]);
221 out[8][1] = sign_[2]*(22.5 - 22.5*in[1] - 45.0*in[0] + 45.0*in[0]*in[1]);
222
223 out[9][0] = sign_[3]*(-1.25 + 7.5*in[0] - 7.5*in[0]*in[0]);
224 out[9][1] = sign_[3]*(-7.5*in[1] + 15.0*in[0]*in[1]);
225
226 out[10][0] = 0.0;
227 out[10][1] = 6.0*in[1];
228
229 out[11][0] = sign_[3]*(-3.75 + 22.5*in[0] - 22.5*in[0]*in[0]);
230 out[11][1] = sign_[3]*(-22.5*in[1] + 45*in[0]*in[1]);
231
232 out[12][0] = 6.0 - 12.0*in[0];
233 out[12][1] = 0.0;
234
235 out[13][0] = 0.0;
236 out[13][1] = 0.0;
237 break;
238 case 1:
239 out[0][0] = sign_[0]*(7.5 - 7.5*in[0] - 15.0*in[1] + 15.0*in[0]*in[1]);
240 out[0][1] = sign_[0]*(-1.25 + 7.5*in[1] - 7.5*in[1]*in[1]);
241
242 out[1][0] = -6.0 + 6.0*in[0];
243 out[1][1] = 0.0;
244
245 out[2][0] = sign_[0]*(22.5 - 22.5*in[0] - 45.0*in[1] + 45.0*in[0]*in[1]);
246 out[2][1] = sign_[0]*(-3.75 + 22.5*in[1] - 22.5*in[1]*in[1]);
247
248 out[3][0] = sign_[1]*(-7.5*in[0] + 15.0*in[0]*in[1]);
249 out[3][1] = sign_[1]*(-1.25 + 7.5*in[1] - 7.5*in[1]*in[1]);
250
251 out[4][0] = -6.0*in[0];
252 out[4][1] = 0.0;
253
254 out[5][0] = sign_[1]*(-22.5*in[0] + 45.0*in[0]*in[1]);
255 out[5][1] = sign_[1]*(-3.75 + 22.5*in[1] - 22.5*in[1]*in[1]);
256
257 out[6][0] = 0.0;
258 out[6][1] = sign_[2]*(5.25 - 7.5*in[0]- 6.0*in[1] + 7.5*in[0]*in[0]);
259
260 out[7][0] = 0.0;
261 out[7][1] = 3.0 - 6.0*in[0];
262
263 out[8][0] = 0.0;
264 out[8][1] = sign_[2]*(3.75 - 22.5*in[0] + 22.5*in[0]*in[0]);
265
266 out[9][0] = 0.0;
267 out[9][1] = sign_[3]*(-0.75 - 7.5*in[0] + 6.0*in[1] + 7.5*in[0]*in[0]);
268
269 out[10][0] = 0.0;
270 out[10][1] = -3.0 + 6.0*in[0];
271
272 out[11][0] = 0.0;
273 out[11][1] = sign_[3]*(3.75 - 22.5*in[0] + 22.5*in[0]*in[0]);
274
275 out[12][0] = 0.0;
276 out[12][1] = 0.0;
277
278 out[13][0] = 0.0;
279 out[13][1] = 6.0 - 12.0*in[1];
280 break;
281 default:
282 DUNE_THROW(RangeError, "Component out of range.");
283 }
284 } else {
285 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
286 }
287 }
288
290 unsigned int order() const
291 {
292 return 3;
293 }
294
295 private:
296 std::array<R,4> sign_;
297 };
298} // end namespace Dune
299#endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_CUBE2D_LOCALBASIS_HH
Dune::BDM2Cube2DLocalBasis
First order Brezzi-Douglas-Marini shape functions on quadrilaterals.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:28
Dune::BDM2Cube2DLocalBasis::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: brezzidouglasmarini2cube2dlocalbasis.hh:183
Dune::BDM2Cube2DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:290
Dune::BDM2Cube2DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:65
Dune::BDM2Cube2DLocalBasis::BDM2Cube2DLocalBasis
BDM2Cube2DLocalBasis()
Standard constructor.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:36
Dune::BDM2Cube2DLocalBasis::BDM2Cube2DLocalBasis
BDM2Cube2DLocalBasis(std::bitset< 4 > s)
Make set number s, where 0 <= s < 16.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:47
Dune::BDM2Cube2DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:54
Dune::BDM2Cube2DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: brezzidouglasmarini2cube2dlocalbasis.hh:106
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)