void Foam::processorFvPatchField<Type>::initInterfaceMatrixUpdate
(
    gpuField<Type>&,
    const gpuField<Type>& psiInternal,
    const scalargpuField&,
    const Pstream::commsTypes commsType
) const
{
    this->patch().patchInternalField(psiInternal, gpuSendBuf_);

    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
    {
        // Fast path.
        if (debug && !this->ready())
        {
            FatalErrorIn
            (
                "processorFvPatchField<Type>::initInterfaceMatrixUpdate(..)"
            )   << "On patch " << procPatch_.name()
                << " outstanding request."
                << abort(FatalError);
        }

        std::streamsize nBytes = gpuSendBuf_.byteSize();

        Type* receive;
        const Type* send;

        if(Pstream::gpuDirectTransfer)
        {
            // Fast path.
            gpuReceiveBuf_.setSize(gpuSendBuf_.size());

            send = gpuSendBuf_.data();
            receive = gpuReceiveBuf_.data();
        }
        else
        {
            sendBuf_.setSize(gpuSendBuf_.size());
            receiveBuf_.setSize(sendBuf_.size());
            thrust::copy
            (
                gpuSendBuf_.begin(),
                gpuSendBuf_.end(),
                sendBuf_.begin()
            );

            send = sendBuf_.begin();
            receive = receiveBuf_.begin();
        }

        outstandingRecvRequest_ = UPstream::nRequests();
        IPstream::read
        (
            Pstream::nonBlocking,
            procPatch_.neighbProcNo(),
            reinterpret_cast<char*>(receive),
            nBytes,
            procPatch_.tag(),
            procPatch_.comm()
        );

        outstandingSendRequest_ = UPstream::nRequests();
        OPstream::write
        (
            Pstream::nonBlocking,
            procPatch_.neighbProcNo(),
            reinterpret_cast<const char*>(send),
            nBytes,
            procPatch_.tag(),
            procPatch_.comm()
        );
    }
    else
    {
        procPatch_.compressedSend(commsType, gpuSendBuf_);
    }

    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = false;
}

