function [ ns, xs, vs ] = paraheat_pwc_1d ( nc, xc, vc )

%*****************************************************************************80
%
%% paraheat_pwc_1d demonstrates the parameterized 1d heat equation solver.
%
%  Discussion:
%
%    paraheat_pwc_1d is a first draft for a project that studies parameterized
%    heat equations for which the solution is sampled at a few locations.
%
%    The steady heat equation is modeled, with piecewise constant diffusivity.
%    The diffusivity parameters are regarded as the input, and
%    the sensor readings as the output, of a sort of function y=f(x).
%
%    We may be interested in attempting to estimate values of x given
%    the values of y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NC, the number of values used to define diffusivity.
%
%    Input, real XC(NC-1), the breakpoints to define diffusivity.
%
%    Input, real VC(NC), the diffusivity values.
%
%    Output, integer NS, the number of sensors.
%
%    Output, real XS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the sensor values.
%
%  Local parameters:
%
%    Local, integer NN, the number of nodes.
%
%    Local, real XN(NN), the coordinates of the nodes.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'paraheat_pwc_1d\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Parameterized heat equation in 1D.\n' );
  fprintf ( 1, '  Steady heat problem with piecewise constant diffusivity.\n' );
  fprintf ( 1, '  Diffusivity parameters are input.\n' );
  fprintf ( 1, '  Sensor readings are output.\n' );
%
%  Define the region.
%
  [ xmin, xmax ] = region_box ( );
%
%  Define the mesh.
%
  nn = 21;
  xn = linspace ( xmin, xmax, nn );
%
%  Define the sensor locations.
%
  ns = 3;
  xs = [ 0.22, 0.43, 0.76 ];
%
%  Solve the heat equation, and return the solution value at sensor locations.
%
  vs = heat1d ( nc, xc, vc, nn, xn, ns, xs );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'paraheat_pwc_1d\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function left = bracketer ( n, t, tval )

%*****************************************************************************80
%
%% bracketer finds the interval containing or nearest a given value.
%
%  Discussion:
%
%    The routine always returns the index LEFT of the sorted array
%    T with the property that either
%    *  T is contained in the interval [ T(LEFT), T(LEFT+1) ], or
%    *  T < T(LEFT) = T(1), or
%    *  T > T(LEFT+1) = T(N).
%
%    The routine is useful for interpolation problems, where
%    the abscissa must be located within an interval of data
%    abscissas for interpolation, or the "nearest" interval
%    to the (extreme) abscissa must be found so that extrapolation
%    can be carried out.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, length of the input array.
%
%    Input, real T(N), an array that has been sorted into ascending order.
%
%    Input, real TVAL, a value to be bracketed by entries of T.
%
%    Output, integer LEFT, is set so that the interval [ T(LEFT), T(LEFT+1) ]
%    is the closest to TVAL; it either contains TVAL, or else TVAL
%    lies outside the interval [ T(1), T(N) ].
%

%
%  Check the input data.
%
  if ( n < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'bracketer - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 2.\n' );
    error ( 'bracketer - Fatal error!' );
  end
%
%  Set LEFT to the middle value.
%
  left = floor ( ( n + 1 ) / 2 );
%
%  CASE 1: TVAL < T(LEFT):
%  Search for TVAL in [T(I), T(I+1)] for intervals I = 1 to LEFT-1.
%
  if ( tval < t(left) )

    if ( left == 1 )
      left = 1;
      return
    elseif ( left == 2 )
      left = 1;
      return
    elseif ( t(left-1) <= tval )
      left = left - 1;
      return
    elseif ( tval <= t(2) )
      left = 1;
      return
    end
%
%  ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = 2 to LEFT-2.
%
    low = 2;
    high = left - 2;

    while ( true )

      if ( low == high )
        left = low;
        return
      end

      mid = floor ( ( low + high + 1 ) / 2 );

      if ( t(mid) <= tval )
        low = mid;
      else
        high = mid - 1;
      end

    end
