point.h

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

#ifndef dealii_point_h
#define dealii_point_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/tensor.h>
#include <cmath>

DEAL_II_NAMESPACE_OPEN

template <int dim, typename Number = double>
class Point : public Tensor<1,dim,Number>
{
public:
  Point ();

  explicit Point (const Tensor<1,dim,Number> &);

  explicit Point (const Number x);

  Point (const Number x,
         const Number y);

  Point (const Number x,
         const Number y,
         const Number z);

  static Point<dim,Number> unit_vector(const unsigned int i);

  Number operator () (const unsigned int index) const;

  Number &operator () (const unsigned int index);

  /*
   * @name Addition and subtraction of points.
   * @{
   */

  Point<dim,Number> operator + (const Tensor<1,dim,Number> &) const;

  Tensor<1,dim,Number> operator - (const Point<dim,Number> &) const;

  Point<dim,Number> operator - (const Tensor<1,dim,Number> &) const;

  Point<dim,Number> operator - () const;

  /*
   * @name Multiplication and scaling of points. Dot products. Norms.
   * @{
   */

  template <typename OtherNumber>
  Point<dim,typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::type>
  operator * (const OtherNumber) const;

  template <typename OtherNumber>
  Point<dim,typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::type>
  operator / (const OtherNumber) const;

  Number operator * (const Tensor<1,dim,Number> &p) const;

  typename numbers::NumberTraits<Number>::real_type square () const;

  typename numbers::NumberTraits<Number>::real_type distance (const Point<dim,Number> &p) const;

  typename numbers::NumberTraits<Number>::real_type distance_square (const Point<dim,Number> &p) const;

  template <class Archive>
  void serialize(Archive &ar, const unsigned int version);
};

/*------------------------------- Inline functions: Point ---------------------------*/

#ifndef DOXYGEN

// At least clang-3.7 requires us to have a user-defined constructor
// and we can't use 'Point<dim,Number>::Point () = default' here.
template <int dim, typename Number>
inline
Point<dim,Number>::Point ()
{}



template <int dim, typename Number>
inline
Point<dim,Number>::Point (const Tensor<1,dim,Number> &t)
  :
  Tensor<1,dim,Number>(t)
{}



template <int dim, typename Number>
inline
Point<dim,Number>::Point (const Number x)
{
  Assert (dim==1,
          ExcMessage ("You can only initialize Point<1> objects using the constructor "
                      "that takes only one argument. Point<dim> objects with dim!=1 "
                      "require initialization with the constructor that takes 'dim' "
                      "arguments."));

  // we can only get here if we pass the assertion. use the switch anyway so
  // as to avoid compiler warnings about uninitialized elements or writing
  // beyond the end of the 'values' array
  switch (dim)
    {
    case 1:
      this->values[0] = x;
      break;

    default:
      ;
    }
}



template <int dim, typename Number>
inline
Point<dim,Number>::Point (const Number x,
                          const Number y)
{
  Assert (dim==2,
          ExcMessage ("You can only initialize Point<2> objects using the constructor "
                      "that takes two arguments. Point<dim> objects with dim!=2 "
                      "require initialization with the constructor that takes 'dim' "
                      "arguments."));
  // we can only get here if we pass the assertion. use the indirection anyway
  // so as to avoid compiler warnings about uninitialized elements or writing
  // beyond the end of the 'values' array
  constexpr unsigned int y_index = (dim<2)?0:1;
  this->values[0] = x;
  this->values[y_index] = y;
}



template <int dim, typename Number>
inline
Point<dim,Number>::Point (const Number x,
                          const Number y,
                          const Number z)
{
  Assert (dim==3,
          ExcMessage ("You can only initialize Point<3> objects using the constructor "
                      "that takes three arguments. Point<dim> objects with dim!=3 "
                      "require initialization with the constructor that takes 'dim' "
                      "arguments."));

  // we can only get here if we pass the assertion. use the indirection anyway
  // so as to avoid compiler warnings about uninitialized elements or writing
  // beyond the end of the 'values' array
  constexpr unsigned int y_index = (dim<2)?0:1;
  constexpr unsigned int z_index = (dim<3)?0:2;
  this->values[0] = x;
  this->values[y_index] = y;
  this->values[z_index] = z;
}


template <int dim, typename Number>
inline
Point<dim,Number>
Point<dim,Number>::unit_vector(unsigned int i)
{
  Point<dim,Number> p;
  p[i] = 1.;
  return p;
}


