Fluent UDF中调用Matlab矩阵运算函数(以二维插值为例)

Fluent UDF中经常需要用到一些常见算法，例如插值、拟合、矩阵运算等等，这些在UDF中是没有现成函数实现的，理论上需要我们自己去写函数。另一方面我们又注意到这些运算恰恰是Matlab的强项，几乎调用一个现成的函数就完成了目的。所以我们有什么办法把Matlab函数给UDF直接调用呢？

这里提供一种借助于VC++ UDF Studio插件实现调用Matlab函数的方法，且编译好以后的UDF库拿到没有安装对应Matlab版本的机器上仍然可以使用。以前硫酸亚铜博客（https://www.cnblogs.com/liusuanyatong/p/12128082.html）提供的方法是依赖于Matlab运行库，编译好的UDF库如果拿到没有安装相应Matlab版本机器上是没法运行的。

         前面博客文章已经介绍过用VC++ UDF Studio插件实现调用误差函数erf的例子（https://www.cnblogs.com/SuperUDF/articles/16114086.html）。但是该函数输入参数是标量，和矩阵输入参数还是有所不同。所以，这里再以调用Matlab中的interp2二维插值函数为例来说明矩阵作为输入参数时的使用过程。

1. 官网下载VC++ UDF Studio插件并安装：VC++ UDF Studio，建议下载学术版（如想进一步采购注册，对高校老师学生比较优惠）

2. 安装Visual Studio（插件2022R2开始支持VS2010~2019社区，专业或旗舰版，建议安装VS2010旗舰版）， C++和C#一起安装，对于64位Fluent还要勾选X64编译器。

 

 3. 安装Matlab 2014a ~ 2021b任一版本，必须勾选Matlab Coder，其它视自己喜好安装。

 

4. 管理员权限打开桌面图标，选择需要的版本并勾选“调用Matlab”后会自动启动Fluent，读入case并点击Fluent嵌入菜单中的“Start Visual Studio”子菜单。

 

5. 把自带的matlab函数文件MatlabAdd.m改名为MatlabInterpolate.m，另外一个自动生成的文件MatlabFunctionTester.m是用来在Matlab中测试调试函数文件用的，下一步会介绍用法，这一步不用管。然后双击打开编辑MatlabInterpolate.m，输入以下自定义Matlab函数体。 

function [resultValue]= MatlabInterpolate(X,Y,V,Xi,Yi)
 resultValue=interp2(X,Y,V,Xi,Yi);
end

  

6. 双击打开MatlabFunctionTester.m，输入如下测试代码，从而便于测试我们前面定义的MatlabInterpolate函数。

clear all
[X,Y]=meshgrid(-2:0.75:2); %6*6的均布网格，每隔0.75一个点
R=sqrt(X.^2+Y.^2)+1E-6; 
V=sin(R)./(R);  %每个网格点的对应函数值
surf(X,Y,V);  %画出曲面
xlim([-2 2]); %限制显示x坐标区域为[-2,2]
ylim([-2 2]); %限制显示y坐标区域为[-2,2]
Xi=0; %要插值点的x坐标
Yi=0; %要插值点的y坐标
Vi=MatlabInterpolate(X,Y,V,Xi,Yi); %Vi就是插值得到的x=0,y=0处的值

 7. 鼠标右键在MatlabFunctionTester.m文件上单击弹出菜单，选择“用Matlab打开”，这样就可以在Matlab里面一步一步调试我们定义的MatlabInterpolate函数，排除错误后关闭Matlab。

 

 画出来的曲面和各矩阵值如下，插值运算得到Vi结果为0.9590。

 

8. 点击工具栏上“将.m文件转为C/C++”按钮，输入参数均设为float类型（双精度fluent时为double），设置矩阵尺寸请点击“矩阵”按钮，然后勾选“动态”并确定。注意：试用版没法开启“矩阵“作为输入参数，必须购买注册后才能使用。

 

 

 其中，Dyn*Dyn代表该矩阵行和列都是“动态尺寸”，例如第一个参数X转换成C/C++代码后就会表示为

const emxArray_real32_T *X

 其中，emxArray_real32_T是一个存放动态数组的结构体，其定义为

struct emxArray_real32_T
{
  float *data;  //指向存放数据的数组
  int *size;    //存放维数的数组，实际就是两个元素size[0]和size[1]
  int allocatedSize; //总的元素个数，等于行数乘以列数
  int numDimensions; //数组的维数，插件中永远是2，代表数组有行和列2个维数
  boolean_T canFreeData; //是否需要调用free来释放由calloc或malloc开辟的data数组
};

 9. 等待片刻，转换完成后，会自动将对应的转换得到的C++头文件MatlabLibrary.h加入到UDF工程中。

 转换后的C/C++函数原型为

float MatlabInterpolate(const emxArray_real32_T *X, const emxArray_real32_T *Y, const emxArray_real32_T *V, float Xi, float Yi);

