MATLAB Examples

Contents

2D Unsteady Navier-Stokes problem

Based on "Finite Element Methods for flow problems" of Jean Donea and Antonio Huerta

% Andrea La Spina
% https://it.linkedin.com/in/andrealaspina
% https://independent.academia.edu/AndreaLaSpina
% https://www.researchgate.net/profile/Andrea_La_Spina

Equation

% u_t-v*div(grad(u))+(u.div)u+grad(p)=b   in W
% div(u)=0                                in W
% u=u_D                                   on W_D

Initialization

initialization

Input

% Model parameters --------------------------------------------------------
v=1;                                    % Kinematic viscosity
b_fun=@(x,y) ...                        % Body force
    [(12-24*y)*x^4+(-24+48*y)*x^3+(-48*y+72*y^2-48*y^3+12)*x^2+...
     (-2+24*y-72*y^2+48*y^3)*x+(1-4*y+12*y^2-8*y^3),...
     (8-48*y+48*y^2)*x^3+(-12+72*y-72*y^2)*x^2+...
     (4-24*y+48*y^2-48*y^3+24*y^4)*x+(-12*y^2+24*y^3-12*y^4)];
% -------------------------------------------------------------------------

% Space -------------------------------------------------------------------
x_i=0;                                  % Initial point (x)
x_f=1;                                  % Final point (x)
y_i=0;                                  % Initial point (y)
y_f=1;                                  % Final point (y)
% -------------------------------------------------------------------------

% Time --------------------------------------------------------------------
t_i=0;                                  % Initial time
t_f=0.1;                                % Final time
dt=t_f/10;                              % Time step
animation_time=10;                      % Animation time
theta=1;                                % Theta for time integration scheme
                                        % ( 0  = Forward Euler)
                                        % (1/2 = Crank-Nicolson)
                                        % (2/3 = Galerkin)
                                        % ( 1  = Backward Euler)
% -------------------------------------------------------------------------

% FEM ---------------------------------------------------------------------
n_el_x=10;                              % Number of finite elements (x)
n_el_y=10;                              % Number of finite elements (y)
n_gauss_side=3;                         % Number of Gauss points per side
body_force_selection='gauss';           % Selection points of body force
                                        % ('centre','gauss','nodes_v')
toll_p=1e-6;                            % Tollerance for pressure
toll_v=1e-6;                            % Tollerance for velocity
% -------------------------------------------------------------------------

% Boundary conditions -----------------------------------------------------
bound_cond_p=0;                         % Boundary conditions (pressure)
bound_cond_v=0;                         % Boundary conditions (velocity)
dof_constrained_string_p=...            % DOFs constrained (pressure)
    ['[1]'];
dof_constrained_string_v=...            % DOFs constrained (velocity)
    ['[2:n_np_x_v-1,',...
     'n_np_v-n_np_x_v+2:n_np_v-1,'...
     '1:n_np_x_v:n_np_v-n_np_x_v+1,'...
     'n_np_x_v:n_np_x_v:n_np_v]'];
% -------------------------------------------------------------------------

% Initial conditions ------------------------------------------------------
radius=0.2;                             % Radius of the hill
x_0=1/2;                                % Centre of the hill (x)
y_0=1/2;                                % Centre of the hill (y)
u_0_fun=@(x,y) [0,0];                   % Initial condition (velocity)

% -------------------------------------------------------------------------

% Plots -------------------------------------------------------------------
plot_geometry='yes';                    % Plot geometry ('yes','no')
plot_shape_functions='yes';             % Plot shape functions ('yes','no')
plot_test_functions='yes';              % Plot test functions ('yes','no')
% -------------------------------------------------------------------------

Derived parameters

% Common parameters -------------------------------------------------------
n_sd=2;                                 % Number of dimensional space
L_x=x_f-x_i;                            % Domain length (x)
L_y=y_f-y_i;                            % Domain length (y)
n_el=n_el_x*n_el_y;                     % Number of finite elements
n_gauss=n_gauss_side^2;                 % Number of Gauss points
L_el_x=L_x/n_el_x;                      % Length of a finite element (x)
L_el_y=L_y/n_el_y;                      % Length of a finite element (y)
% -------------------------------------------------------------------------

% Pressure ----------------------------------------------------------------
n_np_x_p=n_el_x+1;                      % Number of nodal points (x)
n_np_y_p=n_el_y+1;                      % Number of nodal points (y)
n_np_p=n_np_x_p*n_np_y_p;               % Number of nodal points
dof_el_p=4;                             % Number of DOFs per element
dof_p=n_np_p;                           % Total number of DOFs
h_x_p=L_el_x;                           % Spatial step (x)
h_y_p=L_el_y;                           % Spatial step (y)
x_p=x_i:h_x_p:x_f;                      % Space vector (x)
y_p=y_i:h_y_p:y_f;                      % Space vector (y)
% -------------------------------------------------------------------------

