3 #ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
4 #define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
11 #include <dune/common/fvector.hh>
12 #include <dune/common/fmatrix.hh>
19 template<
class D,
class R,
int dim>
25 DUNE_THROW(Dune::NotImplemented,
"HierarchicalSimplexP2LocalBasis not implemented for dim > 3.");
43 template<
class D,
class R>
52 unsigned int size ()
const
59 std::vector<typename Traits::RangeType>& out)
const
65 out[2] = 1-4*(in[0]-0.5)*(in[0]-0.5);
71 std::vector<typename Traits::JacobianType>& out)
const
77 out[2][0][0] = 4-8*in[0];
82 unsigned int order ()
const
109 template<
class D,
class R>
118 unsigned int size ()
const
125 std::vector<typename Traits::RangeType>& out)
const
129 out[0] = 1 - in[0] - in[1];
130 out[1] = 4*in[0]*(1-in[0]-in[1]);
132 out[3] = 4*in[1]*(1-in[0]-in[1]);
133 out[4] = 4*in[0]*in[1];
135 out[6] = 27*in[0]*in[1]*(1-in[0]-in[1]);
142 std::vector<typename Traits::JacobianType>& out)
const
146 out[0][0][0] = -1; out[0][0][1] = -1;
147 out[1][0][0] = 4-8*in[0]-4*in[1]; out[1][0][1] = -4*in[0];
148 out[2][0][0] = 1; out[2][0][1] = 0;
149 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1];
150 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0];
151 out[5][0][0] = 0; out[5][0][1] = 1;
154 out[6][0][0] = 27 * in[1] * (1 - 2*in[0] - in[1]);
155 out[6][0][1] = 27 * in[0] * (1 - 2*in[1] - in[0]);
161 unsigned int order ()
const
192 template<
class D,
class R>
201 unsigned int size ()
const
208 std::vector<typename Traits::RangeType>& out)
const
212 out[0] = 1 - in[0] - in[1] - in[2];
213 out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]);
215 out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]);
216 out[4] = 4 * in[0] * in[1];
218 out[6] = 4 * in[2] * (1 - in[0] - in[1] - in[2]);
219 out[7] = 4 * in[0] * in[2];
220 out[8] = 4 * in[1] * in[2];
224 out[10] = 81*in[0]*in[1]*in[2]*(1-in[0]-in[1]-in[2]);
229 std::vector<typename Traits::JacobianType>& out)
const
233 out[0][0][0] = -1; out[0][0][1] = -1; out[0][0][2] = -1;
234 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];
235 out[2][0][0] = 1; out[2][0][1] = 0; out[2][0][2] = 0;
236 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];
237 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; out[4][0][2] = 0;
238 out[5][0][0] = 0; out[5][0][1] = 1; out[5][0][2] = 0;
239 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];
240 out[7][0][0] = 4*in[2]; out[7][0][1] = 0; out[7][0][2] = 4*in[0];
241 out[8][0][0] = 0; out[8][0][1] = 4*in[2]; out[8][0][2] = 4*in[1];
242 out[9][0][0] = 0; out[9][0][1] = 0; out[9][0][2] = 1;
244 out[10][0][0] = 81 * in[1] * in[2] * (1 - 2*in[0] - in[1] - in[2]);
245 out[10][0][1] = 81 * in[0] * in[2] * (1 - in[0] - 2*in[1] - in[2]);
246 out[10][0][2] = 81 * in[0] * in[1] * (1 - in[0] - in[1] - 2*in[2]);
252 unsigned int order ()
const
289 static const int numVertices = dim+1;
292 static const int numEdges = (dim+1)*dim / 2;
297 : li(numVertices+numEdges + 1)
300 DUNE_THROW(NotImplemented,
"only for 2d");
314 return numVertices+numEdges + 1;
324 std::vector<Dune::LocalKey> li;
333 template<
typename F,
typename C>
336 typename LB::Traits::DomainType x;
337 typename LB::Traits::RangeType y;
342 x[0] = 0.0; x[1] = 0.0; f.evaluate(x,y); out[0] = y;
343 x[0] = 1.0; x[1] = 0.0; f.evaluate(x,y); out[2] = y;
344 x[0] = 0.0; x[1] = 1.0; f.evaluate(x,y); out[5] = y;
347 x[0] = 0.5; x[1] = 0.0; f.evaluate(x,y);
348 out[1] = y - out[0]*(1-x[0]) - out[2]*x[0];
350 x[0] = 0.0; x[1] = 0.5; f.evaluate(x,y);
351 out[3] = y - out[0]*(1-x[1]) - out[5]*x[1];
353 x[0] = 0.5; x[1] = 0.5; f.evaluate(x,y);
354 out[4] = y - out[2]*x[0] - out[5]*x[1];
357 x[0] = 1.0/3; x[1] = 1.0/3; f.evaluate(x,y);
361 std::vector<typename LB::Traits::RangeType> sfValues;
362 shapeFunctions.evaluateFunction(x, sfValues);
365 for (
int i=0; i<6; i++)
366 out[6] -= out[i]*sfValues[i];