function [ u, v, p ] = box_flow ( total_time, nu, Lambda, ... 
  solver_type, save_file, node_file, element_file )

%*****************************************************************************80
%
%% BOX_FLOW is the main program for a fluid flow calculation.
%
%  Discussion:
%
%    This program solves the Navier-Stokes equations (NSE), with the continuity
%    equation replaced by a pressure Poisson equation (PPE).  Appropriate 
%    boundary conditions are developed for the PPE, which allow for a fully
%    decoupled numerical scheme to recover the pressure.  The variational
%    form of the NSE with PPE is used in the Galerkin finite element
%    discretization.  The Galerkin finite element method is then used to
%    solve the NSE with PPE.
%
%  Modified:
%
%    03 December 2018
%
%  Modified:
%
%    John Burkardt
%
%  Author:
%
%    John Cornthwaite
%
%  Reference:
%
%    John Cornthwaite,
%    Pressure Poisson Method for the Incompressible Navier-Stokes Equations
%    Using Galerkin Finite Elements, (2013)
%    Electronic Theses and Dissertations 831,
%    https://digitalcommons.georgiasouthern.edu/etd/831
%
%  Parameters:
%
%    Input, real total time : amount of time to be evaluated.
%
%    Input, real nu: viscosity.
%
%    Input, real Lambda: stability parameter (0-100).
%
%    Input, integer solver type : 1-Backward-Euler , 2-Crank-Nicolson.
%
%    Input, string, save_file, name to save workspace.  Use '' if no save.
%
%    Input, string node_file, the name of the node file.
%
%    Input, string element_file, the name of the element file.
%
%    Output, real u(), v(): horizontal and vertical velocity.
%
%    Output, real p(): pressure.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_FLOW:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Navier-Stokes solver using Pressure Poisson Equation.\n' );
%
%  Gauss quadrature rule ; only 7 is allowed.
%
  quad_num = 7;
%
%  Reynolds number.
%
  Re = 1.0 / nu;

  node_xy = dlmread ( node_file )';
  [ node_m, node_n ] = size ( node_xy );
  fprintf ( 1, '  %d nodes\n', node_n );

  element_node = dlmread ( element_file )';
  [ element_m, element_n ] = size ( element_node );
  fprintf ( 1, '  %d elements of order %d\n', element_n, element_m );
%
%  Determine normal vectors, node_numbering, boundaries, 
%  various element information.
%
  [ pres, pres_num, node_num, element_num, node_boundary, dt, ...
    steps, normal_node, variable_num, node_u_variable, ...
    node_v_variable, node_p_variable ] ....
    = Set_Elements ( total_time, element_node, node_xy );
%
%  Compute mass, stiffness, and gradient matrices.
%
  [ M, D, GP ] = Set_Matrix ( quad_num, element_node, node_xy, ...
    node_num, node_u_variable, node_v_variable, ...
    node_p_variable, element_num, pres_num );
%
%  Invert M.
%
  im = M \ speye ( 2*node_num, 2*node_num );
%
%  Laplace operator for pressure.
%
  Mp = GP' * im * GP;

  pbound = logical ( node_boundary(pres) );
  removebnd = Mp(:,pbound);
  Mp(pbound,:) = 0;
  Mp(:,pbound) = 0;
  Mp(pbound,pbound) = diag ( ones ( nnz(node_boundary(pres)), 1 ) ); 
%
%  Initial conditions
%
%  initial velocity. initial pressure is solved for
%
  [ u, v, ~ ] = external_force ( node_xy, 0, 1, Re );
  node_c = [ u'; v'; zeros(pres_num,1) ];
%
%  Solve for initial pressure
%
  [ Ku, Kv ] = Set_PPE_RHS ( Re, quad_num, element_node, ...
    node_xy, node_num, node_p_variable, node_u_variable, ...
    element_num, pres_num, normal_node, Lambda );
%
%  Compute external force funtion for RHS of momentum equation.
%
  [ ~, ~, F3 ] = rhs ( 0, dt, node_num, node_u_variable, ...
    node_v_variable, node_p_variable, pres_num, element_num, ...
    element_node, quad_num, node_xy, Re );
%
%  Right hand side of PPE.
%
  K1 = Ku * node_c(1:node_num) + Kv * node_c(node_num+1:2*node_num);
%
%  Adjust the boundary.
%
  K2 = removebnd * K1(pbound);
  K = F3 - K2;
  K(pbound) = K1(pbound);
%
%  Set first pressure node to zero.
%
  K(1) = 0;

  P = Mp \ K ;
  node_c(2*node_num+1:variable_num,1) = P;
%
%  END Initial Conditions
%
%  Solve the system.
%
  if ( solver_type == 1 )
    node_c = Backward_Euler ( Re, M, D, GP, Mp, Ku, Kv, node_c, ...
      removebnd, pbound, node_boundary, steps, dt, variable_num, ...
      element_node, node_xy, quad_num, node_num, element_num, ...
      pres_num, node_p_variable, node_u_variable, ...
      node_v_variable );
  else
    node_c = Crank_Nick ( Re, M, D, GP, Mp, Ku, Kv, node_c, ...
      removebnd, pbound, node_boundary, steps, dt, variable_num, ...
      element_node, node_xy, quad_num, node_num, element_num, ... 
      pres_num, node_p_variable, node_u_variable, ...
      node_v_variable );
  end

  u = node_c(1:node_num);
  v = node_c(1+node_num:2*node_num);
  p = node_c(1+2*node_num:end);
%
%  Post-processing
%
  plotdata ( u, v, p, node_xy, node_num, pres, steps, dt, Re );
%
%  Save data if argument exists.
%
  if ( 0 < length ( save_file ) ) 
    save ( save_file );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_FLOW:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function [ pres, pres_num, node_num, element_num, ...
  node_boundary, dt, steps, normal_node, variable_num, ... 
  node_u_variable, node_v_variable, node_p_variable ] ...
  = Set_Elements ( total_time, element_node, node_xy )

