
[Sponsors] 
January 23, 2012, 12:44 
Point interpolation in OF

#1 
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
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 . 

January 23, 2012, 16:18 

#2 
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,558
Blog Entries: 39
Rep Power: 97 
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
__________________


January 24, 2012, 06:07 

#3  
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
Quote:
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. 

January 24, 2012, 06:57 

#4 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21 
check this thread
http://www.cfdonline.com/Forums/ope...rpolation.html 

January 24, 2012, 11:23 

#5  
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
Quote:
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); 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 unqualifiedid 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); 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! 

January 24, 2012, 12:50 

#6 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21 
you need to include interpolation.H.
add the line below at the top of includes #include "interpolation.H" 

January 24, 2012, 13:33 

#7  
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
Quote:
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); 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 

January 24, 2012, 14:34 

#8 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21 
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")); 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) ); 

January 25, 2012, 18:40 

#9 
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
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); } } 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 

January 26, 2012, 02:52 

#10 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21 
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 Newcall 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); } } 

January 26, 2012, 13:22 

#11 
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
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! 

March 30, 2012, 09:45 

#12  
Member
valerio
Join Date: Jun 2009
Location: Nancy
Posts: 30
Rep Power: 9 
Quote:
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) 

December 13, 2013, 14:28 

#13 
Member
Join Date: Feb 2012
Posts: 49
Rep Power: 6 
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. 

December 29, 2013, 21:10 

#14  
Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 9,558
Blog Entries: 39
Rep Power: 97 
Greetings to all!
Quote:
Quote:
Quote:
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
__________________


Thread Tools  
Display Modes  


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 
blockstructured mesh for tjunction  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 