refinedlagrangebasis.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_REFINEDLAGRANGEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_REFINEDLAGRANGEBASIS_HH
 
#include <type_traits>
 
#include <dune/common/exceptions.hh>
#include <dune/common/math.hh>
 
#include <dune/localfunctions/refined.hh>
 
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
#include <dune/functions/functionspacebases/nodes.hh>
 
#include <dune/geometry/type.hh>
 
#include <dune/grid/common/mcmgmapper.hh>
 
 
namespace Dune {
namespace Functions {
 
template<typename GV, int k, typename R>
class RefinedLagrangeNode;
 
template <typename GV, int k, typename R = double>
class RefinedLagrangePreBasis :
  public LeafPreBasisMapperMixin< GV >
{
  using Base = LeafPreBasisMapperMixin< GV >;
 
  static const int dim = GV::dimension;
 
  // refined basis only implemented for P0 and P1
  static_assert(k == 0 || k == 1);
 
  // the layout is defined in terms of a MCMGLayout specialized for k == 0 or 1
  static MCMGLayout dofLayout()
  {
    if constexpr(k == 0)
      // a refined P0 basis assigns each element 2^dim DOFs
      return [](GeometryType gt, int) -> size_t {
        return (gt.dim() == dim) ? (1 << dim) : 0;
      };
    else if constexpr(k == 1)
      // a refined P1 basis has the same layout as a P2 basis
      return [](GeometryType gt, int) -> size_t {
        return Dune::binomial(int(k),int(gt.dim()));
      };
    else
      DUNE_THROW(Dune::NotImplemented,
        "Refined basis not implemented for higher-order Lagrange (k>=2) elements.");
  }
 
public:
 
  using GridView = GV;
 
  using Node = RefinedLagrangeNode<GV, k, R>;
 
  RefinedLagrangePreBasis (const GridView& gv)
    : Base(gv, dofLayout())
  {
    for (auto gt : gv.indexSet().types(0)) {
      if (!gt.isSimplex())
        DUNE_THROW(Dune::NotImplemented,
          "Refined Lagrange basis only implemented for simplex grids.");
    }
  }
 
  Node makeNode () const
  {
    return Node{};
  }
 
  static constexpr unsigned int order()
  {
    return k;
  }
};
 
 
 
template <typename GV, int k, typename R>
class RefinedLagrangeNode
  : public LeafBasisNode
{
  static constexpr int dim = GV::dimension;
 
  // refined basis only implemented for P0 and P1
  static_assert(k == 0 || k == 1);
 
public:
  using Element = typename GV::template Codim<0>::Entity;
 
  using FiniteElement = std::conditional_t<(k==0),
    Dune::RefinedP0LocalFiniteElement<typename GV::ctype,R,dim>,
    Dune::RefinedP1LocalFiniteElement<typename GV::ctype,R,dim>>;
 
  RefinedLagrangeNode ()
    : finiteElement_{}
    , element_(nullptr)
  {}
 
  const Element& element () const
  {
    return *element_;
  }
 
  const FiniteElement& finiteElement () const
  {
    return finiteElement_;
  }
 
  void bind (const Element& e)
  {
    element_ = &e;
    this->setSize(finiteElement_.size());
  }
 
  static constexpr unsigned int order()
  {
    return k;
  }
 
protected:
  const FiniteElement finiteElement_;
  const Element* element_;
};
 
 
 
namespace BasisFactory {
 
template <int k, typename R=double>
auto refinedLagrange ()
{
  return [](const auto& gridView) {
    return RefinedLagrangePreBasis<std::decay_t<decltype(gridView)>, k, R>(gridView);
  };
}
 
} // end namespace BasisFactory
 
 
template <typename GV, int k, typename R=double>
using RefinedLagrangeBasis = DefaultGlobalBasis<RefinedLagrangePreBasis<GV,k,R> >;
 
} // end namespace Functions
} // end namespace Dune
 
 
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_REFINEDLAGRANGEBASIS_HH