%*****************************************************************************80
%
%% SET_ELEMENTS sets up the element information.
%
%  Modified:
%
%    03 December 2018
%
  pres = unique ( element_node(1:3,:) ); %pressure nodes
  pres_num = length(pres); % number of pressure nodes
  node_num = max(max(element_node )); %total number of nodes 
  element_num = size (element_node ,2); %number of elements
  node_boundary = triangulation_order6_boundary_node ( node_num, ...
    node_xy', element_node );
%
%  Determine the time step
%
  del_x = min_dist ( node_xy, element_node );
  dt = 0.95 * del_x ^ 2;
  steps = ceil ( total_time / dt ) - 1;
%
%  Determine the spatial range.
%
  max_x = max ( node_xy(1,:) ); 
  min_x = min ( node_xy(1,:) ); 
  max_y = max ( node_xy(2,:) );
  min_y = min ( node_xy(2,:) );
%
%  Normal vectors.
%
  normal_node = zeros ( 2, node_num );

  normvec = ( node_xy(1,:) == max_x );
  normal_node(1,:) = normvec;
  normvec = ( node_xy(1,:) == min_x );
  normal_node(1,:) = normal_node(1,:) - normvec; 

  normvec = ( node_xy(2,:) == max_y );
  normal_node(2,:) = normvec;
  normvec = ( node_xy(2,:) == min_y );
  normal_node(2,:) = normal_node(2,:) - normvec;

  normvec = abs ( normal_node(1,:) ) + abs ( normal_node(2,:) ); 
  normal_node(:,normvec==2) = 0;
%
%  Handle corners.
%
  normvec = find ( node_xy(1,:) == max_x & node_xy(2,:) == max_y );
  if ( 0 < length ( normvec ) )
    normal_node(:,normvec) = [ +1, +1 ]';
  end
  normvec = find ( node_xy(1,:) == min_x & node_xy(2,:) == max_y );
  if ( 0 < length ( normvec ) )
    normal_node(:,normvec) = [ -1, +1 ]';
  end
  normvec = find ( node_xy(1,:) == min_x & node_xy(2,:) == min_y );
  if ( 0 < length ( normvec ) )
    normal_node(:,normvec) = [ -1, -1 ]';
  end
  normvec = find ( node_xy(1,:) == max_x & node_xy(2,:) == min_y );
  if ( 0 < length ( normvec ) )
    normal_node(:,normvec) = [ +1, -1 ]';
  end
%
%  Number variables.
%
  variable_num = 2 * node_num + pres_num;

  node_u_variable = 1:1:node_num;
  node_v_variable = node_num+1:1:2*node_num;
  node_p_variable = zeros ( 1, node_num );
  node_p_variable(pres) = 1+2*node_num:1:variable_num;

  return 
end
function [ u, v, p ] = external_force ( xy,t, type, Re )

%*****************************************************************************80
%
%% EXTERNAL_FORCE evaluates the external force.
%
%  Modified:
%
%    03 December 2018
%
  x = xy(1,:);
  y = xy(2,:);
  ct = cos(t);
  st = sin(t);
  s2py = sin(2*pi*y);
  spx = sin(pi*x); 
  s2px = sin(2*pi*x);
  spiy = sin(pi*y); 
  cpx = cos(pi*x); 
  c2px = cos(2*pi*x); 
  cpy = cos(pi*y); 
  c2py = cos(2*pi*y);

  if ( type == 1 )

    u =   pi * ct * s2py .* spx .^ 2;
    v = - pi * ct * s2px .* spiy .^ 2; 
    p = - ct * cpx .* spiy;

    dudx =   2 * pi * ct .* s2py .* spx .* cpx;
    dvdy = - 2 * pi * ct .* s2px .* spiy .* cpy;
    divu = norm ( dudx + dvdy, inf );

  else

    ut = - pi * st * s2py .* spx .^ 2;
    vt =   pi * st * s2px .* spiy .^ 2;
    lapu = 2 * pi ^ 3 * ct * s2py .* ( cpx .^ 2 - spx .^ 2 ) ...
      - 4 * pi ^ 3 * ct .* s2py .* spx .^ 2;
    lapv = 4 * pi ^ 3 * ct * s2px .* spiy .^ 2 ...
      - 2 * pi ^ 3 * ct * s2px .* ( cpy .^ 2 - spiy .^ 2 );
    px =   pi * ct * spx .* spiy;
    py = - pi * ct * cpx .* cpy;
    u = ut + px - ( 1 / Re ) * lapu;
    v = vt + py - ( 1 / Re ) * lapv;
    p = 0.0;

  end

  return
end
function [ M, D, GP ] = Set_Matrix ( quad_num, element_node, ...
  node_xy, node_num, node_u_variable, ...
  node_v_variable, node_p_variable, element_num, pres_num )

%*****************************************************************************80
%
%% SET_MATRIX computes the mass, diffusion, and gradient matrices.
%
%  Modified:
%
%    03 December 2018
%
  vec_len = element_num * 36;
  mu = zeros ( vec_len, 3 );
  mv = mu;
  uu = mu;
  vv = mu;
  row_cnt = 1;
  vec_len = element_num * 18;
  pu = zeros ( vec_len, 3 );
  pv = pu;
  row_cnt2 = 1;
  velocity_dof = 2 * node_num;
  [ quad_w, quad_xy ] = quad_rule ( quad_num );

  for element = 1 : element_num
%
%  Make a copy of the triangle.
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
    t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element));
%
%  Map the quadrature points QUADXY to points XY in the physical triangle.
%
    xy(1:2,1:quad_num) = reference_to_physical_t6 ( t6, quad_num, quad_xy );
    area = abs ( triangle_area_2d ( t3 ) ); 
    w(1:quad_num) = area * quad_w(1:quad_num);
