C++ volVectorField::correctBoundaryConditions方法代码示例

void averagingModel::setVectorSumSimple
(
    volVectorField& field,
    double**& value,
    double**& weight,
    int nP
) const
{
    label cellI;
    label subCell=0;
    vector valueVec;
    scalar weightP;

    for(int index=0; index< nP; index++)
    {
        cellI = particleCloud_.cellIDs()[index][subCell];

        if (cellI >= 0)
        {
            for(int i=0;i<3;i++) valueVec[i] = value[index][i];
            weightP = weight[index][subCell];
            field[cellI] += valueVec*weightP;
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

void averagingModel::undoVectorSum
(
    volVectorField& field,
    double** const& value,
    double** const& weight,
    double**const& mask
) const
{
    label cellI;
    vector valueVec;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                //Info << "subCell=" << subCell << endl;
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                    weightP = weight[index][subCell];

                    field[cellI] -= valueVec*weightP;
                }
            }//forAllSubPoints
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

void dilute::setVectorAverage
(
    volVectorField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask,
    double**const& weight2,         //allows the specification of a 2nd weight field
    bool weightWithWeight2    //switch to activate 2nd weight field
) const
{
    label cellI;
    vector valueVec;
    scalar weightP;

    if(weightWithWeight2) //use weight2, e.g., mass-averaged - has no effect, just weight is DIFFERENT!
    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {
            cellI = particleCloud_.cellIDs()[index][subCell];
            if (cellI >= 0)
            {
                for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                weightP = weight[index][subCell]*weight2[index][subCell];
                weightField[cellI] += weightP;
                if(weightP > 0) field[cellI] = valueVec; //field[cellI] = valueVec/weightP;
            }
        }
    }
    else //standard, i.e., volume-averaged - has no effect, just weight is DIFFERENT!
    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {
            //Info << "subCell=" << subCell << endl;
            cellI = particleCloud_.cellIDs()[index][subCell];

            if (cellI >= 0)
            {
                for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                weightP = weight[index][subCell];
                weightField[cellI] += weightP;
                if(weightP > 0) field[cellI] = valueVec; //field[cellI] = valueVec/weightP;
                else Warning << "!!! W A R N I N G --- weightP <= 0" << endl;
            }
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
void averagingModel::undoVectorAverage
(
    volVectorField& fieldPrev,
    volVectorField& fieldNext,
    volScalarField& weightField,
    double** const& value,
    double** const& weight,
    double**const& mask,
    bool single
) const
{
// WARNING - not sure if this is valid for dilute model!!!

    if(!single) fieldPrev.internalField() = fieldNext.internalField();

    label cellI;
    vector valueVec;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                //Info << "subCell=" << subCell << endl;
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                    weightP = weight[index][subCell];

                    if(weightField[cellI] == weightP)
                    {
                        fieldNext[cellI] = vector(0,0,0);
                    }else
                    {
                        fieldNext[cellI] = (fieldNext[cellI]*weightField[cellI]-valueVec*weightP)/(weightField[cellI]-weightP);
                    }
                }
            }
        }
    }

    // correct cell values to patches
    fieldNext.correctBoundaryConditions();
}

void dense::setVectorAverage
(
    volVectorField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask
) const
{
    label cellI;
    vector valueVec;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                    weightP = weight[index][subCell];

                    // first entry in this cell
                    if(weightField[cellI] == 0)
                    {
                        field[cellI] = valueVec;
                        weightField[cellI] = weightP;
                    }
                    else
                    {
                        field[cellI] = (field[cellI]*weightField[cellI]+valueVec*weightP)/(weightField[cellI]+weightP);
                        weightField[cellI] += weightP;
                    }
                }
            }//forAllSubPoints
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

