%% MAIN is the main routine of FEM2D_POISSON.
%
%  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.
%
%  Modified:
%
%    24 March 2005
%
%  Local parameters:
%
%    Local, real A(3*IB+1,NUNK), 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(NUNK), the right hand side.
%
%    Local, integer IB, the half-bandwidth of the matrix.
%
%    Local, integer INDX(NODE_NUM), gives the index of the unknown quantity
%    associated with the given node.
%
%    Local, integer NNODES, the number of nodes used to form one element.
%
%    Local, integer NQ, the number of quadrature points used for assembly.
%
%    Local, integer NQE, the number of quadrature points used for error
%    estimation.
%
%    Local, integer NUNK, the number of unknowns.
%
%    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 NODE_XY(2,NODE_NUM), the nodes.
%
%    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.
%
  clear

  nnodes = 6;
  nq = 3;
  nqe = 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_POISSON (MATLAB version):\n' );
  fprintf ( 1, '  Solution of the Poisson equation on a unit box\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, with piecewise\n' );
  fprintf ( 1, '  quadratic basis functions on 6 node triangular\n' );
  fprintf ( 1, '  elements.\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 );

  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 are boundary nodes and which have a
%  finite element unknown.  Then set the boundary values.
%
  [ indx, nunk ] = indx_set ( nx, ny, node_num );

  fprintf ( 1, '  Number of unknowns =       %d\n', nunk );
%
%  Determine the bandwidth of the coefficient matrix.
%
  ib = bandwidth ( nnodes, element_num, element_node, node_num, indx );

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

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

    nodes_plot ( node_eps_file_name, node_num, node_xy, node_label );

    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_POISSON:\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 = 'fem2d_poisson_nodes.txt';

  nodes_write ( node_num, node_xy, node_txt_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON:\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 = 'fem2d_poisson_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_POISSON:\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 = 'fem2d_poisson_elements.txt';

  element_write ( nnodes, element_num, element_node, ...
    triangulation_txt_file_name );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON:\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' );
%
%  Assemble the coefficient matrix A and the right-hand side F of the
%  finite element equations.
%
  [ a, f ] = assemble ( node_num, node_xy, nnodes, element_num, ...
    element_node, nq, wq, xq, yq, element_area, indx, ib, nunk );

  dgb_print_some ( nunk, nunk, ib, ib, a, 1, 1, 5, 5, ...
    '  Initial 5 x 5 block of matrix A' );

  r8vec_print_some ( nunk, f, 1, 10, '  Part of right hand side vector:' );
%
%  Modify the coefficient matrix and right hand side to account for
%  boundary conditions.
%
  [ a, f ] = boundary ( nx, ny, node_num, node_xy, indx, ib, nunk, a, f );

  dgb_print_some ( nunk, nunk, ib, ib, a, 1, 1, 5, 5, ...
    '  A after boundary adjustments' );

  r8vec_print_some ( nunk, f, 1, 10, '  Part of right hand side vector:' );
%
%  Solve the linear system using a banded solver.
%
  [ a_lu, pivot, info ] = dgb_fa ( nunk, ib, ib, a );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM2D_POISSON - 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_POISSON - Fatal error!' );
  end

  job = 0;
  c = dgb_sl ( nunk, ib, ib, a_lu, pivot, f, job );

  r8vec_print_some ( nunk, f, 1, 10, '  Part of the solution vector:' );
%
%  Calculate error using 13 point quadrature rule.
%
  [ el2, eh1 ] = errors ( element_area, element_node, indx, node_xy, ...
    c, element_num, nnodes, nqe, nunk, node_num );

  compare ( node_num, node_xy, indx, nunk, c );
%
%  Write an ASCII file that can be read into MATLAB.
%
  solution_txt_file_name = 'fem2d_poisson_solution.txt';

  solution_write ( c, indx, node_num, nunk, solution_txt_file_name, ...
    node_xy );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_POISSON:\n' );
  fprintf ( 1, '  Wrote an ASCII solution file\n' );
  fprintf ( 1, '    "%s"\n', solution_txt_file_name );
  fprintf ( 1, '  of the form\n' );
  fprintf ( 1, '    U ( X(I), Y(I) )\n' );
  fprintf ( 1, '  which can be used for plotting.\n' );

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

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