CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

How to get cell indices and manipulate field data

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By ch1

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 23, 2019, 10:17
Default How to get cell indices and manipulate field data
  #1
ch1
New Member
 
Christian
Join Date: Jun 2019
Posts: 13
Rep Power: 2
ch1 is on a distinguished road
Hey everybody,
the following routine written in MATLAB illustrates a basic field manipulation problem, which I am trying to implement in OF:

Code:
 1    %% Initiliaze some example fields
 2    A = rand(4,4,4); %initiliaze a scalar field with random values (mesh1)
 3    B = zeros(2,4,4); %initiliaze an empty scalar field with zeros (mesh2)
 4    x = linspace(0.25,1,4) - 0.125; %get the x-coordinates of the cell centers.
 5    %% The following needs to be implemented in OpenFOAM
 6    idx = x < 0.5; %find the indices of the cells, which are within the desired     range (x < 0.5)
 7    Aidx = A(idx,:,:); %get the corresponding field values of the cells, which are in the desired range
 8    idx2 = Aidx > 0.1; %find the indices of the cells, whose field values are higher than 0.1
 9    B(idx2) = Aidx(idx2) .* 2; %do some operation with these values and write them to the new field in mesh2
10    B(~idx2) = 0; %set the other values of that new field in mesh2 to zero
The first one is similar to the problem of mapping e.g. a volScalarField between two meshes. This is something I was also searching on the forum but I still didn't figure out how to implement it. However, to keep things simple I want to solve my problem differently: Consider two meshes A and B with different dimensions but with the same cell size (so that B could be a sub-mesh of A). I want to "copy" fields between those meshes just by appropriate addressing of the cells, see lines 6 and 7.

After that, the next step would be to do some manipulation with the field data, see lines 8 to 10.

If the MATLAB Code is unclear I can explain in more detail. I just found that it's the most compact way to describe my problem.

Best regards
Christian
ch1 is offline   Reply With Quote

Old   June 24, 2019, 05:00
Default
  #2
ch1
New Member
 
Christian
Join Date: Jun 2019
Posts: 13
Rep Power: 2
ch1 is on a distinguished road
Here is an example of how I am trying to manipulate the field:


Code:
dimensionedScalar TR ("TR", dimensionSet(0,0,0,1,0,0,0), 1000);

forAll (T,idx)
{
    if (T[idx] > TR)
    {
        qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR);
    }
    else
    {
        qR[idx] = 0.0;
    }
}
Which gives me the following compilation error:


qREqn.H: In function int main(int, char**): /home/of/foam/foam-extend-4.0/src/foam/lnInclude/UList.H:410:36: error: overloaded function with no contextual type information for (Foam::label i=0; i<(list).size(); i++) ^ qREqn.H:3:1: note: in expansion of macro forAll forAll (T,idx) ^~~~~~ In file included from rhoEmFoam.C:161:0: qREqn.H:5:14: error: invalid types <unresolved overloaded function type>[Foam::label {aka int}] for array subscript if (T[idx] > TR) ^ qREqn.H:7:34: error: invalid types <unresolved overloaded function type>[Foam::label {aka int}] for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:58: error: invalid types <unresolved overloaded function type>[Foam::label {aka int}] for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:74: error: invalid types <unresolved overloaded function type>[Foam::label {aka int}] for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR);
ch1 is offline   Reply With Quote

Old   June 24, 2019, 08:12
Default
  #3
Super Moderator
 
Tobi's Avatar
 
Tobias Holzmann
Join Date: Oct 2010
Location: Augsburg
Posts: 2,347
Blog Entries: 6
Rep Power: 40
Tobi has a spectacular aura aboutTobi has a spectacular aura about
Send a message via ICQ to Tobi Send a message via Skype™ to Tobi
Hi,

of which kind of class does T belong? I guess it is the temperature field.

Code:
const scalarField& Tcell = T.internalField();

forAll(Tcell, celli)
{

}
By the way. You cannot compare a scalar with a dimensionedScalar which is the error you get. Either you define your limit as a normal scalar or you do extract the value of the dimensionedScalar:
Code:
    if (T[idx] > TR.value())
    {
        qR[idx] = 5600.0 * (T[idx] - TR.value()) + 181.0 * (T[idx] - TR.value()) * (T[idx] - TR.value());
    }
__________________
Keep foaming,
Tobias Holzmann
Tobi is offline   Reply With Quote

Old   June 24, 2019, 12:36
Default
  #4
ch1
New Member
 
Christian
Join Date: Jun 2019
Posts: 13
Rep Power: 2
ch1 is on a distinguished road
Hey Tobi,
thanks for your reply. I made some corrections, this code is working now:


Code:
volScalarField& T = const_cast<volScalarField&> (thermo.T()); //get the temperature field

