[OpenFOAM] 张量的操作

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

1. tensor2D.H

#include "tensor2D.H"
#include "IOstreams.H"

using namespace Foam;

int main()
{   
    /* Define a tenor in the script */

    // define a 2*2 tensor by dot product of two vectors
    vector2D v1(1, 2), v2(3, 4);
    tensor2D t = v1*v2;

    // define a 2*2 tensor by specifying its components
    tensor2D t1(1, 2, 3, 4);

    // the inverse of a 2*2 tensor
    tensor2D t2 = inv(t1);

    Info<< "v1(1, 2)*v2(3, 4) = " << t << endl;

    Info<< "t1 = " << t1 << endl;

    Info << "t.x() = " << t.x() << endl;
    Info << "t.y() = " << t.y() << endl;

    Info<< "The inverse of tensor t1 is: t2 = " << t2 <<endl;

    return 0;
}

./Make/files:

Test-tensor2D.C
EXE = $(FOAM_USER_APPBIN)/Test-tensor2D

./Make/options：空白文件即可。

2. tensor.H

#include "fvCFD.H"
#include "tensor.H"
#include "symmTensor.H"
#include "transform.H"
#include "stringList.H"
#include "IOstreams.H"

using namespace Foam;

int main()
{
    tensor t1(1, 2, 3, 
              4, 5, 6, 
              7, 8, 9);

    tensor t2(1, 2, 3, 
              1, 2, 3, 
              1, 2, 3);

    tensor t3 = t1 + t2;

    Info<< "t1 = " << t1 << endl;
    Info<< "t2 = " << t2 << endl;

    Info<< "t1 + t2 = " << t3 << endl;

    Info<< "t1.x() = " << t1.x() << endl;
    Info<< "t1.y() = " << t1.y() << endl;
    Info<< "t1.z() = " << t1.z() << endl;

    Info << "------------------------------------------------------" << endl;

    // Inverse of a tensor
    tensor t4(3,-2,1,-2,2,0,1, 0, 4);

    Info << "t4 = " << t4 << endl;
    Info << "Inverse of t4 = " << inv(t4) << endl;
    Info << "Dot product t4 & inv(t4) = " << (inv(t4) & t4) << endl; // should be identity, & is dot product
    Info << "t4(0,0) = " << t4(0,0) << endl;  // slice of tensor
    Info << "------------------------------------------------------" << endl;




    tensor t6(1,0,-4,0,5,4,-4,4,3);
    // tensor t6(1,2,0,2,5,0,0,0,0);

    Info<< "tensor " << t6 << endl;
    vector e = eigenValues(t6);
    Info<< "eigenvalues " << e << endl;

    tensor ev = eigenVectors(t6);
    Info<< "eigenvectors " << ev << endl;

    // The determinant of a tensor is the product of its eigenvalues
    Info<< "Check determinant " << e.x()*e.y()*e.z() << " " << det(t6) << endl;

    Info<< "Check eigenvectors "
        << (eigenVectors(t6, e) & t6) << " "
        << (e.x()*eigenVectors(t6, e).x())
        << (e.y()*eigenVectors(t6, e).y())
        << (e.z()*eigenVectors(t6, e).z())
        << endl;

    Info << "------------------------------------------------------" << endl;

    Info<< "Check eigenvalues for symmTensor "
        << eigenValues(symm(t6)) - eigenValues(tensor(symm(t6))) << endl;

    Info<< "Check eigenvectors for symmTensor "
        << eigenVectors(symm(t6)) - eigenVectors(tensor(symm(t6))) << endl;

    tensor t7(1, 2, 3, 2, 4, 5, 3, 5, 6);

    Info<< "Check transformation "
        << (t1 & t7 & t1.T()) << " " << transform(t1, t7) << endl;

    symmTensor st1(1, 2, 3, 4, 5, 6);
    symmTensor st2(7, 8, 9, 10, 11, 12);

    Info << "symmTensor st1 = " << st1 << endl;

    t1 = tensor(st1);

    Info << "tensor t1 = " << t1 << endl;

    Info << "------------------------------------------------------" << endl;

    Info<< "Check symmetric transformation "
        << transform(t1, st1) << endl;

    Info<< "Check dot product of symmetric tensors "
        << (st1 & st2) << endl;

    Info<< "Check inner sqr of a symmetric tensor "
        << innerSqr(st1) << " " << innerSqr(st1) - (st1 & st1) << endl;

    Info<< "Check symmetric part of dot product of symmetric tensors "
        << twoSymm(st1&st2) - ((st1&st2) + (st2&st1)) << endl;

    tensor sk1 = skew(t6);
    Info<< "Check dot product of symmetric and skew tensors "
        << twoSymm(st1&sk1) - ((st1&sk1) - (sk1&st1)) << endl;

    vector v1(1, 2, 3);

    Info<< sqr(v1) << endl;
    Info<< symm(t7) << endl;
    Info<< twoSymm(t7) << endl;
    Info<< magSqr(st1) << endl;
    Info<< mag(st1) << endl;

    Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl;
    Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl;


    // Intrinsic permeability       
    Info << nl << "Reading permeability field K" << endl;
    volTensorField K
    (
        IOobject
        (
            "K",
            "./",
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    // interpolated permeability
    surfaceTensorField Kf(fvc::interpolate(K,"K"));


    return 0;
}

3. Read a tensor field from file

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    // Intrinsic permeability       
    Info << nl << "Reading permeability field K" << endl;

    volTensorField K
    (
        IOobject
        (
            "K",
            "runTime.timeName()",
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    // interpolated permeability
    surfaceTensorField Kf(fvc::interpolate(K,"K"));

    Info << "Kf = " << Kf << endl;

    forAll(Kf, cellI)
    {
        Info << "Kf[" << cellI << "] = " << Kf[cellI] << endl;
    }

    return 0;
}

第三部分代码的测试文献：case