localinterpolation.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_LOCALFUNCTIONS_ENRICHED_CUBEQ1BUBBLE_LOCALINTERPOLATION_HH
#define DUNE_LOCALFUNCTIONS_ENRICHED_CUBEQ1BUBBLE_LOCALINTERPOLATION_HH
 
#include <type_traits>
#include <vector>
 
namespace Dune
{
  template<class LB>
  class CubeQ1BubbleLocalInterpolation
  {
    static const int dim = LB::dimension;
    static const int numVertices = power(2, dim);
 
    using DomainType = typename LB::Traits::DomainType;
    using RangeType = typename LB::Traits::RangeType;
 
  public:
    template<class F, class C,
      class R = std::invoke_result_t<F, DomainType>,
      std::enable_if_t<std::is_convertible_v<R, C>, int> = 0>
    static constexpr void interpolate (const F& f, std::vector<C>& out)
    {
      out.resize(numVertices+1);
 
      // vertices
      DomainType x(0);
      for (int i = 0; i < numVertices; ++i)
      {
        // Generate coordinate of the i-th corner of the reference cube
        for (int j = 0; j < dim; ++j)
          x[j] = (i & (1<<j)) ? 1 : 0;
 
        out[i] = f(x);
      }
 
      // element bubble
      x = 0.5;
      R y = f(x);
 
      // evaluate the other shape functions in x and subtract this value
      std::vector<RangeType> sfValues;
      LB::evaluateFunction(x, sfValues);
 
      out[numVertices] = y;
      for (int i = 0; i < numVertices; ++i)
        out[numVertices] -= out[i]*sfValues[i];
    }
  };
 
} // end namespace Dune
 
#endif // DUNE_LOCALFUNCTIONS_ENRICHED_CUBEQ1BUBBLE_LOCALINTERPOLATION_HH