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/functors/analytical/interpolationF3D.hh | 499 ++++++++++++++++++++++++++++
 1 file changed, 499 insertions(+)
 create mode 100644 src/functors/analytical/interpolationF3D.hh

(limited to 'src/functors/analytical/interpolationF3D.hh')

diff --git a/src/functors/analytical/interpolationF3D.hh b/src/functors/analytical/interpolationF3D.hh
new file mode 100644
index 0000000..5075633
--- /dev/null
+++ b/src/functors/analytical/interpolationF3D.hh
@@ -0,0 +1,499 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2012-2017 Lukas Baron, Tim Dornieden, Mathias J. Krause,
+ *  Albert Mink, Fabian Klemens, Benjamin Förster, Marie-Luise Maier,
+ *  Adrian Kummerlönder
+ *  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.
+ */
+
+#ifndef INTERPOLATION_F_3D_HH
+#define INTERPOLATION_F_3D_HH
+
+#include <algorithm>
+
+#include "interpolationF3D.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+
+namespace olb {
+
+
+/// trilinear interpolation for rectangular lattice with dimensions delta[i];
+/// if the cuboid is a plane (e.g. nZ==1) this functor will convert to bilinear interpolation and to linear interpolation for a "line cuboid" (e.g. nY==nZ==1)
+template <typename T, typename W>
+SpecialAnalyticalFfromBlockF3D<T,W>::SpecialAnalyticalFfromBlockF3D(
+  BlockF3D<W>& f, Cuboid3D<T>& cuboid,
+  Vector<T,3> delta, T scale)
+  : AnalyticalF3D<T,W>(f.getTargetDim()), _f(f), _cuboid(cuboid), _delta(delta), _scale(scale)
+{
+  this->getName() = "fromBlockF";
+}
+
+
+template <typename T, typename W>
+bool SpecialAnalyticalFfromBlockF3D<T,W>::operator()(W output[],
+    const T physC[])
+{
+  Vector<T,3> origin = _cuboid.getOrigin();
+
+  // scale physC in all 3 dimensions
+  Vector<T,3> physCv;
+  for (int i=0; i<3; i++) {
+    physCv[i] = origin[i] + (physC[i] - origin[i]) * ( _cuboid.getDeltaR() / _delta[i] );
+  }
+
+  int latticeR[3];
+  for (int i=0; i<3; i++) {
+    latticeR[i] = std::max((int)floor( (physCv[i] - origin[i])/
+                                       _cuboid.getDeltaR()), 0);
+  }
+  Vector<T,3> physRiC;
+  Vector<W,3> d, e;
+  W output_tmp[3];
+  Vector<T,3> latticeRv;
+
+  for (int i=0; i<3; i++) {
+    latticeRv[i] = (T) latticeR[i];
+  }
+  physRiC = origin + latticeRv * _cuboid.getDeltaR();
+  T dr = 1. / _cuboid.getDeltaR();
+
+  // compute weights
+  d = (physCv - physRiC) * dr;
+  e = 1. - d;
+
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] = W();
+    output_tmp[iD] = W();
+  }
+
+  //0=1=2=
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * e[0]*e[1]*e[2];
+  }
+
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]++;
+  }
+
+  //0=1+2=
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * e[0]*d[1]*e[2];
+  }
+
+  if (_cuboid.getNx() != 1) {
+    latticeR[0]++;
+  }
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]--;
+  }
+  //0+1=2=
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * d[0]*e[1]*e[2];
+  }
+
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]++;
+  }
+  //0+1+2=
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * d[0]*d[1]*e[2];
+  }
+
+  if (_cuboid.getNx() != 1) {
+    latticeR[0]--;
+  }
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]--;
+  }
+  if (_cuboid.getNz() != 1) {
+    latticeR[2]++;
+  }
+  //0=1=2+
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * e[0]*e[1]*d[2];
+  }
+
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]++;
+  }
+  //0=1+2+
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * e[0]*d[1]*d[2];
+  }
+
+  if (_cuboid.getNx() != 1) {
+    latticeR[0]++;
+  }
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]--;
+  }
+  //0+1=2+
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * d[0]*e[1]*d[2];
+  }
+
+  if (_cuboid.getNy() != 1) {
+    latticeR[1]++;
+  }
+  //0+1+2+
+  _f(output_tmp, latticeR);
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] += output_tmp[iD] * d[0]*d[1]*d[2];
+  }
+
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] *= _scale;
+  }
+
+  return true;
+}
+
+
+template <typename T, typename W>
+AnalyticalFfromBlockF3D<T,W>::AnalyticalFfromBlockF3D(
+  BlockF3D<W>& f, Cuboid3D<T>& cuboid, const int overlap)
+  : AnalyticalF3D<T,W>(f.getTargetDim()),
+    _f(f), _cuboid(cuboid), _overlap(overlap)
+{
+  this->getName() = "fromBlockF";
+}
+
+/// trilinear interpolation on cubic lattice
+template <typename T, typename W>
+bool AnalyticalFfromBlockF3D<T,W>::operator()(W output[], const T physC[])
+{
+  int latticeC[3];
+  int latticeR[3];
+  _cuboid.getFloorLatticeR(latticeR, physC);
+
+  if ( latticeR[0] >= -_overlap && latticeR[0] + 1 < _cuboid.getNx() + _overlap &&
+       latticeR[1] >= -_overlap && latticeR[1] + 1 < _cuboid.getNy() + _overlap &&
+       latticeR[2] >= -_overlap && latticeR[2] + 1 < _cuboid.getNz() + _overlap ) {
+    const int& locX = latticeR[0];
+    const int& locY = latticeR[1];
+    const int& locZ = latticeR[2];
+
+    Vector<T,3> physRiC;
+    Vector<T,3> physCv(physC);
+    _cuboid.getPhysR(physRiC.data, locX, locY, locZ);
+
+    // compute weights
+    Vector<W,3> d = (physCv - physRiC) * (1. / _cuboid.getDeltaR());
+    Vector<W,3> e = 1. - d;
+
+    W output_tmp[_f.getTargetDim()];
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX;
+    latticeC[1] = locY;
+    latticeC[2] = locZ;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * e[0] * e[1] * e[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX;
+    latticeC[1] = locY + 1;
+    latticeC[2] = locZ;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * e[0] * d[1] * e[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX + 1;
+    latticeC[1] = locY;
+    latticeC[2] = locZ;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * d[0] * e[1] * e[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX + 1;
+    latticeC[1] = locY + 1;
+    latticeC[2] = locZ;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * d[0] * d[1] * e[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX;
+    latticeC[1] = locY;
+    latticeC[2] = locZ + 1;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * e[0] * e[1] * d[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX;
+    latticeC[1] = locY + 1;
+    latticeC[2] = locZ + 1;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * e[0] * d[1] * d[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX + 1;
+    latticeC[1] = locY;
+    latticeC[2] = locZ + 1;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * d[0] * e[1] * d[2];
+      output_tmp[iD] = W();
+    }
+
+    latticeC[0] = locX + 1;
+    latticeC[1] = locY + 1;
+    latticeC[2] = locZ + 1;
+    _f(output_tmp,latticeC);
+    for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+      output[iD] += output_tmp[iD] * d[0] * d[1] * d[2];
+      output_tmp[iD] = W();
+    }
+
+    return true;
+  } else {
+    return false;
+  }
+}
+
+
+template <typename T, typename W>
+AnalyticalFfromSuperF3D<T,W>::AnalyticalFfromSuperF3D(SuperF3D<T,W>& f,
+    bool communicateToAll, int overlap, bool communicateOverlap)
+  : AnalyticalF3D<T,W>(f.getTargetDim()),
+    _communicateToAll(communicateToAll),
+    _communicateOverlap(communicateOverlap),
+    _f(f),
+    _cuboidGeometry(f.getSuperStructure().getCuboidGeometry()),
+    _overlap(overlap)
+{
+  this->getName() = "fromSuperF";
+
+  if (overlap == -1) {
+    _overlap = _f.getSuperStructure().getOverlap();
+  }
+
+  LoadBalancer<T>& load = _f.getSuperStructure().getLoadBalancer();
+
+  if ( _f.getBlockFSize() == load.size() ) {
+    for (int iC = 0; iC < load.size(); ++iC) {
+      this->_blockF.emplace_back(
+        new AnalyticalFfromBlockF3D<T>(_f.getBlockF(iC),
+                                       _cuboidGeometry.get(load.glob(iC)),
+                                       _overlap)
+      );
+    }
+  }
+}
+
+template <typename T, typename W>
+bool AnalyticalFfromSuperF3D<T,W>::operator()(W output[], const T physC[])
+{
+  for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+    output[iD] = W();
+  }
+
+  int latticeR[4];
+  if (!_cuboidGeometry.getLatticeR(latticeR, physC)) {
+    return false;
+  }
+
+  if (_communicateOverlap) {
+    _f.getSuperStructure().communicate();
+  }
+
+  int dataSize = 0;
+  int dataFound = 0;
+
+  int latticeC[4] = {};
+
+  LoadBalancer<T>& load = _f.getSuperStructure().getLoadBalancer();
+
+  for (int iC = 0; iC < load.size(); ++iC) {
+    latticeC[0] = load.glob(iC);
+    Cuboid3D<T>& cuboid = _cuboidGeometry.get(latticeC[0]);
+    cuboid.getFloorLatticeR(latticeR, physC);
+
+    // latticeR within cuboid extended by overlap
+    if ( latticeR[0] >= -_overlap && latticeR[0] + 1 < cuboid.getNx() + _overlap &&
+         latticeR[1] >= -_overlap && latticeR[1] + 1 < cuboid.getNy() + _overlap &&
+         latticeR[2] >= -_overlap && latticeR[2] + 1 < cuboid.getNz() + _overlap ) {
+      if (_blockF.empty()) {
+        const int& locX = latticeR[0];
+        const int& locY = latticeR[1];
+        const int& locZ = latticeR[2];
+
+        Vector<T,3> physRiC;
+        Vector<T,3> physCv(physC);
+        cuboid.getPhysR(physRiC.data, locX, locY, locZ);
+
+        // compute weights
+        Vector<W,3> d = (physCv - physRiC) * (1. / cuboid.getDeltaR());
+        Vector<W,3> e = 1. - d;
+
+        W output_tmp[_f.getTargetDim()];
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX;
+        latticeC[2] = locY;
+        latticeC[3] = locZ;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * e[0] * e[1] * e[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX;
+        latticeC[2] = locY + 1;
+        latticeC[3] = locZ;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * e[0] * d[1] * e[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX + 1;
+        latticeC[2] = locY;
+        latticeC[3] = locZ;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * d[0] * e[1] * e[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX + 1;
+        latticeC[2] = locY + 1;
+        latticeC[3] = locZ;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * d[0] * d[1] * e[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX;
+        latticeC[2] = locY;
+        latticeC[3] = locZ + 1;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * e[0] * e[1] * d[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX;
+        latticeC[2] = locY + 1;
+        latticeC[3] = locZ + 1;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * e[0] * d[1] * d[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX + 1;
+        latticeC[2] = locY;
+        latticeC[3] = locZ + 1;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * d[0] * e[1] * d[2];
+          output_tmp[iD] = W();
+        }
+
+        latticeC[1] = locX + 1;
+        latticeC[2] = locY + 1;
+        latticeC[3] = locZ + 1;
+        _f(output_tmp,latticeC);
+        for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+          output[iD] += output_tmp[iD] * d[0] * d[1] * d[2];
+          output_tmp[iD] = W();
+        }
+      } else {
+        _blockF[iC]->operator()(output, physC);
+      }
+
+      dataSize += _f.getTargetDim();
+      ++dataFound;
+    }
+  }
+
+  if (_communicateToAll) {
+#ifdef PARALLEL_MODE_MPI
+    singleton::mpi().reduceAndBcast(dataFound, MPI_SUM);
+    singleton::mpi().reduceAndBcast(dataSize, MPI_SUM);
+#endif
+    dataSize /= dataFound;
+#ifdef PARALLEL_MODE_MPI
+    for (int iD = 0; iD < dataSize; ++iD) {
+      singleton::mpi().reduceAndBcast(output[iD], MPI_SUM);
+    }
+#endif
+    for (int iD = 0; iD < dataSize; ++iD) {
+      output[iD]/=dataFound;
+    }
+  } else {
+    if (dataFound!=0) {
+      dataSize /= dataFound;
+      for (int iD = 0; iD < dataSize; ++iD) {
+        output[iD]/=dataFound;
+      }
+    }
+  }
+
+  if (dataFound>0) {
+    return true;
+  }
+  return false;
+}
+
+template <typename T, typename W>
+int AnalyticalFfromSuperF3D<T,W>::getBlockFSize() const
+{
+  OLB_ASSERT(_blockF.size() < INT32_MAX,
+             "it is safe to cast std::size_t to int");
+  return _blockF.size();
+}
+
+template <typename T, typename W>
+AnalyticalFfromBlockF3D<T,W>& AnalyticalFfromSuperF3D<T,W>::getBlockF(int iCloc)
+{
+  OLB_ASSERT(iCloc < int(_blockF.size()) && iCloc >= 0,
+             "block functor index within bounds");
+  return *(_blockF[iCloc]);
+}
+
+
+}  // end namespace olb
+
+#endif
-- 
cgit v1.2.3