3 #ifndef DUNE_PK2DLOCALBASIS_HH 4 #define DUNE_PK2DLOCALBASIS_HH 6 #include <dune/common/fmatrix.hh> 24 template<
class D,
class R,
unsigned int k>
30 enum {
N = (k+1)*(k+2)/2};
38 Dune::FieldMatrix<R,1,2>, 2 >
Traits;
43 for (
unsigned int i=0; i<=k; i++)
44 pos[i] = (1.0*i)/std::max(k,(
unsigned int)1);
55 std::vector<typename Traits::RangeType>& out)
const 65 for (
unsigned int j=0; j<=k; j++)
66 for (
unsigned int i=0; i<=k-j; i++)
69 for (
unsigned int alpha=0; alpha<i; alpha++)
70 out[n] *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
71 for (
unsigned int beta=0; beta<j; beta++)
72 out[n] *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
73 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
74 out[n] *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
82 std::vector<typename Traits::JacobianType>& out)
const 88 out[0][0][0] = 0; out[0][0][1] = 0;
93 for (
unsigned int j=0; j<=k; j++)
94 for (
unsigned int i=0; i<=k-j; i++)
99 for (
unsigned int beta=0; beta<j; beta++)
100 factor *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
101 for (
unsigned int a=0; a<i; a++)
104 for (
unsigned int alpha=0; alpha<i; alpha++)
106 product *= D(1)/(pos[i]-pos[alpha]);
108 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
109 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
110 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
111 out[n][0][0] += product;
113 for (
unsigned int c=i+j+1; c<=k; c++)
116 for (
unsigned int alpha=0; alpha<i; alpha++)
117 product *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
118 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
120 product *= -D(1)/(pos[gamma]-pos[i]-pos[j]);
122 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
123 out[n][0][0] += product;
129 for (
unsigned int alpha=0; alpha<i; alpha++)
130 factor *= (x[0]-pos[alpha])/(pos[i]-pos[alpha]);
131 for (
unsigned int b=0; b<j; b++)
134 for (
unsigned int beta=0; beta<j; beta++)
136 product *= D(1)/(pos[j]-pos[beta]);
138 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
139 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
140 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
141 out[n][0][1] += product;
143 for (
unsigned int c=i+j+1; c<=k; c++)
146 for (
unsigned int beta=0; beta<j; beta++)
147 product *= (x[1]-pos[beta])/(pos[j]-pos[beta]);
148 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
150 product *= -D(1)/(pos[gamma]-pos[i]-pos[j]);
152 product *= (pos[gamma]-x[0]-x[1])/(pos[gamma]-pos[i]-pos[j]);
153 out[n][0][1] += product;
162 template<
unsigned int dOrder>
163 inline void evaluate(
const std::array<int,dOrder>& directions,
165 std::vector<typename Traits::RangeType>& out)
const 170 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
177 for (
unsigned int j=0; j<=k; j++)
178 for (
unsigned int i=0; i<=k-j; i++, n++)
181 for (
unsigned int no1=0; no1 < k; no1++)
183 R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
184 for (
unsigned int no2=0; no2 < k; no2++)
186 factor *= lagrangianFactor(no2, i, j, in);
196 std::fill(out.begin(), out.end(), 0.0);
202 for (
unsigned int j=0; j<=k; j++)
203 for (
unsigned int i=0; i<=k-j; i++, n++)
207 for (
unsigned int no1=0; no1 < k; no1++)
209 R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
210 for (
unsigned int no2=0; no2 < k; no2++)
214 R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
215 for (
unsigned int no3=0; no3 < k; no3++)
217 if (no3 == no1 || no3 == no2)
219 factor2 *= lagrangianFactor(no3, i, j, in);
241 return (x[0]-pos[no])/(pos[i]-pos[no]);
243 return (x[1]-pos[no-i])/(pos[j]-pos[no-i]);
244 return (pos[no+1]-x[0]-x[1])/(pos[no+1]-pos[i]-pos[j]);
253 return (direction == 0) ? 1.0/(pos[i]-pos[no]) : 0;
256 return (direction == 0) ? 0: 1.0/(pos[j]-pos[no-i]);
258 return -1.0/(pos[no+1]-pos[i]-pos[j]);
Definition: pk2dlocalbasis.hh:30
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pk2dlocalbasis.hh:81
number of partial derivatives supported
Definition: localbasis.hh:74
unsigned int order() const
Polynomial order of the shape functions.
Definition: pk2dlocalbasis.hh:231
Lagrange shape functions of arbitrary order on the reference triangle.
Definition: pk2dlocalbasis.hh:25
R RangeType
range type
Definition: localbasis.hh:61
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
Pk2DLocalBasis()
Standard constructor.
Definition: pk2dlocalbasis.hh:41
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:163
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pk2dlocalbasis.hh:54
Definition: pk2dlocalbasis.hh:35
unsigned int size() const
number of shape functions
Definition: pk2dlocalbasis.hh:48
D DomainType
domain type
Definition: localbasis.hh:49
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition: pk2dlocalbasis.hh:38