function [ ns, xs, ys, vs ] = paraheat_gaussian ( xc, yc, sc, vc )

%*****************************************************************************80
%
%% paraheat_gaussian demonstrates the parameterized 2d heat equation solver.
%
%  Discussion:
%
%    The diffusivity parameters xc, yc, sc, vc are considered the "input", 
%    and the sensor readings ns, xs, ys, vs are considered the output.  
%
%    The values ns, xs, ys are usually fixed, but vs can be considered
%    a function of the input parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Cohen, Wolfgang Dahmen, Ron DeVore, Jalal Fadili, Olga Mula,
%    James Nichols,
%    Optimal reduced model algorithms for data-based state estimation,
%    ArXiv:1903.07938v1 [math.NA] 19 Mar 2019
%
%  Parameters:
%
%    Input, real XC, YC, SC, VC, the gaussian diffusivity parameters.
%
%    Output, integer NS, the number of sensors.
%
%    Output, real XS(NS), YS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the sensor readings.
%
%  Local parameters:
%
%    Local, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
%    Local, integer E_NODES(3,NE), the nodes forming each element.
%
%    Local, integer NE, the number of elements.
%
%    Local, integer NN, the number of nodes.
%
%    Local, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Local, integer SENSOR_LOCATION_SEED, a seed for the random number 
%    generator, used if sensor locations will be randomly selected.
%
%    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), YN(NN), the node coordinates.
%
  global xc_save
  global yc_save
  global sc_save
  global vc_save

  sensor_location_seed = 987654321;

  request_comments = false;
  request_diffusivity_plot = false;
  request_diffusivity_print = false;
  request_element_plot = false;
  request_element_report = false;
  request_fem_contour_plot = false;
  request_fem_report = false;
  request_fem_surface_plot = false;
  request_grid_plot = false;
  request_grid_report = false;
  request_rhs_report = false;
  request_sensor_location_plot = false;
  request_sensor_location_report = false;
  request_sensor_value_report = false;

  if ( request_comments )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'paraheat_gaussian\n' );
    fprintf ( 1, '  MATLAB version\n' );
    fprintf ( 1, '  Parameterized heat equation in 2D.\n' );
    fprintf ( 1, '  Steady heat problem with gaussian diffusivity.\n' );
    fprintf ( 1, '  Diffusivity parameters are input.\n' );
    fprintf ( 1, '  Sensor readings are output.\n' );
  end
%
%  Copy the input into the global data.
%
  xc_save = xc;
  yc_save = yc;
  sc_save = sc;
  vc_save = vc;
%
%  Close all previous figures.
%
  close ( 'all' );
%
%  Diffusivity.
%
  if ( request_diffusivity_print )
    diffusivity_gaussian_print ( );
  end

  if ( request_diffusivity_plot )
    diffusivity_gaussian_contour ( );
    diffusivity_gaussian_surface ( );
  end
%
%  Grid.
%
  [ nxn, nyn, nn, xn, yn, bn ] = grid_define ( );

  if ( request_grid_report )
    grid_report ( nxn, nyn, nn, xn, yn, bn );
  end

  if ( request_grid_plot )
    grid_plot ( nn, xn, yn, bn );
  end
%
%  Elements.
%
  [ ne, e_nodes ] = element_define ( nxn, nyn );

  if ( request_element_report )
    element_report ( ne, e_nodes );
  end

  if ( request_element_plot )
    element_plot ( ne, e_nodes, nn, xn, yn );
  end
%
%  FEM.
%
  if ( request_rhs_report )
    rhs_report ( );
  end

  un = fem_define ( nn, nxn, nyn, xn, yn, bn, ne, e_nodes );

  if ( request_fem_report )
    fem_report ( nxn, nyn, nn, xn, yn, bn, un, nonzeros );
  end

  if ( request_fem_contour_plot )
    fem_contour_plot ( nxn, nyn, nn, xn, yn, un );
  end

  if ( request_fem_surface_plot )
    fem_surface_plot ( nn, xn, yn, un, ne, e_nodes );
  end
%
%  Sensor locations.
%
  [ ns, xs, ys ] = sensor_location_define ( sensor_location_seed );

  if ( request_sensor_location_report )
    sensor_location_report ( ns, xs, ys, sensor_location_seed );
  end

  if ( request_sensor_location_plot )
    sensor_location_plot ( ns, xs, ys, nn, xn, yn, bn, sensor_location_seed );
  end
%
%  Sensor values.
%
  if ( false )
    vs = sensor_value_define ( nn, xn, yn, un, ne, e_nodes, ns, xs, ys );
  else
    vs = sensor_value_average_define ( nn, xn, yn, un, ne, e_nodes, ns, xs, ys );
  end

  if ( request_sensor_value_report )
    sensor_value_report ( ns, xs, ys, vs );
  end
%
%  Terminate.
%
  if ( request_comments )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'paraheat_gaussian' );
    fprintf ( 1, '  Normal end of execution.\n' );
  end

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

%*****************************************************************************80
%
%% bracket_left finds the left endpoint of a bracketing interval.
%
%  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, 'bracket_left - Fatal error!\n' );
    fprintf ( 1, '  N must be at least 2.\n' );
    error ( 'bracket_left - 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 diffusivity_gaussian_contour ( )

%*****************************************************************************80
%
%% diffusivity_gaussian_contour plots the gaussian diffusivity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2019
%
%  Author:
%
%    John Burkardt
%
  rng_status = rng ( );
  diffusivity_seed = rng_status.Seed;
  filename = sprintf ( 'diffusivity_gaussian_contour_%d.png', diffusivity_seed );

  figure ( 'Name', 'Diffusivity Contour Plot' );
  clf ( );

  nxp = 101;
  nyp = 101;
  [ xmin, xmax, ymin, ymax ] = region_box ( );
  xp = linspace ( xmin, xmax, nxp );
  yp = linspace ( ymin, ymax, nyp );
  [ Xmat, Ymat ] = ndgrid ( xp, yp );
  np = nxp * nyp;

  Amat = diffusivity_gaussian_func ( np, Xmat, Ymat );
  Amat = reshape ( Amat, nxp, nyp );
%
%  Plot the diffusivity field.
%
  [ M, c ] = contour ( Xmat, Ymat, Amat, 20 );
  c.LineWidth = 3;
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );
  pause ( 5 );
  close ( );

  return
