Re: divergence function

[Javier Enciso]

Jaroslav Hajek wrote:
> On Tue, Sep 29, 2009 at 11:16 AM, Javier Enciso <address@hidden> wrote:
> > Hi All,
> >
> > I' not sure whether this email already reach you guys, I wasn't member of
> > the mailing list, anyway I just resend it.
> >
> > I'm glad to contribute with the 'divergence' function. This function
> > is used to compute divergence of vector field and is part of Matlab's
> > core functions.
> >
> > In doubt, please feel free to contact me.
> >
> > Regards,
> > Javier
>
> Hello Javier,
>
> I have a couple of suggestions:
> 1. Please don't use tabs (or 8 spaces) to indent the code; use two
> spaces. This is done everywhere in Octave sources.
>
> 2. blocks like
> >                        x = varargin {1};
> >                        y = varargin {2};
> >                        z = varargin {3};
> >                        u = varargin {4};
> >                        v = varargin {5};
> >                        w = varargin {6};
> can be abbreviated to
>
> [x,y,z,u,v,w] = deal (varargin{1:6});
>
> this improves readability. Consider doing it.
>
> 3. in check_dim, the for loops can be avoided:
>
> function [] = check_dim (varargin)
>   if (!size_equal (varargin{:})
>     error ("Size of matrices must be equal");
>   endif
>
>   switch (nargin)
>   case {2, 4}
>     if (any (cellfun (@ndims, varargin) != 2))
>       error ("Matrices must have 2 dimensions");
>     endif
>   case {3, 6}                   
>     if (any (cellfun (@ndims, varargin) != 3))
>       error ("Matrices must have 3 dimensions");
>     endif
>   endswitch
> endfunction
>
>
> no big concern about performance in this case, since nargin is small,
> but it also seems to improve readability, at least to me...

Hi Jaroslav,

Thanks for your suggestions, it really helps to improve the readability 
of the code. ;-)

I already made the corrections and I hope the code is ready to be 
shipped. If you (or someone else) have any other concern, please let me 
know.

Regards,
Javier

## Copyright (C) 2009   Javier Enciso   <address@hidden>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see 
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function file} {} divergence (@var{X}, @var{Y}, @var{Z}, 
@var{U}, @var{V}, @var{W})
## @deftypefnx {Function file} {} divergence (@var{U}, @var{V}, @var{W})
## @deftypefnx {Function file} {} divergence (@var{X}, @var{Y}, @var{U}, 
@var{V})
## @deftypefnx {Function file} {} divergence (@var{U}, @var{V})
## Compute divergence of vector field. 
## 
## @code{divergence (X, Y, Z, U, V, W)} computes the divergence of a 
3-dimensional vector field 
## @var{U}, @var{V}, @var{W} at the points @var{X}, @var{Y}, @var{Z}. 
Coordinate points 
## must preserve the given order and define a 3-dimensional grid.
## 
## If @var{X}, @var{Y}, and @var{Z} are omitted, @code{divergence (U, V, W)} 
assumes 
## that the coordinates are given by the expression @code{[X, Y, Z] = meshgrid 
(1:n, 1:m, 1:p)}
## where @code{[m, n, p] = size (U)}.
## 
## @code{divergence (X, Y, U, V)} computes the divergence of a 2-dimensional 
vector field @var{U}, 
## @var{V} at the points @var{X}, @var{Y}. Coordinate points must preserve the 
given order 
## and define a 2-dimensional grid.
## 
## If @var{X} and @var{Y} are omitted, @code{divergence (U, V)} assumes that 
the 
## coordinate points @var{X} and @var{Y} are given by the expression @code{[X, 
Y] = 
## meshgrid (1:n, 1:m)} where @code{[m, n] = size (U)}.
## @end deftypefn

function [div] = divergence (varargin)
  switch (nargin)
    case 2
      [u, v] = deal (varargin{1:2});
      check_dim (u, v);
      [du, dump] = gradient (u);
      [dump, dv] = gradient (v);
      div = du + dv;
    case 3
      [u, v, w] = deal (varargin{1:3});
      check_dim (u, v, w);
      [du, dump, dump] = gradient (u);
      [dump, dv, dump] = gradient (v);
      [dump, dump, dw] = gradient (w);
      div = du + dv + dw;
    case 4
      [x, y, u, v] = deal (varargin{1:4});
      check_dim (x, y, u, v);
      [du, dump] = gradient (u);
      [dump, dv] = gradient (v);
      [dx, dump] = gradient (x);
      [dump, dy] = gradient (y);
      div = du./dx + dv./dy;
    case 6
      [x, y, z, u, v, w] = deal (varargin{1:6});
      check_dim (x, y, z, u, v, w);
      [du, dump, dump] = gradient (u);
      [dump, dv, dump] = gradient (v);
      [dump, dump, dw] = gradient (w);
      [dx, dump, dump] = gradient (x);
      [dump, dy, dump] = gradient (y);
      [dump, dump, dz] = gradient (z);
      div = du./dx + dv./dy + dw./dz;
    otherwise
      print_usage ();
  endswitch
endfunction

function [] = check_dim (varargin)
  if (!size_equal (varargin{:}))
    error ("Size of matrices must be equal");
  endif

  switch (nargin)
  case {2, 4}
    if (any (cellfun (@ndims, varargin) != 2))
      error ("Matrices must have 2 dimensions");
    endif
  case {3, 6}
    if (any (cellfun (@ndims, varargin) != 3))
      error ("Matrices must have 3 dimensions");
    endif
  endswitch
endfunction

%!test
%! a = [1, 3; 1, 2];
%! b = [3, 3; 1, 4];
%! expect = [0, 3; -1, 2];
%! get = divergence (a, b);
%! if (any(size (expect) != size (get)))
%!   error ("wrong size: expected %d,%d but got %d,%d", size (expect), size 
(get));
%! elseif (any (any (expect!=get)))
%!   error ("didn't get what was expected.");
%! endif

%!test
%! a(:,:,1) = [1, 3; 1, 2];
%! a(:,:,2) = [2, 1; 2, 2];
%! b(:,:,1) = [5, 5; 0, 4];
%! b(:,:,2) = [6, 2; 4, 2];
%! c(:,:,1) = [0, 1; 2, 3];
%! c(:,:,2) = [4, 5; 5, 4];
%! expect(:,:,1) = [1, 5; -1, 1];
%! expect(:,:,2) = [1, 3; 1, 1];
%! get = divergence (a, b, c);
%! if (any(size (expect) != size (get)))
%!   error ("wrong size: expected %d,%d but got %d,%d", size (expect), size 
(get));
%! elseif (any (any (expect!=get)))
%!   error ("didn't get what was expected.");
%! endif

%!test
%! a = rand (3, 3, 3);
%! b = rand (3, 3, 3);
%! c = rand (3, 3, 4);
%!error divergence (a, b, c);

%!test
%! a = rand (2, 3);
%! b = rand (2, 3);
%! c = rand (2, 3);
%!error divergence (a, b, c);