%
%  CASE2: T(LEFT+1) < TVAL:
%  Search for TVAL in {T(I),T(I+1)] for intervals I = LEFT+1 to N-1.
%
  elseif ( t(left+1) < tval )

    if ( left == n - 1 )
      return
    elseif ( left == n - 2 )
      left = left + 1;
      return
    elseif ( tval <= t(left+2) )
      left = left + 1;
      return
    elseif ( t(n-1) <= tval )
      left = n - 1;
      return
    end
%
%  ...Binary search for TVAL in [T(I), T(I+1)] for intervals I = LEFT+2 to N-2.
%
    low = left + 2;
    high = n - 2;

    while ( true )

      if ( low == high )
        left = low;
        return
      end

      mid = floor ( ( low + high + 1 ) / 2 );

      if ( t(mid) <= tval )
        low = mid;
      else
        high = mid - 1;
      end

    end
%
%  CASE3: T(LEFT) <= TVAL <= T(LEFT+1):
%  T is in [T(LEFT), T(LEFT+1)], as the user said it might be.
%
  else

  end

  return
end
function vp = diffusivity_pwc_1d ( nc, xc, vc, np, xp )

%*****************************************************************************80
%
%% diffusivity_pwc_1d: piecewise constant diffusivity function in 1D.
%
%  Discussion:
%
%    A piecewise constant function is defined over NC intervals, 
%    with interval IC associated with constant value VC(I),
%    separated by NC-1 ascending sorted breakpoints.
%
%    The function is to be evaluated at NP points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NC, the number of function values.
%
%    Input, real XC(NC-1), the breakpoints, in ascending order.
%
%    Input, real VC(NC), the function values over each interval.
%
%    Input, integer NP, the number of evaluation points.
%
%    Input, real XP(NP), the evaluation points.
%
%    Output, real VP(NP), the function value at the evaluation points.
%
  vp = zeros ( np, 1 );

  for ip = 1 : np
    kc = 1;
    for ic = 1 : nc - 1
      if ( xp(ip) < xc(ic) )
        break
      end
      kc = kc + 1;
    end
    vp(ip) = vc(kc);
  end

  return
end
function value = f ( x )

%*****************************************************************************80
%
%% F evaluates the right hand side function.
%
%  Discussion:
%
%    We are assuming the xn_min = 0.0, xn_max = 1.0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real VALUE, the value of F(X).
%
  value = 1000.0 * x * ( 1.0 - x );

  return
end
function vp = fem_value ( nn, xn, un, np, xp )

%*****************************************************************************80
%
%% fem_value evaluates a piecewise linear finite element function.
%
%  Discussion:
%
%    The finite element function is a piecewise linear function
%    with grid points stored in xn and values in un.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%
%    Input, integer NP, the number of evaluation points.
%
%    Input, real XP(NP), the evaluation coordinates.
%
%    Output, real VP(NP), the evaluation values.
%
  vp = zeros ( np, 1 );

  for i = 1 : np
    j = bracketer ( nn, xn, xp(i) );
    vp(i) = ( un(j)   * ( xn(j+1) - xp(i)         ) ...
            + un(j+1) * (           xp(i) - xn(j) ) ) ...
            /           ( xn(j+1)         - xn(j) );
  end

  return
end
function vs = heat1d ( nc, xc, vc, nn, xn, ns, xs )

%*****************************************************************************80
%
%% heat1d solves a steady 1d heat problem with piecewise constant diffusivity.
%
%  Discussion:
%
%    After meeting on 26 March, final changes to this program were:
%    * Zero Dirichlet boundary conditions;
%    * Source term f(x) = 1000 * x (1-x);
%    * Make number and value of diffusivity coefficients an input;
%    * Make number and location of sensors an input, and value an output.
%
%    After these changes, move on to 2D steady code.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NC, the number of values used to define diffusivity.
%
%    Input, real XC(NC-1), the breakpoints to define diffusivity.
%
%    Input, real VC(NC), the diffusivity values.
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), the node coordinates.
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the sensor values.
%
%  Local parameters:
%
%    Local, real UA, UB, the prescribed value of U at the first and last nodes.
%
%    Local, real UN(NN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%
%    Local, real XN(NN), the node coordinates.
%

%
%  Set up the problem.
%
  [ ua, ub ] = preprocess ( nc, xc, vc, nn, xn );
%
%  Compute the solution.
%
  un = process ( nn, xn, nc, xc, vc, ua, ub );
%
%  Report the results.
%
  vs = postprocess ( nn, xn, un, nc, xc, vc, ns, xs );

  return
end
function vs = postprocess ( nn, xn, un, nc, xc, vc, ns, xs )

%*****************************************************************************80
%
%% postprocess reports the results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%
%    Input, integer NC, the number of values used to define diffusivity.
%
%    Input, real XC(NC-1), the breakpoints to define diffusivity.
%
%    Input, real VC(NC), the values of diffusivity.
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the sensor values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'postprocess:\n' );
  fprintf ( 1, '  Report the results:\n' );
