atmBoundaryLayer::atmBoundaryLayer(const vectorField& p, const dictionary& dict)
:
    flowDir_(dict.lookup("flowDir")),
    zDir_(dict.lookup("zDir")),
    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
    Uref_(readScalar(dict.lookup("Uref"))),
    Zref_(readScalar(dict.lookup("Zref"))),
    z0_("z0", dict, p.size()),
    zGround_("zGround", dict, p.size()),
    Ustar_(p.size())
{
    if (mag(flowDir_) < SMALL || mag(zDir_) < SMALL)
    {
        FatalErrorInFunction
            << "magnitude of n or z must be greater than zero"
            << abort(FatalError);
    }

    // Ensure direction vectors are normalized
    flowDir_ /= mag(flowDir_);
    zDir_ /= mag(zDir_);

    Ustar_ = kappa_*Uref_/(log((Zref_ + z0_)/z0_));
}

Foam::tmp<Foam::vectorField> Foam::sphericalCS::localToGlobal
(
    const vectorField& local,
    bool translate
) const
{
    const scalarField r(local.component(vector::X));
    const scalarField theta
    (
        local.component(vector::Y)
       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
    );
    const scalarField phi
    (
        local.component(vector::Z)
       *(inDegrees_ ? constant::mathematical::pi/180.0 : 1.0)
    );

    vectorField lc(local.size());
    lc.replace(vector::X, r*cos(theta)*sin(phi));
    lc.replace(vector::Y, r*sin(theta)*sin(phi));
    lc.replace(vector::Z, r*cos(phi));

    return coordinateSystem::localToGlobal(lc, translate);
}

void Foam::dx<Type>::writeDXData
(
    const pointField& points,
    const vectorField& values,
    Ostream& os
) const
{
    // Write data
    os  << "object 3 class array type float rank 1 shape 3 items "
        << values.size()
        << " data follows" << nl;

    forAll(values, elemI)
    {
        os  << float(values[elemI].x()) << ' '
            << float(values[elemI].y()) << ' '
            << float(values[elemI].z()) << nl;
    }

    if (values.size() == points.size())
    {
        os  << nl << "attribute \"dep\" string \"positions\""
            << nl << nl;
    }
    else
    {
        os  << nl << "attribute \"dep\" string \"connections\""
            << nl << nl;
    }
}

// Calculate a simple check on a vector field
void record(vectorField& V)
{
    // Calculate the magnitude of a field
    const int x=0, y=1, z=2;
    double magx = sum(pow2(V[x])) / V.numElements();
    double magy = sum(pow2(V[y])) / V.numElements();
    double magz = sum(pow2(V[z])) / V.numElements();

    cout << "norm = [" << magx
        << " " << magy << " " << magz << " ]" << endl;
}

// Update vectorField for all the new cells. Takes over value of
// original cell.
void Foam::multiDirRefinement::update
(
    const Map<label>& splitMap,
    vectorField& field
)
{
    field.setSize(field.size() + splitMap.size());

    forAllConstIter(Map<label>, splitMap, iter)
    {
        field[iter()] = field[iter.key()];
    }
}

Foam::tmp<Foam::Field<Type> >
Foam::basicSymmetryFvPatchField<Type>::snGradTransformDiag() const
{
    const vectorField nHat(this->patch().nf());

    vectorField diag(nHat.size());

    diag.replace(vector::X, mag(nHat.component(vector::X)));
    diag.replace(vector::Y, mag(nHat.component(vector::Y)));
    diag.replace(vector::Z, mag(nHat.component(vector::Z)));

    return transformFieldMask<Type>(pow<vector, pTraits<Type>::rank>(diag));
}

Foam::tmp<Foam::Field<Type> >
Foam::partialSlipFvPatchField<Type>::snGradTransformDiag() const
{
    const vectorField nHat(this->patch().nf());
    vectorField diag(nHat.size());

    diag.replace(vector::X, mag(nHat.component(vector::X)));
    diag.replace(vector::Y, mag(nHat.component(vector::Y)));
    diag.replace(vector::Z, mag(nHat.component(vector::Z)));

    return
        valueFraction_*pTraits<Type>::one
      + (1.0 - valueFraction_)
       *transformFieldMask<Type>(pow<vector, pTraits<Type>::rank>(diag));
}

