Point.h

// Copyright (C) 2006-2014 Anders Logg
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can 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 3 of the License, or
// (at your option) any later version.
//
// DOLFIN is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// Modified by Garth N. Wells 2006
// Modified by Andre Massing 2009
//
// First added:  2006-06-12
// Last changed: 2017-09-28

#ifndef __POINT_H
#define __POINT_H

#include <array>
#include <cmath>
#include <iostream>
#include <dolfin/log/log.h>
#include <dolfin/common/Array.h>

namespace dolfin
{

  class Point
  {
  public:

    explicit Point(const double x=0.0, const double y=0.0, const double z=0.0)
      : _x({{x, y, z}}) {}

    Point(std::size_t dim, const double* x) : _x({{0.0, 0.0, 0.0}})
    {
      for (std::size_t i = 0; i < dim; i++)
        _x[i] = x[i];
    }

    Point(const Array<double>& x) : _x({{0.0, 0.0, 0.0}})
    {
      for (std::size_t i = 0; i < x.size(); i++)
        _x[i] = x[i];
    }

    Point(const Point& p) : _x({{p._x[0], p._x[1], p._x[2]}}) {}

    ~Point() {}

    double& operator[] (std::size_t i)
    { dolfin_assert(i < 3); return _x[i]; }

    double operator[] (std::size_t i) const
    { dolfin_assert(i < 3); return _x[i]; }

    double x() const
    { return _x[0]; }

    double y() const
    { return _x[1]; }

    double z() const
    { return _x[2]; }

    double* coordinates()
    { return _x.data(); }

    const double* coordinates() const
    { return _x.data(); }

    std::array<double, 3> array() const
    {
      return _x;
    }

    Point operator+ (const Point& p) const
    { Point q(_x[0] + p._x[0], _x[1] + p._x[1], _x[2] + p._x[2]); return q; }

    Point operator- (const Point& p) const
    { Point q(_x[0] - p._x[0], _x[1] - p._x[1], _x[2] - p._x[2]); return q; }

    const Point& operator+= (const Point& p)
    { _x[0] += p._x[0]; _x[1] += p._x[1]; _x[2] += p._x[2]; return *this; }

    const Point& operator-= (const Point& p)
    { _x[0] -= p._x[0]; _x[1] -= p._x[1]; _x[2] -= p._x[2]; return *this; }

    Point operator- ()
    { Point p(-_x[0], -_x[1], -_x[2]); return p; }

    Point operator* (double a) const
    { Point p(a*_x[0], a*_x[1], a*_x[2]); return p; }

    const Point& operator*= (double a)
    { _x[0] *= a; _x[1] *= a; _x[2] *= a; return *this; }

    Point operator/ (double a) const
    { Point p(_x[0]/a, _x[1]/a, _x[2]/a); return p; }

    const Point& operator/= (double a)
    { _x[0] /= a; _x[1] /= a; _x[2] /= a; return *this; }

    const Point& operator= (const Point& p)
    { _x = {{p._x[0], p._x[1], p._x[2]}}; return *this; }

    bool operator== (const Point& p) const
    { return _x == p._x; }

    bool operator!= (const Point& p) const
    { return _x != p._x; }

    double squared_distance(const Point& p) const;

    inline double distance(const Point& p) const
    { return sqrt(squared_distance(p)); }

    double norm() const
    { return std::sqrt(_x[0]*_x[0] + _x[1]*_x[1] + _x[2]*_x[2]); }

    double squared_norm() const
    { return _x[0]*_x[0] + _x[1]*_x[1] + _x[2]*_x[2]; }

    const Point cross(const Point& p) const;

    double dot(const Point& p) const;

    Point rotate(const Point& a, double theta) const;

    // Note: Not a subclass of Variable for efficiency!

    std::string str(bool verbose=false) const;

  private:

    std::array<double, 3> _x;

  };

  inline Point operator*(double a, const Point& p)
  { return p*a; }

  inline std::ostream& operator<<(std::ostream& stream, const Point& point)
  { stream << point.str(false); return stream; }

}

#endif