end
function [ xc, yc, sc, vc ] = diffusivity_gaussian_define ( )

%*****************************************************************************80
%
%% diffusivity_gaussian_define defines the gaussian diffusivity.
%
%  Discussion:
%
%    The gaussian diffusivity is modeled by a gaussian function centered
%    at (xc,yc) with standard deviation sc and maximum value vc.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real XC, YC, the center of the Gaussian.
%
%    Output, real SC, the standard deviation of the Gaussian.
%
%    Output, real VC, the maximum value of the Gaussian.
%
  global xc_save
  global yc_save
  global sc_save
  global vc_save

  xc = xc_save;
  yc = yc_save;
  sc = sc_save;
  vc = vc_save;

  return
end
function vp = diffusivity_gaussian_func ( np, xp, yp )

%*****************************************************************************80
%
%% diffusivity_gaussian_func evaluates the gaussian diffusivity function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real NP, the number of evaluation points.
%
%    Input, real XP(NP), Y(NP), the evaluation points.
%
%    Output, real VP(NP), the values of F.
%
  [ xc, yc, sc, vc ] = diffusivity_gaussian_define ( );

  t = ( xp(1:np) - xc ) .^ 2 + ( yp(1:np) - yc ) .^ 2;
  vp = zeros(np,1);
  vp(1:np) = 0.5 + vc * exp ( - 0.5 * t / sc^2 );

  return
end
function diffusivity_gaussian_print ( xc, yc, sc, vc )

%*****************************************************************************80
%
%% diffusivity_gaussian_print: information about the diffusivity.
%
%  Discussion:
%
%    The gaussian test diffusivity is modeled by a gaussian function centered
%    at (xc,yc) with standard deviation sc and maximum value vc.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XC, YC, the center of the Gaussian.
%
%    Input, real SC, the standard deviation of the Gaussian.
%
%    Input, real VC, the maximum value of the Gaussian.
%
  [ xc, yc, sc, vc ] = diffusivity_gaussian_define ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'diffusivity_gaussian_print:\n' );
  fprintf ( 1, '  Diffusivity is a gaussian function:\n' );
  fprintf ( 1, '  d(x,y) = 0.5 + vc * exp ( - 0.5 * ((x-xc)^2+(y-yc)^2)/sc^2 )\n' );
  fprintf ( 1, '  (xc,yc) = {%g,%g)\n', xc, yc );
  fprintf ( 1, '  sc = %g, vc = %g\n', sc, vc );

  return
end
function diffusivity_gaussian_surface ( )

%*****************************************************************************80
%
%% diffusivity_gaussian_surface plots the diffusivity.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2019
%
%  Author:
%
%    John Burkardt
%
  rng_status = rng ( );
  diffusivity_seed = rng_status.Seed;
  filename = sprintf ( 'diffusivity_gaussian_surface_%d.png', diffusivity_seed );

  figure ( 'Name', 'Diffusivity Surface Plot' );
  clf ( );

  nxp = 101;
  nyp = 101;
  [ xmin, xmax, ymin, ymax ] = region_box ( );
  xp = linspace ( xmin, xmax, nxp );
  yp = linspace ( ymin, ymax, nyp );
  [ Xmat, Ymat ] = ndgrid ( xp, yp );
  np = nxp * nyp;

  Amat = diffusivity_gaussian_func ( np, Xmat, Ymat );
  Amat = reshape ( Amat, nxp, nyp );
%
%  Plot the diffusivity field.
%
  surf ( Xmat, Ymat, Amat, 'EdgeColor', 'interp' );
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );
  pause ( 5 );
  close ( );

  return
end
function [ ne, e_nodes ] = element_define ( nxn, nyn )

%*****************************************************************************80
%
%% element_define defines a grid of pairs of 3 node triangles.
%
%  Discussion:
%
%    A rectangular grid is assumed, which can be thought of as a
%    grid of (nelemx+1) x (nelemy+1) evenly spaced points, or as
%    a grid of nelemx x nelemy squares, each square being subdivided
%    into a pair of triangles.
%
%    Node numbering starts at the lower left, moves right first, and
%    then upwards.
%
%    Element numbering goes similarly, square by square, with the
%    lower triangle of each pair numbered first.
%
%    The element_node array records, for each element, the indices of
%    the three nodes that form the vertices of the triangle, in counter
%    clockwise order, starting with the node associated with the 90 degree
%    angle.
%
%  Example:
%
%    Input:
%
%      NXN = 4, NYN = 3
%
%    Output:
%
%      e_nodes =
%         1,  2,  5;
%         6,  5,  2;
%         2,  3,  6;
%         7,  6,  3;
%         3,  4,  7;
%         8,  7,  4;
%         5,  6,  9;
%        10,  9,  6;
%         6,  7, 10;
%        11, 10,  7;
%         7,  8, 11;
%        12, 11,  8.
%
%  Grid:
%
%    9---10---11---12     +----+----+----+
%    |\   |\   |\   |     |\ 8 |\10 |\12 |
%    | \  | \  | \  |     | \  | \  | \  |
%    |  \ |  \ |  \ |     |  \ |  \ |  \ |
%    |   \|   \|   \|     |  7\|  9\| 11\|
%    5----6----7----8     +----+----+----+
%    |\   |\   |\   |     |\ 2 |\ 4 |\ 6 |
%    | \  | \  | \  |     | \  | \  | \  |
%    |  \ |  \ |  \ |     |  \ |  \ |  \ |
%    |   \|   \|   \|     |  1\|  3\|  5\|
%    1----2----3----4     +----+----+----+
%      Node indexing      Element indexing
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer nxn, nyn, the number of nodes in the X and Y directions.
%
%    Output, integer ne, the number of elements.
%
%    Output, integer e_nodes(3,ne), the nodes that form each element.
%
  ne = 2 * ( nxn - 1 ) * ( nyn - 1 );
  e_nodes = zeros ( 3, ne );

  e = 0;

  for j = 1 : nyn - 1
    for i = 1 : nxn - 1

      swn = i     + ( j - 1 ) * nxn;
      sen = i + 1 + ( j - 1 ) * nxn;
      nwn = i     +   j       * nxn;
      nen = i + 1 +   j       * nxn;

      e = e + 1;
      e_nodes(1,e) = swn;
      e_nodes(2,e) = sen;
      e_nodes(3,e) = nwn;

      e = e + 1;
      e_nodes(1,e) = nen;
      e_nodes(2,e) = nwn;
      e_nodes(3,e) = sen;

    end
  end

  return