% Velocity ----------------------------------------------------------------
n_np_x_v=2*n_el_x+1;                    % Number of nodal points (x)
n_np_y_v=2*n_el_y+1;                    % Number of nodal points (y)
n_np_v=n_np_x_v*n_np_y_v;               % Number of nodal points
dof_el_v=9;                             % Number of DOFs per element
dof_v=n_np_v;                           % Total number of DOFs
h_x_v=L_el_x/2;                         % Spatial step (x)
h_y_v=L_el_y/2;                         % Spatial step (y)
x_v=x_i:h_x_v:x_f;                      % Space vector (x)
y_v=y_i:h_y_v:y_f;                      % Space vector (y)
% -------------------------------------------------------------------------

% Body force --------------------------------------------------------------
n_np_x_b=n_el_x;                        % Number of nodal points (x)
n_np_y_b=n_el_y;                        % Number of nodal points (y)
n_np_b=n_np_x_b*n_np_y_b;               % Number of nodal points
h_x_b=L_el_x;                           % Spatial step (x)
h_y_b=L_el_y;                           % Spatial step (y)
x_b=x_i+L_el_x/2:h_x_b:x_f-L_el_x/2;    % Space vector (x)
y_b=y_i+L_el_y/2:h_y_b:y_f-L_el_y/2;    % Space vector (y)
% -------------------------------------------------------------------------

% Time --------------------------------------------------------------------
T=t_i:dt:t_f;                           % Time vector
% -------------------------------------------------------------------------

Gauss parameters

gauss=[];
[gauss]=Gauss_parameters_2D(n_gauss_side,gauss);

% Trasformation of coordinated for the Gauss integration points
for n=1:n_gauss
    x_gauss(n)=L_el_x/2*(1+gauss(n).csi);
    y_gauss(n)=L_el_y/2*(1+gauss(n).eta);
end

% Jacobian of the transformation
J_mat=[L_el_x/2    0
          0    L_el_y/2];
J=det(J_mat);

% Computation of shape and test functions (and derivatives) at Gauss points
% for pressure
[gauss]=shape_functions_Gauss_points_2D_p(gauss);
[gauss]=test_functions_Gauss_points_2D_p(gauss);

% Computation of shape and test functions (and derivatives) at Gauss points
% for velocity
[gauss]=shape_functions_Gauss_points_2D_v(gauss);
[gauss]=test_functions_Gauss_points_2D_v(gauss);

Plot of geometry

if strcmp(plot_geometry,'yes')==1
    if n_el<=100
        plot_geometry_2D(x_p,y_p,x_v,y_v,L_x,L_y,n_gauss,...
                                            L_el_x,L_el_y,x_gauss,y_gauss);
    end
end

Plot of shape and test functions

% Normalized domain
csi_plot=-1:0.1:1;
eta_plot=-1:0.1:1;

if strcmp(plot_shape_functions,'yes')==1
    % Shape functions for pressure
    N_plot_p=f_N_plot_2D_p(csi_plot,eta_plot);
    plot_shape_functions_2D(csi_plot,eta_plot,N_plot_p,dof_el_p);

    % Shape functions for velocity
    N_plot_v=f_N_plot_2D_v(csi_plot,eta_plot);
    plot_shape_functions_2D(csi_plot,eta_plot,N_plot_v,dof_el_v);
end

if strcmp(plot_test_functions,'yes')==1
    % Test functions for pressure
    W_plot_p=f_W_plot_2D_p(csi_plot,eta_plot);
    plot_test_functions_2D(csi_plot,eta_plot,W_plot_p,dof_el_p);

    % Test functions for velocity
    W_plot_v=f_W_plot_2D_v(csi_plot,eta_plot);
    plot_test_functions_2D(csi_plot,eta_plot,W_plot_v,dof_el_v);
end

Evaluate matrices and vectors

% Afference matrix for pressure
[A_p]=afference_matrix_2D_p(n_np_x_p,n_np_y_p,dof_el_p);

% Afference matrix for velocity
[A_v]=afference_matrix_2D_v(n_np_x_v,n_np_y_v,dof_el_v);

% Element mass matrix
for n=1:n_el
    el(n).M=element_mass_matrix_2D(dof_el_v,gauss,J);
end

% Strain rate - velocity matrix
[gauss]=strain_rate_velocity_matrix_2D(dof_el_v,gauss,L_el_x,L_el_y);

% Element viscosity matrix
for n=1:n_el
    el(n).K=element_viscosity_matrix_2D(v,dof_el_v,gauss,J);
end

% Element gradient operator matrix
for n=1:n_el
    el(n).G=element_gradient_operator_matrix_2D(dof_el_v,dof_el_p,gauss,...
                                                          J,L_el_x,L_el_y);
end