10. 在udf_source.cpp文件中输入如下示例源代码，并点击“编译UDF”按钮直到编译通过。有任何错误提示，可以双击提示行直接定位到源码中的错误行。编译通过后按“UDF库加载到Fluent”按钮即可载入到Fluent中。

#include "udf.h"

extern "C"
{
#include "MatlabLibrary.h" 
}

#define MATRIX_SIZE 6

DEFINE_ON_DEMAND(Interpolate)
{
	int X_size[2]={MATRIX_SIZE, MATRIX_SIZE}; //X_size指定X_data的尺寸，本例为6*6矩阵
	int Y_size[2]={MATRIX_SIZE, MATRIX_SIZE}; //Y_size指定Y_data的尺寸，本例为6*6矩阵
	int V_size[2]={MATRIX_SIZE, MATRIX_SIZE}; //V_size指定V_data的尺寸，本例为6*6矩阵

	float X_data[]={   // X_data的具体数据，C++语言是行优先的，所以数据应该是Matlab矩阵数据按照列优先排列
		-2.0,-2.0,-2.0,-2.0,-2.0,-2.0,
		-1.25, -1.25,-1.25,-1.25,-1.25,-1.25,
		-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,
		0.25,0.25,0.25,0.25,0.25,0.25,
		1.0,1.0,1.0,1.0,1.0,1.0,
		1.75,1.75,1.75,1.75,1.75,1.75
	};

	float Y_data[]={   // Y_data的具体数据，C++语言是行优先的，所以数据应该是Matlab矩阵数据按照列优先排列
		-2.0,-1.25,-0.5,0.25,1.0,1.75,
		-2.0,-1.25,-0.5,0.25,1.0,1.75,
		-2.0,-1.25,-0.5,0.25,1.0,1.75,
		-2.0,-1.25,-0.5,0.25,1.0,1.75,
		-2.0,-1.25,-0.5,0.25,1.0,1.75,
		-2.0,-1.25,-0.5,0.25,1.0,1.75
	};

	float V_data[MATRIX_SIZE*MATRIX_SIZE];  // V_data存储坐标x,y处对应的散点值，即Matlab程序中的V

	int RowSize=X_size[0]; // 获得矩阵的行尺寸
	int ColSize=X_size[1]; // 获得矩阵的列尺寸

	for(int i=0;i<RowSize;i++)  // 填充坐标x,y对应的散点值
	{
		for(int j=0;j<ColSize;j++)
		{
			float Radius = sqrt(pow(X_data[i*RowSize+j],2)+pow(Y_data[i*RowSize+j],2))+1E-6; //即Matlab程序中的R
			V_data[i*RowSize+j]=sin(Radius)/Radius;
		}
	}
	
	emxArray_real32_T X, Y, V;
	X.data = X_data; //指向定长一维数组X_data
	X.size = X_size; //指向定长二维数组X_size
	X.allocatedSize = X_size[0] * X_size[1]; //行数乘以列数，MatlabInterpolate函数中实际并未用到，此行也可以不写
	X.numDimensions = 2; //永远为2，代表二维数组，行和列，MatlabInterpolate函数中实际并未用到，此行也可以不写
	X.canFreeData = false;  //因为X_data不是由calloc或malloc开辟的动态数组，所以不需要free，MatlabInterpolate函数中实际并未用到，此行也可以不写

	Y.data = Y_data; //指向定长一维数组Y_data
	Y.size = Y_size; //指向定长二维数组Y_size
	Y.allocatedSize = Y_size[0] * Y_size[1]; //行数乘以列数，MatlabInterpolate函数中实际并未用到，此行也可以不写
	Y.numDimensions = 2; //永远为2，代表二维数组，行和列，MatlabInterpolate函数中实际并未用到，此行也可以不写
	Y.canFreeData = false;  //因为Y_data不是由calloc或malloc开辟的动态数组，所以不需要free，MatlabInterpolate函数中实际并未用到，此行也可以不写

	V.data = V_data; //指向定长一维数组V_data
	V.size = V_size; //指向定长二维数组V_size
	V.allocatedSize = V_size[0] * V_size[1]; //行数乘以列数，MatlabInterpolate函数中实际并未用到，此行也可以不写
	V.numDimensions = 2; //永远为2，代表二维数组，行和列，MatlabInterpolate函数中实际并未用到，此行也可以不写
	V.canFreeData = false;  //因为V_data不是由calloc或malloc开辟的动态数组，所以不需要free，MatlabInterpolate函数中实际并未用到，此行也可以不写

	float Xi=0, Yi=0; //插值的目标x,y值
	float Vi=MatlabInterpolate(&X,&Y,&V, Xi,Yi); //调用Matlab插值函数
	Message0("The interpolated value at coordinate (0,0) is %g\n", Vi);
}

如果出现INFINITY，NAN未声明的标识符的错误，那么请使用较高版本的Visual Studio，例如Visual Studio2015或更高。

11. 执行DEFINE宏，本例由于插值函数放在DEFINE_ON_DEMAND宏中，所以在Execute On Demand对话框里面手动执行。

 

 12. Fluent中运行结果如下，和前面直接Matlab里面的结果是一致的。