end
function element_plot ( ne, e_nodes, nn, xn, yn )

%*****************************************************************************80
%
%% element_plot plots the elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ne, the number of elements.
%
%    Input, integer e_nodes(3,ne), the nodes that form each element.
%
%    Input, integer nn, the number of nodes.
%
%    Input, real xn(nn), yn(nn), the coordinates of the nodes.
%
  figure ( 'Name', 'FEM Elements' );
  clf ( );

  filename = 'elements.png';
%
%  Get the scale.
%
  [ xmin, xmax, ymin, ymax ] = region_box ( );
  margin = 0.05 * max ( xmax - xmin, ymax - ymin );
%
%  Draw the element lines.
%
  hold ( 'on' );

  for e = 1 : ne
    p = [ e_nodes(1:3,e); e_nodes(1,e) ];
    line ( 'XData', xn(p), 'YData', yn(p), 'Color', 'r', 'LineWidth', 2 );
  end
%
%  Label the elements
%
  for e = 1 : ne

    ave_x = 0.0;
    ave_y = 0.0;

    for i = 1 : 3
      n = e_nodes(i,e);
      ave_x = ave_x + xn(n);
      ave_y = ave_y + yn(n);
    end

    ave_x = ave_x / 3;
    ave_y = ave_y / 3;

    text ( ave_x, ave_y, num2str ( e ) );

  end
%
%  Label the plot.
%
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  axis ( [ xmin, xmax, ymin, ymax ] );
  axis ( 'equal' );
  hold ( 'off' );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );

  return
end
function element_report ( ne, e_nodes )

%*****************************************************************************80
%
%% element_report reports the elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ne, the number of elements.
%
%    Input, integer e_nodes(3,ne), the nodes that form each element.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'element_report:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of elements = %d\n', ne );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index   N1  N2  N3\n' );
  fprintf ( 1, '\n' );
  for j = 1 : ne
    fprintf ( 1, '  %3d:  %3d  %3d  %3d\n', j, e_nodes(1:3,j) );
  end

  return
end
function [ un, nonzeros ] = fem_define ( nn, nxn, nyn, xn, yn, bn, ne, e_nodes )

%*****************************************************************************80
%
%% fem_define solves the steady 2D heat equation with finite elements.
%
%  Discussion:
%
%    Grid: specified by nn, nxn, nyn, xn, yn;
%    Elements: specified by ne, e_nodes;
%    Boundary conditions: bn indicates location of zero Dirichlet nodes;
%    Diffusivity: user sets choice in diffusivity_gaussian_func();
%    RHS: user sets choice in rhs_func();
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
%    Input, integer NE, the number of elements.
%
%    Input, integer E_NODES(3,NE), the nodes forming each element.
%
%    Output, real UN(NN), the finite element coefficients, which are also
%    the value of the computed solution at the mesh points.
%
%    Output, integer NONZEROS, the number of nonzero entries in the sparse
%    FEM matrix.
%
%  Local Parameters:
%
%    Local, integer NQ, the number of quadrature points.
%
%    Local, real WQ(NQ), the quadrature weights.
%
%    Local, real XQ(NQ), YQ(NQ), the quadrature point locations.
%
 
%
%  Set the quadrature.
%
  nq = 7;
  [ wq, xq, yq ] = quad_define ( nq );
%
%  Allocate linear system varriables.
%
  nzmax = 8 * nn;
  Amat = sparse ( [], [], [], nn, nn, nzmax );
  bvec = zeros ( nn, 1 );
  uvec = zeros ( nn, 1 );
%
%  Consider element E.
%
  for e = 1 : ne
%
%  Extract element E vertex coordinates.
%
    xt3(1:3) = xn(e_nodes(1:3,e));
    yt3(1:3) = yn(e_nodes(1:3,e));
%
%  Map the quadrature rule to element E.
%
    np = nq;
    [ xp, yp ] = reference_to_physical_t3 ( xt3, yt3, nq, xq, yq );
    area = triangle_area ( xt3, yt3 );
    wp(1:nq) = area * wq(1:nq);
%
%  Evaluate diffusivity and rhs at all quadrature points.
%
    lq = diffusivity_gaussian_func ( np, xp, yp );
    fq = rhs_func ( np, xp, yp );
%
%  Consider the Q-th quadrature point.
%
    for q = 1 : nq
%
%  Consider the TEST-th test function.
%
%  We generate an integral for every node associated with an unknown.
%
      for test = 1 : 3

        tp1 = 1 + mod ( test,     3 );
        tp2 = 1 + mod ( test + 1, 3 );
        bi = ( ( xt3(tp2) - xt3(tp1) ) * ( yp(q) - yt3(tp1) ) ...
             - ( yt3(tp2) - yt3(tp1) ) * ( xp(q) - xt3(tp1) ) ) / area;

        i = e_nodes(test,e);
        bvec(i) = bvec(i) + wp(q) * fq(q) * bi;
        dbidx = - ( yt3(tp2) - yt3(tp1) ) / area;
        dbidy =   ( xt3(tp2) - xt3(tp1) ) / area;
