Initialize at openlb-1-3

1 files changed, 471 insertions, 0 deletions

diff --git a/src/core/util.h b/src/core/util.h
new file mode 100644
index 0000000..7faa31a
--- /dev/null
+++ b/src/core/util.h
@@ -0,0 +1,471 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2006, 2007 Jonas Latt
+ *  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
+ * Set of functions commonly used in LB computations
+ *  -- header file
+ */
+#ifndef UTIL_H
+#define UTIL_H
+
+#include<sstream>
+#include<algorithm>
+#include<utilities/vectorHelpers.h>
+
+#include "dynamics/descriptorFunction.h"
+
+// patch due to ploblems with older compilers
+namespace std {
+template<typename T>
+std::string to_string(const T &n)
+{
+  std::ostringstream s;
+  s << n;
+  return s.str();
+}
+}
+
+namespace olb {
+
+namespace util {
+
+template<typename T> T norm(const std::vector<T>& a);
+
+template <typename T>
+inline int sign(T val)
+{
+  return (0 < val) - (val < 0);
+}
+
+template <typename T>
+inline bool aligned_to_x(const std::vector<T>& vec)
+{
+  return (vec[0]!=0 and vec[1]==0 and vec[2]==0);
+}
+
+template <typename T>
+inline bool aligned_to_y(const std::vector<T>& vec)
+{
+  return (vec[0]==0 and vec[1]!=0 and vec[2]==0);
+}
+
+template <typename T>
+inline bool aligned_to_z(const std::vector<T>& vec)
+{
+  return (vec[0]==0 and vec[1]==0 and vec[2]!=0);
+}
+
+template <typename T>
+inline bool aligned_to_grid(const std::vector<T>& vec)
+{
+  return (aligned_to_x<T>(vec) or
+          aligned_to_y<T>(vec) or
+          aligned_to_z<T>(vec));
+}
+
+
+
+
+
+inline bool intersect (
+  int x0, int x1, int y0, int y1,
+  int x0_, int x1_, int y0_, int y1_,
+  int& newX0, int& newX1, int& newY0, int& newY1 )
+{
+  newX0 = std::max(x0,x0_);
+  newY0 = std::max(y0,y0_);
+
+  newX1 = std::min(x1,x1_);
+  newY1 = std::min(y1,y1_);
+
+  return newX1>=newX0 && newY1>=newY0;
+}
+
+inline bool intersect (
+  int x0, int x1, int y0, int y1, int z0, int z1,
+  int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
+  int& newX0, int& newX1, int& newY0, int& newY1, int& newZ0, int& newZ1 )
+{
+  newX0 = std::max(x0,x0_);
+  newY0 = std::max(y0,y0_);
+  newZ0 = std::max(z0,z0_);
+
+  newX1 = std::min(x1,x1_);
+  newY1 = std::min(y1,y1_);
+  newZ1 = std::min(z1,z1_);
+
+  return newX1>=newX0 && newY1>=newY0 && newZ1>=newZ0;
+}
+
+inline bool contained(int x, int y,
+                      int x0, int x1, int y0, int y1)
+{
+  return x>=x0 && x<=x1 &&
+         y>=y0 && y<=y1;
+}
+
+inline bool contained(int x, int y, int z,
+                      int x0, int x1, int y0, int y1, int z0, int z1)
+{
+  return x>=x0 && x<=x1 &&
+         y>=y0 && y<=y1 &&
+         z>=z0 && z<=z1;
+}
+
+
+template<typename T>
+T sqr(T arg)
+{
+  return arg*arg;
+}
+
+/// Compute norm square of a d-dimensional vector
+template<typename T, unsigned D>
+T normSqr(const T u[D])
+{
+  T uSqr = T();
+  for (unsigned iD=0; iD < D; ++iD) {
+    uSqr += u[iD]*u[iD];
+  }
+  return uSqr;
+}
+
+/// Compute norm square of a d-dimensional vector
+template<typename T, unsigned D>
+T normSqr(const Vector<T,D>& u)
+{
+  T uSqr = T();
+  for (unsigned iD=0; iD < D; ++iD) {
+    uSqr += u[iD]*u[iD];
+  }
+  return uSqr;
+}
+
+template<typename T, int d>
+T scalarProduct(const T u1[d], const T u2[d])
+{
+  T prod = T();
+  for (int iD=0; iD<d; ++iD) {
+    prod += u1[iD]*u2[iD];
+  }
+  return prod;
+}
+
+template<typename T>
+T scalarProduct(const std::vector<T>& u1, const std::vector<T>& u2)
+{
+  T prod = T();
+  if (u1.size() == u2.size()) {
+    for (int iD=0; iD<u1.size(); ++iD) {
+      prod += u1[iD]*u2[iD];
+    }
+  }
+  return prod;
+}
+
+/// Compute number of elements of a symmetric d-dimensional tensor
+template <typename DESCRIPTORBASE> struct TensorVal {
+  static const int n =
+    (DESCRIPTORBASE::d*(DESCRIPTORBASE::d+1))/2; ///< result stored in n
+};
+
+/// Compute the opposite of a given direction
+template <typename DESCRIPTORBASE> inline int opposite(int iPop)
+{
+  return descriptors::opposite<DESCRIPTORBASE>(iPop);
+}
+
+template <typename DESCRIPTORBASE, int index, int value>
+class SubIndex {
+private:
+  SubIndex()
+  {
+    for (int iVel=0; iVel<DESCRIPTORBASE::q; ++iVel) {
+      if (descriptors::c<DESCRIPTORBASE>(iVel,index)==value) {
+        indices.push_back(iVel);
+      }
+    }
+  }
+
+  std::vector<int> indices;
+
+  template <typename DESCRIPTORBASE_, int index_, int value_>
+  friend std::vector<int> const& subIndex();
+};
+
+template <typename DESCRIPTORBASE, int index, int value>
+std::vector<int> const& subIndex()
+{
+  static SubIndex<DESCRIPTORBASE, index, value> subIndexSingleton;
+  return subIndexSingleton.indices;
+}
+
+template <typename DESCRIPTORBASE>
+int findVelocity(const int v[DESCRIPTORBASE::d])
+{
+  for (int iPop=0; iPop<DESCRIPTORBASE::q; ++iPop) {
+    bool fit = true;
+    for (int iD=0; iD<DESCRIPTORBASE::d; ++iD) {
+      if (descriptors::c<DESCRIPTORBASE>(iPop,iD) != v[iD]) {
+        fit = false;
+        break;
+      }
+    }
+    if (fit) {
+      return iPop;
+    }
+  }
+  return DESCRIPTORBASE::q;
+}
+
+/**
+* finds distributions incoming into the wall
+* but we want the ones outgoing from the wall,
+* therefore we have to take the opposite ones.
+*/
+template <typename DESCRIPTORBASE, int direction, int orientation>
+class SubIndexOutgoing {
+private:
+  SubIndexOutgoing()   // finds the indexes outgoing from the walls
+  {
+    indices = util::subIndex<DESCRIPTORBASE,direction,orientation>();
+
+    for (unsigned iPop = 0; iPop < indices.size(); ++iPop) {
+      indices[iPop] = util::opposite<DESCRIPTORBASE>(indices[iPop]);
+    }
+
+  }
+
+  std::vector<int> indices;
+
+  template <typename DESCRIPTORBASE_, int direction_, int orientation_>
+  friend std::vector<int> const& subIndexOutgoing();
+};
+
+template <typename DESCRIPTORBASE, int direction, int orientation>
+std::vector<int> const& subIndexOutgoing()
+{
+  static SubIndexOutgoing<DESCRIPTORBASE, direction, orientation> subIndexOutgoingSingleton;
+  return subIndexOutgoingSingleton.indices;
+}
+
+///finds all the remaining indexes of a lattice given some other indexes
+template <typename DESCRIPTORBASE>
+std::vector<int> remainingIndexes(const std::vector<int> &indices)
+{
+  std::vector<int> remaining;
+  for (int iPop = 0; iPop < DESCRIPTORBASE::q; ++iPop) {
+    bool found = false;
+    for (unsigned jPop = 0; jPop < indices.size(); ++jPop) {
+      if (indices[jPop] == iPop) {
+        found = true;
+      }
+    }
+    if (!found) {
+      remaining.push_back(iPop);
+    }
+  }
+  return remaining;
+}
+
+template <typename DESCRIPTORBASE, int plane, int normal1, int normal2>
+class SubIndexOutgoing3DonEdges {
+private:
+  SubIndexOutgoing3DonEdges()
+  {
+    int normalX,normalY,normalZ;
+    typedef DESCRIPTORBASE L;
+
+    switch (plane) {
+    case 0: {
+      normalX=0;
+      if (normal1==1) {
+        normalY= 1;
+      } else {
+        normalY=-1;
+      }
+      if (normal2==1) {
+        normalZ= 1;
+      } else {
+        normalZ=-1;
+      }
+    }
+    case 1: {
+      normalY=0;
+      if (normal1==1) {
+        normalX= 1;
+      } else {
+        normalX=-1;
+      }
+      if (normal2==1) {
+        normalZ= 1;
+      } else {
+        normalZ=-1;
+      }
+    }
+    case 2: {
+      normalZ=0;
+      if (normal1==1) {
+        normalX= 1;
+      } else {
+        normalX=-1;
+      }
+      if (normal2==1) {
+        normalY= 1;
+      } else {
+        normalY=-1;
+      }
+    }
+    }
+
+    // add zero velocity
+    //knownIndexes.push_back(0);
+    // compute scalar product with boundary normal for all other velocities
+    for (int iP=1; iP<L::q; ++iP) {
+      if (descriptors::c<L>(iP,0)*normalX + descriptors::c<L>(iP,1)*normalY + descriptors::c<L>(iP,2)*normalZ<0) {
+        indices.push_back(iP);
+      }
+    }
+  }
+  std::vector<int> indices;
+
+  template <typename DESCRIPTORBASE_,  int plane_, int normal1_, int normal2_>
+  friend std::vector<int> const& subIndexOutgoing3DonEdges();
+};
+
+template <typename DESCRIPTORBASE,  int plane, int normal1, int normal2>
+std::vector<int> const& subIndexOutgoing3DonEdges()
+{
+  static SubIndexOutgoing3DonEdges<DESCRIPTORBASE,  plane, normal1, normal2> subIndexOutgoing3DonEdgesSingleton;
+  return subIndexOutgoing3DonEdgesSingleton.indices;
+}
+
+// For 2D Corners
+template <typename DESCRIPTORBASE, int normalX, int normalY>
+class SubIndexOutgoing2DonCorners {
+private:
+  SubIndexOutgoing2DonCorners()
+  {
+    typedef DESCRIPTORBASE L;
+
+    // add zero velocity
+    //knownIndexes.push_back(0);
+    // compute scalar product with boundary normal for all other velocities
+    for (int iPop=1; iPop<L::q; ++iPop) {
+      if (descriptors::c<L>(iPop,0)*normalX + descriptors::c<L>(iPop,1)*normalY<0) {
+        indices.push_back(iPop);
+      }
+    }
+  }
+
+  std::vector<int> indices;
+
+  template <typename DESCRIPTORBASE_,  int normalX_, int normalY_>
+  friend std::vector<int> const& subIndexOutgoing2DonCorners();
+};
+
+template <typename DESCRIPTORBASE,  int normalX, int normalY>
+std::vector<int> const& subIndexOutgoing2DonCorners()
+{
+  static SubIndexOutgoing2DonCorners<DESCRIPTORBASE, normalX, normalY> subIndexOutgoing2DonCornersSingleton;
+  return subIndexOutgoing2DonCornersSingleton.indices;
+}
+
+// For 3D Corners
+template <typename DESCRIPTORBASE, int normalX, int normalY, int normalZ>
+class SubIndexOutgoing3DonCorners {
+private:
+  SubIndexOutgoing3DonCorners()
+  {
+    typedef DESCRIPTORBASE L;
+
+    // add zero velocity
+    //knownIndexes.push_back(0);
+    // compute scalar product with boundary normal for all other velocities
+    for (int iP=1; iP<L::q; ++iP) {
+      if (descriptors::c<L>(iP,0)*normalX + descriptors::c<L>(iP,1)*normalY + descriptors::c<L>(iP,2)*normalZ<0) {
+        indices.push_back(iP);
+      }
+    }
+  }
+
+  std::vector<int> indices;
+
+  template <typename DESCRIPTORBASE_,  int normalX_, int normalY_, int normalZ_>
+  friend std::vector<int> const& subIndexOutgoing3DonCorners();
+};
+
+template <typename DESCRIPTORBASE,  int normalX, int normalY, int normalZ>
+std::vector<int> const& subIndexOutgoing3DonCorners()
+{
+  static SubIndexOutgoing3DonCorners<DESCRIPTORBASE, normalX, normalY, normalZ> subIndexOutgoing3DonCornersSingleton;
+  return subIndexOutgoing3DonCornersSingleton.indices;
+}
+
+/// Util Function for Wall Model of Malaspinas
+/// get link with smallest angle to a vector
+template <typename T, typename DESCRIPTOR>
+int get_nearest_link(const std::vector<T>& vec)
+{
+  T max=-1;
+  int max_index = 0;
+  for (int iQ=1; iQ<DESCRIPTOR::q; ++iQ) {
+    std::vector<T> c_i(DESCRIPTOR::c[iQ], DESCRIPTOR::c[iQ]+3);
+    T tmp = util::scalarProduct<T>(c_i, vec)/util::norm(c_i);
+    if (tmp > max) {
+      max = tmp;
+      max_index = iQ;
+    }
+  }
+  return max_index;
+}
+
+namespace tensorIndices2D {
+enum { xx=0, xy=1, yy=2 };
+}
+
+namespace tensorIndices3D {
+enum { xx=0, xy=1, xz=2, yy=3, yz=4, zz=5 };
+}
+
+/// compute lattice density from lattice pressure
+template <typename T, typename DESCRIPTOR>
+T densityFromPressure( T latticePressure )
+{
+  // rho = p / c_s^2 + 1
+  return latticePressure * descriptors::invCs2<T,DESCRIPTOR>() + 1.0;
+}
+
+/// compute lattice pressure from lattice density
+template <typename T, typename DESCRIPTOR>
+T pressureFromDensity( T latticeDensity )
+{
+  // p = (rho - 1) * c_s^2
+  return (latticeDensity - 1.0) / descriptors::invCs2<T,DESCRIPTOR>();
+}
+
+}  // namespace util
+
+}  // namespace olb
+
+#endif