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