%
%  Consider the BASIS-th basis function.
%
        for basis = 1 : 3

          tp1 = 1 + mod ( basis,     3 );
          tp2 = 1 + mod ( basis + 1, 3 );
          bj = ( ( xt3(tp2) - xt3(tp1) ) * ( yp(q) - yt3(tp1) ) ...
               - ( yt3(tp2) - yt3(tp1) ) * ( xp(q) - xt3(tp1) ) ) / area;
          dbjdx = - ( yt3(tp2) - yt3(tp1) ) / area;
          dbjdy =   ( xt3(tp2) - xt3(tp1) ) / area;

          j = e_nodes(basis,e);
          Amat(i,j) = Amat(i,j) ...
            + wp(q) * ( lq(q) * ( dbidx * dbjdx + dbidy * dbjdy ) ...
            );

        end

      end

    end

  end
%
%  Handle the zero Dirichlet boundary conditions.
%
  i = find ( bn == 1 );
  Amat(i,:) = 0;
  Amat(i+(i-1)*nn) = 1.0;
  bvec(i) = 0.0;

  nonzeros = nnz ( Amat );
%
%  Solve the linear system.
%
  un = Amat \ bvec;

  return
end
function fem_contour_plot ( nxn, nyn, nn, xn, yn, un )

%*****************************************************************************80
%
%% fem_contour_plot plots the finite element function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element function at the nodes.
%
  figure ( 'Name', 'FEM Contour Plot' );
  clf ( );

  rng_status = rng ( );
  diffusivity_seed = rng_status.Seed;

  un_max = max ( un );
  un_min = min ( un );
  un_min = min ( un_min, 0.0 );
  un_range = un_max - un_min;
  caxis ( [ un_min, un_max ] );

  filename = sprintf ( 'fem_contour_%d.png', diffusivity_seed );

  levels = linspace ( un_min, un_max, 21 );
  [ M, c ] = contour ( reshape(xn,nxn,nyn), reshape(yn,nxn,nyn), ...
    reshape(un,nxn,nyn), levels );
  c.LineWidth = 3;
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  axis ( 'equal' );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );

  return
end
function fem_surface_plot ( nn, xn, yn, un, ne, e_nodes )

%*****************************************************************************80
%
%% fem_surface_plot plots the finite element function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element function at the nodes.
%
%    Input, integer NE, the number of elements.
%
%    Input, integer E_NODES(3,NE), the nodes forming each element.
%
  figure ( 'Name', 'FEM Surface Plot' );
  clf ( );

  rng_status = rng ( );
  diffusivity_seed = rng_status.Seed;

  un_max = max ( un );
  un_min = min ( un );
  un_min = min ( un_min, 0.0 );
  un_range = un_max - un_min;
  caxis ( [ un_min, un_max ] );

  filename = sprintf ( 'fem_surface_%d.png', diffusivity_seed );

  trisurf ( e_nodes', xn, yn, un, 'Facecolor', 'interp', 'Edgecolor', 'None' );
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  daspect ( [ 1.0, 1.0, un_range ] );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );

  return
end
function fem_report ( nxn, nyn, nn, xn, yn, bn, un, nonzeros )

%*****************************************************************************80
%
%% fem_report reports the finite element function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
%    Input, real UN(NN), the finite element function at the nodes.
%
%    Input, integer NONZEROS, the number of nonzeros in the sparse FEM matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'fem_report:\n' );
  fprintf ( 1, '  Number of nodes NN = %d\n', nn );
  fprintf ( 1, '  FEM sparse matrix uses %d nonzero entries.\n', nonzeros );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    K    I    J    X    Y    B  U\n' );
  k = 0;
  for j = 1 : nyn
    fprintf ( 1, '\n' );
    for i = 1 : nxn
      k = k + 1;
      fprintf ( 1, '  %3d  %3d  %3d  %f  %f  %f  %f\n', ...
      k, i, j, xn(k), yn(k), bn(k), un(k) );
    end
  end

  return
end
function [ nxn, nyn, nn, xn, yn, bn ] = grid_define ( )

%*****************************************************************************80
%
%% grid_define defines the grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Output, integer NN, the total number of nodes.
%
%    Output, real XN(NN), YN(NN), the node coordinates.
%
%    Output, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
  nxn = 129;
  nyn = 129;

  [ xmin, xmax, ymin, ymax ] = region_box ( );
  x = linspace ( xmin, xmax, nxn );
  y = linspace ( ymin, ymax, nyn );

  nn = nxn * nyn;
  xn = zeros ( nn, 1 );
  yn = zeros ( nn, 1 );
  bn = zeros ( nn, 1 );

  k = 0;
  for j = 1 : nyn
    for i = 1 : nxn
      k = k + 1;
      xn(k) = x(i);
      yn(k) = y(j);
      if ( i == 1 || i == nxn || j == 1 || j == nyn )
        bn(k) = 1;
      else
        bn(k) = 0;
      end 
    end
  end

  return
end
function grid_plot ( nn, xn, yn, bn )

%*****************************************************************************80
%
%% grid_plot plots the grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
  figure ( 'Name', 'FEM Grid' );
  clf ( );

  filename = 'grid.png';
%
%  Get the scale.
%
  [ xmin, xmax, ymin, ymax ] = region_box ( );
  margin = 0.05 * max ( xmax - xmin, ymax - ymin );
%
%  Draw the nodes.
%
  hold ( 'on' );
  interior = find ( bn == 0 );
  boundary = find ( bn == 1 );

  plot ( xn(interior), yn(interior), 'bo', 'markersize', 5 );
  plot ( xn(boundary), yn(boundary), 'k.', 'markersize', 20 );
%
%  Label the nodes.
%
  abit = 0.05;
  for j = 1 : nn
    text ( xn(j), yn(j) + abit, num2str(j) );
  end
%
%  Label the plot.
%
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  axis ( [ xmin - margin, xmax + margin, ymin - margin, ymax + margin ] );
  axis ( 'equal' );
  hold ( 'off' );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );

  return
end
function grid_report ( nxn, nyn, nn, xn, yn, bn )

