signaling_nan.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_signaling_nan_h
#define dealii_signaling_nan_h

#include <deal.II/base/config.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/point.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;
        }
      };
    }
  }




  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 ();
  }
}


DEAL_II_NAMESPACE_CLOSE

#endif