navier2d.m

ZEUS is a family of Eulerian (grid based) Magneto-Hydrodynamic codes (MHD) for use in astrophysics,

% Show IC's
set(figure,'Name','Initial Conditions')
subplot(1,3,1), trimesh(data.mesh.t,x,y,U), title('U velocity'), axis square
subplot(1,3,2), trimesh(data.mesh.t,x,y,V), title('V velocity'), axis square
subplot(1,3,3), trimesh(data.mesh.t,x,y,S), title('Tracer')    , axis square

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ok = checkvar(expr)

% Helper for init

var = symvar(expr);     % Get variables in expr
num = numel(var);       % Number of variables

ok = false;
if num==0
    ok = true;
elseif num<=2
    for i = 1:num       % Check that expr only contains 'x' and/or 'y'
        if any(strcmp(var{i},{'x','y'}))
            ok = true;
        end
    end
end
if ~ok
    errordlg('Inputs must contain only constants or functions of [x,y]')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function run(varargin)

% Main run button

% Get GUI data
data = get(figure(1),'UserData');

% Check to see that everything has been specified
if isempty(data.settings)
    errordlg('Integration settings have not been specified','Error')
    return
end
if isempty(data.mesh)
    errordlg('Mesh file has not been loaded','Error')
    return
end
if isempty(data.init)
    errordlg('Initial conditions have not been set')
    return
end

% UVP BC check
e            = data.mesh.e;
be           = data.mesh.be;
boundary     = false(size(e,1),1);
boundary(be) = true;
unassigned   = data.bc(:,1)==-1 & boundary;

% Solver settings
solver = get(data.button(16),'Value');              % UVP, tracer, thermal
if solver>1
    unassigned = (data.bc(:,7)==-1 & boundary) | unassigned;
end
data.solver = solver;

if any(unassigned)
    errordlg('Boundary conditions have not been specified','Error')
    return
end

% Animation settings
data.animation(3) = get(data.button(11),'Value');   % Plot type
data.animation(4) = get(data.button(12),'Value');   % Plot view

% CALL NAVIER-STOKES SOLVER
tvd_rk2(data);      % 2 stage TVD Runge-Kutta solver 

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function help_main(varargin)

% Quick help for the main window

helptext = {'Quick help for Navier2d.'
            ''
            ['The main window is broken into four sections, corresponding ', ...
             'to the four main steps involved in setting up a simulation.']
            ''
            'Step 1: Load a mesh file'
            'Step 2: Specify the boundary conditions'
            'Step 3: Specify the integration settings and initial conditions'
            'Step 4: Specify the animation settings'
            ''
            'The "Run" button will then start the simulation.'
            ''
            ['Mesh files must be created separately. The mesh generator mesh2d ', ...
             'was developed for this purpose.']   
           };

uiwait(msgbox(helptext,'Help','none')); 

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function load_mesh(varargin)

% Prompt user to load in a mesh file.

uiload

% Get current GUI data
data = get(figure(1),'UserData');

% Error checking
if ~exist('p','var')||~exist('t','var')
    errordlg('Not a valid mesh file.','Error')
    return
end
if (size(p,2)~=2)||(size(t,2)~=3)
    errordlg('Incorrect dimensions for p or t.','Error')
    return
end
if (max(t(:))>size(p,1))||(min(t(:))<=0)
    errordlg('t is not a valid triangulation of the nodes in p.','Error')
    return
end

% Run fixmesh to make sure the mesh is CCW and non-duplicate
[p,t] = fixmesh(p,t);

% Setup the mesh based data structures
data.mesh      = data_structure(p,t);
data.bc        = repmat(-1,size(data.mesh.e,1),8);
data.init      = [];
data.settings  = [];
data.animation = [];
data.flag      = false(size(data.mesh.e,1),1);

% Set GUI data
set(figure(1),'UserData',data);

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_mesh(varargin)

% Show the current mesh

data = get(figure(1),'UserData');

if isempty(data.mesh)
    errordlg('No mesh file loaded.','Error')
    return
end

set(figure,'Name','Mesh'); 

patch('faces',data.mesh.t,'vertices',data.mesh.p,'facecolor','w','edgecolor','b'); 

axis equal, axis off

title([num2str(size(data.mesh.p,1)),' Nodes, ', ...
       num2str(size(data.mesh.t,1)),' Triangles'])

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function show_median(varargin)

% Show the median dual mesh used by NS.m

% Get GUI data
data = get(figure(1),'UserData');

if isempty(data.mesh)
    errordlg('No mesh file loaded.','Error')
    return
end

set(figure,'Name','Median Dual Mesh'); 

p    = data.mesh.p;     % Nodes
t    = data.mesh.t;     % Triangulation
pc   = data.mesh.pc;    % Centroids
n2n  = data.mesh.n2n;   % Node to node connectivity
n2e  = data.mesh.n2e;   % Node to edge connectivity
e2t  = data.mesh.e2t;   % Edge to triangle connectivity
numn = size(p,1);

w = waitbar(0,'Building median mesh'); drawnow

