Functions in matlab of general or very specific use:

(the old and original page was lost for a distraction error so that I need to add them back with time)

  • Function to compute the Jacobian matrices of the Euler equations using the symbolic function of matlab. It was written to avoid to make stupid algebraic mistakes by hand since the algebra of it is simple but tedious!






  • Function to read and plot two different types of computational grids: structured and unstructured. This function cannot plot mixed grids yet although it is very easy to modify the algorithm to add this ability. Feel free to do so if you need to; I'd appreciate to have the modification posted as well in case anyone will do that!

airfoil.jpg
mesh.jpg

function drawmesh(CONNfile, COORDSfile, conn_flag)
 
% plotmesh(CONN, COORDS, conn_flag)
%
% Function to draw an unstructed mesh defined by COORDS y CONN
% Input:  COORDSfile:    node coordinates file
%         CONNfile:      connectivity file
%
%     conn_flag = 1 means that the 1st column of the CONN matrix represents
%*                                      the element number (see below);
%*                        conn_flag = 0 means that the 1st column of the CONN matrix is already
%*                                      a node number (see below);
%*                        ex with conn_flag = 1:
%*                        1 2 5 3
%*                        2 6 3 1
%*                        3 3 1 8
%*                        4 9 8 5
%*                        5 ... ... ...
%*                        ...
%*                        ex with conn_flag = 0:
%*                        2 5 3
%*                        6 3 1
%*                        3 1 8
%*                        9 8 5
%*                        ...
%
%************************************************************************************
 
%Load input files:
CONN = load(CONNfile);
COORDS = load(COORDSfile);
 
[nelem,max_nnodes] = size(CONN);
[nnodes,nsd] = size(COORDS);
 
    if conn_flag == 1
    max_nnodes = max_nnodes - 1
        for i = 1:nnodes
            for j = 1:nsd-1
                TEMP(i,j) = COORDS(i,j+1);
            end
        end
    COORDS = TEMP;
    end
    nsd = nsd - 1;
 
    clear TEMP;
    if conn_flag == 1
        for i = 1:nelem
            for j = 1:max_nnodes
                TEMP(i,j) = CONN(i,j+1);
            end
        end
    CONN = TEMP;
    end
    clear TEMP;
 
nen = size(CONN,2);
 
%Plot the mesh:
    elemx = CONN(:,[1:nen,1])';
    xelems = reshape (COORDS(elemx, 1), nen+1, nelem);
    yelems = reshape (COORDS(elemx, 2), nen+1, nelem);
 
    %plot(xelems(:,1:4),yelems(:,1:4),"k", xelems(:,5:13),yelems(:,5:13),"r");
    plot(xelems, yelems, "k");
 


  • Function to count the maximum number of points surrounding points in a mesh graph (see C verision here )
  • This function comes from the algorithm of Lohner explained in his book [1]. I report its Octave and Matlab implementation here, with a couple of explanations on its use and the type of mesh that it accepts.


    [1] Lohner, Rainald (2008),'Applied computational fluid dynamics techniques', John Wiley and Sons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [max_psup] = point_surr_point(coordinates_file, connectivity_file, npoints, index)
% function [max_psup] = point_surr_point(coordinates_file, connectivity_file, npoints, index)
%
%    Input:
%    connectivity_file: text file of the connectivity matrix.
%                        It must contain the conn matrix in the format:
%                        conn(1:nelem, 1:elnnodes)
%                        where nelem is the mesh number of elements and
%                        elnnodes is the max number of nnodes that any element may have.
%                        ex:
%                        2 5 3
%                        6 3 1
%                        3 1 8
%                        9 8 5
%                        ...
%
%                        The input 'index' must be 1 if the connectivity does not have the first column of elements
%                        (as ex. above), and must be 2 if the connectivity has such column. Ex.:
%
%                        1 2 5 3
%                        2 6 3 1
%                        3 3 1 8
%                        4 9 8 5
%                        5 ...
%
%                        NOTE that the first column does NOT contain the numbering of the elements!!!
%
%    Output:
%    max_psup is the output: maximum number of nodal contacts inside the mesh
%
%
%    PLEASE READ THIS:
%
%    REFERENCES
%    [1] Lohner, Rainald (2008),'Applied computational fluid dynamics techniques', John Wiley and Sons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
disp('#Elements surrounding points#');
%Read connectivity matrix:
    inpoel = load(connectivity_file);
    coords = load(coordinates_file);
 
    [nelem, max_nnode] = size(inpoel);
 
    if(index == 1)
        CONN = inpoel(:,1:max_nnode);
    else
        CONN = inpoel(:,2:max_nnode);
    end
 
    inpoel = CONN;
    [nelem, max_nnode] = size(inpoel);
    [npoints, dumb] = size(coords);
 
    nnodes = npoints;
