自学教程：C++ FatalErrorIn函数代码示例

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_ == "none")    {	//Info << "before field intros" << endl;        const fvPatchField<scalar>& rho =            patch().lookupPatchField<volScalarField, scalar>(rhoName_);	//Info << "soundSpeed" << endl;        const fvPatchField<scalar>& soundSpeed =            patch().lookupPatchField<volScalarField, scalar>("soundSpeed");	//Info << "T" << endl;        const fvPatchField<scalar>& T =            patch().lookupPatchField<volScalarField, scalar>("T");        const scalarField Apatch = patch().magSf();	        const scalar R = 8.3144;	//Info << "before Un" << endl;        const scalarField Un = massFlow_/Apatch/rho;	const scalarField M = Un/soundSpeed;	//Info << "before M" << endl;        const scalarField realGamma = soundSpeed*soundSpeed*Mave_/R/T;	//Info << "before pTotCoeff" << endl;        const scalarField pTotCoeff = pow(1.0+0.5*(realGamma-1.0)*M*M,realGamma/(realGamma-1.0));	//Info << "before pExtrap" << endl;        const scalarField pExtrap = rho*R*T/Mave_;        //operator==(p0p - 0.5*rho*(1.0 - pos(phip))*magSqr(Up)); 	operator==(pExtrap/pTotCoeff);    }    else    {        FatalErrorIn        (            "massFlowTotalPressureFvPatchScalarField::updateCoeffs()"        )   << " rho or psi set inconsistently, rho = " << rhoName_            << ", psi = " << psiName_ << "./n"            << "    Set either rho or psi or neither depending on the "               "definition of total pressure." << nl            << "    Set the unused variable(s) to 'none'./n"            << "    on patch " << this->patch().name()            << " of field " << this->dimensionedInternalField().name()            << " in file " << this->dimensionedInternalField().objectPath()            << exit(FatalError);    }    fixedValueFvPatchScalarField::updateCoeffs();}

void Foam::layerARGambit::addZonesAndModifiers(){    // Add the zones and mesh modifiers to operate piston motion    if    (        pointZones().size() > 0     || faceZones().size() > 0     || cellZones().size() > 0    )    {        Info<< "void layerARGambit::addZonesAndModifiers() : "            << "Zones and modifiers already present.  Skipping."            << endl;        if (topoChanger_.size() == 0)        {            FatalErrorIn            (                "void layerARGambit::addZonesAndModifiers()"            )   << "Mesh modifiers not read properly"                << abort(FatalError);        }        setVirtualPistonPosition();        checkAndCalculate();        return;    }    checkAndCalculate();    Info<< "Time = " << engTime().theta() << endl        << "Adding zones to the engine mesh" << endl;    //fz = 1: faces where layer are added/removed    //pz = 2: points below the virtual piston faces and head points    List<pointZone*> pz(2);    List<faceZone*> fz(1);    List<cellZone*> cz(0);    label nPointZones = 0;    label nFaceZones = 0;    // Add the piston zone    if (piston().patchID().active() && offSet() > SMALL)    {        // Piston position        label pistonPatchID = piston().patchID().index();        scalar zPist = max(boundary()[pistonPatchID].patch().localPoints()).z();        scalar zPistV = zPist + offSet();        labelList zone1(faceCentres().size());        boolList flipZone1(faceCentres().size(), false);        label nZoneFaces1 = 0;        bool foundAtLeastOne = false;        scalar zHigher = GREAT;        scalar dh = GREAT;        scalar dl = GREAT;        forAll (faceCentres(), faceI)        {            scalar zc = faceCentres()[faceI].z();            vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);            scalar dd = n & vector(0,0,1);            if (mag(dd) > 0.1)            {                if (zPistV - zc > 0 && zPistV - zc < dl)                {                    dl = zPistV - zc;                }                if (zc - zPistV > 0 && zc - zPistV < dh)                {                    zHigher = zc;                    dh = zc - zHigher;                }                if                (                    zc > zPistV - delta()                    && zc < zPistV + delta()                )                {                    foundAtLeastOne = true;                    if ((faceAreas()[faceI] & vector(0,0,1)) < 0)                    {                        flipZone1[nZoneFaces1] = true;                    }                    zone1[nZoneFaces1] = faceI;                    nZoneFaces1++;//.........这里部分代码省略.........

void Foam::KRR4::solve(    scalar& x,    scalarField& y,    scalarField& dydx,    const scalar eps,    const scalarField& yScale,    const scalar hTry,    scalar& hDid,    scalar& hNext) const{    scalar xTemp = x;    yTemp_ = y;    dydxTemp_ = dydx;    const label nEqns = ode_.nEqns();    ode_.jacobian(xTemp, yTemp_, dfdx_, dfdy_);    scalar h = hTry;    for (register label jtry=0; jtry<maxtry; jtry++)    {        for (register label i=0; i<nEqns; i++)        {            for (register label j=0; j<nEqns; j++)            {                a_[i][j] = -dfdy_[i][j];            }            a_[i][i] += 1.0/(gamma*h);        }        simpleMatrix<scalar>::LUDecompose(a_, pivotIndices_);        for (register label i=0; i<nEqns; i++)        {            g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];        }        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g1_);        for (register label i=0; i<nEqns; i++)        {            y[i] = yTemp_[i] + a21*g1_[i];        }        x = xTemp + a2X*h;        ode_.derivatives(x, y, dydx_);        for (register label i=0; i<nEqns; i++)        {            g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;        }        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g2_);        for (register label i=0; i<nEqns; i++)        {            y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];        }        x = xTemp + a3X*h;        ode_.derivatives(x, y, dydx_);        for (register label i=0; i<nEqns; i++)        {            g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;        }        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g3_);        for (register label i=0; i<nEqns; i++)        {            g4_[i] = dydx_[i] + h*c4X*dfdx_[i]                + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;        }        simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g4_);        for (register label i=0; i<nEqns; i++)        {            y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];            yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];        }        x = xTemp + h;        if (x == xTemp)        {            FatalErrorIn("ODES::KRR4")                << "stepsize not significant"                << exit(FatalError);        }        scalar maxErr = 0.0;        for (register label i=0; i<nEqns; i++)        {            maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));        }//.........这里部分代码省略.........

