multidatafields3d.hh

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2007 Bernd Stahl and Jonas Latt
 *  Address: Battelle Batiment A, Route de Drize 7, 1227 Carouge, Switzerland
 *  E-mail: [e-mail]
 *
 *  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
 * Scalar, vector and tensor fields for 3D data analysis -- generic implementation.
 */

#ifndef MULTI_DATA_FIELDS_3D_HH
#define MULTI_DATA_FIELDS_3D_HH

#include "multiDataFields3D.h"
#include "complexGrids/multiBlockStructure/multiDataGeometry3D.h"
#include <algorithm>
#include <limits>

namespace olb {

/////// Class MultiScalarField3D //////////////////////////////////

template<typename T>
MultiScalarField3D<T>::MultiScalarField3D(MultiDataDistribution3D const& dataDistribution_)
    : locatedBlock(0), constructed(false),
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      xSlice(dataDistribution_.getNy(), dataDistribution_.getNz() ),
      ySlice(dataDistribution_.getNx(), dataDistribution_.getNz() ), 
      zSlice(dataDistribution_.getNx(), dataDistribution_.getNy() )
{
#ifdef PARALLEL_MODE_MPI
    multiDataFieldHandler = new ParallelMultiDataFieldHandler3D<T>(dataDistribution_);
#else
    multiDataFieldHandler = new SerialMultiDataFieldHandler3D<T>(dataDistribution_);
#endif
    allocateFields();
}

template<typename T>
MultiScalarField3D<T>::~MultiScalarField3D() {
    deConstruct();
    deAllocateFields();
    delete multiDataFieldHandler;
    delete serializer; delete unSerializer;
}

template<typename T>
MultiScalarField3D<T>::MultiScalarField3D(MultiScalarField3D<T> const& rhs)
    : locatedBlock(rhs.locatedBlock),
      constructed(false), // must equal false to run this.construct(), below
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      xSlice(rhs.getMultiData().getNy(), rhs.getMultiData().getNz() ),
      ySlice(rhs.getMultiData().getNx(), rhs.getMultiData().getNz() ), 
      zSlice(rhs.getMultiData().getNx(), rhs.getMultiData().getNy() )
{
#ifdef PARALLEL_MODE_MPI
    multiDataFieldHandler = new ParallelMultiDataFieldHandler3D<T>( 
        rhs.multiDataFieldHandler->getMultiDataDistribution() );
#else
    multiDataFieldHandler = new SerialMultiDataFieldHandler3D<T>(
        rhs.multiDataFieldHandler->getMultiDataDistribution() );
#endif
    allocateFields();
    if (rhs.isConstructed()) {
        construct();
    }
    if (isConstructed()) {
        for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
            if (fields[iBlock]) {
                *(fields[iBlock]) = *(rhs.fields[iBlock]);
            }
        }
    }
}

template<typename T>
MultiScalarField3D<T>& MultiScalarField3D<T>::operator=(MultiScalarField3D<T> const& rhs) {
    MultiScalarField3D<T> tmp(rhs);
    swap(tmp);
    return *this;
}

template<typename T>
void MultiScalarField3D<T>::swap(MultiScalarField3D<T>& rhs) {
    std::swap(multiDataFieldHandler, rhs.multiDataFieldHandler);
    fields.swap(rhs.fields);
    std::swap(locatedBlock, rhs.locatedBlock);
    std::swap(constructed, rhs.constructed);
}

template<typename T>
bool MultiScalarField3D<T>::isConstructed() const {
    return constructed; 
}

template<typename T>
void MultiScalarField3D<T>::construct() {
    if (!isConstructed()) {
        for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
            if (fields[iBlock]) {
                fields[iBlock] -> construct();
            }
        }
        constructed = true;
    }
}

template<typename T>
void MultiScalarField3D<T>::deConstruct() {
    if (isConstructed()) {
        for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
            if (fields[iBlock]) {
                fields[iBlock] -> deConstruct();
            }
        }
        constructed = false;
    }
}

template<typename T>
void MultiScalarField3D<T>::reset() {
    OLB_PRECONDITION(isConstructed());
    for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
        if (fields[iBlock]) {
            fields[iBlock] -> reset();
        }
    }
}

template<typename T>
void MultiScalarField3D<T>::allocateFields() 
{
    for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
        int lx=0, ly=0, lz=0;
        if (multiDataFieldHandler->getLocalEnvelope(iBlock, lx, ly, lz)) {
            fields.push_back(new ScalarField3D<T>(lx,ly,lz));
        }
        else {
            fields.push_back( 0 );
        }
    }
}

