Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Minor changes and VTK output #1

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Minor changes and VTK output #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
75fb7b5
Fix y normals in CylBC2D
lcw Apr 3, 2015
1d13459
Correct return format for CylIC2D
lcw Apr 3, 2015
553fc90
Add VTK output functions
lcw Apr 3, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Codes1.1/CFD2D/CylBC2D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,15 @@
% Wall conditions -- reflective, isothermal, i.e., n.u=0, T=T(t=0)
Temp = pin/rhoin/(gamma-1);
rhoW = rho(gmapW); rhouW = rhou(gmapW); rhovW = rhov(gmapW); EnerW = Ener(gmapW);
nxW = gnx(gmapW); nyW = nyin(gmapW);
nxW = gnx(gmapW); nyW = gny(gmapW);

rhou(gmapW) = -rhou(gmapW); rhov(gmapW) = -rhov(gmapW);
Ener(gmapW) = rhoW*Temp + 0.5*(rhou(gmapW).^2 + rhov(gmapW).^2)./rhoW;

% cylinder conditions -- reflective, isothermal, i.e., n.u=0, T=T(t=0)
Temp = pin/rhoin/(gamma-1);
rhoC = rho(gmapC); rhouC = rhou(gmapC); rhovC = rhov(gmapC); EnerC = Ener(gmapC);
nxC = gnx(gmapC); nyC = nyin(gmapC);
nxC = gnx(gmapC); nyC = gny(gmapC);

rhou(gmapC) = -rhou(gmapC); rhov(gmapC) = -rhov(gmapC);
Ener(gmapC) = rhoC*Temp + 0.5*(rhou(gmapC).^2 + rhov(gmapC).^2)./rhoC;
Expand Down
12 changes: 6 additions & 6 deletions Codes1.1/CFD2D/CylIC2D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
function Q = ChannelIC2D(x, y, time);
function [rho,rhou,rhov,Ener] = CylIC2D(x, y, time);

% function [Q] = ChannelIC2D(x, y, time)
% function [rho,rhou,rhov,Ener] = CylIC2D(x, y, time)
% Purpose: Impose uniform plane flow

% Example is Mach ** 0.4 ** flow in wind tunnel
Expand All @@ -12,8 +12,8 @@

% pack modified conserved variables

Q(:,:,1) = rhoin*ones(size(x));
Q(:,:,2)= rhoin*(1/.41)^2*6*(y+.2).*(0.41 - (y+.2));
Q(:,:,3) = 0;
Q(:,:,4) = Ein + 0.5*(Q(:,:,2).^2 + Q(:,:,3).^2)./Q(:,:,1);
rho = rhoin*ones(size(x));
rhou = rhoin*(1/.41)^2*6*(y+.2).*(0.41 - (y+.2));
rhov = 0;
Ener = Ein + 0.5*(rhou.^2 + rhov.^2)./rho;
return
72 changes: 72 additions & 0 deletions Codes1.1/ServiceRoutines/WriteVTK2D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
function WriteVTK2D(filename, Nout, vnames, varargin)
nfields = nargin - 3;

Globals2D;

% build equally spaced grid on reference triangle
Npout = (Nout+1)*(Nout+2)/2;
rout = zeros(Npout,1); sout = zeros(Npout,1);
sk = 1;
for n=1:Nout+1
for m=1:Nout+2-n
rout(sk) = -1 + 2*(m-1)/Nout;
sout(sk) = -1 + 2*(n-1)/Nout;
counter(n,m) = sk; sk = sk+1;
end
end

% build matrix to interpolate field data to equally spaced nodes
interp = InterpMatrix2D(rout, sout);

% build triangulation of equally spaced nodes on reference triangle
tri = [];
for n=1:Nout+1
for m=1:Nout+1-n,
v1 = counter(n,m); v2 = counter(n,m+1);
v3 = counter(n+1,m); v4 = counter(n+1,m+1);
if(v4)
tri = [tri;[[v1 v2 v3];[v2 v4 v3]]];
else
tri = [tri;[[v1 v2 v3]]];
end
end
end

% build triangulation for all equally spaced nodes on all elements
TRI = [];
for k=1:K
TRI = [TRI; tri+(k-1)*Npout];
end

% interpolate node coordinates and field to equally spaced nodes
xout = interp*x;
yout = interp*y;
zout = zeros(size(xout));

uout = cell(1, nfields);
for n=1:nfields
uout{n} = interp*varargin{n};
end

Ntotal = length(xout(:));
[nTRI, nVERT] = size(TRI);

fid = fopen(filename, 'w');
fprintf(fid, '# vtk DataFile Version 2');
fprintf(fid, '\nNUDG simulation');
fprintf(fid, '\nASCII');
fprintf(fid, '\nDATASET UNSTRUCTURED_GRID\n');
fprintf(fid, '\nPOINTS %d double', Ntotal);
fprintf(fid, '\n%25.16e %25.16e %25.16e', [xout(:) yout(:) zout(:)]');
fprintf(fid, '\nCELLS %d %d', nTRI, nTRI*4);
fprintf(fid, '\n3 %10d %10d %10d', (TRI-1)');
fprintf(fid, '\nCELL_TYPES %d', nTRI);
fprintf(fid, '\n%d', repmat([5], nTRI, 1));
fprintf(fid, '\nPOINT_DATA %d', Ntotal);
for n=1:nfields
fprintf(fid, '\nSCALARS %s double 1', vnames{n});
fprintf(fid, '\nLOOKUP_TABLE default');
fprintf(fid, '\n%25.16e', uout{n});
end
fclose(fid);
end
121 changes: 121 additions & 0 deletions Codes1.1/ServiceRoutines/WriteVTK3D.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
function WriteVTK3D(filename, Nout, vnames, varargin)
nfields = nargin - 3;

Globals3D;

Npout = (Nout+1)*(Nout+2)*(Nout+3)/6;
[rout,sout,tout] = EquiNodes3D(Nout);


% build matrix to interpolate field data to equally spaced nodes
Vout = Vandermonde3D(N, rout, sout, tout);
interp = Vout*invV;


% form symmetric tetrahedralization of local nodes
startrow = zeros(Nout+1, Nout+1);

sk = 1;
for i=0:Nout
for j=0:Nout-i
startrow(j+1, i+1) = sk;
sk = sk + Nout+1-i-j;
end
end

% contruct tetrahedralization
tet = zeros(Nout*Nout*Nout,4);

sk = 1;
for i=0:Nout
for j=0:Nout-i
for k=0:Nout-i-j-1
% Add Tet 1
tet(sk,1) = startrow(j+1, i+1)+k;
tet(sk,2) = startrow(j+1, i+1)+k+1;
tet(sk,3) = startrow(j+2, i+1)+k;
tet(sk,4) = startrow(j+1, i+2)+k;
sk = sk+1;

if(k < Nout-i-j-1)
% Add Tet 2
tet(sk,1) = startrow(j+1, i+1)+k+1;
tet(sk,2) = startrow(j+2, i+1)+k;
tet(sk,3) = startrow(j+1, i+2)+k;
tet(sk,4) = startrow(j+1, i+2)+k+1;
sk = sk+1;

% Add Tet 3
tet(sk,1) = startrow(j+1, i+1)+k+1;
tet(sk,2) = startrow(j+2, i+1)+k+1;
tet(sk,3) = startrow(j+2, i+1)+k;
tet(sk,4) = startrow(j+1, i+2)+k+1;
sk = sk+1;

% Add Tet 4
tet(sk,1) = startrow(j+1, i+2)+k;
tet(sk,2) = startrow(j+2, i+2)+k;
tet(sk,3) = startrow(j+1, i+2)+k+1;
tet(sk,4) = startrow(j+2, i+1)+k;
sk = sk+1;

% Add Tet 5
tet(sk,1) = startrow(j+2, i+1)+k;
tet(sk,2) = startrow(j+2, i+1)+k+1;
tet(sk,3) = startrow(j+2, i+2)+k;
tet(sk,4) = startrow(j+1, i+2)+k+1;
sk = sk+1;
end

if(k < Nout-i-j-2)
% Add Tet 6
tet(sk,1) = startrow(j+1, i+2)+k+1;
tet(sk,2) = startrow(j+2, i+1)+k+1;
tet(sk,3) = startrow(j+2, i+2)+k;
tet(sk,4) = startrow(j+2, i+2)+k+1;
sk = sk+1;
end
end
end
end

Ntet = sk-1;

% create global tet mesh
TET = [];
for k=1:K
TET = [TET; tet+(k-1)*Npout];
end

% interpolate node coordinates and field to equally spaced nodes
xout = interp*x;
yout = interp*y;
zout = interp*z;

uout = cell(1, nfields);
for n=1:nfields
uout{n} = interp*varargin{n};
end

Ntotal = length(xout(:));
[nTET, nVERT] = size(TET);

fid = fopen(filename, 'w');
fprintf(fid, '# vtk DataFile Version 2');
fprintf(fid, '\nNUDG simulation');
fprintf(fid, '\nASCII');
fprintf(fid, '\nDATASET UNSTRUCTURED_GRID\n');
fprintf(fid, '\nPOINTS %d double', Ntotal);
fprintf(fid, '\n%25.16e %25.16e %25.16e', [xout(:) yout(:) zout(:)]');
fprintf(fid, '\nCELLS %d %d', nTET, nTET*5);
fprintf(fid, '\n4 %10d %10d %10d %10d', (TET-1)');
fprintf(fid, '\nCELL_TYPES %d', nTET);
fprintf(fid, '\n%d', repmat([10], nTET, 1));
fprintf(fid, '\nPOINT_DATA %d', Ntotal);
for n=1:nfields
fprintf(fid, '\nSCALARS %s double 1', vnames{n});
fprintf(fid, '\nLOOKUP_TABLE default');
fprintf(fid, '\n%25.16e', uout{n});
end
fclose(fid);
end