utilities.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2021 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_utilities_h
#define dealii_utilities_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
 
#include <functional>
#include <string>
#include <tuple>
#include <type_traits>
#include <typeinfo>
#include <utility>
#include <vector>
 
#ifdef DEAL_II_WITH_TRILINOS
#  include <Epetra_Comm.h>
#  include <Epetra_Map.h>
#  include <Teuchos_Comm.hpp>
#  include <Teuchos_RCP.hpp>
#  ifdef DEAL_II_WITH_MPI
#    include <Epetra_MpiComm.h>
#  else
#    include <Epetra_SerialComm.h>
#  endif
#endif
 
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
 
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/core/demangle.hpp>
#include <boost/iostreams/device/array.hpp>
#include <boost/iostreams/device/back_inserter.hpp>
#include <boost/iostreams/filtering_streambuf.hpp>
#include <boost/serialization/array.hpp>
#include <boost/serialization/complex.hpp>
#include <boost/serialization/vector.hpp>
 
#ifdef DEAL_II_WITH_ZLIB
#  include <boost/iostreams/filter/gzip.hpp>
#endif
 
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
DEAL_II_NAMESPACE_OPEN
 
// forward declare Point
#ifndef DOXYGEN
template <int dim, typename Number>
class Point;
#endif
 
namespace Utilities
{
  std::string
  dealii_version_string();
 
  template <int dim, typename Number>
  std::vector<std::array<std::uint64_t, dim>>
  inverse_Hilbert_space_filling_curve(
    const std::vector<Point<dim, Number>> &points,
    const int                              bits_per_dim = 64);
 
  template <int dim>
  std::vector<std::array<std::uint64_t, dim>>
  inverse_Hilbert_space_filling_curve(
    const std::vector<std::array<std::uint64_t, dim>> &points,
    const int                                          bits_per_dim = 64);
 
  template <int dim>
  std::uint64_t
  pack_integers(const std::array<std::uint64_t, dim> &index,
                const int                             bits_per_dim);
 
  std::string
  compress(const std::string &input);
 
  std::string
  decompress(const std::string &compressed_input);
 
  std::string
  encode_base64(const std::vector<unsigned char> &binary_input);
 
  std::vector<unsigned char>
  decode_base64(const std::string &base64_input);
 
  std::string
  int_to_string(const unsigned int value,
                const unsigned int digits = numbers::invalid_unsigned_int);
 
  template <typename number>
  std::string
  to_string(const number       value,
            const unsigned int digits = numbers::invalid_unsigned_int);
 
  unsigned int
  needed_digits(const unsigned int max_number);
 
  template <typename Number>
  Number
  truncate_to_n_digits(const Number number, const unsigned int n_digits);
 
  int
  string_to_int(const std::string &s);
 
  std::string
  dim_string(const int dim, const int spacedim);
 
  std::vector<int>
  string_to_int(const std::vector<std::string> &s);
 
  double
  string_to_double(const std::string &s);
 
 
  std::vector<double>
  string_to_double(const std::vector<std::string> &s);
 
 
  std::vector<std::string>
  split_string_list(const std::string &s, const std::string &delimiter = ",");
 
 
  std::vector<std::string>
  split_string_list(const std::string &s, const char delimiter);
 
 
  std::vector<std::string>
  break_text_into_lines(const std::string &original_text,
                        const unsigned int width,
                        const char         delimiter = ' ');
 
  bool
  match_at_string_start(const std::string &name, const std::string &pattern);
 
  std::pair<int, unsigned int>
  get_integer_at_position(const std::string &name, const unsigned int position);
 
  std::string
  replace_in_string(const std::string &input,
                    const std::string &from,
                    const std::string &to);
 
  std::string
  trim(const std::string &input);
 
  double
  generate_normal_random_number(const double a, const double sigma);
 
  template <class T>
  std::string
  type_to_string(const T &t);
 
  template <int N, typename T>
  T
  fixed_power(const T t);
 
