%% MAIN is the main program for FEM1D_NONLINEAR.
%
%  Discussion:
%
%    FEM1D_NONLINLEAR solves a nonlinear one dimensional boundary value problem.
%
%    The differential equation has the form:
%
%      -d/dx ( p(x) du/dx ) + q(x) * u + u * u' = f(x)
%
%    The finite-element method uses piecewise linear basis functions.
%
%    Here U is an unknown scalar function of X defined on the
%    interval [XL,XR], and P, Q and F are given functions of X.
%
%    The values of U or U' at XL and XR are also specified.
%
%    Sample problem #1:
%
%    u(x) = x, 
%    p(x) = 1, 
%    q(x) = 0, 
%    f(x) = x, 
%    u(0) = 0, 
%    u'(1) = 1.
%
%    The code should solve this problem exactly.
%
%    Sample problem #2:
%
%    u(x)  = 2*(1-cos(0.5*pi*x))/pi, 
%    p(x)  = 1, 
%    q(x)  = 0,
%    f(x)  = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
%    u(0)  = 0, 
%    u'(1) = 1.
%
%  Modified:
%
%    02 November 2006
%
%  Parameters:
%
%    real ADIAG(NU).
%    ADIAG(I) is the "diagonal" coefficient of the I-th
%    equation in the linear system.  That is, ADIAG(I) is
%    the coefficient of the I-th unknown in the I-th equation.
%
%    real ALEFT(NU).
%    ALEFT(I) is the "left hand" coefficient of the I-th
%    equation in the linear system.  That is, ALEFT(I) is the
%    coefficient of the (I-1)-th unknown in the I-th equation.
%    There is no value in ALEFT(1), since the first equation
%    does not refer to a "0-th" unknown.
%
%    real ARITE(NU).
%    ARITE(I) is the "right hand" coefficient of the I-th
%    equation in the linear system.  ARITE(I) is the coefficient
%    of the (I+1)-th unknown in the I-th equation.  There is
%    no value in ARITE(NU) because the NU-th equation does not
%    refer to an "NU+1"-th unknown.
%
%    real H(N), the length of the subintervals.  
%
%    integer IBC.
%    IBC declares what the boundary conditions are.
%    1, at the left endpoint, U has the value UL,
%       at the right endpoint, U' has the value UR.
%    2, at the left endpoint, U' has the value UL,
%       at the right endpoint, U has the value UR.
%    3, at the left endpoint, U has the value UL,
%       and at the right endpoint, U has the value UR.
%    4, at the left endpoint, U' has the value UL,
%       at the right endpoint U' has the value UR.
%
%    integer IMAX.
%    The number of Newton iterations to carry out.
%
%    integer INDX(1:N+1).
%    For a node I, INDX(I) is the index of the unknown
%    associated with node I.
%    If INDX(I) is equal to -1, then no unknown is associated
%    with the node, because a boundary condition fixing the
%    value of U has been applied at the node instead.
%    Unknowns are numbered beginning with 1.
%    If IBC is 2 or 4, then there is an unknown value of U
%    at node 0, which will be unknown number 1.  Otherwise,
%    unknown number 1 will be associated with node 1.
%    If IBC is 1 or 4, then there is an unknown value of U
%    at node N, which will be unknown N or N+1,
%    depending on whether there was an unknown at node 0.
%
%    integer NL.
%    The number of basis functions used in a single
%    subinterval.  (NL-1) is the degree of the polynomials
%    used.  For this code, NL is fixed at 2, meaning that
%    piecewise linear functions are used as the basis.
%
%    integer NODE(NL,N).
%    For each subinterval I:
%    NODE(1,I) is the number of the left node, and
%    NODE(2,I) is the number of the right node.
%
%    integer NPRINT.
%    The number of points at which the computed solution
%    should be printed out when compared to the exact solution.
%
%    integer NQUAD.
%    The number of quadrature points used in a subinterval.
%    This code uses NQUAD = 1.
%
%    integer N.
%    The number of subintervals into which the interval
%    [XL,XR] is broken.
%
%    integer NU.
%    NU is the number of unknowns in the linear system.
%    Depending on the value of IBC, there will be N-1,
%    N, or N+1 unknown values, which are the coefficients
%    of basis functions.
%
%    integer PROBLEM, indicates which problem to be solved.
%    * 1, u = x, p=1, q=0, f=x, u(0)=0, u'(1)=1.
%    * 2, u = 2*(1-cos(0.5*pi*x))/pi, p=1, q=0,
%    f = -0.5*pi*cos(0.5*pi*x) + 2*sin(0.5*pi*x)*(1-cos(0.5*pi*x)/pi
%    u(0) = 0, u'(1)=1.
%
%    real RHS(NU), the right hand side of the linear equations
%    for the Picard or Newton iterations.
%
%    real U(NU).
%    U contains the value of the solution from the previous iteration,
%    and is used in ASSEMBLE to add correction terms to the
%    matrix and right hand side.
%
%    real UL.
%    If IBC is 1 or 3, UL is the value that U is required
%    to have at X = XL.
%    If IBC is 2 or 4, UL is the value that U' is required
%    to have at X = XL.
%
%    real UR.
%    If IBC is 2 or 3, UR is the value that U is required
%    to have at X = XR.
%    If IBC is 1 or 4, UR is the value that U' is required
%    to have at X = XR.
%
%    real XL.
%    XL is the left endpoint of the interval over which the
%    differential equation is being solved.
%
%    real XN(1:N+1).
%    XN(I) is the location of the I-th node.  XN(1) is XL,
%    and XN(N+1) is XR.
%
%    real XNEWT.
%    A factor which is 0 for the first three iterations, and
%    1 thereafter.  It multiplies some terms in ASSEMBLE,
%    and thus is used to "weaken" the nonlinearity for
%    the first few steps.
%
%    real XQUAD(N)
%    XQUAD(I) is the location of the single quadrature point
%    in interval I.
%
%    real XR.
%    XR is the right endpoint of the interval over which the
%    differential equation is being solved.
%
  n = 10;
  nl = 2;

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

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_NONLINEAR\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve a nonlinear boundary value problem:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    -d/dx (p(x) du/dx) + q(x)*u  + u*u'' =  f(x)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  on an interval [xl,xr], with the values of\n' );
  fprintf ( 1, '  u or u'' specified at xl and xr.\n' );
%
%  Initialize variables that define the problem.
%
  [ ibc, imax, nprint, nquad, problem, ul, ur, xl, xr ] = init ( 'DUMMY' );
%
%  Compute the quantities that describe the geometry of the problem.
%
  [ h, indx, node, nu, xn, xquad ] = geometry ( ibc, nl, n, xl, xr );
%
%  The initial solution estimate is 0.
%
  u(1:nu) = 0.0;
%
%  Begin the Newton iteration.
%
  for i = 1 : imax
%
%  Is it time for full nonlinear Newton iteration?
%  For the first 3 steps, we use a weaker iteration that is less
%  likely to diverge.
%
    if ( i <= 5 )
      [ adiag, aleft, arite, rhs ] = assemble_picard ( u, h, indx, n, nl, node, ...
        nquad, nu, problem, ul, ur, xn, xquad );
    else
      [ adiag, aleft, arite, rhs ] = assemble_newton ( u, h, indx, n, nl, node, ...
        nquad, nu, problem, ul, ur, xn, xquad );
    end
%
%  Print out the linear system, just once.
%
    if ( i == 1 || i == 6 )
      prsys ( adiag, aleft, arite, rhs, nu );
    end
%
%  Solve the linear system A * Unew = RHS.
%
    u = solve ( adiag, aleft, arite, rhs, nu );
%
%  Print the current solution.
%
    output ( u, ibc, indx, n, nu, ul, ur, xn );

  end
%
%  Compare the solution to the exact solution.
%
  compare ( u, indx, n, nl, node, nprint, nu, problem, ul, ur, xl, xn, xr );

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

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