void Foam::sixDofTopoMotion::addZonesAndModifiers()
{
    // Add zones and modifiers for motion action

    if (useTopoSliding_)
    {
        if
        (
            pointZones().size() > 0
         || faceZones().size() > 0
         || cellZones().size() > 0
        )
        {
            Info<< "void sixDofTopoMotion::addZonesAndModifiers() : "
                << "Zones and modifiers already present.  Skipping."
                << endl;

            if (topoChanger_.size() == 0)
            {
                FatalErrorIn
                (
                    "void sixDofTopoMotion::addZonesAndModifiers()"
                )   << "Mesh modifiers not read properly"
                    << abort(FatalError);
            }

            return;
        }

        Info<< "Time = " << time().timeName() << endl
            << "Adding zones and modifiers to the mesh" << endl;

        // Add zones
        List<pointZone*> pz(3*bodies_.size());
        List<faceZone*> fz(3*bodies_.size());
        List<cellZone*> cz(0);

        label npz = 0;
        label nfz = 0;
        label nSliders = 0;

        forAll (bodies_, bodyI)
        {
            const floatingBody& curBody = bodies_[bodyI];

            if
            (
                curBody.hullSlider().active()
             && curBody.fixedSlider().active()
            )
            {
                nSliders++;

                // Add an empty zone for cut points
                pz[npz] = new pointZone
                (
                    curBody.name() + "CutPointZone",
                    labelList(0),
                    npz,
                    pointZones()
                );
                npz++;

                // Do face zones for slider

                // Inner slider
                const polyPatch& innerSlider =
                    boundaryMesh()[curBody.hullSlider().index()];

                labelList isf(innerSlider.size());

                forAll (isf, i)
                {
                    isf[i] = innerSlider.start() + i;
                }

                fz[nfz] = new faceZone
                (
                    curBody.name() + "InsideSliderZone",
                    isf,
                    boolList(innerSlider.size(), false),
                    nfz,
                    faceZones()
                );
                nfz++;

                // Outer slider
                const polyPatch& outerSlider =
                    boundaryMesh()[curBody.fixedSlider().index()];

                labelList osf(outerSlider.size());

                forAll (osf, i)
                {
                    osf[i] = outerSlider.start() + i;
                }

                fz[nfz] = new faceZone
                (
                    curBody.name() + "OutsideSliderZone",
//.........这里部分代码省略.........

// Call Metis with options from dictionary.
Foam::label Foam::metisDecomp::decompose
(
    const List<int>& adjncy,
    const List<int>& xadj,
    const scalarField& cWeights,

    List<int>& finalDecomp
)
{
    // C style numbering
    int numFlag = 0;

    // Method of decomposition
    // recursive: multi-level recursive bisection (default)
    // k-way: multi-level k-way
    word method("k-way");

    int numCells = xadj.size()-1;

    // decomposition options. 0 = use defaults
    List<int> options(5, 0);

    // processor weights initialised with no size, only used if specified in
    // a file
    Field<floatScalar> processorWeights;

    // cell weights (so on the vertices of the dual)
    List<int> cellWeights;

    // face weights (so on the edges of the dual)
    List<int> faceWeights;


    // Check for externally provided cellweights and if so initialise weights
    scalar minWeights = gMin(cWeights);
    if (cWeights.size() > 0)
    {
        if (minWeights <= 0)
        {
            WarningIn
            (
                "metisDecomp::decompose"
                "(const pointField&, const scalarField&)"
            )   << "Illegal minimum weight " << minWeights
                << endl;
        }

        if (cWeights.size() != numCells)
        {
            FatalErrorIn
            (
                "metisDecomp::decompose"
                "(const pointField&, const scalarField&)"
            )   << "Number of cell weights " << cWeights.size()
                << " does not equal number of cells " << numCells
                << exit(FatalError);
        }
        // Convert to integers.
        cellWeights.setSize(cWeights.size());
        forAll(cellWeights, i)
        {
            cellWeights[i] = int(cWeights[i]/minWeights);
        }
    }


    // Check for user supplied weights and decomp options
    if (decompositionDict_.found("metisCoeffs"))
    {
        const dictionary& metisCoeffs =
            decompositionDict_.subDict("metisCoeffs");
        word weightsFile;

        if (metisCoeffs.readIfPresent("method", method))
        {
            if (method != "recursive" && method != "k-way")
            {
                FatalErrorIn("metisDecomp::decompose()")
                    << "Method " << method << " in metisCoeffs in dictionary : "
                    << decompositionDict_.name()
                    << " should be 'recursive' or 'k-way'"
                    << exit(FatalError);
            }

            Info<< "metisDecomp : Using Metis method     " << method
                << nl << endl;
        }

        if (metisCoeffs.readIfPresent("options", options))
        {
            if (options.size() != 5)
            {
                FatalErrorIn("metisDecomp::decompose()")
                    << "Number of options in metisCoeffs in dictionary : "
                    << decompositionDict_.name()
                    << " should be 5"
                    << exit(FatalError);
            }

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

void Foam::simpleTwoStroke::calcMovingMasks() const
{
    if (debug)
    {
        Info<< "void movingSquaresTM::calcMovingMasks() const : "
            << "Calculating point and cell masks"
            << endl;
    }

    if (movingPointsMaskPtr_)
    {
        FatalErrorIn("void movingSquaresTM::calcMovingMasks() const")
            << "point mask already calculated"
            << abort(FatalError);
    }

    // Set the point mask
    movingPointsMaskPtr_ = new scalarField(allPoints().size(), 0);
    scalarField& movingPointsMask = *movingPointsMaskPtr_;

    const cellList& c = cells();
    const faceList& f = allFaces();

    const labelList& cellAddr =
        cellZones()[cellZones().findZoneID("movingCells")];

    forAll (cellAddr, cellI)
    {
        const cell& curCell = c[cellAddr[cellI]];

        forAll (curCell, faceI)
        {
            // Mark all the points as moving
            const face& curFace = f[curCell[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 1;
            }
        }
    }


    if(foundScavPorts())
    {

        const word innerScavZoneName
        (
            scavInCylPatchName_  + "Zone"
        );

        const labelList& innerScavAddr =
            faceZones()[faceZones().findZoneID(innerScavZoneName)];

        forAll (innerScavAddr, faceI)
        {
            const face& curFace = f[innerScavAddr[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 1;
            }
        }

        const word outerScavZoneName
        (
            scavInPortPatchName_ + "Zone"
        );

        const labelList& outerScavAddr =
            faceZones()[faceZones().findZoneID(outerScavZoneName)];

        forAll (outerScavAddr, faceI)
        {
            const face& curFace = f[outerScavAddr[faceI]];

            forAll (curFace, pointI)
            {
                movingPointsMask[curFace[pointI]] = 0;
            }
        }

    }

}

Foam::scalar Foam::BisectionRoot<Func>::root
(
    const scalar x0,
    const scalar x1
) const
{
    scalar f, fMid, dx, rtb, xMid;

    f = f_(x0);
    fMid = f_(x1);

    if (f*fMid >= 0)
    {
        FatalErrorIn
        (
            "Foam::scalar Foam::BisectionRoot<Func>::root\n"
            "(\n"
            "    const scalar x0,\n"
            "    const scalar x1\n"
            ") const"
        )   << "Root is not bracketed.  f(x0) = " << f << " f(x1) = " << fMid
            << abort(FatalError);
    }

    // Orient the search such that f > 0 lies at x + dx
    if (f < 0)
    {
        dx = x1 - x0;
        rtb = x0;
    }
    else
    {
        dx = x0 - x1;
        rtb = x1;
    }

    for (label nIter = 0; nIter < maxIter; nIter++)
    {
        dx *= 0.5;
        xMid = rtb + dx;

        fMid = f_(xMid);

        if (fMid <= 0)
        {
            rtb = xMid;
        }

        if (mag(dx) < eps_ || mag(fMid) < SMALL)
        {
            return rtb;
        }
    }

    FatalErrorIn
    (
        "Foam::scalar Foam::BisectionRoot<Func>::root\n"
        "(\n"
        "    const scalar x0,\n"
        "    const scalar x1\n"
        ") const"
    )   << "Maximum number of iterations exceeded" << abort(FatalError);

    // Dummy return to keep compiler happy
    return x0;
}

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

    label nBoundaryFaces=0;
    forAll(cells, celli)
    {
        nBoundaryFaces += cells[celli].nFaces();
    }

    nBoundaryFaces -= 2*bm.nInternalFaces();

    label nDefinedBoundaryFaces=0;
    forAll(patches, patchi)
    {
        nDefinedBoundaryFaces += patches[patchi].size();
    }


    Info<< nl << tab << "Basic statistics" << endl;

    Info<< tab << tab << "Number of internal faces : "
        << bm.nInternalFaces() << endl;

    Info<< tab << tab << "Number of boundary faces : "
        << nBoundaryFaces << endl;

    Info<< tab << tab << "Number of defined boundary faces : "
        << nDefinedBoundaryFaces << endl;

    Info<< tab << tab << "Number of undefined boundary faces : "
        << nBoundaryFaces - nDefinedBoundaryFaces << endl;

    if ((nBoundaryFaces - nDefinedBoundaryFaces) > 0)
    {
        Info<< tab << tab << tab
            << "(Warning : only leave undefined the front and back planes "
            << "of 2D planar geometries!)" << endl;
    }

    Info<< nl << tab << "Checking patch -> block consistency" << endl;


    forAll(patches, patchi)
    {
        const faceList& Patch = patches[patchi];

        forAll(Patch, patchFacei)
        {
            const face& patchFace = Patch[patchFacei];
            bool patchFaceOK = false;

            forAll(cells, celli)
            {
                const labelList& cellFaces = cells[celli];

                forAll(cellFaces, cellFacei)
                {
                    if (patchFace == faces[cellFaces[cellFacei]])
                    {
                        patchFaceOK = true;

                        if
                        (
                            (
                                patchFace.normal(points)
                              & faces[cellFaces[cellFacei]].normal(points)
                            ) < 0.0
                        )
                        {
                            Info<< tab << tab
                                << "Face " << patchFacei
                                << " of patch " << patchi
                                << " (" << patches[patchi].name() << ")"
                                << " points inwards"
                                << endl;

                            blockMeshOK = false;
                        }
                    }
                }
            }

            if (!patchFaceOK)
            {
                Info<< tab << tab
                    << "Face " << patchFacei
                    << " of patch " << patchi
                    << " (" << patches[patchi].name() << ")"
                    << " does not match any block faces" << endl;

                blockMeshOK = false;
            }
        }
    }

    if (!blockMeshOK)
    {
        FatalErrorIn("blockMesh::checkBlockMesh(const polyMesh& bm)")
            << "Block mesh topology incorrect, stopping mesh generation!"
            << exit(FatalError);
    }
}

void Foam::GAMGAgglomeration::agglomerateLduAddressing
(
    const label fineLevelIndex
)
{
    const lduMesh& fineMesh = meshLevel(fineLevelIndex);
    const lduAddressing& fineMeshAddr = fineMesh.lduAddr();

    const unallocLabelList& upperAddr = fineMeshAddr.upperAddr();
    const unallocLabelList& lowerAddr = fineMeshAddr.lowerAddr();

    label nFineFaces = upperAddr.size();

    // Get restriction map for current level
    const labelField& restrictMap = restrictAddressing(fineLevelIndex);

    if (min(restrictMap) == -1)
    {
        FatalErrorIn("GAMGAgglomeration::agglomerateLduAddressing")
            << "min(restrictMap) == -1" << exit(FatalError);
    }

    if (restrictMap.size() != fineMeshAddr.size())
    {
        FatalErrorIn
        (
            "GAMGAgglomeration::agglomerateLduAddressing"
            "(const label fineLevelIndex)"
        )   << "restrict map does not correspond to fine level. " << endl
            << " Sizes: restrictMap: " << restrictMap.size()
            << " nEqns: " << fineMeshAddr.size()
            << abort(FatalError);
    }


    // Get the number of coarse cells
    const label nCoarseCells = nCells_[fineLevelIndex];

    // Storage for coarse cell neighbours and coefficients

    // Guess initial maximum number of neighbours in coarse cell
    label maxNnbrs = 10;

    // Number of faces for each coarse-cell
    labelList cCellnFaces(nCoarseCells, 0);

    // Setup initial packed storage for coarse-cell faces
    labelList cCellFaces(maxNnbrs*nCoarseCells);

    // Create face-restriction addressing
    faceRestrictAddressing_.set(fineLevelIndex, new labelList(nFineFaces));
    labelList& faceRestrictAddr = faceRestrictAddressing_[fineLevelIndex];

    // Initial neighbour array (not in upper-triangle order)
    labelList initCoarseNeighb(nFineFaces);

    // Counter for coarse faces
    label nCoarseFaces = 0;

    // Loop through all fine faces
    forAll (upperAddr, fineFacei)
    {
        label rmUpperAddr = restrictMap[upperAddr[fineFacei]];
        label rmLowerAddr = restrictMap[lowerAddr[fineFacei]];

        if (rmUpperAddr == rmLowerAddr)
        {
            // For each fine face inside of a coarse cell keep the address
            // of the cell corresponding to the face in the faceRestrictAddr
            // as a negative index
            faceRestrictAddr[fineFacei] = -(rmUpperAddr + 1);
        }
        else
        {
            // this face is a part of a coarse face

            label cOwn = rmUpperAddr;
            label cNei = rmLowerAddr;

            // get coarse owner and neighbour
            if (rmUpperAddr > rmLowerAddr)
            {
                cOwn = rmLowerAddr;
                cNei = rmUpperAddr;
            }

            // check the neighbour to see if this face has already been found
            label* ccFaces = &cCellFaces[maxNnbrs*cOwn];

            bool nbrFound = false;
            label& ccnFaces = cCellnFaces[cOwn];

            for (int i=0; i<ccnFaces; i++)
            {
                if (initCoarseNeighb[ccFaces[i]] == cNei)
                {
                    nbrFound = true;
                    faceRestrictAddr[fineFacei] = ccFaces[i];
                    break;
                }
//.........这里部分代码省略.........

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

        // Check that the addressing is valid
        if (min(mfc) < 0 || min(sfc) < 0)
        {
            if (debug)
            {
                forAll (mfc, faceI)
                {
                    if (mfc[faceI] < 0)
                    {
                        Pout << "No cell next to master patch face " << faceI
                            << ".  Global face no: " << mfc[faceI]
                            << " own: " << own[masterPatchFaces[faceI]]
                            << " nei: " << nei[masterPatchFaces[faceI]]
                            << " flip: " << masterFlip[faceI] << endl;
                    }
                }

                forAll (sfc, faceI)
                {
                    if (sfc[faceI] < 0)
                    {
                        Pout << "No cell next to slave patch face " << faceI
                            << ".  Global face no: " << sfc[faceI]
                            << " own: " << own[slavePatchFaces[faceI]]
                            << " nei: " << nei[slavePatchFaces[faceI]]
                            << " flip: " << slaveFlip[faceI] << endl;
                    }
                }
            }

            FatalErrorIn
            (
                "void slidingInterface::calcAttachedAddressing()"
                "const"
            )   << "Error is zone face-cell addressing.  Probable error in "
                << "decoupled mesh or sliding interface definition."
                << abort(FatalError);
        }

        // Calculate stick-out faces
        const labelListList& pointFaces = mesh.pointFaces();

        // Master side
        labelHashSet masterStickOutFaceMap
        (
            primitiveMesh::facesPerCell_*(masterPatch.size())
        );

        const labelList& masterMeshPoints = masterPatch.meshPoints();

        forAll (masterMeshPoints, pointI)
        {
            const labelList& curFaces = pointFaces[masterMeshPoints[pointI]];

            forAll (curFaces, faceI)
            {
                // Check if the face belongs to the master face zone;
                // if not add it
                if
                (
                    faceZones.whichZone(curFaces[faceI])
                 != masterFaceZoneID_.index()
                )

Foam::semiTransSurfaceCoupledFvPatchScalarField::
semiTransSurfaceCoupledFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    mixedFvPatchScalarField(p, iF),
    neighbourFieldName_(dict.lookup("neighbourFieldName")),
{
	
	
    if (!isA<directMappedPatchBase>(this->patch().patch()))
    {
        FatalErrorIn
        (
            "semiTransSurfaceCoupledFvPatchScalarField::"
            "semiTransSurfaceCoupledFvPatchScalarField\n"
            "(\n"
            "    const fvPatch& p,\n"
            "    const DimensionedField<scalar, volMesh>& iF,\n"
            "    const dictionary& dict\n"
            ")\n"
        )   << "\n    patch type '" << p.type()
            << "' not type '" << directMappedPatchBase::typeName << "'"
            << "\n    for patch " << p.name()
            << " of field " << dimensionedInternalField().name()
            << " in file " << dimensionedInternalField().objectPath()
            << exit(FatalError);
    }
    
    
    

    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));

    if (dict.found("refValue"))
    {
		
		
        // Full restart
        refValue() = scalarField("refValue", dict, p.size());
        refGrad()  = scalarField("refGradient", dict, p.size());
        valueFraction() = scalarField("valueFraction", dict, p.size());
        
        
    }
    else
    {
    
        // Start from user entered data. Assume fixedValue.
        refValue() = *this;
        refGrad() = 0.0;
        valueFraction() = 1.0;
        
        
    }
    
    
}

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

                autoPtr<compressible::LESModel> sgsModel
                (
                    compressible::LESModel::New(rho, U, phi, thermo())
                );

                PePtr.set
                (
                    new surfaceScalarField
                    (
                        IOobject
                        (
                            "Pe",
                            runTime.timeName(),
                            mesh,
                            IOobject::NO_READ
                        ),
                        mag(phi)
                        /(
                            mesh.magSf()
                            * mesh.surfaceInterpolation::deltaCoeffs()
                            * fvc::interpolate(sgsModel->muEff())
                        )
                    )
                );
            }
            else
            {
                IOdictionary transportProperties
                (
                    IOobject
                    (
                        "transportProperties",
                        runTime.constant(),
                        mesh,
                        IOobject::MUST_READ,
                        IOobject::NO_WRITE
                    )
                );

                dimensionedScalar mu(transportProperties.lookup("mu"));

                PePtr.set
                (
                    new surfaceScalarField
                    (
                        IOobject
                        (
                            "Pe",
                            runTime.timeName(),
                            mesh,
                            IOobject::NO_READ
                        ),
                        mesh.surfaceInterpolation::deltaCoeffs()
                        * (mag(phi)/(mesh.magSf()))*(runTime.deltaT()/mu)
                    )
                );
            }
        }
        else
        {
            FatalErrorIn(args.executable())
                    << "Incorrect dimensions of phi: " << phi.dimensions()
                    << abort(FatalError);
        }


        // can also check how many cells exceed a particular Pe limit
        /*
        {
            label count = 0;
            label PeLimit = 200;
            forAll(PePtr(), i)
            {
                if (PePtr()[i] > PeLimit)
                {
                    count++;
                }

            }

            Info<< "Fraction > " << PeLimit << " = "
                << scalar(count)/Pe.size() << endl;
        }
        */

        Info << "Pe max : " << max(PePtr()).value() << endl;

        if (writeResults)
        {
            PePtr().write();
        }
    }
    else
    {
        Info<< "    No phi" << endl;
    }

    Info<< "\nEnd\n" << endl;
}

void Pstream::exchange
(
    const List<Container>& sendBufs,
    List<Container>& recvBufs,
    labelListList& sizes,
    const int tag,
    const bool block
)
{
    if (!contiguous<T>())
    {
        FatalErrorIn
        (
            "Pstream::exchange(..)"
        )   << "Continuous data only." << Foam::abort(FatalError);
    }

    if (sendBufs.size() != UPstream::nProcs())
    {
        FatalErrorIn
        (
            "Pstream::exchange(..)"
        )   << "Size of list:" << sendBufs.size()
            << " does not equal the number of processors:"
            << UPstream::nProcs()
            << Foam::abort(FatalError);
    }

    sizes.setSize(UPstream::nProcs());
    labelList& nsTransPs = sizes[UPstream::myProcNo()];
    nsTransPs.setSize(UPstream::nProcs());

    forAll(sendBufs, procI)
    {
        nsTransPs[procI] = sendBufs[procI].size();
    }

    // Send sizes across. Note: blocks.
    combineReduce(sizes, UPstream::listEq(), tag);

    if (Pstream::parRun())
    {
        label startOfRequests = Pstream::nRequests();

        // Set up receives
        // ~~~~~~~~~~~~~~~

        recvBufs.setSize(sendBufs.size());
        forAll(sizes, procI)
        {
            label nRecv = sizes[procI][UPstream::myProcNo()];

            if (procI != Pstream::myProcNo() && nRecv > 0)
            {
                recvBufs[procI].setSize(nRecv);
                UIPstream::read
                (
                    UPstream::nonBlocking,
                    procI,
                    reinterpret_cast<char*>(recvBufs[procI].begin()),
                    nRecv*sizeof(T),
                    tag
                );
            }
        }


        // Set up sends
        // ~~~~~~~~~~~~

        forAll(sendBufs, procI)
        {
            if (procI != Pstream::myProcNo() && sendBufs[procI].size() > 0)
            {
                if
                (
                   !UOPstream::write
                    (
                        UPstream::nonBlocking,
                        procI,
                        reinterpret_cast<const char*>(sendBufs[procI].begin()),
                        sendBufs[procI].size()*sizeof(T),
                        tag
                    )
                )
                {
                    FatalErrorIn("Pstream::exchange(..)")
                        << "Cannot send outgoing message. "
                        << "to:" << procI << " nBytes:"
                        << label(sendBufs[procI].size()*sizeof(T))
                        << Foam::abort(FatalError);
                }
            }
        }


        // Wait for all to finish
        // ~~~~~~~~~~~~~~~~~~~~~~

        if (block)
//.........这里部分代码省略.........

Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
{
    if (obr_.foundObject<compressible::RASModel>("RASProperties"))
    {
        const compressible::RASModel& ras
            = obr_.lookupObject<compressible::RASModel>("RASProperties");

        return ras.devRhoReff();
    }
    else if (obr_.foundObject<incompressible::RASModel>("RASProperties"))
    {
        const incompressible::RASModel& ras
            = obr_.lookupObject<incompressible::RASModel>("RASProperties");

        return rho()*ras.devReff();
    }
    else if (obr_.foundObject<compressible::LESModel>("LESProperties"))
    {
        const compressible::LESModel& les =
        obr_.lookupObject<compressible::LESModel>("LESProperties");

        return les.devRhoReff();
    }
    else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
    {
        const incompressible::LESModel& les
            = obr_.lookupObject<incompressible::LESModel>("LESProperties");

        return rho()*les.devReff();
    }
    else if (obr_.foundObject<fluidThermo>("thermophysicalProperties"))
    {
        const fluidThermo& thermo =
             obr_.lookupObject<fluidThermo>("thermophysicalProperties");

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
    }
    else if
    (
        obr_.foundObject<singlePhaseTransportModel>("transportProperties")
    )
    {
        const singlePhaseTransportModel& laminarT =
            obr_.lookupObject<singlePhaseTransportModel>
            ("transportProperties");

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
    }
    else if (obr_.foundObject<dictionary>("transportProperties"))
    {
        const dictionary& transportProperties =
             obr_.lookupObject<dictionary>("transportProperties");

        dimensionedScalar nu(transportProperties.lookup("nu"));

        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);

        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
    }
    else
    {
        FatalErrorIn("forces::devRhoReff()")
            << "No valid model for viscous stress calculation."
            << exit(FatalError);

        return volSymmTensorField::null();
    }
}

// Check 1D/2Dness of edges. Gets passed the non-empty directions and
// checks all edges in the mesh whether they:
// - have no component in a non-empty direction or
// - are only in a singe non-empty direction.
// Empty direction info is passed in as a vector of labels (synchronised)
// which are 1 if the direction is non-empty, 0 if it is.
bool Foam::polyMesh::checkEdgeAlignment
(
    const pointField& p,
    const bool report,
    const Vector<label>& directions,
    labelHashSet* setPtr
) const
{
    if (debug)
    {
        Info<< "bool polyMesh::checkEdgeAlignment("
            << "const bool, const Vector<label>&, labelHashSet*) const: "
            << "checking edge alignment" << endl;
    }

    label nDirs = 0;
    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
    {
        if (directions[cmpt] == 1)
        {
            nDirs++;
        }
        else if (directions[cmpt] != 0)
        {
            FatalErrorIn
            (
                "polyMesh::checkEdgeAlignment"
                "(const bool, const Vector<label>&, labelHashSet*)"
            )   << "directions should contain 0 or 1 but is now " << directions
                << exit(FatalError);
        }
    }

    if (nDirs == vector::nComponents)
    {
        return false;
    }


    const faceList& fcs = faces();

    EdgeMap<label> edgesInError;

    forAll(fcs, faceI)
    {
        const face& f = fcs[faceI];

        forAll(f, fp)
        {
            label p0 = f[fp];
            label p1 = f.nextLabel(fp);
            if (p0 < p1)
            {
                vector d(p[p1]-p[p0]);
                scalar magD = mag(d);

                if (magD > ROOTVSMALL)
                {
                    d /= magD;

                    // Check how many empty directions are used by the edge.
                    label nEmptyDirs = 0;
                    label nNonEmptyDirs = 0;
                    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
                    {
                        if (mag(d[cmpt]) > 1e-6)
                        {
                            if (directions[cmpt] == 0)
                            {
                                nEmptyDirs++;
                            }
                            else
                            {
                                nNonEmptyDirs++;
                            }
                        }
                    }

                    if (nEmptyDirs == 0)
                    {
                        // Purely in ok directions.
                    }
                    else if (nEmptyDirs == 1)
                    {
                        // Ok if purely in empty directions.
                        if (nNonEmptyDirs > 0)
                        {
                            edgesInError.insert(edge(p0, p1), faceI);
                        }
                    }
                    else if (nEmptyDirs > 1)
                    {
                        // Always an error
                        edgesInError.insert(edge(p0, p1), faceI);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    solutionD_(Vector<label>::zero),
    pointZones_
    (
        IOobject
        (
            "pointZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    faceZones_
    (
        IOobject
        (
            "faceZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    cellZones_
    (
        IOobject
        (
            "cellZones",
            instance(),
            meshSubDir,
            *this,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        *this,
        0
    ),
    globalMeshDataPtr_(NULL),
    moving_(false),
    changing_(false),
    curMotionTimeIndex_(time().timeIndex()),
    oldAllPointsPtr_(NULL),
    oldPointsPtr_(NULL)
{
    // Check if the faces and cells are valid
    forAll (allFaces_, faceI)
    {
        const face& curFace = allFaces_[faceI];

        if (min(curFace) < 0 || max(curFace) > allPoints_.size())
        {
            FatalErrorIn
            (
                "polyMesh::polyMesh\n"
                "(\n"
                "    const IOobject&,\n"
                "    const Xfer<pointField>&,\n"
                "    const Xfer<faceList>&,\n"
                "    const Xfer<cellList>&\n"
                ")\n"
            )   << "Face " << faceI << "contains vertex labels out of range: "
                << curFace << " Max point index = " << allPoints_.size()
                << abort(FatalError);
        }
    }

    // transfer in cell list
    cellList cLst(cells);

    // Check if cells are valid
    forAll (cLst, cellI)
    {
        const cell& curCell = cLst[cellI];

        if (min(curCell) < 0 || max(curCell) > allFaces_.size())
        {
            FatalErrorIn
            (
                "polyMesh::polyMesh\n"
                "(\n"
                "    const IOobject&,\n"
                "    const Xfer<pointField>&,\n"
                "    const Xfer<faceList>&,\n"
                "    const Xfer<cellList>&\n"
                ")\n"
            )   << "Cell " << cellI << "contains face labels out of range: "
                << curCell << " Max face index = " << allFaces_.size()
                << abort(FatalError);
        }
    }

    // Set the primitive mesh
    initMesh(cLst);
}

void decomposeFaces::decomposeMeshFaces(const boolList& decomposeFace)
{
    done_ = false;

    newFacesForFace_.setSize(mesh_.faces().size());
    forAll(newFacesForFace_, fI)
        newFacesForFace_.setRowSize(fI, 0);

    const label nIntFaces = mesh_.nInternalFaces();
    label nFaces(0), nPoints = mesh_.points().size();
    point p;

    pointFieldPMG& points = mesh_.points();
    forAll(decomposeFace, fI)
        if( decomposeFace[fI] )
            ++nFaces;

    points.setSize(nPoints + nFaces);

    polyMeshGenModifier meshModifier(mesh_);
    faceListPMG& faces = meshModifier.facesAccess();

    if( decomposeFace.size() != faces.size() )
        FatalErrorIn
        (
            "void decomposeFaces::decomposeMeshFaces(const boolList&)"
        ) << "Incorrect size of the decomposeFace list!" << abort(FatalError);

    nFaces = 0;
    VRWGraph newFaces;

    //- decompose internal faces
    for(label faceI=0;faceI<nIntFaces;++faceI)
    {
        const face& f = faces[faceI];

        if( decomposeFace[faceI] )
        {
            # ifdef DEBUGDec
            Info << "Decomposing internal face " << faceI << " with nodes "
                << f << endl;
            # endif

            FixedList<label, 3> newF;

            forAll(f, pI)
            {
                newF[0] = f[pI];
                newF[1] = f.nextLabel(pI);
                newF[2] = nPoints;

                # ifdef DEBUGDec
                Info << "Storing face " << newF << " with label "
                    << nFaces << endl;
                # endif

                newFaces.appendList(newF);
                newFacesForFace_.append(faceI, nFaces++);
            }

            p = f.centre(points);
            points[nPoints] = p;
            ++nPoints;
        }
        else
        {

void Foam::massFlowTotalPressureFvPatchScalarField::updateCoeffs
(
    const scalarField& p0p,
    const vectorField& Up
)
{
    //Info << "start of updateCoeff" << endl;
    if (updated())
    {
        return;
    }
    /*
    const fvsPatchField<scalar>& phip =
        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);

    if (psiName_ == "none" && rhoName_ == "none")
    {
        operator==(p0p - 0.5*(1.0 - pos(phip))*magSqr(Up));
    }
    else if (rhoName_ == "none")
    {
        const fvPatchField<scalar>& psip =
            patch().lookupPatchField<volScalarField, scalar>(psiName_);

        if (gamma_ > 1.0)
        {
            scalar gM1ByG = (gamma_ - 1.0)/gamma_;

            operator==
            (
                p0p
               /pow
                (
                    (1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)),
                    1.0/gM1ByG
                )
            );
        }
        else
        {
            operator==(p0p/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
        }
    }*/
    if (psiName