  template <typename T>
  constexpr T
  pow(const T base, const int iexp)
  {
#if defined(DEBUG) && !defined(DEAL_II_CXX14_CONSTEXPR_BUG)
    // Up to __builtin_expect this is the same code as in the 'Assert' macro.
    // The call to __builtin_expect turns out to be problematic.
    if (!(iexp >= 0))
      ::deal_II_exceptions::internals::issue_error_noreturn(
        ::deal_II_exceptions::internals::abort_or_throw_on_exception,
        __FILE__,
        __LINE__,
        __PRETTY_FUNCTION__,
        "iexp>=0",
        "ExcMessage(\"The exponent must not be negative!\")",
        ExcMessage("The exponent must not be negative!"));
#endif
    // The "exponentiation by squaring" algorithm used below has to be
    // compressed to one statement due to C++11's restrictions on constexpr
    // functions. A more descriptive version would be:
    //
    // <code>
    // if (iexp <= 0)
    //   return 1;
    //
    // // avoid overflow of one additional recursion with pow(base * base, 0)
    // if (iexp == 1)
    //   return base;
    //
    // // if the current exponent is not divisible by two,
    // // we need to account for that.
    // const unsigned int prefactor = (iexp % 2 == 1) ? base : 1;
    //
    // // a^b = (a*a)^(b/2)      for b even
    // // a^b = a*(a*a)^((b-1)/2 for b odd
    // return prefactor * ::Utilities::pow(base*base, iexp/2);
    // </code>
 
    static_assert(std::is_integral<T>::value, "Only integral types supported");
 
    return iexp <= 0 ?
             1 :
             (iexp == 1 ? base :
                          (((iexp % 2 == 1) ? base : 1) *
                           ::Utilities::pow(base * base, iexp / 2)));
  }
 
  template <typename Iterator, typename T>
  Iterator
  lower_bound(Iterator first, Iterator last, const T &val);
 
 
  template <typename Iterator, typename T, typename Comp>
  Iterator
  lower_bound(Iterator first, Iterator last, const T &val, const Comp comp);
 
  template <typename Integer>
  std::vector<Integer>
  reverse_permutation(const std::vector<Integer> &permutation);
 
  template <typename Integer>
  std::vector<Integer>
  invert_permutation(const std::vector<Integer> &permutation);
 
  template <typename T>
  size_t
  pack(const T &          object,
       std::vector<char> &dest_buffer,
       const bool         allow_compression = true);
 
  template <typename T>
  std::vector<char>
  pack(const T &object, const bool allow_compression = true);
 
  template <typename T>
  T
  unpack(const std::vector<char> &buffer, const bool allow_compression = true);
 
  template <typename T>
  T
  unpack(const std::vector<char>::const_iterator &cbegin,
         const std::vector<char>::const_iterator &cend,
         const bool                               allow_compression = true);
 
  template <typename T, int N>
  void
  unpack(const std::vector<char> &buffer,
         T (&unpacked_object)[N],
         const bool allow_compression = true);
 
  template <typename T, int N>
  void
  unpack(const std::vector<char>::const_iterator &cbegin,
         const std::vector<char>::const_iterator &cend,
         T (&unpacked_object)[N],
         const bool allow_compression = true);
 
  bool
  get_bit(const unsigned char number, const unsigned int n);
 
 
  void
  set_bit(unsigned char &number, const unsigned int n, const bool x);
 
 
  template <typename To, typename From>
  std::unique_ptr<To>
  dynamic_unique_cast(std::unique_ptr<From> &&p);
 
  template <typename T>
  T &
  get_underlying_value(T &p);
 
  template <typename T>
  T &
  get_underlying_value(std::shared_ptr<T> &p);
 
  template <typename T>
  T &
  get_underlying_value(const std::shared_ptr<T> &p);
 
  template <typename T>
  T &
  get_underlying_value(std::unique_ptr<T> &p);
 
  template <typename T>
  T &
  get_underlying_value(const std::unique_ptr<T> &p);
 
  namespace System
  {
    double
    get_cpu_load();
 
    const std::string
    get_current_vectorization_level();
 
    struct MemoryStats
    {
      unsigned long int VmPeak;
 
      unsigned long int VmSize;
 
      unsigned long int VmHWM;
 