template<typename T>
void MultiScalarField3D<T>::deAllocateFields() 
{
    for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
        delete fields[iBlock];
    }
}
template<typename T>
inline T& MultiScalarField3D<T>::get(int iX, int iY, int iZ) {
    OLB_PRECONDITION(iX>=0 && iX<getNx());
    OLB_PRECONDITION(iY>=0 && iY<getNy());
    OLB_PRECONDITION(iZ>=0 && iZ<getNz());
    OLB_PRECONDITION(isConstructed());

    locatedBlock = getMultiData().locate(iX, iY, iZ, locatedBlock);
    if (locatedBlock == -1) { 
        locatedBlock = 0;
        dummyScalar = std::numeric_limits<T>::signaling_NaN();
        return dummyScalar;
    }
    T* returnScalar = &dummyScalar;
    if (fields[locatedBlock]) { 
        returnScalar = &(fields[locatedBlock] -> get( 
                            getParameters(locatedBlock).toLocalX(iX),
                            getParameters(locatedBlock).toLocalY(iY),
                            getParameters(locatedBlock).toLocalZ(iZ) ));
    }
    multiDataFieldHandler -> broadCastScalar(*returnScalar, locatedBlock);
    return *returnScalar;
}

template<typename T>
inline T const& MultiScalarField3D<T>::get(int iX, int iY, int iZ) const {
    OLB_PRECONDITION(iX>=0 && iX<getNx());
    OLB_PRECONDITION(iY>=0 && iY<getNy());
    OLB_PRECONDITION(iZ>=0 && iZ<getNz());
    OLB_PRECONDITION(isConstructed());

    locatedBlock = getMultiData().locate(iX, iY, iZ, locatedBlock);
    if (locatedBlock == -1) { 
        locatedBlock = 0;
        dummyScalar = std::numeric_limits<T>::signaling_NaN();
        return dummyScalar;
    }
    T* returnScalar = &dummyScalar;
    if (fields[locatedBlock]) { 
        returnScalar = &(fields[locatedBlock] -> get( 
                            getParameters(locatedBlock).toLocalX(iX),
                            getParameters(locatedBlock).toLocalY(iY),
                            getParameters(locatedBlock).toLocalZ(iZ) ));
    }
    multiDataFieldHandler -> broadCastScalar(*returnScalar, locatedBlock);
    return *returnScalar;
}

template<typename T>
ScalarField2D<T> const& MultiScalarField3D<T>::sliceX(int xVal) const {
    xSlice.construct();
    copySerializedData (
            this->getSubSerializer(xVal, xVal, 0, getNy()-1, 0, getNz()-1, IndexOrdering::forward),
            xSlice.getUnSerializer(IndexOrdering::forward) );
    return xSlice;
}

template<typename T>
ScalarField2D<T> const& MultiScalarField3D<T>::sliceY(int yVal) const {
    ySlice.construct();
    copySerializedData (
            this->getSubSerializer(0, getNx()-1, yVal, yVal, 0, getNz()-1, IndexOrdering::forward),
            ySlice.getUnSerializer(IndexOrdering::forward) );
    return ySlice;
}

template<typename T>
ScalarField2D<T> const& MultiScalarField3D<T>::sliceZ(int zVal) const {
    zSlice.construct();
    copySerializedData (
            this->getSubSerializer(0, getNx()-1, 0, getNy()-1, zVal, zVal, IndexOrdering::backward),
            zSlice.getUnSerializer(IndexOrdering::backward) );
    return zSlice;
}

template<typename T>
DataSerializer<T> const& MultiScalarField3D<T>::getSerializer(IndexOrdering::OrderingT ordering) const
{
    delete serializer;
    serializer = new MultiSerializer3D<T>(serializerPolicy, ordering);
    return *serializer;
}

template<typename T>
DataUnSerializer<T>& MultiScalarField3D<T>::getUnSerializer(IndexOrdering::OrderingT ordering)
{
    delete unSerializer;
    unSerializer = new MultiUnSerializer3D<T>(unSerializerPolicy, ordering);
    return *unSerializer;
}

template<typename T>
DataSerializer<T> const& MultiScalarField3D<T>::getSubSerializer (
            int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
            IndexOrdering::OrderingT ordering ) const
{
    delete serializer;
    serializer = new MultiSerializer3D<T> (
            serializerPolicy, x0_, x1_, y0_, y1_, z0_, z1_, ordering );
    return *serializer;
}

template<typename T>
DataUnSerializer<T>& MultiScalarField3D<T>::getSubUnSerializer (
            int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
            IndexOrdering::OrderingT ordering )
{
    delete unSerializer;
    unSerializer = new MultiUnSerializer3D<T> (
            unSerializerPolicy, x0_, x1_, y0_, y1_, z0_, z1_, ordering );
    return *unSerializer;
}


