function fem2d_heat ( )

%*****************************************************************************80
%
%% MAIN is the main routine of FEM2D_HEAT.
%
%  Discussion:
%
%    This program solves
%
%      dUdT - Laplacian U(X,Y,T) = F(X,Y,T)
%
%    in a rectangular region in the plane.
%
%    Along the boundary of the region, Dirichlet conditions
%    are imposed:
%
%      U(X,Y,T) = G(X,Y,T)
%
%    At the initial time, the value of U is given at all points
%    in the region:
%
%      U(X,Y,T_INIT) = H(X,Y,T_INIT)
%
%    The code uses continuous piecewise quadratic basis functions on
%    triangles determined by a uniform grid of NX by NY points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Local parameters:
%
%    Local, real A(3*IB+1,NODE_NUM), the coefficient matrix.
%
%    Local, real EH1, the H1 seminorm error.
%
%    Local, real EL2, the L2 error.
%
%    Local, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
%    Local, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Local, integer ELEMENT_NUM, the number of elements.
%
%    Local, real F(NODE_NUM), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, integer NODE_BOUNDARY(NODE_NUM), is
%    0, if a node is an interior node;
%    1, if a node is a Dirichlet boundary node.
%
%    Local, integer NNODES, the number of nodes used to form one element.
%
%    Local, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Local, integer NX, the number of points in the X direction.
%
%    Local, integer NY, the number of points in the Y direction.
%
%    Local, integer QUAD_NUM, the number of quadrature points used for assembly.
%
%    Local, integer QUAD2_NUM, the number of quadrature points used for error
%    estimation.
%
%    Local, real TIME, the current time.
%
%    Local, real TIME_FINAL, the final time.
%
%    Local, real TIME_INIT, the initial time.
%
%    Local, real TIME_OLD, the time at the previous time step.
%
%    Local, integer TIME_STEP_NUM, the number of time steps to take.
%
%    Local, real TIME_STEP_SIZE, the size of the time steps.
%
%    Local, real U(NODE_NUM), the finite element coefficients defining 
%    the solution at the current time.
%
%    Local, real WQ(QUAD_NUM), quadrature weights.
%
%    Local, real XL, XR, YB, YT, the X coordinates of
%    the left and right sides of the rectangle, and the Y coordinates
%    of the bottom and top of the rectangle.
%
%    Local, real XQ(QUAD_NUM,ELEMENT_NUM), YQ(QUAD_NUM,ELEMENT_NUM), the
%    coordinates of the quadrature points in each element.
%
  nnodes = 6;
  quad_num = 3;
  quad2_num = 13;
  nx = 7;
  ny = 7;

  element_num = ( nx - 1 ) * ( ny - 1 ) * 2;
  node_num = ( 2 * nx - 1 ) * ( 2 * ny - 1 );

  xl = 0.0;
  xr = 1.0;
  yb = 0.0;
  yt = 1.0;

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_HEAT\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution of the time dependent heat equation on a unit\n' );
  fprintf ( 1, '  box in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Ut - Uxx - Uyy = F(x,y,t) in the box\n' );
  fprintf ( 1, '        U(x,y,t) = G(x,y,t) for (x,y) on the boundary.\n' );
  fprintf ( 1, '        U(x,y,t) = H(x,y,t) for t = T_INIT.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method is used, with piecewise\n' );
  fprintf ( 1, '  quadratic basis functions on 6 node triangular\n' );
  fprintf ( 1, '  elements.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The backward Euler formula is used for the time derivative.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The corner nodes of the triangles are generated by an\n' );
  fprintf ( 1, '  underlying grid whose dimensions are\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NX =                         %d\n', nx );
  fprintf ( 1, '  NY =                         %d\n', ny );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes =            %d\n', node_num );
  fprintf ( 1, '  Number of elements =         %d\n', element_num );
%
%  Set the XY coordinates of the nodes.
%
  node_xy = xy_set ( nx, ny, node_num, xl, xr, yb, yt );
%
%  Organize the nodes into a grid of 6-node triangles.
%
  element_node = grid_t6 ( nx, ny, nnodes, element_num );
%
%  Set the quadrature rule for assembly.
%
  [ wq, xq, yq ] = quad_a ( node_xy, element_node, element_num, node_num, ...
    nnodes );
%
%  Determine the areas of the elements.
%
  element_area = area_set ( node_num, node_xy, nnodes, element_num, ...
    element_node );
%
%  Determine which nodes have an associated finite element unknown.
%
  node_boundary = node_boundary_set ( nx, ny, node_num );
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( nnodes, element_num, element_node, node_num );

  fprintf ( 1, '  The matrix half bandwidth is %d\n', ib );
  fprintf ( 1, '  The matrix row size is       %d\n', 3 * ib + 1 );
%
%  Make an EPS picture of the nodes.
%
  if ( nx <= 10 & ny <= 10 )

    node_label = 1;
    node_eps_file_name = 'nodes.eps';

    nodes_plot ( node_eps_file_name, node_num, node_xy, node_label );

    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_HEAT:\n' );
    fprintf ( 1, '  Wrote an EPS node file\n' ); 
    fprintf ( 1, '    "%s".\n', node_eps_file_name );
    fprintf ( 1, '  containing a picture of the nodes.\n' );

  end
%
%  Write an ASCII node file.
%
  node_txt_file_name = 'nodes.txt';

  r8mat_write ( node_txt_file_name, 2, node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_HEAT:\n' );
  fprintf ( 1, '  Wrote an ASCII node file\n' );
  fprintf ( 1, '    "%s"\n', node_txt_file_name );
  fprintf ( 1, '  of the form\n' );
  fprintf ( 1, '    X(I), Y(I)\n' );
  fprintf ( 1, '  which can be used for plotting.\n' );
%
%  Make a picture of the elements.
%
  if ( nx <= 10 & ny <= 10 )

    triangulation_eps_file_name = 'elements.eps';
    node_show = 1;
    triangle_show = 2;

    triangulation_order6_plot ( triangulation_eps_file_name, node_num, ...
      node_xy, element_num, element_node, node_show, triangle_show );

    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_HEAT:\n' );
    fprintf ( 1, '  Wrote an EPS triangulation file\n' ); 
    fprintf ( 1, '    "%s".\n', triangulation_eps_file_name );
    fprintf ( 1, '  containing a picture of the triangulation.\n' );

  end
%
%  Write an ASCII triangle file.
%
  triangulation_txt_file_name = 'elements.txt';

  i4mat_write ( triangulation_txt_file_name, nnodes, element_num, element_node );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_HEAT:\n' );
  fprintf ( 1, '  Wrote an ASCII element file\n' );
  fprintf ( 1, '    "%s"\n', triangulation_txt_file_name );
  fprintf ( 1, '  of the form\n' );
  fprintf ( 1, '    Node(1) Node(2) Node(3) Node(4) Node(5) Node(6)\n' );
  fprintf ( 1, '  which can be used for plotting.\n' );
%
%  Set time stepping quantities.
%
  time_init = 0.0;
  time_final = 0.5;
  time_step_num = 10;
  time_step_size = ( time_final - time_init ) / time_step_num;
%
%  Initialize the names of the time and solution file.
%
  time_file_name = 'time.txt';
  u_file_name = 'u0000.txt';
%
%  Set the value of U at the initial time.
%
  time = time_init;
  [ u_exact, dudx_exact, dudy_exact ] = exact_u ( node_num, node_xy, time );
  u(1:node_num) = u_exact(1:node_num);

  time_unit = fopen ( time_file_name, 'wt');
  fprintf ( time_unit, '  %12f\n', time );
  r8mat_write ( u_file_name, 1, node_num, u );
%
%  Time looping.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '     Time        L2 Error       H1 Error\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  %12f  %12f  %12f\n', time, 0.0, 0.0 );

  for time_step = 1 : time_step_num

    time_old = time;
    u_old(1:node_num) = u(1:node_num);

    time = ( ( time_step_num - time_step ) * time_init    ...
           + (                 time_step ) * time_final ) ...
           / ( time_step_num             );