%
%  Evaluate the basis functions at the quadrature points.
%
    [ Psi, dPsidx, dPsidy ] = basis_mn_t6 ( t6, quad_num, xy );
    [ ~, dPhidx, dPhidy ] = basis_mn_t3 ( t3, quad_num, xy );

    iu(1:6) = node_u_variable(element_node(1:6,element)); 
    iv(1:6) = node_v_variable(element_node(1:6,element)); 
    ip(1:3) = node_p_variable(element_node(1:3,element)) - velocity_dof;

    for i = 1 : 6
      for j = 1 : 6
%
%  Mass matrix.
%
        mu(row_cnt,1:2) = [iu(i),iu(j)];
        mu(row_cnt,3) = sum ( w(1:quad_num) ...
          .* Psi(i,1:quad_num) .* Psi(j,1:quad_num)); 
        mv(row_cnt,1:2) = [ iv(i), iv(j) ];
        mv(row_cnt,3) = mu(row_cnt,3);
%
%  Stiffness matrix.
%
        uu(row_cnt,1:2) = [ iu(i), iu(j) ];
        uu(row_cnt,3) = sum ( w(1:quad_num) ...
          .* ( dPsidx(i,1:quad_num) ...
          .* dPsidx(j,1:quad_num) + dPsidy(i,1:quad_num) ... 
          .* dPsidy(j,1:quad_num) ) );

        vv(row_cnt,1:2) = [ iv(i), iv(j) ]; 
        vv(row_cnt,3) = uu(row_cnt,3);

        if ( i < 4 )
%
%  Gradient of P.
%
          pu(row_cnt2,1:2) = [ iu(j), ip(i) ];
          pu(row_cnt2,3) = sum ( w(1:quad_num) ... 
            .* dPhidx(i,1:quad_num) .* Psi(j,1:quad_num));
          pv(row_cnt2,1:2) = [ iv(j), ip(i) ];
          pv(row_cnt2,3) = sum ( w(1:quad_num) ...
            .* dPhidy(i,1:quad_num) ...
            .* Psi(j,1:quad_num));
          row_cnt2 = row_cnt2 + 1;
        end
        row_cnt = row_cnt + 1; 
      end
    end
  end

  mu = mu ( any ( mu(:,3), 2 ), : ); 
  mv = mv ( any ( mv(:,3), 2 ), : ); 
  mu = [ mu; mv ];

  M = sparse ( mu(:,1), mu(:,2), mu(:,3), 2 * node_num, 2 * node_num );

  uu ( any ( uu(:,3), 2 ), : ); 
  vv ( any ( vv(:,3), 2 ), : );
  uu = [ uu; vv ];

  D = sparse ( uu(:,1), uu(:,2), uu(:,3), 2 * node_num, 2 * node_num );

  pu = pu ( any ( pu(:,3), 2 ), : ); 
  pv = pv ( any ( pv(:,3), 2 ), : ); 
  pu = [ pu; pv ];

  GP = sparse ( pu(:,1), pu(:,2), pu(:,3), 2 * node_num, pres_num );

  return
end
function [ node_c ] = Backward_Euler ( Re, M, D, GP, Mp, Ku, Kv, ...
  node_c, removebnd, pbound, node_boundary, steps, dt, ... 
  variable_num, element_node, node_xy, quad_num, node_num, ... 
  element_num, pres_num, node_p_variable, ...
  node_u_variable, node_v_variable )