%*****************************************************************************80
%
%% grid_report reports the grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXN, NYN, the number of nodes in the X and Y directions.
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'grid_report:\n' );
  fprintf ( 1, '  X grid uses NXN = %d\n', nxn );
  fprintf ( 1, '  Y grid uses NXN = %d\n', nyn );
  fprintf ( 1, '  Number of nodes NN = NXN * NYN = %d\n', nn );
  fprintf ( 1, '  %f <= x <= %f\n', min ( xn ), max ( xn ) );
  fprintf ( 1, '  %f <= y <= %f\n', min ( yn ), max ( yn ) );
  h = min ( ( max ( xn ) - min ( xn ) ) / ( nxn - 1 ), ...
            ( max ( yn ) - min ( yn ) ) / ( nyn - 1 ) );
  fprintf ( 1, '  h = min ( hx, hy ) = %g\n', h );

  fprintf ( 1, '\n' );
  fprintf ( 1, '    K    I    J    X    Y    B\n' );
  k = 0;
  for j = 1 : nyn
    fprintf ( 1, '\n' );
    for i = 1 : nxn
      k = k + 1;
      fprintf ( 1, '  %3d  %3d  %3d  %f  %f  %f\n', ...
      k, i, j, xn(k), yn(k), bn(k) );
    end
  end

  return
end
function pwc_bars_2d ( nxc, nyc, xc, yc, vc, filename )

%*****************************************************************************80
%
%% pwc_bars_2d plots a piecewise constant function of (X,Y) as bars.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXC, NYC, the number of columns and rows of data.
%
%    Input, real XC(NXC+1), YC(NYC+1), the breakpoint coordinates.
%
%    Input, real VC(NXC,NYC), the values.
%
%    Input, string FILENAME, the filename to be used.
%
  figure ( 'Name', 'Diffusivity Bar Plot' );
  clf ( );
%
%  Determine the mapping between colors and values.
%
  zmax = max ( max ( vc ) );
  zmin = min ( min ( vc ) );
  caxis ( [ 0.0, zmax ] );
%
%  Pick the colors that will correspond to the minimum and maximum
%  values of Z.
%
  rmax = 0.8;
  gmax = 0.2;
  bmax = 0.1;

  rmin = 0.1;
  gmin = 0.3;
  bmin = 0.7;
%
%  Start drawing!
%
  hold ( 'on' );

  for j = 1 : nyc

    y1 = yc(j);
    y2 = yc(j);
    y3 = yc(j+1);
    y4 = yc(j+1);

    for i = 1 : nxc

      x1 = xc(i);
      x2 = xc(i+1);
      x3 = xc(i+1);
      x4 = xc(i);

      z = vc(i,j);
%
%  Draw the top of the prism, using a color corresponding to the height.
%
      r = ( ( zmax - z ) * rmin + ( z - zmin ) * rmax ) / ( zmax - zmin );
      g = ( ( zmax - z ) * gmin + ( z - zmin ) * gmax ) / ( zmax - zmin );
      b = ( ( zmax - z ) * bmin + ( z - zmin ) * bmax ) / ( zmax - zmin );
    
      fill3 ( [ x1, x2, x3, x4 ], [ y1, y2, y3, y4 ], [ z, z, z, z ], z );
%
%  Draw the bottom of the prism, using black.
%
      rgb = zmin;
      fill3 ( [ x1, x2, x3, x4 ], [ y1, y2, y3, y4 ], [ 0, 0, 0, 0 ], rgb );
%
%  Draw the sides of the prism, using a lighter shade of the top color.
%
      rgb = zmin + 0.95 * ( z - zmin );
    
      fill3 ( [ x1, x2, x2, x1 ], [ y1, y2, y2, y1 ], [ 0, 0, z, z ], rgb );
      fill3 ( [ x2, x3, x3, x2 ], [ y2, y3, y3, y2 ], [ 0, 0, z, z ], rgb );
      fill3 ( [ x3, x4, x4, x3 ], [ y3, y4, y4, y3 ], [ 0, 0, z, z ], rgb );
      fill3 ( [ x4, x1, x1, x4 ], [ y4, y1, y1, y4 ], [ 0, 0, z, z ], rgb );

    end
  end

  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'Fontsize', 16 );
  daspect ( [ 1.0, 1.0, zmax ] );
  view ( 3 );
  hold ( 'off' );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );
  
  return
end
function pwc_patch_2d ( nxc, nyc, xc, yc, vc, filename )

%*****************************************************************************80
%
%% pwc_patch_2d plots a piecewise constant function as flat 2D patches.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NXC, NYC, the number of values used to define diffusivity.
%
%    Input, real XC(NXC+1), YC(NYC+1), the diffusivity breakpoint coordinates.
%
%    Input, real VC(NXC,NYC), the diffusivity values.
%
%    Input, string FILENAME, the filename to be used.
%
  figure ( 'Name', 'Diffusivity Patch Plot' );
  clf ( );
%
%  Determine the mapping between colors and values.
%
  zmax = max ( max ( vc ) );
  zmin = min ( min ( vc ) );
  caxis ( [ 0.0, zmax ] );

  for j = 1 : nyc
    y = [ yc(j), yc(j), yc(j+1), yc(j+1) ];
    for i = 1 : nxc
      x = [ xc(i), xc(i+1), xc(i+1), xc(i) ];
      patch ( x, y, vc(i,j), 'EdgeColor', 'k' );
    end
  end

  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  zlabel ( '<--Z(X,Y)-->', 'FontSize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'Fontsize', 16 );
  axis ( 'equal' );
  axis ( 'tight' );
  colorbar ( );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );
  pause ( 5 );
  close ( );

  return
end
function [ quad_w, quad_x, quad_y ] = quad_define ( quad_num )