      unsigned long int VmRSS;
    };
 
 
    void
    get_memory_stats(MemoryStats &stats);
 
 
    std::string
    get_hostname();
 
 
    std::string
    get_time();
 
    std::string
    get_date();
 
    void
    posix_memalign(void **memptr, std::size_t alignment, std::size_t size);
  } // namespace System
 
 
#ifdef DEAL_II_WITH_TRILINOS
  namespace Trilinos
  {
    const Epetra_Comm &
    comm_world();
 
    const Epetra_Comm &
    comm_self();
 
    const Teuchos::RCP<const Teuchos::Comm<int>> &
    tpetra_comm_self();
 
    Epetra_Comm *
    duplicate_communicator(const Epetra_Comm &communicator);
 
    void
    destroy_communicator(Epetra_Comm &communicator);
 
    unsigned int
    get_n_mpi_processes(const Epetra_Comm &mpi_communicator);
 
    unsigned int
    get_this_mpi_process(const Epetra_Comm &mpi_communicator);
 
    Epetra_Map
    duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm);
  } // namespace Trilinos
 
#endif
 
 
} // namespace Utilities
 
 
// --------------------- inline functions
 
namespace Utilities
{
  template <int N, typename T>
  inline T
  fixed_power(const T x)
  {
    Assert(
      !std::is_integral<T>::value || (N >= 0),
      ExcMessage(
        "The non-type template parameter N must be a non-negative integer for integral type T"));
 
    if (N == 0)
      return T(1.);
    else if (N < 0)
      return T(1.) / fixed_power<-N>(x);
    else
      // Use exponentiation by squaring:
      return ((N % 2 == 1) ? x * fixed_power<N / 2>(x * x) :
                             fixed_power<N / 2>(x * x));
  }
 
 
 
  template <class T>
  inline std::string
  type_to_string(const T &t)
  {
    return boost::core::demangle(typeid(t).name());
  }
 
 
 
  template <typename Iterator, typename T>
  inline Iterator
  lower_bound(Iterator first, Iterator last, const T &val)
  {
    return Utilities::lower_bound(first, last, val, std::less<T>());
  }
 
 
 
  template <typename Iterator, typename T, typename Comp>
  inline Iterator
  lower_bound(Iterator first, Iterator last, const T &val, const Comp comp)
  {
    // verify that the two iterators are properly ordered. since
    // we need operator- for the iterator type anyway, do the
    // test as follows, rather than via 'last >= first'
    Assert(last - first >= 0,
           ExcMessage(
             "The given iterators do not satisfy the proper ordering."));
 
    unsigned int len = static_cast<unsigned int>(last - first);
 
    if (len == 0)
      return first;
 
    while (true)
      {
        // if length equals 8 or less,
        // then do a rolled out
        // search. use a switch without
        // breaks for that and roll-out
        // the loop somehow
        if (len < 8)
          {
            switch (len)
              {
                case 7:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 6:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 5:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 4:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 3:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 2:
                  if (!comp(*first, val))
                    return first;
                  ++first;
                  DEAL_II_FALLTHROUGH;
                case 1:
                  if (!comp(*first, val))
                    return first;
                  return first + 1;
                default:
                  // indices seem
                  // to not be
                  // sorted
                  // correctly!? or
                  // did len
                  // become==0
                  // somehow? that
                  // shouldn't have
                  // happened
                  Assert(false, ExcInternalError());
              }
          }
 
 
 
        const unsigned int half   = len >> 1;
        const Iterator     middle = first + half;
 
        // if the value is larger than
        // that pointed to by the
        // middle pointer, then the
        // insertion point must be
        // right of it
        if (comp(*middle, val))
          {
            first = middle + 1;
            len -= half + 1;
          }
        else
          len = half;
      }
  }
 
 
  // --------------------- non-inline functions
 
