function fem2d_poisson_rectangle ( )

%*****************************************************************************80
%
%% MAIN is the main routine of FEM2D_POISSON_RECTANGLE.
%
%  Discussion:
%
%    This program solves
%
%      -Laplacian U(X,Y) = F(X,Y)
%
%    in a rectangular region in the plane.
%
%    Along the boundary of the region, homogeneous Dirichlet conditions
%    are imposed:
%
%      U(Xb,Yb) = G(X,Y)
%
%    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:
%
%    13 October 2009
%
%  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(ELEMENT_ORDER,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, integer ELEMENT_ORDER, the number of nodes used to form one element.
%
%    Local, real F(NODE_NUM), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, integer NQ, the number of quadrature points used for assembly.
%
%    Local, integer NODE_NUM, the number of nodes.
%
%    Local, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Local, integer NQE, the number of quadrature points used for error
%    estimation.
%
%    Local, integer NX, the number of points in the X direction.
%
%    Local, integer NY, the number of points in the Y direction.
%
%    Local, real WQ(NQ), 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(NQ,ELEMENT_NUM), YQ(NQ,ELEMENT_NUM), the
%    coordinates of the quadrature points in each element.
%
  element_order = 6;
  nq = 3;
  nqe = 13;
  nx = 97;
  ny = 97;

  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_POISSON_RECTANGLE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Solution of the Poisson equation in the unit square\n' );
  fprintf ( 1, '  in 2 dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - Uxx - Uyy = F(x,y) in the box\n' );
  fprintf ( 1, '       U(x,y) = G(x,y) on the boundary.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The finite element method is used.\n' );
  fprintf ( 1, '  Triangular elements are used, of order %d.\n', element_order );
  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 );

  element_node = grid_t6 ( nx, ny, element_order, element_num );
%
%  Set the quadrature rule for assembly.
%
  [ wq, xq, yq ] = quad_a ( node_xy, element_node, element_num, node_num, ...
    element_order );
%
%  Determine the areas of the elements.
%
  element_area = area_set ( node_xy, element_num, element_node );
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( element_order, element_num, element_node );

  fprintf ( 1, '  The matrix bandwidth is    %d\n', 3 * ib + 1 );
%
%  Assemble the coefficient matrix A and the right-hand side F of the
%  finite element equations.
%
  tic

  [ a, f ] = assemble ( node_num, node_xy, element_order, element_num, ...
    element_node, nq, wq, xq, yq, element_area );

  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix assembly time = %f seconds.\n', wtime );
%
%  Modify the coefficient matrix and right hand side to account for
%  boundary conditions.
%
  [ a, f ] = boundary ( nx, ny, node_num, node_xy, ib, a, f );
%
%  Solve the linear system using a banded solver.
%
  tic;

  c = a \ f';

  wtime = toc;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Linear system factor and solve time = %f seconds.\n', wtime );

  if ( 0 )
    r8vec_print_some ( node_num, f, 1, 5, '  Part of the solution vector:' );
  end
%
%  Calculate error using 13 point quadrature rule.
%
  errors ( element_area, element_node, node_xy, ...
    c, element_num, element_order, nqe, node_num );

  if ( 0 )
    compare ( node_num, node_xy, c );
  end