%*****************************************************************************80
%
%% BACKWARD_EULER carries out the backward Euler method.
%
%  Modified:
%
%    03 December 2018
%
  boundary_size = nnz ( [node_boundary,node_boundary] ); 
  lbound = logical ( [node_boundary,node_boundary]' );
%
%  Pre-solve the LHS momentum equation. 
%
  Lhs = ( M + ( dt / Re ) * D );
%
%  Adjust the boundary.
%
  Lhs(lbound,:) = 0;
  Lhs(:,lbound) = 0;
%
%  Adjust the boundary conditions.
%
  Lhs(lbound,lbound) = diag ( ones(boundary_size,1) );

  for time = 1 : steps

    tic
%
%  Compute external force funtion for RHS of momentum equation.
%
    [ F, ~, F3 ] = rhs ( time * dt, dt, node_num, node_u_variable, ...
      node_v_variable, node_p_variable, pres_num, ...
      element_num, element_node, quad_num, node_xy, Re );

    Rhs = M * node_c(1:2*node_num) - dt * GP ... 
      * node_c(1+2*node_num:end) + dt * F;
%
%  Adjust boundary.
%
    Rhs(lbound) = 0;

    node_c(1:2*node_num) = Lhs \ Rhs;
%
%  Solve the PPE.
%
    K1 = Ku * node_c(1:node_num) + Kv * node_c(node_num+1:2*node_num); 
    K2 = removebnd * K1(pbound);
    K = F3 - K2;
    K(pbound) = K1(pbound);
    K(1) = 0;
    P = Mp \ K ;
    node_c(2*node_num+1:variable_num,1) = P;
    tm = toc;
    fprintf ( 'step %i of %i completed in %3.4f s \n', time, steps, toc ) 

  end

  return
end
function [ node_c ] = Crank_Nick ( Re, M, D, GP, Mp, Ku, Kv, node_c, ...
  removebnd, pbound, node_boundary, steps, dt, variable_num, ...
  element_node, node_xy, quad_num, node_num, element_num, ...
  pres_num , node_p_variable, node_u_variable, node_v_variable )

%*****************************************************************************80
%
%% CRANK_NICK implements the Crank-Nicolson / Adams-Bashforth method.
%
%  Modified:
%
%    03 December 2018
%
  solution = zeros ( variable_num, 2 ); %current and previous
  boundary_size = nnz ( [ node_boundary, node_boundary ] );
  lbound = logical ( [ node_boundary, node boundary ]' );
  solution(:,2) = node_c;
  Lhs2 = ( M + .5 * dt / Re * D ); % pre-solve LHS momentum eq. 
  Lhs2(lbound,:) = 0; %adjust boundary
  Lhs2(:,lbound) = 0; %adjust boundary
  Lhs2(lbound,lbound) = diag ( ones ( boundary_size, 1 ) ); 
  Lhs2 = Lhs2 \ speye ( 2 * node_num, 2 * node_num );

  for time = 1 : steps

    tic

    if ( time == 1 )
      [ F, ~, F3 ] = rhs ( time * dt, dt, node_num, node_u_variable, ...
        node_v_variable, node_p_variable, pres_num, ... 
        element_num, element_node, quad_num, node_xy, Re );
      Rhs = M * node_c(1:2*node_num) - dt * GP ... 
        * node_c(1+2*node_num:end) + dt * F;
      Rhs(lbound) = 0; %adjust boundary
      Lhs = ( M + ( dt / Re ) * D ); % pre-solve LHS momentum eq. 
      Lhs(lbound,:) = 0; %adjust boundary
      Lhs(:,lbound) = 0; %adjust boundary
      Lhs(lbound,lbound) = diag ( ones ( boundary_size, 1) ); 
      node c(1:2*node_num) = Lhs \ Rhs;
    else
      [ F, F2, F3 ] = rhs ( time * dt, dt, node_num, ...
        node_u_variable, node_v_variable, node_p_variable, ...
        pres_num, element_num, element_node, quad_num, node_xy, Re );
      F = .5 * ( F + F2 );
      Rhs = ( M - .5 * dt / Re * D ) * node_c(1:2*node_num,1) ... 
        + dt * GP * ( - 1.5 * node_c(1+2*node_num:end) ...
        + .5 * solution(1+2*node_num:end,1)) + dt * F;
      Rhs(lbound) = 0;
      node_c(1:2*node_num,1) = Lhs2 * Rhs; 
    end
%
%  Solve the PPE.
%
    K1 = Ku * node_c(1:node_num) + Kv * node_c(node_num+1:2*node_num);
    K2 = removebnd * K1(pbound);
    K = F3 - K2;
    K(pbound) = K1(pbound);
    K(1) = 0;
    P = Mp \ K;
    node_c(2*node_num+1:variable_num,1) = P;
    solution = [ solution(:,2), node_c ];
    tm = toc;
    fprintf ( 'step %i of %i completed in %3.4f s \n', time, steps, toc )

  end
 
  return
end
function [ Ku, Kv ] = Set_PPE_RHS ( Re, quad_num, element_node, ...
  node_xy, node_num, node_p_variable, node_u_variable, ... 
  element_num, pres_num, normal_node, Lambda )

%*****************************************************************************80
%
%% SET_PPE_RHS sets the right hand side for the Pressure Poisson equation.
%
%  Modified:
%
%    03 December 2018
%
  vec_len = element_num * 18;
  uu = zeros ( vec_len, 3 );
  vv = uu;
  row_cnt = 1;
  velocity_dof = 2 * node_num;
%
%  Get the quadrature weights and nodes.
%
  [ quad_w, quad_xy ] = quad_rule ( quad_num );
%
%  Consider all quantities associated with a given ELEMENT.
%
  for element = 1 : element_num
%
%  Make a copy of the triangle .
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
    t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element));
%
%  Map the quad points QUADXY to points XY in the physical triangle.
%
    xy(1:2 ,1:quad_num) = reference_to_physical_t6 ( t6, quad_num, quad_xy );
    area = abs ( triangle_area_2d ( t3 ) );
    w(1:quad_num) = area * quad_w(1:quad_num);
%
%  Evaluate the basis functions at the quadrature points. 
%
    [ Psi, dpsidx, dpsidy ] = basis_mn_t6 ( t6, quad_num, xy );
    [ Phi, dphidx, dphidy ] = basis_mn_t3 ( t3, quad_num, xy );
%
%  Extract the node indices for this element.
%
    iu(1:6) = node_u_variable(element_node(1:6,element)); 
    ip(1:3) = node_p_variable(element_node(1:3,element)) - velocity_dof;
%
%  Normal vector for each node in the element.
%
    normal_vec = normal_node(:,element_node(1:6,element));
%
%  Compute <curl(u), n * grad(Phi)> + lambda<n * u,Phi>.
%
    for i = 1 : 3
      for j = 1 : 6

        uu(row_cnt,1:2) = [ ip(i), iu(j) ];
        dphidtau = normal_vec(1,i) * dphidy(i,1:quad_num) ...
          - normal_vec(2,i) * dphidx(i,1:quad_num);
        uu(row_cnt,3) = sum ( w(1:quad_num) .* ( - dphidtau ) ...
          .* dpsidy(j,1:quad_num) + Lambda ... 
          * normal_vec(1,j) * Psi(j,1:quad_num) ...
          .* Phi(i,1:quad_num) );

        vv(row_cnt,1:2) = [ ip(i), iu(j) ];
        vv(row_cnt ,3) = sum ( w(1:quad_num) .* dphidtau ...
          .* dpsidx(j,1:quad_num) + Lambda ... 
          * normal_vec(2,j) * Psi(j,1:quad_num) ...
          .* Phi(i,1:quad_num));

        row_cnt = row_cnt + 1;
      end
    end

  end

  uu = uu ( any ( uu(:,3), 2 ), : );
  vv = vv ( any ( vv(:,3), 2 ), : );

  Ku = ( 1.0 / Re ) * sparse ( uu(:,1), uu(:,2), uu(:,3), pres_num, node_num ); 
  Kv = ( 1.0 / Re ) * sparse ( vv(:,1), vv(:,2), vv(:,3), pres_num, node_num );

  return 