%*****************************************************************************80
%
%% quad_define defines the quadrature rule.
%
%  Discussion:
%
%    The quadrature rule is given for a reference element, points (X,Y) such that
%
%      0 <= X,
%      0 <= Y, and
%      X + Y <= 1.
%
%      ^
%    1 | *
%      | |\
%    Y | | \
%      | |  \
%    0 | *---*
%      +------->
%        0 X 1
%
%    The rules have the following precision:
%
%    QUAD_NUM  Precision
%
%     1        1
%     3        2
%     4        3
%     6        4
%     7        5
%     9        6
%    13        7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer QUAD_NUM, the number of quadrature nodes.
%    Legal values are 1, 3, 4, 6, 7, 9, 13.
%
%    Output, real QUAD_W(QUAD_NUM), the quadrature weights.
%
%    Output, real QUAD_X(QUAD_NUM), QUAD_Y(QUAD_NUM), the quadrature nodes.
%
  if ( quad_num == 1 )

    quad_x = 1.0 / 3.0;
    quad_y = 1.0 / 3.0;
    quad_w = 1.0;

  elseif ( quad_num == 3 )

    quad_x = [ 0.5; 0.5; 0.0 ];
    quad_y = [ 0.0; 0.5; 0.5 ];
    quad_w = [ 1.0; 1.0; 1.0 ] / 3.0;

  elseif ( quad_num == 4 )

    a =   6.0;
    b =  10.0;
    c =  18.0;
    d =  25.0;
    e = -27.0;
    f =  30.0;
    g =  48.0;

    quad_x = [ b; c; a; a ] / f;
    quad_y = [ b; a; c; a ] / f;
    quad_w = [ e; d; d; d ] / g;

  elseif ( quad_num == 6 )

    a = 0.816847572980459;
    b = 0.091576213509771;
    c = 0.108103018168070;
    d = 0.445948490915965;
    v = 0.109951743655322;
    w = 0.223381589678011;

    quad_x = [ a; b; b; c; d; d ];
    quad_y = [ b; a; b; d; c; d ];
    quad_w = [ v, v, v, w, w, w ];

  elseif ( quad_num == 7 )

    a = 1.0 / 3.0;
    b = ( 9.0 + 2.0 * sqrt ( 15.0 ) ) / 21.0;
    c = ( 6.0 -       sqrt ( 15.0 ) ) / 21.0;
    d = ( 9.0 - 2.0 * sqrt ( 15.0 ) ) / 21.0;
    e = ( 6.0 +       sqrt ( 15.0 ) ) / 21.0;
    u = 0.225;
    v = ( 155.0 - sqrt ( 15.0 ) ) / 1200.0;
    w = ( 155.0 + sqrt ( 15.0 ) ) / 1200.0;

    quad_x = [ a; b; c; c; d; e; e ];
    quad_y = [ a; c; b; c; e; d; e ];
    quad_w = [ u; v; v; v; w; w; w ];

  elseif ( quad_num == 9 )

    a = 0.124949503233232;
    b = 0.437525248383384;
    c = 0.797112651860071;
    d = 0.165409927389841;
    e = 0.037477420750088;

    u = 0.205950504760887;
    v = 0.063691414286223;

    quad_x = [ a; b; b; c; c; d; d; e; e ];
    quad_y = [ b; a; b; d; e; c; e; c; d ];
    quad_w = [ u; u; u; v; v; v; v; v; v ];

  elseif ( quad_num == 13 )

    h = 1.0 / 3.0;
    a = 0.479308067841923;
    b = 0.260345966079038;
    c = 0.869739794195568;
    d = 0.065130102902216;
    e = 0.638444188569809;
    f = 0.312865496004875;
    g = 0.048690315425316;

    w = -0.149570044467670;
    t =  0.175615257433204;
    u =  0.053347235608839;
    v =  0.077113760890257;

    quad_x = [ h; a; b; b; c; d; d; e; e; f; f; g; g ];
    quad_y = [ h; b; a; b; d; c; d; f; g; e; g; e; f ];
    quad_w = [ w; t; t; t; u; u; u; v; v; v; v; v; v ];

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'quad_define - Fatal error!\n' );
    fprintf ( 1, '  No rule is available of order QUAD_NUM = %d\n', ...
      quad_num );
    error ( 'quad_define - Fatal error!\n' );

  end

  return
end
function [ xp, yp ] = reference_to_physical_t3 ( xt, yt, n, xr, yr )

%*****************************************************************************80
%
%% REFERENCE_TO_PHYSICAL_T3 maps a reference point to a physical point.
%
%  Discussion:
%
%    Given the vertices of an order 3 physical triangle and a point
%    (XSI,ETA) in the reference triangle, the routine computes the value
%    of the corresponding image point (X,Y) in physical space.
%
%    Note that this routine may also be appropriate for an order 6
%    triangle, if the mapping between reference and physical space
%    is linear.  This implies, in particular, that the sides of the
%    image triangle are straight and that the "midside" nodes in the
%    physical triangle are halfway along the sides of
%    the physical triangle.
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XT(3), YT(3), the coordinates of the vertices.  The vertices 
%    are assumed to be the images of (0,0), (1,0) and (0,1) respectively.
%
%    Input, integer N, the number of points to transform.
%
%    Input, real XR(N), YR(N), the coordinates of points in the reference space.
%
%    Output, real XP(N), YP(N), the coordinates of the corresponding points in the 
%    physical space.
%
  xp(1:n) = xt(1) * ( 1.0 - xr(1:n) - yr(1:n) ) ...
          + xt(2) *         xr(1:n)             ...
          + xt(3) *                   yr(1:n);

  yp(1:n) = yt(1) * ( 1.0 - xr(1:n) - yr(1:n) ) ...
          + yt(2) *         xr(1:n)             ...
          + yt(3) *                   yr(1:n);

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

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

  return
end
function vp = rhs_func ( np, xp, yp )

%*****************************************************************************80
%
%% rhs_func evaluates the right hand side function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real NP, the number of evaluation points.
%
%    Input, real XP(NP), YP(NN), the evaluation points.
%
%    Output, real VP(NP), the values of F.
%
  vp(1:np) = 1000.0;

  return
end
function rhs_report ( )

%*****************************************************************************80
%
%% rhs_report reports the FEM right hand side function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2019
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'rhs_report:\n' );
  fprintf ( 1, '  rhs(x,y) = 1000\n' );

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, the string to be de-escaped.
%
%    Output, string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function [ ns, xs, ys ] = sensor_location_define ( sensor_location_seed )

