summaryrefslogtreecommitdiff
path: root/src/functors/analytical/interpolationF2D.hh
diff options
context:
space:
mode:
Diffstat (limited to 'src/functors/analytical/interpolationF2D.hh')
-rw-r--r--src/functors/analytical/interpolationF2D.hh275
1 files changed, 275 insertions, 0 deletions
diff --git a/src/functors/analytical/interpolationF2D.hh b/src/functors/analytical/interpolationF2D.hh
new file mode 100644
index 0000000..12724df
--- /dev/null
+++ b/src/functors/analytical/interpolationF2D.hh
@@ -0,0 +1,275 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2012-2018 Lukas Baron, Tim Dornieden, Mathias J. Krause,
+ * Albert Mink, Adrian Kummerlaender
+ * 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_2D_HH
+#define INTERPOLATION_F_2D_HH
+
+#include "interpolationF2D.h"
+#include "core/superLattice2D.h"
+#include "dynamics/lbHelpers.h"
+
+namespace olb {
+
+
+template <typename T, typename W>
+AnalyticalFfromBlockF2D<T,W>::AnalyticalFfromBlockF2D(
+ BlockF2D<W>& f, Cuboid2D<T>& cuboid, const int overlap)
+ : AnalyticalF2D<T,W>(f.getTargetDim()),
+ _f(f), _cuboid(cuboid), _overlap(overlap)
+{
+ this->getName() = "fromBlockF";
+}
+
+template <typename T, typename W>
+bool AnalyticalFfromBlockF2D<T,W>::operator()(W output[], const T physC[])
+{
+ int latticeC[2];
+ int latticeR[2];
+ _cuboid.getFloorLatticeR(latticeR, physC);
+
+ if ( latticeR[0] >= -_overlap && latticeR[0] + 1 < _cuboid.getNx() + _overlap &&
+ latticeR[1] >= -_overlap && latticeR[1] + 1 < _cuboid.getNy() + _overlap ) {
+ const int& locX = latticeR[0];
+ const int& locY = latticeR[1];
+
+ Vector<T,2> physRiC;
+ Vector<T,2> physCv(physC);
+ _cuboid.getPhysR(physRiC.data, locX, locY);
+
+ // compute weights
+ Vector<W,2> d = (physCv - physRiC) * (1. / _cuboid.getDeltaR());
+ Vector<W,2> e = 1. - d;
+
+ T output_tmp[_f.getTargetDim()];
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output_tmp[iD] = T();
+ }
+
+ latticeC[0] = locX;
+ latticeC[1] = locY;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * e[0] * e[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[0] = locX;
+ latticeC[1] = locY + 1;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * e[0] * d[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[0] = locX + 1;
+ latticeC[1] = locY;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * d[0] * e[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[0] = locX + 1;
+ latticeC[1] = locY + 1;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * d[0] * d[1];
+ output_tmp[iD] = T();
+ }
+
+ return true;
+ }
+ else {
+ return false;
+ }
+}
+
+template <typename T, typename W>
+AnalyticalFfromSuperF2D<T,W>::AnalyticalFfromSuperF2D(SuperF2D<T>& f,
+ bool communicateToAll, int overlap, bool communicateOverlap)
+ : AnalyticalF2D<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 AnalyticalFfromBlockF2D<T>(_f.getBlockF(iC),
+ _cuboidGeometry.get(load.glob(iC)),
+ _overlap)
+ );
+ }
+ }
+}
+
+template <typename T, typename W>
+bool AnalyticalFfromSuperF2D<T,W>::operator() (T output[], const T physC[])
+{
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] = W();
+ }
+
+ int latticeR[3];
+ if (!_cuboidGeometry.getLatticeR(latticeR, physC)) {
+ return false;
+ }
+
+ if (_communicateOverlap) {
+ _f.getSuperStructure().communicate();
+ }
+
+ int dataSize = 0;
+ int dataFound = 0;
+
+ int latticeC[3] = {};
+
+ LoadBalancer<T>& load = _f.getSuperStructure().getLoadBalancer();
+
+ for (int iC = 0; iC < load.size(); ++iC) {
+ latticeC[0] = load.glob(iC);
+ Cuboid2D<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 ) {
+ if (_blockF.empty()) {
+ const int& locX = latticeR[0];
+ const int& locY = latticeR[1];
+
+ Vector<T,2> physRiC;
+ Vector<T,2> physCv(physC);
+ cuboid.getPhysR(physRiC.data, locX, locY);
+
+ // compute weights
+ Vector<W,2> d = (physCv - physRiC) * (1. / cuboid.getDeltaR());
+ Vector<W,2> e = 1. - d;
+
+ T output_tmp[_f.getTargetDim()];
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output_tmp[iD] = T();
+ }
+
+ latticeC[1] = locX;
+ latticeC[2] = locY;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * e[0] * e[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[1] = locX;
+ latticeC[2] = locY + 1;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * e[0] * d[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[1] = locX + 1;
+ latticeC[2] = locY;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * d[0] * e[1];
+ output_tmp[iD] = T();
+ }
+
+ latticeC[1] = locX + 1;
+ latticeC[2] = locY + 1;
+ _f(output_tmp,latticeC);
+ for (int iD = 0; iD < _f.getTargetDim(); ++iD) {
+ output[iD] += output_tmp[iD] * d[0] * d[1];
+ output_tmp[iD] = T();
+ }
+ }
+ 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 AnalyticalFfromSuperF2D<T,W>::getBlockFSize() const
+{
+ OLB_ASSERT(_blockF.size() < UINT32_MAX,
+ "it is safe to cast std::size_t to int");
+ return _blockF.size();
+}
+
+template <typename T, typename W>
+AnalyticalFfromBlockF2D<T,W>& AnalyticalFfromSuperF2D<T,W>::getBlockF(int iCloc)
+{
+ OLB_ASSERT(size_t(iCloc) < _blockF.size() && iCloc >= 0,
+ "block functor index within bounds");
+ return *(_blockF[iCloc]);
+}
+
+
+} // end namespace olb
+
+#endif