3#ifndef DUNE_PK2DLOCALBASIS_HH
4#define DUNE_PK2DLOCALBASIS_HH
8#include <dune/common/fmatrix.hh>
26 template<
class D,
class R,
unsigned int k>
32 enum {
N = (k+1)*(k+2)/2};
40 Dune::FieldMatrix<R,1,2>, 2 >
Traits;
45 for (
unsigned int i=0; i<=k; i++)
46 pos_[i] = (1.0*i)/std::max(k,(
unsigned int)1);
57 std::vector<typename Traits::RangeType>& out)
const
67 for (
unsigned int j=0; j<=k; j++)
68 for (
unsigned int i=0; i<=k-j; i++)
71 for (
unsigned int alpha=0; alpha<i; alpha++)
72 out[n] *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
73 for (
unsigned int beta=0; beta<j; beta++)
74 out[n] *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
75 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
76 out[n] *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
84 std::vector<typename Traits::JacobianType>& out)
const
90 out[0][0][0] = 0; out[0][0][1] = 0;
95 for (
unsigned int j=0; j<=k; j++)
96 for (
unsigned int i=0; i<=k-j; i++)
101 for (
unsigned int beta=0; beta<j; beta++)
102 factor *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
103 for (
unsigned int a=0; a<i; a++)
106 for (
unsigned int alpha=0; alpha<i; alpha++)
108 product *= D(1)/(pos_[i]-pos_[alpha]);
110 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
111 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
112 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
113 out[n][0][0] += product;
115 for (
unsigned int c=i+j+1; c<=k; c++)
118 for (
unsigned int alpha=0; alpha<i; alpha++)
119 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
120 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
122 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
124 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
125 out[n][0][0] += product;
131 for (
unsigned int alpha=0; alpha<i; alpha++)
132 factor *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
133 for (
unsigned int b=0; b<j; b++)
136 for (
unsigned int beta=0; beta<j; beta++)
138 product *= D(1)/(pos_[j]-pos_[beta]);
140 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
141 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
142 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
143 out[n][0][1] += product;
145 for (
unsigned int c=i+j+1; c<=k; c++)
148 for (
unsigned int beta=0; beta<j; beta++)
149 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
150 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
152 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
154 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
155 out[n][0][1] += product;
170 std::vector<typename Traits::RangeType>& out)
const
172 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
181 std::array<int,1> directions;
182 directions[0] = std::find(
order.begin(),
order.end(), 1)-
order.begin();
183 evaluate<1>(directions, in, out);
188 std::array<int,2> directions;
189 unsigned int counter = 0;
190 auto nonconstOrder =
order;
191 for (
int i=0; i<2; i++)
193 while (nonconstOrder[i])
195 directions[counter++] = i;
200 evaluate<2>(directions, in, out);
204 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
209 template<
unsigned int dOrder>
210 inline void evaluate(
const std::array<int,dOrder>& directions,
212 std::vector<typename Traits::RangeType>& out)
const
217 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
224 for (
unsigned int j=0; j<=k; j++)
225 for (
unsigned int i=0; i<=k-j; i++, n++)
228 for (
unsigned int no1=0; no1 < k; no1++)
230 R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
231 for (
unsigned int no2=0; no2 < k; no2++)
233 factor *= lagrangianFactor(no2, i, j, in);
243 std::fill(out.begin(), out.end(), 0.0);
249 for (
unsigned int j=0; j<=k; j++)
250 for (
unsigned int i=0; i<=k-j; i++, n++)
254 for (
unsigned int no1=0; no1 < k; no1++)
256 R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
257 for (
unsigned int no2=0; no2 < k; no2++)
261 R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
262 for (
unsigned int no3=0; no3 < k; no3++)
264 if (no3 == no1 || no3 == no2)
266 factor2 *= lagrangianFactor(no3, i, j, in);
288 return (x[0]-pos_[no])/(pos_[i]-pos_[no]);
290 return (x[1]-pos_[no-i])/(pos_[j]-pos_[no-i]);
291 return (pos_[no+1]-x[0]-x[1])/(pos_[no+1]-pos_[i]-pos_[j]);
300 return (direction == 0) ? 1.0/(pos_[i]-pos_[no]) : 0;
303 return (direction == 0) ? 0: 1.0/(pos_[j]-pos_[no-i]);
305 return -1.0/(pos_[no+1]-pos_[i]-pos_[j]);
Definition brezzidouglasmarini1cube2d.hh:14
Type traits for LocalBasisVirtualInterface.
Definition localbasis.hh:38
R RangeType
range type
Definition localbasis.hh:61
@ diffOrder
number of partial derivatives supported
Definition localbasis.hh:74
D DomainType
domain type
Definition localbasis.hh:49
Lagrange shape functions of arbitrary order on the reference triangle.
Definition pk2dlocalbasis.hh:28
void evaluate(const std::array< int, dOrder > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate higher derivatives of all shape functions.
Definition pk2dlocalbasis.hh:210
@ O
Definition pk2dlocalbasis.hh:37
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition pk2dlocalbasis.hh:56
unsigned int size() const
number of shape functions
Definition pk2dlocalbasis.hh:50
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 any order of all shape functions.
Definition pk2dlocalbasis.hh:168
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition pk2dlocalbasis.hh:83
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition pk2dlocalbasis.hh:40
@ N
Definition pk2dlocalbasis.hh:32
Pk2DLocalBasis()
Standard constructor.
Definition pk2dlocalbasis.hh:43
unsigned int order() const
Polynomial order of the shape functions.
Definition pk2dlocalbasis.hh:278