DUNE-FEM (unstable)

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