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