5#ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
6#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGECUBE_HH
15#include <dune/geometry/referenceelements.hh>
17#include <dune/localfunctions/common/localbasis.hh>
18#include <dune/localfunctions/common/localfiniteelementtraits.hh>
19#include <dune/localfunctions/common/localkey.hh>
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 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);
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));
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];
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
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());
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());
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);
752 Impl::LagrangeCubeLocalBasis<D,R,dim,k> basis_;
753 Impl::LagrangeCubeLocalCoefficients<dim,k> coefficients_;
754 Impl::LagrangeCubeLocalInterpolation<Impl::LagrangeCubeLocalBasis<D,R,dim,k> > interpolation_;
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
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
Default exception for dummy implementations.
Definition: exceptions.hh:263
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: interpolate.hh:54
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
static const ReferenceElement & cube()
get hypercube reference elements
Definition: referenceelements.hh:168
D DomainType
domain type
Definition: 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