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

Interpolate U.internalField to arbitrary distributed points

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

Like Tree2Likes
  • 1 Post By mtgoncalves
  • 1 Post By mtgoncalves

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 29, 2016, 10:55
Default Interpolate U.internalField to arbitrary distributed points
  #1
New Member
 
Join Date: Aug 2016
Posts: 19
Rep Power: 10
mtgoncalves is on a distinguished road
Given a set of n arbitrary distributed points stored in an array:

Code:
double xyz[3*n] = {x1, y1, z1, ..., xn, yn, zn};
How can I interpolate U.internalField to the points contained in xyz? I want to do this purely programmatically without relying on configuration files like sampleDict.
mtgoncalves is offline   Reply With Quote

Old   September 6, 2016, 04:48
Default
  #2
Member
 
David GISEN
Join Date: Jul 2009
Location: Germany
Posts: 68
Rep Power: 17
David* is on a distinguished road
Hi,

some code snippets to start with:

Code:
#include "volFields.H"          // Diverse fields called by fvCFD.H
#include "fvCFD.H"              // OpenFOAM Finite Volume tools - location in makefile
#include "interpolationCellPoint.H"   // Interpolate within a cell using triangulation

interpolationCellPoint<vector> triangulateCellsU(U);
Foam::point FSO_position(xFSO, yFSO, zFSO);
Foam::label FSO_cell = mesh.findCell(FSO_position);
UatSP = Foam::vector::zero;
UatSP = triangulateCellsU.interpolate(FSO_position, FSO_cell);
(FSO is arbitrary)

Cheers, David
David* is offline   Reply With Quote

Old   November 4, 2016, 07:01
Default
  #3
New Member
 
Join Date: Aug 2016
Posts: 19
Rep Power: 10
mtgoncalves is on a distinguished road
Hi David*,

Sorry for my very late reply. I have been quite busy.

Thank you for your suggestion. I gave it a try this week. It compiles fine, however, I am getting a "Segmentation Violation" error. This can be related to either invalid memory access or using too much memory. Do you know how memory-intensive this method is? I am interpolating from 116k cells to 500k points.
mtgoncalves is offline   Reply With Quote

Old   November 9, 2016, 11:28
Default
  #4
New Member
 
Join Date: Aug 2016
Posts: 19
Rep Power: 10
mtgoncalves is on a distinguished road
It seems that I was interpolating to points which were located outside the source domain. By resolving that I am no longer getting the error. However, the result does include some nan's

How can this be possible when the source volVectorField does not include any nan's?

EDIT:
Never mind there are nan's in the source field.
liquidspoon likes this.
mtgoncalves is offline   Reply With Quote

Old   November 10, 2016, 03:23
Default
  #5
Member
 
David GISEN
Join Date: Jul 2009
Location: Germany
Posts: 68
Rep Power: 17
David* is on a distinguished road
Good to hear you figured it out! You probably know already how to check if a point is within the mesh:
Code:
(continued from above)
if (FSO_cell == -1)   // point out of mesh
{ ... }
Cheers, David
David* is offline   Reply With Quote

Old   November 10, 2016, 04:50
Default
  #6
New Member
 
Join Date: Aug 2016
Posts: 19
Rep Power: 10
mtgoncalves is on a distinguished road
I was not yet aware of that, thanks for the tip!
mtgoncalves is offline   Reply With Quote

Old   November 10, 2016, 09:26
Default
  #7
New Member
 
Join Date: Aug 2016
Posts: 19
Rep Power: 10
mtgoncalves is on a distinguished road
So now I have my code running correctly. It is, however, a bit slow because I interpolate to a lot of points. Maybe you guys have suggestions to make this more efficient.

I converted pisoFoam into an object-oriented solver. One of the methods of this PisoFoamSolver class is used to interpolate the vorticity field to arbitrary points (the number of points exceed 100k). Please see the code below I have come up with thanks to David*.

NOTE: Because I am running 2D simulations I am only interested in the z-component of the vorticity field.

Code:
void PisoFoamSolver::interpolate_velocity(double *xProbes,
                                          double *yProbes,
                                          double *zProbes,
                                          double *wProbes,
                                          int numProbes)
{
    // The attributes of class PisoFoamSolver are not objects.
    // The attributes are defined as pointers to the required objects.
    // Therefore 'mesh' and 'U' are brought into the scope of this method
    // by binding references to dereferenced pointers (i.e. dereferenced class attributes).
    Foam::fvMesh& mesh = *meshPtr;
    volVectorField& vorticity = *vorticityPtr;

    // Declarations
    Foam::point probeXYZ(Foam::vector::zero);
    Foam::label probeCell = 0;
    Foam::vector probeW(Foam::vector::zero);
    interpolationCellPoint<vector> interpolationOperator(vorticity);
    int i = 0;

    for (i = 0; i < numProbes; i++)
    {
        // Define coordinate of probe
        probeXYZ[0] = *(xProbes+i);
        probeXYZ[1] = *(yProbes+i);
        probeXYZ[2] = *(zProbes+i);

        // Find cell containing probe
        probeCell = mesh.findCell(probeXYZ);

        // Compute target vorticity
        probeW = interpolationOperator.interpolate(probeXYZ,
                                                   probeCell);

        // Store z-component of vorticity
        *(wProbes+i) = probeW[2];
    }
}
rajibroy likes this.
mtgoncalves is offline   Reply With Quote

Reply

Tags
interpolate, probe, sample

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


Similar Threads
Thread Thread Starter Forum Replies Last Post
[snappyHexMesh] No layers in a small gap bobburnquist OpenFOAM Meshing & Mesh Conversion 6 August 26, 2015 09:38
Interpolate face velocitys to points Tobi OpenFOAM Programming & Development 1 December 16, 2014 08:19
[snappyHexMesh] crash sHM H25E OpenFOAM Meshing & Mesh Conversion 11 November 10, 2014 11:27
[DesignModeler] DM's JScript: FPoint()'s GetPoint(i) function picks points backwards? ANT ANSYS Meshing & Geometry 2 July 23, 2012 15:25
[blockMesh] BlockMeshmergePatchPairs hjasak OpenFOAM Meshing & Mesh Conversion 11 August 15, 2008 07:36


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