function phi = basis_brick20 ( n, p )

%*****************************************************************************80
%
%% basis_brick20(): BRICK20 basis functions at natural coordinates.
%
%  Discussion:
%
%        8----19---7
%       /|        /|
%     20 |      18 |        t   s
%     /  16     /  15       |  /
%    5----17---6   |        | /
%    |   |     |   |        |/
%    |   4--11-|---3        0---------r
%   13  /     14  /
%    | 12      | 10
%    |/        |/
%    1----9----2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the number of evaluation points.
%
%    real P(3,N), natural coordinates of evaluation points.
%
%  Output:
%
%    real PHI(20,N), the basis function values.
%
  phi(1,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( - p(1,1:n) - p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(2,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( + p(1,1:n) - p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(3,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( + p(1,1:n) + p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(4,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) ...
    .* ( - p(1,1:n) + p(2,1:n) - p(3,1:n) - 2.0 ) / 8.0;
  phi(5,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( - p(1,1:n) - p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(6,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( + p(1,1:n) - p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(7,1:n) = ...
    ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( + p(1,1:n) + p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;
  phi(8,1:n) = ...
    ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) ...
    .* ( - p(1,1:n) + p(2,1:n) + p(3,1:n) - 2.0 ) / 8.0;

  phi(9,1:n) =  ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(10,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(11,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 + p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(12,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(13,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(14,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(15,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(16,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 + p(3,1:n) ) .* ( 1.0 - p(3,1:n) ) / 4.0;
  phi(17,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(18,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(19,1:n) = ( 1.0 + p(1,1:n) ) .* ( 1.0 - p(1,1:n) ) ...
              .* ( 1.0 + p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;
  phi(20,1:n) = ( 1.0 - p(1,1:n) ) .* ( 1.0 + p(2,1:n) ) ...
              .* ( 1.0 - p(2,1:n) ) .* ( 1.0 + p(3,1:n) ) / 4.0;

  return
end
