path: root/src/boundary/boundaryPostProcessors3D.hh

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006, 2007 Jonas Latt
 *  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 BOUNDARY_POST_PROCESSORS_3D_HH
#define BOUNDARY_POST_PROCESSORS_3D_HH

#include "boundaryPostProcessors3D.h"
#include "core/finiteDifference3D.h"
#include "core/blockLattice3D.h"
#include "dynamics/firstOrderLbHelpers.h"
#include "core/util.h"
#include "utilities/vectorHelpers.h"

namespace olb {

////////  PlaneFdBoundaryProcessor3D ///////////////////////////////////

template<typename T, typename DESCRIPTOR, int direction, int orientation>
PlaneFdBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::
PlaneFdBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
  : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
{
  OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
}

template<typename T, typename DESCRIPTOR, int direction, int orientation>
void PlaneFdBoundaryProcessor3D<T,DESCRIPTOR,direction,orientation>::
processSubDomain(BlockLattice3D<T,DESCRIPTOR>& blockLattice,
                 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
  using namespace olb::util::tensorIndices3D;

  int newX0, newX1, newY0, newY1, newZ0, newZ1;
  if ( util::intersect (
         x0, x1, y0, y1, z0, z1,
         x0_, x1_, y0_, y1_, z0_, z1_,
         newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
    int iX;

#ifdef PARALLEL_MODE_OMP
    #pragma omp parallel for
#endif
    for (iX=newX0; iX<=newX1; ++iX) {
      T dx_u[DESCRIPTOR::d], dy_u[DESCRIPTOR::d], dz_u[DESCRIPTOR::d];
      for (int iY=newY0; iY<=newY1; ++iY) {
        for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
          Cell<T,DESCRIPTOR>& cell = blockLattice.get(iX,iY,iZ);
          Dynamics<T,DESCRIPTOR>* dynamics = blockLattice.getDynamics(iX, iY, iZ);
          T rho, u[DESCRIPTOR::d];
          cell.computeRhoU(rho,u);

          interpolateGradients<0> ( blockLattice, dx_u, iX, iY, iZ );
          interpolateGradients<1> ( blockLattice, dy_u, iX, iY, iZ );
          interpolateGradients<2> ( blockLattice, dz_u, iX, iY, iZ );
          T dx_ux = dx_u[0];
          T dy_ux = dy_u[0];
          T dz_ux = dz_u[0];
          T dx_uy = dx_u[1];
          T dy_uy = dy_u[1];
          T dz_uy = dz_u[1];
          T dx_uz = dx_u[2];
          T dy_uz = dy_u[2];
          T dz_uz = dz_u[2];
          T omega = dynamics->getOmega();
          T sToPi = - rho / descriptors::invCs2<T,DESCRIPTOR>() / omega;
          T pi[util::TensorVal<DESCRIPTOR >::n];
          pi[xx] = (T)2 * dx_ux * sToPi;
          pi[yy] = (T)2 * dy_uy * sToPi;
          pi[zz] = (T)2 * dz_uz * sToPi;
          pi[xy] = (dx_uy + dy_ux) * sToPi;
          pi[xz] = (dx_uz + dz_ux) * sToPi;
          pi[yz] = (dy_uz + dz_uy) * sToPi;

          // Computa