%
%  Assemble the finite element contributions to the coefficient matrix A 
%  and the right-hand side F.
%
    [ a, f ] = assemble ( node_num, node_xy, nnodes, element_num, ...
      element_node, quad_num, wq, xq, yq, element_area, ib, time );

    if ( 0 ) 
      dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, 5, 5, ...
        '  Initial 5 x 5 block of matrix A:' );

      r8vec_print_some ( node_num, f, 1, 10, '  Part of right hand side F:' );
    end
%
%  Modify the coefficient matrix and right hand side to account for the dU/dt
%  term, which we are treating using the backward Euler formula.
%
    [ a, f ] = adjust_backward_euler ( node_num, node_boundary, ib, ...
      time_step_size, u_old, a, f );

    if ( 0 )
      dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, 5, 5, ...
        '  A after DT adjustment:' );

      r8vec_print_some ( node_num, f, 1, 10, '  F after DT adjustment:' );
    end 
%
%  Modify the coefficient matrix and right hand side to account for
%  boundary conditions.
%
    [ a, f ] = adjust_boundary ( node_num, node_xy, node_boundary, ib, ...
      time, a, f );

    if ( 0 )
      dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, 5, 5, ...
        '  Matrix A after BC adjustment:'  );

      r8vec_print_some ( node_num, f, 1, 10, '  F after BC adjustment:' );
    end
%
%  Solve the linear system using a banded solver.
%
    [ a_lu, pivot, info ] = dgb_fa ( node_num, ib, ib, a );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'FEM2D_HEAT - Fatal error!\n' );
      fprintf ( 1, '  DGB_FA returned an error condition.\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  The linear system was not factored, and the\n' );
      fprintf ( 1, '  algorithm cannot proceed.\n' );
      error ( 'FEM2D_HEAT - Fatal error!' );
    end

    job = 0;
    u = dgb_sl ( node_num, ib, ib, a_lu, pivot, f, job );

    if ( 0 )
      r8vec_print_some ( node_num, u, 1, 10, '  Part of the solution vector:' );
    end
%
%  Calculate error using 13 point quadrature rule.
%
    [ el2, eh1 ] = errors ( element_area, element_node, node_xy, u, ...
      element_num, nnodes, quad2_num, node_num, time );

    if ( 0 )
      compare ( node_num, node_xy, time, u );
    end
%
%  Write an ASCII file that can be read into MATLAB.
%
    file_name_inc ( u_file_name );

    fprintf ( time_unit, '  %12f\n', time );

    r8mat_write ( u_file_name, 1, node_num, u );

  end

  fclose ( time_unit );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_HEAT:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ a, f ] = adjust_backward_euler ( node_num, node_boundary, ib, ...
  time_step_size, u_old, a, f )

%*****************************************************************************80
%
%% ADJUST_BACKWARD_EULER modifies the system for the backward Euler time term.
%
%  Discussion:
%
%    The term
%
%      dU/dT
%
%    is approximated by
%
%      ( U_NEW - U_OLD ) / DELTA T
%
%    which is handled by:
%
%      adding 1/DELTA T to each diagonal element of the system matrix, 
%      adding U_OLD / DELTA T to the right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer NODE_BOUNDARY(NODE_NUM), is
%    0, if a node is an interior node;
%    1, if a node is a Dirichlet boundary node.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real TIME_STEP_SIZE, the current time step.
%
%    Input, real U_OLD(NODE_NUM), the solution at the previous time.
%
%    Input, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM 
%    coefficient matrix, stored in a compressed format.  
%
%    Input, real F(NODE_NUM), the right hand side.
%
%    Output, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM
%    coefficient matrix, adjusted for boundary conditions.
%
%    Output, real F(NODE_NUM), the right hand side, adjusted for boundary 
%    conditions.
%

%
%  Consider each node.
%
  for node = 1 : node_num

    if ( node_boundary(node) == 0 )

      a(node-node+2*ib+1,node) = a(node-node+2*ib+1,node) ...
        + 1.0 / time_step_size;
    
      f(node) = f(node) + u_old(node) / time_step_size;

    end

  end

  return
end
function [ a, f ] = adjust_boundary ( node_num, node_xy, node_boundary, ib, ...
  time, a, f )

%*****************************************************************************80
%
%% ADJUST_BOUNDARY modifies the linear system for boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Input, integer NODE_BOUNDARY(NODE_NUM), is
%    0, if a node is an interior node;
%    1, if a node is a Dirichlet boundary node.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real TIME, the current time.
%
%    Input, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM
%    coefficient matrix, stored in a compressed format.
%
%    Input, real F(NODE_NUM), the right hand side.
%
%    Output, real A(3*IB+1,NODE_NUM), has been adjusted for boundary conditions.
%
%    Output, real F(NODE_NUM), has been adjusted for boundary conditions.
%

%
%  Get the exact solution at every node.
%
  u_exact = exact_u ( node_num, node_xy, time );

  for node = 1 : node_num

    if ( node_boundary(node) ~= 0 )

      jlo = max ( node - ib, 1 );
      jhi = min ( node + ib, node_num );

      for j = jlo : jhi
        a(node-j+2*ib+1,j) = 0.0;
      end

      a(node-node+2*ib+1,node) = 1.0;
    
      f(node) = u_exact(node);

    end

  end

  return
end
function element_area = area_set ( node_num, node_xy, nnodes, ...
  element_num, element_node )

%*****************************************************************************80
%
%% AREA_SET sets the area of each element.
%
%  Discussion:
%
%    The areas of the elements are needed in order to adjust
%    the integral estimates produced by the quadrature formulas.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer NNODES, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Output, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
  for element = 1 : element_num

    i1 = element_node(1,element);
    x1 = node_xy(1,i1);
    y1 = node_xy(2,i1);

    i2 = element_node(2,element);
    x2 = node_xy(1,i2);
    y2 = node_xy(2,i2);

    i3 = element_node(3,element);
    x3 = node_xy(1,i3);
    y3 = node_xy(2,i3);

    element_area(element) = 0.5E+00 * abs ...
      ( y1 * ( x2 - x3 ) ...
      + y2 * ( x3 - x1 ) ...
      + y3 * ( x1 - x2 ) );

  end

  return
end
function [ a, f ] = assemble ( node_num, node_xy, nnodes, element_num, ...
  element_node, quad_num, wq, xq, yq, element_area, ib, time )

%*****************************************************************************80
%
%% ASSEMBLE assembles the coefficient matrix A and right hand side F.
%
%  Discussion:
%
%    The matrix is known to be banded.  A special matrix storage format
%    is used to reduce the space required.  Details of this format are
%    discussed in the routine DGB_FA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer QUAD_NUM, the number of quadrature points used in assembly.
%
%    Input, real WQ(QUAD_NUM), quadrature weights.
%
%    Input, real XQ(QUAD_NUM,ELEMENT_NUM), YQ(QUAD_NUM,ELEMENT_NUM), the
%    coordinates of the quadrature points in each element.
%
%    Input, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, real TIME, the current time.
%
%    Output, real A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM
%    coefficient matrix, stored in a compressed format.
%
%    Output, real F(NODE_NUM), the right hand side.
%
%  Local parameters:
%
%    Local, real BB, BX, BY, the value of some basis function
%    and its first derivatives at a quadrature point.
%
%    Local, real BBB, BBX, BBY, the value of another basis
%    function and its first derivatives at a quadrature point.
%

%
%  Initialize the arrays to zero.
%
  f(1:node_num) = 0.0;
  a(1:3*ib+1,1:node_num) = 0.0;
%
%  The actual values of A and F are determined by summing up
%  contributions from all the elements.
%
  for element = 1 : element_num
%
%  Consider a quadrature point QUAD, with coordinates (X,Y).
%
    for quad = 1 : quad_num

      x = xq(quad,element);
      y = yq(quad,element);
      w = element_area(element) * wq(quad);
