BoundaryFunctions.cpp

Go to the documentation of this file.

// ==============================================================================================
// Microflow 3D, http://www.microflow.pwr.edu.pl/
// Created by Roman Szafran on 18.05.19.
// Modified by Roman Szafran, last on 07.09.19.
// Copyright (c) 2019 Wroclaw University of Science and Technology.
// Distributed under the Apache License, Version 2.0. You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0 or see accompanying file license.txt.
// Redistributions of source code must retain the above copyright and license notice.
// ==============================================================================================

#include "BoundaryFunctions.h"
#define f(x) (pNode->FQ19[x])
#define XVelocityLB_ux (Vector4.x)
#define YVelocityLB_uy (Vector4.y)
#define ZVelocityLB_uz (Vector4.z)
#define DensityLB_rho (Vector4.rho)
#define BoundaryDensityLB_rhoB (Thread_Ptr->m_BoundaryValue_RhoLB)
#define BoundaryXVelocityLB_uxB (Thread_Ptr->m_BoundaryValue_ULB.x)
#define BoundaryYVelocityLB_uyB (Thread_Ptr->m_BoundaryValue_ULB.y)
#define BoundaryZVelocityLB_uzB (Thread_Ptr->m_BoundaryValue_ULB.z)
#define ReverseDirectionVector_Rc(x) (MF::Solver_CPU::BoundaryFunctions::m_LatticeParameters_Ptr->ReverseDirectionVector_Rc[x])