void snapshot(vectorField& P)
{
    static int snapshotNum = 0;

    ++snapshotNum;
    char filename[128];
    sprintf(filename, "velocity%03d.m", snapshotNum);

    ofstream ofs(filename);
    int N = P.length(firstDim);

    int k = N/2;

    ofs << "P" << snapshotNum << " = [ ";
    for (int i=0; i < N; ++i)
    {
        for (int j=0; j < N; ++j)
        {
            float value = norm(P(k,j,N-i-1));
            ofs << value << " ";
        }
        if (i < N-1)
            ofs << ";" << endl;
    }
    ofs << "];" << endl;
}

// Update vectorField for all the new cells. Takes over value of
// original cell.
void Foam::multiDirRefinement::update
(
    const Map<label>& splitMap,
    vectorField& field
)
{
    field.setSize(field.size() + splitMap.size());

    for
    (
        Map<label>::const_iterator iter = splitMap.begin();
        iter != splitMap.end();
        ++iter
    )
    {
        field[iter()] = field[iter.key()];
    }
}

Foam::tmp<Foam::vectorField> Foam::transform
(
    const quaternion& q,
    const vectorField& tf
)
{
    tmp<vectorField > tranf(new vectorField(tf.size()));
    transform(tranf.ref(), q, tf);
    return tranf;
}

Foam::tmp<Foam::vectorField> Foam::waveSuperposition::velocity
(
    const scalar t,
    const vectorField& xyz
) const
{
    vectorField result(xyz.size(), vector::zero);

    forAll(waveModels_, wavei)
    {
        const vector2D d(cos(waveAngles_[wavei]), sin(waveAngles_[wavei]));
        const vector2DField xz
        (
            zip
            (
                d & zip(xyz.component(0), xyz.component(1)),
                tmp<scalarField>(xyz.component(2))
            )
        );
        const vector2DField uw
        (
            waveModels_[wavei].velocity(t, xz)
        );
        result += zip
        (
            d.x()*uw.component(0),
            d.y()*uw.component(0),
            uw.component(1)
        );
    }

    tmp<scalarField> s = scale(zip(xyz.component(0), xyz.component(1)));

    return s*result;
}

/*
 * Adjust the time step according to the CFL stability criterion
 */
void adjustTimeStep(vectorField& V)
{
    // Find maximum velocity magnitude
    double maxV = 0.0;

    // NEEDS_WORK: Blitz should provide a norm(vectorField) function.
    // This is ugly.

    for (int i=V.lbound(0); i <= V.ubound(0); ++i)
      for (int j=V.lbound(1); j <= V.ubound(1); ++j)
        for (int k=V.lbound(2); k <= V.ubound(2); ++k)
        {
            double normV = norm(V(i,j,k));
            if (normV > maxV)
                maxV = normV;
        }

    cout << "Maximum velocity is " << maxV << " m/s" << endl;

    maxV += 1e-10;   // Avoid divide-by-zero

    // Steve K: need to have spatialStep^2 / diffusion constant
    // diffusion constant = eta * recip_rho

    delta_t = 0.3 * spatialStep / maxV;
    const double maxTimeStep = 0.01;

    if (delta_t > maxTimeStep)
        delta_t = maxTimeStep;

    cout << "Set time step to " << delta_t << " s" << endl;
}

