Dune Core Modules (2.6.0)

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 
15 namespace 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: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.80.0 (Apr 27, 22:29, 2024)