CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   Point interpolation in OF (http://www.cfd-online.com/Forums/openfoam/96450-point-interpolation.html)

Alucard January 23, 2012 12: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 16: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 06:07

Quote:

Originally Posted by wyldckat (Post 340716)
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

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 06:57

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

Alucard January 24, 2012 11:23

Quote:

Originally Posted by niklas (Post 340810)

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 12:50

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

#include "interpolation.H"

Alucard January 24, 2012 13: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 14: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 18: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 02: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 13: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)
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);

    }
}


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 14: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 21: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


All times are GMT -4. The time now is 16:23.