%
%  Consider one of the basis functions, which will play the
%  role of test function in the integral.
%
%  We generate an integral for every node associated with an unknown.
%  But if a node is associated with a boundary condition, we do nothing.
%
      for test = 1 : nnodes

        i = element_node(test,element);

        [ bi, dbidx, dbidy ] = qbf ( x, y, element, test, node_xy, ...
          element_node, element_num, nnodes, node_num );

        f(i) = f(i) + w * rhs ( x, y, time ) * bi;
%
%  Consider a basis function used to form the value of the solution function.
%
%  If this basis function is associated with a boundary condition,
%  then subtract the term from the right hand side.
%
%  Otherwise, this term is included in the system matrix.
%
%  Logically, this term goes in entry A(I,J).  Because of the
%  band matrix storage, entry (I,J) is actually stored in
%  A(I-J+2*NHBA+1,J).
%
        for basis = 1 : nnodes

          j = element_node(basis,element);

          [ bj, dbjdx, dbjdy ] = qbf ( x, y, element, basis, node_xy, ...
            element_node, element_num, nnodes, node_num );

          aij = dbidx * dbjdx + dbidy * dbjdy;

          a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) + w * aij;

        end

      end

    end

  end

  return
end
function nhba = bandwidth ( nnodes, element_num, element_node, node_num )

%*****************************************************************************80
%
%% BANDWIDTH determines the bandwidth of the coefficient matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NNODES, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, integer NHBA, the half bandwidth of the matrix.
%
  nhba = 0;
  for element = 1 : element_num
    for iln = 1 : nnodes
      i = element_node(iln,element);
      for jln = 1 : nnodes
        j = element_node(jln,element);
        nhba = max ( nhba, j - i );
      end
    end
  end

  return
end
function compare ( node_num, node_xy, time, u )

