function fem1d_approximate ( mesh_file, data_file, weight_a, weight_d, weight_b )

%*****************************************************************************80
%
%% fem1d_approximate approximates data with a finite element function.
%
%  Discussion:
%
%    Given a pre-existing finite element mesh (node_x), and a set of data pairs
%    (data_x,data_y), we seek finite element coefficients (node_c), which
%    define a finite element function fem(x), so that we minimize the 
%    approximation error:
%
%      e^2 = sum ( fem(data_x) - data_y )^2 )
%
%    But doing only this, the finite element function won't satisfy boundary
%    conditions, and will be zero in intervals where no data is supplied.
%
%    We add a derivative condition, that the second derivative be zero,
%    and a boundary condition.
%
%    By multiplying the approximation, derivative, and boundary conditions
%    by appropriate weights, we can manipulate the finite element function
%    so that it does a reasonable job of approximating the data, and
%    at the same time tries to minimize the second derivative - which means
%    the function doesn't flatline to zero in data-free intervals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string mesh_file, the file containing the mesh node coordinates.
%
%    Input, string data_file, the file containing data coordinates and values.
%
%    Input, real weight_a, weight_d, weight_b, the weights for approximation,
%    derivative minimization, and boundary condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fem1d_approximate\n' );
  fprintf ( 1, '  MATLAb version\n' );
  fprintf ( 1, '  Given pairs of data,\n' );
  fprintf ( 1, '  construct a finite element approximating function on a\n' );
  fprintf ( 1, '  pre-existing mesh.\n' );
%
%  Read the mesh.
%
  node_x = load ( mesh_file );
  node_x = node_x(:);
  node_num = length ( node_x );
  r8vec_print ( node_num, node_x, '  Mesh node coordinates:' );
%
%  Read the data.
%
  data = load ( data_file );
  data_x = data(:,1);
  data_y = data(:,2);
  data_num = length ( data_x );
  r8vec2_print ( data_num, data_x, data_y, '  Data (x,y) pairs:' );
%
%  Bracket the data.
%
  [ data_l, data_r ] = data_bracket ( node_num, node_x, data_num, data_x );
  is_legal = ( node_x(1) <= data_x & data_x <= node_x(node_num) );
  legal_index = find ( ( node_x(1) <= data_x & data_x <= node_x(node_num) ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Bracketing:\n' );
  for i = 1 : data_num
    fprintf ( 1, '\n' );
    if ( ~ is_legal(i) )
      fprintf ( 1, '               data_x(%d) = %6.2f  Outside the interval\n', i, data_x(i) );
    else
      fprintf ( 1, '  node_x(%d) <= data_x(%d) <= node_x(%d)\n', data_l(i), i, data_r(i) );
      fprintf ( 1, '  %f  <= %f <= %f\n', node_x(data_l(i)), data_x(i), node_x(data_r(i)) );
    end
  end
%
%  Count the equations.
%
  eq_num = sum ( data_l ~= -1 ) + node_num;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  There will be %d conditions:\n', eq_num );
  fprintf ( 1, '  Approximation conditions weighted by %g\n', weight_a );
  fprintf ( 1, '  Derivative conditions weighted by    %g\n', weight_d );
  fprintf ( 1, '  Boundary conditions weighted by      %g\n', weight_b );
%
%  Write the error equations.
%
  A = zeros ( eq_num, node_num );
  b = zeros ( eq_num, 1 );

  eq_i = 0;
%
%  Approximation conditions.
%
  for i = 1 : data_num
    if ( is_legal(i) )
      eq_i = eq_i + 1;
      fem_l = data_l(i);
      fem_r = data_r(i);
%
%  fem = fem_c(fem_l);
%
      if ( data_l(i) == data_r(i) )
        A(eq_i,fem_l) = weight_a * 1.0;
%
%  fem = ( ( node_x(fem_r) - data_x(i)                 ) * node_c(fem_l) 
%        + (                 data_x(i) - node_x(fem_l) ) * node_c(fem_r) ) 
%        / ( node_x(fem_r)             - node_x(fem_l) );
%
      else
        A(eq_i,fem_l) = weight_a * ( node_x(fem_r) - data_x(i) ) ...
                                / ( node_x(fem_r) - node_x(fem_l) );
        A(eq_i,fem_r) = weight_a *  ( data_x(i)     - node_x(fem_l) ) ...
                                / ( node_x(fem_r) - node_x(fem_l) );
      end
      b(eq_i) = weight_a * data_y(i);
    end
  end
%
%  Derivative conditions.
%
  for i = 2 : node_num - 1
    eq_i = eq_i + 1;
    h_left = node_x(i) - node_x(i-1);
    h_right = node_x(i+1) - node_x(i);
    A(eq_i,i-1) = weight_d / h_left;
    A(eq_i,i) = - weight_d / h_right - weight_d / h_left ;
    A(eq_i,i+1) = weight_d / h_right;
  end
%
%  Boundary conditions.
%
  eq_i = eq_i + 1;
  A(eq_i,1) = weight_b * 1.0;
  b(eq_i) = weight_b * 0.0;

  eq_i = eq_i + 1;
  A(eq_i,node_num) = weight_b * 1.0;
  b(eq_i) = weight_b * 0.0;
%
%  For debugging, you might want to see A, b, node_c.
%
  node_c = A \ b;

  if ( false )
    A
    b
    node_c
  end
%
%  Plot FEM and data.
%
  clf ( );
  hold ( 'on' );
  plot ( node_x, node_c, 'b-', 'linewidth', 3 );
  plot ( node_x, node_c, 'ko', 'markersize', 15 );
  plot ( data_x(legal_index), data_y(legal_index), 'r.', 'markersize', 30 );
  grid ( 'on' );
  xlabel ( '<--- X --->', 'fontsize', 16 );
  ylabel ( '<--- Y --->', 'fontsize', 16 );
  label = s_escape_tex ( data_file );
  title ( label, 'fontsize', 16 );
  hold ( 'off' );

  filename = 'fem1d_approximate_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"', filename );
  pause ( 5 );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fem1d_approximate\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end