end
function [ F1, F2, F3 ] = rhs ( t, dt, node_num, node_u_variable, ... 
  node_v_variable, node_p_variable, pres_num, element_num, ... 
  element_node, quad_num, node_xy, Re )

%*****************************************************************************80
%
%% RHS evaluates the right hand side.
%
%  Modified:
%
%    03 December 2018
%
  F1 = zeros(2*node_num,1); % F at time = t
  F2 = F1; % F at time = t + .5dt
  F3 = zeros ( pres_num,1); % F at time = t + dt
  [ quad_w, quad_xy ] = quad_rule ( quad_num );
%
%  Consider all quantities associated with a given ELEMENT. 
%
  for element = 1 : element_num
%
%  Make a copy of the triangle.
%
    t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element));
    t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element));
%
%  Map the quadrature points QUADXY to points XY in the physical triangle.
%
    xy(1:2,1:quad_num) = reference_to_physical_t6 ( t6, quad_num, quad_xy );
    area = abs ( triangle_area_2d ( t3 ) );
    w(1:quad_num) = area * quad_w(1:quad_num);
    iu(1:6) = node_u_variable(element_node(1:6,element));
    iv(1:6) = node_v_variable(element_node(1:6,element));
    ip(1:3) = node_p_variable(element_node(1:3,element)) - 2 * node_num;
%
%  Evaluate the basis functions at the quadrature points. 
%
    [ Psi, ~, ~ ] = basis_mn_t6 ( t6, quad_num, xy );
    [ ~, dphidx, dphidy ] = basis_mn_t3 ( t3, quad_num, xy );
%
%  Evaluate the forcing function.
%
    [ u_rhs,  v_rhs,  ~ ] = external_force ( xy, t,       2, Re );
    [ u_rhs2, v_rhs2, ~ ] = external_force ( xy, t+.5*dt, 2, Re ); 
    [ u_rhs3, v_rhs3, ~ ] = external_force ( xy, t+dt,    2, Re );
%
%  Weak form.
%
    for i = 1 : 6

      F1(iu(i)) = F1(iu(i)) + sum ( w .* u_rhs(1:quad_num) .* Psi(i,1:quad_num));
      F1(iv(i)) = F1(iv(i)) + sum ( w .* v_rhs(1:quad_num) .* Psi(i,1:quad_num));
      F2(iu(i)) = F2(iu(i)) + sum ( w .* u_rhs2(1:quad_num) .* Psi(i,1:quad_num));
      F2(iv(i)) = F2(iv(i)) + sum ( w .* v_rhs2(1:quad_num) .* Psi(i,1:quad_num));

      if ( i < 4 )
        F3(ip(i)) = F3(ip(i)) + sum ( w .* ( u_rhs3(1:quad_num) ... 
          .* dphidx(i,1:quad_num) + v_rhs3(1:quad_num) ... 
          .* dphidy(i,1:quad_num)));
      end

    end

  end

  return 
end
function [ quad_w, quad_xy ] = quad_rule ( quad_num )