%*****************************************************************************80
%
%% sensor_location_define defines the sensor locations.
%
%  Discussion:
%
%    The sensors are modeled as NS points.
%
%    If sensor_location_option is 0, then the sensors are randomly located 
%    within a subrectangle of the region that avoids an outer band of width 
%    1/10 of the relative width and height.
%
%    Otherwise, the sensors are randomly located throughout the region.
%
%    This function calls the random number generator with a specific seed.
%    In order to avoid affecting other random number calculations, this
%    function saves the current seed, uses its seed, and then restores
%    the saved seed.
%    
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SENSOR_LOCATION_SEED, a seed for the random number
%    generator, to be used if the sensor locations will be randomly chosen.
%
%    Output, integer NS, the number of sensors.
%
%    Output, real XS(NS), YS(NS), the sensor coordinates.
%
  sensor_location_option = 1;
%
%  Save the current random number generator status, 
%  and temporarily use the sensor seed.
%
  rng_status = rng ( );
  rng ( sensor_location_seed );

  [ xmin, xmax, ymin, ymax ] = region_box ( );
  ns = 50;

  if ( sensor_location_option == 0 )
    xs = xmin + 0.1 * ( xmax - xmin ) + 0.8 * ( xmax - xmin ) * rand ( ns, 1 );
    ys = ymin + 0.1 * ( ymax - ymin ) + 0.8 * ( ymax - ymin ) * rand ( ns, 1 );
  else
    xs = xmin + ( xmax - xmin ) * rand ( ns, 1 );
    ys = ymin + ( ymax - ymin ) * rand ( ns, 1 );
  end
%
%  Restore the random number generator status.
%
  rng ( rng_status );

  return
end
function sensor_location_plot ( ns, xs, ys, nn, xn, yn, bn, ...
  sensor_location_seed )

%*****************************************************************************80
%
%% sensor_location_plot plots the sensor locations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Input, integer NN, the total number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, integer BN(NN), an indicator for nodal boundary conditions.
%    0, this is an interior node.
%    1, this is a boundary node with zero Dirichlet condition.
%
%    Input, integer sensor_location_seed, the random number seed used
%    when generating the sensor locations.
%
  figure ( 'Name', 'Sensor Locations' );
  clf ( );

  filename = sprintf ( 'sensor_location_%d.png', sensor_location_seed );
%
%  Get the scale.
%
  [ xmin, xmax, ymin, ymax ] = region_box ( );
  margin = 0.05 * max ( xmax - xmin, ymax - ymin );
%
%  Draw the nodes.
%
  hold ( 'on' );
  if ( nn < 150 )
    interior = find ( bn == 0 );
    plot ( xn(interior), yn(interior), 'bo', 'markersize', 5 );
  end
  boundary = find ( bn == 1 );
  plot ( xn(boundary), yn(boundary), 'k.', 'markersize', 20 );
  plot ( xs, ys, 'r.', 'markersize', 20 );
%
%  Label the sensors.
%
  abit = 0.025;
  for j = 1 : ns
    text ( xs(j) + abit, ys(j), num2str(j) );
  end
%
%  Label the plot.
%
  xlabel ( '<---X--->', 'Fontsize', 16 );
  ylabel ( '<---Y--->', 'Fontsize', 16 );
  filename2 = s_escape_tex ( filename );
  datestring = datestr ( datevec ( now ), 0 );
  title ( { filename2; datestring }, 'fontsize', 16 );
  axis ( [ xmin - margin, xmax + margin, ymin - margin, ymax + margin ] );
  axis ( 'equal' );
  grid ( 'on' );
  hold ( 'off' );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  pause ( 5 );
  close ( );

  return
end
function sensor_location_report ( ns, xs, ys, sensor_location_seed )

%*****************************************************************************80
%
%% sensor_location_report reports sensor locations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Input, integer SENSOR_LOCATION_SEED, a seed used for the random number
%    generator if the sensor locations are randomized.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'sensor_location_report:\n' );
  fprintf ( 1, '  sensor_location_seed = %d\n', sensor_location_seed );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    S (    X(s),    Y(s))\n' );
  fprintf ( 1, '\n' );

  for s = 1 : ns
    fprintf ( 1, '  %3d (%8.4f,%8.4f)\n', s, xs(s), ys(s) );
  end

  return
end
function vs = sensor_value_define ( nn, xn, yn, un, ne, e_nodes, ns, xs, ys )

%*****************************************************************************80
%
%% sensor_value_define evaluates the solution at the sensor locations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element solution at the nodes.
%
%    Input, integer NE, the number of elements.
%
%    Input, integer E_NODES(3,NE), the nodes forming each element.
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the finite element solution at the sensor locations.
%
  vs = zeros ( ns, 1 );

  for s = 1 : ns
%
%  Assuming a Delaunay triangulation, and the existence of a
%  element_neighbor() array, we could search efficiently.
%
    e = -1;

    for e_test = 1 : ne
      xt3(1:3) = xn(e_nodes(1:3,e_test));
      yt3(1:3) = yn(e_nodes(1:3,e_test));
      if ( triangle_contains_point ( xt3, yt3, xs(s), ys(s) ) )
        e = e_test;
        break
      end
    end

    if ( e == -1 )
      continue
    end

    if ( false )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Sensor %d is in element %d\n', s, e );
      fprintf ( 1, '  Sensor = (%f,%f)\n', xs(s), ys(s) );
      fprintf ( 1, '  V1 = (%f,%f)\n', xt3(1), yt3(1) );
      fprintf ( 1, '  V2 = (%f,%f)\n', xt3(2), yt3(2) );
      fprintf ( 1, '  V3 = (%f,%f)\n', xt3(2), yt3(3) );
    end

    area = triangle_area ( xt3, yt3 );

    for basis = 1 : 3
      i = e_nodes(basis,e);
      tp1 = 1 + mod ( basis,     3 );
      tp2 = 1 + mod ( basis + 1, 3 );
      bi = ( ( xt3(tp2) - xt3(tp1) ) * ( ys(s) - yt3(tp1) ) ...
           - ( yt3(tp2) - yt3(tp1) ) * ( xs(s) - xt3(tp1) ) ) / area;
      vs(s) = vs(s) + un(i) * bi;
    end

  end

  return
end
function vs = sensor_value_average_define ( nn, xn, yn, un, ne, e_nodes, ...
  ns, xs, ys )

