src/gromacs/math/vectypes.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifndef GMX_MATH_VECTYPES_H
  38 #define GMX_MATH_VECTYPES_H
  39
  40 #include <cmath>
  41
  42 #include <algorithm>
  43 #include <type_traits>
  44
  45 #include "gromacs/utility/gmxassert.h"
  46 #include "gromacs/utility/real.h"
  47
  48 #define XX      0 /* Defines for indexing in */
  49 #define YY      1 /* vectors                 */
  50 #define ZZ      2
  51 #define DIM     3 /* Dimension of vectors    */
  52
  53 typedef real    rvec[DIM];
  54
  55 typedef double  dvec[DIM];
  56
  57 typedef real    matrix[DIM][DIM];
  58
  59 typedef real    tensor[DIM][DIM];
  60
  61 typedef int     ivec[DIM];
  62
  63 namespace gmx
  64 {
  65
  66 /*! \brief
  67  * C++ class for 3D vectors.
  68  *
  69  * \tparam ValueType  Type
  70  *
  71  * This class provides a C++ version of rvec/dvec/ivec that can be put into STL
  72  * containers etc.  It is more or less a drop-in replacement for `rvec` and
  73  * friends: it can be used in most contexts that accept the equivalent C type.
  74  * However, there is one case where explicit conversion is necessary:
  75  *  - An array of these objects needs to be converted with as_vec_array() (or
  76  *    convenience methods like as_rvec_array()).
  77  *
  78  * For the array conversion to work, the compiler should not add any extra
  79  * alignment/padding in the layout of this class;  that this actually works as
  80  * intended is tested in the unit tests.
  81  *
  82  * \inpublicapi
  83  */
  84 template <typename ValueType>
  85 class BasicVector
  86 {
  87     public:
  88         //! Underlying raw C array type (rvec/dvec/ivec).
  89         using RawArray = ValueType[DIM];
  90
  91         // The code here assumes ValueType has been deduced as a data type like int
  92         // and not a pointer like int*. If there is a use case for a 3-element array
  93         // of pointers, the implementation will be different enough that the whole
  94         // template class should have a separate partial specialization. We try to avoid
  95         // accidental matching to pointers, but this assertion is a no-cost extra check.
  96         static_assert(!std::is_pointer<typename std::remove_cv<ValueType>::type>::value,
  97                       "BasicVector value type must not be a pointer.");
  98
  99         //! Constructs default (uninitialized) vector.
 100         BasicVector() {}
 101         //! Constructs a vector from given values.
 102         BasicVector(ValueType x, ValueType y, ValueType z) : x_ {x, y, z}
 103         {}
 104         /*! \brief
 105          * Constructs a vector from given values.
 106          *
 107          * This constructor is not explicit to support implicit conversions
 108          * that allow, e.g., calling `std::vector<RVec>:``:push_back()` directly
 109          * with an `rvec` parameter.
 110          */
 111         BasicVector(const RawArray x) : x_ {x[XX], x[YY], x[ZZ]}
 112         {}
 113         //! Default copy constructor.
 114         BasicVector(const BasicVector &src) = default;
 115         //! Default copy assignment operator.
 116         BasicVector &operator=(const BasicVector &v) = default;
 117         //! Default move constructor.
 118         BasicVector(BasicVector &&src) noexcept = default;
 119         //! Default move assignment operator.
 120         BasicVector &operator=(BasicVector &&v) noexcept = default;
 121         //! Indexing operator to make the class work as the raw array.
 122         ValueType &operator[](int i) { return x_[i]; }
 123         //! Indexing operator to make the class work as the raw array.
 124         ValueType operator[](int i) const { return x_[i]; }
 125         //! Allow inplace addition for BasicVector
 126         BasicVector<ValueType> &operator+=(const BasicVector<ValueType> &right)
 127         {
 128             return *this = *this + right;
 129         }
 130         //! Allow inplace subtraction for BasicVector
 131         BasicVector<ValueType> &operator-=(const BasicVector<ValueType> &right)
 132         {
 133             return *this = *this - right;
 134         }
 135         //! Allow vector addition
 136         BasicVector<ValueType> operator+(const BasicVector<ValueType> &right) const
 137         {
 138             return {x_[0] + right[0], x_[1] + right[1], x_[2] + right[2]};
 139         }
 140         //! Allow vector subtraction
 141         BasicVector<ValueType> operator-(const BasicVector<ValueType> &right) const
 142         {
 143             return {x_[0] - right[0], x_[1] - right[1], x_[2] - right[2]};
 144         }
 145         //! Allow vector scalar division
 146         BasicVector<ValueType> operator/(const ValueType &right) const
 147         {
 148             GMX_ASSERT(right != 0, "Cannot divide by zero");
 149
 150             return *this*(1/right);
 151         }
 152         //! Scale vector by a scalar
 153         BasicVector<ValueType> &operator*=(const ValueType &right)
 154         {
 155             x_[0] *= right;
 156             x_[1] *= right;
 157             x_[2] *= right;
 158
 159             return *this;
 160         }
 161         //! Divide vector by a scalar
 162         BasicVector<ValueType> &operator/=(const ValueType &right)
 163         {
 164             GMX_ASSERT(right != 0, "Cannot divide by zero");
 165
 166             return *this *= 1/right;
 167         }
 168         //! Return dot product
 169         ValueType dot(const BasicVector<ValueType> &right) const
 170         {
 171             return x_[0]*right[0] + x_[1]*right[1] + x_[2]*right[2];
 172         }
 173
 174         //! Allow vector vector multiplication (cross product)
 175         BasicVector<ValueType> cross(const BasicVector<ValueType> &right) const
 176         {
 177             return {
 178                        x_[YY]*right.x_[ZZ]-x_[ZZ]*right.x_[YY],
 179                        x_[ZZ]*right.x_[XX]-x_[XX]*right.x_[ZZ],
 180                        x_[XX]*right.x_[YY]-x_[YY]*right.x_[XX]
 181             };
 182         }
 183
 184         //! Return normalized to unit vector
 185         BasicVector<ValueType> unitVector() const
 186         {
 187             const ValueType vectorNorm = norm();
 188             GMX_ASSERT(vectorNorm != 0, "unitVector() should not be called with a zero vector");
 189
 190             return *this/vectorNorm;
 191         }
 192
 193         //! Length^2 of vector
 194         ValueType norm2() const
 195         {
 196             return dot(*this);
 197         }
 198
 199         //! Norm or length of vector
 200         ValueType norm() const
 201         {
 202             return std::sqrt(norm2());
 203         }
 204
 205         //! cast to RVec
 206         BasicVector<real> toRVec() const
 207         {
 208             return {real(x_[0]), real(x_[1]), real(x_[2])};
 209         }
 210
 211         //! cast to IVec
 212         BasicVector<int> toIVec() const
 213         {
 214             return { static_cast<int>(x_[0]), static_cast<int>(x_[1]), static_cast<int>(x_[2]) };
 215         }
 216
 217         //! cast to DVec
 218         BasicVector<double> toDVec() const
 219         {
 220             return {double(x_[0]), double(x_[1]), double(x_[2])};
 221         }
 222
 223         //! Converts to a raw C array where implicit conversion does not work.
 224         RawArray &as_vec() { return x_; }
 225         //! Converts to a raw C array where implicit conversion does not work.
 226         const RawArray &as_vec() const { return x_; }
 227         //! Makes BasicVector usable in contexts where a raw C array is expected.
 228         operator RawArray &() { return x_; }
 229         //! Makes BasicVector usable in contexts where a raw C array is expected.
 230         operator const RawArray &() const { return x_; }
 231     private:
 232         RawArray x_;
 233 };
 234
 235 //! Allow vector scalar multiplication
 236 template<typename ValueType>
 237 BasicVector<ValueType> operator*(const BasicVector<ValueType> &basicVector,
 238                                  const ValueType              &scalar)
 239 {
 240     return {
 241                basicVector[0]*scalar, basicVector[1]*scalar, basicVector[2]*scalar
 242     };
 243 }
 244
 245 //! Allow scalar vector multiplication
 246 template<typename ValueType>
 247 BasicVector<ValueType> operator*(const ValueType              &scalar,
 248                                  const BasicVector<ValueType> &basicVector)
 249 {
 250     return {
 251                scalar*basicVector[0], scalar*basicVector[1], scalar*basicVector[2]
 252     };
 253 }
 254
 255 /*! \brief
 256  * unitv for gmx::BasicVector
 257  */
 258 template <typename VectorType> static inline
 259 VectorType unitVector(const VectorType &v)
 260 {
 261     return v.unitVector();
 262 }
 263
 264 /*! \brief
 265  * norm for gmx::BasicVector
 266  */
 267 template <typename ValueType> static inline
 268 ValueType norm(BasicVector<ValueType> v)
 269 {
 270     return v.norm();
 271 }
 272
 273 /*! \brief
 274  * Square of the vector norm for gmx::BasicVector
 275  */
 276 template <typename ValueType> static inline
 277 ValueType norm2(BasicVector<ValueType> v)
 278 {
 279     return v.norm2();
 280 }
 281
 282 /*! \brief
 283  * cross product for gmx::BasicVector
 284  */
 285 template <typename VectorType> static inline
 286 VectorType cross(const VectorType &a, const VectorType &b)
 287 {
 288     return a.cross(b);
 289 }
 290
 291 /*! \brief
 292  * dot product for gmx::BasicVector
 293  */
 294 template <typename ValueType> static inline
 295 ValueType dot(BasicVector<ValueType> a, BasicVector<ValueType> b)
 296 {
 297     return a.dot(b);
 298 }
 299
 300 /*! \brief
 301  * Multiply two vectors element by element and return the result.
 302  */
 303 template <typename VectorType>
 304 static inline VectorType scaleByVector(const VectorType &a, const VectorType &b)
 305 {
 306     return {a[0] * b[0], a[1] * b[1], a[2] * b[2]};
 307 }
 308
 309 /*! \brief
 310  * Return the element-wise minimum of two vectors.
 311  */
 312 template <typename VectorType>
 313 static inline VectorType elementWiseMin(const VectorType &a, const VectorType &b)
 314 {
 315     return {std::min(a[0], b[0]), std::min(a[1], b[1]), std::min(a[2], b[2])};
 316 }
 317
 318 /*! \brief
 319  * Return the element-wise maximum of two vectors.
 320  */
 321 template <typename VectorType>
 322 static inline VectorType elementWiseMax(const VectorType &a, const VectorType &b)
 323 {
 324     return {std::max(a[0], b[0]), std::max(a[1], b[1]), std::max(a[2], b[2])};
 325 }
 326
 327 /*! \brief
 328  * Casts a gmx::BasicVector array into an equivalent raw C array.
 329  */
 330 template <typename ValueType> static inline
 331 typename BasicVector<ValueType>::RawArray *
 332 as_vec_array(BasicVector<ValueType> *x)
 333 {
 334     return reinterpret_cast<typename BasicVector<ValueType>::RawArray *>(x);
 335 }
 336
 337 /*! \brief
 338  * Casts a gmx::BasicVector array into an equivalent raw C array.
 339  */
 340 template <typename ValueType> static inline
 341 const typename BasicVector<ValueType>::RawArray *
 342 as_vec_array(const BasicVector<ValueType> *x)
 343 {
 344     return reinterpret_cast<const typename BasicVector<ValueType>::RawArray *>(x);
 345 }
 346
 347 //! Shorthand for C++ `rvec`-equivalent type.
 348 typedef BasicVector<real> RVec;
 349 //! Shorthand for C++ `dvec`-equivalent type.
 350 typedef BasicVector<double> DVec;
 351 //! Shorthand for C++ `ivec`-equivalent type.
 352 typedef BasicVector<int> IVec;
 353 //! Casts a gmx::RVec array into an `rvec` array.
 354 static inline rvec *as_rvec_array(RVec *x)
 355 {
 356     return as_vec_array(x);
 357 }
 358 //! Casts a gmx::RVec array into an `rvec` array.
 359 static inline const rvec *as_rvec_array(const RVec *x)
 360 {
 361     return as_vec_array(x);
 362 }
 363 //! Casts a gmx::DVec array into an `Dvec` array.
 364 static inline dvec *as_dvec_array(DVec *x)
 365 {
 366     return as_vec_array(x);
 367 }
 368 //! Casts a gmx::IVec array into an `ivec` array.
 369 static inline ivec *as_ivec_array(IVec *x)
 370 {
 371     return as_vec_array(x);
 372 }
 373
 374
 375 //! Casts a gmx::DVec array into an `dvec` array.
 376 static inline const dvec *as_dvec_array(const DVec *x)
 377 {
 378     return as_vec_array(x);
 379 }
 380 //! Casts a gmx::IVec array into an `ivec` array.
 381 static inline const ivec *as_ivec_array(const IVec *x)
 382 {
 383     return as_vec_array(x);
 384 }
 385
 386 //! Shorthand for C++ `ivec`-equivalent type.
 387 typedef BasicVector<int> IVec;
 388
 389 }      // namespace gmx
 390
 391 #endif // include guard