quadraturerules.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_GEOMETRY_QUADRATURERULES_HH
#define DUNE_GEOMETRY_QUADRATURERULES_HH

#include <algorithm>
#include <iostream>
#include <limits>
#include <mutex>
#include <utility>
#include <vector>

#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/stdstreams.hh>
#include <dune/common/stdthread.hh>
#include <dune/common/visibility.hh>

#include <dune/geometry/type.hh>
#include <dune/geometry/typeindex.hh>

namespace Dune {
  // forward declaration
  template<typename ct, int dim>
  class QuadraturePoint;
}

// class specialization of standard classes that allow to use structured bindings on QuadraturePoint
namespace std {
  template<typename ct, int dim>
  struct tuple_size<Dune::QuadraturePoint<ct,dim>> : public std::integral_constant<std::size_t,2> {};

  template<typename ct, int dim>
  struct tuple_element<0, Dune::QuadraturePoint<ct,dim>> { using type = Dune::FieldVector<ct, dim>; };

  template<typename ct, int dim>
  struct tuple_element<1, Dune::QuadraturePoint<ct,dim>> { using type = ct; };
}

namespace Dune {

  class QuadratureOrderOutOfRange : public NotImplemented {};

  template<typename ct, int dim>
  class QuadraturePoint {
  public:
    constexpr static int dimension = dim;

    typedef ct Field;

    typedef Dune::FieldVector<ct,dim> Vector;

    QuadraturePoint (const Vector& x, ct w) : local(x), weight_(w)
    {}

    const Vector& position () const
    {
      return local;
    }

    const ct &weight () const
    {
      return weight_;
    }

    template<std::size_t index, std::enable_if_t<(index<=1), int> = 0>
    std::tuple_element_t<index, QuadraturePoint<ct, dim>> get() const
    {
      if constexpr (index == 0) {
        return local;
      }
      else {
        return weight_;
      }
    }

  protected:
    FieldVector<ct, dim> local;
    ct weight_;
  };

  namespace QuadratureType {
    enum Enum {
      GaussLegendre = 0,

      GaussJacobi_1_0 = 1,

      GaussJacobi_2_0 = 2,

      GaussJacobi_n_0 = 3,

      GaussLobatto = 4,

      GaussRadauLeft = 5,

      GaussRadauRight = 6,
      size
    };
  }

  template<typename ct, int dim>
  class QuadratureRule : public std::vector<QuadraturePoint<ct,dim> >
  {
  public:
    QuadratureRule() : delivered_order(-1) {}

  protected:
    QuadratureRule(GeometryType t) : geometry_type(t), delivered_order(-1) {}

    QuadratureRule(GeometryType t, int order) : geometry_type(t), delivered_order(order) {}
  public:
    constexpr static int d = dim;

    typedef ct CoordType;

    virtual int order () const { return delivered_order; }

    virtual GeometryType type () const { return geometry_type; }
    virtual ~QuadratureRule(){}

    typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator iterator;

  protected:
    GeometryType geometry_type;
    int delivered_order;
  };

  // Forward declaration of the factory class,
  // needed internally by the QuadratureRules container class.
  template<typename ctype, int dim> class QuadratureRuleFactory;

  template<typename ctype, int dim>
  class QuadratureRules {

    using QuadratureRule = Dune::QuadratureRule<ctype, dim>;

    // indexed by quadrature order
    using QuadratureOrderVector = std::vector<std::pair<std::once_flag, QuadratureRule> >;

    // indexed by geometry type
    using GeometryTypeVector = std::vector<std::pair<std::once_flag, QuadratureOrderVector> >;

    // indexed by quadrature type enum
    using QuadratureCacheVector = std::vector<std::pair<std::once_flag, GeometryTypeVector> >;