xline = repmat(0,2*size(n2n,2)*numn,2);
yline = xline;
next  = 1;
lim   = 0;
for k = 1:numn
    m = 1;
    n = k/numn;
    if (n>lim)||(n==1)
        waitbar(n,w); lim = n+0.1;
    end
    while n2n(k,m)>0    % Loop around neighbours of node k
        % Triangles associated with current edge
        t1 = e2t(n2e(k,m),1);
        t2 = e2t(n2e(k,m),2);
        % Edge midpoint
        mx = 0.5*(p(k,1)+p(n2n(k,m),1));
        my = 0.5*(p(k,2)+p(n2n(k,m),2));
        % Median edge 1 (joins pc1 and pm)
        xline(next,1) = mx;
        xline(next,2) = pc(t1,1);
        yline(next,1) = my;
        yline(next,2) = pc(t1,2);
        next          = next+1;
        % Median edge 2 (joins pc2 and pm)
        if t2>0
            xline(next,1) = mx;
            xline(next,2) = pc(t2,1);
            yline(next,1) = my;
            yline(next,2) = pc(t2,2);
            next          = next+1; 
        end
        % Counter
        m = m+1;
    end
end
close(w)

% Triangles
patch('faces',t,'vertices',p,'facecolor','none','edgecolor','w'); axis equal, axis off, hold on
% Median CV's
patch('xdata',xline','ydata',yline','facecolor','none','edgecolor','b'); 

title([num2str(size(data.mesh.p,1)),' Nodes, ', ...
       num2str(size(data.mesh.t,1)),' Triangles'])

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function set_bc(varargin)

% Set the boundary conditions

data = get(figure(1),'UserData');

if isempty(data.mesh)
    errordlg('No mesh file loaded.','Error')
    return
end

set(figure, ...
          'Name'        ,'Boundary Conditions', ...
          'DoubleBuffer','On'                 , ...
          'Units'       ,'Normalized');
      
axes('Units'   ,'Normalized', ...
     'Position',[0.25,0.1,0.7,0.8]);      
     
btnx = 0.1;
btny = 0.05;

% Buttons
b1 = uicontrol('Style','PushButton'           , ...
               'Units','Normalized'           , ...
               'Position',[0.05,0.3,btnx,btny], ...
               'String','Select'              , ...
               'Callback',@select);
           
b2 = uicontrol('Style','PushButton'           , ...
               'Units','Normalized'           , ...
               'Position',[0.05,0.1,btnx,btny], ...
               'String','Set'                 , ...
               'Callback',@bc_type);   
           
b3 = uicontrol('Style','PushButton'           , ...
               'Units','Normalized'           , ...
               'Position',[0.05,0.2,btnx,btny], ...
               'String','Clear'               , ...
               'Callback',@clear_sel);  
           
b4 = uicontrol('Style'   ,'PushButton'        , ...
               'Units'   ,'Normalized'        , ...
               'Position',[0.05,0.4,btnx,btny], ...
               'String'  ,'Help'              , ...
               'Callback',@help_bc);              
           
% Headers           
uicontrol('Style'          ,'Text'              , ...
          'Units'          ,'Normalized'        , ...
          'Position'       ,[0.025,0.8,0.2,0.05], ...
          'String'         ,'Black = Unassigned', ...
          'BackgroundColor',[0.8,0.8,0.8]);  
uicontrol('Style'          ,'Text'               , ...
          'Units'          ,'Normalized'         , ...
          'Position'       ,[0.025,0.75,0.2,0.05], ...
          'String'         ,'Blue = Velocity'    , ...
          'BackgroundColor',[0.8,0.8,0.8]);  
uicontrol('Style'          ,'Text'              , ...
          'Units'          ,'Normalized'        , ...
          'Position'       ,[0.025,0.7,0.2,0.05], ...
          'String'         ,'Green = Pressure'  , ...
          'BackgroundColor',[0.8,0.8,0.8]);  
uicontrol('Style'          ,'Text'               , ...
          'Units'          ,'Normalized'         , ...
          'Position'       ,[0.025,0.65,0.2,0.05], ...
          'String'         ,'Yellow = Gradient'  , ...
          'BackgroundColor',[0.8,0.8,0.8]);  

% Boundary edge geometry
e  = data.mesh.e;
be = data.mesh.be;
p  = data.mesh.p;
pe = data.mesh.pe;
pm = 0.5*(p(e(be,1),:)+p(e(be,2),:));

boundary     = false(size(e,1),1);
boundary(be) = true;
unassigned   = data.bc(:,1)==-1 & boundary;
velocity     = data.bc(:,1)== 1;
pressure     = data.bc(:,5) > 0;
gradient     = xor(data.bc(:,1),data.bc(:,3));

% Plot midpoints
plot(pe(unassigned,1),pe(unassigned,2),'k.', ...
     pe(velocity,1)  ,pe(velocity,2)  ,'b.', ...
     pe(pressure,1)  ,pe(pressure,2)  ,'g.', ...
     pe(gradient,1)  ,pe(gradient,2)  ,'y.'), axis equal, axis off, hold on

% Plot edges
patch('faces',e(unassigned,:),'vertices',data.mesh.p,'facecolor','none','edgecolor','k'); 
patch('faces',e(velocity,:)  ,'vertices',data.mesh.p,'facecolor','none','edgecolor','b'); 
patch('faces',e(pressure,:)  ,'vertices',data.mesh.p,'facecolor','none','edgecolor','g'); 
patch('faces',e(gradient,:)  ,'vertices',data.mesh.p,'facecolor','none','edgecolor','y'); 

% Plot arrows for velocity type
if any(velocity)
    quiver(pe(velocity,1),pe(velocity,2),data.bc(velocity,2),data.bc(velocity,4))
end


% New GUI data just for the "set bc" window
bcdata = struct('p'         ,p                  , ...
                'pe'        ,pe                 , ...
                'e'         ,e                  , ...
                'be'        ,be                 , ...
                'in'        ,false(size(pm,1),1), ...
                'unassigned',unassigned         , ...
                'velocity'  ,velocity           , ...
                'pressure'  ,pressure           , ...