CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   How to rotate a pointScalarField? (http://www.cfd-online.com/Forums/openfoam-programming-development/112270-how-rotate-pointscalarfield.html)

carowjp January 25, 2013 01:53

How to rotate a pointScalarField?
 
Hello.

Can anyone please tell me how to rotate a pointScalarField?

I need to keep the mesh fixed and rotate the scalar contents of a field. So for example if you use setFields with boxToCell to set the value of 1 in the box and 0 elsewhere , I would like to then change the orientation of the volScalarField 'box' keeping the mesh fixed.

My idea was to use volPointInterpolation to interpolate cell center values to a pointScalarField, rotate the pointScalarField, and finally interpolate back to the cell centers using pointVolInterpolation.

Thanks for any suggestions on how to accomplish this.

Jim

kmooney January 25, 2013 18:24

Hi Jim,

Just out of curiosity what is it you are trying to do?

One option for rotating the 'box' would be to advect it with an applied rotational velocity field. This would be pretty analogous to the "zalesak's disk" test which is common in multiphase alg. testing. Give it a google and you should see what I mean.

Cheers!

Kyle

carowjp January 25, 2013 18:44

How to rotate a pointScalarField?
 
Kyle,

Thanks for your reply.

I am trying to accomplish rotation of a scalar field at constant rate about a fixed point.

I basically want to calculate the rotated positions of the cell centers and interpolate the values of those grid points back onto the original non-rotated grid.

I would like to do accomplish the result if possible without advecting the field to avoid numerical diffusion, especially since the transformation is constant and known in advance.

Any ideas would be appreciated.

The application I am working on is flow induced by a moving rigid body using a non-body fitted grid and a volume penalization method. The body is represented by a phase fraction and a forcing term is added to the momentum equation.

thanks and regards,

Jim

kmooney January 25, 2013 19:09

Hi Jim,

It sounds like this might be done with some kind of mapping approach.

http://www.openfoam.org/docs/user/mapFields.php

The link shows a console level example but what you would want is something that occurs during runtime. You could probably figure it out from the mapFields source.

Off the top of my head this is what I'm thinking (don't hold me to it):

Start your run with your initial box'd field. Instantiate a second, new mesh that is a copy of the first. Rotate the new mesh. Map the new mesh's field onto the old mesh. Delete the new mesh.
Rinse and repeat.

Just an idea!

Cheers!

Kyle

carowjp January 27, 2013 14:42

Ok, this is as far I have gotten. If there are better ways I would appreciate learning them My doubts are in the copying of a mesh instance and the meshToMesh interpolation.

Code:


    Foam::fvMesh mesh
    (
        Foam::IOobject
        (
            Foam::fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            Foam::IOobject::MUST_READ
        )
    );

    Foam::fvMesh meshR
    (
        Foam::IOobject
        (
            Foam::fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            Foam::IOobject::MUST_READ
        )
    );

    pointField meshPoints(meshR.points());

  // volume fraction of embedded body
    volScalarField alpha
    (
        IOobject
        (
            "alpha",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

  // not sure how to properly 'copy' an instance of a mesh and its fields

    volScalarField alphaR  // Volume fraction of body in mesh to be rotated
    (
        IOobject
        (
            "alphaR",
            runTime.timeName(),
            meshR,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        alpha                                   
    );

    const scalar theta = 45;
    const scalar angle = Foam::cos(theta/2);
    vector axis = vector(0, 0, 1);

    quaternion Rotation(axis,angle);
    tensor R = Rotation.R();

    meshPoints = transform(R, meshPoints);
    meshR.movePoints(meshPoints);            // not sure if this line is necessary

    meshToMesh  meshToMeshInterp(meshR, mesh);

    Info << "interpolate from meshR to mesh." << endl;
    meshToMeshInterp.interpolate(alphaR, alpha, meshToMesh::INTERPOLATE);

Regardless of the order of alpha & alphaR in the arguments I have a crash at the last line of code:

Code:

--> FOAM FATAL ERROR:
the argument field does not correspond to the right mesh. Field size: 529200 mesh size: 529200
From function meshToMesh::interpolateInternalField(Field<Type>& toF, const GeometricField<Type, fvPatchField, volMesh>& fromVf, meshToMesh::order ord) const
in file /home/carowjp/OpenFOAM/OpenFOAM-2.1.1/src/sampling/lnInclude/meshToMeshInterpolate.C at line 129.

This problem is mentioned in the thread: http://www.cfd-online.com/Forums/ope...time-step.html but no solution was posted.

thanks for any tips,

Jim

carowjp January 27, 2013 20:49

how to "Instantiate a second, new mesh that is a copy of the first"?
 
If I change the initial values of the field on the mesh to be rotated as:

Code:

    volScalarField alphaR  // Volume fraction of body in rotated mesh
    (
        IOobject
        (
            "alphaR",
            runTime.timeName(),
            meshR,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        meshR,
        scalar(1)  // was alpha
    );

and use this field order in the arguments:
Code:

meshToMeshInterp.interpolate(alpha, alphaR, meshToMesh::INTERPOLATE);
The code in the previous post will run as a silly test case. :D Although it produces incorrect results if run in parallel. :(

What I am missing though is how to "Instantiate a second, new mesh that is a copy of the first" which also has the values of the field alpha on the original mesh?

thanks,

Jim


All times are GMT -4. The time now is 20:24.