%
%  Write ASCII files.
%
  node_filename = 'rectangle_nodes.txt';

  r8mat_write ( node_filename, 2, node_num, node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Wrote the node file "%s"\n', node_filename );

  triangulation_filename = 'rectangle_elements.txt';

  i4mat_write ( triangulation_filename, element_order, element_num, element_node );

  fprintf ( 1, '  Wrote the element file "%s"\n', triangulation_filename );

  solution_filename = 'rectangle_solution.txt';

  r8mat_write ( solution_filename, 1, node_num, c );

  fprintf ( 1, '  Wrote the solution file "%s"\n', solution_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON_RECTANGLE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

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

  return
end
function element_area = area_set ( node_xy, 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 ELEMENT_ORDER, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,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.
%
  element_area = zeros ( element_num, 1 );

  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.5 * abs ...
      ( y1 * ( x2 - x3 ) ...
      + y2 * ( x3 - x1 ) ...
      + y3 * ( x1 - x2 ) );

  end

  return
end
function [ a, f ] = assemble ( node_num, node_xy, element_order, ...
  element_num, element_node, nq, wq, xq, yq, element_area )

%*****************************************************************************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:
%
%    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 ELEMENT_ORDER, the number of nodes used to form one element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer NQ, the number of quadrature points used in assembly.
%
%    Input, real WQ(NQ), quadrature weights.
%
%    Input, real XQ(NQ,ELEMENT_NUM), YQ(NQ,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.
%
%    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.
%

%
%  Use MATLAB's sparse matrix storage format.
%
  if ( element_order == 3 )
    nz_max = 7 * node_num;
  else
    nz_max = 19 * node_num;
  end

  a = sparse ( [], [], [], node_num, node_num, nz_max );
%
%  Initialize the arrays to zero.
%
  f(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 : nq

      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.
%
      for test = 1 : element_order

        i = element_node(test,element);

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

        f(i) = f(i) + w * rhs ( x, y ) * bi;
%
%  Consider a basis function used to form the value of the solution function.
%
%  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 : element_order

          j = element_node(basis,element);

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

          aij = dbidx * dbjdx + dbidy * dbjdy;

          a(i,j) = a(i,j) + w * aij;

        end

      end

    end

  end

  return
end
function nhba = bandwidth ( element_order, element_num, element_node )

%*****************************************************************************80
%
%% BANDWIDTH determines the bandwidth of the coefficient matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_ORDER, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,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 : element_order
      i = element_node(iln,element);
      if ( 0 < i )
        for jln = 1 : element_order
          j = element_node(jln,element);
          nhba = max ( nhba, j - i );
        end
      end
    end
  end

  return
end
function [ a, f ] = boundary ( nx, ny, node_num, node_xy, ib, a, f )

%*****************************************************************************80
%
%% BOUNDARY modifies the linear system for boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2005
%
%  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 NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    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.
%

%
%  Consider each node.
%
  node = 0;

  for row = 1 : 2 * ny - 1

    for col = 1 : 2 * nx - 1

      node = node + 1;

      if ( row == 1 || row == 2 * ny - 1 || col == 1 || col == 2 * nx - 1 )

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

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

        for j = jlo : jhi
          a(node,j) = 0.0;
        end

        a(node,node) = 1.0;

        f(node) = u;

      end

    end
  end

  return
end
function compare ( node_num, node_xy, f )

%*****************************************************************************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:
%
%    17 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, real F(NODE_NUM), the solution vector of the finite
%    element system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMPARE:\n' );
  fprintf ( 1, '  Compare computed and exact solutions at the nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         X           Y          U           U\n' );
  fprintf ( 1, '                             computed     exact\n' );
  fprintf ( 1, '\n' );

  for node = 1 : node_num

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

    u = exact ( x, y );

    uh = f(node);

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

  end

  return
end
function [ el2, eh1 ] = errors ( element_area, element_node, node_xy, ...
  f, element_num, element_order, nqe, node_num )

%*****************************************************************************80
%
%% ERRORS calculates the L2 and H1-seminorm errors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ELEMENT_AREA(ELEMENT_NUM), the area of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,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 F(NODE_NUM), the coefficients of the solution.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_ORDER, the number of nodes used to form one element.
%
%    Input, integer NQE, the number of points in the quadrature rule.
%    This is actually fixed at 13.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    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.0E+00;
  eh1 = 0.0E+00;

  for element = 1 : element_num

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

    for iq = 1 : nqe

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

      uh = 0.0E+00;
      dudxh = 0.0E+00;
      dudyh = 0.0E+00;

      for in1 = 1 : element_order

        i = element_node(in1,element);

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

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

      end

      [ u, dudx, dudy ] = exact ( x, y );

      el2 = el2 + ( uh - u )^2 * ar;
      eh1 = eh1 + ( ( dudxh - dudx )^2 + ( dudyh - dudy )^2 ) * ar;

    end

  end

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '*********************************************\n' );
  fprintf ( 1, '*                                           *\n' );
  fprintf ( 1, '*  ERRORS:                                  *\n' );
  fprintf ( 1, '*    L2 error =          %14f     *\n', el2 );
  fprintf ( 1, '*    H1-seminorm error = %14f     *\n', eh1 );
  fprintf ( 1, '*                                           *\n' );
  fprintf ( 1, '*********************************************\n' );

  return
end
function [ u, dudx, dudy ] = exact ( x, y )

%*****************************************************************************80
%
%% EXACT calculates the exact solution and its first derivatives.
%
%  Discussion:
%
%    The function specified here depends on the problem being
%    solved.  The user must make sure that EXACT and RHS are consistent.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the coordinates of a point
%    in the region, at which the exact solution is to be evaluated.
%
%    Output, real U, DUDX, DUDY, the value of
%    the exact solution U and its derivatives dUdX
%    and dUdY at the point (X,Y).
%
  u    =      sin ( pi * x ) * sin ( pi * y ) + x;
  dudx = pi * cos ( pi * x ) * sin ( pi * y ) + 1.0;
  dudy = pi * sin ( pi * x ) * cos ( pi * y );

  return
end
function element_node = grid_t6 ( nx, ny, element_order, 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 ELEMENT_ORDER, the number of local nodes per element.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the index of the I-th node of the J-th element.
%
  element_node = zeros ( element_order, element_num );

  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!' );
  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 [ b, dbdx, dbdy ] = qbf ( x, y, element, inode, node_xy, ...
  element_node, element_num, element_order, 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(ELEMENT_ORDER,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 ELEMENT_ORDER, 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.0 * ( 1.0 - r - s ) * ( 0.5 - r - s );
    dbdr = - 3.0 + 4.0 * r + 4.0 * s;
    dbds = - 3.0 + 4.0 * r + 4.0 * s;

  elseif ( inode == 2 )

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

  elseif ( inode == 3 )

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

  elseif ( inode == 4 )

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

  elseif ( inode == 5 )

    b    =   4.0 * r * s;
    dbdr =   4.0 * s;
    dbds =   4.0 * r;

  elseif ( inode == 6 )

    b    =   4.0 * s * ( 1.0 - r - s );
    dbdr = - 4.0 * s;
    dbds =   4.0 - 4.0 * r - 8.0 * 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, element_order )

%*****************************************************************************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(ELEMENT_ORDER,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 ELEMENT_ORDER, 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 = zeros ( 3, 1 );
  xq = zeros ( 3, element_num );
  yq = zeros ( 3, element_num );

  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, element_order, 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(ELEMENT_ORDER,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 ELEMENT_ORDER, 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.
%
  wqe = zeros ( 13, 1 );

  for i = 1 : 3
    wqe(i) = 0.175615257433204;
    ii = i + 3;
    wqe(ii) = 0.053347235608839;
    ii = i + 6;
    iii = ii + 3;
    wqe(ii) = 0.077113760890257;
    wqe(iii) = wqe(ii);
  end

  wqe(13) = -0.14957004446767;

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

  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.0;
  yqe(13) = ( y1 + y2 + y3 ) / 3.0;

  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!' );
  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, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  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 )

%*****************************************************************************80
%
%% RHS gives the right-hand side of the differential equation.
%
%  Discussion:
%
%    The function specified here depends on the problem being
%    solved.  This is one of the routines that a user will
%    normally want to change.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  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.
%
%    Output, real VALUE, the value of the right
%    hand side of the differential equation at (X,Y).
%
  value = 2.0 * pi * pi * sin ( pi * x ) * sin ( pi * y );

  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 [ 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
%
%  Author:
%
%    John Burkardt
%
%  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.
%
  node_xy = zeros ( 2, node_num );

  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
