Scalar interpolation


This UDF replaces Solution 586, which was found to have problems in parallel.

Interpolate p, u, v, and w to point data specified in input file. File format is as follows

x1, y1, z1
x2, y2, z2
...
xN, yN, zN

In 2D, no z-coordinate is specified.
The UDF reads line separated data in an input file and interpolates the scalar values of p, u, v, and w to the locations
specified by the coordinates.

First the cell containing the coordinate is located. Then the scalar value at the cell center, y0, is obtained from the solver. Then the gradient of the scalar is gotten. Then by the follwing formula:

y = y0 + dx . Grad(y),

we obtain the value of the scalar at the coordinate. In this formula dx is the vector from the cell centroid to the coordinate.

The UDF follows:
/**************************************/
/* UDF To interpolat P, U, V, and W, using linear extrapolation.

Data points are specified in (x,y) or (x,y,z), line-separated format in an input file default-named "input".

This file is read, and the cells containing each point are identified. If a point is not found, an error
is produced and no results can be generated.

The scalar quantities at the data-point-containing cell-center are obtained from the solver.

The scalar quantites at the data-point are computed from the following:

y = y0 + Grad(y) DOT dx,

where dx is the displacement from the cell center to the data point
and Grad(y) is the gradient of the scalar at the cell center, which is obtained from the solver.

NOTE: You must type "solve set expert" int the FLUENT TUI, and answer "yes" to the question
about whether memory should not be freed.

The UDF consists of two EOD functions:
read_points: reads the data file and finds the cells containig the points. Determines if points cannot be located.
interpolate: computes the value of the scalar using the above formulae for all scalars (p,u,v,w)

Data are written to a file in the following format


POINT X Y [Z] P U V [W]


AUTHOR: William Wangard, Ph.D.
EMAIL: ww@fluent.com

Note: Fluent does not guarantee the rigor of this user-defined function.

Last Modified on June 3, 2003

*/





#include "udf.h"
#include "surf.h"

#define MAXPOINTS 2000

#define NSCALARS (1+ND_ND)

real coords[MAXPOINTS][ND_ND];
int total_count;

#if !RP_HOST

int total_points_found;

struct interpolation_point{
#if PARALLEL
int partition;
int highest_partition;
#endif
cell_t c;
Thread* t;
boolean found;
};

struct interpolation_point point_list[MAXPOINTS];
boolean interpolation_initialized;

/* Function Prototypes */
boolean is_point_in_cell(cell_t c, Thread* t, real* x);

#endif



DEFINE_ON_DEMAND(interpolate)
{

#if !RP_HOST
Thread* t;
cell_t c;
real NV_VEC(dy), NV_VEC(grady), NV_VEC(centroid);
real y0, y;
int point, i, n;
real values[MAXPOINTS][NSCALARS+ND_ND];

#if PARALLEL
real work[(NSCALARS+ND_ND)*MAXPOINTS];
#endif


/* Initialize values */
for(i=0; i<MAXPOINTS; i++)
for(n=0; n<NSCALARS+ND_ND; n++)
values[i][n] = 0.0;

#if PARALLEL
for(i=0; i<(NSCALARS+ND_ND)*MAXPOINTS; i++)
work[i] = 0.0;
#endif

if(!interpolation_initialized)
{
Message0("nnERROR: Problem did not initialize properly.n");
}
else
{


Message0("Computing interpolated scalar data and writing to file...n");


for(point=0; point<total_count; point++)
{

if (point_list[point].found
#if PARALLEL
&& (myid == point_list[point].highest_partition)
#endif
)
{
c = point_list[point].c;
t = point_list[point].t;
C_CENTROID(centroid,c,t);

for(i=0; i<ND_ND; i++)
values[point][i] = coords[point][i];

/* Displacement vector */
y0 = 0.0;
NV_VV(dy, =, coords[point], -, centroid);


for(n=0; n<NSCALARS; n++)
{
switch (n)
{
case 0:
NV_V(grady, =, C_P_G(c,t));
y0 = C_P(c,t);
break;

case 1:
NV_V(grady, =, C_U_G(c,t));
y0 = C_U(c,t);
break;

case 2:
NV_V(grady, =, C_V_G(c,t));
y0 = C_V(c,t);
break;

case 3:
NV_V(grady, =, C_W_G(c,t));
y0 = C_W(c,t);
break;

default:
break;
}

y = y0 + NV_DOT(grady, dy);

values[point][ND_ND + n] = y;
}
}
}

#if PARALLEL
PRF_GRSUM(&values[0][0], (NSCALARS+ND_ND)*MAXPOINTS, work);
#endif


#if PARALLEL
if (I_AM_NODE_ZERO_P)
#endif
{
FILE* output;
output = fopen("interp.out","a");
if (!output)
{
Message("nnERROR: Could not open interpolation output file.n");
return;
}

for(point=0; point<total_points_found; point++)
{
fprintf(output, "%5i: ", point+1);
for(i=0; i< NSCALARS+ND_ND; i++)
{
fprintf(output,"%12.5e ", values[point][i]);
}
fprintf(output,"n");
}
fprintf(output,"n");
fclose(output);
}

}

Message0("Done.nn");

#endif

}