%
%  Print the solution.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    XN         UN\n' );
  fprintf ( 1, '\n' );

  for i = 1 : nn
    fprintf ( 1, '  %4d  %8f  %8f\n', i, xn(i), un(i) );
  end
%
%  Plot the solution.
%
  un_max = max ( un );
  un_max = 1.05 * un_max;

  figure ( 2 )
  clf
  hold on
  plot ( xn, un, 'b-', 'Linewidth', 3 );
  plot ( [ xc(1), xc(1) ], [ 0, un_max ], 'r:', 'linewidth', 3 );
  plot ( [ xc(2), xc(2) ], [ 0, un_max ], 'r:', 'linewidth', 3 );
  hold off
  axis ( [ xn(1), xn(nn), 0.0, un_max ] )
  grid ( 'on' );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- U(X) -->' );
  datestring = datestr ( datevec ( now ), 0 );
  title ( {'Solution of steady heat equation', datestring }, 'fontsize', 16 );
  filename = 'solution.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  pause ( 5 );
%
%  Compute and print sensor readings.
%
  vs = fem_value ( nn, xn, un, ns, xs );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sensor readings:\n' );
  fprintf ( 1, '     I    XS         VS\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ns
    fprintf ( 1, '  %4d  %8f  %8f\n', i, xs(i), vs(i) );
  end

  return
end
function [ ua, ub ] = preprocess ( nc, xc, vc, nn, xn )

%*****************************************************************************80
%
%% preprocess sets up the problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NC, the number of values used to define diffusivity.
%
%    Input, real XC(NC-1), the breakpoints to define diffusivity.
%
%    Input, real VC(NC), the values of diffusivity.
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), the node coordinates.
%
%    Output, real UA, UB, the prescribed value of U at the first and last nodes.
%

%
%  Geometry definitions.
%
  ua = 0.0;
  ub = 0.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'preprocess:\n' );
  fprintf ( 1, '  Number of nodes NN = %d\n', nn );
  fprintf ( 1, '  Leftmost node = %f\n', xn(1) );
  fprintf ( 1, '  Rightmost node = %f\n', xn(nn) );
  fprintf ( 1, '  Prescribed U(A) = %f\n', ua );
  fprintf ( 1, '  Prescribed U(B) = %f\n', ub );
%
%  Report the piecewise constant diffusivity:
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of diffusivity values NC = %d\n', nc );
  fprintf ( 1, '\n' );
  for i = 1 : nc
    fprintf ( 1, '  Diffusivity = %g\n', vc(i) );
    if ( i < nc )
      fprintf ( 1, '    until X = %g, then\n', xc(i) );
    end
  end
%
%  Plot the diffusivity.
%
  vc_max = max ( vc );
  vc_max = 1.05 * vc_max;

  figure ( 1 )
  clf
  hold on
  plot ( [ xn(1), xc(1), xc(1), xc(2), xc(2), xn(nn) ], ...
         [ vc(1),  vc(1), vc(2), vc(2), vc(3), vc(3)], 'b-', 'Linewidth', 3 );
  plot ( [ xc(1), xc(1) ], [ 0, vc_max ], 'r:', 'linewidth', 3 );
  plot ( [ xc(2), xc(2) ], [ 0, vc_max ], 'r:', 'linewidth', 3 );
  hold off
  axis ( [ xn(1), xn(nn), 0.0, vc_max ] )
  grid ( 'on' );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- rho(X) -->' );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { 'Diffusivity'; datestring }, 'fontsize', 16 );
  filename = 'diffusivity.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  pause ( 5 );

  return