%*****************************************************************************80
%
%% QUAD_RULE sets the quadrature rule for assembly.
%
%  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:
%
%    17 July 2005
%
%  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_XY(2,QUAD_NUM), the quadrature nodes.
%
  if ( quad_num == 1 )

    quad_xy(1:2,1:quad_num) = [ 1.0 / 3.0, 1.0 / 3.0 ]';

    quad_w(1:quad_num) = 1.0;

  elseif ( quad_num == 3 )

    quad_xy(1:2,1:quad_num) = [ ...
      0.5, 0.0; ...
      0.5, 0.5; ...
      0.0, 0.5 ]';

    quad_w(1:quad_num) = 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_xy(1:2,1:quad_num) = [ ...
      b, b; ...
      c, a; ...
      a, c; ...
      a, a ]' / f;

    quad_w(1:quad_num) = [ 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_xy(1:2,1:quad_num) = [
      a, b; ...
      b, a; ...
      b, b; ...
      c, d; ...
      d, c; ...
      d, d ]';

    quad_w(1:6) = [ 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_xy(1:2,1:quad_num) = [ ...
      a, a; ...
      b, c; ...
      c, b; ...
      c, c; ...
      d, e; ...
      e, d; ...
      e, e ]';

    quad_w(1:quad_num) = [ 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_xy(1:2,1:quad_num) = [ ...
      a, b; ...
      b, a; ...
      b, b; ...
      c, d; ...
      c, e; ...
      d, c; ...
      d, e; ...
      e, c; ...
      e, d ]';

    quad_w(1:quad_num) = [ 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_xy(1:2,1:quad_num) = [
      h, h; ...
      a, b; ...
      b, a; ...
      b, b; ...
      c, d; ...
      d, c; ...
      d, d; ...
      e, f; ...
      e, g; ...
      f, e; ...
      f, g; ...
      g, e; ...
      g, f ]';

    quad_w(1:quad_num) = [ w, t, t, t, u, u, u, v, v, v, v, v, v ];

  else

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

  end

  return
end
function phy = reference_to_physical_t6 ( t, n, ref )

%*****************************************************************************80
%
%% REFERENCE_TO_PHYSICAL_T6 maps T6 reference points to physical points.
%
%  Discussion:
%
%    Given the vertices of an order 6 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.
%
%    The mapping from (XSI,ETA) to (X,Y) has the form:
%
%      X(ETA,XSI) = A1 * XSI**2 + B1 * XSI*ETA + C1 * ETA**2
%                 + D1 * XSI    + E1 * ETA     + F1
%
%      Y(ETA,XSI) = A2 * XSI**2 + B2 * XSI*ETA + C2 * ETA**2
%                 + D2 * XSI    + E2 * ETA     + F2
%
%  Reference Element T6:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,6), the coordinates of the vertices.
%    The vertices are assumed to be the images of (0,0), (1,0),
%    (0,1),(1/2,0), (1/2,1/2) and (0,1/2) respectively.
%
%    Input, integer N, the number of objects to transform.
%
%    Input, real REF(2,N), points in the reference triangle.
%
%    Output, real PHY(2,N), corresponding points in the
%    physical triangle.
%
  for i = 1 : 2

    a(i) =   2.0 * t(i,1) + 2.0 * t(i,2)                ...
           - 4.0 * t(i,4);

    b(i) =   4.0 * t(i,1)                               ...
           - 4.0 * t(i,4) + 4.0 * t(i,5) - 4.0 * t(i,6);

    c(i) =   2.0 * t(i,1)                + 2.0 * t(i,3) ...
                                         - 4.0 * t(i,6);

    d(i) = - 3.0 * t(i,1) -       t(i,2)                ...
           + 4.0 * t(i,4);

    e(i) = - 3.0 * t(i,1)                -       t(i,3) ...
                                         + 4.0 * t(i,6);
    f(i) =         t(i,1);

  end

  for i = 1 : 2
    phy(i,1:n) = a(i) * ref(1,1:n) .* ref(1,1:n) ...
               + b(i) * ref(1,1:n) .* ref(2,1:n) ...
               + c(i) * ref(2,1:n) .* ref(2,1:n) ...
               + d(i) * ref(1,1:n) ...
               + e(i) * ref(2,1:n) ...
               + f(i);
  end

  return
end
function area = triangle_area_2d ( t )

%*****************************************************************************80
%
%% TRIANGLE_AREA_2D computes the area of a triangle in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Output, real AREA, the absolute area of the triangle.
%
  area = 0.5 * abs ( ...
      t(1,1) * ( t(2,2) - t(2,3) ) ...
    + t(1,2) * ( t(2,3) - t(2,1) ) ...
    + t(1,3) * ( t(2,1) - t(2,2) ) );

  return
end
function [ phi, dphidx, dphidy ] = basis_mn_t6 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_T6: all bases for N points in a T6 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices and midside
%    nodes of a triangle.  It works directly with these coordinates, and does
%    not refer to a reference element.
%
%    This routine requires that the midside nodes be 'in line'
%    with the vertices, that is, that the sides of the triangle be
%    straight.  However, the midside nodes do not actually have to
%    be halfway along the side of the triangle.
%
%  Physical element T6:
%
%    This picture indicates the assumed ordering of the six nodes
%    of the triangle.
%
%             3
%            / \
%           /   \
%          6     5
%         /       \
%        /         \
%       1-----4-----2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,6), the nodal oordinates of the element.
%    It is common to list these points in counter clockwise order.
%
%    Input, real P(2,N), the evaluation points.
%
%    Output, real PHI(6,N), the basis functions at the evaluation points.
%
%    Output, real DPHIDX(6,N), DPHIDY(6,N), the basis derivatives at the
%    evaluation points.
%
%  Local Parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%

%
%  Basis function 1: PHI(X,Y) = G(3,2) * H(6,4) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) ...
          - ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) ) ...
          - ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,4) ) * ( t(2,6)   - t(2,4) ) ...
          - ( t(1,6)   - t(1,4) ) * ( p(2,1:n) - t(2,4) );

  hn(1:n) = ( t(1,1) - t(1,4) ) * ( t(2,6) - t(2,4) ) ...
          - ( t(1,6) - t(1,4) ) * ( t(2,1) - t(2,4) );

  phi(1,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(1,1:n) =  ( ( t(2,3) - t(2,2) ) * hx(1:n) ...
                    + gx(1:n) * ( t(2,6) - t(2,4) ) ) ...
                  ./ ( gn(1:n) .* hn(1:n) );

  dphidy(1,1:n) = -( ( t(1,3) - t(1,2) ) * hx(1:n) ...
                   + gx * ( t(1,6) - t(1,4) ) ) ...
                  ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 2: PHI(X,Y) = G(3,1) * H(4,5) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) ...
          - ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );

  gn(1:n) = ( t(1,2) - t(1,1) ) * ( t(2,3) - t(2,1) ) ...
          - ( t(1,3) - t(1,1) ) * ( t(2,2) - t(2,1) );

  hx(1:n) = ( p(1,1:n) - t(1,5) ) * ( t(2,4)   - t(2,5) ) ...
          - ( t(1,4)   - t(1,5) ) * ( p(2,1:n) - t(2,5) );

  hn(1:n) = ( t(1,2) - t(1,5) ) * ( t(2,4) - t(2,5) ) ...
          - ( t(1,4) - t(1,5) ) * ( t(2,2) - t(2,5) );

  phi(2,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(2,1:n) =  ( ( t(2,3) - t(2,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,4) - t(2,5) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(2,1:n) = -( ( t(1,3) - t(1,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,4) - t(1,5) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 3: PHI(X,Y) = G(1,2) * H(5,6) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) ...
          - ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) ) ...
          - ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,6) ) * ( t(2,5)   - t(2,6) ) ...
          - ( t(1,5)   - t(1,6) ) * ( p(2,1:n) - t(2,6) );

  hn(1:n) = ( t(1,3) - t(1,6) ) * ( t(2,5) - t(2,6) ) ...
          - ( t(1,5) - t(1,6) ) * ( t(2,3) - t(2,6) );

  phi(3,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(3,1:n) =  ( ( t(2,1) - t(2,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,5) - t(2,6) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(3,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,5) - t(1,6) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 4: PHI(X,Y) = G(1,3) * H(2,3) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,1)   - t(2,3) ) ...
          - ( t(1,1)   - t(1,3) ) * ( p(2,1:n) - t(2,3) );

  gn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,1) - t(2,3) ) ...
          - ( t(1,1) - t(1,3) ) * ( t(2,4) - t(2,3) );

  hx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,2)   - t(2,3) ) ...
          - ( t(1,2)   - t(1,3) ) * ( p(2,1:n) - t(2,3) );

  hn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,2) - t(2,3) ) ...
          - ( t(1,2) - t(1,3) ) * ( t(2,4) - t(2,3) );

  phi(4,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(4,1:n) =  ( ( t(2,1) - t(2,3) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,2) - t(2,3) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(4,1:n) = -( ( t(1,1) - t(1,3) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,2) - t(1,3) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 5: PHI(X,Y) = G(2,1) * H(3,1) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,2)   - t(2,1) ) ...
          - ( t(1,2)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );;

  gn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,2) - t(2,1) ) ...
          - ( t(1,2) - t(1,1) ) * ( t(2,5) - t(2,1) );

  hx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3)   - t(2,1) ) ...
          - ( t(1,3)   - t(1,1) ) * ( p(2,1:n) - t(2,1) );

  hn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,3) - t(2,1) ) ...
          - ( t(1,3) - t(1,1) ) * ( t(2,5) - t(2,1) );

  phi(5,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(5,1:n) =  ( ( t(2,2) - t(2,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,3) - t(2,1) ) ) ...
               ./ ( gn(1:n) .* hn(1:n) );

  dphidy(5,1:n) = -( ( t(1,2) - t(1,1) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,3) - t(1,1) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );
%
%  Basis function 6: PHI(X,Y) = G(1,2) * H(3,2) / normalization.
%
  gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1)   - t(2,2) ) ...
          - ( t(1,1)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  gn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,1) - t(2,2) ) ...
          - ( t(1,1) - t(1,2) ) * ( t(2,6) - t(2,2) );

  hx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3)   - t(2,2) ) ...
          - ( t(1,3)   - t(1,2) ) * ( p(2,1:n) - t(2,2) );

  hn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,3) - t(2,2) ) ...
          - ( t(1,3) - t(1,2) ) * ( t(2,6) - t(2,2) );

  phi(6,1:n) = ( gx(1:n) .* hx(1:n) ) ./ ( gn(1:n) .* hn(1:n) );

  dphidx(6,1:n) =  ( ( t(2,1) - t(2,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(2,3) - t(2,2) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  dphidy(6,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) ...
              + gx(1:n) * ( t(1,3) - t(1,2) ) ) ...
              ./ ( gn(1:n) .* hn(1:n) );

  return
end
function [ phi, dphidx, dphidy ] = basis_mn_t3 ( t, n, p )

%*****************************************************************************80
%
%% BASIS_MN_T3: all bases functions at N points for a T3 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices of a triangle.
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%    The sides of the triangle DO NOT have to lie along a coordinate
%    axis.
%
%    The routine evaluates the basis functions associated with each vertex,
%    and their derivatives with respect to X and Y.
%
%  Physical Element T3:
%
%            3
%           / \
%          /   \
%         /     \
%        /       \
%       1---------2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the vertices of the triangle.  It is common to list
%    these points in counter clockwise order.
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real P(2,N), the coordinates of the evaluation points.
%
%    Output, real PHI(3,N), the basis functions at the evaluation points.
%
%    Output, real DPHIDX(3,N), DPHIDY(3,N), the basis derivatives
%    at the evaluation points.
%
%  Local parameters:
%
%    Local, real AREA, is (twice) the area of the triangle.
%
  area = t(1,1) * ( t(2,2) - t(2,3) ) ...
       + t(1,2) * ( t(2,3) - t(2,1) ) ...
       + t(1,3) * ( t(2,1) - t(2,2) );

  phi(1,1:n) =     (  ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) )     ...
                    - ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) );
  dphidx(1,1:n) =   - ( t(2,3) - t(2,2) );
  dphidy(1,1:n) =     ( t(1,3) - t(1,2) );

  phi(2,1:n) =     (  ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) )     ...
                    - ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) );
  dphidx(2,1:n) =   - ( t(2,1) - t(2,3) );
  dphidy(2,1:n) =     ( t(1,1) - t(1,3) );

  phi(3,1:n) =     (  ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) )     ...
                    - ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) );
  dphidx(3,1:n) =   - ( t(2,2) - t(2,1) );
  dphidy(3,1:n) =     ( t(1,2) - t(1,1) );
