matlab数据转换为tecplot格式[ASCII码格式下的plt文件]

hyhhyh21

已于 2024-08-27 21:16:22 修改

阅读量1.6w

点赞数 57

文章标签： matlab tecplot txt plt ASCII码

于 2020-02-08 10:31:56 首次发布

本文链接：https://blog.csdn.net/weixin_42943114/article/details/104204172

版权

本文详细介绍了如何将MATLAB数据转换为tecplot支持的ASCII格式plt文件，包括文件格式、文件头、zone相关格式，以及一维、二维、三维数据的转换方法。 tecplot能读取ASCII格式的plt文件，支持多种数据类型和网格结构，如非矩形区域、多块数据和非结构网格。通过示例代码和tecplot的使用说明，帮助读者理解数据转换过程。

摘要由CSDN通过智能技术生成

2021.06更新。（1）修改了之前1.1节演示代码，之前ZONE前面忘记加注释（2）新增了结构网格的多块数据的写入说明。感谢评论区的老铁们。
2021.08更新。（1）添加了离散点数据以及非结构数据的导入。
2022.08更新。（1）偷不了懒，ZONETYPE还差两个多边形导入的没有提到，果然评论区就有人指出来了，O(∩_∩)O哈哈~，更新这两个FEPolygon和FEPolyhedron。

1 tecplot能读取的的常见ASCII格式

1.1 文件格式与文件头

tecplot默认的数据读取与保存格式为二进制格式。与ASCII码格式相比，tecplot对二进制格式的读取数据速度更快。但是为了更方便的文件交互与人为校验，对于数据量不太大的ASCII格式也是必须的。

本文输出的文件后缀为*.plt，由于是ASCII码格式，所以可以用记事本等软件打开。文件命名或者路径可以包含中文。对于单个文件导入时，可以采用load data，或者直接拖拽的方式读取文件。多个文件（比如非定常）不支持拖拽，需要load data。

在文件中，以#号开头的行会作为注释而忽略。字符串需要采用双引号"X"进行括起来，如果字符串内包含双引号，则可以采用反斜杠的方式输入（"）。

文件对换行不是很敏感，多余的换行会被忽略掉（除了注释）。双引号内的字符串不能换行。文件也对大写小写不太敏感，ZONE和zone都可以被相同的识别，除双引号内的字符串。如果遇到并列，采用空格还是逗号也同样无所谓，亦或是混合搭配，这两种方式都可以被识别。

文件中的数字可以采用整数 (101325)，小数(101325.0)，或科学计数法 (1.01325E+05)。数字和数字之间可以用空格、逗号、Tab、回车、换行进行分隔。如果有连续的重复数字，例如：37个120.5，后面跟着100个0.0，可以简写作:37120.5,1000.0。

一个tecplot文件，可以分为文件头和数据两部分。文件头需要标明：TITLE名称、FILETYPE文件类型和VARIABLES变量名称。

TITLE，后面内容为字符串，格式为 TITLE = “Example File”
FILETYPE，可以选择为GRID（网格）、SOLUTION（结果）、FULL（全部）。这一项可以省略，默认为FULL。
VARIABLES，后面接数据名称。同样为字符串格式，如果有多个数据，用空格或逗号隔开。

数据则分为 ZONE、TEXT、GEOMETRY、CUSTOMLABELS、 DATASETAUXDATA 和 VARAUXDATA。本文以zone区域进行演示，这也是最主要被用到的。

其中zone分为 Block和Point两种数据格式，这在之后的章节再解释。

一个典型的tecplot下ASCII格式的plt文件，如下所示：

#保存为Example.plt
#标题
TITLE = "Example File"
#变量名
VARIABLES = "X","F"
#定义ZONE
#ZONE
#定义IJK，一维只需要写I
zone I = 11,t="F=x*x"
DATAPACKING=point
0 0
0 .1 0.01
0.2 0.04
0.3 0.09
0.4 0.16
0.5 0.25
0.6 0.36
0.7 0.49
0.8 0.64
0.9 0.81
1 1

zone I = 11,t="F=x"
DATAPACKING=block
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

导入到tecplot，勾选symbol，和mapping style 中的show map选项，显示的图像如下：
在这里插入图片描述
更多相关信息，可以搜索tecplot的帮助文档，我的版本对应的pdf文档名称为：360_data_format_guide.pdf。

1.2 zone相关的格式

数据以zone开头，标志着后面所有内容均为其下的内容，直到另一个数据类型标识为止。

必须指定的变量有：

变量	含义
I,J,K	数据点的个数。如果是二维，只需要定义IJ。如果是1维，只需定义I

可以不指定的变量有：

变量	含义
T	zone的名称
ZONETYPE	可选ordered（默认）和 finite element两种。有限元需要额外的nodes、elements、faces 等信息
DATAPACKING	可选block（默认）和point两种。point可认为是每一行代表一个点，block可认为是每一列代表一个点（虽然tecplot不在乎换行）
DT	储存精度，默认Single。可选Double,single,longint,shortint,byte,bit。其中bit只支持Block格式。
VARLOCATION	可选 NODAL（默认）和CELLCENTERED。其中Cellcentered只支持Block格式，允许在中心定义变量，此时数据点个数应为(I-1)×(J-1)×(K-1)
SOLUTIONTIME	求解时间

2 一维数据

一维数据，以抛物线为示例。主要代码为plt_Head和plt_Zone。

clear
fclose('all');


t=0:0.1:1;
F=t.^2;

%1准备数据
filename='example_1D.plt';
title=[];
variables={'X','F'};
Mat_Data=[t(:),F(:)];
IJK=length(t);

%2创建文件
if exist(filename,'file') 
    delete(filename)
end
f_id=fopen(filename,'a');
fclose(f_id);

%3创建表头
plt_Head(filename,'',variables)
%4创建zone（point）格式
plt_Zone(filename,'1',IJK,Mat_Data)



function plt_Head(filename,title,variables)
%创建表头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end
%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone(filename,zone_title,IJK,Mat_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);

Dim=numel(IJK);
if Dim==1
    s=['zone I=',num2str( IJK(1) )];