DEFINE_ON_DEMAND(read_points)
{

#if !RP_HOST
Domain* d = Get_Domain(1);
Thread* t;
cell_t c;
int i, point;
#endif


#if !RP_NODE

char* input_file_name = "interp.inp";

FILE* input;
int n, m;

for(n=0; n<MAXPOINTS; n++)
for(m=0; m<ND_ND; m++)
coords[n][m] = 0.0;

n = 0;

input = fopen(input_file_name, "r");

while(!feof(input))
{
#if RP_DOUBLE
#if RP_3D
fscanf(input,"%lg %lg %lgn", &coords[n][0], &coords[n][1], &coords[n][2]);
#else
fscanf(input,"%lg %lgn", &coords[n][0], &coords[n][1]);
#endif
#else
#if RP_3D
fscanf(input,"%g %g %gn", &coords[n][0], &coords[n][1], &coords[n][2]);
#else
fscanf(input,"%g %gn", &coords[n][0], &coords[n][1]);
Message("%g %gn", coords[n][0], coords[n][1]);
#endif
#endif

n++;

if (n == MAXPOINTS)
{
Message("nnWARNING: Number of points in input file has exceeded MAXPOINTS, which is set to %in", MAXPOINTS);
Message(" Recompile UDF with MAXPOINTS >= number of data points in input file.n");
Message(" ... only %i points will be processed...nn", MAXPOINTS);
break;
}


}



total_count = n;

Message("nnThere are %i sets of coordinates read from input file.n",total_count);

fclose(input);

#endif


/* Initialize coordinates on COMPUTE NODES */
host_to_node_int_1(total_count);
host_to_node_real(&coords[0][0],ND_ND*MAXPOINTS);

#if !RP_HOST

interpolation_initialized = FALSE;

for(point=0; point<total_count; point++)
{
point_list[point].found = FALSE;
}


/* Search over all cells */
thread_loop_c(t,d)
{
begin_c_loop_int(c,t)
{

for(point=0; point<total_count; point++)
{
if (point_list[point].found == FALSE)
{
if (is_point_in_cell(c,t,coords[point]))
{
#if PARALLEL
point_list[point].partition = myid;
#endif
point_list[point].c = c;
point_list[point].t = t;
point_list[point].found = TRUE;
}
}
}
}
end_c_loop_int(c,t);
}



total_points_found = 0;


#if PARALLEL
for(point=0; point<total_count; point++)
{
point_list[point].highest_partition = PRF_GIHIGH1(point_list[point].partition);
point_list[point].found = PRF_GIHIGH1(point_list[point].found);
}
#endif

for(point=0; point<total_count; point++)
{
if (point_list[point].found)
{
#if PARALLEL
if (myid == point_list[point].highest_partition)
#endif
total_points_found++;
}
else
{
Message0("Could not find point %i: (", point+1);
for(i=0; i<ND_ND; i++)
Message0(" %12.5e ", coords[point][i]);
Message0(") in the domain.n");
}
}

#if PARALLEL
total_points_found = PRF_GRSUM1(total_points_found);
#endif

if (total_points_found == total_count)
{
Message0("nAll points were successfully located.n");
interpolation_initialized = TRUE;
}
else
{
Message0("ERROR: All points have not been located.n");
interpolation_initialized = FALSE;
}

#if PARALLEL
interpolation_initialized = PRF_GRLOW1(interpolation_initialized);
#endif

#endif

return;

}



#if !RP_HOST

boolean is_point_in_cell(cell_t c, Thread* t, real* x)
{

int i;
face_t f;
Thread* tf;
real A[ND_ND], n[ND_ND], v[ND_ND], face_centroid[ND_ND], z[ND_ND], cell_centroid[ND_ND];
real Amag;

/* Center of cell */
C_CENTROID(cell_centroid,c,t);

/* Loop over all faces of cell */
c_face_loop(c,t,i)
{
/* Face i of cell*/
f = C_FACE(c,t,i);
tf = C_FACE_THREAD(c,t,i);

/* Face normal */
F_AREA(A,f,tf);

/* Probably sufficient to work with A, but normalizing it
will reduce truncation error */
Amag = NV_MAG(A);
NV_VS(n, =, A, /, Amag);

/* Centroid on face i */
F_CENTROID(face_centroid,f,tf);

/* Vector from face centroid to point x */
NV_VV(v, =, x, -, face_centroid);

/* Vector from cell centroid to face centroid */
NV_VV(z, =, face_centroid, -, cell_centroid);

/* Perform test to make sure face normal points outwards */
if (NV_DOT(n,z) < 0.0) NV_S(n,*=,-1.0);

/* If n.v > 0, then point is beyond "plane" that defines the face
and it cannot be inside the cell */
if (NV_DOT(n,v) > 0.0) return FALSE;

}

/* Otherwise it must be in the cell */
return TRUE;
}

#endif
/**************************************/





Show Form
No comments yet. Be the first to add a comment!