signaling_nan.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2018 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_signaling_nan_h
#define dealii_signaling_nan_h

#include <deal.II/base/config.h>

#include <deal.II/base/derivative_form.h>
#include <deal.II/base/point.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>

#include <limits>


DEAL_II_NAMESPACE_OPEN

namespace numbers
{
  namespace internal
  {
    namespace SignalingNaN
    {
      template <typename T>
      struct NaNInitializer;


      template <>
      struct NaNInitializer<float>
      {
        static float
        invalid_element()
        {
          return std::numeric_limits<float>::signaling_NaN();
        }
      };


      template <>
      struct NaNInitializer<double>
      {
        static double
        invalid_element()
        {
          return std::numeric_limits<double>::signaling_NaN();
        }
      };


      template <int dim, typename T>
      struct NaNInitializer<Tensor<1, dim, T>>
      {
        static Tensor<1, dim, T>
        invalid_element()
        {
          Tensor<1, dim, T> nan_tensor;

          for (unsigned int i = 0; i < dim; ++i)
            nan_tensor[i] = NaNInitializer<T>::invalid_element();

          return nan_tensor;
        }
      };



      template <int rank, int dim, typename T>
      struct NaNInitializer<Tensor<rank, dim, T>>
      {
        static Tensor<rank, dim, T>
        invalid_element()
        {
          Tensor<rank, dim, T> nan_tensor;

          // recursively initialize sub-tensors with invalid elements
          for (unsigned int i = 0; i < dim; ++i)
            nan_tensor[i] =
              NaNInitializer<Tensor<rank - 1, dim, T>>::invalid_element();

          return nan_tensor;
        }
      };



      template <int dim, typename T>
      struct NaNInitializer<Point<dim, T>>
      {
        static Point<dim, T>
        invalid_element()
        {
          Point<dim, T> nan_point;

          for (unsigned int i = 0; i < dim; ++i)
            nan_point[i] = NaNInitializer<T>::invalid_element();

          return nan_point;
        }
      };



      template <int rank, int dim, typename T>
      struct NaNInitializer<SymmetricTensor<rank, dim, T>>
      {
        static SymmetricTensor<rank, dim, T>
        invalid_element()
        {
          // initialize symmetric tensors via the unrolled list of elements
          T initializers
            [SymmetricTensor<rank, dim, T>::n_independent_components];
          for (unsigned int i = 0;
               i < SymmetricTensor<rank, dim, T>::n_independent_components;
               ++i)
            initializers[i] = NaNInitializer<T>::invalid_element();

          return SymmetricTensor<rank, dim, T>(initializers);
        }
      };



      template <int order, int dim, int spacedim, typename T>
      struct NaNInitializer<DerivativeForm<order, dim, spacedim, T>>
      {
        static DerivativeForm<order, dim, spacedim, T>
        invalid_element()
        {
          DerivativeForm<order, dim, spacedim, T> form;

          // recursively initialize sub-tensors with invalid elements
          for (unsigned int i = 0; i < spacedim; ++i)
            form[i] = NaNInitializer<Tensor<order, dim, T>>::invalid_element();

          return form;
        }
      };
    } // namespace SignalingNaN
  }   // namespace internal



  template <class T>
  T
  signaling_nan()
  {
    // dispatch to the classes in the internal namespace because there
    // we can do partial specializations, which is not possible for
    // template functions such as the current one
    return internal::SignalingNaN::NaNInitializer<T>::invalid_element();
  }
} // namespace numbers


DEAL_II_NAMESPACE_CLOSE

#endif