elseif Dim==2
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) )];
elseif Dim==3
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) ),',K=',num2str( IJK(3) )];
end

%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%格式为point
s='DATAPACKING=point';
fprintf(f_id,'%s \r\n',s);

%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
fclose(f_id);
end

导入到tecplot的结果为
在这里插入图片描述

3 二维数据

二维数据，采用标准的xy网格。包括X、Y、U、V、P5个变量。主要代码为plt_Head和plt_Zone。
代码和一维相似，如下所示。

这里要注意IJ的排列顺序，Tecplot里的IJ等于matlab里size(meshgrid(**))生成的尺寸。Tecplot里数据的读取顺序，等效于meshgrid生成的数据，再进行x(:)操作的读取顺序。tecplot里如果这里出错会造成读取错误，但两种软件相似的数据读取方式也减少了不必要的麻烦。

clear
fclose('all');

%0数据
[x,y]=meshgrid(-5:0.5:5,-4:0.5:4);
u=-sin(x);
v=cos(y);
F=u.^2+v.^2;

%1准备数据
filename='example_2D.plt';
title='';%无标题
variables={'X','Y','U','V','F'};
zone_title='';%无标题
Mat_Data=[x(:),y(:),u(:),v(:),F(:)];
IJK=size(x);
time=[];%无时间
 
%2创建文件
if exist(filename,'file') 
    delete(filename)
end
f_id=fopen(filename,'a');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone(filename,zone_title,IJK,time,Mat_Data)




function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone(filename,zone_title,IJK,time,Mat_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);

Dim=numel(IJK);
if Dim==1
    s=['zone I=',num2str( IJK(1) )];
elseif Dim==2
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) )];
elseif Dim==3
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) ),',K=',num2str( IJK(3) )];
end


%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%格式为point
s='DATAPACKING=point';
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
fclose(f_id);
end

将生成的plt文件，导入tecplot生成结果为：
在这里插入图片描述
对于二维非定常问题，需要添加SOLUTIONTIME选项，生成程序如下。自定义函数plt_Head和plt_Zone与上一个程序内的相同。

clear
fclose('all');

k=0;
for t=0.2:0.2:3
    k=k+1;
    %0数据
    [x,y]=meshgrid(-5:0.5:5,-4:0.5:4);
    u=-sin(x+t);
    v=cos(y);
    F=u.^2+v.^2;
    
    %1准备数据
    filename=['example_2D',num2str(k,'%02.f'),'.plt'];
    title='';%无标题
    variables={'X','Y','U','V','F'};
    zone_title='';%无标题
    Mat_Data=[x(:),y(:),u(:),v(:),F(:)];
    IJK=size(x);
    time=t;%非定常时间
    
    %2创建文件
    if exist(filename,'file')
        delete(filename)
    end
    f_id=fopen(filename,'a');
    fclose(f_id);
    
    %3创建文件头
    plt_Head(filename,title,variables)
    %4创建zone（point）格式
    plt_Zone(filename,zone_title,IJK,time,Mat_Data)

end



function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end

function plt_Zone(filename,zone_title,IJK,time,Mat_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);

Dim=numel(IJK);
if Dim==1
    s=['zone I=',num2str( IJK(1) )];
elseif Dim==2
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) )];
elseif Dim==3
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) ),',K=',num2str( IJK(3) )];
end

%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%格式为point
s='DATAPACKING=point';
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end

%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
fclose(f_id);
end

生成的plt文件有15个，每个文件代表某一时刻的流场。打开Tecplot，选择File-Load Data，就可以导入所有时刻的流场。选择Animation-Time就可以进行非定常展示。可以添加文本框，输入 t=&(solutiontime)s即可以显示文字版时间。
在这里插入图片描述

4 三维数据

三维数据采用matlab自带的wind数据。输入load wind，便可得到包含meshgrid格式的x,y,z,u,v,w流场数据。

程序如下所示，其中自定义函数plt_Head和plt_Zone与第3节程序内的相同。

clear
fclose('all');
load wind

V2=sqrt(u.^2+v.^2+w.^2);

%1准备数据
filename=['example_3D.plt'];
title='';%无标题
variables={'X','Y','Z','U','V','W','V2'};
zone_title='';%无标题
Mat_Data=[x(:),y(:),z(:),u(:),v(:),w(:),V2(:)];
IJK=size(x);
time=[];

%2创建文件
if exist(filename,'file')
    delete(filename)
end
f_id=fopen(filename,'a');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone(filename,zone_title,IJK,time,Mat_Data)

function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone(filename,zone_title,IJK,time,Mat_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);

Dim=numel(IJK);
if Dim==1
    s=['zone I=',num2str( IJK(1) )];
elseif Dim==2
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) )];
elseif Dim==3
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) ),',K=',num2str( IJK(3) )];
end


%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%格式为point
s='DATAPACKING=point';
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
fclose(f_id);
end

生成example_3D.plt文件，导入tecplot内，可以生成如下的流场示意图。
在这里插入图片描述

4.1 非矩形区域的数据导入

非矩形数据也可以用ijk的形式导入，只要是结构网格。一般的思路是：IJK结构矩形网格→uvw中间参数网格→xyz最终网格。当然网格简单的话，可以直接跳过ijk网格的构造，或跳过uvw网格的构造。复杂的网格建议利用其它结构网格划分软件进行划分。

下面以一个圆柱坐标系为例。

clear
fclose('all');
%非矩形的结构网格构造方法
%以柱面坐标为例
r0=1:0.2:10;
z0=0:0.1:2;
theta0=0:0.05:pi;

[r,theta,z]=meshgrid(r0,theta0,z0);
%转换为xyz坐标系
x=r.*cos(theta);
y=r.*sin(theta);

p=cos(r).*sin(4*theta);

%1准备数据
filename=['example_3D.plt'];
title='';%无标题
variables={'X','Y','Z','p'};
zone_title='';%无标题
Mat_Data=[x(:),y(:),z(:),p(:)];
IJK=size(x);
time=[];

%2创建文件
if exist(filename,'file')
    delete(filename)
end
f_id=fopen(filename,'a');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone(filename,zone_title,IJK,time,Mat_Data)

