│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │ -
5#ifndef DUNE_REFINED_P0_LOCALFINITEELEMENT_HH
│ │ │ │ -
6#define DUNE_REFINED_P0_LOCALFINITEELEMENT_HH
│ │ │ │ +
5#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
│ │ │ │ +
6#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
│ │ │ │
│ │ │ │ -
8#include <dune/geometry/type.hh>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
11#include <dune/common/fmatrix.hh>
│ │ │ │ +
12#include <dune/common/fvector.hh>
│ │ │ │ +
13#include <dune/common/math.hh>
│ │ │ │ +
│ │ │ │ +
15#include <dune/geometry/referenceelements.hh>
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
27 template<
class D,
class R,
int dim>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
42 template<
class D,
class R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
76 return interpolation_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
87 static constexpr GeometryType
type ()
│ │ │ │ -
│ │ │ │ -
89 return GeometryTypes::line;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
21namespace Dune {
namespace Impl
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
24 template<
class LocalBasis>
│ │ │ │ +
25 class LagrangeCubeLocalInterpolation;
│ │ │ │ +
│ │ │ │ +
37 template<
class D,
class R,
unsigned int dim,
unsigned int k>
│ │ │ │ +
38 class LagrangeCubeLocalBasis
│ │ │ │ +
│ │ │ │ +
40 friend class LagrangeCubeLocalInterpolation<LagrangeCubeLocalBasis<D,R,dim,k> >;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
43 static R p(
unsigned int i, D x)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
46 for (
unsigned int j=0; j<=k; j++)
│ │ │ │ +
47 if (j!=i) result *= (k*x-j)/((
int)i-(int)j);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
52 static R dp(
unsigned int i, D x)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
56 for (
unsigned int j=0; j<=k; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
60 R prod( (k*1.0)/((
int)i-(
int)j) );
│ │ │ │ +
61 for (
unsigned int l=0; l<=k; l++)
│ │ │ │ +
│ │ │ │ +
63 prod *= (k*x-l)/((
int)i-(int)l);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
72 static R ddp(
unsigned int j, D x)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
76 for (
unsigned int i=0; i<=k; i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
83 for (
unsigned int m=0; m<=k; m++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
88 R prod( (k*1.0)/((
int)j-(
int)m) );
│ │ │ │ +
89 for (
unsigned int l=0; l<=k; l++)
│ │ │ │ +
90 if (l!=i && l!=j && l!=m)
│ │ │ │ +
91 prod *= (k*x-l)/((
int)j-(int)l);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
95 result += sum * ( (k*1.0)/((
int)j-(int)i) );
│ │ │ │ +
│ │ │ │
│ │ │ │ -
102 template<
class D,
class R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
102 static std::array<unsigned int,dim> multiindex (
unsigned int i)
│ │ │ │ +
│ │ │ │ +
104 std::array<unsigned int,dim> alpha;
│ │ │ │ +
105 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
107 alpha[j] = i % (k+1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
129 return coefficients_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
136 return interpolation_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
142 return basis_.size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
147 static constexpr GeometryType
type ()
│ │ │ │ -
│ │ │ │ -
149 return GeometryTypes::triangle;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
114 using Traits = LocalBasisTraits<D,dim,FieldVector<D,dim>,R,1,FieldVector<R,1>,FieldMatrix<R,1,dim> >;
│ │ │ │ +
│ │ │ │ +
118 static constexpr unsigned int size ()
│ │ │ │ +
│ │ │ │ +
120 return power(k+1, dim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
125 std::vector<typename Traits::RangeType>& out)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
138 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
142 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
144 out[i] *= (i & (1<<j)) ? x[j] : 1-x[j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
150 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
153 std::array<unsigned int,dim> alpha(multiindex(i));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
162 template<
class D,
class R>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
189 return coefficients_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
196 return interpolation_;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
202 return basis_.size();
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
159 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
160 out[i] *= p(alpha[j],x[j]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
170 std::vector<typename Traits::JacobianType>& out)
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
177 std::fill(out[0][0].begin(), out[0][0].end(), 0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
185 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
188 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
192 out[i][0][j] = (i & (1<<j)) ? 1 : -1;
│ │ │ │ +
│ │ │ │ +
194 for (
unsigned int l=0; l<dim; l++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
198 out[i][0][j] *= (i & (1<<l)) ? x[l] : 1-x[l];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
207 static constexpr GeometryType
type ()
│ │ │ │ -
│ │ │ │ -
209 return GeometryTypes::tetrahedron;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
208 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
211 std::array<unsigned int,dim> alpha(multiindex(i));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
214 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
218 out[i][0][j] = dp(alpha[j],x[j]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
221 for (
unsigned int l=0; l<dim; l++)
│ │ │ │ +
│ │ │ │ +
223 out[i][0][j] *= p(alpha[l],x[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
234 void partial(
const std::array<unsigned int,dim>& order,
│ │ │ │ +
│ │ │ │ +
236 std::vector<typename Traits::RangeType>& out)
const
│ │ │ │ +
│ │ │ │ +
238 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
244 out[0] = (totalOrder==0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
252 evaluateFunction(in, out);
│ │ │ │ +
│ │ │ │ +
254 else if (totalOrder == 1)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
258 auto direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
│ │ │ │ +
259 if (direction >= dim)
│ │ │ │ +
260 DUNE_THROW(RangeError,
"Direction of partial derivative not found!");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
263 for (std::size_t i = 0; i < size(); ++i)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
267 out[i] = (i & (1<<direction)) ? 1 : -1;
│ │ │ │ +
│ │ │ │ +
269 for (
unsigned int j = 0; j < dim; ++j)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
273 out[i] *= (i & (1<<j)) ? in[j] : 1-in[j];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
277 else if (totalOrder == 2)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
280 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
283 std::array<unsigned int,dim> alpha(multiindex(i));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
289 for (std::size_t l=0; l<dim; l++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
294 out[i][0] *= p(alpha[l],in[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
298 out[i][0] *= dp(alpha[l],in[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
304 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
310 DUNE_THROW(NotImplemented,
"Partial derivative of order " << totalOrder <<
" is not implemented!");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
318 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
321 std::array<unsigned int,dim> alpha(multiindex(i));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
327 for (std::size_t l=0; l<dim; l++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
332 out[i][0] *= p(alpha[l],in[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
335 out[i][0] *= dp(alpha[l],in[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
338 out[i][0] *= ddp(alpha[l],in[l]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
341 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
348 static constexpr unsigned int order ()
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
359 template<
unsigned int dim,
unsigned int k>
│ │ │ │ +
360 class LagrangeCubeLocalCoefficients
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
363 static std::array<unsigned int,dim> multiindex (
unsigned int i)
│ │ │ │ +
│ │ │ │ +
365 std::array<unsigned int,dim> alpha;
│ │ │ │ +
366 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
│ │ │ │ +
368 alpha[j] = i % (k+1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
375 void setup1d(std::vector<unsigned int>& subEntity)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
379 unsigned lastIndex=0;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
386 subEntity[lastIndex++] = 0;
│ │ │ │ +
387 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
388 subEntity[lastIndex++] = 0;
│ │ │ │ +
│ │ │ │ +
390 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
392 assert(power(k+1,dim)==lastIndex);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
395 void setup2d(std::vector<unsigned int>& subEntity)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
399 unsigned lastIndex=0;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
413 subEntity[lastIndex++] = 0;
│ │ │ │ +
414 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
415 subEntity[lastIndex++] = 2;
│ │ │ │ +
│ │ │ │ +
417 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
420 for (
unsigned e = 0; e < k - 1; ++e) {
│ │ │ │ +
421 subEntity[lastIndex++] = 0;
│ │ │ │ +
422 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
423 subEntity[lastIndex++] = 0;
│ │ │ │ +
424 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
428 subEntity[lastIndex++] = 2;
│ │ │ │ +
429 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
430 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
432 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
434 assert(power(k+1,dim)==lastIndex);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
437 void setup3d(std::vector<unsigned int>& subEntity)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
441 unsigned lastIndex=0;
│ │ │ │ +
│ │ │ │ +
443 const unsigned numIndices = power(k+1,dim);
│ │ │ │ +
444 const unsigned numFaceIndices = power(k+1,dim-1);
│ │ │ │ +
│ │ │ │ +
446 const unsigned numInnerEdgeDofs = k-1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
467 subEntity[lastIndex++] = 0;
│ │ │ │ +
468 for (
unsigned i = 0; i < numInnerEdgeDofs; ++i)
│ │ │ │ +
469 subEntity[lastIndex++] = 6;
│ │ │ │ +
│ │ │ │ +
471 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
474 for (
unsigned e = 0; e < numInnerEdgeDofs; ++e) {
│ │ │ │ +
475 subEntity[lastIndex++] = 4;
│ │ │ │ +
476 for (
unsigned i = 0; i < numInnerEdgeDofs; ++i)
│ │ │ │ +
477 subEntity[lastIndex++] = 4;
│ │ │ │ +
478 subEntity[lastIndex++] = 5;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
482 subEntity[lastIndex++] = 2;
│ │ │ │ +
483 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
484 subEntity[lastIndex++] = 7;
│ │ │ │ +
485 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
487 assert(numFaceIndices==lastIndex);
│ │ │ │ +
│ │ │ │ +
491 for(
unsigned f = 0; f < numInnerEdgeDofs; ++f) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
494 subEntity[lastIndex++] = 0;
│ │ │ │ +
495 for (
unsigned i = 0; i < numInnerEdgeDofs; ++i)
│ │ │ │ +
496 subEntity[lastIndex++] = 2;
│ │ │ │ +
497 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
500 for (
unsigned e = 0; e < numInnerEdgeDofs; ++e) {
│ │ │ │ +
501 subEntity[lastIndex++] = 0;
│ │ │ │ +
502 for (
unsigned i = 0; i < numInnerEdgeDofs; ++i)
│ │ │ │ +
503 subEntity[lastIndex++] = 0;
│ │ │ │ +
504 subEntity[lastIndex++] = 1;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
508 subEntity[lastIndex++] = 2;
│ │ │ │ +
509 for (
unsigned i = 0; i < numInnerEdgeDofs; ++i)
│ │ │ │ +
510 subEntity[lastIndex++] = 3;
│ │ │ │ +
511 subEntity[lastIndex++] = 3;
│ │ │ │ +
│ │ │ │ +
513 assert(lastIndex==(f+1+1)*numFaceIndices);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
518 subEntity[lastIndex++] = 4;
│ │ │ │ +
519 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
520 subEntity[lastIndex++] = 10;
│ │ │ │ +
521 subEntity[lastIndex++] = 5;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
524 for (
unsigned e = 0; e < k - 1; ++e) {
│ │ │ │ +
525 subEntity[lastIndex++] = 8;
│ │ │ │ +
526 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
527 subEntity[lastIndex++] = 5;
│ │ │ │ +
528 subEntity[lastIndex++] = 9;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
532 subEntity[lastIndex++] = 6;
│ │ │ │ +
533 for (
unsigned i = 0; i < k - 1; ++i)
│ │ │ │ +
534 subEntity[lastIndex++] = 11;
│ │ │ │ +
535 subEntity[lastIndex++] = 7;
│ │ │ │ +
│ │ │ │ +
537 assert(numIndices==lastIndex);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
542 LagrangeCubeLocalCoefficients ()
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
547 localKeys_[0] = LocalKey(0,0,0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
553 for (std::size_t i=0; i<size(); i++)
│ │ │ │ +
554 localKeys_[i] = LocalKey(i,dim,0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
561 std::vector<unsigned int> codim(size());
│ │ │ │ +
│ │ │ │ +
563 for (std::size_t i=0; i<codim.size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
569 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ +
570 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
571 if (mIdx[j]==0 or mIdx[j]==k)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
580 std::vector<unsigned int> index(size());
│ │ │ │ +
│ │ │ │ +
582 for (std::size_t i=0; i<size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
586 std::array<unsigned int,dim> mIdx = multiindex(i);
│ │ │ │ +
│ │ │ │ +
588 for (
int j=dim-1; j>=0; j--)
│ │ │ │ +
589 if (mIdx[j]>0 && mIdx[j]<k)
│ │ │ │ +
590 index[i] = (k-1)*index[i] + (mIdx[j]-1);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
594 std::vector<unsigned int> subEntity(size());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
609 DUNE_THROW(Dune::NotImplemented,
"LagrangeCubeLocalCoefficients for order " << k <<
" and dim == " << dim);
│ │ │ │ +
│ │ │ │ +
611 for (
size_t i=0; i<size(); i++)
│ │ │ │ +
612 localKeys_[i] = LocalKey(subEntity[i], codim[i], index[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
616 static constexpr std::size_t size ()
│ │ │ │ +
│ │ │ │ +
618 return power(k+1,dim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
622 const LocalKey& localKey (std::size_t i)
const
│ │ │ │ +
│ │ │ │ +
624 return localKeys_[i];
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
628 std::vector<LocalKey> localKeys_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
635 template<
class LocalBasis>
│ │ │ │ +
636 class LagrangeCubeLocalInterpolation
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
647 template<
typename F,
typename C>
│ │ │ │ +
648 void interpolate (
const F& f, std::vector<C>& out)
const
│ │ │ │ +
│ │ │ │ +
650 constexpr auto dim = LocalBasis::Traits::dimDomain;
│ │ │ │ +
651 constexpr auto k = LocalBasis::order();
│ │ │ │ +
652 using D =
typename LocalBasis::Traits::DomainFieldType;
│ │ │ │ +
│ │ │ │ +
654 typename LocalBasis::Traits::DomainType x;
│ │ │ │ +
│ │ │ │ +
656 out.resize(LocalBasis::size());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
661 auto center = ReferenceElements<D,dim>::cube().position(0,0);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
669 for (
unsigned int i=0; i<LocalBasis::size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
672 for (
int j=0; j<dim; j++)
│ │ │ │ +
673 x[j] = (i & (1<<j)) ? 1.0 : 0.0;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
681 for (
unsigned int i=0; i<LocalBasis::size(); i++)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
684 std::array<unsigned int,dim> alpha(LocalBasis::multiindex(i));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
687 for (
unsigned int j=0; j<dim; j++)
│ │ │ │ +
688 x[j] = (1.0*alpha[j])/k;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
707 template<
class D,
class R,
int dim,
int k>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
714 Impl::LagrangeCubeLocalCoefficients<dim,k>,
│ │ │ │ +
715 Impl::LagrangeCubeLocalInterpolation<Impl::LagrangeCubeLocalBasis<D,R,dim,k> > >;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
728 return coefficients_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
735 return interpolation_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
739 static constexpr std::size_t
size ()
│ │ │ │ +
│ │ │ │ +
741 return power(k+1,dim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
746 static constexpr GeometryType
type ()
│ │ │ │ +
│ │ │ │ +
748 return GeometryTypes::cube(dim);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
752 Impl::LagrangeCubeLocalBasis<D,R,dim,k> basis_;
│ │ │ │ +
753 Impl::LagrangeCubeLocalCoefficients<dim,k> coefficients_;
│ │ │ │ +
754 Impl::LagrangeCubeLocalInterpolation<Impl::LagrangeCubeLocalBasis<D,R,dim,k> > interpolation_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
Definition bdfmcube.hh:18
│ │ │ │ +
D DomainType
domain type
Definition common/localbasis.hh:43
│ │ │ │
traits helper struct
Definition localfiniteelementtraits.hh:13
│ │ │ │
LB LocalBasisType
Definition localfiniteelementtraits.hh:16
│ │ │ │
LC LocalCoefficientsType
Definition localfiniteelementtraits.hh:20
│ │ │ │
LI LocalInterpolationType
Definition localfiniteelementtraits.hh:24
│ │ │ │ -
Local finite element that is piecewise P0 on a once uniformly refined reference geometry.
Definition refinedp0.hh:29
│ │ │ │ -
P0LocalFiniteElement< D, R, dim >::Traits Traits
Definition refinedp0.hh:35
│ │ │ │ -
const Traits::LocalBasisType & localBasis() const
Definition refinedp0.hh:60
│ │ │ │ -
static constexpr GeometryType type()
Definition refinedp0.hh:87
│ │ │ │ -
RefinedP0LocalFiniteElement()
Definition refinedp0.hh:55
│ │ │ │ -
const Traits::LocalInterpolationType & localInterpolation() const
Definition refinedp0.hh:74
│ │ │ │ -
const Traits::LocalCoefficientsType & localCoefficients() const
Definition refinedp0.hh:67
│ │ │ │ -
LocalFiniteElementTraits< RefinedP0LocalBasis< D, R, 1 >, RefinedP0LocalCoefficients< 1 >, RefinedP0LocalInterpolation< RefinedP0LocalBasis< D, R, 1 > > > Traits
Definition refinedp0.hh:51
│ │ │ │ -
unsigned int size() const
Number of shape functions in this finite element.
Definition refinedp0.hh:80
│ │ │ │ -
unsigned int size() const
Number of shape functions in this finite element.
Definition refinedp0.hh:140
│ │ │ │ -
RefinedP0LocalFiniteElement()
Definition refinedp0.hh:115
│ │ │ │ -
static constexpr GeometryType type()
Definition refinedp0.hh:147
│ │ │ │ -
LocalFiniteElementTraits< RefinedP0LocalBasis< D, R, 2 >, RefinedP0LocalCoefficients< 2 >, RefinedP0LocalInterpolation< RefinedP0LocalBasis< D, R, 2 > > > Traits
Definition refinedp0.hh:111
│ │ │ │ -
const Traits::LocalCoefficientsType & localCoefficients() const
Definition refinedp0.hh:127
│ │ │ │ -
const Traits::LocalInterpolationType & localInterpolation() const
Definition refinedp0.hh:134
│ │ │ │ -
const Traits::LocalBasisType & localBasis() const
Definition refinedp0.hh:120
│ │ │ │ -
RefinedP0LocalFiniteElement()
Definition refinedp0.hh:175
│ │ │ │ -
const Traits::LocalBasisType & localBasis() const
Definition refinedp0.hh:180
│ │ │ │ -
const Traits::LocalInterpolationType & localInterpolation() const
Definition refinedp0.hh:194
│ │ │ │ -
static constexpr GeometryType type()
Definition refinedp0.hh:207
│ │ │ │ -
unsigned int size() const
Number of shape functions in this finite element.
Definition refinedp0.hh:200
│ │ │ │ -
const Traits::LocalCoefficientsType & localCoefficients() const
Definition refinedp0.hh:187
│ │ │ │ -
LocalFiniteElementTraits< RefinedP0LocalBasis< D, R, 3 >, RefinedP0LocalCoefficients< 3 >, RefinedP0LocalInterpolation< RefinedP0LocalBasis< D, R, 3 > > > Traits
Definition refinedp0.hh:171
│ │ │ │ -
Uniformly refined constant shape functions on a unit simplex in R^dim.
Definition refinedp0localbasis.hh:41
│ │ │ │ -
Layout map for RefinedP0 elements.
Definition refinedp0localcoefficients.hh:26
│ │ │ │ -
Definition refinedp0localinterpolation.hh:17
│ │ │ │ +
Lagrange finite element for cubes with arbitrary compile-time dimension and polynomial order.
Definition lagrangecube.hh:709
│ │ │ │ +
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition lagrangecube.hh:719
│ │ │ │ +
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition lagrangecube.hh:733
│ │ │ │ +
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition lagrangecube.hh:746
│ │ │ │ +
static constexpr std::size_t size()
The number of shape functions.
Definition lagrangecube.hh:739
│ │ │ │ +
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition lagrangecube.hh:726
│ │ │ │ +
│ │ │ │
1.9.8
│ │ │ │
│ │ │ │
│ │ │ │ ├── html2text {}
│ │ │ │ │ @@ -1,36 +1,29 @@
│ │ │ │ │ dune-localfunctions 2.10
│ │ │ │ │ Loading...
│ │ │ │ │ Searching...
│ │ │ │ │ No Matches
│ │ │ │ │ * _d_u_n_e
│ │ │ │ │ * _l_o_c_a_l_f_u_n_c_t_i_o_n_s
│ │ │ │ │ - * _o_r_t_h_o_n_o_r_m_a_l
│ │ │ │ │ -_C_l_a_s_s_e_s | _N_a_m_e_s_p_a_c_e_s | _F_u_n_c_t_i_o_n_s
│ │ │ │ │ -orthonormalcompute.hh File Reference
│ │ │ │ │ -#include