hierarchicalsimplexp2withelementbubble.hh

Go to the documentation of this file.

// -*- 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_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
#define DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH
 
#include <array>
#include <cassert>
#include <numeric>
#include <stdexcept>
#include <vector>
 
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/math.hh>
 
#include <dune/geometry/referenceelement.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/common/rangecast.hh>
 
namespace Dune
{
  template<class D, class R, int dim>
  class HierarchicalSimplexP2WithElementBubbleLocalBasis
  {
    template<class,int> friend class HierarchicalSimplexP2WithElementBubbleLocalInterpolation;
 
    using DomainType = FieldVector<D,dim>;
 
    using RangeType = FieldVector<R,1>;
 
    using JacobianType = FieldMatrix<R,1,dim>;
 
    // Number of vertices
    static constexpr int numVertices = dim+1;
 
    // Number of edges (or zero for dim==1)
    static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
 
    // helper function to evaluate the vertex basis functions
    template <class It>
    static constexpr It evaluateVertexFunctions (const DomainType& in, It outIt)
    {
      *outIt = 1;
      for (int i = 0; i < dim; ++i)
        *outIt -= in[i];
      ++outIt;
      for (int i = 0; i < dim; ++i)
        *outIt++ = in[i];
      return outIt;
    }
 
    // helper function to evaluate the basis functions
    template <class It>
    static constexpr It evaluateAllFunctions (const DomainType& in, It outIt)
    {
      It vertexValues = outIt;
      outIt = evaluateVertexFunctions(in, outIt);
 
      if constexpr(dim > 1) {
        auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
        for (int i = 0; i < numEdges; ++i) {
          const int v0 = refElem.subEntity(i,dim-1,0,dim);
          const int v1 = refElem.subEntity(i,dim-1,1,dim);
          *outIt++ = 4 * vertexValues[v0] * vertexValues[v1];
        }
      }
 
      // element bubble function
      *outIt = power(dim+1, dim+1);
      for (int i = 0; i < numVertices; ++i)
        *outIt *= Impl::rangeCast(vertexValues[i]);
      return outIt;
    }
 
  public:
    using Traits = LocalBasisTraits<D,dim,DomainType,R,1,RangeType,JacobianType>;
 
    static constexpr std::size_t size () noexcept
    {
      return numVertices + numEdges + 1;
    }
 
    static constexpr void evaluateFunction (const DomainType& in,
                                            std::vector<RangeType>& out)
    {
      out.resize(size());
      evaluateAllFunctions(in,out.begin());
    }
 
    static constexpr void evaluateJacobian (const DomainType& in,
                                            std::vector<JacobianType>& out)
    {
      out.resize(size());
 
      // vertex basis functions
      R tmp = 1;
      for (int i = 0; i < dim; ++i) {
        out[0][0][i] = -1;
        for (int j = 0; j < dim; ++j)
          out[j+1][0][i] = (i == j);
        tmp -= in[i];
      }
 
      int n = numVertices;
      std::array<RangeType,numVertices> shapeValues;
      evaluateVertexFunctions(in, shapeValues.begin());
 
      // edge basis functions
      if constexpr(dim > 1) {
        auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
        for (int i = 0; i < numEdges; ++i,++n) {
          const int v0 = refElem.subEntity(i,dim-1,0,dim);
          const int v1 = refElem.subEntity(i,dim-1,1,dim);
          for (int j = 0; j < dim; ++j)
            out[n][0][j] = 4 * (out[v0][0][j] * shapeValues[v1][0] + shapeValues[v0][0] * out[v1][0][j]);
        }
      }
 
      // element bubble function
      for (int i = 0; i < dim; ++i) {
        out[n][0][i] = power(dim+1, dim+1) * (tmp - in[i]);
        for (int j = i+1; j < dim+i; ++j)
          out[n][0][i] *= in[j % dim];
      }
    }
 