在这里插入图片描述
当然如果遇到复杂的结构，单一的结构网格就不太够用了。需要引入多块结构网格，或者非结构网格。

5 多块数据

当block存在多个的时候，比如结构网格是由多个block组成的时候，也可以用plt文件进行导入。

这里只需要保持总的标题不变，其余各个block的按照之前的设置依次输入就行。

同样代码如下。
其中自定义函数plt_Head和plt_Zone与第3节和第4节程序内的相同。

%生成多块网格
clear
clc
fclose('all');
%0生成数据
%block1
load wind
x_1=x;y_1=y;z_1=z;
u_1=u;v_1=v;w_1=w;
V2_1=sqrt(u_1.^2+v_1.^2+w_1.^2);
%block2
x_2=x/1.2;y_2=y;z_2=z-min(z_1,[],'all')+max(z_1,[],'all');%这里稍微更改一下网格
u_2=u;v_2=v;w_2=w;
V2_2=sqrt(u_2.^2+v_2.^2+w_2.^2);

%1准备数据
filename=['example_3D_2Blocks.plt'];
title='';%多block只对应1个文件标题。这里为了省事，设置无标题。
variables={'X','Y','Z','U','V','W','V2'};
zone_title='';%无标题
Mat_Data_1=[x_1(:),y_1(:),z_1(:),u_1(:),v_1(:),w_1(:),V2_1(:)];
Mat_Data_2=[x_2(:),y_2(:),z_2(:),u_2(:),v_2(:),w_2(:),V2_2(:)];
IJK_1=size(x_1);
IJK_2=size(x_2);
time=[];

%2创建文件
if exist(filename,'file')
    delete(filename)
end
f_id=fopen(filename,'a');
fclose(f_id);

%Block1
%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone(filename,zone_title,IJK_1,time,Mat_Data_1)
%Block2
%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone(filename,zone_title,IJK_2,time,Mat_Data_2)

function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone(filename,zone_title,IJK,time,Mat_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);

Dim=numel(IJK);
if Dim==1
    s=['zone I=',num2str( IJK(1) )];
elseif Dim==2
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) )];
elseif Dim==3
    s=['zone I=',num2str( IJK(1) ),',J=',num2str( IJK(2) ),',K=',num2str( IJK(3) )];
end


%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%格式为point
s='DATAPACKING=point';
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
fclose(f_id);
end

在这里插入图片描述

6 非结构网格

6.1 二维平面网格

tecplot的非结构数据形式也比较简单。分为两个部分，第一个部分为每个节点处的信息，第二部分为非结构网格的每个网格的由哪几个点构成。

下面为一个简单的二维非结构的plt

VARIABLES = "X", "Y", "P"
ZONE T="TRIANGLES", NODES=4, ELEMENTS=2, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
#每个点的信息（ NODES=4）
0.0 0.0 0.0
0.0 1.0 1.2
1.0 1.0 -2.0
1.0 0.0 0.0
#每个网格由哪几个点组成（ELEMENTS=2）
1 2 3
1 3 4

在这里插入图片描述

非结构网格的ZONETYPE分为

ZONETYPE	网格	点数
FETRIANGLE	二维三角网格	3点
FEQUADRILATERAL	二维四边形网格	4点
FETETRAHEDRON	三维四面体网格	4点
FEBRICK	三维六面体网格	8点

二维或者曲面的，使用FETRIANGLE和FEQUADRILATERAL就可以。

对于混合了三角形和四边形的网格，按照四边形网格计算就行。但是三角形区域需要通过重复，变成ABCC的形式，凑齐4个点的格式。
还是举个例子：

VARIABLES = "X", "Y", "P"
ZONE T="TRIANGLES", NODES=6, ELEMENTS=3, DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL
#每个点的信息
0.0 0.0 0.0
0.0 1.0 0.0
1.0 1.0 -1.0
1.0 0.0 -1.0
-1.0 0.0 1.0
-1.0 1.0 1.0
#每个网格由哪几个点组成
1 2 3 4
5 6 2 2
1 2 5 5

在这里插入图片描述
左侧两个三角形网格，需要通过5 6 2 2的形式，重复一个节点，来凑齐ZONETYPE = FEQUADRILATERAL 所要求的每个网格4个点的条件。

复杂的非结构网格需要自己依据的模型来构造，也可以借用其他网格划分软件进行网格定义。这里不是本文的重点。

下面给出一个已知很多离散点的信息，在tecplot里绘图的例子：

clear
fclose('all');
%非结构面网格的处理(二维)
%这里是生成的数据
x=gallery('uniformdata',[300,1],0);
y=gallery('uniformdata',[300,1],1);
p1=sin((x+y)/2*4*pi);
p2=cos(sqrt(x.*y)*4*pi);
%1准备数据
filename=['example_2D_Surface_U.plt'];
title='';%无标题
variables={'X','Y','P1','P2'};
zone_title='';%无标题
%准备每个点的信息
Mat_Data=[x(:),y(:),p1(:),p2(:)];
%准备每个非结构网格的信息
Tri_Data=delaunay(x,y);%这里默认是实心的不规则离散点。其它特殊形状的不适用。
time=[];

%2创建文件
f_id=fopen(filename,'w');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone_U(filename,zone_title,time,Mat_Data,Tri_Data)


function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone_U(filename,zone_title,time,Mat_Data,Tri_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);
E=size(Tri_Data,1);

s=['ZONE NODES=',num2str(N),', ELEMENTS=',num2str(E),', DATAPACKING=POINT, ZONETYPE=FETRIANGLE'];

%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
%导入网格
for k=1:E
    fprintf(f_id,'%s \r\n',num2str(Tri_Data(k,:)));
end
fclose(f_id);
end

在这里插入图片描述

6.2 简单的曲面演示

曲面的构造形式和二维一样，这里也简单做一个小演示。
为了演示，这里用的是四边形的非结构网格，ZONETYPE=FEQUADRILATERAL。

