navier2d.m

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

% 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(gradient,:)  ,'vertices',data.mesh.p,'facecolor','none','edgecolor','y'); 

% 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           , ...
                'gradient'  ,gradient);

% Set GUI data
set(gcf,'UserData',bcdata);

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bc_type_tracer(varargin)

% Set the BC type and value

% Get GUI data from the "set bc" window
fig    = gcf;
bcdata = get(fig,'UserData');

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

if sum(in)==0
    errordlg('No edges selected.')
    return
end

% Get main GUI data
data = get(figure(1),'Userdata'); figure(fig);

% Prompt with listbox
[i,ok] = listdlg('PromptString' , 'BC type:'            , ...
                 'Name'         , 'Boundary Conditions' , ...
                 'SelectionMode', 'Single'              , ...
                 'ListString'   , {'S','dS/dn = 0'});
             
if ok
   if i==1      % Value
       % Prompt for value
       value = inputdlg('Tracer value','Set tracer value',1,{'0'});
       if ~isempty(value)
           % Assign to main GUI data
           data.bc(be(in),7:8) = repmat([1,str2double(value)],sum(in),1);
       else
           return
       end
   else         % Gradient
       % Assign to main GUI data
       data.bc(be(in),7:8) = repmat([0,0],sum(in),1); 
   end
end

% Re-evaluate
boundary     = false(size(e,1),1);
pressure     = boundary;            % "Pressure" type doesn't exist for tracer. Set = false
boundary(be) = true;
unassigned   = data.bc(:,7)==-1 & boundary;
velocity     = data.bc(:,7)== 1;    % Re-use "velocity" as value so that the other
gradient     = data.bc(:,7)== 0;    % sub-functions will work...  

% Clear selection
in                = false(size(pm,1),1);
bcdata.in         = in;
bcdata.unassigned = unassigned;
bcdata.velocity   = velocity;
bcdata.pressure   = pressure;
bcdata.gradient   = gradient;

% Plot midpoints
hold off

% Plot midpoints
plot(pe(unassigned,1),pe(unassigned,2),'k.', ...
     pe(velocity,1)  ,pe(velocity,2)  ,'b.', ...
     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(gradient,:)  ,'vertices',data.mesh.p,'facecolor','none','edgecolor','y');  
 
set(fig,'UserData',bcdata);

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

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function set_anim(varargin)

% Prompt for animation variables

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

list = {'U velocity','V velocity','Total velocity','Pressure','Vorticity','Tracer'};

% Prompt with listbox
[i,ok] = listdlg('PromptString', 'Select variables: (CTRL+click)', ...
                 'Name'        , 'Animation'                     , ...
                 'ListString'  , list);
       
if ok
    if length(i)>2
        errordlg('Can''t select more than 2 variables.')
        return
    else
        if length(i)==0
            data.animation(1) = 0;
        else
            data.animation(1) = i(1);
        end
        if length(i)==2
            data.animation(2) = i(2);
        else
            data.animation(2) = 0;
        end
    end
end

% Save GUI data
set(figure(1),'Userdata',data);
            
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function residuals(varargin)

% Select residual type

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

list = {'[U,V] residuals','[x,y] forces'};

% Prompt with listbox
[i,ok] = listdlg('PromptString' , 'Select type', ...
                 'Name'         , 'Monitors'   , ...
                 'ListString'   , list         , ...
                 'SelectionMode','Single');

if ok
    if i==1     % Residuals
        % Un-mark any edges for lift/drag
        data.flag = false(size(data.flag));
        % Reset main GUI data
        set(figure(1),'UserData',data);
    else        % Lift/drag
        lift_drag;
    end
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function lift_drag(varargin)

% Mark edges at which the lift and drag forces will be calculated

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

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

set(figure, ...
          'Name'        ,'Lift/drag Forces', ...
          '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','Mark'                , ...
               'Callback',@set_lift_drag);   
           
b3 = uicontrol('Style','PushButton'           , ...
               'Units','Normalized'           , ...
               'Position',[0.05,0.2,btnx,btny], ...
               'String','Clear'               , ...
               'Callback',@clear_sel);           
           
% 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 = Marked'      , ...
          '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);
pressure     = boundary;            % "Pressure" type doesn't exist. Set = false
gradient     = boundary;            % "Gradient" type doesn't exist. Set = false
boundary(be) = true;
unassigned   = data.flag==0 & boundary;
velocity     = data.flag==1;        % Re-use "velocity" as value so that the other
                                    % sub-functions will work...  

% Plot midpoints
plot(pe(unassigned,1),pe(unassigned,2),'k.', ...
     pe(velocity,1)  ,pe(velocity,2)  ,'b.'), 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'); 

% 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           , ...
                'gradient'  ,gradient);

% Set GUI data
set(gcf,'UserData',bcdata);

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function set_lift_drag(varargin)

% Record marked edges (Lift/drag calc) in the GUI data

% Get GUI data from the "set bc" window
fig    = gcf;
bcdata = get(fig,'UserData');

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

if sum(in)==0
    errordlg('No edges selected.')
    return
end

% Get main GUI data
data = get(figure(1),'Userdata'); figure(fig);

% Set main GUI data
data.flag(be) = in;

boundary     = false(size(e,1),1);
pressure     = boundary;            % "Pressure" type doesn't exist. Set = false
gradient     = boundary;            % "Gradient" type doesn't exist. Set = false
boundary(be) = true;
unassigned   = data.flag==0 & boundary;
velocity     = data.flag==1;        % Re-use "velocity" as value so that the other
                                    % sub-functions will work... 
% Clear selection
in                = false(size(pm,1),1);
bcdata.in         = in;
bcdata.unassigned = unassigned;
bcdata.velocity   = velocity;
bcdata.pressure   = pressure;
bcdata.gradient   = gradient;

% Plot midpoints
hold off

% Plot midpoints
plot(pe(unassigned,1),pe(unassigned,2),'k.', ...
     pe(velocity,1)  ,pe(velocity,2)  ,'b.'), 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'); 
 
set(fig,'UserData',bcdata);

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

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,t] = fixmesh(p,t)

% Remove duplicate nodes & setup CCW (Counter-clockwise) node ordering.

% Darren Engwirda - 2006. Adapted from Distmesh.

% Remove duplicate nodes
snap        = max(max(p,[],1)-min(p,[],1),[],2)*1024*eps;
[foo,ix,jx] = unique(round(p/snap)*snap,'rows');
p           = p(ix,:);
t           = jx(t);
[pix,ix,jx] = unique(t);
t           = reshape(jx,size(t));
p           = p(pix,:);

% CCW ordering
p1  = p(t(:,1),:); 
d12 = p(t(:,2),:)-p1;
d13 = p(t(:,3),:)-p1;

% Negative area?
flip          = (d12(:,1).*d13(:,2)-d12(:,2).*d13(:,1))<0;
t(flip,[1,2]) = t(flip,[2,1]);

return