CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Point interpolation in OF

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

Like Tree3Likes
  • 1 Post By wyldckat
  • 1 Post By niklas
  • 1 Post By wyldckat

Reply
 
LinkBack Thread Tools Display Modes
Old   January 23, 2012, 12:44
Default Point interpolation in OF
  #1
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
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 .
Alucard is offline   Reply With Quote

Old   January 23, 2012, 16:18
Default
  #2
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,301
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
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
elvis likes this.
wyldckat is offline   Reply With Quote

Old   January 24, 2012, 06:07
Default
  #3
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
Quote:
Originally Posted by wyldckat View Post
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.
Alucard is offline   Reply With Quote

Old   January 24, 2012, 06:57
Default
  #4
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
check this thread
http://www.cfd-online.com/Forums/ope...rpolation.html
niklas is offline   Reply With Quote

Old   January 24, 2012, 11:23
Default
  #5
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
Quote:
Originally Posted by niklas View Post
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:
interpolationCellPoint not consistent?


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.

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!
Alucard is offline   Reply With Quote

Old   January 24, 2012, 12:50
Default
  #6
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
you need to include interpolation.H.
add the line below at the top of includes

#include "interpolation.H"
niklas is offline   Reply With Quote

Old   January 24, 2012, 13:33
Default
  #7
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
Quote:
Originally Posted by niklas View Post
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
Alucard is offline   Reply With Quote

Old   January 24, 2012, 14:34
Default
  #8
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
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.
niklas is offline   Reply With Quote

Old   January 25, 2012, 18:40
Default
  #9
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
Quote:
Originally Posted by niklas View Post
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
Alucard is offline   Reply With Quote

Old   January 26, 2012, 02:52
Default
  #10
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
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);

    } 
}
Tobi likes this.
niklas is offline   Reply With Quote

Old   January 26, 2012, 13:22
Default
  #11
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
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 is offline   Reply With Quote

Old   March 30, 2012, 09:45
Default
  #12
Member
 
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 8
Alucard is on a distinguished road
Quote:
Originally Posted by niklas View Post
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 Extracting cell index from x y z coordinates )

(Ref.2

Looping over a volScalarField in a parallel run)
Alucard is offline   Reply With Quote

Old   December 13, 2013, 14:28
Default
  #13
ooo
Member
 
Join Date: Feb 2012
Posts: 49
Rep Power: 5
ooo is on a distinguished road
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.

Last edited by ooo; December 15, 2013 at 13:59.
ooo is offline   Reply With Quote

Old   December 29, 2013, 21:10
Default
  #14
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 8,301
Blog Entries: 34
Rep Power: 84
wyldckat is just really nicewyldckat is just really nicewyldckat is just really nicewyldckat is just really nice
Greetings to all!

Quote:
Originally Posted by ooo View Post
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 View Post
(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: pointVolInterpolation exists anymore? post #5 - the quick answer:
Quote:
Originally Posted by wyldckat View Post
The correct class name should be "volPointInterpolation".
Quote:
Originally Posted by ooo View Post
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
ooo likes this.
wyldckat is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Face and point ordering during interpolation Arnoldinho OpenFOAM Programming & Development 1 November 29, 2011 05:54
block-structured mesh for t-junction Robert@cfd ANSYS Meshing & Geometry 20 November 11, 2011 05:59
Problem with UDF compiling for kTkLW model Wantami FLUENT 0 July 18, 2011 05:11
Gmsh and samplesurface touf Open Source Meshers: Gmsh, Netgen, CGNS, ... 2 December 10, 2007 03:27
CFX4.3 -build analysis form Chie Min CFX 5 July 12, 2001 23:19


All times are GMT -4. The time now is 18:35.