%
%  Normalize.
%
  phi(1:3,1:n)    = phi(1:3,1:n) / area;
  dphidx(1:3,1:n) = dphidx(1:3,1:n) / area;
  dphidy(1:3,1:n) = dphidy(1:3,1:n) / area;

  return
end
function plotdata ( u, v, p, node_xy, node_num, pres, steps, dt, Re )

%*****************************************************************************80
%
%% PLOTDATA plots the results.
%
%  Modified:
%
%    03 December 2018
%

%
%  Plot the velocity magnitude.
%
  figure ( 1 );
  xlin = linspace ( min ( node_xy(1,:) ), max ( node_xy(1,:) ), ...
    min ( 500, ceil ( .5 * node_num ) ) );
  ylin = linspace ( min ( node_xy(2,:) ), max ( node_xy(2,:) ), ...
    min ( 500, ceil ( .5 * node_num ) ) ); 
  [ xx, yy ] = meshgrid ( xlin, ylin );
  vel = sqrt ( u.^2 + v.^2 );
  F = TriScatteredInterp ( node_xy(1,:)', node_xy(2,:)', vel ); 
  z = F(xx,yy);
  colormap ( jet );
  imagesc ( xlin, ylin, z );
  hold on;
%
%  Plot the streamlines.
%
  U = TriScatteredInterp ( node_xy(1,:)', node_xy(2,:)', u );
  V = TriScatteredInterp ( node_xy(1,:)', node_xy(2,:)', v );
  uu = U(xx,yy);
  vv = V(xx,yy);
  U = streamslice ( xx, yy, uu, vv );
  axis ( 'tight' );
  set ( U, 'Color', 'white' )
  xlabel ( 'x', 'fontsize', 10, 'fontweight', 'b' );
  ylabel ( 'y', 'fontsize', 10, 'fontweight', 'b' );
  set ( gca, 'XTick', [0 0.2 0.4 0.6 0.8 1], ...
    'YTick', [0 0.2 0.4 0.6 0.8 1], 'fontsize', 10, 'fontweight', 'b', ...
    'ydir', 'normal' );
  hold off;
  title ( 'Velocity magnitude', 'Fontsize', 16 );
  print ( '-dpng', 'box_flow_velocity_magnitude.png' );