% Element load vector for velocity
for n=1:n_el
    [r,c]=row_column(n,n_el_x);
    x_c=x_i+L_el_x/2+L_el_x*(c-1);
    y_c=y_i+L_el_y*(n_el_y-1/2)-L_el_y*(r-1);
    if strcmp(body_force_selection,'centre')==1
        b_c=feval(b_fun,x_c,y_c);
        el(n).f=element_load_vector_centre_2D(b_c,dof_el_v,gauss,J);
    elseif strcmp(body_force_selection,'gauss')==1
        for k=1:n_gauss
            x_g(k)=x_c-L_el_x/2+x_gauss(k);
            y_g(k)=y_c-L_el_y/2+y_gauss(k);
            b_g(k,:)=feval(b_fun,x_g(k),y_g(k));
        end
        el(n).f=element_load_vector_gauss_2D(b_g,dof_el_v,gauss,J);
    elseif strcmp(body_force_selection,'nodes_v')==1
        x_n(1)=x_c-L_el_x/2; y_n(1)=y_c-L_el_y/2;
        x_n(2)=x_c+L_el_x/2; y_n(2)=y_c-L_el_y/2;
        x_n(3)=x_c+L_el_x/2; y_n(3)=y_c+L_el_y/2;
        x_n(4)=x_c-L_el_x/2; y_n(4)=y_c+L_el_y/2;
        x_n(5)=x_c;          y_n(5)=y_c-L_el_y/2;
        x_n(6)=x_c+L_el_x/2; y_n(6)=y_c;
        x_n(7)=x_c;          y_n(7)=y_c+L_el_y/2;
        x_n(8)=x_c-L_el_x/2; y_n(8)=y_c;
        x_n(9)=x_c;          y_n(9)=y_c;
        for k=1:dof_el_v
            b_n(k,:)=feval(b_fun,x_n(k),y_n(k));
        end
        el(n).f=element_load_vector_nodes_v_2D(b_n,dof_el_v,gauss,J);
    end
end

% Element load vector for pressure
for n=1:n_el
    el(n).h=zeros(dof_el_p,1);
end

Assemblate matrices and vectors

% Assemblage of mass matrix
[M]=assemble_mass_matrix(el,dof_v,n_el,dof_el_v,A_v);

% Assemblage of viscosity matrix
[K]=assemble_viscosity_matrix(el,dof_v,n_el,dof_el_v,A_v);

% Assemblage of gradient operator matrix
[G]=assemble_gradient_operator_matrix(el,dof_v,dof_p,n_el,dof_el_v,...
                                                         dof_el_p,A_v,A_p);

% Divergence operator matrix
[G_T]=G';

% Zero matrix
ZERO=zeros(dof_p,dof_p);

% Zero_1 matrix
ZERO_1=zeros(dof_v*n_sd,dof_p);

% Zero_2 matrix
ZERO_2=zeros(dof_p,dof_v*n_sd);

% Assemblage of load vector for velocity
[f]=assemble_load_vector_v(el,dof_v,n_el,dof_el_v,A_v);

% Assemblage of load vector for pressure
[h]=assemble_load_vector_p(el,dof_p,n_el,dof_el_p,A_p);

Conversion of the matrices from full to sparse

M=sparse(M);
K=sparse(K);
G=sparse(G);
G_T=sparse(G_T);
f=sparse(f);
h=sparse(h);
ZERO=sparse(ZERO);
ZERO_1=sparse(ZERO_1);
ZERO_2=sparse(ZERO_2);

Boundary conditions

% Definition of the constrained DOFs for the pressure
dof_constrained_p=eval(dof_constrained_string_p);
dof_free_p=dof_p-length(dof_constrained_p);
dof_constrained_p=sort(dof_constrained_p);
n_dof_constrained_p=length(dof_constrained_p);

% Definition of the constrained DOFs for the velocity
dof_constrained_v=eval(dof_constrained_string_v);
dof_free_v=dof_v-length(dof_constrained_v);
dof_constrained_v=sort(dof_constrained_v);
n_dof_constrained_v=length(dof_constrained_v);
dof_constrained_v_X_Y=sort([dof_constrained_v.*2-1,dof_constrained_v.*2]);
dof_free_v_X_Y=dof_v*n_sd-length(dof_constrained_v_X_Y);
n_dof_constrained_v_X_Y=length(dof_constrained_v_X_Y);

% Evaluation of the derivative of the boundary conditions for the pressure
bound_cond_p_der=0;

% Evaluation of the derivative of the boundary conditions for the velocity
bound_cond_v_der=0;

% Evaluation of boundary conditions over time for pressure
for m=1:length(T)
    time(m).t=T(m);
    for n=1:n_dof_constrained_p
       constrain_p(n)=bound_cond_p;
       constrain_p_der(n)=bound_cond_p_der;
    end
    time(m).p_p(:,1)=constrain_p';
    time(m).p_der_p(:,1)=constrain_p_der';
end

% Evaluation of boundary conditions over time for velocity
for m=1:length(T)
    time(m).t=T(m);
    for n=1:n_dof_constrained_v_X_Y
       constrain_v(n)=bound_cond_v;
       constrain_v_der(n)=bound_cond_v_der;
    end
    time(m).u_p(:,1)=constrain_v';
    time(m).u_der_p(:,1)=constrain_v_der';
end

% Mass matrix
[M_ff,M_fp,M_pf,M_pp]=constrain_matrix(M,...
                              dof_constrained_v_X_Y,dof_constrained_v_X_Y);

