CFD Online Discussion Forums - Point interpolation in OF

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)

- - Point interpolation in OF (https://www.cfd-online.com/Forums/openfoam/96450-point-interpolation.html)

Alucard

January 23, 2012 11:44

Point interpolation in OF

Goodmorning
I need to perform some interpolations in Openfoam.

Imagine I've a point P(x,y) ( I know x and y) I need to find at which cell volume this point belongs and,obviously I need he cell center coordinates via the .x() and .y() functions.
All this because I need to know the vector position difference between the point P and the cell center
Xdiff=X(c)-X(P)

(X is a vector)

that I'll use later for some interpolations.

Does some clever function exist in OpenFOAM??
I'm surfing in the *.C,*.H openfoam files and I'm reading the forum since four days without success :(.

wyldckat

January 23, 2012 15:18

Hi Valério,

I'm not sure how this can be used, but "Foam::interpolation" seems to be the core of what you need: http://foam.sourceforge.net/docs/cpp/a00905.html

For even more information: http://openfoamwiki.net/index.php/In...%28OpenFOAM%29

Best regards and good luck!
Bruno

Alucard

January 24, 2012 05:07

Quote:

Originally Posted by wyldckat (Post 340716)

Hi Bruno
thank you very much I'm reading all of it!

I've to be honest for the moment I didn't figure out how to do it (I understand what OF inteprolation functions do but I didn't understand how to pass them my P point yet)
I have a square domain and a regular mesh so as I know how to obtain the center of a cell I could do it writing a "function" by myself (I come from Fortran and I was used to code everything by myself).
But I think in my lab they will continue to use the software I'm developping also after i'll leave that's why I'm try to get a "gereral" and efficient way that could works for more general domain shapes and unstructured meshes.

niklas

January 24, 2012 05:57

check this thread
http://www.cfd-online.com/Forums/ope...rpolation.html

Alucard

January 24, 2012 10:23

Quote:

Originally Posted by niklas (Post 340810)

check this thread
http://www.cfd-online.com/Forums/ope...rpolation.html

hi Niklas thank you
I tried to copy/paste the part:

Code:

   dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

    autoPtr<interpolation<vector> > Uinterp = interpolation<vector>::New(interpolationDict, U);

    vector pos(0,0,0);

    label cellI = mesh.findCell(pos);

    vector Ui = Uinterp->interpolate(pos, cellI);

just to see if it is compiled but no. The first line it's ok but for the second one i have:

In file included from interFoam.C:88:
miha.H: In function ‘int main(int, char**)’:
miha.H:160: error: ‘interpolation’ was not declared in this scope
miha.H:160: error: template argument 1 is invalid
miha.H:160: error: expected unqualified-id before ‘>’ token

If I have well understood:
-the first line look in the dictionary to see what kind of interpolation scheme I've chosen.
-the second one just says that I'll interpolate U and Uinterp is the pointer
-in the third one there is the vector position X
-in the fourth one I find the index (label) of the cell that belongs X
-in the fifth I obtain the interpolated value of the U field in the X point.
So if I just tailor these lines and I copy them in my "main" i have:

Code:

dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

    autoPtr<interpolation<vector> > aplhainterp = interpolation<vector>::New(interpolationDict, alpha1);

    vector pos(0,0,0);

    label cellI = mesh.findCell(pos);

    scalar alphaint = alphainterp->interpolate(pos, cellI);

I was also looking to:
http://www.cfd-online.com/Forums/ope...onsistent.html

interpolationCellPoint<vector> UInt(U);
point p1(0.0550,0.0048,0.0001);
Info << "p1=" << UInt.interpolate(p1,meshPhys.findCell(p1)) << endl;

I tried to but I have compilation errors.:confused:

I'm sorry for the stupid questions but I'm still learning C++ and undesrtand in few months how Openfoam works is quite difficult for me cause there is no a clear programmer's guide about those things and I also have some problems in beeing confident with the C++ syntax.

Thank you!

niklas

January 24, 2012 11:50

you need to include interpolation.H.
add the line below at the top of includes

#include "interpolation.H"

Alucard

January 24, 2012 12:33

Quote:

Originally Posted by niklas (Post 340907)

you need to include interpolation.H.
add the line below at the top of includes

#include "interpolation.H"

Thank you, i was thinking about the same thing "I've to include somehow it or it's already included in some other .H files?" you give me the right answer.

so now this part

Code:

//   dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

autoPtr<interpolation<scalar> > alphaint =interpolation<scalar>::New(interpolationDict,

alpha1);

vector pos(0,0,0);

label cellI = mesh.findCell(pos);

scalar alphai = alphaint->interpolate(pos, cellI);

gives me back at each time step the value of alphai in the point (0,0,0).
It's a good start thank you very much.
Now the second part (i've to develop) is to store several "position vectors" at each time step and for each of the points of them perform the interpolation step.
I hope I could ask all of you some support if I still have some problem!
thank you very much
Sincerely
Valerio

niklas

January 24, 2012 13:34

easiest way is to read them from a dictionary

stick this at the end of createFields.H

Code:

    IOdictionary myDict

    (

        IOobject

        (

            "myOwnDict",

            runTime.constant(),

            mesh,

            IOobject::MUST_READ,

            IOobject::NO_WRITE

        )

    );



    vectorList positions(myDict.lookup("positions"));

then you create the file myOwnDict and put it in the constant directory.
Easiest way is to copy one of the other files in the constant directory and
then delete everything below the FoamFile{...} header

then you just add the list of vectors to the positions, like this

Code:

positions

(

    (0 0 0)

    (1 1 1)

);

done.

Alucard

January 25, 2012 17:40

Quote:

Originally Posted by niklas (Post 340929)

easiest way is to read them from a dictionary
.....
[/code]done.

Hi again
thank you for the idea but in my case points change during time:( I'll show you the part of my code because still have a problem:

Code:

   

forAll (mesh.cellCentres(),cellI)

{ 



    if(alpha1[cellI]<0.95 && alpha1[cellI]>0.5)  

    { 

     double X_p=mesh.cellCentres()[cellI].x();

     double Y_p=mesh.cellCentres()[cellI].y();

 

     double N_x=nHat[cellI].x();

     double N_y=nHat[cellI].y();



     double X_ext=X_p+N_x;

     double Y_ext=Y_p+N_y; 



     double X_int=X_p+N_x;

     double Y_int=Y_p+N_y;

     

     vector pos_int(X_int,Y_int,0.);

     vector pos_ext(X_ext,Y_ext,0.);



     dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

  

     autoPtr<interpolation<scalar>> Cl_int = interpolation<scalar>::New(interpolationDict, Cl);

     label cellJ = mesh.findCell(pos_int);

     scalar Clint = Cl_int->interpolate(pos_int, cellJ);





     autoPtr<interpolation<scalar>> Cl_ext = interpolation<scalar>::New(interpolationDict, Cl);

     label cellK = mesh.findCell(pos_ext);

     scalar Clext = Cl_ext->interpolate(pos_ext, cellK);







    }    

    }

The idea is to obtain for each point selected by the "if" statement, the difference "Clext-Clint". The part of the code I wrote is just a "structure" the way I obtain the Xext and Yint is not exactly like this (you add to the X_p a vector proportional to the local normal, in the code for te moment it's just "n" without the scale factor but it's just a simple test.
Cl is the liquid concentration and it's just a scalar field (like apha1)
So the problem is that if I write something like
"
autoPtr<interpolation<scalar>> Cl_int = interpolation<scalar>::New(interpolationDict, Cl);
label cellJ = mesh.findCell(pos_int);
scalar Clint = Cl_int->interpolate(pos_int, cellJ);
"

for only one point outside the "forAll" loop I compile it if I put (this is the case) In a Loop I've compilation errors.
Do you perhaps know where I fail? I'm so close to the end: if this part works I'm happy cause all the rest of the code is already working :)
Thank you very much!
Valerio

niklas

January 26, 2012 01:52

Why do you recreate the dictionary and interpolator every time.
I think you have misunderstood how the interpolator works.
You only need to create it once and then give the position and cell number to get the interpolated value
of the field used in the second argument of the New-call below.

Code:

   

dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

// interpolator for the alphaField

autoPtr<interpolation<scalar>> alphaInterpolator = interpolation<scalar>::New(interpolationDict, alpha1);



forAll (mesh.cellCentres(),cellI)

{ 

    if(alpha1[cellI]<0.95 && alpha1[cellI]>0.5)  

    {

         //double X_p=mesh.cellCentres()[cellI].x();

         //double Y_p=mesh.cellCentres()[cellI].y();

         vector p = mesh.cellCentres()[celli];



         //double N_x=nHat[cellI].x();

         //double N_y=nHat[cellI].y();

         vector N = nHat[celli];



         //double X_ext=X_p+N_x;

         //double Y_ext=Y_p+N_y; 

         vector Pext = p + N;



         //double X_int=X_p+N_x;

         //double Y_int=Y_p+N_y;

         vector Pint = p + N; // this seems to be the same as Pext?!?!



         vector pos_int(Pint.x(),Pint.y(),0.);

         vector pos_ext(Pext.x(),Pext.y(),0.);



         label cellJ = mesh.findCell(pos_int);

         scalar Clint = alphaInterpolator->interpolate(pos_int, cellJ);



         label cellK = mesh.findCell(pos_ext);

         scalar Clext = alphaInterpolator->interpolate(pos_ext, cellK);



    } 

}

Alucard

January 26, 2012 12:22

Thank you very much it works now!
Your help was really important to me!
As I'll finish all the modifications and the code will work I'll show you the results!

Alucard

March 30, 2012 09:45

Quote:

Originally Posted by niklas (Post 341218)

Code:

   

dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");

// interpolator for the alphaField

autoPtr<interpolation<scalar>> alphaInterpolator = interpolation<scalar>::New(interpolationDict, alpha1);



forAll (mesh.cellCentres(),cellI)

{ 

    if(alpha1[cellI]<0.95 && alpha1[cellI]>0.5)  

    {

         //double X_p=mesh.cellCentres()[cellI].x();

         //double Y_p=mesh.cellCentres()[cellI].y();

         vector p = mesh.cellCentres()[celli];



         //double N_x=nHat[cellI].x();

         //double N_y=nHat[cellI].y();

         vector N = nHat[celli];



         //double X_ext=X_p+N_x;

         //double Y_ext=Y_p+N_y; 

         vector Pext = p + N;



         //double X_int=X_p+N_x;

         //double Y_int=Y_p+N_y;

         vector Pint = p + N; // this seems to be the same as Pext?!?!



         vector pos_int(Pint.x(),Pint.y(),0.);

         vector pos_ext(Pext.x(),Pext.y(),0.);



         label cellJ = mesh.findCell(pos_int);

         scalar Clint = alphaInterpolator->interpolate(pos_int, cellJ);



         label cellK = mesh.findCell(pos_ext);

         scalar Clext = alphaInterpolator->interpolate(pos_ext, cellK);



    } 

}

Goodmorning again to everybody.
I have now a "new" problem unfortunately.
The subroutine I wrote (thanks a lot Niklas for your precious help), actually works fine in serial mode.
As I try to run it in parallel i have a "bug" when, moving from a "p" point located in one subdomain, I obtain a Pext (or Pint) point that lies outside.

I read (Ref1) that mesh.findCell doesn't "connect" different subdomains and i really don't find the "way" to solve this problem cause I've not good ideas about how can I "connect" them in parallel. i tried reading some infos about parallel programming in openFOAM (Ref.2) cause I'm really new at it.
I hope Niklas you can help me (it's you again in ref.2)
Thank you very much.

(Ref.1 http://www.cfd-online.com/Forums/ope...ordinates.html)

(Ref.2
https://mail.google.com/mail/images/cleardot.gif
http://www.cfd-online.com/Forums/ope...allel-run.html)

ooo	December 13, 2013 13:28

Hi,

Do you know how we can interpolate in a reverse way, i.e interpolate(distribute) a vector, to a volVectorField?
e.g consider we have a volVectorField U. Now i have some vectors and i would like to interpolate(distribute) their values to the existed volVectorField.
(I tried to use pointVolInterpolation but it seems that it does not exist in OpenFOAM anymore. The file cannot be found in http://www.openfoam.org/docs/cpp/. Also using it in code, gives me error (there isn"t any pointVolInterpolation.H , also it has not been declared in Interpolation.H as well)

Also do you know the difference between the accuracy of cellPoint and cellPointFace scheme?

Thanks in advance.

wyldckat

December 29, 2013 20:10

Greetings to all!

Quote:

Originally Posted by ooo (Post 466318)

Do you know how we can interpolate in a reverse way, i.e interpolate(distribute) a vector, to a volVectorField?
e.g consider we have a volVectorField U. Now i have some vectors and i would like to interpolate(distribute) their values to the existed volVectorField.

Mmm... isn't opposite operation named "reconstruct" in OpenFOAM? You can check the following project for ideas on how this is done: https://github.com/wyldckat/reconstr...te-fields/wiki

Quote:

Originally Posted by ooo (Post 466318)

(I tried to use pointVolInterpolation but it seems that it does not exist in OpenFOAM anymore. The file cannot be found in http://www.openfoam.org/docs/cpp/. Also using it in code, gives me error (there isn"t any pointVolInterpolation.H , also it has not been declared in Interpolation.H as well)

Answered here: http://www.cfd-online.com/Forums/ope...tml#post468049 post #5 - the quick answer:

Quote:

Originally Posted by wyldckat (Post 468049)

The correct class name should be "volPointInterpolation".

Quote:

Originally Posted by ooo (Post 466318)

Also do you know the difference between the accuracy of cellPoint and cellPointFace scheme?

Well... I think it all really depends on where you want to interpolate and what field.
Because the "cellPointFace" scheme will take into account the values at the centre of the cell, the nearest points and the nearest centres of cell faces; the problem is whether the values at the face centres are already interpolated results or are they values in the surface mesh (where there is no other cell on the other side of the face). In other words, it really depends on whether there are already too many interpolated values in the neighbourhood.

Best regards,
Bruno

praveensrikanth91

October 5, 2016 13:40

cellPoint interpolation at boundaries

Hello,

I am trying to define a new boundary condition for a fluid-fluid interface in chtMultiRegionFoam to which I have added phase change. I want to get the temperature at the interface and to do that I want to interpolate the temperatures from the cell center values to the interface from both sides and average them. I know that the cellPointFace interpolation uses prescribed boundary values at these faces which turns out to be a default value. I was wondering How does the cellPoint interpolation work at boundaries if I give it the face center of the boundary faces as the position vector? I tried to go through the code for clues but I do not see any specific treatment at boundaries. Does that mean that it would interpolate values at the boundary?

Any information would be appreciated

Thanks
Praveen

All times are GMT -4. The time now is 11:04.