//const scalar TR = 1000.0;
scalar TR(readScalar(physicalProperties.lookup("TR"))); 

const scalar c1 = 5600.0;
const scalar c2 = 181.0;

forAll(T, idx)
{
    const scalar& cellT = T[idx];
    if (cellT > TR)
    {
        qR[idx] = c1 * (cellT - TR) + c2 * sqr(cellT - TR);
    }
    else
    {
        qR[idx] = 0.0;
    }
}
The first problem was that the temperature field is by default not directly available after the enthalpy Equation is solved (I am using a modified rhoPimpleFoam solver). Now, I found how to access it, see the first line above. Also, I am converting now the cell values of the dimensionedScalar temperature field to scalar values, see the first line inside the loop, in order to avoid the illegal comparison, which you mentioned.

The last task I need to solve is the following: Find all cells of field A in meshA, which are in the range of e.g. (x < xmax), and copy the field values to a field B in meshB. Or, since I know the sizes of the meshes it might be even easier to do:
Get all values of field A but reduced by e.g. 50 cells in x-direction while keeping the number of cellValues in y and z direction. Then copy the result to the fieldB in meshB, which has the matching dimensions.


PS: I assume that you can see my second post because of your moderator status. However, I think it needs to be set to public. At least I cannot see it myself.
ch1 is offline   Reply With Quote

Old   June 26, 2019, 04:24
Default
  #5
ch1
New Member
 
Christian
Join Date: Jun 2019
Posts: 13
Rep Power: 2
ch1 is on a distinguished road
Hello,
I solved my problem by using the meshToMesh.H functions, similar to the code described in runTime mapFields from fineMesh to coarseMesh

However, I had to modify the code because I am using foam-extend 4.0.
Here is the working code:


Code:
// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
    IOobject
    (
        "mapFieldsDict",
        runTime.system(),
        runTime,
        IOobject::MUST_READ,
        IOobject::NO_WRITE,
        false
    )
);

HashTable<word> patchMap;
wordList cuttingPatches;
mapFieldsDict.lookup("patchMap") >> patchMap;
mapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

Info<< "Target mesh size: " << meshTarget.nCells() << endl;

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// specifiy source and target field 
Info<< "    interpolating qJ"      << endl;

volScalarField fieldSource(qJ);

volScalarField fieldTarget
(
    IOobject 
    (
        "qJ",
        runTime.timeName(),//meshTarget.time().timeName(),
        meshTarget,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    meshTarget,
    dimensionedScalar("qJ",dimensionSet(1,-1,-3,0,0,0,0), 0.0)
);

// Interpolate field
meshToMeshInterp.interpolate
(
    fieldTarget,
    fieldSource,
    meshToMesh::INTERPOLATE
);

// Write field
 fieldTarget.write();
Also, make sure that

Code:
#include "meshTools.H"
#include "meshToMesh.H"
is included and linked to the solver.


--------------------------------------------------

Best regards
Christian
frobaux likes this.
ch1 is offline   Reply With Quote

Old   September 25, 2019, 11:17
Default
  #6
Senior Member
 
Join Date: Jan 2012
Posts: 191
Rep Power: 10
itsme_kit is on a distinguished road
// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh();
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();


Hi Christian

I'm wondering how to define meshSource and meshTarget into emMesh and mesh respectively.

Thanks.



Quote:
Originally Posted by ch1 View Post
Hello,
I solved my problem by using the meshToMesh.H functions, similar to the code described in runTime mapFields from fineMesh to coarseMesh

However, I had to modify the code because I am using foam-extend 4.0.
Here is the working code:


Code:
// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
    IOobject
    (
        "mapFieldsDict",
        runTime.system(),
        runTime,
        IOobject::MUST_READ,
        IOobject::NO_WRITE,
        false
    )
);

HashTable<word> patchMap;
wordList cuttingPatches;
mapFieldsDict.lookup("patchMap") >> patchMap;
mapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

Info<< "Target mesh size: " << meshTarget.nCells() << endl;

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// specifiy source and target field 
Info<< "    interpolating qJ"      << endl;

volScalarField fieldSource(qJ);

volScalarField fieldTarget
(
    IOobject 
    (
        "qJ",
        runTime.timeName(),//meshTarget.time().timeName(),
        meshTarget,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    meshTarget,
    dimensionedScalar("qJ",dimensionSet(1,-1,-3,0,0,0,0), 0.0)
);

// Interpolate field
meshToMeshInterp.interpolate
(
    fieldTarget,
    fieldSource,
    meshToMesh::INTERPOLATE
);

// Write field
 fieldTarget.write();
Also, make sure that

Code:
#include "meshTools.H"
#include "meshToMesh.H"
is included and linked to the solver.


--------------------------------------------------

Best regards
Christian
itsme_kit is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 01:02.