% Viscosity matrix
[K_ff,K_fp,K_pf,K_pp]=constrain_matrix(K,...
                              dof_constrained_v_X_Y,dof_constrained_v_X_Y);

% Gradient operator matrix
[G_ff,G_fp,G_pf,G_pp]=constrain_matrix(G,...
                                  dof_constrained_v_X_Y,dof_constrained_p);

% Divergence operator matrix
[G_T_ff,G_T_fp,G_T_pf,G_T_pp]=constrain_matrix(G_T,...
                                  dof_constrained_p,dof_constrained_v_X_Y);

% Zero matrix
[ZERO_ff,ZERO_fp,ZERO_pf,ZERO_pp]=constrain_matrix(ZERO,...
                                      dof_constrained_p,dof_constrained_p);

% Zero_1 matrix
[ZERO_1_ff,ZERO_1_fp,ZERO_1_pf,ZERO_1_pp]=constrain_matrix(ZERO_1,...
                                  dof_constrained_v_X_Y,dof_constrained_p);

% Zero_2 matrix
[ZERO_2_ff,ZERO_2_fp,ZERO_2_pf,ZERO_2_pp]=constrain_matrix(ZERO_2,...
                                  dof_constrained_p,dof_constrained_v_X_Y);

% Load vector for velocity
[f_f,f_p]=constrain_vector(f,dof_constrained_v_X_Y);

% Load vector for pressure
[h_f,h_p]=constrain_vector(h,dof_constrained_p);

Initial conditions

% Initial conditions for the velocity
time(1).iter(1).u=zeros(dof_v*n_sd,1,1);
for n=1:n_np_v
    [r,c]=row_column(n,n_np_x_v);
    xp=x_i+L_el_x/2*(c-1);
    yp=y_i+L_el_y*n_el_y-L_el_y/2*(r-1);
    u_0_xp_yp=feval(u_0_fun,xp,yp);
    time(1).iter(1).u(n_sd*n-1,1)=u_0_xp_yp(1);
    time(1).iter(1).u(n_sd*n,1)=u_0_xp_yp(2);
end
time(1).iter(1).u_f=constrain_vector(time(1).iter(1).u,...
                                                    dof_constrained_v_X_Y);

% Pressure from the static problem ----------------------------------------
K_ff_inv=inv(K_ff);
time(1).iter(1).p_f=(G_T_ff*K_ff_inv*G_ff)\...
                 (...
                 G_T_ff*K_ff_inv*(f_f-K_fp*time(1).u_p-G_fp*time(1).p_p)...
                 +G_T_fp*time(1).u_p-h_f...
                 );
[time(1).iter(1).p]=data_all_dof(time(1).iter(1).p_f,time(1).p_p,...
                                                        dof_constrained_p);
% -------------------------------------------------------------------------
time(1).iter(1).z=[time(1).iter(1).u_f;time(1).iter(1).p_f];

% Conversion of data from vector to matrix
for j=1:n_np_y_v
    [time(1).iter(1).u_x_matrix(:,j)]=time(1).iter(1).u(...
                                     ((n_np_y_v-j)*n_np_x_v+1)*2-1:2:...
                                     ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-1);
    [time(1).iter(1).u_y_matrix(:,j)]=time(1).iter(1).u(...
                                     ((n_np_y_v-j)*n_np_x_v+1)*2-0:2:...
                                     ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-0);
end
time(1).iter(1).u_matrix=sqrt(time(1).iter(1).u_x_matrix.^2+...
                              time(1).iter(1).u_y_matrix.^2);
for j=1:n_np_y_p
    [time(1).iter(1).p_matrix(:,j)]=time(1).iter(1).p(...
                                           (n_np_y_p-j)*n_np_x_p+1:...
                                           (n_np_y_p-j)*n_np_x_p+n_np_x_p);
end

Unteady Navier-Stokes problem

time(1).iter(1).D=K;
[time(1).iter(1).D_ff,time(1).iter(1).D_fp,...
 time(1).iter(1).D_pf,time(1).iter(1).D_pp]=...
 constrain_matrix(time(1).iter(1).D,...
                              dof_constrained_v_X_Y,dof_constrained_v_X_Y);