%
%  Plot the pressure.
%
  figure ( 2 );
  F = TriScatteredInterp ( node_xy(1,pres)', node_xy(2,pres)', p );
  pp = F(xx,yy);
  surf ( xx, yy, pp, 'FaceColor', 'interp', 'FaceLighting', 'phong' );
  xlabel ( 'x', 'fontsize', 10, 'fontweight', 'b' );
  ylabel ( 'y', 'fontsize', 10, 'fontweight', 'b' );
  zlabel ( 'pressure', 'fontsize', 10, 'fontweight', 'b' );
  set ( gca, 'XTick', [0 0.2 0.4 0.6 0.8 1], ...
    'YTick', [0 0.2 0.4 0.6 0.8 1], 'fontsize', 10, 'fontweight', 'b' );
  title ( 'Pressure', 'Fontsize', 16 );
  print ( '-dpng', 'box_flow_pressure.png' );
%
%  Plot the velocity error.
%
  [ ua, va, pa] = external_force ( node_xy, ( steps ) * dt, 1, Re );

  fprintf ( 1, '  Loo norm of UA - U = %g\n', norm ( ua' - u, inf ) );
  fprintf ( 1, '  L2  norm of UA - U = %g\n', norm ( ua' - u, 2 ) );

  xlin = linspace ( min ( node_xy(1,:) ), max ( node_xy(1,:) ), ...
    min ( 50, ceil ( .5 * node_num) ) );
  ylin = linspace ( min ( node_xy(2,:) ), max ( node_xy(2,:) ), ...
    min ( 50, ceil ( .5 * node_num) ) );
  [ xx, yy ] = meshgrid ( xlin, ylin );

  figure ( 3 );
  true_v = sqrt ( ua .^ 2 + va .^ 2 )';
  F = TriScatteredInterp ( node_xy(1,:)', node_xy(2,:)', ( vel - true_v ) );
  z = F(xx,yy);
  colormap ( jet );
  surf ( xx, yy, z, 'FaceColor', 'interp', 'EdgeColor', 'none', ...
    'FaceLighting', 'phong ' ); 
  title ( 'Velocity error', 'Fontsize', 16 );
  colorbar;
  print ( '-dpng', 'box_flow_velocity_error.png' );
%
%  Plot the pressure error.
%
  figure ( 4 );
  F = TriScatteredInterp ( node_xy(1,pres)', node_xy(2,pres)', p - pa(pres)' );
  pp = F(xx,yy);
  colormap( jet );
  surf ( xx, yy, pp, 'FaceColor', 'interp', 'EdgeColor', 'none', ...
    'FaceLighting', 'phong' );
  title ( 'Pressure error', 'Fontsize', 16 );
  colorbar;
  print ( '-dpng', 'box_flow_pressure_error.png' );

  return
end
function bnd = triangulation_order6_boundary_node ( num_node, xy, t )

%*****************************************************************************80
%
%% TRIANGULATION ORDER6 BOUNDARY NODE indicates which nodes are on the boundary.
%
%  p contains coordinates
%  t contains elements
%
%  Modified:
%
%    03 December 2018
%
  t = t';

  t2 = [ ...
    t(:,4); ...
    t(:,5); ...
    t(:,6) ];

  t = [ ...
    t(:,1:2); ...
    t(:,2:3); ...
    [ t(:,3), t(:,1) ] ];

  t = sort ( t, 2 );
  t = [ t, t2 ];
  edges = unique ( t, 'rows' );
  tl = length ( edges );
%
%  Identify boundary nodes
%
  bnd = zeros ( 1, num_node );
  for i = 1 : tl
    indx2 = find ( t(:,1) == edges(i,1) & t(:,2) == edges(i,2) ); 
    if ( length(indx2) == 1 )
      bnd(1,edges(i,1)) = 1;
      bnd(1,edges(i,2)) = 1;
      bnd(1,edges(i,3)) = 1;
    end
  end
%
%  Interior boundaries on rectangle domain.
%
  max_x = max ( xy(:,1) ) - eps;
  min_x = min ( xy(:,1) ) + eps;
  max_y = max ( xy(:,2) ) - eps;
  min_y = min ( xy(:,2) ) + eps;

  for i = 1 : num_node
    if ( bnd(i) == 1 )
      if ( xy(i,1) < max_x & xy(i,1) > min_x & ...
        xy(i,2) < max_y & xy(i,2) > min_y )
        bnd(i) = 2;
      end
    end
  end

  return 
end
function del_x = min_dist ( coords, nodes ) 

%*****************************************************************************80
%
%% MIN_DIST returns a minimum distance between nodes.
%
%  Modified:
%
%    03 December 2018
%
  [ rows cols ] = size ( nodes );

  del_x = 9999;

  for element = 1 : 1 : cols
    side1 = norm ( coords(:,nodes(1,element)) - coords(:,nodes(4,element)) );
    side3 = norm ( coords(:,nodes(1,element)) - coords(:,nodes(6,element)) );
    side2 = norm ( coords(:,nodes(2,element)) - coords(:,nodes(5,element)) );
    del_x = min ( [ del_x, side1, side2, side3 ] );
  end

  return
end

