如何将自编程有限元程序求解的场结果转化为vtk文件,以便在Paraview中显示一直是困扰我的难题。
前两天恶补vtk的语法知识,终于攻克了这个难题,之前一直就是套用已有代码,云里雾里,懂了语法结构后,就可以随心所欲的使用各种代码语言实现它,就以最简单的Matlab语言为例,今天分享一下使用心得,希望对该话题感兴趣的小伙伴有帮助!
官方帮助文档:
https://docs.vtk.org/en/latest/
头文件:包含了 VTK 文件的版本信息和描述
# vtk DataFile Version x.x
Description
ASCII
DATASET data_type
数据集描述:此部分指定了数据集的类型和大小。不同类型的数据集需要不同的描述。
% n 节点数, dataType为double float等数据类型
POINTS [n] [dataType]
[x_0] [y_0] [z_0]
[x_1] [y_1] [z_1]
[x_0] [y_0] [z_0]
CELLS [n_cells] [n_list] % 单元个数、数据行占的格子
[单元0上的顶点总数目] [单元0顶点0的编号] [单元0顶点1的编号] [单元0顶点2的编号]
[单元1上的顶点总数目] [单元1顶点0的编号] [单元1顶点1的编号] [单元1顶点2的编号]
CELL_TYPES [n_cells] % 单元个数
[单元0类型] % 表示几何类型的整数
数据属性:这部分包括与数据相关的属性,例如,标量、矢量、张量等。这些属性可以随节点或单元关联。
POINT_DATA [n] % n表示节点数
CELL_DATA [n_cells] % n_cells表示单元数
% dataName数据名字,dataType是double之类的,numComp表示每个节点有两个分量
SCALARS [dataName] [dataType] [numComp]
load('out.mat');
vtk_file = fopen('output.vtk', 'w');
node_coords =Node;
element_info =Element;
displacement_info =U;
AxialForce_info = AxialForce';
% 写入 VTK 文件头
fprintf(vtk_file, '# vtk DataFile Version 3.0\n');
fprintf(vtk_file, 'Finite Element Results\n');
fprintf(vtk_file, 'ASCII\n');
fprintf(vtk_file, 'DATASET UNSTRUCTURED_GRID\n');
% 写入节点坐标
num_nodes = size(node_coords, 1);
fprintf(vtk_file, 'POINTS %d float\n', num_nodes);
for i = 1:num_nodes
fprintf(vtk_file, '%f %f 0.0\n', node_coords(i, 1), node_coords(i, 2));
end
% 写入单元信息
num_elements = size(element_info, 1);
fprintf(vtk_file, 'CELLS %d %d\n', num_elements, num_elements * 3);
for i = 1:num_elements
fprintf(vtk_file, '2 %d %d\n', element_info(i, 2) - 1, element_info(i, 3) - 1);
end
% 写入单元类型
fprintf(vtk_file, 'CELL_TYPES %d\n', num_elements);
for i = 1:num_elements
fprintf(vtk_file, '3\n'); % 2 节点线元素
end
%%%%%%%%%%%%%%%%%%%%%%%%POINT_DAT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 写入位移数据
fprintf(vtk_file, 'POINT_DATA %d\n', num_nodes);
fprintf(vtk_file, 'SCALARS Displacement float 2\n');
fprintf(vtk_file, 'LOOKUP_TABLE default\n');
for i = 1:num_nodes
% fprintf(vtk_file, '%f %f\n', displacement_info(i, 2), displacement_info(i, 3));
fprintf(vtk_file, '%f %f \n', displacement_info(i * 2 - 1), displacement_info(i * 2));
end
% 写入节点轴力数据
% fprintf(vtk_file, 'POINT_DATA %d\n', num_nodes);
fprintf(vtk_file, 'SCALARS AxialForce float 2\n');
fprintf(vtk_file, 'LOOKUP_TABLE default\n');
for i = 1:num_nodes
fprintf(vtk_file, '%f %f \n', AxialForce_info(i * 2 - 1), AxialForce_info(i * 2));
end
fclose(vtk_file);
Paraview显示效果如下: