From 75fb7b55fd90a61b32b86ba3d50bb1390622ea24 Mon Sep 17 00:00:00 2001 From: Lucas C Wilcox Date: Fri, 3 Apr 2015 11:41:33 -0700 Subject: [PATCH 1/3] Fix y normals in CylBC2D --- Codes1.1/CFD2D/CylBC2D.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Codes1.1/CFD2D/CylBC2D.m b/Codes1.1/CFD2D/CylBC2D.m index 599bc44..3a69632 100644 --- a/Codes1.1/CFD2D/CylBC2D.m +++ b/Codes1.1/CFD2D/CylBC2D.m @@ -26,7 +26,7 @@ % 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; @@ -34,7 +34,7 @@ % 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; From 1d1345907d7f17be04e65b508ff8a52f35f3554c Mon Sep 17 00:00:00 2001 From: Lucas C Wilcox Date: Fri, 3 Apr 2015 11:44:32 -0700 Subject: [PATCH 2/3] Correct return format for CylIC2D --- Codes1.1/CFD2D/CylIC2D.m | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Codes1.1/CFD2D/CylIC2D.m b/Codes1.1/CFD2D/CylIC2D.m index 59a3b82..6b0cd9b 100644 --- a/Codes1.1/CFD2D/CylIC2D.m +++ b/Codes1.1/CFD2D/CylIC2D.m @@ -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 @@ -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 From 553fc90a58cb734271a2546955ba3ea59ca9067c Mon Sep 17 00:00:00 2001 From: Lucas C Wilcox Date: Fri, 3 Apr 2015 11:46:45 -0700 Subject: [PATCH 3/3] Add VTK output functions --- Codes1.1/ServiceRoutines/WriteVTK2D.m | 72 +++++++++++++++ Codes1.1/ServiceRoutines/WriteVTK3D.m | 121 ++++++++++++++++++++++++++ 2 files changed, 193 insertions(+) create mode 100644 Codes1.1/ServiceRoutines/WriteVTK2D.m create mode 100644 Codes1.1/ServiceRoutines/WriteVTK3D.m diff --git a/Codes1.1/ServiceRoutines/WriteVTK2D.m b/Codes1.1/ServiceRoutines/WriteVTK2D.m new file mode 100644 index 0000000..34d938b --- /dev/null +++ b/Codes1.1/ServiceRoutines/WriteVTK2D.m @@ -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 diff --git a/Codes1.1/ServiceRoutines/WriteVTK3D.m b/Codes1.1/ServiceRoutines/WriteVTK3D.m new file mode 100644 index 0000000..eea8587 --- /dev/null +++ b/Codes1.1/ServiceRoutines/WriteVTK3D.m @@ -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