offLatticeBoundaryCondition3D.hh

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OFF_LATTICE_BOUNDARY_CONDITION_3D_HH
#define OFF_LATTICE_BOUNDARY_CONDITION_3D_HH

#include "core/globalDefs.h"
#include "offLattice/offLatticeBoundaryCondition3D.h"
#include "offLattice/triangularSurfaceMesh.h"
#include "offLattice/offLatticeBoundaryProfiles3D.h"
#include "triangleToDef.h"

namespace plb {


/* ********** OffLatticeBoundaryCondition3D ********************************** */

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::
    OffLatticeBoundaryCondition3D (
        OffLatticeModel3D<T,BoundaryType>* offLatticeModel_,
        VoxelizedDomain3D<T>& voxelizedDomain_,
        MultiBlockLattice3D<T,Descriptor>& lattice_ )
    : voxelizedDomain(voxelizedDomain_),
      lattice(lattice_),
      boundaryShapeArg(lattice_),
      offLatticeModel(offLatticeModel_),
      offLatticePattern(lattice)
{
    std::vector<MultiBlock3D*> offLatticeIniArg;
    // First argument for compute-off-lattice-pattern.
    offLatticeIniArg.push_back(&offLatticePattern);
    // Remaining arguments for inner-flow-shape.
    offLatticeIniArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeIniArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeIniArg.push_back(&boundaryShapeArg);
    applyProcessingFunctional (
            new OffLatticePatternFunctional3D<T,BoundaryType> (
                offLatticeModel->clone() ),
            offLatticePattern.getBoundingBox(), offLatticeIniArg );
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::
    OffLatticeBoundaryCondition3D (
        OffLatticeModel3D<T,BoundaryType>* offLatticeModel_,
        VoxelizedDomain3D<T>& voxelizedDomain_,
        MultiBlockLattice3D<T,Descriptor>& lattice_,
        MultiParticleField3D<DenseParticleField3D<T,Descriptor> >& particleField_ )
    : offLatticeModel(offLatticeModel_),
      voxelizedDomain(voxelizedDomain_),
      lattice(lattice_),
      boundaryShapeArg(particleField_),
      offLatticePattern(lattice)
{
    std::vector<MultiBlock3D*> offLatticeIniArg;
    // First argument for compute-off-lattice-pattern.
    offLatticeIniArg.push_back(&offLatticePattern);
    // Remaining arguments for inner-flow-shape.
    offLatticeIniArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeIniArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeIniArg.push_back(&boundaryShapeArg);
    applyProcessingFunctional (
            new OffLatticePatternFunctional3D<T,BoundaryType> (
                offLatticeModel->clone() ),
            offLatticePattern.getBoundingBox(), offLatticeIniArg );
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::OffLatticeBoundaryCondition3D (
        OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType> const& rhs )
    : offLatticeModel(rhs.offLatticeModel.clone()),
      voxelizedDomain(rhs.voxelizedDomain),
      lattice(rhs.lattice),
      boundaryShapeArg(rhs.boundaryShapeArg),
      offLatticePattern(rhs.offLatticePattern)
{ }


template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::~OffLatticeBoundaryCondition3D()
{
    delete offLatticeModel;
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
void OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::insert()
{
    std::vector<MultiBlock3D*> offLatticeArg;
    // First three arguments for Guo algorithm.
    offLatticeArg.push_back(&lattice);
    offLatticeArg.push_back(&offLatticePattern);
    // Remaining arguments for inner-flow-shape.
    offLatticeArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeArg.push_back(&boundaryShapeArg);
    plint processorLevel = 1;
    plint numShapeArgs = 3;
    plint numCompletionArgs = 0;
    integrateProcessingFunctional (
            new OffLatticeCompletionFunctional3D<T,Descriptor,BoundaryType> (
                offLatticeModel->clone(), numShapeArgs, numCompletionArgs ),
            boundaryShapeArg.getBoundingBox(), offLatticeArg, processorLevel );
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
void OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::insert (
        std::vector<MultiBlock3D*> const& completionArg )
{
    std::vector<MultiBlock3D*> offLatticeArg;
    // First three arguments for Guo algorithm.
    offLatticeArg.push_back(&lattice);
    offLatticeArg.push_back(&offLatticePattern);
    // Next arguments for inner-flow-shape.
    offLatticeArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeArg.push_back(&boundaryShapeArg);
    // Remaining are optional arguments for completion algorithm.
    plint numCompletionArgs = (plint)completionArg.size();
    for (plint i=0; i<numCompletionArgs; ++i) {
        offLatticeArg.push_back(completionArg[i]);
    }
    plint processorLevel = 1;
    plint numShapeArgs = 3;
    integrateProcessingFunctional (
            new OffLatticeCompletionFunctional3D<T,Descriptor,BoundaryType> (
                offLatticeModel->clone(), numShapeArgs, numCompletionArgs ),
            boundaryShapeArg.getBoundingBox(), offLatticeArg, processorLevel );
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
void OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::apply()
{
    std::vector<MultiBlock3D*> offLatticeArg;
    // First three arguments for Guo algorithm.
    offLatticeArg.push_back(&lattice);
    offLatticeArg.push_back(&offLatticePattern);
    // Remaining arguments for inner-flow-shape.
    offLatticeArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeArg.push_back(&boundaryShapeArg);
    plint numShapeArgs = 3;
    plint numCompletionArgs = 0;
    applyProcessingFunctional (
            new OffLatticeCompletionFunctional3D<T,Descriptor,BoundaryType> (
                offLatticeModel->clone(), numShapeArgs, numCompletionArgs ),
            boundaryShapeArg.getBoundingBox(), offLatticeArg );
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
void OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::apply (
        std::vector<MultiBlock3D*> const& completionArg )
{
    std::vector<MultiBlock3D*> offLatticeArg;
    // First three arguments for Guo algorithm.
    offLatticeArg.push_back(&lattice);
    offLatticeArg.push_back(&offLatticePattern);
    // Next arguments for inner-flow-shape.
    offLatticeArg.push_back(&voxelizedDomain.getVoxelMatrix());
    offLatticeArg.push_back(&voxelizedDomain.getTriangleHash());
    offLatticeArg.push_back(&boundaryShapeArg);
    // Remaining are optional arguments for completion algorithm.
    plint numCompletionArgs = (plint)completionArg.size();
    for (plint i=0; i<numCompletionArgs; ++i) {
        offLatticeArg.push_back(completionArg[i]);
    }
    plint numShapeArgs = 3;
    applyProcessingFunctional (
            new OffLatticeCompletionFunctional3D<T,Descriptor,BoundaryType> (
                offLatticeModel->clone(), numShapeArgs, numCompletionArgs ),
            boundaryShapeArg.getBoundingBox(), offLatticeArg );
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
Array<T,3> OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::getForceOnObject()
{
    std::vector<MultiBlock3D*> arg;
    arg.push_back(&offLatticePattern);
    GetForceOnObjectFunctional3D<T,BoundaryType> functional(offLatticeModel->clone());
    applyProcessingFunctional (
            functional, boundaryShapeArg.getBoundingBox(), arg );
    return functional.getForce();
}

    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiTensorField3D<T,3> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocity(Box3D domain)
{
    std::auto_ptr<MultiTensorField3D<T,3> > velocity(plb::computeVelocity(lattice,domain));
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant<T,3>(*velocity, voxelizedDomain.getVoxelMatrix(), solidFlag,
                       domain, Array<T,3>(T(),T(),T()));
    setToConstant<T,3>(*velocity, voxelizedDomain.getVoxelMatrix(), solidBorderFlag,
                       domain, Array<T,3>(T(),T(),T()));
    return velocity;
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiTensorField3D<T,3> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocity()
{
    return computeVelocity(lattice.getBoundingBox());
}

    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiTensorField3D<T,3> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVorticity(Box3D domain)
{
    std::auto_ptr<MultiTensorField3D<T,3> > vorticity (
            plb::computeBulkVorticity(*plb::computeVelocity(lattice,domain), domain) );
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant<T,3>(*vorticity, voxelizedDomain.getVoxelMatrix(), solidFlag,
                       domain, Array<T,3>(T(),T(),T()));
    setToConstant<T,3>(*vorticity, voxelizedDomain.getVoxelMatrix(), solidBorderFlag,
                       domain, Array<T,3>(T(),T(),T()));
    return vorticity;
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiTensorField3D<T,3> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVorticity()
{
    return computeVorticity(lattice.getBoundingBox());
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocityNorm(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> > velocityNorm(plb::computeVelocityNorm(lattice,domain));
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant(*velocityNorm, voxelizedDomain.getVoxelMatrix(),
                  solidFlag, domain, (T)0);
    setToConstant(*velocityNorm, voxelizedDomain.getVoxelMatrix(),
                  solidBorderFlag, domain, (T)0);
    return velocityNorm;
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocityNorm()
{
    return computeVelocityNorm(lattice.getBoundingBox());
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocityComponent (
            Box3D domain, plint iComp )
{
    std::auto_ptr<MultiScalarField3D<T> > velocityComponent (
            plb::computeVelocityComponent(lattice,domain, iComp));
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant(*velocityComponent, voxelizedDomain.getVoxelMatrix(),
                  solidFlag, domain, (T)0);
    setToConstant(*velocityComponent, voxelizedDomain.getVoxelMatrix(),
                  solidBorderFlag, domain, (T)0);
    return velocityComponent;
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeVelocityComponent(plint iComp)
{
    return computeVelocityComponent(lattice.getBoundingBox(), iComp);
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computePressure(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> > pressure(plb::computeDensity(lattice,domain));
    T averageDensity = computeAverageDensity(domain);
    subtractInPlace(*pressure, averageDensity, domain);
    multiplyInPlace(*pressure, Descriptor<T>::cs2, domain);
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant(*pressure, voxelizedDomain.getVoxelMatrix(),
                  solidFlag, domain, (T)0);
    setToConstant(*pressure, voxelizedDomain.getVoxelMatrix(),
                  solidBorderFlag, domain, (T)0);
    return pressure;
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computePressure()
{
    return computePressure(lattice.getBoundingBox());
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeDensity(Box3D domain, T solidDensity)
{
    std::auto_ptr<MultiScalarField3D<T> > density(plb::computeDensity(lattice,domain));
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant(*density, voxelizedDomain.getVoxelMatrix(),
                  solidFlag, domain, solidDensity);
    setToConstant(*density, voxelizedDomain.getVoxelMatrix(),
                  solidBorderFlag, domain, solidDensity);
    return density;
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeDensity(T solidDensity)
{
    return computeDensity(lattice.getBoundingBox(), solidDensity);
}

    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeStrainRateNorm(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> >
         strainRateNorm(computeSymmetricTensorNorm(*plb::computeStrainRateFromStress(lattice,domain)));
    int flowType = voxelizedDomain.getFlowType();
    int solidFlag = voxelFlag::invert(flowType);
    int solidBorderFlag = voxelFlag::borderFlag(solidFlag);
    setToConstant(*strainRateNorm, voxelizedDomain.getVoxelMatrix(),
                  solidFlag, domain, (T)0);
    setToConstant(*strainRateNorm, voxelizedDomain.getVoxelMatrix(),
                  solidBorderFlag, domain, (T)0);
    return strainRateNorm;
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
std::auto_ptr<MultiScalarField3D<T> >
    OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeStrainRateNorm()
{
    return computeStrainRateNorm(lattice.getBoundingBox());
}


template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeAverageVelocityComponent(Box3D domain, plint iComponent)
{
    std::auto_ptr<MultiScalarField3D<T> > density (
            plb::computeVelocityComponent(lattice,domain, iComponent) );
    MultiScalarField3D<int> flagMatrix((MultiBlock3D&)voxelizedDomain.getVoxelMatrix());
    int flowType = voxelizedDomain.getFlowType();
    int fluidBorderFlag = voxelFlag::borderFlag(flowType);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  flowType, domain, 1);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  fluidBorderFlag, domain, 1);
    return computeAverage(*density, flagMatrix, 1, domain);
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeAverageDensity()
{
    return computeAverageDensity(voxelizedDomain.getVoxelMatrix().getBoundingBox());
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeAverageDensity(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> > density (
            plb::computeDensity(lattice,domain) );
    std::auto_ptr<MultiScalarField3D<T> > density2 (
            plb::computeDensity(lattice,domain) );
    MultiScalarField3D<int> flagMatrix((MultiBlock3D&)voxelizedDomain.getVoxelMatrix());
    int flowType = voxelizedDomain.getFlowType();
    int fluidBorderFlag = voxelFlag::borderFlag(flowType);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  flowType, domain, 1);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  fluidBorderFlag, domain, 1);
    return computeAverage(*density, flagMatrix, 1, domain);
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeAverageEnergy()
{
    return computeAverageEnergy(voxelizedDomain.getVoxelMatrix().getBoundingBox());
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeAverageEnergy(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> > energy (
            plb::computeKineticEnergy(lattice,domain) );
    MultiScalarField3D<int> flagMatrix((MultiBlock3D&)voxelizedDomain.getVoxelMatrix());
    int flowType = voxelizedDomain.getFlowType();
    int fluidBorderFlag = voxelFlag::borderFlag(flowType);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  flowType, domain, 1);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  fluidBorderFlag, domain, 1);
    return computeAverage(*energy, flagMatrix, 1, domain);
}
    
template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeRMSvorticity()
{
    return computeRMSvorticity(voxelizedDomain.getVoxelMatrix().getBoundingBox());
}

template< typename T,
          template<typename U> class Descriptor,
          class BoundaryType >
T OffLatticeBoundaryCondition3D<T,Descriptor,BoundaryType>::computeRMSvorticity(Box3D domain)
{
    std::auto_ptr<MultiScalarField3D<T> > vorticityNormSqr (
            plb::computeNormSqr(*plb::computeBulkVorticity(*plb::computeVelocity(lattice,domain)) ) );
    MultiScalarField3D<int> flagMatrix((MultiBlock3D&)voxelizedDomain.getVoxelMatrix());
    int flowType = voxelizedDomain.getFlowType();
    int fluidBorderFlag = voxelFlag::borderFlag(flowType);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  flowType, domain, 1);
    setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
                  fluidBorderFlag, domain, 1);
    return sqrt(computeAverage(*vorticityNormSqr, flagMatrix, 1, domain));
}
    
}  // namespace plb

#endif  // OFF_LATTICE_BOUNDARY_CONDITION_3D_HH