clear
fclose('all');
%非结构面网格的处理(二维)
%这里是生成的数据
[x,y]=meshgrid(-3:0.2:3,-3:0.2:3);
z = 0.5*peaks(x,y);
p1=sin((x+y)/2*pi);
p2=cos(sqrt(x.*y)*4*pi);
%1准备数据
filename=['example_2D_Surface_U4.plt'];
title='';%无标题
variables={'X','Y','Z','P1','P2'};
zone_title='';%无标题
%准备每个点的信息
Mat_Data=[x(:),y(:),z(:),p1(:),p2(:)];
%准备每个非结构网格的信息(这里用的是投影在xy平面上的点，划分二维网格)
Tri_Data=delaunay(x,y);%这里默认是实心的不规则离散点。其它特殊形状的不适用。
Nx=size(x,1);
Ny=size(y,2);
Quad_Data=zeros((Nx-1)*(Ny-1),4);
n=1;
for k=1:Nx-1
    for j=1:Ny-1
        Quad_Data(n,1)=sub2ind(size(x),k,j);
        Quad_Data(n,2)=sub2ind(size(x),k+1,j);
        Quad_Data(n,3)=sub2ind(size(x),k+1,j+1);
        Quad_Data(n,4)=sub2ind(size(x),k,j+1);
        n=n+1;
    end
end

time=[];
%2创建文件
f_id=fopen(filename,'w');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone_U4(filename,zone_title,time,Mat_Data,Quad_Data)


function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone_U4(filename,zone_title,time,Mat_Data,Tri_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);
E=size(Tri_Data,1);

s=['ZONE NODES=',num2str(N),', ELEMENTS=',num2str(E),', DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL'];

%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
%导入网格
%四点网格
for k=1:E
    fprintf(f_id,'%s \r\n',num2str(Tri_Data(k,:)));
end
fclose(f_id);
end

效果如下（后面的背景是利用tecplot上面菜单栏的Plot-Axis，点击选中Show-X-Axis、Show-Y-Axis、Show-Z-Axis画的）。
在这里插入图片描述

6.3 三维非结构网格

之前介绍了二维非结构网格。三维的网格和二维网格思路上一样，就是zonetype不一样。

三维这里使用FETETRAHEDRON 和 FEBRICK ，分别对应三维四面体网格和三维六面体网格。

如果遇到了混合网格，或不规则的三棱柱、四棱锥等网格，处理方法也和二维的思路一样，统一采用FEBRICK格式，然后将那些没有凑齐8个点的网格，都利用重复补齐至8个点。

下面同样给出一个已知很多离散点的信息，在tecplot里绘图的例子：

clear
fclose('all');
%非结构体网格的处理(三维)
%这里是生成的数据
x=gallery('uniformdata',[800,1],0);
y=gallery('uniformdata',[800,1],1);
z=gallery('uniformdata',[800,1],2);
p1=sin((x+y+z)/3*4*pi);
p2=cos((x.*y.*z).^(1/3)*4*pi);
%1准备数据
filename=['example_3D_Block_U.plt'];
title='';%无标题
variables={'X','Y','Z','P1','P2'};
zone_title='';%无标题
%准备每个点的信息
Mat_Data=[x(:),y(:),z(:),p1(:),p2(:)];
%准备每个非结构网格的信息
Tri_Data=delaunay(x,y,z);%这里默认是实心的不规则离散点。其它特殊形状的不适用。
time=[];

%2创建文件
f_id=fopen(filename,'w');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone_U3(filename,zone_title,time,Mat_Data,Tri_Data)


function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end

%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);

end


function plt_Zone_U3(filename,zone_title,time,Mat_Data,Tri_Data)
%创建zone,point格式
f_id=fopen(filename,'a');
N=size(Mat_Data,1);
E=size(Tri_Data,1);

s=['ZONE NODES=',num2str(N),', ELEMENTS=',num2str(E),', DATAPACKING=POINT, ZONETYPE=FETETRAHEDRON'];

%标题
if ~isempty(zone_title)
    s=[s,',t="',zone_title,'"'];
end
fprintf(f_id,'%s \r\n',s);
%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end


%导入数据
for k=1:N
    fprintf(f_id,'%s \r\n',num2str(Mat_Data(k,:)));
end
%导入网格
for k=1:E
    fprintf(f_id,'%s \r\n',num2str(Tri_Data(k,:)));
end
fclose(f_id);
end

在这里插入图片描述

7 多面体网格

多面体网格顾名思义，就是不限制每个网格的形状，对于二维可以是三角形四边形，也可以是更多边形状；对于三维，则不一定是四面体，六面体，也可以是其它异型图形。

Tecplot官方也给出了多面体网格的构建例子，如examples文件夹中的Sphere.lpk和Pyramid.plt文件。

多面体网格对应的ZONETYPE为下面两种

ZONETYPE	网格	点数
FEPolygon	二维多边形网格	不限
FEPolyhedron	三维多面体网格	不限

7.1 二维多边形网格

具体格式为：

VARIABLES=....这里填变量名称
ZONE 
NODES= 节点数，和后面block里的数量要对应
FACES= 边的数量，和后面face里的数量要对应
ELEMENTS= 网格（或者说多边形）的数量，和后面left和right的最大值要对应
DATAPACKING =BLOCK 这里必须是block不能是point
ZONETYPE=FEPOLYGON 这是二维多边形网格的zonetype

TotalNumFaceNodes= 等于后面face里的元素总数量。因为对于2维，每个边只有两个点，所以等于2*FACES
NumConnectedBoundaryFaces=0 涉及多个zone交互，本文暂不考虑
TotalNumBoundaryConnections=0 涉及多个zone交互，本文暂不考虑

后面是具体的数据
1 Block数据
2 Faces 多边形的所有边，格式为N行2列的形式，代表每个边是由哪两个点组成的。
3 left elements 一行，每条边左手的那个多边形编号，如果没有就是0
4 right elements 一行，每条边右手的那个多边形编号，如果没有就是0

下面举一个例子，只有左右两个正方形网格。这样共有6个点，点的编号见下图。共有7条边，每个边由两个点组成。最后是每条边的左边和右边网格编号。