int main(int argc, char *argv[]){#   include "addRegionOption.H"#   include "addLoadFunctionPlugins.H"    argList::validOptions.insert("overwrite", "");    argList::validOptions.insert("expression","expression to write");    argList::validOptions.insert("dictExt","Extension to the dictionary");    argList::validOptions.insert("relative", "");#   include "setRootCase.H"    printSwakVersion();#   include "createTime.H"    runTime.functionObjects().off();    Foam::instantList timeDirs = Foam::timeSelector::select0(runTime, args);#   include "createNamedMesh.H"#   include "loadFunctionPlugins.H"    const word oldInstance = mesh.pointsInstance();    bool overwrite    = args.optionFound("overwrite");    bool relative     = false;    pointField newPoints;    if (!overwrite)    {        runTime++;    }    if(args.optionFound("expression")) {        if(args.optionFound("dictExt")) {            FatalErrorIn(args.executable())                << "Can't specify 'dictExt' and 'expression' at the same time"                    << endl                    << exit(FatalError);        }        relative=args.optionFound("relative");        exprString expression(            args.options()["expression"],            dictionary::null        );        FieldValueExpressionDriver driver(            runTime.timeName(),            runTime,            mesh        );        // no clearVariables needed here        driver.parse(expression);        if(!driver.resultIsTyp<pointVectorField>()) {            FatalErrorIn(args.executable())                << "Expression " << expression                    << " does not evaluate to a pointVectorField but a "                    << driver.typ()                    << endl                    << exit(FatalError);        }        newPoints=driver.getResult<pointVectorField>().internalField();    } else {        Info << "Dictionary mode" << nl << endl;        if(args.optionFound("relative")) {            FatalErrorIn(args.executable())                << "Option 'relative' not allowed in dictionary-mode"                    << endl                    << exit(FatalError);        }        word dictName("funkyWarpMeshDict");        if(args.optionFound("dictExt")) {            dictName+="."+word(args.options()["dictExt"]);        }        IOdictionary warpDict            (                IOobject                (                    dictName,                    runTime.system(),                    mesh,                    IOobject::MUST_READ,                    IOobject::NO_WRITE                )            );        const word mode(warpDict.lookup("mode"));        if(mode=="set") {            relative=readBool(warpDict.lookup("relative"));            exprString expression(                warpDict.lookup("expression"),                warpDict            );            FieldValueExpressionDriver driver(                runTime.timeName(),                runTime,                mesh            );//.........这里部分代码省略.........

const Foam::Tuple2<Foam::scalar, Type>&Foam::interpolationTable<Type>::operator[](const label i) const{    label ii = i;    label n  = this->size();    if (n <= 1)    {        ii = 0;    }    else if (ii < 0)    {        switch (boundsHandling_)        {            case interpolationTable::ERROR:            {                FatalErrorIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "index (" << ii << ") underflow" << nl                    << exit(FatalError);                break;            }            case interpolationTable::WARN:            {                WarningIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "index (" << ii << ") underflow" << nl                    << "    Continuing with the first entry"                    << endl;                // fall-through to 'CLAMP'            }            case interpolationTable::CLAMP:            {                ii = 0;                break;            }            case interpolationTable::REPEAT:            {                while (ii < 0)                {                    ii += n;                }                break;            }        }    }    else if (ii >= n)    {        switch (boundsHandling_)        {            case interpolationTable::ERROR:            {                FatalErrorIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "index (" << ii << ") overflow" << nl                    << exit(FatalError);                break;            }            case interpolationTable::WARN:            {                WarningIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "index (" << ii << ") overflow" << nl                    << "    Continuing with the last entry"                    << endl;                // fall-through to 'CLAMP'            }            case interpolationTable::CLAMP:            {                ii = n - 1;                break;            }            case interpolationTable::REPEAT:            {                while (ii >= n)                {                    ii -= n;                }                break;            }        }    }    return List<Tuple2<scalar, Type> >::operator[](ii);}

void Foam::thoboisSliding::checkAndCalculate(){        label pistonIndex = -1;    bool foundPiston = false;    label linerIndex = -1;    bool foundLiner = false;    label cylinderHeadIndex = -1;    bool foundCylinderHead = false;            forAll(boundary(), i)    {        if (boundary()[i].name() == "piston")        {            pistonIndex = i;            foundPiston = true;        }        else if (boundary()[i].name() == "liner")        {            linerIndex = i;            foundLiner = true;        }        else if (boundary()[i].name() == "cylinderHead")        {            cylinderHeadIndex = i;            foundCylinderHead = true;        }    }        reduce(foundPiston, orOp<bool>());    reduce(foundLiner, orOp<bool>());    reduce(foundCylinderHead, orOp<bool>());    if (!foundPiston)    {        FatalErrorIn("Foam::verticalValvesGambit::checkAndCalculate()")            << " : cannot find piston patch"            << abort(FatalError);    }    if (!foundLiner)    {         FatalErrorIn("Foam::verticalValvesGambit::checkAndCalculate()")            << " : cannot find liner patch"            << abort(FatalError);    }    if (!foundCylinderHead)    {         FatalErrorIn("Foam::verticalValvesGambit::checkAndCalculate()")            << " : cannot find cylinderHead patch"            << exit(FatalError);    }    {        if (linerIndex != -1)        {            pistonPosition() =                max(boundary()[pistonIndex].patch().localPoints()).z();        }        reduce(pistonPosition(), minOp<scalar>());        if (cylinderHeadIndex != -1)        {            deckHeight() = min            (                boundary()[cylinderHeadIndex].patch().localPoints()            ).z();        }        reduce(deckHeight(), minOp<scalar>());        Info<< "deckHeight: " << deckHeight() << nl            << "piston position: " << pistonPosition() << endl;    }            } 

Type Foam::interpolationTable<Type>::operator()(const scalar value) const{    label n = this->size();    if (n <= 1)    {        return List<Tuple2<scalar, Type> >::operator[](0).second();    }    scalar minLimit = List<Tuple2<scalar, Type> >::operator[](0).first();    scalar maxLimit = List<Tuple2<scalar, Type> >::operator[](n-1).first();    scalar lookupValue = value;    if (lookupValue < minLimit)    {        switch (boundsHandling_)        {            case interpolationTable::ERROR:            {                FatalErrorIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const scalar) const"                )   << "value (" << lookupValue << ") underflow" << nl                    << exit(FatalError);                break;            }            case interpolationTable::WARN:            {                WarningIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const scalar) const"                )   << "value (" << lookupValue << ") underflow" << nl                    << "    Continuing with the first entry"                    << endl;                // fall-through to 'CLAMP'            }            case interpolationTable::CLAMP:            {                return List<Tuple2<scalar, Type> >::operator[](0).second();                break;            }            case interpolationTable::REPEAT:            {                // adjust lookupValue to >= minLimit                scalar span = maxLimit-minLimit;                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;                break;            }        }    }    else if (lookupValue >= maxLimit)    {        switch (boundsHandling_)        {            case interpolationTable::ERROR:            {                FatalErrorIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "value (" << lookupValue << ") overflow" << nl                    << exit(FatalError);                break;            }            case interpolationTable::WARN:            {                WarningIn                (                    "Foam::interpolationTable<Type>::operator[]"                    "(const label) const"                )   << "value (" << lookupValue << ") overflow" << nl                    << "    Continuing with the last entry"                    << endl;                // fall-through to 'CLAMP'            }            case interpolationTable::CLAMP:            {                return List<Tuple2<scalar, Type> >::operator[](n-1).second();                break;            }            case interpolationTable::REPEAT:            {                // adjust lookupValue <= maxLimit                scalar span = maxLimit-minLimit;                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;                break;            }        }    }    label lo = 0;    label hi = 0;    // look for the correct range    for (label i = 0; i < n; ++i)    {        if (lookupValue >= List<Tuple2<scalar, Type> >::operator[](i).first())        {//.........这里部分代码省略.........

void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs(){    if (this->updated())    {        return;    }    const scalarField& Tp =        patch().lookupPatchField<volScalarField, scalar>(TName_);    const radiationModel& radiation =        db().lookupObject<radiationModel>("radiationProperties");    const fvDOM& dom(refCast<const fvDOM>(radiation));    label rayId = -1;    label lambdaId = -1;    dom.setRayIdLambdaId(dimensionedInternalField().name(), rayId, lambdaId);    const label patchI = patch().index();    if (dom.nLambda() != 1)    {        FatalErrorIn        (            "Foam::radiation::"            "greyDiffusiveRadiationMixedFvPatchScalarField::updateCoeffs"        )   << " a grey boundary condition is used with a non-grey "            << "absorption model" << nl << exit(FatalError);    }    scalarField& Iw = *this;    vectorField n = patch().Sf()/patch().magSf();    radiativeIntensityRay& ray =        const_cast<radiativeIntensityRay&>(dom.IRay(rayId));    ray.Qr().boundaryField()[patchI] += Iw*(n & ray.dAve());    forAll(Iw, faceI)    {        scalar Ir = 0.0;        for (label rayI=0; rayI < dom.nRay(); rayI++)        {            const vector& d = dom.IRay(rayI).d();            const scalarField& IFace =                dom.IRay(rayI).ILambda(lambdaId).boundaryField()[patchI];            if ((-n[faceI] & d) < 0.0)            {                // q into the wall                const vector& dAve = dom.IRay(rayI).dAve();                Ir += IFace[faceI]*mag(n[faceI] & dAve);            }        }        const vector& d = dom.IRay(rayId).d();        if ((-n[faceI] & d) > 0.0)        {            // direction out of the wall            refGrad()[faceI] = 0.0;            valueFraction()[faceI] = 1.0;            refValue()[faceI] =                (                    Ir*(1.0 - emissivity_)                    + emissivity_*radiation::sigmaSB.value()*pow4(Tp[faceI])                )                /mathematicalConstant::pi;            // Emmited heat flux from this ray direction            ray.Qem().boundaryField()[patchI][faceI] =                refValue()[faceI]*(n[faceI] & ray.dAve());        }        else        {            // direction into the wall            valueFraction()[faceI] = 0.0;            refGrad()[faceI] = 0.0;            refValue()[faceI] = 0.0; //not used            // Incident heat flux on this ray direction            ray.Qin().boundaryField()[patchI][faceI] =                Iw[faceI]*(n[faceI] & ray.dAve());        }    }

void Foam::multiSolver::setInitialSolverDomain(const word& solverDomainName){    if (!solverDomains_.found(solverDomainName))    {        FatalErrorIn("multiSolver::setInitialSolverDomain")            << "Initial solverDomainName '" << solverDomainName << "' does"            << " not exist in multiSolver dictionary.  Found entries are: "            << solverDomains_.toc()            << abort(FatalError);    }    currentSolverDomain_ = solverDomainName;    setSolverDomainControls(currentSolverDomain_);    // Purge all time directories from case directory root    purgeTimeDirs(multiDictRegistry_.path());    // Purge any constant/superLoopData/    fileName superLoopDataPath    (        multiDictRegistry_.path()/multiDictRegistry_.constant()            /"superLoopData"    );    if (exists(superLoopDataPath))    {        rmDir(superLoopDataPath);    }    // Read initial settings and determine data source (from which path the    // initial data is copied, the starting superLoop_, and the current    // globalTime (used to determine globalOffset).  Rules that are applied:    //    //   1. superLoop_ = data source superLoop    //       a. unless data source solverDomain != currentSolverDomain_, in    //          which case, superLoop_ = data source superLoop + 1    //   2. globalTime = data source globalTime.  globalTime does not increment    //      when swapping solver domains.    //   3. startTime = data source local time    //       a. unless data source solverDomain != currentSolverDomain_, in    //          which case, startTime is dictated by the solverDomains    //          subdictionary.    //   4. endTime is determined by the solverDomains subdictionary    //       a. unless the finalStopAt trumps it    // Find initial data source    timeCluster tcSource(initialDataSource());    fileName sourcePath(findInstancePath(tcSource, tcSource.size() - 1));    superLoop_ = tcSource.superLoop();    globalIndex_ = tcSource.globalIndex();    // If starting from initial conditions, superLoop_ = -1    if (superLoop_ < 0) superLoop_ = 0;    scalar globalTime(tcSource.globalValue(tcSource.size() - 1));    scalar localStartTime(tcSource.localValue(tcSource.size() -1));    // Now to apply the exceptions if currentSolverDomain_ != data source    // solverDomain (see long comment above).    if (sourcePath.path().path().name() != currentSolverDomain_)    {        superLoop_++;        globalIndex_++;        switch (startFrom_)        {            case mtsFirstTime:                localStartTime = 0;                break;            case mtsStartTime:                localStartTime = startTime_;                break;            case mtsLatestTimeThisDomain:                {                    timeCluster tcTemp                    (                        findLatestLocalTime                        (                            readSolverDomainTimes(currentSolverDomain_)                        )                    );                    localStartTime = tcTemp.localValue(0);                }                break;            case mtsLatestTimeAllDomains:                localStartTime = globalTime;                break;        }    }    startTime_ = localStartTime;    globalTimeOffset_ = globalTime - startTime_;    // Give multiDictRegistry a time value (required for regIOobject::write()    // to case/[timeValue]    multiDictRegistry_.setTime(startTime_, 0);    // Copy the source data and any previous time directories to//.........这里部分代码省略.........

void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&){    if (axis_ != vector::rootMax)    {        return;    }    if (returnReduce(size(), sumOp<label>()))    {        const vectorField& nf(faceNormals());        n_ = gAverage(nf);        if (debug)        {            Info<< "Patch " << name() << " calculated average normal "                << n_ << endl;        }        // Check the wedge is planar        forAll(nf, faceI)        {            if (magSqr(n_ - nf[faceI]) > SMALL)            {                // only issue warning instead of error so that the case can                // still be read for post-processing                WarningIn                (                    "wedgePolyPatch::calcGeometry(PstreamBuffers&)"                )                    << "Wedge patch '" << name() << "' is not planar." << nl                    << "At local face at "                    << primitivePatch::faceCentres()[faceI]                    << " the normal " << nf[faceI]                    << " differs from the average normal " << n_                    << " by " << magSqr(n_ - nf[faceI]) << nl                    << "Either correct the patch or split it into planar parts"                    << endl;            }        }        centreNormal_ =            vector            (                sign(n_.x())*(max(mag(n_.x()), 0.5) - 0.5),                sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),                sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)            );        centreNormal_ /= mag(centreNormal_);        cosAngle_ = centreNormal_ & n_;        const scalar cnCmptSum =            centreNormal_.x() + centreNormal_.y() + centreNormal_.z();        if (mag(cnCmptSum) < (1 - SMALL))        {            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")                << "wedge " << name()                << " centre plane does not align with a coordinate plane by "                << 1 - mag(cnCmptSum)                << exit(FatalError);        }        axis_ = centreNormal_ ^ n_;        scalar magAxis = mag(axis_);        if (magAxis < SMALL)        {            FatalErrorIn("wedgePolyPatch::calcGeometry(PstreamBuffers&)")                << "wedge " << name()                << " plane aligns with a coordinate plane." << nl                << "    The wedge plane should make a small angle (~2.5deg)"                   " with the coordinate plane" << nl                << "    and the the pair of wedge planes should be symmetric"                << " about the coordinate plane." << nl                << "    Normal of wedge plane is " << n_                << " , implied coordinate plane direction is " << centreNormal_                << exit(FatalError);        }        axis_ /= magAxis;        faceT_ = rotationTensor(centreNormal_, n_);        cellT_ = faceT_ & faceT_;    }}

Foam::label Foam::UIPstream::read(    const commsTypes commsType,    const int fromProcNo,    char* buf,    const std::streamsize bufSize,    const int tag,    const label communicator){    if (debug)    {        Pout<< "UIPstream::read : starting read from:" << fromProcNo            << " tag:" << tag << " comm:" << communicator            << " wanted size:" << label(bufSize)            << " commsType:" << UPstream::commsTypeNames[commsType]            << Foam::endl;    }    if (UPstream::warnComm != -1 && communicator != UPstream::warnComm)    {        Pout<< "UIPstream::read : starting read from:" << fromProcNo            << " tag:" << tag << " comm:" << communicator            << " wanted size:" << label(bufSize)            << " commsType:" << UPstream::commsTypeNames[commsType]            << " warnComm:" << UPstream::warnComm            << Foam::endl;        error::printStack(Pout);    }    if (commsType == blocking || commsType == scheduled)    {        MPI_Status status;        if        (            MPI_Recv            (                buf,                bufSize,                MPI_BYTE,                fromProcNo,                tag,                PstreamGlobals::MPICommunicators_[communicator],                &status            )        )        {            FatalErrorIn            (                "UIPstream::read"                "(const int fromProcNo, char* buf, std::streamsize bufSize)"            )   << "MPI_Recv cannot receive incomming message"                << Foam::abort(FatalError);            return 0;        }        // Check size of message read        int messageSize;        MPI_Get_count(&status, MPI_BYTE, &messageSize);        if (debug)        {            Pout<< "UIPstream::read : finished read from:" << fromProcNo                << " tag:" << tag << " read size:" << label(bufSize)                << " commsType:" << UPstream::commsTypeNames[commsType]                << Foam::endl;        }        if (messageSize > bufSize)        {            FatalErrorIn            (                "UIPstream::read"                "(const int fromProcNo, char* buf, std::streamsize bufSize)"            )   << "buffer (" << label(bufSize)                << ") not large enough for incomming message ("                << messageSize << ')'                << Foam::abort(FatalError);        }        return messageSize;    }    else if (commsType == nonBlocking)    {        MPI_Request request;        if        (            MPI_Irecv            (                buf,                bufSize,                MPI_BYTE,                fromProcNo,                tag,                PstreamGlobals::MPICommunicators_[communicator],                &request//.........这里部分代码省略.........