%Initialize esup2:
    esup2(1:nnodes+1) = 0;
 
 
%Element pass 1: count the number of elements connected to each point:
    for ielem = 1:nelem
        for inode = 1:max_nnode
        %Update storage counter, sotring ahead:
            ipoi1 = inpoel(ielem, inode) + 1;
            esup2(ipoi1) = esup2(ipoi1) + 1;
            %fprintf('PASS 1 esup2[%d] = %d\n', ipoi1, esup2(ipoi1));
        end
    end
 
%Storage/reshugffling pass 1:
    for ipoin = 2:nnodes+1
    %Update storage counter and store:
        esup2(ipoin) = esup2(ipoin) + esup2(ipoin-1);
        %fprintf('PASS 1 esup2[%d] = %d\n', ipoin, esup2(ipoin ));
    end
    %fprintf('PASS 1 esup2[1] = %d\n', esup2(1));
    szesup1 = esup2(nnodes);
    %printf("reallocated size of _esup1 = %d\n",szesup1);
 
%Element pass 2: store the elements in esup1:
    for ielem = 1:nelem
        for inode = 1:max_nnode
        %Update storage counter, storing in esup1:
            ipoin = inpoel(ielem, inode);
%            fprintf('ipoin = %d\n', ipoin);
            istor = esup2(ipoin) + 1;
%            fprintf('istor = %d\n', istor);
            esup2(ipoin) = istor;
%            fprintf('ielem = %d\n', ielem);
            esup1(istor) = ielem;
%            fprintf('PASS 2 esup2[%d] = %d\n', ipoin, esup2(ipoin ));
%            fprintf('PASS 2 esup1[%d] = %d\n', istor, esup1(istor ));
        end
    end
%    size(esup1);
 
%Storage/reshuffling pass 2 and count the number of elements sorrounding point ipoin:
    for ipoin = nnodes+1:-1:2
        nesup(ipoin) = esup2(ipoin) - esup2(ipoin-1);
        esup2(ipoin) = esup2(ipoin-1);
    end
        esup2(1) = 0;
        nesup(1) = esup2(2) - esup2(1);
 
    %for i = 1:nnodes+1
    %    fprintf("nesup2[%d] = %d\n", i, nesup(i))
    %    fprintf('esup2[%d] = %d\n', i, esup2(i))
    %end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Points surrounding points:
disp('\n###############################\n# Points surrounding points #\n###############################\n\n');
 
    lpoin(1:nnodes) = 0;
    psup2(1) = 0;
    istor = 0;
 
    for ipoin = 1:nnodes
    %Loop over the elements sorrounding the point:
        for iesup = esup2(ipoin)+1 : esup2(ipoin+1)
            ielem = esup1(iesup);
            %fprintf("iesup = %d  ielem = %d\n", iesup, ielem)
 
            for inode = 1:max_nnode    %Loop over the nodes of the element
                %disp('CONN(%d,%d) %d\n', ielem, inode, inpoel(ielem,inode));
                jpoin = inpoel(ielem, inode);
                %fprintf("jpoin[%d] = %d\n", inode, jpoin);
 
                if(jpoin ~= ipoin && lpoin(jpoin) ~= ipoin)
                    %Update storage counter, storing ahead, and mark lpoin:
                    istor = istor + 1;
                    %fprintf(" ISTOR = %d\n ", istor);
                    psup1(istor) = jpoin;
                    %fprintf(" psup1[%d] = %d\n ",istor, psup1(istor));
                    lpoin(jpoin) = ipoin;
                    %fprintf(" lpoin[%d] = %d\n ",jpoin, lpoin(jpoin));
                end
            end
        end
 
        %Update storage counters
        psup2(ipoin+1) = istor;
        %fprintf('psup2[%d] = %d\n', ipoin+1,psup2(ipoin+1));
        npsup(ipoin) = psup2(ipoin+1)-psup2(ipoin);
        %fprintf('npsup[%d] = %d\n', ipoin,npsup(ipoin));
        printf(' # Node %d is sorrounded by %d nodes. X[%d] = %f  Z[%d] = %f \n', ipoin,npsup(ipoin), ipoin, coords(ipoin,2), ipoin, coords(ipoin,3));
    end
    %fprintf(" ISTOR = %d\n ", istor)
 
    %Get the maximum number of points surrounding points:
    % This can be used as the maximum number of nonzero elements
    % appearing among all the rows of the stiffness matrix
    max_psup = 0;
    for ipoin = 1:nnodes
        if(npsup(ipoin) > max_psup)
        %printf(" _npsup[%d] = %d\n", ipoin, npsup(ipoin));
            max_psup = npsup(ipoin);
        end
    end
    fprintf("#\n# Maximum number of points surrounding points: %d\n###############################\n\n", max_psup);
 
 
