Dune Core Modules (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 
17 namespace 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:95
Dune::NotImplemented
Default exception for dummy implementations.
Definition: exceptions.hh:263
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:254
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:218
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
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  |  generated with Hugo v0.80.0 (Apr 29, 22:29, 2024)