3 #ifndef DUNE_PK2DLOCALBASIS_HH 4 #define DUNE_PK2DLOCALBASIS_HH 8 #include <dune/common/fmatrix.hh> 9 #include <dune/common/deprecated.hh> 27 template<
class D,
class R,
unsigned int k>
33 enum {
N = (k+1)*(k+2)/2};
41 Dune::FieldMatrix<R,1,2>, 2 >
Traits;
46 for (
unsigned int i=0; i<=k; i++)
47 pos_[i] = (1.0*i)/std::max(k,(
unsigned int)1);
58 std::vector<typename Traits::RangeType>& out)
const 68 for (
unsigned int j=0; j<=k; j++)
69 for (
unsigned int i=0; i<=k-j; i++)
72 for (
unsigned int alpha=0; alpha<i; alpha++)
73 out[n] *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
74 for (
unsigned int beta=0; beta<j; beta++)
75 out[n] *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
76 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
77 out[n] *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
85 std::vector<typename Traits::JacobianType>& out)
const 91 out[0][0][0] = 0; out[0][0][1] = 0;
96 for (
unsigned int j=0; j<=k; j++)
97 for (
unsigned int i=0; i<=k-j; i++)
102 for (
unsigned int beta=0; beta<j; beta++)
103 factor *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
104 for (
unsigned int a=0; a<i; a++)
107 for (
unsigned int alpha=0; alpha<i; alpha++)
109 product *= D(1)/(pos_[i]-pos_[alpha]);
111 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
112 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
113 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
114 out[n][0][0] += product;
116 for (
unsigned int c=i+j+1; c<=k; c++)
119 for (
unsigned int alpha=0; alpha<i; alpha++)
120 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
121 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
123 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
125 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
126 out[n][0][0] += product;
132 for (
unsigned int alpha=0; alpha<i; alpha++)
133 factor *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
134 for (
unsigned int b=0; b<j; b++)
137 for (
unsigned int beta=0; beta<j; beta++)
139 product *= D(1)/(pos_[j]-pos_[beta]);
141 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
142 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
143 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
144 out[n][0][1] += product;
146 for (
unsigned int c=i+j+1; c<=k; c++)
149 for (
unsigned int beta=0; beta<j; beta++)
150 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
151 for (
unsigned int gamma=i+j+1; gamma<=k; gamma++)
153 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
155 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
156 out[n][0][1] += product;
171 std::vector<typename Traits::RangeType>& out)
const 173 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
182 std::array<int,1> directions;
183 directions[0] = std::find(order.begin(), order.end(), 1)-order.begin();
184 DUNE_NO_DEPRECATED_BEGIN
185 evaluate<1>(directions, in, out);
186 DUNE_NO_DEPRECATED_END
191 std::array<int,2> directions;
192 unsigned int counter = 0;
193 auto nonconstOrder =
order;
194 for (
int i=0; i<2; i++)
196 while (nonconstOrder[i])
198 directions[counter++] = i;
203 DUNE_NO_DEPRECATED_BEGIN
204 evaluate<2>(directions, in, out);
205 DUNE_NO_DEPRECATED_END
209 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
214 template<
unsigned int dOrder>
215 inline void DUNE_DEPRECATED_MSG(
"Use method 'partial' instead!")
217 const typename Traits::DomainType& in,
218 std::vector<typename Traits::RangeType>& out) const
223 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
230 for (
unsigned int j=0; j<=k; j++)
231 for (
unsigned int i=0; i<=k-j; i++, n++)
234 for (
unsigned int no1=0; no1 < k; no1++)
236 R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
237 for (
unsigned int no2=0; no2 < k; no2++)
239 factor *= lagrangianFactor(no2, i, j, in);
249 std::fill(out.begin(), out.end(), 0.0);
255 for (
unsigned int j=0; j<=k; j++)
256 for (
unsigned int i=0; i<=k-j; i++, n++)
260 for (
unsigned int no1=0; no1 < k; no1++)
262 R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
263 for (
unsigned int no2=0; no2 < k; no2++)
267 R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
268 for (
unsigned int no3=0; no3 < k; no3++)
270 if (no3 == no1 || no3 == no2)
272 factor2 *= lagrangianFactor(no3, i, j, in);
294 return (x[0]-pos_[no])/(pos_[i]-pos_[no]);
296 return (x[1]-pos_[no-i])/(pos_[j]-pos_[no-i]);
297 return (pos_[no+1]-x[0]-x[1])/(pos_[no+1]-pos_[i]-pos_[j]);
306 return (direction == 0) ? 1.0/(pos_[i]-pos_[no]) : 0;
309 return (direction == 0) ? 0: 1.0/(pos_[j]-pos_[no-i]);
311 return -1.0/(pos_[no+1]-pos_[i]-pos_[j]);
R RangeType
range type
Definition: localbasis.hh:61
unsigned int size() const
number of shape functions
Definition: pk2dlocalbasis.hh:51
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition: pk2dlocalbasis.hh:41
Definition: pk2dlocalbasis.hh:33
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
number of partial derivatives supported
Definition: localbasis.hh:74
D DomainType
domain type
Definition: localbasis.hh:49
Definition: pk2dlocalbasis.hh:38
Lagrange shape functions of arbitrary order on the reference triangle.
Definition: pk2dlocalbasis.hh:28
Pk2DLocalBasis()
Standard constructor.
Definition: pk2dlocalbasis.hh:44
unsigned int order() const
Polynomial order of the shape functions.
Definition: pk2dlocalbasis.hh:284
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pk2dlocalbasis.hh:84
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:169
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:216
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:37
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pk2dlocalbasis.hh:57