template<typename T>
T MultiScalarField3D<T>::computeReduction(DataReduction<T>& reduction) const 
{
    OLB_PRECONDITION(isConstructed());
    reduction.reset();
    for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
        if (fields[iBlock]) {
            // data must be taken on the bulk only, not on the envelope
            BlockParameters3D blockParameters = getMultiData().getBlockParameters(iBlock);
            BlockCoordinates3D bulkLocal = blockParameters.toLocal( blockParameters.getBulk() );
            for (int iX=bulkLocal.x0; iX<=bulkLocal.x1; iX++) {
                for (int iY=bulkLocal.y0; iY<=bulkLocal.y1; iY++) {
                    for (int iZ=bulkLocal.z0; iZ<=bulkLocal.z1; iZ++) {
                        reduction.takeElement( fields[iBlock]->get(iX,iY,iZ) );
                    }
                }
            }
        }
    }
    reduction.reduceParallel();
    return reduction.getResult();
}

template<typename T>
MultiDataDistribution3D const& MultiScalarField3D<T>::getMultiData() const {
    return multiDataFieldHandler -> getMultiDataDistribution();
}

template<typename T>
BlockParameters3D const& MultiScalarField3D<T>::getParameters(int iParam) const {
    return getMultiData().getBlockParameters(iParam);
}

template<typename T>
int MultiScalarField3D<T>::getNumBlocks() const {
    return getMultiData().getNumBlocks();
}

template<typename T>
MultiDataDistribution3D MultiScalarField3D<T>::getDataDistribution() const {
    return getMultiData();
}

template<typename T>
SpatiallyExtendedObject3D* MultiScalarField3D<T>::getComponent(int iBlock) {
    OLB_PRECONDITION( iBlock<getScalarFields().size() );
    return getScalarFields()[iBlock];
}

template<typename T>
SpatiallyExtendedObject3D const* MultiScalarField3D<T>::getComponent(int iBlock) const {
    OLB_PRECONDITION( iBlock<getScalarFields().size() );
    return getScalarFields()[iBlock];
}

template<typename T>
multiPhysics::MultiPhysicsId MultiScalarField3D<T>::getMultiPhysicsId() const {
    return multiPhysics::getMultiPhysicsScalarId<T>();
}


//////// Class MultiTensorField3D //////////////////////////////////

template<typename T, int nDim>
MultiTensorField3D<T,nDim>::MultiTensorField3D(MultiDataDistribution3D const& dataDistribution_)
    : locatedBlock(0), constructed(false),
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      xSlice(dataDistribution_.getNy(), dataDistribution_.getNz() ),
      ySlice(dataDistribution_.getNx(), dataDistribution_.getNz() ), 
      zSlice(dataDistribution_.getNx(), dataDistribution_.getNy() )
{
#ifdef PARALLEL_MODE_MPI
    multiDataFieldHandler = new ParallelMultiDataFieldHandler3D<T>(dataDistribution_);
#else
    multiDataFieldHandler = new SerialMultiDataFieldHandler3D<T>(dataDistribution_);
#endif
    allocateFields();
    allocateComponents();
}

template<typename T, int nDim>
MultiTensorField3D<T,nDim>::~MultiTensorField3D() {
    deAllocateComponents();
    deConstruct();
    deAllocateFields();
    delete multiDataFieldHandler;
    delete serializer; delete unSerializer;
}

template<typename T, int nDim>
MultiTensorField3D<T,nDim>::MultiTensorField3D(MultiTensorField3D<T,nDim> const& rhs)
    : locatedBlock(rhs.locatedBlock), constructed(false), // must equal false to  run this.construct(), below
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      xSlice(rhs.getMultiData().getNy(), rhs.getMultiData().getNz() ),
      ySlice(rhs.getMultiData().getNx(), rhs.getMultiData().getNz() ), 
      zSlice(rhs.getMultiData().getNx(), rhs.getMultiData().getNy() )
{
#ifdef PARALLEL_MODE_MPI
    multiDataFieldHandler = new ParallelMultiDataFieldHandler3D<T>( 
        rhs.multiDataFieldHandler->getMultiDataDistribution() );
#else
    multiDataFieldHandler = new SerialMultiDataFieldHandler3D<T>(
        rhs.multiDataFieldHandler->getMultiDataDistribution() );
#endif
    allocateFields();
    allocateComponents();
    if (rhs.isConstructed()) {
        construct();
    }
    if (isConstructed()) {
        for (int iBlock=0; iBlock<getNumBlocks(); ++iBlock) {
            if (fields[iBlock]) {
                *(fields[iBlock]) = *(rhs.fields[iBlock]);
            }
        }
    }
}

template<typename T, int nDim>
MultiTensorField3D<T,nDim>& MultiTensorField3D<T,nDim>::operator=(MultiTensorField3D<T,nDim> const& rhs) {
    MultiTensorField3D<T,nDim> tmp(rhs);