  template <typename T>
  size_t
  pack(const T &          object,
       std::vector<char> &dest_buffer,
       const bool         allow_compression)
  {
    std::size_t size = 0;
 
    // see if the object is small and copyable via memcpy. if so, use
    // this fast path. otherwise, we have to go through the BOOST
    // serialization machinery
#ifdef DEAL_II_HAVE_CXX17
    if constexpr (std::is_trivially_copyable<T>() && sizeof(T) < 256)
#else
    if (std::is_trivially_copyable<T>() && sizeof(T) < 256)
#endif
      {
        (void)allow_compression;
        const std::size_t previous_size = dest_buffer.size();
        dest_buffer.resize(previous_size + sizeof(T));
 
        std::memcpy(dest_buffer.data() + previous_size, &object, sizeof(T));
 
        size = sizeof(T);
      }
    else
      {
        // use buffer as the target of a compressing
        // stream into which we serialize the current object
        const std::size_t previous_size = dest_buffer.size();
        {
          boost::iostreams::filtering_ostreambuf fosb;
#ifdef DEAL_II_WITH_ZLIB
          if (allow_compression)
            fosb.push(boost::iostreams::gzip_compressor());
#else
          (void)allow_compression;
#endif
          fosb.push(boost::iostreams::back_inserter(dest_buffer));
 
          boost::archive::binary_oarchive boa(fosb);
          boa << object;
          // the stream object has to be destroyed before the return statement
          // to ensure that all data has been written in the buffer
        }
        size = dest_buffer.size() - previous_size;
      }
 
    return size;
  }
 
 
  template <typename T>
  std::vector<char>
  pack(const T &object, const bool allow_compression)
  {
    std::vector<char> buffer;
    pack<T>(object, buffer, allow_compression);
    return buffer;
  }
 
 
  template <typename T>
  T
  unpack(const std::vector<char>::const_iterator &cbegin,
         const std::vector<char>::const_iterator &cend,
         const bool                               allow_compression)
  {
    T object;
 
    // see if the object is small and copyable via memcpy. if so, use
    // this fast path. otherwise, we have to go through the BOOST
    // serialization machinery
#ifdef DEAL_II_HAVE_CXX17
    if constexpr (std::is_trivially_copyable<T>() && sizeof(T) < 256)
#else
    if (std::is_trivially_copyable<T>() && sizeof(T) < 256)
#endif
      {
        (void)allow_compression;
        Assert(std::distance(cbegin, cend) == sizeof(T), ExcInternalError());
        std::memcpy(&object, &*cbegin, sizeof(T));
      }
    else
      {
        // decompress the buffer section into the object
        boost::iostreams::filtering_istreambuf fisb;
#ifdef DEAL_II_WITH_ZLIB
        if (allow_compression)
          fisb.push(boost::iostreams::gzip_decompressor());
#else
        (void)allow_compression;
#endif
        fisb.push(boost::iostreams::array_source(&*cbegin, &*cend));
 
        boost::archive::binary_iarchive bia(fisb);
        bia >> object;
      }
 
    return object;
  }
 
 
  template <typename T>
  T
  unpack(const std::vector<char> &buffer, const bool allow_compression)
  {
    return unpack<T>(buffer.cbegin(), buffer.cend(), allow_compression);
  }
 
 
  template <typename T, int N>
  void
  unpack(const std::vector<char>::const_iterator &cbegin,
         const std::vector<char>::const_iterator &cend,
         T (&unpacked_object)[N],
         const bool allow_compression)
  {
    // see if the object is small and copyable via memcpy. if so, use
    // this fast path. otherwise, we have to go through the BOOST
    // serialization machinery
    if (std::is_trivially_copyable<T>() && sizeof(T) * N < 256)
      {
        Assert(std::distance(cbegin, cend) == sizeof(T) * N,
               ExcInternalError());
        std::memcpy(unpacked_object, &*cbegin, sizeof(T) * N);
      }
    else
      {
        // decompress the buffer section into the object
        boost::iostreams::filtering_istreambuf fisb;
#ifdef DEAL_II_WITH_ZLIB
        if (allow_compression)
          fisb.push(boost::iostreams::gzip_decompressor());
#else
        (void)allow_compression;
#endif
        fisb.push(boost::iostreams::array_source(&*cbegin, &*cend));
 
        boost::archive::binary_iarchive bia(fisb);
        bia >> unpacked_object;
      }
  }
 
 
  template <typename T, int N>
  void
  unpack(const std::vector<char> &buffer,
         T (&unpacked_object)[N],
         const bool allow_compression)
  {
    unpack<T, N>(buffer.cbegin(),
                 buffer.cend(),
                 unpacked_object,
                 allow_compression);
  }
 
 
 