template <int dim, typename Number>
inline
Number
Point<dim,Number>::operator () (const unsigned int index) const
{
  AssertIndexRange((int) index, dim);
  return this->values[index];
}



template <int dim, typename Number>
inline
Number &
Point<dim,Number>::operator () (const unsigned int index)
{
  AssertIndexRange((int) index, dim);
  return this->values[index];
}



template <int dim, typename Number>
inline
Point<dim,Number>
Point<dim,Number>::operator + (const Tensor<1,dim,Number> &p) const
{
  Point<dim,Number> tmp = *this;
  tmp += p;
  return tmp;
}



template <int dim, typename Number>
inline
Tensor<1,dim,Number>
Point<dim,Number>::operator - (const Point<dim,Number> &p) const
{
  return (Tensor<1,dim,Number>(*this) -= p);
}



template <int dim, typename Number>
inline
Point<dim,Number>
Point<dim,Number>::operator - (const Tensor<1,dim,Number> &p) const
{
  Point<dim,Number> tmp = *this;
  tmp -= p;
  return tmp;
}



template <int dim, typename Number>
inline
Point<dim,Number>
Point<dim,Number>::operator - () const
{
  Point<dim,Number> result;
  for (unsigned int i=0; i<dim; ++i)
    result.values[i] = -this->values[i];
  return result;
}



template <int dim, typename Number>
template <typename OtherNumber>
inline
Point<dim,typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::type>
Point<dim,Number>::operator * (const OtherNumber factor) const
{
  Point<dim,typename ProductType<Number, OtherNumber>::type> tmp;
  for (unsigned int i=0; i<dim; ++i)
    tmp[i] = this->operator[](i) * factor;
  return tmp;
}



template <int dim, typename Number>
template <typename OtherNumber>
inline
Point<dim,typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::type>
Point<dim,Number>::operator / (const OtherNumber factor) const
{
  Point<dim,typename ProductType<Number, OtherNumber>::type> tmp;
  for (unsigned int i=0; i<dim; ++i)
    tmp[i] = this->operator[](i) / factor;
  return tmp;
}



template <int dim, typename Number>
inline
Number
Point<dim,Number>::operator * (const Tensor<1,dim,Number> &p) const
{
  Number res = Number();
  for (unsigned int i=0; i<dim; ++i)
    res += this->operator[](i) * p[i];
  return res;
}


template <int dim, typename Number>
inline
typename numbers::NumberTraits<Number>::real_type
Point<dim,Number>::square () const
{
  return this->norm_square();
}



template <int dim, typename Number>
inline
typename numbers::NumberTraits<Number>::real_type
Point<dim,Number>::distance (const Point<dim,Number> &p) const
{
  return std::sqrt(distance_square(p));
}



template <int dim, typename Number>
inline
typename numbers::NumberTraits<Number>::real_type
Point<dim,Number>::distance_square (const Point<dim,Number> &p) const
{
  Number sum = internal::NumberType<Number>::value(0.0);
  for (unsigned int i=0; i<dim; ++i)
    {
      const Number diff=static_cast<Number>(this->values[i])-p(i);
      sum += numbers::NumberTraits<Number>::abs_square (diff);
    }

  return sum;
}



template <int dim, typename Number>
template <class Archive>
inline
void
Point<dim,Number>::serialize(Archive &ar, const unsigned int)
{
  // forward to serialization
  // function in the base class
  ar   &static_cast<Tensor<1,dim,Number> &>(*this);
}

#endif // DOXYGEN


/*------------------------------- Global functions: Point ---------------------------*/


template <int dim, typename Number, typename OtherNumber>
inline
Point<dim,typename ProductType<Number, typename EnableIfScalar<OtherNumber>::type>::type>
operator * (const OtherNumber        factor,
            const Point<dim,Number> &p)
{
  return p * factor;
}



template <int dim, typename Number>
inline
std::ostream &operator << (std::ostream            &out,
                           const Point<dim,Number> &p)
{
  for (unsigned int i=0; i<dim-1; ++i)
    out << p[i] << ' ';
  out << p[dim-1];

  return out;
}



template <int dim, typename Number>
inline
std::istream &operator >> (std::istream      &in,
                           Point<dim,Number> &p)
{
  for (unsigned int i=0; i<dim; ++i)
    in >> p[i];

  return in;
}


#ifndef DOXYGEN

template <typename Number>
inline
std::ostream &operator << (std::ostream &out,
                           const Point<1,Number> &p)
{
  out << p[0];

  return out;
}

#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE

#endif