%% TRIANGULATION_ORDER6_CONTOUR contour plots data at the nodes of a six node triangle.
%
%  Discussion:
%
%    This program can read data files defining a set of nodes, the triangulation of
%    those nodes using 6 node triangles, and the value of a scalar at the nodes.  
%    It then rearranges the data to set up calls to MATLAB that display the
%    triangulation, and color contour plots of the scalar.
%
%    The program reads three data files:
%
%    * nodes.txt, containing 
%      X, Y coordinates of points;
%      
%    * elements.txt, containing a list of sets of six nodes
%      making up each triangular element;
%      
%    * solution.txt, containing the value of a scalar function
%      U(X,Y) evaluated at each node.
%
%    The MATLAB plotting routines require a 3-node triangulation.  This program
%    shows two ways to turn a 6-node triangulation into a 3-node triangulation.
%    The default MATLAB plot of data on a triangular mesh uses constant color
%    over each triangle.  After displaying these plots, the program shows how to
%    get a much nicer plot using interpolated color.
%
%  Modified:
%
%    24 March 2005
%
%  Author:
%
%    John Burkardt
%
  timestamp;
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_ORDER6_CONTOUR:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot a scalar defined on a 6-node triangle triangulation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program expects to find three files to read:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * "nodes.txt",     the node file,\n' );
  fprintf ( 1, '  * "elements.txt",  the element file,\n' );
  fprintf ( 1, '  * "solution.txt",  the solution file,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  It reads the files, and makes two plots of the solution.\n' );
%
%  Read the data.
%
  p = load ( 'nodes.txt' );
  t = load ( 'elements.txt' );
  u = load ( 'solution.txt' );
%
%  Make a 3-node triangulation by discarding the midside nodes.
%
  t2 = t(:,1:3);
%
%  Display the mesh.
%
  trimesh ( t2, p(:,1), p(:,2) );

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Dropping midside nodes', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Display the solution on this crude mesh.
%
  trisurf ( t2, p(:,1), p(:,2), u )
  
  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'FEM solution U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Make a 3-node triangulation by breaking each 6 node triangle 
%  into 4 3-node triangles.
%
  tri_num1 = size ( t, 1 );
  t3 = triangulation_order6_to_order3 ( tri_num1, t' );
  t3 = t3';
%
%  Display the mesh.
%
  trimesh ( t3, p(:,1), p(:,2) );
  
  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Splitting 6-Node Triangles', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Display the solution on this finer mesh.
%
  trisurf ( t3, p(:,1), p(:,2), u );
  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'FEM solution U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause
%
%  Make a nicer plot on the finer mesh by using color interpolation.
%
  trisurf ( t3, p(:,1), p(:,2), u, 'FaceColor', 'interp', ...
    'EdgeColor', 'interp' )

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  zlabel ( 'U', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'FEM solution U(X,Y)', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return...\n' );
  pause

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

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