    static constexpr void partial (const std::array<unsigned int, dim>& order,
                                   const DomainType& in,
                                   std::vector<RangeType>& out)
    {
      unsigned int totalOrder = 0;
      for (int i = 0; i < dim; ++i)
        totalOrder += order[i];
 
      switch (totalOrder) {
      case 0:
        evaluateFunction(in,out);
        break;
      case 1: {
        out.resize(size());
        int d = 0; // the direction of differentiation
        for (int i = 0; i < dim; ++i)
          d += i * order[i];
 
        // vertex basis functions
        R tmp = 1;
        for (int i = 0; i < dim; ++i) {
          out[0] = -1;
          for (int j = 0; j < dim; ++j)
            out[j+1] = (dim == j);
          tmp -= in[i];
        }
 
        int n = numVertices;
        std::array<RangeType,numVertices> shapeValues;
        evaluateVertexFunctions(in, shapeValues.begin());
 
        // edge basis functions
        if constexpr(dim > 1) {
          auto refElem = referenceElement<D,dim>(GeometryTypes::simplex(dim));
          for (int i = 0; i < numEdges; ++i,++n) {
            const int v0 = refElem.subEntity(i,dim-1,0,dim);
            const int v1 = refElem.subEntity(i,dim-1,1,dim);
            out[n] = 4 * (out[v0] * shapeValues[v1] + shapeValues[v0] * out[v1]);
          }
        }
 
        // element bubble function
        out[n] = power(dim+1, dim+1) * (tmp - in[d]);
        for (int j = d+1; j < dim+d; ++j)
          out[n] *= in[j % dim];
      } break;
      default:
        throw std::runtime_error("Desired derivative order is not implemented");
      }
    }
 
    static constexpr unsigned int order () noexcept
    {
      return dim+1;
    }
  };
 
 
  template <int dim>
  class HierarchicalSimplexP2WithElementBubbleLocalCoefficients
  {
    // Number of vertices
    static constexpr int numVertices = dim+1;
 
    // Number of edges (or zero for dim==1)
    static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
 
  public:
    HierarchicalSimplexP2WithElementBubbleLocalCoefficients () noexcept
    {
      int n = 0;
      for (int i = 0; i < numVertices; ++i)
        li_[n++] = LocalKey(i,dim,0);       // Vertex
      if constexpr(dim > 1) {
        for (int i = 0; i < numEdges; ++i)
          li_[n++] = LocalKey(i,dim-1,0);   // Edges
      }
      li_[n++] = LocalKey(0,0,0);           // Element
    }
 
    static constexpr std::size_t size () noexcept
    {
      return numVertices + numEdges + 1;
    }
 
    const LocalKey& localKey (std::size_t i) const noexcept
    {
      return li_[i];
    }
 
  private:
    std::array<LocalKey, numVertices+numEdges+1> li_;
  };
 
  template<class LB, int dim>
  class HierarchicalSimplexP2WithElementBubbleLocalInterpolation
  {
    using LocalBasis = LB;
    using DomainType = typename LB::Traits::DomainType;
    using RangeType = typename LB::Traits::RangeType;
 
    // Number of vertices
    static constexpr int numVertices = dim+1;
 
    // Number of edges (or zero for dim==1)
    static constexpr int numEdges = (dim > 1 ? ((dim+1)*dim / 2) : 0);
 
  public:
    template<class F, class C>
    requires ( requires (const F& f, const DomainType& x, C& y){ y = Impl::rangeCast(f(x)); } )
    static constexpr void interpolate (const F& f, std::vector<C>& out)
    {
      auto refElem = referenceElement<typename LB::Traits::DomainFieldType,dim>(GeometryTypes::simplex(dim));
 
      out.resize(LB::size());
      int n = 0;
 
      // vertices
      assert(numVertices == refElem.size(dim));
      for (int i = 0; i < numVertices; ++i)
        out[n++] = Impl::rangeCast(f(refElem.position(i,dim)));
 
      std::array<RangeType,LB::size()> shapeValues;
 
      // edge bubbles
      if constexpr(dim > 1) {
        assert(numEdges == refElem.size(dim-1));
        for (int i = 0; i < numEdges; ++i) {
          auto y = Impl::rangeCast(f(refElem.position(i,dim-1)));
          LB::evaluateVertexFunctions(refElem.position(i,dim-1), shapeValues.begin());
          for (int j = 0; j < numVertices; ++j)
            y -= out[j]*shapeValues[j][0];
          out[n++] = y;
        }
      }
 
      // element bubble
      auto y = Impl::rangeCast(f(refElem.position(0,0)));
      LB::evaluateAllFunctions(refElem.position(0,0), shapeValues.begin());
      for (int j = 0; j < numVertices+numEdges; ++j)
        y -= out[j]*shapeValues[j][0];
      out[n++] = y;
    }
  };
 
} // end namespace Dune
 
#endif // DUNE_HIERARCHICAL_SIMPLEX_P2_WITH_ELEMENT_BUBBLE_LOCALBASIS_HH