Fluent UDF中沿指定方向获取邻接网格

有朋友在VC++UDF Studio插件群里问Fluent UDF有没有宏可以沿着指定方向搜索到下一个邻接的网格？答案是No，没有现成的宏可以实现这个目的，但是我们可以写一个函数来手动实现。这个思路如下：对于当前的网格，我们可以利用c_face_loop宏配合C_FACE及C_FACE_THREAD宏，对当前网格的所有face进行循环，例如下图中的当前网格循环会有0,1,2,3四个face。然后我们再利用F_C0分别取这四个face的c0，将c0与当前网格比较，如果c0就是当前网格，那么邻接网格就是c1，反之，如果c0不是当前网格，那么邻接网格就是c0。这样，我们就能得到四个邻接网格，至于需要指定方向的邻接网格，可以将当前网格与这四个邻接网格的中心连成矢量，然后分别求这四个矢量与指定的方向矢量的夹角，哪个夹角最小就是我们需要的邻接网格。如果夹角对应的face恰好为边界面则没有邻接网格，函数返回c=-1，表示没有找到邻接网格。

 例子代码如下（VC++ UDF Studio 2022R1学术版上编译通过）

#include "udf.h"

extern "C"
{
	#include "sg.h" //INTERIOR_FACE_GEOMETRY宏和BOUNDARY_FACE_GEOMETRY宏需要此头文件
};

cat_t GetNextCellByCell(cat_t cur_cat,real NV_VEC(direction))  
{
	int i;
	real angle_min=180;  //较大的初始值
	real ds,A_by_es,NV_VEC(A),NV_VEC(es),NV_VEC(dr0),NV_VEC(dr1),NV_VEC(vector);
	cat_t NextCellAndThread={-1, NULL}; //用于存贮最终邻接网格

	c_face_loop(cur_cat.c,cur_cat.t,i) //对当前cell的所有face循环
	{
		face_t f=C_FACE(cur_cat.c,cur_cat.t,i);  //face的f
		Thread* tf=C_FACE_THREAD(cur_cat.c,cur_cat.t,i);  //face的Thread
		if(!BOUNDARY_FACE_THREAD_P(tf)) //如果是内部face
		{
			cell_t c0=F_C0(f,tf); //取face的c0
			Thread*t0 = THREAD_T0(tf);
			INTERIOR_FACE_GEOMETRY(f,tf,A,ds,es,A_by_es,dr0,dr1);
			if(cur_cat.c==c0)  //当前网格是face的c0，那么邻接网格就是c1
			{
				NV_V(vector, =, es); //vector为c0中心到c1中心的单位矢量
				real cos_angle=NV_DOT(direction,vector)/(NV_MAG(direction)*NV_MAG(vector)); //点积公式获取两个矢量夹角的余弦值
				cos_angle=MIN(1.0,MAX(-1.0,cos_angle)); //将范围限制在[-1,1]，以免acos函数出错
				real angle=acos(cos_angle)/M_PI*180;  // 根据余弦值计算角度
				if(angle<angle_min)  //如果小于其它夹角，则就是我们要的cell
				{
					angle_min=angle;
					NextCellAndThread.c=F_C1(f,tf);
					NextCellAndThread.t=THREAD_T1(tf);
				}
			}
			else  //当前网格是face的c1，那么邻接网格就是c0
			{
				NV_VS(vector, =, es,*,-1); //vector为c1中心到c0中心的单位矢量
				real cos_angle=NV_DOT(direction,vector)/(NV_MAG(direction)*NV_MAG(vector)); //点积公式获取两个矢量夹角的余弦值
				cos_angle=MIN(1.0,MAX(-1.0,cos_angle));  //将范围限制在[-1,1]，以免acos函数出错
				real angle=acos(cos_angle)/M_PI*180; // 根据余弦值计算角度
				if(angle<angle_min)  //如果小于其它夹角，则就是我们要的cell
				{
					angle_min=angle;
					NextCellAndThread.c=F_C0(f,tf);
					NextCellAndThread.t=THREAD_T0(tf);
				}
			}
		}
		else //如果是边界face
		{			
			cell_t c0=F_C0(f,tf); //取face的c0
			Thread*t0 = THREAD_T0(tf);
			BOUNDARY_FACE_GEOMETRY(f,tf,A,ds,es,A_by_es,dr0);
			NV_V(vector, =, es);  //vector为c0中心到边界face的单位矢量
			real cos_angle=NV_DOT(direction,vector)/(NV_MAG(direction)*NV_MAG(vector)); //点积公式获取两个矢量夹角的余弦值
			cos_angle=MIN(1.0,MAX(-1.0,cos_angle)); //将范围限制在[-1,1]，以免acos函数出错
			real angle=acos(cos_angle)/M_PI*180; // 根据余弦值计算角度
			if(angle<angle_min)  //如果小于其它夹角，则没有下一个邻接网格
			{
				angle_min=angle;
				NextCellAndThread.c=-1;
				NextCellAndThread.t=NULL;
			}
		}
	}

	return NextCellAndThread;
}

DEFINE_SOURCE(test_source, c,t, dS, eqn)
{
	cat_t current_cat={c,t};
	real NV_VEC(search_dir);

	NV_D(search_dir, =, 1.0, 0, 0); //假设要搜索x正方向的邻接网格
	cat_t NextCellAndThread=GetNextCellByCell(current_cat,search_dir);

	dS[eqn] = 0.0;
	return 0.0;  //只是为了说明GetNextCellByCell函数用法，所以不加源项
}

以上GetNextCellByCell函数即实现了针对cur_cat网格沿着指定矢量方向direction搜索邻接网格的功能。其中cat_t为Fluent内置类型，定义如下，cat就是cell and thread的缩写，实际就是cell的序号和Thread合成一个结构体类型。

typedef struct cell_and_thread_struct
{
  cell_t c;
  struct thread_struct *t;
} cat_t;

注意：以上代码目前仅考虑了真串行，如果是假串行或者并行，还要考虑并行架构，以及跨计算节点去搜索邻接网格，更加复杂一些。

说到真串行假串行问题，这里补充一个小要点，现在新版本的串行都是1个并行节点实现的假串行，UDF编程时候需要考虑并行架构。如果要不考虑并行架构，必须利用VC++ UDF Studio插件中的串行启动来实现真串行。具体可以参考另外一篇博文《Fluent真串行假串行的问题 (学UDF前须知)》。

Fluent真串行假串行的问题 (学UDF前须知)_SuperUDF的博客-CSDN博客