dune-localfunctions  2.2.1
hierarchicalsimplexp2withelementbubble.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
3 #define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
4 
9 #include <vector>
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
12 
15 
16 namespace Dune
17 {
18  template<class D, class R, int dim>
20  {
21  public:
23  {
24  DUNE_THROW(Dune::NotImplemented,"HierarchicalSimplexP2LocalBasis not implemented for dim > 3.");
25  }
26  };
27 
42  template<class D, class R>
44  {
45  public:
47  typedef LocalBasisTraits<D,1,Dune::FieldVector<D,1>,R,1,Dune::FieldVector<R,1>,
48  Dune::FieldMatrix<R,1,1> > Traits;
49 
51  unsigned int size () const
52  {
53  return 3;
54  }
55 
57  inline void evaluateFunction (const typename Traits::DomainType& in,
58  std::vector<typename Traits::RangeType>& out) const
59  {
60  out.resize(3);
61 
62  out[0] = 1-in[0];
63  out[1] = in[0];
64  out[2] = 1-4*(in[0]-0.5)*(in[0]-0.5);
65  }
66 
68  inline void
69  evaluateJacobian (const typename Traits::DomainType& in, // position
70  std::vector<typename Traits::JacobianType>& out) const // return value
71  {
72  out.resize(3);
73 
74  out[0][0][0] = -1;
75  out[1][0][0] = 1;
76  out[2][0][0] = 4-8*in[0];
77  }
78 
81  unsigned int order () const
82  {
83  return 2;
84  }
85 
86  };
87 
108  template<class D, class R>
110  {
111  public:
113  typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
114  Dune::FieldMatrix<R,1,2> > Traits;
115 
117  unsigned int size () const
118  {
119  return 7;
120  }
121 
123  inline void evaluateFunction (const typename Traits::DomainType& in,
124  std::vector<typename Traits::RangeType>& out) const
125  {
126  out.resize(7);
127 
128  out[0] = 1 - in[0] - in[1];
129  out[1] = 4*in[0]*(1-in[0]-in[1]);
130  out[2] = in[0];
131  out[3] = 4*in[1]*(1-in[0]-in[1]);
132  out[4] = 4*in[0]*in[1];
133  out[5] = in[1];
134  out[6] = 27*in[0]*in[1]*(1-in[0]-in[1]);
135 
136  }
137 
139  inline void
140  evaluateJacobian (const typename Traits::DomainType& in, // position
141  std::vector<typename Traits::JacobianType>& out) const // return value
142  {
143  out.resize(7);
144 
145  out[0][0][0] = -1; out[0][0][1] = -1;
146  out[1][0][0] = 4-8*in[0]-4*in[1]; out[1][0][1] = -4*in[0];
147  out[2][0][0] = 1; out[2][0][1] = 0;
148  out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1];
149  out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0];
150  out[5][0][0] = 0; out[5][0][1] = 1;
151 
152  // Cubic bubble
153  out[6][0][0] = 27 * in[1] * (1 - 2*in[0] - in[1]);
154  out[6][0][1] = 27 * in[0] * (1 - 2*in[1] - in[0]);
155 
156  }
157 
160  unsigned int order () const
161  {
162  return 3;
163  }
164 
165  };
166 
191  template<class D, class R>
193  {
194  public:
196  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
197  Dune::FieldMatrix<R,1,3> > Traits;
198 
200  unsigned int size () const
201  {
202  return 11;
203  }
204 
206  void evaluateFunction (const typename Traits::DomainType& in,
207  std::vector<typename Traits::RangeType>& out) const
208  {
209  out.resize(10);
210 
211  out[0] = 1 - in[0] - in[1] - in[2];
212  out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]);
213  out[2] = in[0];
214  out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]);
215  out[4] = 4 * in[0] * in[1];
216  out[5] = in[1];
217  out[6] = 4 * in[2] * (1 - in[0] - in[1] - in[2]);
218  out[7] = 4 * in[0] * in[2];
219  out[8] = 4 * in[1] * in[2];
220  out[9] = in[2];
221 
222  // quartic element bubble
223  out[10] = 81*in[0]*in[1]*in[2]*(1-in[0]-in[1]-in[2]);
224  }
225 
227  void evaluateJacobian (const typename Traits::DomainType& in, // position
228  std::vector<typename Traits::JacobianType>& out) const // return value
229  {
230  out.resize(10);
231 
232  out[0][0][0] = -1; out[0][0][1] = -1; out[0][0][2] = -1;
233  out[1][0][0] = 4-8*in[0]-4*in[1]-4*in[2]; out[1][0][1] = -4*in[0]; out[1][0][2] = -4*in[0];
234  out[2][0][0] = 1; out[2][0][1] = 0; out[2][0][2] = 0;
235  out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1]-4*in[2]; out[3][0][2] = -4*in[1];
236  out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; out[4][0][2] = 0;
237  out[5][0][0] = 0; out[5][0][1] = 1; out[5][0][2] = 0;
238  out[6][0][0] = -4*in[2]; out[6][0][1] = -4*in[2]; out[6][0][2] = 4-4*in[0]-4*in[1]-8*in[2];
239  out[7][0][0] = 4*in[2]; out[7][0][1] = 0; out[7][0][2] = 4*in[0];
240  out[8][0][0] = 0; out[8][0][1] = 4*in[2]; out[8][0][2] = 4*in[1];
241  out[9][0][0] = 0; out[9][0][1] = 0; out[9][0][2] = 1;
242 
243  out[10][0][0] = 81 * in[1] * in[2] * (1 - 2*in[0] - in[1] - in[2]);
244  out[10][0][1] = 81 * in[0] * in[2] * (1 - in[0] - 2*in[1] - in[2]);
245  out[10][0][2] = 81 * in[0] * in[1] * (1 - in[0] - in[1] - 2*in[2]);
246  }
247 
248 
251  unsigned int order () const
252  {
253  return 4;
254  }
255 
256  };
257 
258 
284 template <int dim>
286 {
287  // The binomial coefficient: dim+1 over 1
288  static const int numVertices = dim+1;
289 
290  // The binomial coefficient: dim+1 over 2
291  static const int numEdges = (dim+1)*dim / 2;
292 
293 public:
296  : li(numVertices+numEdges + 1)
297  {
298  if (dim!=2)
299  DUNE_THROW(NotImplemented, "only for 2d");
300 
301  li[0] = Dune::LocalKey(0,2,0); // Vertex (0,0)
302  li[1] = Dune::LocalKey(0,1,0); // Edge (0.5, 0)
303  li[2] = Dune::LocalKey(1,2,0); // Vertex (1,0)
304  li[3] = Dune::LocalKey(1,1,0); // Edge (0, 0.5)
305  li[4] = Dune::LocalKey(2,1,0); // Edge (0.5, 0.5)
306  li[5] = Dune::LocalKey(2,2,0); // Vertex (0,1)
307  li[6] = Dune::LocalKey(0,0,0); // Element (1/3, 1/3)
308  }
309 
311  size_t size () const
312  {
313  return numVertices+numEdges + 1;
314  }
315 
317  const Dune::LocalKey& localKey (size_t i) const
318  {
319  return li[i];
320  }
321 
322 private:
323  std::vector<Dune::LocalKey> li;
324 };
325 
326 template<class LB>
328 {
329 public:
330 
332  template<typename F, typename C>
333  void interpolate (const F& f, std::vector<C>& out) const
334  {
335  typename LB::Traits::DomainType x;
336  typename LB::Traits::RangeType y;
337 
338  out.resize(7);
339 
340  // vertices
341  x[0] = 0.0; x[1] = 0.0; f.evaluate(x,y); out[0] = y;
342  x[0] = 1.0; x[1] = 0.0; f.evaluate(x,y); out[2] = y;
343  x[0] = 0.0; x[1] = 1.0; f.evaluate(x,y); out[5] = y;
344 
345  // edge bubbles
346  x[0] = 0.5; x[1] = 0.0; f.evaluate(x,y);
347  out[1] = y - out[0]*(1-x[0]) - out[2]*x[0];
348 
349  x[0] = 0.0; x[1] = 0.5; f.evaluate(x,y);
350  out[3] = y - out[0]*(1-x[1]) - out[5]*x[1];
351 
352  x[0] = 0.5; x[1] = 0.5; f.evaluate(x,y);
353  out[4] = y - out[2]*x[0] - out[5]*x[1];
354 
355  // element bubble
356  x[0] = 1.0/3; x[1] = 1.0/3; f.evaluate(x,y);
357 
360  std::vector<typename LB::Traits::RangeType> sfValues;
361  shapeFunctions.evaluateFunction(x, sfValues);
362 
363  out[6] = y;
364  for (int i=0; i<6; i++)
365  out[6] -= out[i]*sfValues[i];
366 
367  }
368 
369 };
370 
371 
372 }
373 #endif