boundaryTemplates.h

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 BOUNDARY_TEMPLATES_H
#define BOUNDARY_TEMPLATES_H

#include "boundaryCondition/regularizedBoundaryDynamics.h"
#include "core/cell.h"
#include "latticeBoltzmann/indexTemplates.h"

namespace plb {

template<typename T, template<typename U> class Descriptor, int direction, int orientation>
struct boundaryTemplates {
    static void compute_PiNeq (
         Dynamics<T,Descriptor> const& dynamics,
         Cell<T,Descriptor> const& cell, T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr,
         Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq )
    {
        typedef Descriptor<T> L;

        std::vector<plint> const& onWallIndices = indexTemplates::subIndex<L, direction, 0>();
        std::vector<plint> const& normalIndices = indexTemplates::subIndex<L, direction, orientation>();

        // Compute off-equilibrium for known particle populations.
        Array<T,Descriptor<T>::q> fNeq;
        for (pluint fIndex=0; fIndex<onWallIndices.size(); ++fIndex) {
            plint iPop = onWallIndices[fIndex];
            fNeq[iPop] = cell[iPop] - dynamics.computeEquilibrium(iPop, rhoBar, j, jSqr);
        }
        for (pluint fIndex=0; fIndex<normalIndices.size(); ++fIndex) {
            plint iPop = normalIndices[fIndex];
            if (iPop == 0) {
                fNeq[iPop] = T();  // fNeq[0] will not be used anyway
            }
            else {
                fNeq[iPop] = cell[iPop] - dynamics.computeEquilibrium(iPop, rhoBar, j, jSqr);
            }
        }

        // Compute PiNeq from fNeq, by using "bounce-back of off-equilibrium part" rule.
        int iPi = 0;
        for (int iAlpha=0; iAlpha<L::d; ++iAlpha) {
            for (int iBeta=iAlpha; iBeta<L::d; ++iBeta) {
                PiNeq[iPi] = T();
                for (pluint fIndex=0; fIndex<onWallIndices.size(); ++fIndex)
                {
                    const plint iPop = onWallIndices[fIndex];
                    PiNeq[iPi] += L::c[iPop][iAlpha]*L::c[iPop][iBeta]*fNeq[iPop];
                }
                for (pluint fIndex=0; fIndex<normalIndices.size(); ++fIndex)
                {
                    const plint iPop = normalIndices[fIndex];
                    PiNeq[iPi] += (T)2 * L::c[iPop][iAlpha]*L::c[iPop][iBeta]* fNeq[iPop];
                }
                ++iPi;
            }
        }
    }

};  // struct boundaryTemplates

}  // namespace plb

#endif  // BOUNDARY_TEMPLATES_H