end
function un = process ( nn, xn, nc, xc, vc, ua, ub )

%*****************************************************************************80
%
%% process solves the steady 1D heat equation with finite elements.
%
%  Discussion:
%
%    The program uses the finite element method, with piecewise linear basis
%    functions to solve the steady state heat equation in one dimension.
%
%    The problem is defined on the region x_min <= x <= x_max.
%
%    The following differential equation is imposed between A and B:
%
%      - d/dx k(x) du/dx = f(x)
%
%    where k(x) and f(x) are given functions.
%
%    At the boundaries, the following conditions are applied:
%
%      u(x_min) = UA
%      u(x_max) = UB
%
%    A set of NN equally spaced nodes XN is defined on this
%    interval, with x_min = XN(1) < XN(2) < ... < XN(NN) = x_max.
%
%    At each node I, we associate a piecewise linear basis function V(I,X),
%    which is 0 at all nodes except node I.  This implies that V(I,X) is
%    everywhere 0 except that
%
%    for XN(I-1) <= X <= XN(I):
%
%      V(I,X) = ( X - XN(I-1) ) / ( XN(I) - XN(I-1) ) 
%
%    for XN(I) <= X <= XN(I+1):
%
%      V(I,X) = ( XN(I+1) - X ) / ( XN(I+1) - XN(I) )
%
%    We now assume that the solution U(X) can be written as a linear
%    sum of these basis functions:
%
%      U(X) = sum ( 1 <= J <= NN ) UN(J) * V(J,X)
%
%    where U(X) is the function of X, and UN(J)
%    is a coefficient of a basis functions.
%
%    To determine the coefficients UN, we multiply the original
%    differential equation by the basis function V(J,X), and use
%    integration by parts, to arrive at the I-th finite element equation:
%
%        Integral K(X) * U'(X) * V'(I,X) dx = Integral F(X) * V(I,X) dx
%
%    We note that the functions U(X) and U'(X) can be replaced by
%    the finite element form involving the linear sum of basis functions,
%    but we also note that the resulting integrand will only be nonzero
%    for terms where J = I - 1, I, or I + 1.
%
%    By writing this equation for basis functions I = 2 through NN - 1,
%    and using the boundary conditions, we have NN linear equations
%    for the NN unknown coefficients UN(1) through UN(NN).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), the node coordinates.
%
%    Input, integer NC, the number of values used to define diffusivity.
%
%    Input, real XC(NC-1), the breakpoints to define diffusivity.
%
%    Input, real VC(NC), the values of diffusivity.
%
%    Input, real UA, UB, the prescribed value of U at the first and last nodes.
%
%    Output, real UN(NN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%

%
%  Quadrature definitions.
%
  nq = 2;
  xq(1) = -0.577350269189625764509148780502;
  xq(2) = +0.577350269189625764509148780502;
  wq(1) = 1.0;
  wq(2) = 1.0;
%
%  Make room for the matrix and right hand side.
%
  nzmax = 3 * nn;
  Amat = sparse ( [], [], [], nn, nn, nzmax );
  bvec = zeros ( nn, 1 );
%
%  Equation 1 is the left boundary condition, U(x_min) = UA;
%
  Amat(1,1) = 1.0;
  bvec(1) = ua;
%
%  Equation I involves the basis function at node I.
%  This basis function is nonzero from X(I-1) to X(I+1).
%  Equation I looks like this:
%
%    Integral A(X) U'(X) V'(I,X) dx = Integral F(X) V(I,X) dx
%
%  Then, we realize that U(X) = sum ( 1 <= J <= NN ) U(J) * V(J,X), 
%  (U(X) means the function; U(J) is the coefficient of V(J,X) ).
%
%  The only V functions that are nonzero when V(I,X) is nonzero are
%  V(I-1,X) and V(I+1,X). 
%
%  Let's use the shorthand 
%
%    VL(X) = V(I-1,X)
%    VM(X) = V(I,X)
%    VR(X) = V(I+1,X)
%
%  So our equation becomes
%
%    Integral A(X) [ VL'(X) U(I-1) + VM'(X) U(I) + VR'(X) U(I+1) ] * VM'(X) dx
%  = Integral F(X) VM(X) dx.
%
%  
%
%  This is actually a set of NN-2 linear equations for the NN coefficients UN.
%
%  Now gather the multipliers of U(I-1) to get the matrix entry A(I,I-1), and so on.
%
  for i = 2 : nn - 1
%
%  Get the left, right and middle coordinates.
%
    xl = xn(i-1);
    xm = xn(i);
    xr = xn(i+1);
%
%  Make temporary variables for A(I,I-1), A(I,I), A(I,I+1) and B(I).
%
    al = 0.0;
    am = 0.0;
    ar = 0.0;
    bm = 0.0;
%
%  We approximate the integrals by using a weighted sum of
%  the integrand values at quadrature points.
%
    for q = 1 : nq
%
%  Integrate over the LEFT interval, between XL and XM, where:
%
%  VL(X) = ( XM - X       ) / ( XM - XL )
%  VM(X) = (      X  - XL ) / ( XM - XL )
%  VR(X) = 0
%
%  VL'(X) =             - 1 / ( XM - XL )
%  VM'(X) =             + 1 / ( XM - XL ) 
%  VR'(X) = 0
%
      xql = ( ( 1.0 - xq(q) ) * xl   ...
            + ( 1.0 + xq(q) ) * xm ) ...
            / 2.0;

      wql = wq(q) * ( xm - xl ) / 2.0;

      vl = ( xm - xql ) / ( xm - xl ); 
      vlp =     - 1.0   / ( xm - xl );

      vm = ( xql - xl ) / ( xm - xl );
      vmp = +1.0        / ( xm - xl );

      vr =  0.0;
      vrp = 0.0;

      kxql = diffusivity_pwc_1d ( nc, xc, vc, 1, xql );
      fxql = f ( xql );

      al = al + wql * ( kxql * vlp * vmp );
      am = am + wql * ( kxql * vmp * vmp );
      ar = ar + wql * ( kxql * vrp * vmp );
      bm = bm + wql * ( fxql * vm );
%
%  Integrate over the RIGHT interval, between XM and XR, where:
%
%  VL(X) = 0
%  VM(X) = ( XR - X       ) / ( XR - XM )
%  VR(X) = (      X  - XM ) / ( XR - XM )
%
%  VL'(X) = 0
%  VM'(X) =             - 1 / ( XR - XM )
%  VR'(X) =             + 1 / ( XR - XM ) 
%
      xqr = ( ( 1.0 - xq(q) ) * xm   ...
            + ( 1.0 + xq(q) ) * xr ) ...
            / 2.0;

      wqr = wq(q) * ( xr - xm ) / 2.0;

      vl = 0.0;
      vlp = 0.0;

      vm = ( xr - xqr ) / ( xr - xm );
      vmp =     - 1.0   / ( xr - xm );

      vr = ( xqr - xm ) / ( xr - xm );
      vrp =  1.0        / ( xr - xm );

      kxqr = diffusivity_pwc_1d ( nc, xc, vc, 1, xqr );
      fxqr = f ( xqr );

      al = al + wqr * ( kxqr * vlp * vmp );
      am = am + wqr * ( kxqr * vmp * vmp );
      ar = ar + wqr * ( kxqr * vrp * vmp );
      bm = bm + wqr * ( fxqr * vm );

    end

    Amat(i,i-1) = al;
    Amat(i,i)   = am;
    Amat(i,i+1) = ar;
    bvec(i)     = bm;

  end
%
%  Equation NN is the right boundary condition, U(x_max) = UB;
%
  Amat(nn,nn) = 1.0;
  bvec(nn) = ub;
%
%  Solve the linear system.
%
  un = Amat \ bvec;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'process:\n' );
  fprintf ( 1, '  The solution un has been computed.\n' );

  return
end
function [ xmin, xmax ] = region_box ( )

%*****************************************************************************80
%
%% region_box returns the minimum and maximum x values of the region.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real XMIN, XMAX, the argument ranges.
%
  xmin = 0.0;
  xmax = 1.0;

  return
end