disp('TIME INTEGRATION')
disp('*******************************************************************')
time(1).n_iter=1;
tic
for m=1:length(T)-1

    % Initialization of variables
    time(m+1).iter(1).u=time(m).iter(end).u;
    time(m+1).iter(1).p=time(m).iter(end).p;

    % Initialization of unknowns
    [time(m+1).iter(1).u_f,time(m+1).iter(1).u_p]=...
               constrain_vector(time(m+1).iter(1).u,dof_constrained_v_X_Y);
    [time(m+1).iter(1).p_f,time(m+1).iter(1).p_p]=...
               constrain_vector(time(m+1).iter(1).p,dof_constrained_p);
    time(m+1).iter(1).z=[time(m+1).iter(1).u_f;time(m+1).iter(1).p_f];

    % Initialization of Viscosity+Convection matrix
    time(m+1).iter(1).D=time(m).iter(time(m).n_iter).D;
    [time(m+1).iter(1).D_ff,time(m+1).iter(1).D_fp,...
     time(m+1).iter(1).D_pf,time(m+1).iter(1).D_pp]=...
     constrain_matrix(time(m+1).iter(1).D,...
                              dof_constrained_v_X_Y,dof_constrained_v_X_Y);

    % Initialization of errors (>tollerances)
    time(m+1).iter(1).norm_err_p=toll_p+1;
    time(m+1).iter(1).norm_err_v=toll_v+1;

    % Convergence plots
    figure('Color',[1 1 1])
    axes('FontSize',14)
    subplot(1,2,1,'YScale','log')
    hold on
    title(['Time step ',num2str(m),'/',num2str(length(T)-1)],'FontSize',14)
    xlabel('Iteration','FontSize',14)
    ylabel('Norm((p_k-p_k_-_1)/p_k)','FontSize',14)
    grid on
    grid minor
    subplot(1,2,2,'YScale','log')
    hold on
    title(['Time step ',num2str(m),'/',num2str(length(T)-1)],'FontSize',14)
    xlabel('Iteration','FontSize',14)
    ylabel('Norm((u_k-u_k_-_1)/u_k)','FontSize',14)
    grid on
    grid minor
    pause(eps)

    fprintf('\nTIME STEP %d/%d\n',m,length(T)-1)
disp('-------------------------------------------------------------------')
    clear time_analysis time_inversion time_solution_p time_solution_v
    k=1;
    while time(m+1).iter(end).norm_err_p>toll_p ||...
          time(m+1).iter(end).norm_err_v>toll_v

        % Previous evaluation
        p_old=time(m+1).iter(k).p;
        u_old=time(m+1).iter(k).u;

        % Element convection matrix
        for n=1:n_el
            u_el=zeros(dof_el_v*n_sd,1);
            for a=1:dof_el_v
                for i=1:n_sd
                    r=n_sd*(a-1)+i;
                    u_el(r,1)=time(m+1).iter(k).u(A_v(n,r),1);
                end
            end
            el(n).C=element_convection_matrix_2D(u_el,dof_el_v,gauss,...
                                                          J,L_el_x,L_el_y);
        end

        % Assemblage of convection matrix
        [C]=assemble_convection_matrix(el,dof_v,n_el,dof_el_v,A_v);

        % Conversion of the convection matrix from full to sparse
        C=sparse(C);

        % Viscosity+Convection matrix
        time(m+1).iter(k).D=K+C;

        % Constrain viscosity+convection matrix
        [time(m+1).iter(k).D_ff,time(m+1).iter(k).D_fp,...
         time(m+1).iter(k).D_pf,time(m+1).iter(k).D_pp]=...
         constrain_matrix(time(m+1).iter(k).D,...
                              dof_constrained_v_X_Y,dof_constrained_v_X_Y);

        % -----------------------------------------------------------------
        % Definition of matrices
        H=[M_ff,ZERO_1_ff;ZERO_2_ff,ZERO_ff];
        Q2=[time(m+1).iter(k).D_ff,G_ff;G_T_ff,ZERO_ff];
        Q1=[time(m).iter(time(m).n_iter).D_ff,G_ff;G_T_ff,ZERO_ff];
        S2=[f_f;h_f]-...
           [M_fp,ZERO_1_fp;ZERO_2_fp,ZERO_fp]*...
           [time(m+1).u_der_p;time(m+1).p_der_p]-...
           [time(m+1).iter(k).D_fp,G_fp;G_T_fp,ZERO_fp]*...
           [time(m+1).u_p;time(m+1).p_p];
        S1=[f_f;h_f]-...
            [M_fp,ZERO_1_fp;ZERO_2_fp,ZERO_fp]*...
            [time(m).u_der_p;time(m).p_der_p]-...
            [time(m).iter(time(m).n_iter).D_fp,G_fp;G_T_fp,ZERO_fp]*...
            [time(m).u_p;time(m).p_p];

        % Evaluation of the new unknowns
        time(m+1).iter(k+1).z=(H+dt*  theta  *Q2)\...
                     (...
                     (H-dt*(1-theta)*Q1)*time(m).iter(time(m).n_iter).z+...
                        dt*  theta  *S2+...
                        dt*(1-theta)*S1...
                     );

        time(m+1).iter(k+1).u_f=...
                                 time(m+1).iter(k+1).z(1:dof_free_v_X_Y,1);
        time(m+1).iter(k+1).p_f=...
                             time(m+1).iter(k+1).z(dof_free_v_X_Y+1:end,1);
        % -----------------------------------------------------------------

        % Data for all dof
        [time(m+1).iter(k+1).u]=...
                     data_all_dof(time(m+1).iter(k+1).u_f,time(m+1).u_p,...
                                                    dof_constrained_v_X_Y);
        [time(m+1).iter(k+1).p]=...
                     data_all_dof(time(m+1).iter(k+1).p_f,time(m+1).p_p,...
                                                        dof_constrained_p);

        % Current evaluation
        u_new=time(m+1).iter(k+1).u;
        p_new=time(m+1).iter(k+1).p;

        % Evaluation of errors
        time(m+1).iter(k+1).norm_err_p=norm((p_new-p_old)/p_new);
        time(m+1).iter(k+1).norm_err_v=norm((u_new-u_old)/u_new);

        % Display results
        fprintf('\nIteration = %d',k-1)
        fprintf('\tNorm_p = %.2e',time(m+1).iter(k+1).norm_err_p)
        fprintf('\tNorm_v = %.2e',time(m+1).iter(k+1).norm_err_v)
        fprintf('\tElapsed time %.2f sec',toc)
        fprintf('\tReynolds = %.1f',...
                                max(time(m+1).iter(k+1).u)*mean(L_x,L_y)/v)

        % Convergence plot update
        subplot(1,2,1)
        semilogy(k-1,time(m+1).iter(k+1).norm_err_p,'k+','LineWidth',2)
        plot([0,k+1],[toll_p,toll_p],'r','LineWidth',2)
        xlim([0,k+1])
        ylim([min(toll_p/10,min(time(m+1).iter(k+1).norm_err_p)),...
              max(10,max(time(m+1).iter(k+1).norm_err_p)*10)])
        subplot(1,2,2)
        semilogy(k-1,time(m+1).iter(k+1).norm_err_v,'k+','LineWidth',2)
        plot([0,k+1],[toll_v,toll_v],'r','LineWidth',2)
        xlim([0,k+1])
        ylim([min(toll_v/10,min(time(m+1).iter(k+1).norm_err_v)),...
              max(10,max(time(m+1).iter(k+1).norm_err_v)*10)])
        pause(eps)

        k=k+1;
    end
    disp(' ')
disp('-------------------------------------------------------------------')

    time(m+1).n_iter=k-1;
end
disp(' ')
disp('*******************************************************************')
disp(' ')

% Conversion of data from vector to matrix
for m=1:length(T)
    for j=1:n_np_y_v
        [time(m).u_x_matrix(:,j)]=time(m).iter(1).u(...
                                     ((n_np_y_v-j)*n_np_x_v+1)*2-1:2:...
                                     ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-1);
        [time(m).u_y_matrix(:,j)]=time(m).iter(1).u(...
                                     ((n_np_y_v-j)*n_np_x_v+1)*2-0:2:...
                                     ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-0);
    end
    time(m).u_matrix=sqrt(time(m).u_x_matrix.^2+time(m).u_y_matrix.^2);

    for j=1:n_np_y_p
        [time(m).p_matrix(:,j)]=time(m).iter(1).p(...
                                           (n_np_y_p-j)*n_np_x_p+1:...
                                           (n_np_y_p-j)*n_np_x_p+n_np_x_p);
    end
end
TIME INTEGRATION
*******************************************************************

TIME STEP 1/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 1.60e+00	Norm_v = 1.32e+01	Elapsed time 1.90 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 1.43e-06	Norm_v = 1.11e-06	Elapsed time 3.54 sec	Reynolds = 0.0
Iteration = 2	Norm_p = 1.50e-13	Norm_v = 1.15e-12	Elapsed time 5.15 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 2/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 7.36e-01	Norm_v = 4.80e+00	Elapsed time 7.02 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 7.74e-08	Norm_v = 3.72e-07	Elapsed time 8.65 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 3/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 4.32e-01	Norm_v = 2.47e+00	Elapsed time 10.51 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 8.00e-08	Norm_v = 2.08e-07	Elapsed time 12.14 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 4/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 2.76e-01	Norm_v = 1.45e+00	Elapsed time 13.99 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 4.32e-08	Norm_v = 1.29e-07	Elapsed time 15.62 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 5/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 1.83e-01	Norm_v = 9.07e-01	Elapsed time 17.47 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 2.49e-08	Norm_v = 8.39e-08	Elapsed time 19.11 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 6/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 1.23e-01	Norm_v = 5.87e-01	Elapsed time 21.01 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 1.56e-08	Norm_v = 5.58e-08	Elapsed time 22.63 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 7/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 8.32e-02	Norm_v = 3.87e-01	Elapsed time 24.52 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 1.02e-08	Norm_v = 3.75e-08	Elapsed time 26.16 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 8/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 5.66e-02	Norm_v = 2.59e-01	Elapsed time 28.05 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 6.85e-09	Norm_v = 2.55e-08	Elapsed time 29.70 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 9/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 3.86e-02	Norm_v = 1.75e-01	Elapsed time 31.64 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 4.65e-09	Norm_v = 1.73e-08	Elapsed time 33.28 sec	Reynolds = 0.0 
-------------------------------------------------------------------

TIME STEP 10/10
-------------------------------------------------------------------

Iteration = 0	Norm_p = 2.64e-02	Norm_v = 1.19e-01	Elapsed time 35.68 sec	Reynolds = 0.0
Iteration = 1	Norm_p = 3.17e-09	Norm_v = 1.18e-08	Elapsed time 37.80 sec	Reynolds = 0.0 
-------------------------------------------------------------------
 
*******************************************************************
 

Plots

