2 #ifndef DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
3 #define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
10 #include <dune/common/fvector.hh>
11 #include <dune/common/fmatrix.hh>
18 template<
class D,
class R,
int dim>
24 DUNE_THROW(Dune::NotImplemented,
"HierarchicalSimplexP2LocalBasis not implemented for dim > 3.");
42 template<
class D,
class R>
51 unsigned int size ()
const
58 std::vector<typename Traits::RangeType>& out)
const
64 out[2] = 1-4*(in[0]-0.5)*(in[0]-0.5);
70 std::vector<typename Traits::JacobianType>& out)
const
76 out[2][0][0] = 4-8*in[0];
81 unsigned int order ()
const
108 template<
class D,
class R>
117 unsigned int size ()
const
124 std::vector<typename Traits::RangeType>& out)
const
128 out[0] = 1 - in[0] - in[1];
129 out[1] = 4*in[0]*(1-in[0]-in[1]);
131 out[3] = 4*in[1]*(1-in[0]-in[1]);
132 out[4] = 4*in[0]*in[1];
134 out[6] = 27*in[0]*in[1]*(1-in[0]-in[1]);
141 std::vector<typename Traits::JacobianType>& out)
const
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;
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]);
160 unsigned int order ()
const
191 template<
class D,
class R>
200 unsigned int size ()
const
207 std::vector<typename Traits::RangeType>& out)
const
211 out[0] = 1 - in[0] - in[1] - in[2];
212 out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]);
214 out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]);
215 out[4] = 4 * in[0] * 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];
223 out[10] = 81*in[0]*in[1]*in[2]*(1-in[0]-in[1]-in[2]);
228 std::vector<typename Traits::JacobianType>& out)
const
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;
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]);
251 unsigned int order ()
const
288 static const int numVertices = dim+1;
291 static const int numEdges = (dim+1)*dim / 2;
296 : li(numVertices+numEdges + 1)
299 DUNE_THROW(NotImplemented,
"only for 2d");
313 return numVertices+numEdges + 1;
323 std::vector<Dune::LocalKey> li;
332 template<
typename F,
typename C>
335 typename LB::Traits::DomainType x;
336 typename LB::Traits::RangeType y;
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;
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];
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];
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];
356 x[0] = 1.0/3; x[1] = 1.0/3; f.evaluate(x,y);
360 std::vector<typename LB::Traits::RangeType> sfValues;
361 shapeFunctions.evaluateFunction(x, sfValues);
364 for (
int i=0; i<6; i++)
365 out[6] -= out[i]*sfValues[i];