utility.hh

#ifndef DUNE_PDELAB_FINITEELEMENTMAP_UTILITY_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_UTILITY_HH
 
#include <cstddef>
 
#include <dune/common/std/type_traits.hh>
#include <dune/common/keywords.hh>
 
#include <dune/geometry/type.hh>
 
namespace Dune {
  namespace PDELab {
 
 
    template<typename FEM>
    using StaticFEMSize = decltype(FEM::size(GeometryTypes::vertex));
 
#ifndef DOXYGEN
 
    namespace Impl {
 
      // This function iterates over all geometry types up to the dimension of the finite element map
      // and returns the value of FEM::size(gt) iff that number is constant for all geometry types for
      // which the returned size is > 0. Otherwise it returns 0. As this only works if FEM::size() is
      // static, we use the additional argument to provide a separate overload if it is not.
      // Note that as there is no way to easily construct the set of "valid" geometry types for a
      // given dimension, we manually iterate over all possible topology ids. This creates weird
      // geometry types, but we just assume that FEM::size() will return 0 for invalid ones.
      template<typename FEM>
      constexpr std::size_t _femBlockSize(std::true_type)
      {
        constexpr int dim = FEM::dimension;
        std::size_t size = 0;
        for (int d = 0 ; d <= dim ; ++d)
          {
            std::size_t gt_size = FEM::size(GeometryTypes::none(d));
            if (gt_size > 0)
            {
              if (size > 0 and size != gt_size)
                return 0;
              else
                size = gt_size;
            }
            for (unsigned int topology_id = 0 ; topology_id < (1 << dim) ; ++topology_id)
              {
                std::size_t gt_size = FEM::size(GeometryType(topology_id,d));
                if (gt_size > 0)
                {
                  if (size > 0 and size != gt_size)
                    return 0;
                  else
                    size = gt_size;
                }
              }
          }
        return size;
      }
 
      // fallback version if `FEM::size()` is an instance method.
      template<typename FEM>
      constexpr std::size_t _femBlockSize(std::false_type)
      {
        return 0;
      }
 
    } // namespace Impl
 
#endif // DOXYGEN
 
 
    template<typename FEM>
    constexpr std::size_t finiteElementMapBlockSize()
    {
      return Impl::_femBlockSize<FEM>(Std::is_detected<StaticFEMSize,FEM>());
    }
 
    template<typename FEM>
    using FiniteElementMapBlockSize = std::integral_constant<std::size_t,finiteElementMapBlockSize<FEM>()>;
 
  } // namespace PDELab
} //namespace Dune
 
#endif // DUNE_PDELAB_FINITEELEMENTMAP_UTILITY_HH