% Grid matrices
[X_v,Y_v]=meshgrid(x_v,y_v);
[X_p,Y_p]=meshgrid(x_p,y_p);
[X_b,Y_b]=meshgrid(x_b,y_b);

% Evaluation of the body force field
for n=1:n_np_b
    [r,c]=row_column(n,n_np_x_b);
    xp=x_i+L_el_x*(c-1)+L_el_x/2;
    yp=y_i+L_el_y*n_el_y-L_el_y*(r-1)-L_el_y/2;
    b(n,:)=feval(b_fun,xp,yp);
end
b_x=b(:,1);
b_y=b(:,2);

for j=1:n_np_y_b
    [b_x_matrix(:,j)]=b_x((n_np_y_b-j)*n_np_x_b+1:...
                          (n_np_y_b-j)*n_np_x_b+n_np_x_b);
    [b_y_matrix(:,j)]=b_y((n_np_y_b-j)*n_np_x_b+1:...
                          (n_np_y_b-j)*n_np_x_b+n_np_x_b);
end
b_matrix=sqrt(b_x_matrix.^2+b_y_matrix.^2);

% Limit values
u_0_x_min=min(min(time(1).iter(1).u_x_matrix));
u_0_x_max=max(max(time(1).iter(1).u_x_matrix));
u_0_y_min=min(min(time(1).iter(1).u_y_matrix));
u_0_y_max=max(max(time(1).iter(1).u_y_matrix));
p_0_min=min(min(time(1).iter(1).p_matrix));
p_0_max=max(max(time(1).iter(1).p_matrix));
u_f_x_min=min(min(time(length(T)).u_x_matrix));
u_f_x_max=max(max(time(length(T)).u_x_matrix));
u_f_y_min=min(min(time(length(T)).u_y_matrix));
u_f_y_max=max(max(time(length(T)).u_y_matrix));
p_f_min=min(min(time(length(T)).p_matrix));
p_f_max=max(max(time(length(T)).p_matrix));

