div

 div - divergence operator

    fd = div(Px,Py, options);
    fd = div(P, options);

   options.bound = 'per' or 'sym'
   options.order = 1 (backward differences)
                 = 2 (centered differences)

   Note that the -div and grad operator are adjoint
   of each other such that 
       <grad(f),g>=<f,-div(g)>

   See also: grad.

    Copyright (c) 2007 Gabriel Peyre

0001 function fd = div(Px,Py, options)
0002 
0003 % div - divergence operator
0004 %
0005 %    fd = div(Px,Py, options);
0006 %    fd = div(P, options);
0007 %
0008 %   options.bound = 'per' or 'sym'
0009 %   options.order = 1 (backward differences)
0010 %                 = 2 (centered differences)
0011 %
0012 %   Note that the -div and grad operator are adjoint
0013 %   of each other such that
0014 %       <grad(f),g>=<f,-div(g)>
0015 %
0016 %   See also: grad.
0017 %
0018 %    Copyright (c) 2007 Gabriel Peyre
0019 
0020 % retrieve number of dimensions
0021 nbdims = 2;
0022 if size(Px,1)==1 || size(Px,2)==1
0023     nbdims = 1;
0024 end
0025 if size(Px,3)>1
0026     if nargin>1
0027         options = Py;
0028         clear Py;
0029     end
0030     if size(Px,4)<=1
0031         Py = Px(:,:,2);
0032         Px = Px(:,:,1);
0033     else
0034         Pz = Px(:,:,:,3);
0035         Py = Px(:,:,:,2);
0036         Px = Px(:,:,:,1); 
0037         nbdims = 3;
0038     end
0039 end
0040 
0041 options.null = 0;
0042 bound = getoptions(options, 'bound', 'sym');
0043 order = getoptions(options, 'order', 1);
0044 
0045 if strcmp(bound, 'sym')
0046     if order==1
0047         fx = Px-Px([1 1:end-1],:,:);         
0048         fx(1,:,:)   = Px(1,:,:);        % boundary
0049         fx(end,:,:) = -Px(end-1,:,:);        
0050         if nbdims>=2
0051             fy = Py-Py(:,[1 1:end-1],:);
0052             fy(:,1,:)   = Py(:,1,:);    % boundary
0053             fy(:,end,:) = -Py(:,end-1,:);
0054         end        
0055         if nbdims>=3
0056             fz = Pz-Pz(:,:,[1 1:end-1]);  
0057             fz(:,:,1)   = Pz(:,:,1);    % boundary
0058             fz(:,:,end) = -Pz(:,:,end-1);
0059         end
0060     else
0061         fx = (Px([2:end end],:,:)-Px([1 1:end-1],:,:))/2;
0062         fx(1,:,:)   = +Px(2,:,:)/2+Px(1,:,:);   % boundary
0063         fx(2,:,:)   = +Px(3,:,:)/2-Px(1,:,:);
0064         fx(end,:,:) = -Px(end,:,:)-Px(end-1,:,:)/2;
0065         fx(end-1,:,:) = Px(end,:,:)-Px(end-2,:,:)/2;
0066         if nbdims>=2
0067             fy = (Py(:,[2:end end],:)-Py(:,[1 1:end-1],:))/2;
0068             fy(:,1,:)   = +Py(:,2,:)/2+Py(:,1,:);
0069             fy(:,2,:)   = +Py(:,3,:)/2-Py(:,1,:);
0070             fy(:,end,:) = -Py(:,end,:)-Py(:,end-1,:)/2;
0071             fy(:,end-1,:) = Py(:,end,:)-Py(:,end-2,:)/2;
0072         end
0073         if nbdims>=3
0074             fz = (Pz(:,:,[2:end end])-Pz(:,:,[1 1:end-1]))/2;            
0075             fz(:,:,1)   = +Pz(:,:,2)/2+Pz(:,:,1);   % boundary
0076             fz(:,:,2)   = +Pz(:,:,3)/2-Pz(:,:,1);
0077             fz(:,:,end) = -Pz(:,:,end)-Pz(:,:,end-1)/2;
0078             fz(:,:,end-1) = Pz(:,:,end)-Pz(:,:,end-2)/2;
0079         end
0080     end 
0081 else
0082     if order==1
0083         fx = Px-Px([end 1:end-1],:,:);
0084         if nbdims>=2
0085             fy = Py-Py(:,[end 1:end-1],:);
0086         end
0087         if nbdims>=3
0088             fz = Pz-Pz(:,:,[end 1:end-1]);
0089         end
0090     else
0091         fx = (Px([2:end 1],:,:)-Px([end 1:end-1],:,:))/2;
0092         if nbdims>=2
0093             fy = (Py(:,[2:end 1],:)-Py(:,[end 1:end-1],:))/2;
0094         end
0095         if nbdims>=3
0096             fz = (Pz(:,:,[2:end 1])-Pz(:,:,[end 1:end-1]))/2;
0097         end
0098     end
0099 end
0100 
0101 % gather result
0102 if nbdims==3
0103     fd = fx+fy+fz;
0104 elseif nbdims==2
0105     fd = fx+fy;
0106 else
0107     fd = fx;
0108 end