%*****************************************************************************80
%
%% sensor_value_average_define evaluates the solution at the sensor locations.
%
%  Discussion:
%
%    sensor_value_average_define() returns the FEM solution at the
%    sensor location.
%
%    This function, instead, returns the average of several nearby values. 
%    Ideally, it would return an average value based on an integral over
%    a small region.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NN, the number of nodes.
%
%    Input, real XN(NN), YN(NN), the node coordinates.
%
%    Input, real UN(NN), the finite element solution at the nodes.
%
%    Input, integer NE, the number of elements.
%
%    Input, integer E_NODES(3,NE), the nodes forming each element.
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Output, real VS(NS), the finite element solution at the sensor locations.
%
  [ xmin, xmax, ymin, ymax ] = region_box ( );

  nk = 4;
  dk = 0.001 * max ( xmax - xmin, ymax - ymin );

  vs = zeros ( ns, 1 );

  for s = 1 : ns

    xk(1) = ( 1.0 - dk ) * xs(s) + dk * xmin;
    xk(2) = ( 1.0 - dk ) * xs(s) + dk * xmax;
    xk(3) = ( 1.0 - dk ) * xs(s) + dk * xmax;
    xk(4) = ( 1.0 - dk ) * xs(s) + dk * xmin;
    yk(1) = ( 1.0 - dk ) * ys(s) + dk * ymin;
    yk(2) = ( 1.0 - dk ) * ys(s) + dk * ymin;
    yk(3) = ( 1.0 - dk ) * ys(s) + dk * ymax;
    yk(4) = ( 1.0 - dk ) * ys(s) + dk * ymax;

    for k = 1 : 4
%
%  Find E, the index of the element containing the sensor sample point.
%
      e = -1;

      for e_test = 1 : ne
        xt3(1:3) = xn(e_nodes(1:3,e_test));
        yt3(1:3) = yn(e_nodes(1:3,e_test));
        if ( triangle_contains_point ( xt3, yt3, xk(k), yk(k) ) )
          e = e_test;
          break
        end
      end

      if ( e == -1 )
        continue
      end

      if ( false )
        fprintf ( 1, '\n' );
        fprintf ( 1, '  Sensor %d sample point %d is in element %d\n', s, k, e );
        fprintf ( 1, '  Sensor = (%f,%f)\n', xs(s), ys(s) );
        fprintf ( 1, '  Sensor sample point = (%f,%f)\n', xk(k), yk(k) );
        fprintf ( 1, '  V1 = (%f,%f)\n', xt3(1), yt3(1) );
        fprintf ( 1, '  V2 = (%f,%f)\n', xt3(2), yt3(2) );
        fprintf ( 1, '  V3 = (%f,%f)\n', xt3(2), yt3(3) );
      end
%
%  Evaluate the finite element solution in element E at the sensor sample point.
%
      area = triangle_area ( xt3, yt3 );

      for basis = 1 : 3
        i = e_nodes(basis,e);
        tp1 = 1 + mod ( basis,     3 );
        tp2 = 1 + mod ( basis + 1, 3 );
        bi = ( ( xt3(tp2) - xt3(tp1) ) * ( yk(k) - yt3(tp1) ) ...
             - ( yt3(tp2) - yt3(tp1) ) * ( xk(k) - xt3(tp1) ) ) / area;
        vs(s) = vs(s) + un(i) * bi;
      end

    end
%
%  Divide total by number of samples.
%
    vs(s) = vs(s) / nk;

  end

  return
end
function sensor_value_report ( ns, xs, ys, vs )

%*****************************************************************************80
%
%% sensor_value_report reports sensor values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of sensors.
%
%    Input, real XS(NS), YS(NS), the sensor coordinates.
%
%    Input, real VS(NS), the sensor values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'sensor_value_report:\n' );
  fprintf ( 1, '    S (    X(s),    Y(s))      V(s)\n' );
  fprintf ( 1, '\n' );

  for s = 1 : ns
    fprintf ( 1, '  %3d (%8.4f,%8.4f)  %8.4f\n', s, xs(s), ys(s), vs(s) );
  end

  return
end
function area = triangle_area ( xt, yt )

%*****************************************************************************80
%
%% TRIANGLE_AREA computes the area of a triangle in 2D.
%
%  Discussion:
%
%    If the triangle's vertices are given in counterclockwise order,
%    the area will be positive.  If the triangle's vertices are given
%    in clockwise order, the area will be negative!
%
%    An earlier version of this routine always returned the absolute
%    value of the computed area.  I am convinced now that that is
%    a less useful result!  For instance, by returning the signed 
%    area of a triangle, it is possible to easily compute the area 
%    of a nonconvex polygon as the sum of the (possibly negative) 
%    areas of triangles formed by node 1 and successive pairs of vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XT(3), YT(3), the triangle vertices.
%
%    Output, real AREA, the area of the triangle.
%
  area = 0.5 * ( ...
      xt(1) * ( yt(2) - yt(3) ) ...
    + xt(2) * ( yt(3) - yt(1) ) ...
    + xt(3) * ( yt(1) - yt(2) ) );

  return
end
function inside = triangle_contains_point ( xt, yt, xp, yp )

%*****************************************************************************80
%
%% TRIANGLE_CONTAINS_POINT finds if a point is inside a triangle in 2D.
%
%  Discussion:
%
%    The routine assumes that the vertices are given in counter-clockwise
%    order.  If the triangle vertices are actually given in clockwise
%    order, this routine will behave as though the triangle contains
%    no points whatsoever!
%
%    The routine determines if a point P is "to the right of" each of the lines
%    that bound the triangle.  It does this by computing the cross product
%    of vectors from a vertex to its next vertex, and to P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XT(3), YT(3), the triangle vertices.
%    The vertices should be given in counter clockwise order.
%
%    Input, real X, Y, the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside
%    the triangle or on its boundary.
%
  for j = 1 : 3

    k = mod ( j, 3 ) + 1;

    if ( 0.0 < ( xp - xt(j) ) * ( yt(k) - yt(j) ) ...
             - ( yp - yt(j) ) * ( xt(k) - xt(j) ) )
      inside = false;
      return
    end

  end

  inside = true;

  return
end