void MF::Solver_CPU::BoundaryFunctions::SetBoundaryNodePointerToFunc(const std::shared_ptr<MF::Database::ThreadArray>& ThreadArray_Ptr) {
    for (auto & Thread_Ptr : *ThreadArray_Ptr->m_ThreadsTable_Ptr) {
        if (MF::Solver_CPU::BoundaryFunctions::m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::QuasiCompressible) { //Quasi compressible fluid flow model
            switch (Thread_Ptr->m_NodeType) {
                default:
                case 0: {break;} //solid node
                case 1: {break;} //fluid node
                    //------------------------------------------------------------------------------------
                    //Periodic
                    //------------------------------------------------------------------------------------
                case 4: {break;}
                    //------------------------------------------------------------------------------------
                    //Bounce Back 2
                    //------------------------------------------------------------------------------------
                case 61: {  //BB typ 2 full way standard
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_BounceBack2_61;
                    Thread_Ptr->m_DoPreCollision = false;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_BounceBack2_61;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Dirichlet Velocity
                    //------------------------------------------------------------------------------------
                case 21:    //North velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityNorth_21;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityNorth_21;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 22:    //South velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocitySouth_22;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocitySouth_22;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 23:    //East velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityEast_23;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityEast_23;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 24:    // West velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityWest_24;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityWest_24;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 25:    // Bottom velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityBottom_25;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityBottom_25;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 26:    // Top velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityTop_26;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityTop_26;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Pressure, Neumann
                    //------------------------------------------------------------------------------------
                case 31:    //North pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureNorth_31;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureNorth_31;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 32:    //South pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureSouth_32;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureSouth_32;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 33:    //East pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureEast_33;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureEast_33;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 34:    //West pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureWest_34;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureWest_34;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 35:    // Bottom pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureBottom_35;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureBottom_35;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 36:    // Top pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFQC_PressureTop_36;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureTop_36;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Dirichlet Velocity 0
                    //------------------------------------------------------------------------------------
                case 41:    //North velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0North_41;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0North_41;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 42:    //South velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0South_42;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0South_42;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 43:    //East velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0East_43;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0East_43;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 44:    // West velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0West_44;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0West_44;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 45:    // Bottom velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0Bottom_45;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0Bottom_45;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 46:    // Top velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0Top_46;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0Top_46;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }

            }
        }
        else if (MF::Solver_CPU::BoundaryFunctions::m_CaseParameters_Ptr->FluidFlowModel_MT == MF::Solver_CPU::FlowTypeClass::Incompressible) { //Incompressible fluid flow model
            switch (Thread_Ptr->m_NodeType) {
                default:
                case 0: {break;} //solid node
                case 1: {break;} //fluid node
                    //------------------------------------------------------------------------------------
                    //Periodic
                    //------------------------------------------------------------------------------------
                case 4: {break;}
                    //------------------------------------------------------------------------------------
                    //Bounce Back 2
                    //------------------------------------------------------------------------------------
                case 61:    //BB typ 2 full way standard
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_BounceBack2_61;
                    Thread_Ptr->m_DoPreCollision = false;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_BounceBack2_61;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Dirichlet Velocity
                    //------------------------------------------------------------------------------------
                case 21:    //North velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityNorth_21;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityNorth_21;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 22:    //South velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocitySouth_22;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocitySouth_22;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 23:    //East velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityEast_23;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityEast_23;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 24:    // West velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityWest_24;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityWest_24;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 25:    // Bottom velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityBottom_25;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityBottom_25;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 26:    // Top velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityTop_26;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityTop_26;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Pressure, Neumann
                    //------------------------------------------------------------------------------------
                case 31:    //North pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureNorth_31;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureNorth_31;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 32:    //South pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureSouth_32;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureSouth_32;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 33:    //East pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureEast_33;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureEast_33;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 34:    //West pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureWest_34;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureWest_34;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 35:    // Bottom pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureBottom_35;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureBottom_35;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 36:    // Top pressure
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BFIC_PressureTop_36;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureTop_36;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                    //------------------------------------------------------------------------------------
                    //Dirichlet Velocity
                    //------------------------------------------------------------------------------------
                case 41:    //North velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0North_41;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0North_41;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 42:    //South velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0South_42;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0South_42;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 43:    //East velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0East_43;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0East_43;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 44:    // West velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0West_44;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0West_44;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 45:    // Bottom velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0Bottom_45;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0Bottom_45;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
                case 46:    // Top velocity
                {
                    Thread_Ptr->m_pBoundaryFunction = MF::Solver_CPU::BoundaryFunctions::BF_Velocity0Top_46;
                    Thread_Ptr->m_pVRLBFunction = MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0Top_46;
                    Thread_Ptr->m_DoPreCollision = true;
                    break;
                }
            }
        }
    }
}

inline double MF::Solver_CPU::BoundaryFunctions::MeanRHOAdjFluidNodes(const MF::Database::Node *pNode) { // Returns mean rho LB of values from adjacent nodes to pNode.
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    auto nodeNr = Thread_Ptr->getNodeNr(pNode);
    double rho {0};
    uint32_t nNode {0};

    for (auto f = 1; f < MFQ19; f++) {
        if (Thread_Ptr->getLinkedNodePointer(nodeNr, f) != nullptr) { // not a solid node
            if (Thread_Ptr->getLinkedNodePointer(nodeNr, f)-> pMyThread_Ptr->operator*().m_NodeType == 1 ) { //Linked node is fluid(1) node.
                nNode++;
                const MF::Database::Node * neighboringNode = Thread_Ptr->getLinkedNodePointer(nodeNr, f);
                const std::shared_ptr<MF::Database::Thread>& Thread_Ptr5 = *Thread_Ptr->getLinkedNodePointer(nodeNr, f)->pMyThread_Ptr;
                const MF::Database::Vec4<double> Vector5 = (*Thread_Ptr5->m_pVRLBFunction)(neighboringNode);
                rho += Vector5.rho;
            }
        }
    }
    if (nNode != 0)
        return (rho/nNode);
    else
        return -1;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_BounceBack2_61(const MF::Database::Node *pNode) { // Rho LB for BB 61, quasi compressible flow model
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = MF::Solver_CPU::BoundaryFunctions::MeanRHOAdjFluidNodes(pNode);

    if (DensityLB_rho == -1) {
        const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
        auto nodeNr = Thread_Ptr->getNodeNr(pNode);
        double rho {0};
        uint32_t nNode {0};
        for (auto f = 1; f < MFQ19; f++) {
            if (Thread_Ptr->getLinkedNodePointer(nodeNr, f) != nullptr) { // not a solid node
                const MF::Database::Node *neighboringNode = Thread_Ptr->getLinkedNodePointer(nodeNr, f);
                if (MF::Solver_CPU::BoundaryFunctions::MeanRHOAdjFluidNodes(neighboringNode) != -1) {
                    rho += MF::Solver_CPU::BoundaryFunctions::MeanRHOAdjFluidNodes(neighboringNode);
                    nNode++;
                }
            }
        }
        if (nNode != 0)
            DensityLB_rho = rho/nNode;
        else
            DensityLB_rho = NAN;
    }
    return Vector4;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_BounceBack2_61(const MF::Database::Node *pNode) { // Rho LB for BB 61, incompressible compressible flow model
    return MF::Solver_CPU::BoundaryFunctions::VRBFQC_BounceBack2_61(pNode); // The same formula as for quasi compressible flow.
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0North_41(const MF::Database::Node *pNode) {       // Rho LB for North Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(15) + f(16) + f(2) + f(7) + f(8)));
    return Vector4;
}

// Velocity 0, Dirichlet type, own derivation
void MF::Solver_CPU::BoundaryFunctions::BF_Velocity0North_41(MF::Database::Node * pNode){        // North Velocity0
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3));
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6));

    f(4) = f(2);
    f(9) = f(7) + N_x;
    f(10) = f(8) - N_x;
    f(17) = f(15) + N_z;
    f(18) = f(16) - N_z;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0South_42(const MF::Database::Node *pNode) {       // Rho LB for South Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(10) + f(17) + f(18) + f(4) + f(9));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BF_Velocity0South_42(MF::Database::Node * pNode){        // South Velocity0
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3));
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6));

    f(2) = f(4);
    f(7) = f(9) - N_x;
    f(8) = f(10) + N_x;
    f(15) = f(17) - N_z;
    f(16) = f(18) + N_z;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0East_43(const MF::Database::Node *pNode){        // Rho LB for East Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(1) + f(10) + f(11) + f(12) + f(7));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BF_Velocity0East_43(MF::Database::Node * pNode){        // East Velocity0
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4));
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6));

    f(3) = f(1);
    f(13) = f(11) + N_z;
    f(14) = f(12) - N_z;
    f(9) = f(7) + N_y;
    f(8) = f(10) - N_y;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0West_44(const MF::Database::Node *pNode){        // Rho LB for West Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(13) + f(14) + f(3) + f(8) + f(9));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BF_Velocity0West_44(MF::Database::Node * pNode){        // West Velocity0
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4));
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6));

    f(1) = f(3);
    f(11) = f(13) - N_z;
    f(12) = f(14) + N_z;
    f(7) = f(9) - N_y;
    f(10) = f(8) + N_y;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::BoundaryFunctions::VRBF_Velocity0Bottom_45(const MF::Database::Node *pNode) {       // Rho LB for Bottom Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(12) + f(13) + f(16) + f(17) + f(6));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BoundaryFunctions::BF_Velocity0Bottom_45(MF::Database::Node * pNode){        // Bottom Velocity0
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9));
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9));

    f(5) = f(6);
    f(11) = f(13) - N_x;
    f(18) = f(16) + N_y;
    f(15) = f(17) - N_y;
    f(14) = f(12) + N_x;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBF_Velocity0Top_46(const MF::Database::Node *pNode){        // Rho LB for Top Velocity0
    MF::Database::Vec4<double> Vector4 {0,0,0,0};
    DensityLB_rho = f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(11) + f(14) + f(15) + f(18) + f(5));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BF_Velocity0Top_46(MF::Database::Node * pNode){        // Top Velocity0
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9));
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9));

    f(6) = f(5);
    f(13) = f(11) + N_x;
    f(16) = f(18) - N_y;
    f(17) = f(15) + N_y;
    f(12) = f(14) - N_x;
}