    DUNE_EXPORT const QuadratureRule& _rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre) const
    {
      assert(t.dim()==dim);

      DUNE_ASSERT_CALL_ONCE();

      static QuadratureCacheVector quadratureCache(QuadratureType::size);

      auto& [ onceFlagQuadratureType, geometryTypes ] = quadratureCache[qt];
      // initialize geometry types for this quadrature type once
      std::call_once(onceFlagQuadratureType, [&types = geometryTypes]{
        types = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
      });

      auto& [ onceFlagGeometryType, quadratureOrders ] = geometryTypes[LocalGeometryTypeIndex::index(t)];
      // initialize quadrature orders for this geometry type and quadrature type once
      std::call_once(onceFlagGeometryType, [&, &orders = quadratureOrders]{
        // we only need one quadrature rule for points, not maxint
        const auto numRules = dim == 0 ? 1 : QuadratureRuleFactory<ctype,dim>::maxOrder(t, qt)+1;
        orders = QuadratureOrderVector(numRules);
      });

      // we only have one quadrature rule for points
      auto& [ onceFlagQuadratureOrder, quadratureRule ] = quadratureOrders[dim == 0 ? 0 : p];
      // initialize quadrature rule once
      std::call_once(onceFlagQuadratureOrder, [&, &rule = quadratureRule]{
        rule = QuadratureRuleFactory<ctype,dim>::rule(t, p, qt);
      });

      return quadratureRule;
    }

    DUNE_EXPORT static QuadratureRules& instance()
    {
      static QuadratureRules instance;
      return instance;
    }

    QuadratureRules () = default;
  public:
    static unsigned
    maxOrder(const GeometryType& t,
             QuadratureType::Enum qt=QuadratureType::GaussLegendre)
    {
      return QuadratureRuleFactory<ctype,dim>::maxOrder(t,qt);
    }

    static const QuadratureRule& rule(const GeometryType& t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
    {
      return instance()._rule(t,p,qt);
    }

    static const QuadratureRule& rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
    {
      GeometryType gt(t,dim);
      return instance()._rule(gt,p,qt);
    }
  };

} // end namespace Dune

#define DUNE_INCLUDING_IMPLEMENTATION

// 0d rules
#include "quadraturerules/pointquadrature.hh"
// 1d rules
#include "quadraturerules/gausslobattoquadrature.hh"
#include "quadraturerules/gaussquadrature.hh"
#include "quadraturerules/gaussradauleftquadrature.hh"
#include "quadraturerules/gaussradaurightquadrature.hh"
#include "quadraturerules/jacobi1quadrature.hh"
#include "quadraturerules/jacobi2quadrature.hh"
#include "quadraturerules/jacobiNquadrature.hh"
// 3d rules
#include "quadraturerules/prismquadrature.hh"
// general rules
#include "quadraturerules/simplexquadrature.hh"
#include "quadraturerules/tensorproductquadrature.hh"

#undef DUNE_INCLUDING_IMPLEMENTATION

namespace Dune {

  template<typename ctype, int dim>
  class QuadratureRuleFactory {
  private:
    friend class QuadratureRules<ctype, dim>;
    static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
    {
      return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
    }
    static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
    {
      return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
    }
  };

