function elliptic_wavefront_domain ( nc )

%*****************************************************************************80
%
%% ELLIPTIC_WAVEFRONT_DOMAIN defines the domain for the wavefront problem.
%
%  Modified:
%
%    05 November 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NC, the grid parameter; this is the log base 2 of the
%    number of elements in each direction.  Thus NC = 3 generates an 8x8 grid.
%    2 <= NC.
%
%  Local Parameters:
%
%    Local, integer BOUND(4*N) or BOUND(4*2*N), the boundary vertex vector 
%    containing a list of nodes at which Dirichlet boundary conditions 
%    are imposed.
%
%    Local, integer GRID_TYPE, 1 for uniform, 2 for stretched.
%
%    Local, integer MBOUND(?,4), the boundary edge matrix.  The first entry
%    is an element index, and the second is 1, 2, 3 or 4, indicating that
%    the element touches the boundary on its bottom, right, top or left side.
%
%    Local, integer MV(E_NUM,9), the macroelement mapping.  MV(E,LN) contains 
%    the global node number GN of local node LN in element E.
%
%    Local, integer N, the number of elements along either the X or Y
%    direction.
%
%    Local, integer OUTBC, set this to 1.
%
%    Local, real X(N+1), the X grid coordinates.
%
%    Local, real XY((N+1)*(N+1),2), the X and Y coordinates of every node.
%
%    Local, real Y(N+1), the Y grid coordinates.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPTIC_WAVEFRONT_GRID:\n' );
  fprintf ( 1, '  Generate the grid information for the elliptic equation\n' );
  fprintf ( 1, '  modeling a wavefront on the unit square.\n' );

  if ( nargin < 1 )
    nc = 4;
  end

  n = 2 ^ nc;
  np = n / 2;
  nq = n / 4;
  grid_type = 1;
%
%  X and Y coordinates.
%
  x = linspace ( 0.0, 1.0, n + 1 );
  y = linspace ( 0.0, 1.0, n + 1 );
%
%  Compute biquadratic element coordinates.
%
  nvtx = ( n + 1 ) * ( n + 1 );
  [ X, Y ] = meshgrid ( x, y );
  xx = reshape ( X', nvtx, 1 );
  yy = reshape ( Y', nvtx, 1 );
  xy = [ xx(:), yy(:) ];

  kx = 1;
  ky = 1;
  mel = 0;
  for j = 1 : np
    for i = 1 : np
      mref = ( n + 1 ) * ( ky - 1 ) + kx;
      mel = mel + 1;
      nvv(1) = mref;
      nvv(2) = mref + 2;
      nvv(3) = mref + 2 * n + 4;
      nvv(4) = mref + 2 * n + 2;
      nvv(5) = mref + 1;
      nvv(6) = mref + n + 3; 
      nvv(7) = mref + 2 * n + 3; 
      nvv(8)=  mref + n + 1;
      nvv(9)=  mref + n + 2; 
      mv(mel,1:9) = nvv(1:9);
      kx = kx + 2;
    end
    ky = ky + 2; 
    kx = 1;
  end
%
%  Compute the boundary vertices and edges.
%
%  Four boundary edges.
%
  k1 = find ( xy(:,2) == 0.0 );
  e1 = [];
  for k = 1 : mel
    if any ( mv(k,5) == k1 )
      e1 = [ e1, k ]; 
    end
  end
  ef1 = ones ( size ( e1 ) );

  k2 = find ( xy(:,1) == 1.0 & xy(:,2) <= 1.0 & xy(:,2) > 0.0 );
  e2 = [];
  for k = 1 : mel
    if any ( mv(k,6) == k2 )
      e2 = [ e2, k ]; 
    end
  end
  ef2 = 2.0 * ones ( size ( e2 ) );

  k3 = find ( xy(:,2) == 1.0 & xy(:,1) < 1.0 & xy(:,1) > 0.0 );
  e3 = [];
  for k = 1 : mel
    if any ( mv(k,7) == k3 )
      e3 = [ e3, k ]; 
    end
  end
  ef3 = 3 * ones ( size(e3) );

  k4 = find ( xy(:,1) == 0.0 & xy(:,2) <= 1.0 & xy(:,2) > 0.0 );
  e4 = [];
  for k = 1 : mel
    if any ( mv(k,8) == k4 )
      e4 = [ e4, k ];
    end
  end
  ef4 = 4 * ones ( size ( e4 ) );

  bound = sort ( [ k1; k2; k3; k4 ] );

  mbound = [ e1', ef1'; ...
             e2', ef2'; ...
             e3', ef3'; ...
             e4', ef4' ];

  outbc = 1;
%
%  Move to the IFISS DATAFILES directory to place the datafile
%  in a standard location.
%
  gohome
  cd datafiles
%
%  Store the grid information.
%
  filename = 'elliptic_wavefront_grid.mat';

  save ( filename, 'mv', 'xy', 'bound', 'mbound', 'grid_type', 'outbc', ...
    'x', 'y' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The grid information has been stored in "datafiles/%s"\n', ...
    filename );

  return
end