VARIABLES =  "X", "Y"
ZONE
NODES=6
FACES=7
ELEMENTS=2
DATAPACKING=BLOCK
ZONETYPE=FEPOLYGON 
#TotalNumFaceNodes是FACES那一区域的数量，对于2D，就是2*FACES
TotalNumFaceNodes=14
NumConnectedBoundaryFaces=0
TotalNumBoundaryConnections=0
#NODES（坐标）
0.0 1.0 2.0 0.0 1.0 2.0
0.0 0.0 0.0 1.0 1.0 1.0
# node count per face（下面face读取的时候会用到）（2D不需要，默认都是2）
#2 2 2 2 2 2 2 （相加需要等于14）
#FACES 7行2列，共14个
1 2
2 3
3 6
6 5
5 4
4 1
2 5
# left elements（在这个face左侧的单元）
1 2 2 2 1 1 1
# right elements
0 0 0 0 0 0 2

请添加图片描述
小注释：
1 由于plt文件对于换行不敏感，所以实际有些文件会非常混乱，不像上面例子分着这几部分这么清晰
2 plt文件每一行字数似乎有限制，不建议太长，建议10个数字就换一行
3 最后的左和右如果判断不清也没关系，tecplot不太注意这个。比如“左1右2”写成了“左2右1”，tecplot也能正常导入不会报错。可能会影响光照渲染，但不影响网格结构。

知道了这个简单的例子后，就可以扩展，导入和整理较为复杂的图形。

下面用一个近六边形的网格，同6.2节，在tecplot绘制matlab图形，代码见下。

这里面核心代码只有第二部分（2数据输出部分），前面0和1都是在生成数据和整理数据。由于不同数据格式不一样，需要整理的方法也不一样，所以本文只是做一个示范。具体自己的数据，如果格式不同，还需要按照自己的方法整理。

clear
clc
close all
%最终决定用这个输出
%% 0数据准备，这个程序里的数据是自己生成创建的
%绘制六边形格式的散点
x=0;y=0;

N=6;
th1=linspace(0,2*pi,N+1);

