C++ vectorField类代码示例

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;
}