  template<typename ctype>
  class QuadratureRuleFactory<ctype, 0> {
  private:
    constexpr static int dim = 0;
    friend class QuadratureRules<ctype, dim>;
    static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum)
    {
      if (t.isVertex())
      {
        return std::numeric_limits<int>::max();
      }
      DUNE_THROW(Exception, "Unknown GeometryType");
    }
    static QuadratureRule<ctype, dim> rule(const GeometryType& t, int , QuadratureType::Enum)
    {
      if (t.isVertex())
      {
        return PointQuadratureRule<ctype>();
      }
      DUNE_THROW(Exception, "Unknown GeometryType");
    }
  };

  template<typename ctype>
  class QuadratureRuleFactory<ctype, 1> {
  private:
    constexpr static int dim = 1;
    friend class QuadratureRules<ctype, dim>;
    static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
    {
      if (t.isLine())
      {
        switch (qt) {
        case QuadratureType::GaussLegendre :
          return GaussQuadratureRule1D<ctype>::highest_order;
        case QuadratureType::GaussJacobi_1_0 :
          return Jacobi1QuadratureRule1D<ctype>::highest_order;
        case QuadratureType::GaussJacobi_2_0 :
          return Jacobi2QuadratureRule1D<ctype>::highest_order;
        case QuadratureType::GaussLobatto :
          return GaussLobattoQuadratureRule1D<ctype>::highest_order;
        case QuadratureType::GaussJacobi_n_0 :
          return JacobiNQuadratureRule1D<ctype>::maxOrder();
        case QuadratureType::GaussRadauLeft :
          return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
        case QuadratureType::GaussRadauRight :
          return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
        default :
          DUNE_THROW(Exception, "Unknown QuadratureType");
        }
      }
      DUNE_THROW(Exception, "Unknown GeometryType");
    }
    static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
    {
      if (t.isLine())
      {
        switch (qt) {
        case QuadratureType::GaussLegendre :
          return GaussQuadratureRule1D<ctype>(p);
        case QuadratureType::GaussJacobi_1_0 :
          return Jacobi1QuadratureRule1D<ctype>(p);
        case QuadratureType::GaussJacobi_2_0 :
          return Jacobi2QuadratureRule1D<ctype>(p);
        case QuadratureType::GaussLobatto :
          return GaussLobattoQuadratureRule1D<ctype>(p);
        case QuadratureType::GaussJacobi_n_0 :
          return JacobiNQuadratureRule1D<ctype>(p);
        case QuadratureType::GaussRadauLeft :
          return GaussRadauLeftQuadratureRule1D<ctype>(p);
        case QuadratureType::GaussRadauRight :
          return GaussRadauRightQuadratureRule1D<ctype>(p);
        default :
          DUNE_THROW(Exception, "Unknown QuadratureType");
        }
      }
      DUNE_THROW(Exception, "Unknown GeometryType");
    }
  };

  template<typename ctype>
  class QuadratureRuleFactory<ctype, 2> {
  private:
    constexpr static int dim = 2;
    friend class QuadratureRules<ctype, dim>;
    static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
    {
      unsigned order =
        TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
      if (t.isSimplex())
        order = std::max
          (order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
      return order;
    }
    static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
    {
      if (t.isSimplex()
        && ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
        && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
      {
        return SimplexQuadratureRule<ctype,dim>(p);
      }
      return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
    }
  };

  template<typename ctype>
  class QuadratureRuleFactory<ctype, 3> {
  private:
    constexpr static int dim = 3;
    friend class QuadratureRules<ctype, dim>;
    static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt)
    {
      unsigned order =
        TensorProductQuadratureRule<ctype,dim>::maxOrder(t.id(), qt);
      if (t.isSimplex())
        order = std::max
          (order, static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
      if (t.isPrism())
        order = std::max
          (order, static_cast<unsigned>(PrismQuadratureRule<ctype,dim>::highest_order));
      return order;
    }
    static QuadratureRule<ctype, dim> rule(const GeometryType& t, int p, QuadratureType::Enum qt)
    {

      if (t.isSimplex()
        && ( qt == QuadratureType::GaussLegendre || qt == QuadratureType::GaussJacobi_n_0 )
        && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
      {
        return SimplexQuadratureRule<ctype,dim>(p);
      }
      if (t.isPrism()
        && qt == QuadratureType::GaussLegendre
        && p <= PrismQuadratureRule<ctype,dim>::highest_order)
      {
        return PrismQuadratureRule<ctype,dim>(p);
      }
      return TensorProductQuadratureRule<ctype,dim>(t.id(), p, qt);
    }
  };

#ifndef DUNE_NO_EXTERN_QUADRATURERULES
  extern template class GaussLobattoQuadratureRule<double, 1>;
  extern template class GaussQuadratureRule<double, 1>;
  extern template class GaussRadauLeftQuadratureRule<double, 1>;
  extern template class GaussRadauRightQuadratureRule<double, 1>;
  extern template class Jacobi1QuadratureRule<double, 1>;
  extern template class Jacobi2QuadratureRule<double, 1>;
  extern template class JacobiNQuadratureRule<double, 1>;
  extern template class PrismQuadratureRule<double, 3>;
  extern template class SimplexQuadratureRule<double, 2>;
  extern template class SimplexQuadratureRule<double, 3>;
#endif // !DUNE_NO_EXTERN_QUADRATURERULES

} // end namespace

#endif // DUNE_GEOMETRY_QUADRATURERULES_HH