CFD Online Discussion Forums - Interpolate U.internalField to arbitrary distributed points

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

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Interpolate U.internalField to arbitrary distributed points (https://www.cfd-online.com/Forums/openfoam-programming-development/176865-interpolate-u-internalfield-arbitrary-distributed-points.html)

Interpolate U.internalField to arbitrary distributed points

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.

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

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.

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.

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

I was not yet aware of that, thanks for the tip!

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];

    }

}