for k=1:20
    R=(k^1.2)*sqrt(3);
    x2=R*sin(th1);
    y2=R*cos(th1);
    %然后第几层六边形，就扩充几个
    N1=0:N*k;
    N2=0:k:N*k;
    x3=interp1(N2,x2,N1);
    y3=interp1(N2,y2,N1);
    %删除最后一个重复点
    x3=x3(2:end);y3=y3(2:end);
    %合并，删除重复的点
    x=[x,x3];
    y=[y,y3];
    [~,ia,~]=unique([x;y]','stable','rows');
    x=x(ia);
    y=y(ia);
end

x=x+0.5*gallery('normaldata',size(x),1);
y=y+0.5*gallery('normaldata',size(x),2);


MaxXY=k^1.2*sqrt(3);
XYShow=8;%缩放到[-8,8]这个区间
dN=MaxXY/XYShow;
x=x/dN;
y=y/dN;


%绘制voronoin网格
[v,c]=voronoin([x;y]');
%删除不用的网格
x2=XYShow*sin(th1);
y2=XYShow*cos(th1);%划分界限
IfInBoundary=@(x) all(inpolygon( v(x,1),v(x,2),x2,y2 ));
c=c(cellfun(IfInBoundary,c));%只用在界限范围内的网格
%有些点没有用上，所以删掉
Node_Sum=[];
for k=1:size(c,1)
    Node_Sum=[Node_Sum,c{k}];
end
Node_unique=unique(Node_Sum);%最终只输出实际用到的点
indx=true(size(v,1),1);
indx(Node_unique)=false;
v(indx,:)=0;%用不到的点赋值为0
%% 1数据整理，这里只是演示一下，怎么样转换数据格式。
%不同数据的数据格式转换方法不一样，所以可能需要重新编程整理
%把[v,c]数据格式转换成FACES_Data和FElement_Data格式
%FACES_Data是所有边的数据，N行2列，N是网格边的总数量
%FElement_Data是单元网格的数据，分为Left和Right两个。分别为每个边左边网格和右边网格信息。

%构建输出的Mat_Data
vx=v(:,1);
vy=v(:,2);
z=peaks(vx,vy)/3;
p1=sin(vx)+cos(vy);
Mat_Data=[vx(:),vy(:),z(:),p1(:)];
%先把c的数据按照逆时针顺序排序
for k=1:size(c,1)
    c_k=c{k};
    x_k=v(c_k,1);
    y_k=v(c_k,2);
    if trapz([x_k;x_k(1)],[y_k;y_k(1)])>0
        c{k}=fliplr(c_k);
    end
end
%这个Line数据好整理，就是读取所有边
FACES_Data=[];
for k=1:size(c,1)
    c_k=c{k};
    Line_k=zeros(length(c_k),2);
    %按顺序读取边
    c_k=[c_k,c_k(1)];
    Line_k(:,1)=c_k(1:end-1);
    Line_k(:,2)=c_k(2:end);
    %合并
    FACES_Data=[FACES_Data;Line_k];
end
%删除重复边
FACES_Data=sort(FACES_Data,2);
FACES_Data=unique(FACES_Data,'rows','stable');
%然后是判断左右。因为前面每个element已经按照逆时针整理过了，所以边序号相同就是左，相反就是右。
NF=size(FACES_Data,1);
Left_Elements=zeros(1,NF);
Right_Elements=zeros(1,NF);
for k=1:NF
    FACES_k=FACES_Data(k,:);%当前边
    FACES_kF=FACES_Data(k,:);%当前边反
    %对于2维网格，每个边只对应2个element，一个左，一个右
    for kc=1:size(c,1) %对每个element都找一下（这个方法慢，但直观，先不优化）
        c_k=c{kc};
        IfInList1=Find_Line(FACES_k,c_k);
        if IfInList1==1
            Left_Elements(k)=kc;
        elseif IfInList1==2
            Right_Elements(k)=kc;
        end
    end
end
FElement_Data.N_Element=size(c,1);
FElement_Data.Left=Left_Elements;
FElement_Data.Right=Right_Elements;

%% 2数据输出
%1准备数据
filename=['example_2D_Surface_UN.plt'];
title='';%无标题
variables={'X','Y','Z','P1'};%和前面Mat_Data对应
zone_title='';%无标题

time=[];
%2创建文件
f_id=fopen(filename,'w');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone_UN(filename,zone_title,time,Mat_Data,FACES_Data,FElement_Data)

%% 展示网格
figure(1)
hold on
for i=1:length(c)
    x_plot=v(c{i},1);
    y_plot=v(c{i},2);
    patch(x_plot,y_plot,i);
end
hold off
axis equal


%% 后置函数
function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end
%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);
end


function plt_Zone_UN(filename,zone_title,time,Nodes_Data,FACES_Data,FElement_Data)
%创建zone,block格式（ZONETYPE=FEPOLYGON时必须采用这个形式）
f_id=fopen(filename,'a');
N=size(Nodes_Data,1);
F=size(FACES_Data,1);
E=FElement_Data.N_Element;
Left_Elements=FElement_Data.Left;
Right_Elements=FElement_Data.Right;
s='ZONE';
fprintf(f_id,'%s \r\n',s);
s=['NODES=',num2str(N),', FACES=',num2str(F),', ELEMENTS=',num2str(E),', DATAPACKING=BLOCK, ZONETYPE=FEPOLYGON'];
fprintf(f_id,'%s \r\n',s);
s=['TotalNumFaceNodes=',num2str(2*F)];
fprintf(f_id,'%s \r\n',s);
s=['NumConnectedBoundaryFaces=0'];
fprintf(f_id,'%s \r\n',s);
s=['TotalNumBoundaryConnections=0'];
fprintf(f_id,'%s \r\n',s);
%标题
if ~isempty(zone_title)
    s=['t="',zone_title,'"'];
    fprintf(f_id,'%s \r\n',s); 
end

%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end

%变成BLOCk格式
Nodes_Data=Nodes_Data';
%导入数据 
for k=1:size(Nodes_Data,1)
    fprintf(f_id,'%s \r\n',['# 变量',num2str(k)]);
    Mat_Data_k=Nodes_Data(k,:);
    if N>10
        %为了防止太长导致读取失败，把数据再拆成每10个一行
        N10=fix(numel(Mat_Data_k)/10);
        Mat_2=zeros(10,N10);
        Mat_2(:)=Mat_Data_k(1:N10*10);
        Mat_2=Mat_2';
        for k2=1:N10
            fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
        end
        fprintf(f_id,'%s \r\n',num2str(Mat_Data_k(N10*10+1:end)));
    else
        fprintf(f_id,'%s \r\n',num2str(Mat_Data_k));
    end
end

%导入Face边
fprintf(f_id,'%s \r\n','# Faces 多边形的所有边');
for k=1:size(FACES_Data,1)
    fprintf(f_id,'%s \r\n',num2str(FACES_Data(k,:)));
end

%导入left elements（在这个face左侧的单元）
fprintf(f_id,'%s \r\n','# left elements（在这个face左侧的单元）');
%为了防止太长导致读取失败，把数据再拆成每10个一行
if numel(Left_Elements)>10
    N10=fix(numel(Left_Elements)/10);
    Mat_2=zeros(10,N10);
    Mat_2(:)=Left_Elements(1:N10*10);
    Mat_2=Mat_2';
    for k2=1:N10
        fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
    end
    fprintf(f_id,'%s \r\n',num2str(Left_Elements(N10*10+1:end)));
else
    fprintf(f_id,'%s \r\n',num2str(Left_Elements));
end

%导入right elements（在这个face右侧的单元）
fprintf(f_id,'%s \r\n','# right elements（在这个face右侧的单元）');
%为了防止太长导致读取失败，把数据再拆成每10个一行
if numel(Right_Elements)>10
    N10=fix(numel(Right_Elements)/10);
    Mat_2=zeros(10,N10);
    Mat_2(:)=Right_Elements(1:N10*10);
    Mat_2=Mat_2';
    for k2=1:N10
        fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
    end
    fprintf(f_id,'%s \r\n',num2str(Right_Elements(N10*10+1:end)));
else
    fprintf(f_id,'%s \r\n',num2str(Right_Elements));
end
fclose(f_id);
end

function IfInList=Find_Line(AB,List)
%查找AB直线是否在List中，如果在，输出在左或在右
%0不在，1在左，2在右
List=[0,List,List(1),0];
indx1=find(List==AB(1));
IfInList=0;%默认是0
if ~isempty(indx1)
    %左侧（顺序相同）
    for k=1:numel(indx1)
        if List(indx1(k)+1)==AB(2)
            IfInList=1;%找到相同AB，输出
            break
        end
    end
    %右侧（顺序相反）
    for k=1:numel(indx1)
        if List(indx1(k)-1)==AB(2)
            IfInList=2;%找到BA，输出
            break
        end
    end
end
end

最终在Tecplot里的效果如下：
请添加图片描述

7.2 三维多面网格

下面举一个三维多面体网格

VARIABLES =  "X", "Y", "Z",
ZONE
NODES=9
FACES=10
ELEMENTS=2
DATAPACKING=BLOCK
ZONETYPE=FEPolyhedron
#node count的数量为36（4*6+3*4）
TotalNumFaceNodes=36 
NumConnectedBoundaryFaces=0
TotalNumBoundaryConnections=0
#NODES
0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.5
0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0 0.5
0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 2.0
# node count per face（下面face读取的时候会用到，每个面点的数量）
4 4 4 4 4 4 3 3 3 3
#FACES 10行，每个面一行
1 2 4 3
6 8 4 2
1 3 7 5
3 4 8 7
5 6 2 1
7 8 6 5
9 8 6
9 6 5
9 5 7
9 7 8
# left elements（在这个face左侧的单元）
0 0 0 0 0 2 0 0 0 0
# right elements 按照上面那个面，给出的点的顺序，按照右手定则，正方向的在右，反方向则在左。
1 1 1 1 1 1 2 2 2 2

请添加图片描述
同样，这个如果有了上面二维的基础下，理解3维的数据结构也很容易。就是多了一个node count per face这一块，因为二维一条边只有两个点组成，但三维一个面却不一定有多少点组成，所以需要单独列出这一块。

之后以一个三维多边形来展示：

clear
clc
close all
%最终决定用这个输出
%3D输出
%% 0数据准备，这个程序里的数据是自己生成创建的
[x,y,z]=ndgrid(0:0.25:1,0:0.25:1,0:0.25:1);
P0=[x(:),y(:),z(:)];
P0=P0+0.1*(gallery('uniformdata',size(P0),1)-0.5);
%gallery('uniformdata',[100,3],1)

[v,c]=voronoin(P0);
%删除不用的网格
xmin=-0.1;
xmax=1.1;
IfInBox=@(x) all(v(x,1)<xmax) && all(v(x,2)<xmax) && all(v(x,3)<xmax) &&...
    all(v(x,1)>xmin) && all(v(x,2)>xmin) && all(v(x,3)>xmin);
c=c(cellfun(IfInBox,c));%只用在[0,1]界限范围内的网格
%有些点没有用上，所以删掉
Node_Sum=[];
for k=1:size(c,1)
    Node_Sum=[Node_Sum,c{k}];
end
Node_unique=unique(Node_Sum);%最终只输出实际用到的点
indx=true(size(v,1),1);
indx(Node_unique)=false;
v(indx,:)=0;%用不到的点赋值为0
%% 1数据整理，这里只是演示一下，怎么样转换数据格式。
%不同数据的数据格式转换方法不一样，所以可能需要重新编程整理
%把[v,c]数据格式转换成FACES_Data和FElement_Data格式
%FACES_Data是所有面的数据，N行，N是网格面的总数量，每行由各个面的点按照一定顺序输出
%FElement_Data是单元网格的数据，分为Left和Right两个。分别为每个边左边网格和右边网格信息。
%FElement_Data的左右似乎不是很严格
%构建输出的Mat_Data
vx=v(:,1);
vy=v(:,2);
vz=v(:,3);
p1=sin(2*pi*vx)+cos(2*pi*vy)+cos(2*pi*vz);
Mat_Data=[vx(:),vy(:),vz(:),p1(:)];
%读取所有的面（这里用convhulln函数全部转化成三角形面了）
FACES_Data=[];
Faces_Ele_Sum=[];
for k=1:size(c,1)
    Node_k=c{k};
    xyzDraw=v(Node_k,:);
    CH_k = convhulln(xyzDraw) ;%convhulln输出的三角形方向均向内
    %根据CH_k索引，找到对应的点
    FACES_Data_k=Node_k(CH_k);
    Faces_Ele=ones(size(CH_k,1),1)*k;%标记这个面属于哪个多边形
    FACES_Data=[FACES_Data;FACES_Data_k];%储存所有面
    Faces_Ele_Sum=[Faces_Ele_Sum;Faces_Ele];
end
%删除重复面
FACES_Data2=sort(FACES_Data,2);%因为是三角形，所以即使改变顺序也是一圈方向。多边形不行。
[FACES_Data3,ind_F1,~]=unique(FACES_Data2,'rows','stable');%排序且去重复的Data
FACES_Data2=FACES_Data(ind_F1,:);%再按照不sort的顺序，输出一个去重复但不排序的Data
%每个面的Node数量
FACES_NODE_NUM=3*ones(1,size(FACES_Data3,1));%都是三角形，所以每个面上点的数量都是3
%然后是判断左右。三维需要按照FACES_Data给出的顺序，然后右手定则判断，先不判断。
%不判断左右也可以，能输出，但是会影响光照渲染。
NF=size(FACES_Data3,1);
Left_Elements=zeros(1,NF);
Right_Elements=Faces_Ele_Sum(ind_F1)';%先默认都在右

%再看看这个面是否在其它多面体中出现
for k=1:NF
    ChangeLR=0;
    FACES_k=FACES_Data2(k,:);
    indLR=find(all([any(FACES_Data==FACES_k(1),2),any(FACES_Data==FACES_k(2),2),any(FACES_Data==FACES_k(3),2)],2));
    %indLR=find(and(and(FACES_Data(:,1)==FACES_k(1),FACES_Data(:,2)==FACES_k(2)),FACES_Data(:,3)==FACES_k(3)));
    if numel(indLR)==1
        %说明这个面在外面，只对应一个体
        te=indLR(1);%这个单元体
        if Find_Line3(FACES_k,FACES_Data(te,:))==2
            ChangeLR=1;%如果当前FACES_k与实际FACES_Data反向，说明待会需要左右面交换一下
        end
        continue
    elseif numel(indLR)==2
        te=indLR(2);
        Left_Elements(k)=Faces_Ele_Sum(te);
        te=indLR(1);%
        if Find_Line3(FACES_k,FACES_Data(te,:))==2
            ChangeLR=1;%如果当前FACES_k与实际FACES_Data反向，说明待会需要左右面交换一下
        end
    end

    if ChangeLR
        Lt=Left_Elements(k);
        Rt=Right_Elements(k);
        Left_Elements(k)=Rt;
        Right_Elements(k)=Lt;
    end
end
FElement_Data.N_Element=size(c,1);
FElement_Data.Left=Left_Elements;
FElement_Data.Right=Right_Elements;


%保存简化的FACES_Data
FACES_Data=FACES_Data2;
%整理一下，最好第一列上下不要重复（玄学）
for k=2:size(FACES_Data,1)
    FACES_Data_k=FACES_Data(k,:);
    if FACES_Data_k(1)==FACES_Data(k-1,1)
        FACES_Data_k=[FACES_Data_k(2:end),FACES_Data_k(1)];
        FACES_Data(k,:)=FACES_Data_k;
    end
end
%% 2数据输出
%1准备数据
filename=['example_3D_Surface_UN.plt'];
title='';%无标题
variables={'X','Y','Z','P1'};%和前面Mat_Data对应
zone_title='';%无标题

time=[];
%2创建文件
f_id=fopen(filename,'w');
fclose(f_id);

%3创建文件头
plt_Head(filename,title,variables)
%4创建zone（point）格式
plt_Zone_UN3(filename,zone_title,time,Mat_Data,FACES_Data,FACES_NODE_NUM,FElement_Data)




%绘图
hold on
for k=1:size(c,1)
    xyzDraw=v(c{k},:);
    CH_k = convhulln(xyzDraw) ;
    trisurf(CH_k,xyzDraw(:,1),xyzDraw(:,2),xyzDraw(:,3));
end
hold off
view(3)
lightangle(gca,-45,30)
lighting gouraud

function plt_Head(filename,title,variables)
%创建文件头
f_id=fopen(filename,'a');
%名称
if ~isempty(title)
    s=['TITLE = "',title,'"'];
    fprintf(f_id,'%s \r\n',s);
end
%变量
v=numel(variables);
s='VARIABLES = ';
for k=1:v
    if k~=1
        s=[s,','];
    end
    s=[s,' "',variables{k},'"'];
end
fprintf(f_id,'%s \r\n',s);
fclose(f_id);
end


function plt_Zone_UN3(filename,zone_title,time,Nodes_Data,FACES_Data,FACES_NODE_NUM,FElement_Data)
%创建zone,block格式（ZONETYPE=FEPolyhedron时必须采用这个形式）
f_id=fopen(filename,'a');
N=size(Nodes_Data,1);
F=size(FACES_Data,1);
E=FElement_Data.N_Element;
Left_Elements=FElement_Data.Left;
Right_Elements=FElement_Data.Right;
s='ZONE';
fprintf(f_id,'%s \r\n',s);
s=['NODES=',num2str(N),', FACES=',num2str(F),', ELEMENTS=',num2str(E),', DATAPACKING=BLOCK, ZONETYPE=FEPolyhedron'];
fprintf(f_id,'%s \r\n',s);
s=['TotalNumFaceNodes=',num2str(numel(FACES_Data))];%所有FACES_Data包含的点的数量。也等于sum(FACES_NODE_NUM)
fprintf(f_id,'%s \r\n',s);
s=['NumConnectedBoundaryFaces=0'];
fprintf(f_id,'%s \r\n',s);
s=['TotalNumBoundaryConnections=0'];
fprintf(f_id,'%s \r\n',s);
%标题
if ~isempty(zone_title)
    s=['t="',zone_title,'"'];
    fprintf(f_id,'%s \r\n',s); 
end

%定义非定常时间
if ~isempty(time)
    s=['SOLUTIONTIME=',num2str(time)];
    fprintf(f_id,'%s \r\n',s);
end

%变成BLOCk格式
Nodes_Data=Nodes_Data';
%导入数据 
for k=1:size(Nodes_Data,1)
    fprintf(f_id,'%s \r\n',['# 变量',num2str(k)]);
    Mat_Data_k=Nodes_Data(k,:);
    if N>10
        %为了防止太长导致读取失败，把数据再拆成每10个一行
        N10=fix(numel(Mat_Data_k)/10);
        Mat_2=zeros(10,N10);
        Mat_2(:)=Mat_Data_k(1:N10*10);
        Mat_2=Mat_2';
        for k2=1:N10
            fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
        end
        fprintf(f_id,'%s \r\n',num2str(Mat_Data_k(N10*10+1:end)));
    else
        fprintf(f_id,'%s \r\n',num2str(Mat_Data_k));
    end
end

%导入每个面分别有多少点组成
fprintf(f_id,'%s \r\n','# node count per face 每个面分别有多少点组成');
%为了防止太长导致读取失败，把数据再拆成每10个一行
if numel(FACES_NODE_NUM)>10
    N10=fix(numel(FACES_NODE_NUM)/10);
    Mat_2=zeros(10,N10);
    Mat_2(:)=FACES_NODE_NUM(1:N10*10);
    Mat_2=Mat_2';
    for k2=1:N10
        fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
    end
    fprintf(f_id,'%s \r\n',num2str(FACES_NODE_NUM(N10*10+1:end)));
else
    fprintf(f_id,'%s \r\n',num2str(FACES_NODE_NUM));
end

%导入Face
fprintf(f_id,'%s \r\n','# Faces 多面体所有面');
for k=1:size(FACES_Data,1)
    fprintf(f_id,'%s \r\n',num2str(FACES_Data(k,:)));
end

%导入left elements（在这个face左侧的单元）
fprintf(f_id,'%s \r\n','# left elements（在这个face左侧的单元）');
%为了防止太长导致读取失败，把数据再拆成每10个一行
if numel(Left_Elements)>10
    N10=fix(numel(Left_Elements)/10);
    Mat_2=zeros(10,N10);
    Mat_2(:)=Left_Elements(1:N10*10);
    Mat_2=Mat_2';
    for k2=1:N10
        fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
    end
    fprintf(f_id,'%s \r\n',num2str(Left_Elements(N10*10+1:end)));
else
    fprintf(f_id,'%s \r\n',num2str(Left_Elements));
end

%导入right elements（在这个face右侧的单元）
fprintf(f_id,'%s \r\n','# right elements（在这个face右侧的单元）');
%为了防止太长导致读取失败，把数据再拆成每10个一行
if numel(Right_Elements)>10
    N10=fix(numel(Right_Elements)/10);
    Mat_2=zeros(10,N10);
    Mat_2(:)=Right_Elements(1:N10*10);
    Mat_2=Mat_2';
    for k2=1:N10
        fprintf(f_id,'%s \r\n',num2str(Mat_2(k2,:)));
    end
    fprintf(f_id,'%s \r\n',num2str(Right_Elements(N10*10+1:end)));
else
    fprintf(f_id,'%s \r\n',num2str(Right_Elements));
end
fclose(f_id);
end

function IfInList=Find_Line3(ABC,List)
%查找ABC三角形是否在List中，如果在，输出在左或在右
%0不在，1在左，2在右
IfInList=0;%默认是0
List=[0,List,List(1),List(2),0];
List2=fliplr(List);

indx1=find(List==ABC(1));
if ~isempty(indx1)
    %方向相同
    for k=1:numel(indx1)
        if List(indx1(k)+1)==ABC(2) && List(indx1(k)+2)==ABC(3)
            IfInList=1;%找到相同ABC，输出
            return
        end
    end
end

indx2=find(List2==ABC(1));
if ~isempty(indx2)
    %方向相反
    for k=1:numel(indx2)
        if List2(indx2(k)+1)==ABC(2) && List2(indx2(k)+2)==ABC(3)
            IfInList=2;%找到CBA，输出
            return
        end
    end
end
end