Foam::tmp<Foam::symmTensorField> Foam::STARCDCoordinateRotation::
transformVector
(
    const vectorField& st
) const
{
    tmp<symmTensorField> tfld(new symmTensorField(st.size()));
    symmTensorField& fld = tfld();

    forAll(fld, i)
    {
        fld[i] = transformPrincipal(R_, st[i]);
    }

Foam::tmp<Foam::vectorField> Foam::toroidalCS::localToGlobal
(
    const vectorField& local,
    bool translate
) const
{
    const scalarField r = local.component(vector::X);

    const scalarField theta =
        local.component(vector::Y)*mathematicalConstant::pi/180.0;

    const scalarField phi =
        local.component(vector::Z)*mathematicalConstant::pi/180.0;

    const scalarField rprime = radius_ + r*sin(phi);

    vectorField lc(local.size());
    lc.replace(vector::X, rprime*cos(theta));
    lc.replace(vector::Y, rprime*sin(theta));
    lc.replace(vector::Z, r*cos(phi));

    return coordinateSystem::localToGlobal(lc, translate);
}

Foam::tmp<Foam::vectorField> Foam::sphericalCS::globalToLocal
(
    const vectorField& global,
    bool translate
) const
{
    const vectorField lc = coordinateSystem::globalToLocal(global, translate);
    const scalarField r = mag(lc);

    tmp<vectorField> tresult(new vectorField(lc.size()));
    vectorField& result = tresult();

    result.replace
    (
        vector::X, r

    );

    result.replace
    (
        vector::Y,
        atan2
        (
            lc.component(vector::Y),
            lc.component(vector::X)
        )*180.0/mathematicalConstant::pi
    );

    result.replace
    (
        vector::Z,
        acos(lc.component(vector::Z)/(r + SMALL))*180.0/mathematicalConstant::pi
    );

    return tresult;
}

void Foam::rotateSearchableSurface::getNormal
(
    const List<pointIndexHit>& info,
    vectorField& normal
) const
{
    vectorField iNormal;
    
    transformationSearchableSurface::getNormal
        (
            info,
            iNormal
        );

    normal.setSize(iNormal.size());

    forAll(normal,i) {
        normal[i]=transform(iNormal[i]);
    }

void Foam::scaleSearchableSurface::getNormal
(
    const List<pointIndexHit>& info,
    vectorField& normal
) const
{
    vectorField iNormal;

    transformationSearchableSurface::getNormal
        (
            info,
            iNormal
        );

    normal.setSize(iNormal.size());

    forAll(normal,i) {
        normal[i]=inverseTransform(iNormal[i]);
        scalar len=mag(normal[i]);
        if(len>SMALL) {
            normal[i]/=len;
        }
    }

Foam::tmp<Foam::vectorField> Foam::parabolicCylindricalCS::localToGlobal
(
    const vectorField& local,
    bool translate
) const
{
    if (min(local.component(vector::Y)) < 0.0)
    {
        FatalErrorIn
        (
            "parabolicCylindricalCS::localToGlobal"
            "(const vectorField&, bool) const"
        )   << "parabolic cylindrical coordinates v < 0"
            << abort(FatalError);
    }

    vectorField lc(local.size());
    lc.replace
    (
        vector::X,
        0.5*
        (
            sqr(local.component(vector::X))
          - sqr(local.component(vector::Y))
        )
    );

    lc.replace
    (
        vector::Y,
        local.component(vector::X) * local.component(vector::Y)
    );

    lc.replace
    (
        vector::Z,
        local.component(vector::Z)
    );

    return coordinateSystem::localToGlobal(lc, translate);
}

void Foam::meshRefinement::findNearest
(
    const labelList& meshFaces,
    List<pointIndexHit>& nearestInfo,
    labelList& nearestSurface,
    labelList& nearestRegion,
    vectorField& nearestNormal
) const
{
    pointField fc(meshFaces.size());
    forAll(meshFaces, i)
    {
        fc[i] = mesh_.faceCentres()[meshFaces[i]];
    }

    const labelList allSurfaces(identity(surfaces_.surfaces().size()));

    surfaces_.findNearest
    (
        allSurfaces,
        fc,
        scalarField(fc.size(), sqr(GREAT)),    // sqr of attraction
        nearestSurface,
        nearestInfo
    );

    // Do normal testing per surface.
    nearestNormal.setSize(nearestInfo.size());
    nearestRegion.setSize(nearestInfo.size());

    forAll(allSurfaces, surfI)
    {
        DynamicList<pointIndexHit> localHits;

        forAll(nearestSurface, i)
        {
            if (nearestSurface[i] == surfI)
            {
                localHits.append(nearestInfo[i]);
            }
        }

        label geomI = surfaces_.surfaces()[surfI];

        pointField localNormals;
        surfaces_.geometry()[geomI].getNormal(localHits, localNormals);

        labelList localRegion;
        surfaces_.geometry()[geomI].getRegion(localHits, localRegion);

        label localI = 0;
        forAll(nearestSurface, i)
        {
            if (nearestSurface[i] == surfI)
            {
                nearestNormal[i] = localNormals[localI];
                nearestRegion[i] = localRegion[localI];
                localI++;
            }
        }
    }

void Foam::dirichletNeumannFriction::correct
(
    const vectorField& slavePressure,
    const PrimitivePatch<face, List, pointField>& masterFaceZonePatch,
    const PrimitivePatch<face, List, pointField>& slaveFaceZonePatch,
    const intersection::algorithm alg,
    const intersection::direction dir,
    const word interpolationMethod,
    const word fieldName,
    const Switch orthotropic,
    const word nonLinear,
    const vectorField& slaveFaceNormals
)
{
    const fvMesh& mesh = mesh_;
    const label slavePatchIndex = slavePatchID();
    const label masterPatchIndex = masterPatchID();
    contactIterNum_++;

    // we have local masterDU and we want to interpolate it to the slave
    // to get local masterDUInterpToSlave (i.e. masterDU interpolated
    // to the slave)
    // so the method is:
    // create global masterDU field
    // interpolate global masterDU from master global face zone to
    // slave global zone
    // then find local masterDUInterpToSlave from the global interpolated field

    vectorField masterDUInterpToSlave
    (mesh.boundaryMesh()[slavePatchIndex].size(), vector::zero);

    // global master DU
    vectorField globalMasterDU(masterFaceZonePatch.size(), vector::zero);

    // lookup current displacement field
    const volVectorField& dispField =
        mesh.objectRegistry::lookupObject<volVectorField>(fieldName);

    // local master and slave DU increment
    vectorField masterDU = dispField.boundaryField()[masterPatchIndex];
    vectorField slaveDU = dispField.boundaryField()[slavePatchIndex];

    if (fieldName == "U")
    {
        // lookup old U
        const volVectorField& dispOldField =
            mesh.objectRegistry::lookupObject<volVectorField>(fieldName+"_0");

        // subtract old U
        masterDU -= dispOldField.boundaryField()[masterPatchIndex];
        slaveDU -= dispOldField.boundaryField()[slavePatchIndex];
    }
    else if (fieldName != "DU")
    {
        FatalError << "dirichletNeumannFriction::correct()\n"
                   " The displacement field must be called U or DU"
                   << exit(FatalError);
    }

    // put local masterDU into globalMasterDU
    const label masterPatchStart
        = mesh.boundaryMesh()[masterPatchIndex].start();
    forAll(masterDU, i)
    {
        globalMasterDU[
            mesh.faceZones()[masterFaceZoneID()
                            ].whichFace(masterPatchStart + i)] =
                                masterDU[i];
    }
    //- exchange parallel data
    // sum because each face is only on one proc
    reduce(globalMasterDU, sumOp<vectorField>());

    // globalMasterDU is interpolated to the slave
    vectorField globalMasterDUInterpToSlave
    (slaveFaceZonePatch.size(), vector::zero);

    // interpolate DU from master to slave using inverse distance or ggi
    if (interpolationMethod == "patchToPatch")
    {
        PatchToPatchInterpolation<
        PrimitivePatch<
        face, List, pointField
        >, PrimitivePatch<face, List, pointField>
        > masterToSlavePatchToPatchInterpolator
        (
            masterFaceZonePatch, // from zone
            slaveFaceZonePatch, // to zone
            alg,
            dir
        );
        globalMasterDUInterpToSlave =
            masterToSlavePatchToPatchInterpolator.faceInterpolate<vector>
            (
                globalMasterDU
            );
    }
    else if (interpolationMethod == "ggi")
    {
        GGIInterpolation<
//.........这里部分代码省略.........