//.........这里部分代码省略.........
            clockM().start(19,"setvoidFraction");
            if(verbose_) Info << "- setvoidFraction()" << endl;
            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
            if(verbose_) Info << "setvoidFraction done." << endl;
            clockM().stop("setvoidFraction");

            // set particles velocity field
            clockM().start(20,"setVectorAverage");
            setVectorAverages();
            clockM().stop("setVectorAverage");

            // set particles forces
            clockM().start(21,"setForce");
            if(verbose_) Info << "- setForce(forces_)" << endl;
            setForces();
            if(verbose_) Info << "setForce done." << endl;
            clockM().stop("setForce");

            // get next force field
            clockM().start(22,"setParticleForceField");
            if(verbose_) Info << "- setParticleForceField()" << endl;
            averagingM().setVectorSum
            (
                forceM(0).impParticleForces(),
                impForces_,
                particleWeights_,
                NULL //mask
            );
            averagingM().setVectorSum
            (
                forceM(0).expParticleForces(),
                expForces_,
                particleWeights_,
                NULL //mask
            );
            if(verbose_) Info << "- setParticleForceField done." << endl;
            clockM().stop("setParticleForceField");

            // write DEM data
            if(verbose_) Info << " -giveDEMdata()" << endl;
            clockM().start(23,"giveDEMdata");
            giveDEMdata();
            clockM().stop("giveDEMdata");
        }//end dataExchangeM().couple()
        Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;

        clockM().start(24,"interpolateEulerFields");
        // update smoothing model
        smoothingM().dSmoothing();

        //============================================
        // update voidFractionField V1
        alpha = voidFractionM().voidFractionInterp();
        smoothingM().smoothen(alpha);
        if(dataExchangeM().couplingStep() < 2)
        {
            alpha.oldTime() = alpha; // supress volume src
            alpha.oldTime().correctBoundaryConditions();
        }
        alpha.correctBoundaryConditions();

        // calc ddt(voidfraction)
        //calcDdtVoidfraction(voidFractionM().voidFractionNext());
        calcDdtVoidfraction(alpha);

        // update particle velocity Field
        Us = averagingM().UsInterp();
        //smoothingM().smoothenReferenceField(Us);
        Us.correctBoundaryConditions();
        /*//============================================
        // update voidFractionField
        volScalarField oldAlpha = alpha.oldTime(); //save old (smooth) alpha field
        alpha.oldTime().internalField() = voidFractionM().voidFractionInterp();
        smoothingM().smoothen(alpha);
        alpha.correctBoundaryConditions();
        alpha.oldTime() = oldAlpha; //set old (smooth) alpha field to allow correct computation of ddt

        // calc ddt(voidfraction)
        if (doCouple) calcDdtVoidfraction(alpha);
        //calcDdtVoidfraction(alpha); // alternative with scale=1! (does not see change in alpha?)

        // update particle velocity Field
        Us.oldTime().internalField() = averagingM().UsInterp();
        smoothingM().smoothenReferenceField(Us);
        Us.correctBoundaryConditions();
        //============================================*/
        clockM().stop("interpolateEulerFields");

        if(verbose_){
            #include "debugInfo.H"
        }

        clockM().start(25,"dumpDEMdata");
        // do particle IO
        IOM().dumpDEMdata();
        clockM().stop("dumpDEMdata");

    }//end ignore
    return doCouple;
}

bool Foam::cfdemCloud::evolve
(
    volScalarField& alpha,
    volVectorField& Us,
    volVectorField& U
)
{
    numberOfParticlesChanged_ = false;
    arraysReallocated_=false;
    bool doCouple=false;

    if(!ignore())
    {
        if (dataExchangeM().doCoupleNow())
        {
            Info << "\n Coupling..." << endl;
            dataExchangeM().couple(0);
            doCouple=true;

            // reset vol Fields
            clockM().start(16,"resetVolFields");
            if(verbose_)
            {
                Info << "couplingStep:" << dataExchangeM().couplingStep() 
                     << "\n- resetVolFields()" << endl;
            }
            averagingM().resetVectorAverage(averagingM().UsPrev(),averagingM().UsNext(),false);
            resetVoidFraction();
            averagingM().resetVectorAverage(forceM(0).impParticleForces(),forceM(0).impParticleForces(),true);
            averagingM().resetVectorAverage(forceM(0).expParticleForces(),forceM(0).expParticleForces(),true);
            averagingM().resetWeightFields();
            for (int i=0;i<momCoupleModels_.size(); i++)
                momCoupleM(i).resetMomSourceField();
            if(verbose_) Info << "resetVolFields done." << endl;
            clockM().stop("resetVolFields");

            if(verbose_) Info << "- getDEMdata()" << endl;
            clockM().start(17,"getDEMdata");
            getDEMdata();
            clockM().stop("getDEMdata");
            if(verbose_) Info << "- getDEMdata done." << endl;

            // search cellID of particles
            clockM().start(18,"findCell");
            if(verbose_) Info << "- findCell()" << endl;
            findCells();
            if(verbose_) Info << "findCell done." << endl;
            clockM().stop("findCell");

            // set void fraction field
            clockM().start(19,"setvoidFraction");
            if(verbose_) Info << "- setvoidFraction()" << endl;
            setVoidFraction();
            if(verbose_) Info << "setvoidFraction done." << endl;
            clockM().stop("setvoidFraction");

            // set average particles velocity field
            clockM().start(20,"setVectorAverage");
            setVectorAverages();


            //Smoothen "next" fields            
            smoothingM().dSmoothing();
            smoothingM().smoothen(voidFractionM().voidFractionNext());

            //only smoothen if we use implicit force coupling in cells void of particles
            //because we need unsmoothened Us field to detect cells for explicit 
            //force coupling
            if(!treatVoidCellsAsExplicitForce())
                smoothingM().smoothenReferenceField(averagingM().UsNext());
            
            clockM().stop("setVectorAverage");
        }
        
        //============================================
        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD 
        //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
        //      QUANTITIES AT THE GRID!
        Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
        clockM().start(24,"interpolateEulerFields");

        // update voidFractionField
        setAlpha(alpha);
        if(dataExchangeM().couplingStep() < 2)
        {
            alpha.oldTime() = alpha; // supress volume src
            alpha.oldTime().correctBoundaryConditions();
        }
        alpha.correctBoundaryConditions();

        // calc ddt(voidfraction)
        calcDdtVoidfraction(alpha,Us);

        // update mean particle velocity Field
        Us = averagingM().UsInterp();
        Us.correctBoundaryConditions();

        clockM().stop("interpolateEulerFields");
        //============================================

//.........这里部分代码省略.........