point.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2020 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_point_h
#define dealii_point_h


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

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

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/geometries/point.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

#include <cmath>

DEAL_II_NAMESPACE_OPEN

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

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

  explicit DEAL_II_CUDA_HOST_DEV
  Point(const Number x);

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

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

  template <std::size_t dummy_dim,
            typename std::enable_if<(dim == dummy_dim) && (dummy_dim != 0),
                                    int>::type = 0>
  Point(const boost::geometry::model::
          point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt);

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

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

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

  template <typename OtherNumber>
  Point<dim, Number> &
  operator=(const Tensor<1, dim, OtherNumber> &p);

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

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

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

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

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

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

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

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

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

  DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV
Point<dim, Number>::Point() // NOLINT
{}



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



template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV
Point<dim, Number>::Point(const Number x)
{
#  ifndef __CUDA_ARCH__
  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."));
#  endif

  // 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 DEAL_II_CUDA_HOST_DEV
Point<dim, Number>::Point(const Number x, const Number y)
{
#  ifndef __CUDA_ARCH__
  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."));
#  endif

  // 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 DEAL_II_CUDA_HOST_DEV
Point<dim, Number>::Point(const Number x, const Number y, const Number z)
{
#  ifndef __CUDA_ARCH__
  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."));
#  endif

  // 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>
template <
  std::size_t dummy_dim,
  typename std::enable_if<(dim == dummy_dim) && (dummy_dim != 0), int>::type>
inline Point<dim, Number>::Point(
  const boost::geometry::model::
    point<Number, dummy_dim, boost::geometry::cs::cartesian> &boost_pt)
{
  Assert(dim <= 3, ExcNotImplemented());
  this->values[0]                = boost::geometry::get<0>(boost_pt);
  constexpr unsigned int y_index = (dim < 2) ? 0 : 1;
  constexpr unsigned int z_index = (dim < 3) ? 0 : 2;

  if (dim >= 2)
    this->values[y_index] = boost::geometry::get<y_index>(boost_pt);

  if (dim >= 3)
    this->values[z_index] = boost::geometry::get<z_index>(boost_pt);
}



template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV Number
Point<dim, Number>::operator()(const unsigned int index) const
{
#  ifndef __CUDA_ARCH__
  AssertIndexRange(index, dim);
#  endif
  return this->values[index];
}



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



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



template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV
    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 DEAL_II_CUDA_HOST_DEV
  Point<dim,
        typename ProductType<Number,
                             typename EnableIfScalar<OtherNumber>::type>::type>
  Point<dim, Number>::operator/(const OtherNumber factor) const
{
  const Tensor<1, dim, Number> &base_object = *this;
  return Point<
    dim,
    typename ProductType<Number,
                         typename EnableIfScalar<OtherNumber>::type>::type>(
    ::operator/(base_object, factor));
}



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



template <int dim, typename Number>
inline DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV 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 DEAL_II_CUDA_HOST_DEV
  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