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 following 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 interpolate 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 quantities 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" into 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 containing 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 rigour 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 /**************************************/ |
||
|