Dune Core Modules (unstable)

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