openMA_modes_interp.m

function [ux,uy] = openMA_modes_interp( cc, nm, dm, bm )
% OPENMA_MODES_INTERP - interpolates modes at a specific set of points
% and returns either the modes themselves or the velocities at those
% points. 
%
% Usage: u = openMA_modes_interp( cc, vfm, dfm, bm )
%        [ux,uy] = openMA_modes_interp( cc, vfm, dfm, bm )
%
% The first form returns the modes themselves at the coordinates in cc,
% which should be a two column matrix of x and y coordinates in the same
% units as the triangular grids of the modes themselves.
%
% The second form returns the velocities of the modes at the coordinates.
%
% NOTE that the number of output arguments is very important to the result.
%
% vfm, dfm, and bm are the vorticity-free modes (formerly known as the
% Neumann modes), divergence-free modes (formerly the Dirichlet modes) and
% the boundary modes (or harmonic modes), respectively.  These should be in
% the format generated by the openMA_pdetool_neumann_modes_solve,
% openMA_pdetool_dirichlet_modes_solve and
% openMA_pdetool_boundary_modes_solve.  If any of them is empty or absent,
% those modes will be ignored.
%
% The format of u, ux and uy will be an array with one row for each point
% in cc and one column for each mode in the order given.
%
% Interpolation of the modes is done via linear interpolation on the
% enclosing triangles.  Interpolation of the velocities is done by
% selecting the velocity of the enclosing triangle.  Both of these
% basically assume that the triangles are small enough that the linear
% approximation is valid over the triangle.  If this assumption isn't
% valid, you probably didn't do a very good job of solving your original
% PDE.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 	$Id: openMA_modes_interp.m 79 2007-03-05 21:51:20Z dmk $	
%
% Copyright (C) 2005 David M. Kaplan
% Licence: GPL (Gnu Public License)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('dm','var')
  dm = [];
end

if ~exist('bm','var')
  bm = [];
end

% Number of modes total.
nnn = length(nm) + length(dm) + length(bm);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do interpolation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout == 1
  % Want modes themselves interpolated
  
  ux = zeros(size(cc,1),nnn);
  
  nn = 0; % Start counter

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over neumann modes and do interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(nm)
    nn = nn + 1;
    ux(:,nn) = pdeintrp_arbitrary( cc, nm(k).p, nm(k).t, ...
				     nm(k).u );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over dirichlet modes and do interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(dm)
    nn = nn + 1;
    ux(:,nn) = pdeintrp_arbitrary( cc, dm(k).p, dm(k).t, ...
				     dm(k).u );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over boundary modes and do interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(bm)
    nn = nn + 1;
    ux(:,nn) = pdeintrp_arbitrary( cc, bm(k).p, bm(k).t, ...
				   bm(k).u );
  end
  
else 
  % Want vector velocities.
  % Will use pdeintrp_arbitrary.
  % This implicitly assumes that velocity is constant over each
  % triangle, as the pde solver would do.  Probably an acceptable
  % assumption in this application.
  
  [ux,uy] = deal( zeros(size(cc,1),nnn) );
  
  nn = 0; % Start counter

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over neumann modes and do interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(nm)
    nn = nn + 1;
    [ux(:,nn),uy(:,nn)] = pdeintrp_arbitrary( cc, nm(k).p, nm(k).t, ...
						  nm(k).u );
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over dirichlet modes and do interpolation
  %
  % Sign in definition of vorticity modes agrees with k x grad(psi)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(dm)
    nn = nn + 1;
    [uy(:,nn),ux(:,nn)] = pdeintrp_arbitrary( cc, dm(k).p, dm(k).t, ...
						  dm(k).u );
    
    % Fix sign of x term
    ux(:,nn) = -ux(:,nn);
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Loop over boundary modes and do interpolation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for k = 1:length(bm)
    nn = nn + 1;
    [ux(:,nn),uy(:,nn)] = pdeintrp_arbitrary( cc, bm(k).p, bm(k).t, ...
						  bm(k).u );
  end
  
end