//----------------------------------------------------------------------------------------------------------------------------
// Quasi compressible
//----------------------------------------------------------------------------------------------------------------------------

// Velocity, Dirichlet type, Hecht and Harting (2010) dla D3Q19, own derivation, (Zou,He 1997,Dirichlet boundry)

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityNorth_21(const MF::Database::Node *pNode){    // Rho LB for North Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(15) + f(16) + f(2) + f(7) + f(8))) / (1.0 + YVelocityLB_uy);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityNorth_21(MF::Database::Node * pNode){    // North Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityNorth_21(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(4) = f(2) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    f(9) = f(7) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(10) = f(8) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy) * DensityLB_rho;
    f(17) = f(15) + N_z - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(18) = f(16) - N_z + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocitySouth_22(const MF::Database::Node *pNode){    // Rho LB for South Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(10) + f(17) + f(18) + f(4) + f(9))) / (1.0 - YVelocityLB_uy);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocitySouth_22(MF::Database::Node * pNode){    // South Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocitySouth_22(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(2) = f(4) + 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    f(7) = f(9) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(8) = f(10) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(15) = f(17) - N_z + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(16) = f(18) + N_z + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityEast_23(const MF::Database::Node *pNode){    // Rho LB for East Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(1) + f(10) + f(11) + f(12) + f(7))) / (1.0 + XVelocityLB_ux);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityEast_23(MF::Database::Node * pNode){    // East Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityEast_23(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(3) = f(1) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    f(13) = f(11) + N_z - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(14) = f(12) - N_z + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(9) = f(7) + N_y - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(8) = f(10) - N_y + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityWest_24(const MF::Database::Node *pNode){    // Rho LB for West Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(13) + f(14) + f(3) + f(8) + f(9))) / (1.0 - XVelocityLB_ux);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityWest_24(MF::Database::Node * pNode){    // West Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityWest_24(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(1) = f(3) + 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    f(11) = f(13) - N_z + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(12) = f(14) + N_z + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz) * DensityLB_rho;
    f(7) = f(9) - N_y + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(10) = f(8) + N_y + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityBottom_25(const MF::Database::Node *pNode){    // Rho LB for Bottom Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(12) + f(13) + f(16) + f(17) + f(6))) / (1.0 - ZVelocityLB_uz);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityBottom_25(MF::Database::Node * pNode){    // Bottom Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityBottom_25(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;

    f(5) = f(6) + 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;
    f(11) = f(13) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(18) = f(16) + N_y + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(15) = f(17) - N_y + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(14) = f(12) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityTop_26(const MF::Database::Node *pNode){        // Rho LB for Top Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(11) + f(14) + f(15) + f(18) + f(5))) / (1.0 + ZVelocityLB_uz);
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_VelocityTop_26(MF::Database::Node * pNode){        // Top Velocity
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_VelocityTop_26(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;

    f(6) = f(5) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;
    f(13) = f(11) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(16) = f(18) - N_y + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz) * DensityLB_rho;
    f(17) = f(15) + N_y - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(12) = f(14) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz) * DensityLB_rho;
}


// Pressure, Neumann type, Hecht and Harting (2010) for D3Q19, own derivation

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureNorth_31(const MF::Database::Node *pNode){    // Velocity LB for North Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    YVelocityLB_uy = -1.0 + (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(15) + f(16) + f(2) + f(7) + f(8))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureNorth_31(MF::Database::Node * pNode){    // North Pressure
    auto Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureNorth_31(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(4) = f(2) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    f(9) = f(7) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(10) = f(8) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy) * DensityLB_rho;
    f(17) = f(15) + N_z - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(18) = f(16) - N_z + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureSouth_32(const MF::Database::Node *pNode){    // Velocity LB for South Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    YVelocityLB_uy = 1.0 - (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(10) + f(17) + f(18) + f(4) + f(9))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureSouth_32(MF::Database::Node * pNode){    // South Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureSouth_32(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(2) = f(4) + 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    f(7) = f(9) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(8) = f(10) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(15) = f(17) - N_z + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(16) = f(18) + N_z + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureEast_33(const MF::Database::Node *pNode){    // Velocity LB for East Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    XVelocityLB_ux = -1.0 + (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(1) + f(10) + f(11) + f(12) + f(7))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureEast_33(MF::Database::Node * pNode){    // East Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureEast_33(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(3) = f(1) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    f(13) = f(11) + N_z - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(14) = f(12) - N_z + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(9) = f(7) + N_y - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(8) = f(10) - N_y + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureWest_34(const MF::Database::Node *pNode){    // Velocity LB for West Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    XVelocityLB_ux = 1.0 - (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(13) + f(14) + f(3) + f(8) + f(9))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureWest_34(MF::Database::Node * pNode){    // West Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureWest_34(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;

    f(1) = f(3) + 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    f(11) = f(13) - N_z + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(12) = f(14) + N_z + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz) * DensityLB_rho;
    f(7) = f(9) - N_y + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy) * DensityLB_rho;
    f(10) = f(8) + N_y + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureBottom_35(const MF::Database::Node *pNode){    // Velocity LB for Bottom Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    ZVelocityLB_uz = 1.0 - (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(12) + f(13) + f(16) + f(17) + f(6))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureBottom_35(MF::Database::Node * pNode){    // Bottom Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureBottom_35(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;

    f(5) = f(6) + 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;
    f(11) = f(13) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(18) = f(16) + N_y + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(15) = f(17) - N_y + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(14) = f(12) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureTop_36(const MF::Database::Node *pNode){    // Velocity LB for Top Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    ZVelocityLB_uz = -1.0 + (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(11) + f(14) + f(15) + f(18) + f(5))) / DensityLB_rho;
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFQC_PressureTop_36(MF::Database::Node * pNode){        // Top pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFQC_PressureTop_36(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux * DensityLB_rho;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy * DensityLB_rho;

    f(6) = f(5) - 1.0 / 3.0 * ZVelocityLB_uz * DensityLB_rho;
    f(13) = f(11) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz) * DensityLB_rho;
    f(16) = f(18) - N_y + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz) * DensityLB_rho;
    f(17) = f(15) + N_y - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz) * DensityLB_rho;
    f(12) = f(14) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz) * DensityLB_rho;
}

//----------------------------------------------------------------------------------------------------------------------------
// Incompressible Zou et al. 1995
//----------------------------------------------------------------------------------------------------------------------------

// Velocity, Dirichlet type, Hecht and Harting (2010) for D3Q19, own derivation. (Zou,He 1997,Dirichlet boundary)

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityNorth_21(const MF::Database::Node *pNode){    // Rho for North Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (-BoundaryYVelocityLB_uyB) + (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(15) + f(16) + f(2) + f(7) + f(8)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityNorth_21(MF::Database::Node * pNode){    // North Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * BoundaryZVelocityLB_uzB;

    f(4) = f(2) - 1.0 / 3.0 * BoundaryYVelocityLB_uyB;
    f(9) = f(7) + N_x - 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
    f(10) = f(8) - N_x + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB - BoundaryYVelocityLB_uyB);
    f(17) = f(15) + N_z - 1.0 / 6.0 * (BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
    f(18) = f(16) - N_z + 1.0 / 6.0 * (-BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocitySouth_22(const MF::Database::Node *pNode){    // Rho for South Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = BoundaryYVelocityLB_uyB + (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(10) + f(17) + f(18) + f(4) + f(9)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocitySouth_22(MF::Database::Node * pNode){    // South Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * BoundaryZVelocityLB_uzB;

    f(2) = f(4) + 1.0 / 3.0 * BoundaryYVelocityLB_uyB;
    f(7) = f(9) - N_x + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
    f(8) = f(10) + N_x + 1.0 / 6.0 * (-BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
    f(15) = f(17) - N_z + 1.0 / 6.0 * (BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
    f(16) = f(18) + N_z + 1.0 / 6.0 * (BoundaryYVelocityLB_uyB - BoundaryZVelocityLB_uzB);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityEast_23(const MF::Database::Node *pNode){    // Rho for East Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (-BoundaryXVelocityLB_uxB) + (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(1) + f(10) + f(11) + f(12) + f(7)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityEast_23(MF::Database::Node * pNode){    // East Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * BoundaryYVelocityLB_uyB;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * BoundaryZVelocityLB_uzB;

    f(3) = f(1) - 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    f(13) = f(11) + N_z - 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
    f(14) = f(12) - N_z + 1.0 / 6.0 * (-BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
    f(9) = f(7) + N_y - 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
    f(8) = f(10) - N_y + 1.0 / 6.0 * (-BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityWest_24(const MF::Database::Node *pNode){    // Rho for West Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = BoundaryXVelocityLB_uxB + (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(13) + f(14) + f(3) + f(8) + f(9)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityWest_24(MF::Database::Node * pNode){    // West Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * BoundaryYVelocityLB_uyB;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * BoundaryZVelocityLB_uzB;

    f(1) = f(3) + 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    f(11) = f(13) - N_z + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
    f(12) = f(14) + N_z + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB - BoundaryZVelocityLB_uzB);
    f(7) = f(9) - N_y + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryYVelocityLB_uyB);
    f(10) = f(8) + N_y + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB - BoundaryYVelocityLB_uyB);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityBottom_25(const MF::Database::Node *pNode){    // Rho for Bottom Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = BoundaryZVelocityLB_uzB + (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(12) + f(13) + f(16) + f(17) + f(6)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityBottom_25(MF::Database::Node * pNode){    // Bottom Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * BoundaryYVelocityLB_uyB;

    f(5) = f(6) + 1.0 / 3.0 * BoundaryZVelocityLB_uzB;
    f(11) = f(13) - N_x + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
    f(18) = f(16) + N_y + 1.0 / 6.0 * (-BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
    f(15) = f(17) - N_y + 1.0 / 6.0 * (BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
    f(14) = f(12) + N_x + 1.0 / 6.0 * (-BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_VelocityTop_26(const MF::Database::Node *pNode){        // Rho for Top Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,0};
    DensityLB_rho = (-BoundaryZVelocityLB_uzB) + (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(11) + f(14) + f(15) + f(18) + f(5)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_VelocityTop_26(MF::Database::Node * pNode){        // Top Velocity
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * BoundaryXVelocityLB_uxB;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * BoundaryYVelocityLB_uyB;

    f(6) = f(5) - 1.0 / 3.0 * BoundaryZVelocityLB_uzB;
    f(13) = f(11) + N_x - 1.0 / 6.0 * (BoundaryXVelocityLB_uxB + BoundaryZVelocityLB_uzB);
    f(16) = f(18) - N_y + 1.0 / 6.0 * (BoundaryYVelocityLB_uyB - BoundaryZVelocityLB_uzB);
    f(17) = f(15) + N_y - 1.0 / 6.0 * (BoundaryYVelocityLB_uyB + BoundaryZVelocityLB_uzB);
    f(12) = f(14) - N_x + 1.0 / 6.0 * (BoundaryXVelocityLB_uxB - BoundaryZVelocityLB_uzB);
}


// Pressure, Neumann type, Hecht and Harting (2010) for D3Q19, own derivation

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureNorth_31(const MF::Database::Node *pNode){    // Velocity LB for North Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    YVelocityLB_uy = (-DensityLB_rho) + (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(15) + f(16) + f(2) + f(7) + f(8)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureNorth_31(MF::Database::Node * pNode){    // North Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureNorth_31(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz;

    f(4) = f(2) - 1.0 / 3.0 * YVelocityLB_uy;
    f(9) = f(7) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy);
    f(10) = f(8) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy);
    f(17) = f(15) + N_z - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz);
    f(18) = f(16) - N_z + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureSouth_32(const MF::Database::Node *pNode){    // Velocity LB for South Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    YVelocityLB_uy = DensityLB_rho - (f(0) + f(1) + f(11) + f(12) + f(13) + f(14) + f(3) + f(5) + f(6) + 2.0 * (f(10) + f(17) + f(18) + f(4) + f(9)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureSouth_32(MF::Database::Node * pNode){    // South Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureSouth_32(pNode);
    const double N_x = 0.5 * (f(1) + f(11) + f(12) - f(13) - f(14) - f(3)) - 1.0 / 3.0 * XVelocityLB_ux;
    const double N_z = 0.5 * (f(11) - f(12) - f(13) + f(14) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz;

    f(2) = f(4) + 1.0 / 3.0 * YVelocityLB_uy;
    f(7) = f(9) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy);
    f(8) = f(10) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy);
    f(15) = f(17) - N_z + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz);
    f(16) = f(18) + N_z + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureEast_33(const MF::Database::Node *pNode){    // Velocity LB for East Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    XVelocityLB_ux = (-DensityLB_rho) + (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(1) + f(10) + f(11) + f(12) + f(7)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureEast_33(MF::Database::Node * pNode){    // East Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureEast_33(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz;

    f(3) = f(1) - 1.0 / 3.0 * XVelocityLB_ux;
    f(13) = f(11) + N_z - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz);
    f(14) = f(12) - N_z + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz);
    f(9) = f(7) + N_y - 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy);
    f(8) = f(10) - N_y + 1.0 / 6.0 * (-XVelocityLB_ux + YVelocityLB_uy);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureWest_34(const MF::Database::Node *pNode){    // Velocity LB for WestPressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    XVelocityLB_ux = DensityLB_rho - (f(0) + f(15) + f(16) + f(17) + f(18) + f(2) + f(4) + f(5) + f(6) + 2.0 * (f(13) + f(14) + f(3) + f(8) + f(9)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureWest_34(MF::Database::Node * pNode){    // West Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureWest_34(pNode);
    const double N_y = 0.5 * (f(15) + f(16) - f(17) - f(18) + f(2) - f(4)) - 1.0 / 3.0 * YVelocityLB_uy;
    const double N_z = 0.5 * (f(15) - f(16) - f(17) + f(18) + f(5) - f(6)) - 1.0 / 3.0 * ZVelocityLB_uz;

    f(1) = f(3) + 1.0 / 3.0 * XVelocityLB_ux;
    f(11) = f(13) - N_z + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz);
    f(12) = f(14) + N_z + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz);
    f(7) = f(9) - N_y + 1.0 / 6.0 * (XVelocityLB_ux + YVelocityLB_uy);
    f(10) = f(8) + N_y + 1.0 / 6.0 * (XVelocityLB_ux - YVelocityLB_uy);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureBottom_35(const MF::Database::Node *pNode){    // Velocity LB for Bottom Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    ZVelocityLB_uz = DensityLB_rho - (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(12) + f(13) + f(16) + f(17) + f(6)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureBottom_35(MF::Database::Node * pNode){    // Bottom Pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureBottom_35(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy;

    f(5) = f(6) + 1.0 / 3.0 * ZVelocityLB_uz;
    f(11) = f(13) - N_x + 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz);
    f(18) = f(16) + N_y + 1.0 / 6.0 * (-YVelocityLB_uy + ZVelocityLB_uz);
    f(15) = f(17) - N_y + 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz);
    f(14) = f(12) + N_x + 1.0 / 6.0 * (-XVelocityLB_ux + ZVelocityLB_uz);
}

inline MF::Database::Vec4<double> MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureTop_36(const MF::Database::Node *pNode){    // Velocity LB for Top Pressure
    const std::shared_ptr<MF::Database::Thread>& Thread_Ptr = *pNode->pMyThread_Ptr;
    MF::Database::Vec4<double> Vector4 = {BoundaryXVelocityLB_uxB,BoundaryYVelocityLB_uyB,BoundaryZVelocityLB_uzB,BoundaryDensityLB_rhoB};
    ZVelocityLB_uz = -DensityLB_rho + (f(0) + f(1) + f(10) + f(2) + f(3) + f(4) + f(7) + f(8) + f(9) + 2.0 * (f(11) + f(14) + f(15) + f(18) + f(5)));
    return Vector4;
}

void MF::Solver_CPU::BoundaryFunctions::BFIC_PressureTop_36(MF::Database::Node * pNode){        // Top pressure
    MF::Database::Vec4<double> Vector4 = MF::Solver_CPU::BoundaryFunctions::VRBFIC_PressureTop_36(pNode);
    const double N_x = 0.5 * (f(1) + f(10) - f(3) + f(7) - f(8) - f(9)) - 1.0 / 3.0 * XVelocityLB_ux;
    const double N_y = 0.5 * (-f(10) + f(2) - f(4) + f(7) + f(8) - f(9)) - 1.0 / 3.0 * YVelocityLB_uy;

    f(6) = f(5) - 1.0 / 3.0 * ZVelocityLB_uz;
    f(13) = f(11) + N_x - 1.0 / 6.0 * (XVelocityLB_ux + ZVelocityLB_uz);
    f(16) = f(18) - N_y + 1.0 / 6.0 * (YVelocityLB_uy - ZVelocityLB_uz);
    f(17) = f(15) + N_y - 1.0 / 6.0 * (YVelocityLB_uy + ZVelocityLB_uz);
    f(12) = f(14) - N_x + 1.0 / 6.0 * (XVelocityLB_ux - ZVelocityLB_uz);
}