IPCS-driven-cavity-/write_vtu_hex_vector.m at main · JD63021/IPCS-driven-cavity- · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
function write_vtu_hex_vector(xv, yv, zv, Ux, Uy, Uz, P, filename)
nx = numel(xv); ny = numel(yv); nz = numel(zv);
Nn = nx*ny*nz; Nx = nx-1; Ny = ny-1; Nz = nz-1; Nc = Nx*Ny*Nz;
[XI,YI,ZI] = ndgrid(xv,yv,zv); Xp=XI(:); Yp=YI(:); Zp=ZI(:);
Uxv=Ux(:); Uyv=Uy(:); Uzv=Uz(:); Pv=P(:); Umag=sqrt(Uxv.^2+Uyv.^2+Uzv.^2);
pid = @(i,j,k) i + (j-1)*nx + (k-1)*nx*ny;
conn=zeros(Nc,8,'int32'); offsets=int32((1:Nc).'*8); types=repmat(uint8(12),Nc,1);
c=0; for k=1:Nz, for j=1:Ny, for i=1:Nx
    p000=pid(i,j,k); p100=pid(i+1,j,k); p110=pid(i+1,j+1,k); p010=pid(i,j+1,k);
    p001=pid(i,j,k+1); p101=pid(i+1,j,k+1); p111=pid(i+1,j+1,k+1); p011=pid(i,j+1,k+1);
    c=c+1; conn(c,:)=int32([p000 p100 p110 p010 p001 p101 p111 p011]-1);
end, end, end
fid=fopen(filename,'w'); if fid<0, error('Cannot open %s',filename); end
fprintf(fid,'<?xml version="1.0"?>\n<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">\n');
fprintf(fid,'  <UnstructuredGrid>\n    <Piece NumberOfPoints="%d" NumberOfCells="%d">\n',Nn,Nc);
fprintf(fid,'      <Points>\n        <DataArray type="Float64" NumberOfComponents="3" format="ascii">\n');
for p=1:Nn, fprintf(fid,'          %.16g %.16g %.16g\n',Xp(p),Yp(p),Zp(p)); end
fprintf(fid,'        </DataArray>\n      </Points>\n');
fprintf(fid,'      <Cells>\n        <DataArray type="Int32" Name="connectivity" format="ascii">\n');
for c=1:Nc, fprintf(fid,'          %d %d %d %d %d %d %d %d\n',conn(c,:)); end
fprintf(fid,'        </DataArray>\n        <DataArray type="Int32" Name="offsets" format="ascii">\n');
for c=1:Nc, fprintf(fid,'          %d\n',offsets(c)); end
fprintf(fid,'        </DataArray>\n        <DataArray type="UInt8" Name="types" format="ascii">\n');
for c=1:Nc, fprintf(fid,'          %u\n',types(c)); end
fprintf(fid,'        </DataArray>\n      </Cells>\n');
fprintf(fid,'      <PointData Scalars="Umag" Vectors="U">\n');
fprintf(fid,'        <DataArray type="Float64" Name="U" NumberOfComponents="3" format="ascii">\n');
for p=1:Nn, fprintf(fid,'          %.16g %.16g %.16g\n',Uxv(p),Uyv(p),Uzv(p)); end
fprintf(fid,'        </DataArray>\n        <DataArray type="Float64" Name="Umag" format="ascii">\n');
for p=1:Nn, fprintf(fid,'          %.16g\n',Umag(p)); end
fprintf(fid,'        </DataArray>\n        <DataArray type="Float64" Name="P" format="ascii">\n');
for p=1:Nn, fprintf(fid,'          %.16g\n',Pv(p)); end
fprintf(fid,'        </DataArray>\n      </PointData>\n');
fprintf(fid,'    </Piece>\n  </UnstructuredGrid>\n</VTKFile>\n'); fclose(fid);
end