int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    #include "addRegionOption.H"
    Foam::argList::addOption
    (
        "dir",
        "word",
        "Direction to remove in TKE. -dir x, or y, or z."
    );
    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"
    
    if (args.optionFound("dir")==true)
    {
        word dir(args.optionRead<word>("dir"));
        
        forAll(timeDirs, timeI)
        {
            runTime.setTime(timeDirs[timeI], timeI);
            Info<< "Time = " << runTime.timeName() << endl;
            
            
            IOobject UPrime2MeanHeader
            (
                "UPrime2Mean",
                runTime.timeName(),
                mesh,
                IOobject::MUST_READ
            );
            
            volScalarField TKE2D
            (
                IOobject
                (
                    "TKE2D",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedScalar
                (
                    "dimTKE2D",
                    pow(dimLength,2)/pow(dimTime,2),
                    0
                ),
                "zeroGradient"
            );
            
            if (UPrime2MeanHeader.headerOk())
            {
                Info<< "    Reading average field UPrime2Mean" << endl;
                volSymmTensorField UPrime2Mean(UPrime2MeanHeader, mesh);         
                
                if (dir=="x")
                {
                    Info << "    Calculating 2D turbulent kinetic energy TKE2D without the " 
                         << dir << "component." << endl;
                    TKE2D = 0.5 * (
//                               UPrime2Mean.component(0)
                              UPrime2Mean.component(3)
                            + UPrime2Mean.component(5)
                            );
                }
                else if (dir=="y")
                {
                    Info << "    Calculating 2D turbulent kinetic energy TKE2D without the " 
                         << dir << "component." << endl;
                    TKE2D = 0.5 * (
                              UPrime2Mean.component(0)
//                             + UPrime2Mean.component(3)
                            + UPrime2Mean.component(5)
                            );
                }
                else if (dir=="z")
                {
                    Info << "    Calculating 2D turbulent kinetic energy TKE2D without the " 
                         << dir << "component." << endl;
                    TKE2D = 0.5 * (
                              UPrime2Mean.component(0)
                            + UPrime2Mean.component(3)
//                             + UPrime2Mean.component(5)
                            );
                }
                else
                {
                    Info << "    dir must be x, y or z" << endl;
                }
            
            TKE2D.write();


            }
            else
            {
                Info << "No UPrime2Mean." << endl;
//.........这里部分代码省略.........

// ***********************************************************
// Construct from components for scalarGeneralExchangePhaseChange
scalarGeneralExchange::scalarGeneralExchange
(
    const dictionary& dict,
    cfdemCloud& sm,
    word        dictName
)
:
    forceModel(dict,sm),
    propsDict_(dict.subDict(dictName + "Props")),
    scalarTransportProperties_                  //this is clumsy, but effective
    (
        IOobject
        (
            "scalarTransportProperties",
            sm.mesh().time().constant(),
            sm.mesh(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),
    generalPropsDict_(scalarTransportProperties_.subDict("generalManualProps")),
    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),             //common names/data
    velFieldName_(propsDict_.lookup("velFieldName")),
    tempFieldName_(propsDict_.lookup("tempFieldName")),                             //temperature names/data
    partTempName_(propsDict_.lookup("partTempName")),
    partHeatFluxName_(propsDict_.lookupOrDefault<word>(      "partHeatFluxName", "none")),
    partHeatTransCoeffName_(propsDict_.lookupOrDefault<word>("partHeatTransCoeffName", "none")),
    partHeatFluidName_(propsDict_.lookupOrDefault<word>(     "partHeatFluidName", "none")),
    partDat_(NULL),
    partDatFlux_(NULL),
    partDatTransCoeff_(NULL),
    partDatFluid_(NULL),
    partDatTmpExpl_(NULL),
    partDatTmpImpl_(NULL),
    validPartFlux_(false),
    validPartTransCoeff_(false),
    validPartFluid_(false),
    haveTemperatureEqn_(false),
    useLiMason_(false),
    useGeneralCorrelation_(false),
    lambda_(readScalar(propsDict_.lookup("lambda"))),
    Prandtl_(readScalar(propsDict_.lookup("Prandtl"))),
    eulerianFieldNames_( generalPropsDict_.lookup("eulerianFields")), 
    partSpeciesNames_(propsDict_.lookup("partSpeciesNames")),                       
    partSpeciesFluxNames_(propsDict_.lookup("partSpeciesFluxNames")),
    partSpeciesTransCoeffNames_(propsDict_.lookup("partSpeciesTransCoeffNames")),
    partSpeciesFluidNames_(propsDict_.lookup("partSpeciesFluidNames")),
    DMolecular_(propsDict_.lookup("DMolecular")),
    partHeatFluxPositionInRegister_(-1),
    partHeatTransCoeffPositionInRegister_(-1),
    partHeatFluidPositionInRegister_(-1),
    maxSource_(1e30),
    scaleDia_(1.)

{
    setForceSubModels(propsDict_);
    setupModel();
    if (probeIt_ && propsDict_.found("suppressProbe"))
        probeIt_=!Switch(propsDict_.lookup("suppressProbe"));
    if(probeIt_)
    {
      forAll(eulerianFieldNames_, fieldIt)
      {
        particleCloud_.probeM().initialize(dictName, dictName + "_" + eulerianFieldNames_[fieldIt] + ".logDat");
        particleCloud_.probeM().vectorFields_.append("Urel");               //first entry must the be the vector to probe
        if(eulerianFieldNames_[fieldIt]==tempFieldName_) //this is the temperature
        {
            particleCloud_.probeM().scalarFields_.append("Rep");            
            particleCloud_.probeM().scalarFields_.append("Nu");                 //other are debug
        }
        else                
            particleCloud_.probeM().scalarFields_.append("Sh"); 
        particleCloud_.probeM().scalarFields_.append("exchangeRate");
        particleCloud_.probeM().writeHeader();
      }
    }

bool Foam::fileFormats::VTKsurfaceFormat<Face>::read
(
    const fileName& filename
)
{
    const bool mustTriangulate = this->isTri();
    this->clear();

    IFstream is(filename);
    if (!is.good())
    {
        FatalErrorInFunction
            << "Cannot read file " << filename
            << exit(FatalError);
    }

    // assume that the groups are not intermixed
    bool sorted = true;


    // Construct dummy time so we have something to create an objectRegistry
    // from
    Time dummyTime
    (
        "dummyRoot",
        "dummyCase",
        "system",
        "constant",
        false           // enableFunctionObjects
    );

    // Make dummy object registry
    objectRegistry obr
    (
        IOobject
        (
            "dummy",
            dummyTime,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        )
    );

    // Read all
    vtkUnstructuredReader reader(obr, is);
    const faceList& faces = reader.faces();

    // Assume all faces in zone0 unless a region field is present
    labelList zones(faces.size(), 0);
    if (reader.cellData().foundObject<scalarIOField>("region"))
    {
        const scalarIOField& region =
            reader.cellData().lookupObject<scalarIOField>
            (
                "region"
            );
        forAll(region, i)
        {
            zones[i] = label(region[i]);
        }
    }

void Foam::fileControl::initialVertices
(
    pointField& pts,
    scalarField& sizes,
    triadField& alignments
) const
{
    Info<< "    Reading points from file     : " << pointsFile_ << endl;

    pointIOField pointsTmp
    (
        IOobject
        (
            pointsFile_,
            runTime_.constant(),
            runTime_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );

    pts.transfer(pointsTmp);

    Info<< "    Reading sizes from file      : " << sizesFile_ << endl;

    scalarIOField sizesTmp
    (
        IOobject
        (
            sizesFile_,
            runTime_.constant(),
            runTime_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );

    sizes.transfer(sizesTmp);

    Info<< "    Reading alignments from file : " << alignmentsFile_ << endl;

    triadIOField alignmentsTmp
    (
        IOobject
        (
            alignmentsFile_,
            runTime_.constant(),
            runTime_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );

    alignments.transfer(alignmentsTmp);

    if ((pts.size() != sizes.size()) || (pts.size() != alignments.size()))
    {
        FatalErrorInFunction
            << "Size of list of points, sizes and alignments do not match:"
            << nl
            << "Points size     = " << pts.size() << nl
            << "Sizes size      = " << sizes.size() << nl
            << "Alignments size = " << alignments.size()
            << abort(FatalError);
    }
}

Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New(const IOobject& io)
{
    wordList libNames(1);
    libNames[0]=word("mesquiteMotionSolver");

    forAll(libNames,i) {
        const word libName("lib"+libNames[i]+".so");

        bool ok=dlLibraryTable::open(libName);
        if(!ok) {
            WarningIn("dynamicFvMesh::New(const IOobject& io)")
                << "Loading of dynamic mesh library " << libName
                    << " unsuccesful. Some dynamic mesh  methods may not be "
                    << " available"
                    << endl;
        }
    }

    // Enclose the creation of the dynamicMesh to ensure it is
    // deleted before the dynamicFvMesh is created otherwise the dictionary
    // is entered in the database twice
    IOdictionary dynamicMeshDict
    (
        IOobject
        (
            "dynamicMeshDict",
            io.time().constant(),
            (io.name() == dynamicFvMesh::defaultRegion ? "" : io.name() ),
            io.db(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );

    word dynamicFvMeshTypeName(dynamicMeshDict.lookup("dynamicFvMesh"));

    Info<< "Selecting dynamicFvMesh " << dynamicFvMeshTypeName << endl;

    dlLibraryTable::open
    (
        dynamicMeshDict,
        "dynamicFvMeshLibs",
        IOobjectConstructorTablePtr_
    );

    if (!IOobjectConstructorTablePtr_)
    {
        FatalErrorIn
        (
            "dynamicFvMesh::New(const IOobject&)"
        )   << "dynamicFvMesh table is empty"
            << exit(FatalError);
    }

    IOobjectConstructorTable::iterator cstrIter =
        IOobjectConstructorTablePtr_->find(dynamicFvMeshTypeName);

    if (cstrIter == IOobjectConstructorTablePtr_->end())
    {
        FatalErrorIn
        (
            "dynamicFvMesh::New(const IOobject&)"
        )   << "Unknown dynamicFvMesh type " << dynamicFvMeshTypeName
            << endl << endl
            << "Valid dynamicFvMesh types are :" << endl
            << IOobjectConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<dynamicFvMesh>(cstrIter()(io));
}

void Foam::calcTypes::fieldMap2d::preCalc
(
  const argList& args,
  const Time& runTime,
  const fvMesh& mesh
)
{
  
  const word dictName("fieldMap2dDict");
  IOdictionary dict
  (
    IOobject
    (
      dictName,
      runTime.system(),
      mesh,
      IOobject::MUST_READ,
      IOobject::NO_WRITE
    )
  );

  if( !dict.readIfPresent<word>("patchName", patchName) )
  {
    SeriousErrorIn("preCalc")
          << "There is no patchName parameter in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  
  if( !dict.readIfPresent<word>("geometry", geometry) )
  {
    SeriousErrorIn("preCalc")
          << "There is no geometry parameter in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  
  point minPoint;
  if( !dict.readIfPresent<point>("minPoint", minPoint) )
  {
    SeriousErrorIn("preCalc")
          << "There is no minPoint parameter in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  point maxPoint;
  if( !dict.readIfPresent<point>("maxPoint", maxPoint) )
  {
    SeriousErrorIn("preCalc")
          << "There is no maxPoint parameter in fieldMap2dDict dictionary"
          << exit(FatalError);
  }

  int integrationDir;
  if( !dict.readIfPresent<int>("integrationDirection", integrationDir) )
  {
    SeriousErrorIn("preCalc")
          << "There is no integrationDirection parameter "
             "in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  intDir = integrationDir;
  
  int flowDir;
  if( !dict.readIfPresent<int>("flowDirection", flowDir) )
  {
    SeriousErrorIn("preCalc")
          << "There is no flowDirection parameter "
             "in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  majDir = flowDir;

  int expectedNumberOfIntersections;
  if
  (
    !dict.readIfPresent<int>
          (
            "expectedNumberOfIntersections", 
            expectedNumberOfIntersections
          ) 
  )
  {
    SeriousErrorIn("preCalc")
          << "There is no expectedNumberOfIntersections parameter "
             "in fieldMap2dDict dictionary"
          << exit(FatalError);
  }
  expNI = expectedNumberOfIntersections;
  
  int numberOfIntegrationPoints;
  if
  (
    !dict.readIfPresent<int>
          (
            "numberOfIntegrationPoints",
            numberOfIntegrationPoints
          ) 
  )
  {
    SeriousErrorIn("preCalc")
          << "There is no numberOfIntegrationPoints parameter "
             "in fieldMap2dDict dictionary"
//.........这里部分代码省略.........

autoPtr<DMDModel<Type> > 
DMDModel<Type>::New
(
   const fvMesh& mesh
)
{
    // get model name, but do not register the dictionary
    // otherwise it is registered in the database twice
    const word DMDModelName
    (
        IOdictionary
        (
            IOobject
            (
                "DMDProperties",
                mesh.time().constant(),
                mesh,
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).lookup("DMDModel")
    );

    Info<< "Selecting DMD model " << DMDModelName << endl;

    if (!dictionaryConstructorTablePtr_)
    {
        FatalErrorIn
        (
            "DMDModel::New"
            "("
                "const fvMesh& mesh,"
                "const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fields"
            ")"
        )   << "DMD model table is empty"
            << exit(FatalError);
    }

    typename dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(DMDModelName);

    if (cstrIter == dictionaryConstructorTablePtr_->end())
    {
        FatalErrorIn
        (
            "DMDModel::New"
            "("
                "const fvMesh& mesh,"
                "const PtrList<GeometricField<Type, fvPatchField, volMesh> >& fields"
            ")"
        )   << "Unknown DMDModel type "
            << DMDModelName << nl << nl
            << "Valid DMDModel types:" << endl
            << dictionaryConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<Foam::DMDModel<Type> >
    (
        cstrIter()(mesh, DMDModelName)
    );
}

Foam::tmp<Foam::surfaceScalarField> Foam::harmonic<Type>::weights
(
    const GeometricField<Type, fvPatchField, volMesh>& phi
) const
{
    // Implement weights-based stabilised harmonic interpolation using
    // magnitude of type
    // Algorithm:
    // 1) calculate magnitude of internal field and neighbour field
    // 2) calculate harmonic mean magnitude
    // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei
    // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei)
    // 5) Use weights to interpolate values

    tmp<surfaceScalarField> tw
    (
        new surfaceScalarField
        (
            IOobject
            (
                "harmonicWeightingFactors" + phi.name(),
                this->mesh().time().timeName(),
                this->mesh()
            ),
            this->mesh() ,
            dimless
        )
    );

    const surfaceScalarField& deltaCoeffs = this->mesh().deltaCoeffs();
    const surfaceScalarField& weights = this->mesh().weights();

    const magLongDelta& mld = magLongDelta::New(this->mesh());
    const surfaceScalarField& longDelta = mld.magDelta();

    surfaceScalarField& w = tw();

    const unallocLabelList& owner = this->mesh().owner();
    const unallocLabelList& neighbour = this->mesh().neighbour();

    scalarField magPhi = mag(phi);

    scalarField& wIn = w.internalField();

    // Larger small for complex arithmetic accuracy
    const scalar kSmall = 1000*SMALL;

    // Calculate internal weights using field magnitude
    forAll (owner, faceI)
    {
        label own = owner[faceI];
        label nei = neighbour[faceI];

        scalar mOwn = magPhi[own]/(1 - weights[faceI]);
        scalar mNei = magPhi[nei]/weights[faceI];

        scalar den = magPhi[nei] - magPhi[own];

        scalar mean = mOwn*mNei/
            ((mOwn + mNei)*longDelta[faceI]*deltaCoeffs[faceI]);

        // Note: complex arithmetic requires extra accuracy
        // This is a division of two close subtractions
        // HJ, 28/Sep/2011
        if (mag(den) > kSmall)
        {
            // Limit weights for round-off safety
            wIn[faceI] =
                Foam::max(0, Foam::min((magPhi[nei] - mean)/den, 1));
        }
        else
        {
            wIn[faceI] = 0.5;
        }
    }

Foam::quadraticFitSnGradData::quadraticFitSnGradData
(
    const fvMesh& mesh,
    const scalar cWeight
)
:
    MeshObject<fvMesh, quadraticFitSnGradData>(mesh),
    centralWeight_(cWeight),
    #ifdef SPHERICAL_GEOMETRY
        dim_(2),
    #else
        dim_(mesh.nGeometricD()),
    #endif
    minSize_
    (
        dim_ == 1 ? 3 :
        dim_ == 2 ? 6 :
        dim_ == 3 ? 9 : 0
    ),
    stencil_(mesh),
    fit_(mesh.nInternalFaces())
{
    if (debug)
    {
        Info<< "Contructing quadraticFitSnGradData" << endl;
    }

    // check input
    if (centralWeight_ < 1 - SMALL)
    {
        FatalErrorIn("quadraticFitSnGradData::quadraticFitSnGradData")
            << "centralWeight requested = " << centralWeight_
            << " should not be less than one"
            << exit(FatalError);
    }

    if (minSize_ == 0)
    {
        FatalErrorIn("quadraticFitSnGradData")
            << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
    }

    // store the polynomial size for each face to write out
    surfaceScalarField snGradPolySize
    (
        IOobject
        (
            "quadraticFitSnGradPolySize",
            "constant",
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("quadraticFitSnGradPolySize", dimless, scalar(0))
    );

    // Get the cell/face centres in stencil order.
    // Centred face stencils no good for triangles of tets. Need bigger stencils
    List<List<point> > stencilPoints(stencil_.stencil().size());
    stencil_.collectData
    (
        mesh.C(),
        stencilPoints
    );

    // find the fit coefficients for every face in the mesh

    for (label faci = 0; faci < mesh.nInternalFaces(); faci++)
    {
        snGradPolySize[faci] = calcFit(stencilPoints[faci], faci);
    }

    if (debug)
    {
        snGradPolySize.write();
        Info<< "quadraticFitSnGradData::quadraticFitSnGradData() :"
            << "Finished constructing polynomialFit data"
            << endl;
    }
}

            );

            // What to do with exposed internal faces if put into this patch?
        }
    }


    // Create the complete field from the pieces
    tmp<GeometricField<Type, fvPatchField, volMesh> > tresF
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "subset"+vf.name(),
                sMesh.time().timeName(),
                sMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            sMesh,
            vf.dimensions(),
            internalField,
            patchFields
        )
    );

    return tresF;
}


template<class Type>

void Foam::scalarTransportPOD::calcDerivativeCoeffs() const
{
    if (derivativeMatrixPtr_)
    {
        FatalErrorIn
        (
            "void scalarTransportPOD::calcDerivativeCoeffs() const"
        )   << "Derivative matrix already calculated"
            << abort(FatalError);
    }

    // Calculate coefficients for differential equation
    // Get times list
    Time& runTime = const_cast<Time&>(this->mesh().time());

    // Remember time index to restore it
    label origTimeIndex = runTime.timeIndex();

    // Read diffusivity

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            this->mesh(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

    Info<< "Reading diffusivity D\n" << endl;

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

    // Find velocity field

    word Uname(this->dict().lookup("velocity"));

    instantList Times = runTime.times();

    volVectorField* Uptr = NULL;

    forAll (Times, i)
    {
        runTime.setTime(Times[i], i);

        Info<< "Time = " << runTime.timeName() << endl;

        IOobject Uheader
        (
            Uname,
            runTime.timeName(),
            this->mesh(),
            IOobject::MUST_READ
        );

        if (Uheader.headerOk())
        {
            Info<< "    Reading " << Uname << endl;

            Uptr = new volVectorField(Uheader, this->mesh());
            break;
        }
        else
        {
            Info<< "    No " << Uname << endl;
        }
    }

Foam::multiSolidBodyMotionFvMesh::multiSolidBodyMotionFvMesh(const IOobject& io)
:
    dynamicFvMesh(io),
    dynamicMeshCoeffs_
    (
        IOdictionary
        (
            IOobject
            (
                "dynamicMeshDict",
                io.time().constant(),
                *this,
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).subDict(typeName + "Coeffs")
    ),
    undisplacedPoints_
    (
        IOobject
        (
            "points",
            io.time().constant(),
            meshSubDir,
            *this,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    )
{
    if (undisplacedPoints_.size() != nPoints())
    {
        FatalIOErrorInFunction
        (
            dynamicMeshCoeffs_
        )   << "Read " << undisplacedPoints_.size()
            << " undisplaced points from " << undisplacedPoints_.objectPath()
            << " but the current mesh has " << nPoints()
            << exit(FatalIOError);
    }


    zoneIDs_.setSize(dynamicMeshCoeffs_.size());
    SBMFs_.setSize(dynamicMeshCoeffs_.size());
    pointIDs_.setSize(dynamicMeshCoeffs_.size());
    label zoneI = 0;

    forAllConstIter(dictionary, dynamicMeshCoeffs_, iter)
    {
        if (iter().isDict())
        {
            zoneIDs_[zoneI] = cellZones().findZoneID(iter().keyword());

            if (zoneIDs_[zoneI] == -1)
            {
                FatalIOErrorInFunction
                (
                    dynamicMeshCoeffs_
                )   << "Cannot find cellZone named " << iter().keyword()
                    << ". Valid zones are " << cellZones().names()
                    << exit(FatalIOError);
            }

            const dictionary& subDict = iter().dict();

            SBMFs_.set
            (
                zoneI,
                solidBodyMotionFunction::New(subDict, io.time())
            );

            // Collect points of cell zone.
            const cellZone& cz = cellZones()[zoneIDs_[zoneI]];

            boolList movePts(nPoints(), false);

            forAll(cz, i)
            {
                label celli = cz[i];
                const cell& c = cells()[celli];
                forAll(c, j)
                {
                    const face& f = faces()[c[j]];
                    forAll(f, k)
                    {
                        label pointi = f[k];
                        movePts[pointi] = true;
                    }
                }
            }

    forAll(databases_, proci)
    {
        meshes_.set
        (
            proci,
            new fvMesh
            (
                IOobject
                (
                    meshName_,
                    databases_[proci].timeName(),
                    databases_[proci]
                )
            )
        );

        pointProcAddressing_.set
        (
            proci,
            new labelIOList
            (
                IOobject
                (
                    "pointProcAddressing",
                    meshes_[proci].facesInstance(),
                    meshes_[proci].meshSubDir,
                    meshes_[proci],
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );

        faceProcAddressing_.set
        (
            proci,
            new labelIOList
            (
                IOobject
                (
                    "faceProcAddressing",
                    meshes_[proci].facesInstance(),
                    meshes_[proci].meshSubDir,
                    meshes_[proci],
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );

        cellProcAddressing_.set
        (
            proci,
            new labelIOList
            (
                IOobject
                (
                    "cellProcAddressing",
                    meshes_[proci].facesInstance(),
                    meshes_[proci].meshSubDir,
                    meshes_[proci],
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );

        boundaryProcAddressing_.set
        (
            proci,
            new labelIOList
            (
                IOobject
                (
                    "boundaryProcAddressing",
                    meshes_[proci].facesInstance(),
                    meshes_[proci].meshSubDir,
                    meshes_[proci],
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );
    }

Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
(
    const fvDOM& dom,
    const fvMesh& mesh,
    const scalar phi,
    const scalar theta,
    const scalar deltaPhi,
    const scalar deltaTheta,
    const label nLambda,
    const absorptionEmissionModel& absorptionEmission,
    const blackBodyEmission& blackBody
)
:
    dom_(dom),
    mesh_(mesh),
    absorptionEmission_(absorptionEmission),
    blackBody_(blackBody),
    I_
    (
        IOobject
        (
            "I" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("I", dimMass/pow3(dimTime), 0.0)
    ),
    Qr_
    (
        IOobject
        (
            "Qr" + name(rayId),
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0)
    ),
    d_(vector::zero),
    dAve_(vector::zero),
    theta_(theta),
    phi_(phi),
    omega_(0.0),
    nLambda_(nLambda),
    ILambda_(nLambda)
{
    scalar sinTheta = Foam::sin(theta);
    scalar cosTheta = Foam::cos(theta);
    scalar sinPhi = Foam::sin(phi);
    scalar cosPhi = Foam::cos(phi);

    omega_ = 2.0*sinTheta*Foam::sin(deltaTheta/2.0)*deltaPhi;
    d_ = vector(sinTheta*sinPhi, sinTheta*cosPhi, cosTheta);
    dAve_ = vector
    (
        sinPhi
       *Foam::sin(0.5*deltaPhi)
       *(deltaTheta - Foam::cos(2.0*theta)
       *Foam::sin(deltaTheta)),
        cosPhi
       *Foam::sin(0.5*deltaPhi)
       *(deltaTheta - Foam::cos(2.0*theta)
       *Foam::sin(deltaTheta)),
        0.5*deltaPhi*Foam::sin(2.0*theta)*Foam::sin(deltaTheta)
    );


    autoPtr<volScalarField> IDefaultPtr;

    forAll(ILambda_, lambdaI)
    {
        IOobject IHeader
        (
            intensityPrefix + "_" + name(rayId) + "_" + name(lambdaI),
            mesh_.time().timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        );

        // check if field exists and can be read
        if (IHeader.headerOk())
        {
            ILambda_.set
            (
                lambdaI,
                new volScalarField(IHeader, mesh_)
            );
        }
        else
        {
            // Demand driven load the IDefault field
            if (!IDefaultPtr.valid())
            {
                IDefaultPtr.reset
//.........这里部分代码省略.........

multivariateSelectionScheme<Type>::multivariateSelectionScheme
(
    const fvMesh& mesh,
    const typename multivariateSurfaceInterpolationScheme<Type>::
        fieldTable& fields,
    const surfaceScalarField& faceFlux,
    Istream& schemeData
)
:
    multivariateSurfaceInterpolationScheme<Type>
    (
        mesh,
        fields,
        faceFlux,
        schemeData
    ),
    schemes_(schemeData),
    faceFlux_(faceFlux),
    weights_
    (
        IOobject
        (
            "multivariateWeights",
            mesh.time().timeName(),
            mesh
        ),
        mesh,
        dimless
    )
{
    typename multivariateSurfaceInterpolationScheme<Type>::
        fieldTable::const_iterator iter = this->fields().begin();

    surfaceScalarField limiter =
    (
        limitedSurfaceInterpolationScheme<Type>::New
        (
            mesh,
            faceFlux_,
            schemes_.lookup(iter()->name())
        )().limiter(*iter())
    );

    for (++iter; iter != this->fields().end(); ++iter)
    {
        limiter = min
        (
            limiter,
            limitedSurfaceInterpolationScheme<Type>::New
            (
                mesh,
                faceFlux_,
                schemes_.lookup(iter()->name())
            )().limiter(*iter())
        );
    }

    weights_ =
        limiter*mesh.surfaceInterpolation::weights()
      + (scalar(1) - limiter)*upwind<Type>(mesh, faceFlux_).weights();
}

Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
:
    IOdictionary
    (
        IOobject
        (
            "fvSchemes",
            obr.time().system(),
            obr,
//             IOobject::MUST_READ,
            IOobject::READ_IF_PRESENT,  // Allow default dictionary creation
            IOobject::NO_WRITE
        )
    ),
    ddtSchemes_
    (
        ITstream
        (
            objectPath() + "::ddtSchemes",
            tokenList()
        )()
    ),
    defaultDdtScheme_
    (
        ddtSchemes_.name() + "::default",
        tokenList()
    ),
    d2dt2Schemes_
    (
        ITstream
        (
            objectPath() + "::d2dt2Schemes",
            tokenList()
        )()
    ),
    defaultD2dt2Scheme_
    (
        d2dt2Schemes_.name() + "::default",
        tokenList()
    ),
    interpolationSchemes_
    (
        ITstream
        (
            objectPath() + "::interpolationSchemes",
            tokenList()
        )()
    ),
    defaultInterpolationScheme_
    (
        interpolationSchemes_.name() + "::default",
        tokenList()
    ),
    divSchemes_
    (
        ITstream
        (
            objectPath() + "::divSchemes",
            tokenList()
        )()
    ),
    defaultDivScheme_
    (
        divSchemes_.name() + "::default",
        tokenList()
    ),
    gradSchemes_
    (
        ITstream
        (
            objectPath() + "::gradSchemes",
            tokenList()
        )()
    ),
    defaultGradScheme_
    (
        gradSchemes_.name() + "::default",
        tokenList()
    ),
    snGradSchemes_
    (
        ITstream
        (
            objectPath() + "::snGradSchemes",
            tokenList()
        )()
    ),
    defaultSnGradScheme_
    (
        snGradSchemes_.name() + "::default",
        tokenList()
    ),
    laplacianSchemes_
    (
        ITstream
        (
            objectPath() + "::laplacianSchemes",
            tokenList()
        )()
    ),
//.........这里部分代码省略.........

void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable()
{
    // Initialise
    if (startSampleTime_ == -1 && endSampleTime_ == -1)
    {
        const polyMesh& pMesh = this->patch().boundaryMesh().mesh()();

        // Read the initial point position
        pointField meshPts;

        if (pMesh.pointsInstance() == pMesh.facesInstance())
        {
            meshPts = pointField(pMesh.points(), this->patch().meshPoints());
        }
        else
        {
            // Load points from facesInstance
            if (debug)
            {
                Info<< "Reloading points0 from " << pMesh.facesInstance()
                    << endl;
            }

            pointIOField points0
            (
                IOobject
                (
                    "points",
                    pMesh.facesInstance(),
                    polyMesh::meshSubDir,
                    pMesh,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE,
                    false
                )
            );
            meshPts = pointField(points0, this->patch().meshPoints());
        }

        pointIOField samplePoints
        (
            IOobject
            (
                "points",
                this->db().time().constant(),
                "boundaryData"/this->patch().name(),
                this->db(),
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE,
                false
            )
        );

        mapperPtr_.reset
        (
            new pointToPointPlanarInterpolation
            (
                samplePoints,
                meshPts,
                perturb_
            )
        );

        // Read the times for which data is available

        const fileName samplePointsFile = samplePoints.filePath();
        const fileName samplePointsDir = samplePointsFile.path();
        sampleTimes_ = Time::findTimes(samplePointsDir);

        if (debug)
        {
            Info<< "timeVaryingMappedFixedValuePointPatchField : In directory "
                << samplePointsDir << " found times "
                << pointToPointPlanarInterpolation::timeNames(sampleTimes_)
                << endl;
        }
    }

    // Find current time in sampleTimes
    label lo = -1;
    label hi = -1;

    bool foundTime = mapperPtr_().findTime
    (
        sampleTimes_,
        startSampleTime_,
        this->db().time().value(),
        lo,
        hi
    );

    if (!foundTime)
    {
        FatalErrorIn
        (
            "timeVaryingMappedFixedValuePointPatchField<Type>::checkTable"
        )   << "Cannot find starting sampling values for current time "
            << this->db().time().value() << nl
            << "Have sampling values for times "
            << pointToPointPlanarInterpolation::timeNames(sampleTimes_) << nl
//.........这里部分代码省略.........