max_nnode
%END code that computes the max number of POINTS sorrounding the points of an unstructured mesh.
 



  • Function to read a sequential text file built as the example below:
# curve0
139.286 0.00677944
155.994 0.0127935
89 89
3 45
5 6
# curve1
957.604 0.020563
1158 0.0212545
1358.3 0.0216704
# curve2
2159.9 0.0146284
2360.3 0.00982097
# curve3
2.9 0.0146284
0.3 0.00982097
1 3

  • Code to read it:
function [data_x, data_y, numbcurves] = read_cuts(file_name, count_curves, write2f, variable_name)
 
% Function to read the file of the cut-off curves from VisIt.
%
% The files have the form:
%
%    # curve0
%    x1 y1
%    x2 y2
%    ... ...
%    # curvei
%    x... y...
%    ... ...
%    xn yn
%
%    The argument count_curves must be 1 or 0 is you want only to count the curves in the file or if you want to read all the data as well.
%    The argument write2f must be 1 or 0 is you want to write to file the output or not respectively.
%    The argument variable_name is needed to write the output file when wrtie2f == 1.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Auxiliary arrays used throughout:
%    Ic(i=1:ncurves)                        It contains the initial absolute line index of every new curve
%    Jc(i=1:ncurves) = Ic(i) - i            It contains the initial line index as if the #curve headers did not appear in the file
%    Lc(i=1:ncurves) = Jc(i) + nel(i) - 1    It contains the final   line index as if the #curve headers did not appear in the file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
    file_name_id = fopen(file_name, "r");
    if( file_name_id == -1)
        printf("!!! Error, could not open file '%s'\n", file_name);
        matrix=[];
        return;
    end
 
        incurves = 0;
        line_cntr = 0;
        nel = 0;
        while( !feof(file_name_id))
 
            string = fgets(file_name_id);
            line_cntr = line_cntr + 1;
            line_cntr_array(line_cntr) = line_cntr;
            %fprintf(" %s ", string);
 
 
            % 1) Count the number of curves in the file by identifying the headers # curveN:
                if( strncmp(string,"# curve",7) == 1 )
                    %fprintf(" %s ", string);
                    incurves = incurves + 1;
                    curve_line(incurves) = line_cntr;
                    Ic(incurves) = line_cntr + 1;
                    Jc(incurves) = Ic(incurves) - incurves;
 
                end
        end %WHILE
 
        % Total number of lines in the file (inclusive of the strings # curveN
            line_cntr = line_cntr - 1;
 
        % Total number of curves: ncurves.
            curve_cntr = incurves;
            numbcurves = curve_cntr;
 
        %if(count_curves == 1)
        %    data_x = -1;
        %    data_y = -1;
        %return;
        %end
 
        nel(1) = curve_line(2) - 2;
        Lc(1) = Jc(1) + nel(1) - 1;
        for i = 2:curve_cntr-1
 
            nel(i) = curve_line(i+1) - curve_line(i) - 1;
            Lc(i) = Jc(i) + nel(i) - 1;
        end
        nel(curve_cntr) = line_cntr - curve_line(curve_cntr) + 1;
        Lc(curve_cntr) = Jc(curve_cntr) + nel(curve_cntr) - 1;
 
        string_last = strcat("# curve", int2str(incurves));
 
 
    % END COUNTING THE LINES AND THE ELEMENTS OF EACH CURVE STORED IN THE COORDINATE FILE
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    % Now we can store the elements in proper arrays:
    frewind(file_name_id);
 
    i = 0;
    x_curvei = 0.0;
    y_curvei = 0.0;
 
    data = [];
 
    while( !feof(file_name_id))
 
        string = fgets(file_name_id);  % reads in a line
        [tmp_data, n_items] = sscanf(string, "%f %f\n");  % scans the line for 2 numbers
 
        if (n_items == 2)      % if 2 numbers were found on the line
            data(end+1, :) = tmp_data;  % add that to the data matrix
        end
    end
 
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Store the arrays for each curves:
 
    l = 0;
    for i = 1:curve_cntr
    out_file = variable_name;
 
        %Set the name of the output file for every curve:
            curveid = int2str(i);
            out_file = strcat(out_file, "_curveCut_");
            out_file = strcat(out_file, curveid);
            out_file = strcat(out_file, ".dat")
 
        %Open file for writing:
            file_id = fopen(out_file, "w");
 
        %Start storing the values of the original file into proper arrays:
        k = 0;
 
        for j = Jc(i):Lc(i)
            k = k + 1;
            l = l + 1;
 
            data_x(k,i) = data(l,1);
            data_y(k,i) = data(l,2);
 
            %Write them to file: file of the form: xi yi
                fprintf(file_id, "%f %f\n", data_x(k,i), data_y(k,i));
                %printf(" %f %f\n", data_x(k,i), data_y(k,i));
        end
        %printf("\n");
 
    fclose(file_id);
    end
 
 
    %Data is of size: [nel(icurve)]x[ncurves]
    %data_x(:,1)
 
    %Only takes elements from 2nd to the one before last for plotting reasons:
    data_x = data_x(2:nel(incurves),:);
    data_y = data_y(2:nel(incurves),:);
 
    [rows,colmn] = size(data_x);
 
    %plot(data_x(:,1:incurves),data_y(:,1:incurves));
 
    %Close file
    fclose(file_name_id);
 
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Miscellaneous:
    %Write data to a new file: only if you selected write2f == 1
    if(write2f == 1)
 
        variable_name = strcat(variable_name, "_cutCurves.dat");
 
        file_name_id = fopen(variable_name, "w");
        if( file_name_id == -1)
            printf("!!! Error, could not open file '%s'\n", file_name);
            matrix=[];
            return;
        else
            printf("! WARNING: printing to file not working YET for this function\n The function will return now\n\n!");
            return;
            %for i = 1:incurves
            %    fprintf(file_name_id, "%f %f",data_x(:,), data_y(:,1:incurves));
            %end
        end
 
        fclose(file_name_id);
    end
 
 
end
 



  • Linear regression: function that takes in the coordinates (X,Y) of a set of scattered points and computes the coefficients of the regression curve Yreg = aX + b

function [YY] = linearRegression(X, Y, plot_flag)
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input:  X, Y -> coordinates of the scattered values
% Output: YY   -> linear function after computiong the regression of Y
%
%    Simone Marras
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    [m,n] = size(X);
    X_mean = sum(X)/m;
    Y_mean = sum(Y)/m;
 
    Xdiff = X-X_mean;
    Ydiff = Y-Y_mean;
 
    b = sum(Xdiff.*Ydiff)/sum(Xdiff.*Xdiff);
    a = Y_mean - b*X_mean;
 
    YY = a + b*X;
 
    if(plot_flag == 1)
        plot(YY,X, Y, X);
    end
 
return;
 





  • curve_length.m
  • Function that returns the length of a curve given on the cartesian space (x,y):
  • % Function that computes the lenght of a curve
    % between two points xa and xb given its definition by arrays x and y:
    function [curveLength] = curve_length(x,y, xa, xb, n)
        % Interpolate to the inermediate points: if the original curve y(x) has more points thatn n, 
        % then use as many points as the original y-curve already has:
        if(length(y) > n)
            n = length(y);
        end
        %Intermediate epoints where to compute the lenghts of each single segment in which the curve is divided
            xi= xa:(xb-xa)/n:xb;
     
        % Interpolate to find the yi corresponding to the intermediate xi points between xa and xb:
            yi=interp1(x,y,xi,'spline');
        curveLength = 0.0;
        dxi = diff(xi);
        dyixi = diff(yi)./diff(xi);
     
        for i = 1:n
            integrand = dxi(i)*sqrt(1.0 + dyixi(i)^2);
            curveLength = curveLength + integrand;
        end





  • Plot horizontal and vertical lines in Matlab:

%Horizontal from x=1.0 to x = 3.0 at z = 4.0
plot([1, 3],[4, 4])
 
 
%Vertical from z=0.0 to z = 2.0 at x = 1.0
plot([1, 1],[0, 2])
 
 





  • Plotting in Matlab: a good linkby Didier Gonze