%*****************************************************************************80
%
%% COMPARE compares the exact and computed solution at the nodes.
%
%  Discussion:
%
%    This is a rough comparison, done only at the nodes.  Such a pointwise
%    comparison is easy, because the value of the finite element
%    solution is exactly the value of the finite element coefficient
%    associated with that node.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, real TIME, the current time.
%
%    Input, real U(NODE_NUM), the solution vector of the finite
%    element system.
%
  [ u_exact, dudx_exact, dudy_exact ] = exact_u ( node_num, node_xy, time );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMPARE:\n' );
  fprintf ( 1, '  Compare computed and exact solutions at the nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  TIME = %f\n', time );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         X           Y          U           U\n' );
  fprintf ( 1, '                              exact      computed\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num

    x = node_xy(1,node);
    y = node_xy(2,node);

    fprintf ( 1, '%12f  %12f  %12f  %12f\n', x, y, u_exact(node), u(node) );

  end

  return
end
function [ alu, pivot, info ] = dgb_fa ( n, ml, mu, a )

%*****************************************************************************80
%
%% DGB_FA performs a LINPACK-style PLU factorization of a DGB matrix.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "DGB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%    The following program segment will set up the input.
%
%      m = ml + mu + 1
%      do j = 1, n
%        i1 = max ( 1, j-mu )
%        i2 = min ( n, j+ml )
%        do i = i1, i2
%          k = i - j + m
%          a(k,j) = afull(i,j)
%        end do
%      end do
%
%    This uses rows ML+1 through 2*ML+MU+1 of the array A.
%    In addition, the first ML rows in the array are used for
%    elements generated during the triangularization.
%    The total number of rows needed in A is 2*ML+MU+1.
%    The ML+MU by ML+MU upper left triangle and the
%    ML by ML lower right triangle are not referenced.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2003
%
%  Author:
%
%    Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N-1.
%
%    Input, real A(2*ML+MU+1,N), the matrix in band storage.  The
%    columns of the matrix are stored in the columns of the array,
%    and the diagonals of the matrix are stored in rows ML+1 through
%    2*ML+MU+1.
%
%    Output, real ALU(2*ML+MU+1,N), the LU factors in band storage.  
%    The L and U matrices are stored in a single array.
%
%    Output, integer PIVOT(N), the pivot vector.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
  alu = a(1:2*ml+mu+1,1:n);

  m = ml + mu + 1;
  info = 0;
%
%  Zero out the initial fill-in columns.
%
  j0 = mu + 2;
  j1 = min ( n, m ) - 1;

  for jz = j0 : j1
    i0 = m + 1 - jz;
    alu(i0:ml,jz) = 0.0E+00;
  end

  jz = j1;
  ju = 0;

  for k = 1 : n - 1
%
%  Zero out the next fill-in column.
%
    jz = jz + 1;
    if ( jz <= n )
      alu(1:ml,jz) = 0.0E+00;
    end
%
%  Find L = pivot index.
%
    lm = min ( ml, n-k );

    l = m;

    for j = m+1 : m+lm
      if ( abs ( alu(l,k) ) < abs ( alu(j,k) ) )
        l = j;
      end
    end

    pivot(k) = l + k - m;
%
%  Zero pivot implies this column already triangularized.
%
    if ( alu(l,k) == 0.0E+00 )
      info = k;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'DGB_FA - Fatal error!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return;
    end
%
%  Interchange if necessary.
%
    t = alu(l,k);
    alu(l,k) = alu(m,k);
    alu(m,k) = t;
%
%  Compute multipliers.
%
    alu(m+1:m+lm,k) = - alu(m+1:m+lm,k) / alu(m,k);
%
%  Row elimination with column indexing.
%
    ju = max ( ju, mu+pivot(k) );
    ju = min ( ju, n );
    mm = m;

    for j = k+1 : ju

      l = l - 1;
      mm = mm - 1;

      if ( l ~= mm )
        t = alu(l,j);
        alu(l,j) = alu(mm,j);
        alu(mm,j) = t;
      end

      alu(mm+1:mm+lm,j) = alu(mm+1:mm+lm,j) + alu(mm,j) * alu(m+1:m+lm,k);

    end

  end

  pivot(n) = n;

  if ( alu(m,n) == 0.0E+00 )
    info = n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DGB_FA - Fatal error!\n' );
    fprintf ( 1, '  Zero pivot on step %d\n', info );
  end

  return
end
function dgb_print_some ( m, n, ml, mu, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% DGB_PRINT_SOME prints some of a DGB matrix.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "DGB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%    Only entries in rows ILO to IHI, columns JLO to JHI are considered.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than min(M,N)-1..
%
%    Input, real A(2*ML+MU+1,N), the M by N band matrix, stored in LINPACK
%    or LAPACK general band storage mode.
%
%    Input, integer ILO, JLO, IHI, JHI, designate the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  incx = 5;
%
%  Print the columns of the matrix, in strips of 5.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2lo = max ( i2lo, j2lo - mu );

    i2hi = min ( ihi, m );
    i2hi = min ( i2hi, j2hi + ml );

    for i = i2lo : i2hi
%
%  Print out (up to) 5 entries in row I, that lie in the current strip.
%
      fprintf ( 1, '%6i  ', i );
      for j = j2lo: j2hi

        if ( ml < i-j | mu < j-i )
          fprintf ( 1, '              ' );
        else
          fprintf ( 1, '%14f', a(i-j+ml+mu+1,j) );
        end

      end

      fprintf ( 1, '\n' );

    end

  end

  fprintf ( 1, '\n' );

  return
end
function x = dgb_sl ( n, ml, mu, a, pivot, b, job )

%*****************************************************************************80
%
%% DGB_SL solves a system factored by DGB_FA.
%
%  Discussion:
%
%    An M by N banded matrix A with lower bandwidth ML and upper bandwidth MU
%    is assumed to be entirely zero, except for the main diagonal, and
%    entries in the ML nearest subdiagonals, and MU nearest superdiagonals.
%
%    LINPACK and LAPACK "DGB" storage for such a matrix generally includes
%    room for ML extra superdiagonals, which may be required to store
%    nonzero entries generated during Gaussian elimination.
%
%    The original M by N matrix is "collapsed" downward, so that diagonals
%    become rows of the storage array, while columns are preserved.  The
%    collapsed array is logically 2*ML+MU+1 by N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2003
%
%  Author:
%
%    Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N-1.
%
%    Input, real A(2*ML+MU+1,N), the LU factors from DGB_FA.
%
%    Input, integer PIVOT(N), the pivot vector from DGB_FA.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, integer JOB.
%    0, solve A * x = b.
%    nonzero, solve A' * x = b.
%
%    Output, real X(N), the solution.
%
  x(1:n) = b(1:n);

  m = mu + ml + 1;
%
%  Solve A * x = b.
%
  if ( job == 0 )
%
%  Solve L * Y = B.
%
    if ( 1 <= ml )

      for k = 1 : n-1

        lm = min ( ml, n-k );
        l = pivot(k);

        if ( l ~= k )
          t = x(l);
          x(l) = x(k);
          x(k) = t;
        end
        
        for i = 1 : lm
          x(k+i) = x(k+i) + x(k) * a(m+i,k);
        end 
      
      end
    end
%
%  Solve U * X = Y.
%
    for k = n : -1 : 1

      x(k) = x(k) / a(m,k);
      lm = min ( k, m ) - 1;
      la = m - lm;
      lb = k - lm;

      for i = 0: lm-1
        x(lb+i) = x(lb+i) - x(k) * a(la+i,k);
      end

    end
%
%  Solve A' * X = B.
%
  else
%
%  Solve U' * Y = B.
%
    for k = 1 : n
      lm = min ( k, m ) - 1;
      la = m - lm;
      lb = k - lm;
      for i = 0 : lm - 1
        x(k) = x(k) - a(la+i,k) * x(lb+i);
      end
      x(k) = x(k) / a(m,k);
    end
%
%  Solve L' * X = Y.
%
    if ( 1 <= ml )

      for k = n-1: -1 : 1

        lm = min ( ml, n-k );
        for i = 1 : lm
          x(k) = x(k) + a(m+i,k) * x(k+i);
        end
        l = pivot(k);

        if ( l ~= k )
          t = x(l);
          x(l) = x(k);
          x(k) = t;
        end

      end

    end

  end

  return
end
function [ el2, eh1 ] = errors ( element_area, element_node, node_xy, u, ...
  element_num, nnodes, quad2_num, node_num, time )

%*****************************************************************************80
%
%% ERRORS calculates the L2 and H1-seminorm errors.
%
%  Discussion:
%
%    This routine uses a 13 point quadrature rule in each element,
%    in order to estimate the values of
%
%      EL2(t) = Sqrt ( Integral ( U(x,y,t) - Uh(x,y,t) )**2 dx dy )
%
%      EH1(t) = Sqrt ( Integral ( Ux(x,y,t) - Uhx(x,y,t) )**2 + 
%                               ( Uy(x,y,t) - Uhy(x,y,t) )**2 dx dy )
%
%    Here U is the exact solution, and Ux and Uy its spatial derivatives,
%    as evaluated by a user-supplied routine.
%
%    Uh, Uhx and Uhy are the computed solution and its spatial derivatives,
%    as specified by the computed finite element solution.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, real U(NODE_NUM), the coefficients of the solution.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Input, integer QUAD2_NUM, the number of points in the quadrature rule.
%    This is actually fixed at 13.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real TIME, the current time.
%
%    Output, real EL2, the L2 error.
%
%    Output, real EH1, the H1 seminorm error.
%
%  Local Parameters:
%
%    Local, real AR, the weight for a given quadrature point
%    in a given element.
%
%    Local, real BB, BX, BY, a basis function and its first
%    derivatives evaluated at a particular quadrature point.
%
%    Local, real EH1, the H1 seminorm error.
%
%    Local, real EL2, the L2 error.
%
%    Local, real UEX, UEXX, UEXY, the exact solution and its first
%    derivatives evaluated at a particular quadrature point.
%
%    Local, real UH, UHX, UHY, the computed solution and its first
%    derivatives evaluated at a particular quadrature point.
%
%    Local, real WQE(NQE), stores the quadrature weights.
%
%    Local, real X, Y, the coordinates of a particular
%    quadrature point.
%
%    Local, real XQE(NQE), YQE(NQE), stores the location
%    of quadrature points in a given element.
%
  el2 = 0.0;
  eh1 = 0.0;

  for element = 1 : element_num

    [ wqe, xqe, yqe ] = quad_e ( node_xy, element_node, element, ...
      element_num, nnodes, node_num, quad2_num );

    for quad = 1 : quad2_num

      ar = element_area(element) * wqe(quad);
      x = xqe(quad);
      y = yqe(quad);

      uh = 0.0;
      dudxh = 0.0;
      dudyh = 0.0;

      for in1 = 1 : nnodes

        i = element_node(in1,element);

        [ bi, dbidx, dbidy ] = qbf ( x, y, element, in1, node_xy, ...
          element_node, element_num, nnodes, node_num );

        uh    = uh    + bi    * u(i);
        dudxh = dudxh + dbidx * u(i);
        dudyh = dudyh + dbidy * u(i);

      end

      xy(1,1) = x;
      xy(2,1) = y;

      [ u_exact, dudx_exact, dudy_exact ] = exact_u ( 1, xy, time );

      el2 = el2 + ar * ( uh - u_exact(1) )^2 ;

      eh1 = eh1 + ar * ( ( dudxh - dudx_exact(1) )^2 ...
                       + ( dudyh - dudy_exact(1) )^2 );

    end

  end

  el2 = sqrt ( el2 );
  eh1 = sqrt ( eh1 );

  fprintf ( 1, '  %12f  %12f  %12f\n', time, el2, eh1 );

  return
end
function [ u, dudx, dudy ] = exact_u ( node_num, node_xy, time )

%*****************************************************************************80
%
%% EXACT_U calculates the exact solution.
%
%  Discussion:
%
%    It is assumed that the user knows the exact solution and its
%    derivatives.  This, of course, is NOT true for a real computation.
%    But for this code, we are interested in studying the convergence
%    behavior of the approximations, and so we really need to assume
%    we know the correct solution.
%
%    As a convenience, this single routine is used for several purposes:
%
%    * it supplies the initial value function H(X,Y,T);
%    * it supplies the boundary value function G(X,Y,T);
%    * it is used by the COMPARE routine to make a node-wise comparison
%      of the exact and approximate solutions.
%    * it is used by the ERRORS routine to estimate the integrals of
%      the L2 and H1 errors of approximation.
%
%    DUDX and DUDY are only needed for the ERRORS calculation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes at which
%    a value is desired.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of 
%    the points where a value is desired.
%
%    Input, real TIME, the current time.
%
%    Output, real U(NODE_NUM), the exact solution.
%
%    Output, real DUDX(NODE_NUM), DUDY(NODE_NUM),
%    the X and Y derivatives of the exact solution.
%
  u(1:node_num) =         sin ( pi * node_xy(1,1:node_num)' ) ...
                       .* sin ( pi * node_xy(2,1:node_num)' ) ...
                        * exp ( -time );
 
  dudx(1:node_num) = pi * cos ( pi * node_xy(1,1:node_num)' ) ...
                       .* sin ( pi * node_xy(2,1:node_num)' ) ...
                        * exp ( -time );

  dudy(1:node_num) = pi * sin ( pi * node_xy(1,1:node_num)' ) ...
                       .* cos ( pi * node_xy(2,1:node_num)' ) ...
                        * exp ( -time );

  return
end
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME(*), the string to be incremented.
%
%    Output, character FILE_NAME_NEW(*), the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;
      
      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';
      
      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end
function element_node = grid_t6 ( nx, ny, nnodes, element_num )

%*****************************************************************************80
%
%% GRID_T6 produces a grid of pairs of 6 node triangles.
%
%  Example:
%
%    Input:
%
%      NX = 4, NY = 3
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  3, 15,  2,  9,  8;
%        17, 15,  3, 16,  9, 10;
%         3,  5, 17,  4, 11, 10;
%        19, 17,  5, 18, 11, 12;
%         5,  7, 19,  6, 13, 12;
%        21, 19,  7, 20, 13, 14;
%        15, 17, 29, 16, 23, 22;
%        31, 29, 17, 30, 23, 24;
%        17, 19, 31, 18, 25, 24;
%        33, 31, 19, 32, 25, 26;
%        19, 21, 33, 20, 27, 26;
%        35, 33, 21, 34, 27, 28.
%
%  Diagram:
%
%   29-30-31-32-33-34-35
%    |\ 8  |\10  |\12  |
%    | \   | \   | \   |
%   22 23 24 25 26 27 28
%    |   \ |   \ |   \ |
%    |  7 \|  9 \| 11 \|
%   15-16-17-18-19-20-21
%    |\ 2  |\ 4  |\ 6  |
%    | \   | \   | \   |
%    8  9 10 11 12 13 14
%    |   \ |   \ |   \ |
%    |  1 \|  3 \|  5 \|
%    1--2--3--4--5--6--7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, controls the number of elements along the
%    X and Y directions.  The number of elements will be
%    2 * ( NX - 1 ) * ( NY - 1 ).
%
%    Input, integer NNODES, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Output, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the index of the I-th node of the J-th element.
%
  element = 0;

  for j = 1: ny - 1
    for i = 1 : nx - 1

      sw = ( j - 1 ) * 2 * ( 2 * nx - 1 ) + 2 * i - 1;
      w  = sw + 1;
      nw = sw + 2;

      s  = sw + 2 * nx - 1;
      c  = s  + 1;
      n  = s  + 2;

      se = s  + 2 * nx - 1;
      e  = se + 1;
      ne = se + 2;

      element = element + 1;
      element_node(1,element) = sw;
      element_node(2,element) = se;
      element_node(3,element) = nw;
      element_node(4,element) = s;
      element_node(5,element) = c;
      element_node(6,element) = w;

      element = element + 1;
      element_node(1,element) = ne;
      element_node(2,element) = nw;
      element_node(3,element) = se;
      element_node(4,element) = n;
      element_node(5,element) = c;
      element_node(6,element) = e;

    end
  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 ) 
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function [ indx, nu ] = indx_set ( nx, ny, node_num )

%*****************************************************************************80
%
%% INDX_SET assigns a boundary value index or unknown value index at each node.
%
%  Discussion:
%
%    Every finite element node will is assigned an index which
%    indicates the finite element basis function and its coefficient
%    which are associated with that node.
%
%  Example:
%
%    On a simple 5 by 5 grid, where the nodes are numbered starting 
%    at the lower left, and increasing in X first, we would have the 
%    following values of INDX:
%
%       21  22  23  24  25
%       16  17  18  19  20
%       11  12  13  14  15
%        6   7   8   9  10
%        1   2   3   4   5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of elements in the X and Y directions.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, integer INDX(NODE_NUM), the index of the unknown in the finite 
%    element linear system.
%
%    Output, integer NU, the number of unknowns.
%
  nu = 0;
  in = 0;

  for j = 1 : 2 * ny - 1

    for i = 1 : 2 * nx - 1

      in = in + 1;
      nu = nu + 1;
      indx(in) = nu;

    end
  end

  return
end
function node_boundary = node_boundary_set ( nx, ny, node_num )

%*****************************************************************************80
%
%% NODE_BOUNDARY_SET stores the boundary status of each node in BOUNDARY_NODE.
%
%  Discussion:
%
%    Every finite element node will is assigned an index which
%    indicates the finite element basis function and its coefficient
%    which are associated with that node.
%
%  Example:
%
%    On a simple 5 by 5 grid, where the nodes are numbered starting 
%    at the lower left, and increasing in X first, we would have the 
%    following values of NODE_BOUNDARY:
%
%       1  1  1  1  1
%       1  0  0  0  1
%       1  0  0  0  1
%       1  0  0  0  1
%       1  1  1  1  1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, the number of elements in the X and Y directions.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, integer NODE_BOUNDARY(NODE_NUM), is
%    0, if a node is an interior node;
%    1, if a node is a Dirichlet boundary node.
%
  node = 0;

  for j = 1 : 2 * ny - 1

    for i = 1 : 2 * nx - 1

      node = node + 1;

      if ( j == 1 | ...
           j == 2 * ny - 1 | ...
           i == 1 | ...
           i == 2 * nx - 1 )
        node_boundary(node) = 1;
      else
        node_boundary(node) = 0;
      end

    end
  end

  return
end
function nodes_plot ( file_name, node_num, node_xy, node_label )

%*****************************************************************************80
%
%% NODES_PLOT plots the nodes.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the file to create.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, logical NODE_LABEL, is TRUE if the nodes should be labeled.
%
  x_ps_max = 576;
  x_ps_max_clip = 594;
  x_ps_min = 36;
  x_ps_min_clip = 18;
  y_ps_max = 666;
  y_ps_max_clip = 684;
  y_ps_min = 126;
  y_ps_min_clip = 108;
%
%  We need to do some figuring here, so that we can determine
%  the range of the data, and hence the height and width
%  of the piece of paper.
%
  x_max = max ( node_xy(1,1:node_num) );
  x_min = min ( node_xy(1,1:node_num) );
  x_scale = x_max - x_min;

  x_max = x_max + 0.05 * x_scale;
  x_min = x_min - 0.05 * x_scale;
  x_scale = x_max - x_min;

  y_max = max ( node_xy(2,1:node_num) );
  y_min = min ( node_xy(2,1:node_num) );
  y_scale = y_max - y_min;

  y_max = y_max + 0.05 * y_scale;
  y_min = y_min - 0.05 * y_scale;
  y_scale = y_max - y_min;

  if ( x_scale < y_scale ) 

    delta = round ( ( x_ps_max - x_ps_min ) ...
      * ( y_scale - x_scale ) / ( 2.0 * y_scale ) );

    x_ps_max = x_ps_max - delta;
    x_ps_min = x_ps_min + delta;

    x_ps_max_clip = x_ps_max_clip - delta;
    x_ps_min_clip = x_ps_min_clip + delta;

    x_scale = y_scale;

  elseif ( y_scale < x_scale )

    delta = round ( ( y_ps_max - y_ps_min ) ...
      * ( x_scale - y_scale ) / ( 2.0 * x_scale ) );

    y_ps_max      = y_ps_max - delta;
    y_ps_min      = y_ps_min + delta;

    y_ps_max_clip = y_ps_max_clip - delta;
    y_ps_min_clip = y_ps_min_clip + delta;

    y_scale = x_scale;

  end

  file_unit = fopen ( file_name, 'wt' );

  if ( file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NODES_PLOT - Fatal error!\n');
    fprintf ( 1, '  Could not open the output EPS file.\n' );
    error ( 'NODES_PLOT - Fatal error!' );
  end

  fprintf ( file_unit, '%!PS-Adobe-3.0 EPSF-3.0\n');
  fprintf ( file_unit, '%%Creator: nodes_plot.m\n');
  fprintf ( file_unit, '%%Title: %s\n', file_name );
  fprintf ( file_unit, '%%Pages: 1\n');
  fprintf ( file_unit, '%%BoundingBox:  %d  %d  %d  %d\n', ... 
    x_ps_min, y_ps_min, x_ps_max, y_ps_max );
  fprintf ( file_unit, '%%Document-Fonts: Times-Roman\n');
  fprintf ( file_unit, '%%LanguageLevel: 1\n');
  fprintf ( file_unit, '%%EndComments\n');
  fprintf ( file_unit, '%%BeginProlog\n');
  fprintf ( file_unit, '/inch {72 mul} def\n');
  fprintf ( file_unit, '%%EndProlog\n');
  fprintf ( file_unit, '%%Page:      1     1\n');
  fprintf ( file_unit, 'save\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '% Set RGB line color to very light gray.\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, ' 0.9000 0.9000 0.9000 setrgbcolor\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '% Draw a gray border around the page.\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, 'newpath\n');
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, 'stroke\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '% Set RGB line color to black.\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, ' 0.0000 0.0000 0.0000 setrgbcolor\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '% Set the font and its size.\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '/Times-Roman findfont\n');
  fprintf ( file_unit, '0.50 inch scalefont\n');
  fprintf ( file_unit, 'setfont\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Print a title:\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  210  702  moveto\n');
  fprintf ( file_unit, '%  (Pointset)  show\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Define a clipping polygon\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%newpath\n');
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, 'clip newpath\n');
%
%  Draw the nodes.
%
  if ( node_num <= 200 )
    circle_size = 5;
  elseif ( node_num <= 500 )
    circle_size = 4;
  elseif ( node_num <= 1000 )
    circle_size = 3;
  elseif ( node_num <= 5000 )
    circle_size = 2;
  else
    circle_size = 1;
  end

  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Draw filled dots at each node:\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Set the RGB color to blue:\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, ' 0.000 0.150 0.750 setrgbcolor\n');

  for node = 1 : node_num

    x_ps = floor ( ...
      ( ( x_max - node_xy(1,node)         ) * x_ps_min   ...
      + (       + node_xy(1,node) - x_min ) * x_ps_max ) ...
      / ( x_max                   - x_min ) );

    y_ps = floor ( ...
      ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
      + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
      / ( y_max                   - y_min ) );

    fprintf ( file_unit, ...
      'newpath  %d  %d  %d  0 360 arc closepath fill\n', ...
      x_ps, y_ps, circle_size );

  end
%
%  Label the nodes.
%
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Label the nodes:\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%  Set the RGB color to darker blue:\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, ' 0.000 0.250 0.850 setrgbcolor\n');

  fprintf ( file_unit, ' 0.0000 0.0000 0.0000 setrgbcolor\n');
  fprintf ( file_unit, '/Times-Roman findfont\n');
  fprintf ( file_unit, '0.20 inch scalefont\n');
  fprintf ( file_unit, 'setfont\n');

  for node = 1 : node_num

    x_ps = floor ( ...
      ( ( x_max - node_xy(1,node)         ) * x_ps_min   ...
      + (       + node_xy(1,node) - x_min ) * x_ps_max ) ...
      / ( x_max                   - x_min ) );

    y_ps = floor ( ...
      ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
      + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
      / ( y_max                   - y_min ) );

    fprintf ( file_unit, '%d  %d  moveto (%d) show\n', x_ps, y_ps+5, node );

  end
 
  fprintf ( file_unit, '%\n' );
  fprintf ( file_unit, 'restore showpage\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '% End of page\n');
  fprintf ( file_unit, '%\n');
  fprintf ( file_unit, '%%Trailer\n');
  fprintf ( file_unit, '%%EOF\n');

  fclose ( file_unit );

  return
end
function [ b, dbdx, dbdy ] = qbf ( x, y, element, inode, node_xy, ...
  element_node, element_num, nnodes, node_num )

%*****************************************************************************80
%
%% QBF evaluates the quadratic basis functions.
%
%  Discussion:
%
%    This routine assumes that the "midpoint" nodes are, in fact,
%    exactly the average of the two extreme nodes.  This is NOT true
%    for a general quadratic triangular element.
%
%    Assuming this property of the midpoint nodes makes it easy to
%    determine the values of (R,S) in the reference element that
%    correspond to (X,Y) in the physical element.
%
%    Once we know the (R,S) coordinates, it's easy to evaluate the
%    basis functions and derivatives.
%
%  The physical element T6:
%
%    In this picture, we don't mean to suggest that the bottom of
%    the physical triangle is horizontal.  However, we do assume that
%    each of the sides is a straight line, and that the intermediate
%    points are exactly halfway on each side.
%
%    |
%    |
%    |        3
%    |       / \
%    |      /   \
%    Y     6     5
%    |    /       \
%    |   /         \
%    |  1-----4-----2
%    |
%    +--------X-------->
%
%  Reference element T6:
%
%    In this picture of the reference element, we really do assume
%    that one side is vertical, one horizontal, of length 1.
%
%    |
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-------->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the (global) coordinates of the point
%    at which the basis function is to be evaluated.
%
%    Input, integer ELEMENT, the index of the element which contains the point.
%
%    Input, integer INODE, the local index (between 1 and 6) that
%    specifies which basis function is to be evaluated.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, real B, DBDX, DBDY, the value of the basis function
%    and its X and Y derivatives at (X,Y).
%
  xn(1:6) = node_xy(1,element_node(1:6,element));
  yn(1:6) = node_xy(2,element_node(1:6,element));
%
%  Determine the (R,S) coordinates corresponding to (X,Y).
%
%  What is happening here is that we are solving the linear system:
%
%    ( X2-X1  X3-X1 ) * ( R ) = ( X - X1 )
%    ( Y2-Y1  Y3-Y1 )   ( S )   ( Y - Y1 )
%
%  by computing the inverse of the coefficient matrix and multiplying
%  it by the right hand side to get R and S.
%
%  The values of dRdX, dRdY, dSdX and dSdY are easily from the formulas
%  for R and S.
%
  det =   ( xn(2) - xn(1) ) * ( yn(3) - yn(1) ) ...
        - ( xn(3) - xn(1) ) * ( yn(2) - yn(1) );

  r = ( ( yn(3) - yn(1) ) * ( x     - xn(1) ) ...
      + ( xn(1) - xn(3) ) * ( y     - yn(1) ) ) / det;

  drdx = ( yn(3) - yn(1) ) / det;
  drdy = ( xn(1) - xn(3) ) / det;

  s = ( ( yn(1) - yn(2) ) * ( x     - xn(1) ) ...
      + ( xn(2) - xn(1) ) * ( y     - yn(1) ) ) / det;

  dsdx = ( yn(1) - yn(2) ) / det;
  dsdy = ( xn(2) - xn(1) ) / det;
%
%  The basis functions can now be evaluated in terms of the
%  reference coordinates R and S.  It's also easy to determine
%  the values of the derivatives with respect to R and S.
%
  if ( inode == 1 )

    b    =   2.0E+00 *     ( 1.0E+00 - r - s ) * ( 0.5E+00 - r - s );
    dbdr = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s;
    dbds = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s;

  elseif ( inode == 2 )

    b    =   2.0E+00 * r * ( r - 0.5E+00 );
    dbdr = - 1.0E+00 + 4.0E+00 * r;
    dbds =   0.0E+00;

  elseif ( inode == 3 )

    b    =   2.0E+00 * s * ( s - 0.5E+00 );
    dbdr =   0.0E+00;
    dbds = - 1.0E+00               + 4.0E+00 * s;

  elseif ( inode == 4 )

    b    =   4.0E+00 * r * ( 1.0E+00 - r - s );
    dbdr =   4.0E+00 - 8.0E+00 * r - 4.0E+00 * s;
    dbds =           - 4.0E+00 * r;

  elseif ( inode == 5 )

    b    =   4.0E+00 * r * s;
    dbdr =                           4.0E+00 * s;
    dbds =             4.0E+00 * r;

  elseif ( inode == 6 )

    b    =   4.0E+00 * s * ( 1.0E+00 - r - s );
    dbdr =                         - 4.0E+00 * s;
    dbds =   4.0E+00 - 4.0E+00 * r - 8.0E+00 * s;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'QBF - Fatal error!\n' );
    fprintf ( 1, '  Request for local basis function INODE = %d\n', inode );
    error ( 'QBF - Fatal error!' );

  end
%
%  We need to convert the derivative information from (R(X,Y),S(X,Y))
%  to (X,Y) using the chain rule.
%
  dbdx = dbdr * drdx + dbds * dsdx;
  dbdy = dbdr * drdy + dbds * dsdy;

  return
end
function [ wq, xq, yq ] = quad_a ( node_xy, element_node, ...
  element_num, node_num, nnodes )

%*****************************************************************************80
%
%% QUAD_A sets the quadrature rule for assembly.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Output, real WQ(3), quadrature weights.
%
%    Output, real XQ(3,ELEMENT_NUM), YQ(3,ELEMENT_NUM), the X and Y
%    coordinates of the quadrature points in each element.
%
  wq(1) = 1.0 / 3.0;
  wq(2) = wq(1);
  wq(3) = wq(1);

  for element = 1 : element_num

    ip1 = element_node(1,element);
    ip2 = element_node(2,element);
    ip3 = element_node(3,element);

    x1 = node_xy(1,ip1);
    x2 = node_xy(1,ip2);
    x3 = node_xy(1,ip3);

    y1 = node_xy(2,ip1);
    y2 = node_xy(2,ip2);
    y3 = node_xy(2,ip3);

    xq(1,element) = 0.5E+00 * ( x1 + x2 );
    xq(2,element) = 0.5E+00 * ( x2 + x3 );
    xq(3,element) = 0.5E+00 * ( x1 + x3 );

    yq(1,element) = 0.5E+00 * ( y1 + y2 );
    yq(2,element) = 0.5E+00 * ( y2 + y3 );
    yq(3,element) = 0.5E+00 * ( y1 + y3 );

  end

  return
end
function [ wqe, xqe, yqe ] =  quad_e ( node_xy, element_node, ...
  element, element_num, nnodes, node_num, nqe )

%*****************************************************************************80
%
%% QUAD_E sets up quadrature information for a 13-point rule in a given element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer ELEMENT, the index of the element for which the quadrature
%    points are to be computed.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer NNODES, the number of nodes used to form one element.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer NQE, the number of points in the quadrature rule.
%    This is actually fixed at 13.
%
%    Output, real WQE(13), the quadrature weights.
%
%    Output, real XQE(NQE), YQE(NQE), the X and Y coordinates
%    of the quadrature points.
%
  for i = 1 : 3
    wqe(i) = 0.175615257433204E+00;
    ii = i + 3;
    wqe(ii) = 0.053347235608839E+00;
    ii = i + 6;
    iii = ii + 3;
    wqe(ii) = 0.077113760890257E+00;
    wqe(iii) = wqe(ii);
  end

  wqe(13) = -0.14957004446767E+00;

  z1 = 0.479308067841923E+00;
  z2 = 0.260345966079038E+00;
  z3 = 0.869739794195568E+00;
  z4 = 0.065130102902216E+00;
  z5 = 0.638444188569809E+00;
  z6 = 0.312865496004875E+00;
  z7 = 0.048690315425316E+00;

  ip1 = element_node(1,element);
  ip2 = element_node(2,element);
  ip3 = element_node(3,element);
  x1 = node_xy(1,ip1);
  x2 = node_xy(1,ip2);
  x3 = node_xy(1,ip3);
  y1 = node_xy(2,ip1);
  y2 = node_xy(2,ip2);
  y3 = node_xy(2,ip3);

  xqe( 1) = z1 * x1 + z2 * x2 + z2 * x3;
  yqe( 1) = z1 * y1 + z2 * y2 + z2 * y3;
  xqe( 2) = z2 * x1 + z1 * x2 + z2 * x3;
  yqe( 2) = z2 * y1 + z1 * y2 + z2 * y3;
  xqe( 3) = z2 * x1 + z2 * x2 + z1 * x3;
  yqe( 3) = z2 * y1 + z2 * y2 + z1 * y3;
  xqe( 4) = z3 * x1 + z4 * x2 + z4 * x3;
  yqe( 4) = z3 * y1 + z4 * y2 + z4 * y3;
  xqe( 5) = z4 * x1 + z3 * x2 + z4 * x3;
  yqe( 5) = z4 * y1 + z3 * y2 + z4 * y3;
  xqe( 6) = z4 * x1 + z4 * x2 + z3 * x3;
  yqe( 6) = z4 * y1 + z4 * y2 + z3 * y3;
  xqe( 7) = z5 * x1 + z6 * x2 + z7 * x3;
  yqe( 7) = z5 * y1 + z6 * y2 + z7 * y3;
  xqe( 8) = z5 * x1 + z7 * x2 + z6 * x3;
  yqe( 8) = z5 * y1 + z7 * y2 + z6 * y3;
  xqe( 9) = z6 * x1 + z5 * x2 + z7 * x3;
  yqe( 9) = z6 * y1 + z5 * y2 + z7 * y3;
  xqe(10) = z6 * x1 + z7 * x2 + z5 * x3;
  yqe(10) = z6 * y1 + z7 * y2 + z5 * y3;
  xqe(11) = z7 * x1 + z5 * x2 + z6 * x3;
  yqe(11) = z7 * y1 + z5 * y2 + z6 * y3;
  xqe(12) = z7 * x1 + z6 * x2 + z5 * x3;
  yqe(12) = z7 * y1 + z6 * y2 + z5 * y3;
  xqe(13) = ( x1 + x2 + x3 ) / 3.0E+00;
  yqe(13) = ( y1 + y2 + y3 ) / 3.0E+00;

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 ) 
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
%  For greater precision, try:
%
%     fprintf ( output_unit, '  %24,16f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %14f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function r8vec_print_some ( n, a, i_lo, i_hi, title )

%*****************************************************************************80
%
%% R8VEC_PRINT_SOME prints "some" of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines to print.
%
%    Input, string TITLE, an optional title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );

  for i = max ( 1, i_lo ) : min ( n, i_hi )
    fprintf ( 1, '  %8d  %12f\n', i, a(i) );
  end

  return
end
function value = rhs ( x, y, time )

%*****************************************************************************80
%
%% RHS gives the right-hand side of the differential equation.
%
%  Discussion:
%
%    The function specified here depends on the problem being
%    solved.  The user must be sure that RHS and EXACT_U are consistent.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2006
%
%  Parameters:
%
%    Input, real X, Y, the coordinates of a point
%    in the region, at which the right hand side of the
%    differential equation is to be evaluated.
%
%    Input, real TIME, the current time.
%
%    Output, real VALUE, the value of the right
%    hand side of the differential equation at (X,Y).
%
  ut =            - sin ( pi * x ) .* sin ( pi * y ) * exp ( - time );
  uxx = - pi * pi * sin ( pi * x ) .* sin ( pi * y ) * exp ( - time );
  uyy = - pi * pi * sin ( pi * x ) .* sin ( pi * y ) * exp ( - time );

  value = ut - uxx - uyy;

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function triangulation_order6_plot ( file_name, node_num, node_xy, tri_num, ...
  triangle_node, node_show, triangle_show )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_PLOT plots a 6-node triangulation of a pointset.
%
%  Discussion:
%
%    The triangulation is most usually a Delaunay triangulation,
%    but this is not necessary.
%
%    In a six node triangulation, it is assumed that nodes 1, 2, and 3
%    are the vertices of the triangles, and that nodes 4, 5, and 6
%    lie between 1 and 2, 2 and 3, and 3 and 1 respectively.
%
%    This routine has been specialized to deal correctly ONLY with
%    a mesh of 6 node elements, with the property that starting
%    from local node 1 and traversing the edges of the element will
%    result in encountering local nodes 1, 4, 2, 5, 3, 6 in that order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character FILE_NAME(*), the name of the output file.
%
%    Input, integer NODE_NUM, the number of points.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer TRI_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(6,TRI_NUM), lists, for each triangle,
%    the indices of the points that form the vertices and midsides 
%    of the triangle.
%
%    Input, logical NODE_SHOW:
%    0, do not show the nodes.
%    1, show the nodes.
%    2, show the nodes, and label them.
%
%    Input, logical TRIANGLE_SHOW, 
%    0, do not show the trangles.
%    1, show the triangles.
%    2, show the triangles, and label them.
%
  order = [ 1, 4, 2, 5, 3, 6 ];
  x_ps_max = 576;
  x_ps_max_clip = 594;
  x_ps_min = 36;
  x_ps_min_clip = 18;
  y_ps_max = 666;
  y_ps_max_clip = 684;
  y_ps_min = 126;
  y_ps_min_clip = 108;
%
%  We need to do some figuring here, so that we can determine
%  the range of the data, and hence the height and width
%  of the piece of paper.
%
  x_max = max ( node_xy(1,1:node_num) );
  x_min = min ( node_xy(1,1:node_num) );
  x_scale = x_max - x_min;
  
  x_max = x_max + 0.05 * x_scale;
  x_min = x_min - 0.05 * x_scale;
  x_scale = x_max - x_min;

  y_max = max ( node_xy(2,1:node_num) );
  y_min = min ( node_xy(2,1:node_num) );
  y_scale = y_max - y_min;

  y_max = y_max + 0.05 * y_scale;
  y_min = y_min - 0.05 * y_scale;
  y_scale = y_max - y_min;

  if ( x_scale < y_scale )

    delta = round ( ( x_ps_max - x_ps_min ) ...
      * ( y_scale - x_scale ) / ( 2.0 * y_scale ) );

    x_ps_max = x_ps_max - delta;
    x_ps_min = x_ps_min + delta;

    x_ps_max_clip = x_ps_max_clip - delta;
    x_ps_min_clip = x_ps_min_clip + delta;

    x_scale = y_scale;

  elseif ( y_scale < x_scale )

    delta = round ( ( y_ps_max - y_ps_min ) ...
      * ( x_scale - y_scale ) / ( 2.0 * x_scale ) );

    y_ps_max      = y_ps_max - delta;
    y_ps_min      = y_ps_min + delta;

    y_ps_max_clip = y_ps_max_clip - delta;
    y_ps_min_clip = y_ps_min_clip + delta;

    y_scale = x_scale;

  end
%
%  Plot the triangulation.
%
  file_unit = fopen ( file_name, 'wt' );

  if ( file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_ORDER3_PLOT - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'TRIANGULATION_ORDER3_PLOT - Fatal error!' );
  end

  fprintf ( file_unit, '%%!PS-Adobe-3.0 EPSF-3.0\n' );
  fprintf ( file_unit, '%%%%Creator: triangulation_order6_plot.m\n' );
  fprintf ( file_unit, '%%%%Title: %s\n', file_name );
  fprintf ( file_unit, '%%%%Pages: 1\n' );
  fprintf ( file_unit, '%%%%BoundingBox:  %d  %d  %d  %d\n', ...
    x_ps_min, y_ps_min, x_ps_max, y_ps_max );
  fprintf ( file_unit, '%%%%Document-Fonts: Times-Roman\n' );
  fprintf ( file_unit, '%%%%LanguageLevel: 1\n' );
  fprintf ( file_unit, '%%%%EndComments\n' );
  fprintf ( file_unit, '%%%%BeginProlog\n' );
  fprintf ( file_unit, '/inch {72 mul} def\n' );
  fprintf ( file_unit, '%%%%EndProlog\n' );
  fprintf ( file_unit, '%%%%Page: 1 1\n' );
  fprintf ( file_unit, 'save\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the RGB color to very light gray.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '0.900  0.900  0.900 setrgbcolor\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Draw a gray border around the page.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, 'stroke\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the RGB color to black.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '0.000  0.000  0.000 setrgbcolor\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the font and its size.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '/Times-Roman findfont\n' );
  fprintf ( file_unit, '0.50 inch scalefont\n' );
  fprintf ( file_unit, 'setfont\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Print a title.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%210  702  moveto\n' );
  fprintf ( file_unit, '%%(Triangulation)  show\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Define a clipping polygon.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, 'clip newpath\n' );
%
%  Draw the nodes.
%
  if ( node_num <= 200 )
  {
    circle_size = 5;
  }
  else if ( node_num <= 500 )
  {
    circle_size = 4;
  }
  else if ( node_num <= 1000 )
  {
    circle_size = 3;
  }
  else if ( node_num <= 5000 )
  {
    circle_size = 2;
  }
  else
  {
    circle_size = 1;
  }
  if ( 1 <= node_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Draw filled dots at the nodes.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to blue.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.000  0.150  0.750 setrgbcolor\n' );
    fprintf ( file_unit, '%%\n' );

    for node = 1 : node_num

      x_ps = floor ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, ...
        '  newpath  %d  %d  %d 0 360 arc closepath fill\n', ...
        x_ps, y_ps, circle_size );

    end

  end 
%
%  Label the nodes.
%
  if ( 2 <= node_show )
      
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Label the nodes.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to darker blue.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.000  0.250  0.850 setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.20 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );
    fprintf ( file_unit, '%%\n' );

    for node = 1 : node_num

      x_ps = floor ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, '  %d  %d  moveto (%d) show\n', ...
        x_ps, y_ps+5, node );

    end

  end
%
%  Draw the triangles.
%
  if ( 1 <= triangle_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to red.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.900  0.200  0.100 setrgbcolor\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%  Draw the triangles.\n' );
    fprintf ( file_unit, '%%\n' );

    for  triangle = 1 : tri_num

      fprintf ( file_unit, 'newpath\n' );

      node = triangle_node(order(1),triangle);

      x_ps = floor ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, '%d  %d  moveto\n', x_ps, y_ps );
  
      for i = 1 : 6

        ip1 = mod ( i, 6 ) + 1;
        node = triangle_node(order(ip1),triangle);

        x_ps = floor ( ...
          ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
          + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
          / ( x_max                   - x_min ) );

        y_ps = floor ( ...
          ( ( y_max - node_xy(2,node)         ) * y_ps_min ...
          + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
          / ( y_max                   - y_min ) );

        fprintf ( file_unit, '%d  %d  lineto\n', x_ps, y_ps );

      end

      fprintf ( file_unit, 'stroke\n' );

    end

  end
%
%  Label the triangles.
%
  if ( 2 <= triangle_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Label the triangles:\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to darker red:\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, ' 0.950  0.250  0.150  setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.30 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );

    for triangle = 1 : tri_num

      ave_x = 0.0;
      ave_y = 0.0;

      for i = 1 : 6
        node = triangle_node(i,triangle);
        ave_x = ave_x + node_xy(1,node);
        ave_y = ave_y + node_xy(2,node);
      end

      ave_x = ave_x / 6.0;
      ave_y = ave_y / 6.0;

      x_ps = floor ( ...
        ( ( x_max - ave_x         ) * x_ps_min ...
        + (         ave_x - x_min ) * x_ps_max ) ...
        / ( x_max         - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - ave_y         ) * y_ps_min ...
        + (         ave_y - y_min ) * y_ps_max ) ...
        / ( y_max         - y_min ) );

      fprintf ( file_unit, '%d  %d  moveto\n', x_ps, y_ps );
      fprintf ( file_unit, '(%d) show\n', triangle );

    end

  end

  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'restore showpage\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%% End of page\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%%%Trailer\n' );
  fprintf ( file_unit, '%%%%EOF\n' );

  fclose ( file_unit );

  return
end
function [ node_xy ] = xy_set ( nx, ny, node_num, xl, xr, yb, yt  )

%*****************************************************************************80
%
%% XY_SET sets the XY coordinates of the nodes.
%
%  Discussion:
%
%    The nodes are laid out in an evenly spaced grid, in the unit square.
%
%    The first node is at the origin.  More nodes are created to the
%    right until the value of X = 1 is reached, at which point
%    the next layer is generated starting back at X = 0, and an
%    increased value of Y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2005
%
%  Parameters:
%
%    Input, integer NX, NY, the number of elements in the X and
%    Y direction.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real XL, XR, YB, YT, the X coordinates of
%    the left and right sides of the rectangle, and the Y coordinates
%    of the bottom and top of the rectangle.
%
%    Output, real NODE_XY(2,NODE_NUM), the nodes.
%
  for j = 1 : 2*ny-1
    for i = 1 : 2*nx - 1

      node_xy(1,i+(j-1)*(2*nx-1)) =    ...
        ( ( 2 * nx - i - 1 ) * xl   ...
        + (          i - 1 ) * xr ) ...
        / ( 2 * nx     - 2 );

      node_xy(2,i+(j-1)*(2*nx-1)) =    ...
        ( ( 2 * ny - j - 1 ) * yb   ...
        + (          j - 1 ) * yt ) ...
        / ( 2 * ny     - 2 );

    end
  end

  return
end