  inline bool
  get_bit(const unsigned char number, const unsigned int n)
  {
    AssertIndexRange(n, 8);
 
    // source:
    // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
    // "Checking a bit"
    return (number >> n) & 1U;
  }
 
 
 
  inline void
  set_bit(unsigned char &number, const unsigned int n, const bool x)
  {
    AssertIndexRange(n, 8);
 
    // source:
    // https://stackoverflow.com/questions/47981/how-do-you-set-clear-and-toggle-a-single-bit
    // "Changing the nth bit to x"
    number ^= (-static_cast<unsigned char>(x) ^ number) & (1U << n);
  }
 
 
 
  template <typename To, typename From>
  inline std::unique_ptr<To>
  dynamic_unique_cast(std::unique_ptr<From> &&p)
  {
    // Let's see if we can cast from 'From' to 'To'. If so, do the cast,
    // and then release the pointer from the old
    // owner
    if (To *cast = dynamic_cast<To *>(p.get()))
      {
        std::unique_ptr<To> result(cast);
        p.release();
        return result;
      }
    else
      throw std::bad_cast();
  }
 
 
 
  template <typename T>
  inline T &
  get_underlying_value(T &p)
  {
    return p;
  }
 
 
 
  template <typename T>
  inline T &
  get_underlying_value(std::shared_ptr<T> &p)
  {
    return *p;
  }
 
 
 
  template <typename T>
  inline T &
  get_underlying_value(const std::shared_ptr<T> &p)
  {
    return *p;
  }
 
 
 
  template <typename T>
  inline T &
  get_underlying_value(std::unique_ptr<T> &p)
  {
    return *p;
  }
 
 
 
  template <typename T>
  inline T &
  get_underlying_value(const std::unique_ptr<T> &p)
  {
    return *p;
  }
 
 
 
  template <typename Integer>
  std::vector<Integer>
  reverse_permutation(const std::vector<Integer> &permutation)
  {
    const std::size_t n = permutation.size();
 
    std::vector<Integer> out(n);
    for (std::size_t i = 0; i < n; ++i)
      out[i] = n - 1 - permutation[i];
 
    return out;
  }
 
 
 
  template <typename Integer>
  std::vector<Integer>
  invert_permutation(const std::vector<Integer> &permutation)
  {
    const std::size_t n = permutation.size();
 
    std::vector<Integer> out(n, numbers::invalid_unsigned_int);
 
    for (std::size_t i = 0; i < n; ++i)
      {
        AssertIndexRange(permutation[i], n);
        out[permutation[i]] = i;
      }
 
    // check that we have actually reached
    // all indices
    for (std::size_t i = 0; i < n; ++i)
      Assert(out[i] != numbers::invalid_unsigned_int,
             ExcMessage("The given input permutation had duplicate entries!"));
 
    return out;
  }
} // namespace Utilities
 
 
DEAL_II_NAMESPACE_CLOSE
 
#ifndef DOXYGEN
namespace boost
{
  namespace serialization
  {
    // Provides boost and c++11 with a way to serialize tuples and pairs
    // automatically.
    template <int N>
    struct Serialize
    {
      template <class Archive, typename... Args>
      static void
      serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version)
      {
        ar &std::get<N - 1>(t);
        Serialize<N - 1>::serialize(ar, t, version);
      }
    };
 
    template <>
    struct Serialize<0>
    {
      template <class Archive, typename... Args>
      static void
      serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version)
      {
        (void)ar;
        (void)t;
        (void)version;
      }
    };
 
    template <class Archive, typename... Args>
    void
    serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version)
    {
      Serialize<sizeof...(Args)>::serialize(ar, t, version);
    }
  } // namespace serialization
} // namespace boost
#endif
 
#endif