Initialize at openlb-1-3

1 files changed, 482 insertions, 0 deletions

diff --git a/src/core/vector.h b/src/core/vector.h
new file mode 100644
index 0000000..d9d8605
--- /dev/null
+++ b/src/core/vector.h
@@ -0,0 +1,482 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2015 Asher Zarth, Mathias J. Krause, Albert Mink
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program 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 General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+*/
+
+
+/** \file
+ * efficient implementation of a vector class
+ */
+
+#ifndef VECTOR_H
+#define VECTOR_H
+
+#include <iostream>
+#include <cstring>
+#include <type_traits>
+#include "util.h"
+#include "math.h"
+#include "olbDebug.h"
+
+typedef double S;
+
+// Use BaseType double as the default BaseType.
+#ifndef BaseType
+///TODO
+//check OLD_PRECONDITION
+//change this typedef, double for testing only
+//#define BaseType S
+#define BaseType S
+#endif
+
+namespace olb {
+
+/** Efficient implementation of a vector class.
+ * Static size.
+ * Default value of vector element is T(0)
+ */
+template <typename T, unsigned Size>
+class Vector {
+public:
+  /// number of dimensions
+  T data[Size];
+public:
+  /// Constructor from array
+  Vector(const T ar[Size]);
+  /// Constructor from std::vector<T>
+  Vector(const std::vector<T>& v);
+  /// Construct with all entries as explicit scalar
+  explicit Vector(const T& s = T());
+  /// Construct 2D, 2 explicit arguments
+  Vector(const T& a, const T& b);
+  /// Construct 3D, 3 explicit arguments
+  Vector(const T& a, const T& b, const T& c);
+  /// copy constructor
+  Vector(const Vector<T,Size>& v);
+  /// copy assignment operator
+  Vector& operator = (const Vector<T,Size>& v);
+  /// destructor
+  ~Vector() = default;
+  /// element access
+  T& operator[](unsigned n);
+  /// element access (read only)
+  const T& operator[](const unsigned n) const;
+  /// cumulative add vector
+  Vector& operator += (const Vector<T,Size>& v);
+  /// cumulatively subtract vector
+  Vector& operator -= (const Vector<T,Size>& v);
+  /// add scalar to each entry
+  Vector& operator += (const T& s);
+  /// subtract scalar from each entry
+  Vector& operator -= (const T& s);
+  /// multiply by scalar
+  Vector& operator *= (const T& s);
+  /// divide by scalar
+  Vector& operator /= (const T& s);
+  /// equal operator returns true if all components match
+  bool operator == (const Vector<T, Size>& v);
+  /// unequal operator returns true if at least one components differs
+  bool operator != (const Vector<T, Size>& v);
+  /// get l2 norm of vector
+  T norm();
+  /// get normalized vector
+  void normalize(const T& scale = T(1));
+  /// check if vector has zero length
+  bool closeToZero();
+  //operator BaseType();
+  /// convert to std::vector
+  std::vector<T> toStdVector() const;
+};
+
+
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const T ar[Size])
+{
+  //static_assert(std::is_trivially_copyable<T>::value, "Only trivially copyable objects can safely be copied by memcpy");
+  std::memcpy( data, ar, Size*sizeof(T));
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const std::vector<T>& v)
+{
+  OLB_PRECONDITION(v.size() == Size);
+  //static_assert(std::is_trivially_copyable<T>::value, "Only trivially copyable objects can safely be copied by memcpy");
+  std::memcpy( data, v.data(), Size*sizeof(T));
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const T& s)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] = s;
+  }
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const T& a, const T& b)
+{
+  OLB_PRECONDITION(Size == 2);
+  data[0] = a;
+  data[1] = b;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const T& a, const T& b, const T& c)
+{
+  OLB_PRECONDITION(Size == 3);
+  data[0] = a;
+  data[1] = b;
+  data[2] = c;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>::Vector(const Vector<T,Size>& v)
+{
+  //static_assert(std::is_trivially_copyable<T>::value, "Only trivially copyable objects can safely be copied by memcpy");
+  std::memcpy( data, v.data, Size*sizeof(T));
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator= (const Vector<T,Size>& v)
+{
+  //static_assert(std::is_trivially_copyable<T>::value, "Only trivially copyable objects can safely be copied by memcpy");
+  std::memcpy( data, v.data, Size*sizeof(T));
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator+= (const Vector<T,Size>& v)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] += v[i];
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator+= (const T& s)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] += s;
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator-= (const Vector<T,Size>& v)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] -= v[i];
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator*= (const T& s)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] *= s;
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator/= (const T& s)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] /= s;
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline Vector<T, Size>& Vector<T, Size>::operator-= (const T& s)
+{
+  for (unsigned i = 0; i < Size; ++i) {
+    data[i] -= s;
+  }
+  return *this;
+}
+
+template <typename T, unsigned Size>
+inline bool Vector<T, Size>::operator== (const Vector<T, Size>& v)
+{
+  bool theSame = true;
+  for (unsigned i = 0; i < Size; ++i) {
+    theSame &= data[i] == v.data[i];
+  }
+  return theSame;
+}
+
+template <typename T, unsigned Size>
+inline bool Vector<T, Size>::operator!= (const Vector<T, Size>& v)
+{
+  return !(*this == v);
+}
+
+template <typename T, unsigned Size>
+inline T& Vector<T, Size>::operator[](unsigned n)
+
+{
+  OLB_PRECONDITION(n < Size);
+  return data[n];
+}
+
+template <typename T, unsigned Size>
+inline const T& Vector<T, Size>::operator[](const unsigned n) const
+{
+  OLB_PRECONDITION(n < Size);
+  return data[n];
+}
+
+template <typename T, unsigned Size>
+inline bool Vector<T, Size>::closeToZero()
+{
+  bool zero = true;
+  T EPSILON = std::numeric_limits<T>::epsilon();
+
+  for (unsigned i = 0; i < Size; ++i) {
+    if (fabs(data[i]) > EPSILON) {
+      zero = false;
+    }
+  }
+  return zero;
+}
+
+template <typename T, unsigned Size>
+inline T Vector<T, Size>::norm()
+{
+  T v(0);
+  for (unsigned iD = 0; iD < Size; ++iD) {
+    v += data[iD]*data[iD];
+  }
+  v = sqrt(v);
+  return v;
+}
+
+template <typename T, unsigned Size>
+inline void Vector<T, Size>::normalize(const T& scale)
+{
+  T invScale = scale/this->norm();
+  //assert(scale > 0);
+  for (unsigned int iDim = 0; iDim < Size; ++iDim) {
+    data[iDim] *= invScale;
+  }
+}
+
+template <typename T, unsigned Size>
+inline std::vector<T> Vector<T, Size>::toStdVector() const
+{
+  return std::vector<T> (data, data+Size);
+}
+
+template <class T, unsigned Size>
+inline std::ostream& operator << (std::ostream& os, const Vector<T,Size>& o)
+{
+  if (Size != 0) {
+    os << "[";
+    for (unsigned i = 0; i < Size-1; ++i) {
+      os << o[i] << " ";
+    }
+    os << o[Size-1]<<"]";
+  }
+  else {
+    os << "[empty]";
+  }
+  return os;
+}
+
+template <typename T>
+inline Vector<T, 3> crossProduct3D(const Vector<T, 3>& a, const Vector<T,3>& b)
+{
+  Vector<T, 3> v(a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]);
+  return v;
+}
+
+template <typename T, unsigned Size>
+inline T norm(const Vector<T, Size>& a)
+{
+  T v(0);
+  for (unsigned iD = 0; iD < Size; ++iD) {
+    v += a[iD]*a[iD];
+  }
+  v = sqrt(v);
+  return v;
+}
+
+template <typename T, unsigned Size>
+inline void normalize(Vector<T, Size>& a)
+{
+  T scale = norm(a);
+  //assert(scale > 0);
+  for (unsigned iDim = 0; iDim < Size; ++iDim) {
+    a[iDim] /= scale;
+  }
+}
+
+
+
+
+// indirect helper functions, providing general arithmetic
+// overloading of operators for interaction with BaseType and itself
+// addition, subtraction, multiplication and division
+
+// PLUS
+template <class T, unsigned Size>
+inline Vector<T, Size> operator+ (const BaseType& a, const Vector<T, Size>& b)
+{
+  return Vector<T, Size>(b)+=a;
+}
+
+template <class T, unsigned Size>
+inline Vector<T, Size> operator+ (const Vector<T, Size>& a, const BaseType& b)
+{
+  return Vector<T, Size>(a)+=b;
+}
+
+template <class T, unsigned Size>
+inline Vector<T, Size> operator+ (const Vector<T, Size>& a, const Vector<T, Size>& b)
+{
+  return Vector<T, Size>(a)+=b;
+}
+
+template <typename T, typename W, unsigned Size>
+inline Vector<decltype(T{}+W{}), Size> operator+ (const Vector<T, Size>& a, const Vector<W, Size>& b)
+{
+  Vector<decltype(T{}+W{}), Size> result;
+  for (unsigned i = 0; i < Size; ++i) {
+    result[i] = a[i] + b[i];
+  }
+  return result;
+}
+
+// MINUS
+template <class T, unsigned Size>
+inline Vector<T, Size> operator- (const BaseType& a, const Vector<T, Size>& b)
+{
+  return Vector<T, Size>(a)-=b;
+}
+
+template <class T, unsigned Size>
+inline Vector<T, Size> operator- (const Vector<T, Size>& a, const BaseType& b)
+{
+  return Vector<T, Size>(a)-=b;
+}
+
+template <class T, unsigned Size>
+inline Vector<T, Size> operator- (const Vector<T, Size>& a, const Vector<T, Size>& b)
+{
+  return Vector<T, Size>(a)-=b;
+}
+
+template <typename T, typename W, unsigned Size>
+inline Vector<decltype(T{}-W{}), Size> operator- (const Vector<T, Size>& a, const Vector<W, Size>& b)
+{
+  Vector<decltype(T{}-W{}), Size> result;
+  for (unsigned i = 0; i < Size; ++i) {
+    result[i] = a[i] - b[i];
+  }
+  return result;
+}
+
+// MULTIPLY
+template <class T, unsigned Size>
+inline Vector<T, Size> operator* (const BaseType& a, const Vector<T, Size>& b)
+{
+  return Vector<T, Size>(b)*=a;
+}
+
+template <class T, unsigned Size>
+inline Vector<T, Size> operator* (const Vector<T, Size>& a, const BaseType& b)
+{
+  return Vector<T, Size>(a)*=b;
+}
+
+template <class T, unsigned Size>
+inline T operator* (const Vector<T, Size>& a, const Vector<T, Size>& b)
+{
+  T scalarProduct = T();
+  for (unsigned iD = 0; iD < Size; ++iD) {
+    scalarProduct += a[iD]*b[iD];
+  }
+  return scalarProduct;
+}
+
+// DIVIDE
+template <class T, unsigned Size>
+inline Vector<T, Size> operator/ (const Vector<T, Size>& a, const BaseType& b)
+{
+  return Vector<T, Size>(a)/=b;
+}
+
+/// < Operator returns true if all components are smaller
+template <class T, unsigned Size>
+inline bool operator< (const Vector<T, Size>& lhs, const Vector<T, Size>& rhs)
+{
+  bool smaller = true;
+  for (unsigned i = 0; i < Size; ++i) {
+    smaller &= (lhs[i] < rhs[i]);
+  }
+  return smaller;
+}
+
+/// > Operator returns true if all components are greater
+template <class T, unsigned Size>
+inline bool operator> (const Vector<T, Size>& lhs, const Vector<T, Size>& rhs)
+{
+  return rhs < lhs;
+}
+
+/// <= Operator returns true if all components are smaller or equal
+template <class T, unsigned Size>
+inline bool operator<=(const Vector<T, Size>& lhs, const Vector<T, Size>& rhs)
+{
+  bool smallerEq = true;
+  for (unsigned i = 0; i < Size; ++i) {
+    smallerEq &= (lhs[i] <= rhs[i]);
+  }
+  return smallerEq;
+}
+
+/// >= Operator returns true if all components are smaller or equal
+template <class T, unsigned Size>
+inline bool operator>=(const Vector<T, Size>& lhs, const Vector<T, Size>& rhs)
+{
+  return rhs <= lhs;
+}
+
+
+template <typename T>
+using Scalar   = Vector<T,1>;
+
+template <typename T>
+using Vector2D = Vector<T,2>;
+
+template <typename T>
+using Vector3D = Vector<T,3>;
+
+template <typename T>
+using Vector4D = Vector<T,4>;
+
+} // end namespace olb
+
+#endif