% Body force vector field
figure('Color',[1 1 1])
axes('FontSize',14)
quiver(X_b',Y_b',b_x_matrix,b_y_matrix,...
       'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5)
title('Body force vector field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Body force field
figure('Color',[1 1 1])
axes('FontSize',14)
contourf(X_b',Y_b',b_matrix,'LineWidth',1)
title('Body force field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Initial condition of velocity in x
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_v',Y_v',time(1).iter(1).u_x_matrix)
title('Initial condition of u_x','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('u_0_x(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
if u_0_x_min~=u_0_x_max
    zlim([u_0_x_min,u_0_x_max])
end

% Initial condition of velocity in y
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_v',Y_v',time(1).iter(1).u_y_matrix)
title('Initial condition of u_y','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('u_0_y(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
if u_0_y_min~=u_0_y_max
    zlim([u_0_y_min,u_0_y_max])
end

% Initial condition of velocity vector field
figure('Color',[1 1 1])
axes('FontSize',14)
quiver(X_v',Y_v',time(1).iter(1).u_x_matrix,time(1).iter(1).u_y_matrix,...
       'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5)
title('Initial condition of velocity vector field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Initial condition of velocity field
figure('Color',[1 1 1])
axes('FontSize',14)
contourf(X_v',Y_v',time(1).iter(1).u_matrix,'LineWidth',1)
title('Initial condition of velocity field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Initial condition of pressure
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_p',Y_p',time(1).iter(1).p_matrix)
title('Initial condition of p','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('p(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
if p_0_min~=p_0_max
    zlim([p_0_min,p_0_max])
end

% Solution of velocity in x
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_v',Y_v',time(length(T)).u_x_matrix)
title('Solution of u_x','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('u_x(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
zlim([u_f_x_min,u_f_x_max])

% Solution of velocity in y
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_v',Y_v',time(length(T)).u_y_matrix)
title('Solution of u_y','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('u_y(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
zlim([u_f_y_min,u_f_y_max])

% Solution of velocity vector field
figure('Color',[1 1 1])
axes('FontSize',14)
quiver(X_v',Y_v',time(length(T)).u_x_matrix,time(length(T)).u_y_matrix,...
       'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5)
title('Solution of velocity vector field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Solution of velocity field
figure('Color',[1 1 1])
axes('FontSize',14)
contourf(X_v',Y_v',time(length(T)).u_matrix,'LineWidth',1)
title('Solution of velocity field','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])

% Solution of pressure
figure('Color',[1 1 1])
axes('FontSize',14)
surf(X_p',Y_p',time(length(T)).p_matrix)
title('Solution of p','FontSize',14)
xlabel('x','FontSize',14)
ylabel('y','FontSize',14)
zlabel('p(x,y)','FontSize',14)
grid on
grid minor
xlim([x_i,x_f])
ylim([y_i,y_f])
zlim([p_f_min,p_f_max])
Warning: Contour not rendered for constant ZData 

Display in command window

disp('MODEL PARAMETERS')
disp('-------------------------------------------------------------------')
fprintf('Length (x)\t\t\t=\t%.1f\n',L_x)
fprintf('Length (y)\t\t\t=\t%.1f\n',L_y)
fprintf('Body force\t\t\t=\t%s\n',char(b_fun))
fprintf('Kinematic viscosity\t\t=\t%.2f\n',v)
fprintf('DOFs constrained (velocity)\t=\t%s\n',dof_constrained_string_v)
fprintf('DOFs constrained (pressure)\t=\t%s\n',dof_constrained_string_p)
fprintf('Boundary conditions (velocity)\t=\t%.1f\n',bound_cond_v)
fprintf('Boundary conditions (pressure)\t=\t%.1f\n',bound_cond_p)
fprintf('Initial condition\t\t=\t%s\n',char(u_0_fun))
fprintf('Initial time\t\t\t=\t%.2f\n',t_i)
fprintf('Final time\t\t\t=\t%.2f\n',t_f)
fprintf('Time step\t\t\t=\t%.2e\n',dt)
fprintf('Theta\t\t\t\t=\t%.2f\n',theta)
disp('-------------------------------------------------------------------')

disp(' ')

disp('FEM PARAMETERS')
disp('-------------------------------------------------------------------')
fprintf('Number of finite elements (x)\t=\t%d\n',n_el_x)
fprintf('Number of finite elements (y)\t=\t%d\n',n_el_y)
fprintf('Length of a finite element (x)\t=\t%.2f\n',L_el_x)
fprintf('Length of a finite element (y)\t=\t%.2f\n',L_el_y)
fprintf('Gauss integration points\t=\t%d\n',n_gauss)
fprintf('Number of nodes (velocity)\t=\t%d\n',n_np_v)
fprintf('Number of nodes (pressure)\t=\t%d\n',n_np_p)
fprintf('DOFs per element (velocity)\t=\t%d\n',dof_el_v)
fprintf('DOFs per element (pressure)\t=\t%d\n',dof_el_p)
fprintf('Total number of DOF (velocity)\t=\t%d\n',dof_v)
fprintf('Total number of DOF (pressure)\t=\t%d\n',dof_p)
disp('-------------------------------------------------------------------')
MODEL PARAMETERS
-------------------------------------------------------------------
Length (x)			=	1.0
Length (y)			=	1.0
Body force			=	@(x,y)[(12-24*y)*x^4+(-24+48*y)*x^3+(-48*y+72*y^2-48*y^3+12)*x^2+(-2+24*y-72*y^2+48*y^3)*x+(1-4*y+12*y^2-8*y^3),(8-48*y+48*y^2)*x^3+(-12+72*y-72*y^2)*x^2+(4-24*y+48*y^2-48*y^3+24*y^4)*x+(-12*y^2+24*y^3-12*y^4)]
Kinematic viscosity		=	1.00
DOFs constrained (velocity)	=	[2:n_np_x_v-1,n_np_v-n_np_x_v+2:n_np_v-1,1:n_np_x_v:n_np_v-n_np_x_v+1,n_np_x_v:n_np_x_v:n_np_v]
DOFs constrained (pressure)	=	[1]
Boundary conditions (velocity)	=	0.0
Boundary conditions (pressure)	=	0.0
Initial condition		=	@(x,y)[0,0]
Initial time			=	0.00
Final time			=	0.10
Time step			=	1.00e-02
Theta				=	1.00
-------------------------------------------------------------------
 
FEM PARAMETERS
-------------------------------------------------------------------
Number of finite elements (x)	=	10
Number of finite elements (y)	=	10
Length of a finite element (x)	=	0.10
Length of a finite element (y)	=	0.10
Gauss integration points	=	9
Number of nodes (velocity)	=	441
Number of nodes (pressure)	=	121
DOFs per element (velocity)	=	9
DOFs per element (pressure)	=	4
Total number of DOF (velocity)	=	441
Total number of DOF (pressure)	=	121
-------------------------------------------------------------------

Animation

% Screen dimensions
scrsz=get(0,'ScreenSize');
bar=64;

% Simulation parameters
fact_ampl_sim=1;
fact_vel_sim=(t_f-t_i)/animation_time;
interval=0.0001;

i=1;
j=1;
tic

figure('Color',[1 1 1],'Position',[0 0 scrsz(3) (scrsz(4)-bar)])
while j<=length(T)

    % Plot
    quiver(X_v',Y_v',time(j).u_x_matrix,time(j).u_y_matrix,...
        'LineWidth',1,'Color',[0 0 0],'AutoScaleFactor',1.5)
    title(['Solution - t = ',num2str(round(T(j)*100)/100),' sec'],...
        'FontSize',14)
    xlabel('x','FontSize',14)
    ylabel('y','FontSize',14)
    zlabel('u(x,y)','FontSize',14)
    grid on
    grid minor
    xlim([x_i,x_f])
    ylim([y_i,y_f])

    % Time calibration
    i=i+1;
    time_real=toc;
    time_simulation=T(j)/fact_vel_sim;
    j=ceil(time_real/dt*fact_vel_sim+eps);

    pause(interval);
end
relative_error=abs(time_simulation-time_real)/abs(time_real);