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
---
.../lattice/integral/superPlaneIntegralF3D.hh | 322 +++++++++++++++++++++
1 file changed, 322 insertions(+)
create mode 100644 src/functors/lattice/integral/superPlaneIntegralF3D.hh
(limited to 'src/functors/lattice/integral/superPlaneIntegralF3D.hh')
diff --git a/src/functors/lattice/integral/superPlaneIntegralF3D.hh b/src/functors/lattice/integral/superPlaneIntegralF3D.hh
new file mode 100644
index 0000000..a99a064
--- /dev/null
+++ b/src/functors/lattice/integral/superPlaneIntegralF3D.hh
@@ -0,0 +1,322 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2017 Adrian Kummerlaender
+ * 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.
+*/
+
+#ifndef SUPER_PLANE_INTEGRAL_F_3D_HH
+#define SUPER_PLANE_INTEGRAL_F_3D_HH
+
+#include "superPlaneIntegralF3D.h"
+#include "utilities/vectorHelpers.h"
+#include "functors/analytical/indicator/indicator2D.hh"
+
+namespace olb {
+
+
+template
+bool SuperPlaneIntegralF3D::isToBeIntegrated(const Vector& physR, int iC)
+{
+ Vector latticeR;
+ //get nearest lattice point
+ if ( _geometry.getCuboidGeometry().getFloorLatticeR(physR, latticeR) ) {
+ const int& iX = latticeR[1];
+ const int& iY = latticeR[2];
+ const int& iZ = latticeR[3];
+
+ // interpolation is possible iff all neighbours are within the indicated subset
+ return _integrationIndicatorF->operator()( iC, iX, iY, iZ )
+ && _integrationIndicatorF->operator()(iC, iX, iY, iZ+1)
+ && _integrationIndicatorF->operator()(iC, iX, iY+1, iZ )
+ && _integrationIndicatorF->operator()(iC, iX, iY+1, iZ+1)
+ && _integrationIndicatorF->operator()(iC, iX+1, iY, iZ )
+ && _integrationIndicatorF->operator()(iC, iX+1, iY, iZ+1)
+ && _integrationIndicatorF->operator()(iC, iX+1, iY+1, iZ )
+ && _integrationIndicatorF->operator()(iC, iX+1, iY+1, iZ+1);
+ }
+ else {
+ return false;
+ }
+}
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const HyperplaneLattice3D& hyperplaneLattice,
+ FunctorPtr>&& integrationIndicator,
+ FunctorPtr>&& subplaneIndicator,
+ BlockDataReductionMode mode)
+ : SuperF3D(f->getSuperStructure(), 2 + f->getTargetDim()),
+ _geometry(geometry),
+ _f(std::move(f)),
+ _integrationIndicatorF(std::move(integrationIndicator)),
+ _subplaneIndicatorF(std::move(subplaneIndicator)),
+ _reductionF(*_f,
+ hyperplaneLattice,
+ BlockDataSyncMode::None,
+ mode),
+ _origin(hyperplaneLattice.getHyperplane().origin),
+ _u(hyperplaneLattice.getVectorU()),
+ _v(hyperplaneLattice.getVectorV()),
+ _normal(hyperplaneLattice.getHyperplane().normal)
+{
+ this->getName() = "SuperPlaneIntegralF3D";
+
+ _normal.normalize();
+ _u.normalize();
+ _v.normalize();
+
+ for ( const std::tuple& pos : _reductionF.getRankLocalSubplane() ) {
+ const int& i = std::get<0>(pos);
+ const int& j = std::get<1>(pos);
+ const int& iC = std::get<2>(pos);
+ const Vector physR = _reductionF.getPhysR(i, j);
+ if (isToBeIntegrated(physR, iC)) {
+ // check if interpolated hyperplane is to be restricted further
+ // e.g. using IndicatorCircle2D
+ if ( _subplaneIndicatorF ) {
+ // determine physical coordinates relative to original hyperplane origin
+ // [!] different from _reductionF._origin in the general case.
+ const Vector physRelativeToOrigin = physR - _origin;
+ const T physOnHyperplane[2] {
+ physRelativeToOrigin * _u,
+ physRelativeToOrigin * _v
+ };
+
+ if ( _subplaneIndicatorF->operator()(physOnHyperplane) ) {
+ _rankLocalSubplane.emplace_back(i, j);
+ }
+ }
+ else {
+ // plane is not restricted further
+ _rankLocalSubplane.emplace_back(i, j);
+ }
+ }
+ }
+}
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Hyperplane3D& hyperplane,
+ FunctorPtr>&& integrationIndicator,
+ FunctorPtr>&& subplaneIndicator,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ HyperplaneLattice3D(geometry.getCuboidGeometry(), hyperplane),
+ std::forward(integrationIndicator),
+ std::forward(subplaneIndicator),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Hyperplane3D& hyperplane,
+ FunctorPtr>&& integrationIndicator,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ hyperplane,
+ std::forward(integrationIndicator),
+ nullptr,
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& origin, const Vector& u, const Vector& v,
+ std::vector materials,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ Hyperplane3D().originAt(origin).spannedBy(u, v),
+ geometry.getMaterialIndicator(std::forward(materials)),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& origin, const Vector& u, const Vector& v,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ origin, u, v,
+ std::vector(1,1),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& origin, const Vector& normal,
+ std::vector materials,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ Hyperplane3D().originAt(origin).normalTo(normal),
+ geometry.getMaterialIndicator(std::forward(materials)),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& origin, const Vector& normal,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ origin, normal,
+ std::vector(1,1),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& normal,
+ std::vector materials,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ Hyperplane3D()
+ .centeredIn(geometry.getCuboidGeometry().getMotherCuboid())
+ .normalTo(normal),
+ geometry.getMaterialIndicator(std::forward(materials)),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const Vector& normal,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ normal,
+ std::vector(1,1),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const IndicatorCircle3D& circle, std::vector materials,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ Hyperplane3D().originAt(circle.getCenter()).normalTo(circle.getNormal()),
+ geometry.getMaterialIndicator(std::forward>(materials)),
+ std::unique_ptr>(new IndicatorCircle2D({0,0}, circle.getRadius())),
+ mode)
+{ }
+
+template
+SuperPlaneIntegralF3D::SuperPlaneIntegralF3D(
+ FunctorPtr>&& f,
+ SuperGeometry3D& geometry,
+ const IndicatorCircle3D& circle,
+ BlockDataReductionMode mode)
+ : SuperPlaneIntegralF3D(
+ std::forward(f),
+ geometry,
+ circle,
+ std::vector(1,1),
+ mode)
+{ }
+
+template
+bool SuperPlaneIntegralF3D::operator()(T output[], const int input[])
+{
+ this->getSuperStructure().communicate();
+
+ _reductionF.update();
+
+ const int flowDim = _reductionF.getTargetDim();
+
+ std::vector flow(flowDim,0.);
+
+ for ( std::tuple& pos : _rankLocalSubplane ) {
+ T outputTmp[flowDim];
+ const int inputTmp[2] { std::get<0>(pos), std::get<1>(pos) };
+
+ _reductionF(outputTmp, inputTmp);
+
+ for ( int j = 0; j < flowDim; j++ ) {
+ flow[j] += outputTmp[j];
+ }
+ }
+
+ int vox = _rankLocalSubplane.size();
+
+#ifdef PARALLEL_MODE_MPI
+ for ( int j = 0; j < flowDim; j++ ) {
+ singleton::mpi().reduceAndBcast(flow[j], MPI_SUM);
+ }
+ singleton::mpi().reduceAndBcast(vox, MPI_SUM);
+#endif
+
+ const T h = _reductionF.getPhysSpacing();
+
+ switch ( flowDim ) {
+ case 1: {
+ output[0] = flow[0] * h * h;
+ break;
+ }
+ case 3: {
+ output[0] = (h*h * Vector(flow)) * _normal;
+ break;
+ }
+ }
+
+ // area
+ output[1] = vox * h * h;
+ // write flow to output[2..]
+ std::copy_n(flow.cbegin(), flowDim, &output[2]);
+
+ return true;
+}
+
+
+}
+
+#endif
--
cgit v1.2.3