From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/core/util.h | 471 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 471 insertions(+)
create mode 100644 src/core/util.h
(limited to 'src/core/util.h')
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
+ *
+ *
+ * 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
+#include
+#include
+
+#include "dynamics/descriptorFunction.h"
+
+// patch due to ploblems with older compilers
+namespace std {
+template
+std::string to_string(const T &n)
+{
+ std::ostringstream s;
+ s << n;
+ return s.str();
+}
+}
+
+namespace olb {
+
+namespace util {
+
+template T norm(const std::vector& a);
+
+template
+inline int sign(T val)
+{
+ return (0 < val) - (val < 0);
+}
+
+template
+inline bool aligned_to_x(const std::vector& vec)
+{
+ return (vec[0]!=0 and vec[1]==0 and vec[2]==0);
+}
+
+template
+inline bool aligned_to_y(const std::vector& vec)
+{
+ return (vec[0]==0 and vec[1]!=0 and vec[2]==0);
+}
+
+template
+inline bool aligned_to_z(const std::vector& vec)
+{
+ return (vec[0]==0 and vec[1]==0 and vec[2]!=0);
+}
+
+template
+inline bool aligned_to_grid(const std::vector& vec)
+{
+ return (aligned_to_x(vec) or
+ aligned_to_y(vec) or
+ aligned_to_z(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
+T sqr(T arg)
+{
+ return arg*arg;
+}
+
+/// Compute norm square of a d-dimensional vector
+template
+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
+T normSqr(const Vector& u)
+{
+ T uSqr = T();
+ for (unsigned iD=0; iD < D; ++iD) {
+ uSqr += u[iD]*u[iD];
+ }
+ return uSqr;
+}
+
+template
+T scalarProduct(const T u1[d], const T u2[d])
+{
+ T prod = T();
+ for (int iD=0; iD
+T scalarProduct(const std::vector& u1, const std::vector& u2)
+{
+ T prod = T();
+ if (u1.size() == u2.size()) {
+ for (int iD=0; iD struct TensorVal {
+ static const int n =
+ (DESCRIPTORBASE::d*(DESCRIPTORBASE::d+1))/2; ///< result stored in n
+};
+
+/// Compute the opposite of a given direction
+template inline int opposite(int iPop)
+{
+ return descriptors::opposite(iPop);
+}
+
+template
+class SubIndex {
+private:
+ SubIndex()
+ {
+ for (int iVel=0; iVel(iVel,index)==value) {
+ indices.push_back(iVel);
+ }
+ }
+ }
+
+ std::vector indices;
+
+ template
+ friend std::vector const& subIndex();
+};
+
+template
+std::vector const& subIndex()
+{
+ static SubIndex subIndexSingleton;
+ return subIndexSingleton.indices;
+}
+
+template
+int findVelocity(const int v[DESCRIPTORBASE::d])
+{
+ for (int iPop=0; iPop(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
+class SubIndexOutgoing {
+private:
+ SubIndexOutgoing() // finds the indexes outgoing from the walls
+ {
+ indices = util::subIndex();
+
+ for (unsigned iPop = 0; iPop < indices.size(); ++iPop) {
+ indices[iPop] = util::opposite(indices[iPop]);
+ }
+
+ }
+
+ std::vector indices;
+
+ template
+ friend std::vector const& subIndexOutgoing();
+};
+
+template
+std::vector const& subIndexOutgoing()
+{
+ static SubIndexOutgoing subIndexOutgoingSingleton;
+ return subIndexOutgoingSingleton.indices;
+}
+
+///finds all the remaining indexes of a lattice given some other indexes
+template
+std::vector remainingIndexes(const std::vector &indices)
+{
+ std::vector 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
+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(iP,0)*normalX + descriptors::c(iP,1)*normalY + descriptors::c(iP,2)*normalZ<0) {
+ indices.push_back(iP);
+ }
+ }
+ }
+ std::vector indices;
+
+ template
+ friend std::vector const& subIndexOutgoing3DonEdges();
+};
+
+template
+std::vector const& subIndexOutgoing3DonEdges()
+{
+ static SubIndexOutgoing3DonEdges subIndexOutgoing3DonEdgesSingleton;
+ return subIndexOutgoing3DonEdgesSingleton.indices;
+}
+
+// For 2D Corners
+template
+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(iPop,0)*normalX + descriptors::c(iPop,1)*normalY<0) {
+ indices.push_back(iPop);
+ }
+ }
+ }
+
+ std::vector indices;
+
+ template
+ friend std::vector const& subIndexOutgoing2DonCorners();
+};
+
+template
+std::vector const& subIndexOutgoing2DonCorners()
+{
+ static SubIndexOutgoing2DonCorners subIndexOutgoing2DonCornersSingleton;
+ return subIndexOutgoing2DonCornersSingleton.indices;
+}
+
+// For 3D Corners
+template
+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(iP,0)*normalX + descriptors::c(iP,1)*normalY + descriptors::c(iP,2)*normalZ<0) {
+ indices.push_back(iP);
+ }
+ }
+ }
+
+ std::vector indices;
+
+ template
+ friend std::vector const& subIndexOutgoing3DonCorners();
+};
+
+template
+std::vector const& subIndexOutgoing3DonCorners()
+{
+ static SubIndexOutgoing3DonCorners subIndexOutgoing3DonCornersSingleton;
+ return subIndexOutgoing3DonCornersSingleton.indices;
+}
+
+/// Util Function for Wall Model of Malaspinas
+/// get link with smallest angle to a vector
+template
+int get_nearest_link(const std::vector& vec)
+{
+ T max=-1;
+ int max_index = 0;
+ for (int iQ=1; iQ c_i(DESCRIPTOR::c[iQ], DESCRIPTOR::c[iQ]+3);
+ T tmp = util::scalarProduct(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
+T densityFromPressure( T latticePressure )
+{
+ // rho = p / c_s^2 + 1
+ return latticePressure * descriptors::invCs2() + 1.0;
+}
+
+/// compute lattice pressure from lattice density
+template
+T pressureFromDensity( T latticeDensity )
+{
+ // p = (rho - 1) * c_s^2
+ return (latticeDensity - 1.0) / descriptors::invCs2();
+}
+
+} // namespace util
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3