triangulation_corner.m.txt

a MATLAB program which tries to correct situations in which a triangulation includes corner triangle

function triangulation_corner ( input_node_filename, input_triangle_filename )

%% MAIN is the main program for TRIANGULATION_CORNER
%
%  Discussion:
%
%    TRIANGULATION_TRIANGLE_CORNER tries to alter a triangular mesh in
%    cases where two sides of a triangle lie on the boundary.
%
%    The program reads the node and triangle data, computes the triangle
%    neighbor information, and writes it to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Usage:
%
%    triangulation_corner ( 'node_file', 'triangle_file' )
%
  debug = 0;

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGULATION_CORNER\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read a node dataset of NODE_NUM points in 2 dimensions.\n' );
  fprintf ( 1, '  Read an associated triangle file of TRIANGLE_NUM \n' );
  fprintf ( 1, '  triangles using 3 or 6 nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Any triangle which has exactly two sides on the boundary\n' );
  fprintf ( 1, '  is a corner triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If there are any corner triangles this program tries to\n' );
  fprintf ( 1, '  eliminate them.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1,'  The "repaired" triangle file is written out.\n' );
%
%  Get the number of command line arguments.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
    input_node_filename = input ( '  Please enter the name of the node file.' );

  end
%
%  If at least two command line arguments, the second is the triangulation file.
%
  if ( nargin < 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
    input_triangle_filename = input ( ...
      '  Please enter the name of the triangulation file.' );

  end
%
%  Read the node data.
%
  [ dim_num, node_num ] = dtable_header_read (  input_node_filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', input_node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  node_xy(1:dim_num,1:node_num) = dtable_data_read ( input_node_filename, ...
    dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', input_node_filename );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, dim_num, 5, ...
    '  First 5 nodes:' );
%
%  Read the triangulation data.
%
  [ triangle_order, triangle_num ] = ...
    itable_header_read ( input_triangle_filename );

  if ( triangle_order ~= 3 & triangle_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_CORNER - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 3-node or 6-node triangulation.\n' );
    error ( 'TRIANGULATION_L2Q - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of ""%s".\n', input_triangle_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Triangle order = %d\n', triangle_order );
  fprintf ( 1, '  Number of triangles TRIANGLE_NUM  = %d\n', triangle_num );

  triangle_node(1:triangle_order,1:triangle_num) = itable_data_read ( ...
    input_triangle_filename, triangle_order, triangle_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in ""%s".\n', input_triangle_filename );

  i4mat_transpose_print_some ( triangle_order, triangle_num, triangle_node, ...
    1, 1, triangle_order, 5, '  First 5 triangles:' );
%
%  Create the triangle neighbor array.
%
  if ( triangle_order == 3 )

    triangle_neighbor = triangulation_order3_neighbor_triangles ( ...
      triangle_num, triangle_node );

  elseif ( triangle_order == 6 )

    triangle_neighbor = triangulation_order6_neighbor_triangles ( ...
      triangle_num, triangle_node );

  end
%
%  Examine the triangle neighbor array.
%
  OFFSET = 1;
  negative_total(OFFSET+0:OFFSET+3) = 0;

  for triangle = 1 : triangle_num

    negative = 0;
    for neighbor = 1 : 3
      if ( triangle_neighbor(neighbor,triangle) < 0 )
        negative = negative + 1;
      end
    end

    negative_total(OFFSET+negative) = negative_total(OFFSET+negative) + 1;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of boundary sides     Number of triangles\n' );
  fprintf ( 1, '\n' );
  for i = 0 : 3
    fprintf ( 1, '            %8d            %8d\n', i, negative_total(OFFSET+i) );
  end
%
%  Try to patch problems.
%
  if ( 0 < negative_total(OFFSET+3) )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_CORNER - Fatal error!\n' );
    fprintf ( 1, '  There is at least one triangle with all sides\n' );
    fprintf ( 1, '  on the boundary.\n' );
    return

  elseif ( 0 == negative_total(OFFSET+2) )

    fprintf ( 1,  '\n' );
    fprintf ( 1,  'TRIANGULATION_CORNER:\n' );
    fprintf ( 1,  '  No corner triangles were found.\n' );
    fprintf ( 1,  '  No corrections need to be made.\n' );
    return

  else
%
%  We need the triangles to be oriented properly.
%
    negative = 0;

    for triangle = 1 : triangle_num

      t3(1:2,1:3) = node_xy(1:2,triangle_node(1:3,triangle));

      area = triangle_area_2d ( t3 );

      if ( area < 0.0 )

        negative = negative + 1;

        node                      = triangle_node(2,triangle);
        triangle_node(2,triangle) = triangle_node(3,triangle);
        triangle_node(3,triangle) = node;

        if ( triangle_order == 6 )
          node                      = triangle_node(4,triangle);
          triangle_node(4,triangle) = triangle_node(6,triangle);
          triangle_node(6,triangle) = node;
        end

        neighbor                      = triangle_neighbor(1,triangle);
        triangle_neighbor(1,triangle) = triangle_neighbor(3,triangle);
        triangle_neighbor(3,triangle) = neighbor;

      end

    end

    if ( 0 < negative )
      fprintf ( 1,  '\n' );
      fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
      fprintf ( 1, '  Reoriented %d triangles.', negative );
      fprintf ( 1,  '\n' );
    else
      fprintf ( 1,  '\n' );
      fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
      fprintf ( 1, '  Triangles were already properly oriented.\n' );
      fprintf ( 1,  '\n' );
    end
%
%  Now consider each triangle that has exactly two boundary sides.
%
    for triangle1 = 1 : triangle_num 

      negative = 0;
      for neighbor = 1 : 3
        if ( triangle_neighbor(neighbor,triangle1) < 0 )
          negative = negative + 1;
        end
      end

     if ( negative == 2 )

       triangle2 = -1;

       for neighbor = 1 : 3
         if ( 0 < triangle_neighbor(neighbor,triangle1) )
           triangle2 = triangle_neighbor(neighbor,triangle1);
           t1_to_t2 = neighbor;
         end
       end

       fprintf ( 1, '  Adjusting triangle %d using triangle %d\n', ...
         triangle1, triangle2 );

       t2_to_t1 = -1;
       for neighbor = 1 : 3
         if ( triangle_neighbor(neighbor,triangle2) == triangle1 )
           t2_to_t1 = neighbor;
         end
       end

       n1_index = i4_wrap ( t1_to_t2 - 1, 1, 3 );
       node = triangle_node(n1_index,triangle1);;

       n1_index = i4_wrap ( t1_to_t2 + 1, 1, 3 );
       n2_index = i4_wrap ( t2_to_t1 - 1, 1, 3 );
       triangle_node(n1_index,triangle1) = triangle_node(n2_index,triangle2);

       n1_index = i4_wrap ( t1_to_t2 - 1, 1, 3 );
       n2_index = i4_wrap ( t2_to_t1 + 1, 1, 3 );
       triangle_node(n2_index,triangle2) = node;

       if ( triangle_order == 6 )
%
%  Adjust coordinates of the new midside node.
%
         j1 = i4_wrap ( t1_to_t2 - 1, 1, 3 );
         j2 = i4_wrap ( t1_to_t2 + 3, 4, 6 );
         j3 = i4_wrap ( t2_to_t1 - 1, 1, 3 );

         node1 = triangle_node(j1,triangle1);
         node2 = triangle_node(j2,triangle1);
         node3 = triangle_node(j3,triangle2);

         node_xy(1,node2) = 0.5D+00 * ( node_xy(1,node1) + node_xy(1,node3) );
         node_xy(2,node2) = 0.5D+00 * ( node_xy(2,node1) + node_xy(2,node3) );
%
%  Update the triangle array.
%
         j1 = i4_wrap ( t1_to_t2 + 4, 4, 6 );
         j2 = i4_wrap ( t1_to_t2 + 3, 4, 6 );
         j3 = i4_wrap ( t2_to_t1 + 4, 4, 6 );
         j4 = i4_wrap ( t2_to_t1 + 3, 4, 6 );

         node                        = triangle_node(j1,triangle1);
         triangle_node(j1,triangle1) = triangle_node(j2,triangle1);
         triangle_node(j2,triangle1) = triangle_node(j3,triangle2);
         triangle_node(j3,triangle2) = triangle_node(j4,triangle2);
         triangle_node(j4,triangle2) = node;

       end
%
%  Update the neighbor array.
%
       n2_index = i4_wrap ( t2_to_t1 + 1, 1, 3 );
       triangle = triangle_neighbor(n2_index,triangle2);
       triangle_neighbor(n2_index,triangle2) = triangle1;
       triangle_neighbor(t2_to_t1,triangle2) = -1;

       n1_index = i4_wrap ( t1_to_t2 + 1, 1, 3 );
       triangle_neighbor(n1_index,triangle1) = triangle2;
       triangle_neighbor(t1_to_t2,triangle1) = triangle;

     end

    end
%
%  Write out the corrected triangle file.
%
    output_triangle_filename = file_name_ext_swap ( input_triangle_filename, ...
      'corner.txt' );

    header = 0;

    itable_write ( output_triangle_filename, triangle_order, triangle_num, ...
      triangle_node, header );

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
    fprintf ( 1, '  New triangle file with repaired corners written to \n' );
    fprintf ( 1, '  "%s".\n', output_triangle_filename );
%
%  For quadratic triangles, write out the corrected node coordinate file.
%
    if ( triangle_order == 6 )

      output_node_filename = file_name_ext_swap ( input_node_filename, ...
        'corner.txt' );

      header = 0;

      dtable_write ( output_node_filename, dim_num, node_num, node_xy, header );

      fprintf ( 1, '\n' );
      fprintf ( 1, 'TRIANGULATION_CORNER:\n' );
      fprintf ( 1, '  New node coordinate file with adjusted midside nodes\n' );
      fprintf ( 1, '  written to "%s".\n', output_node_filename );

    end

  end

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

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

  return
end