Can I vectorize one of the loops ?

Hello everbody,

I have a performance problem with the two loops below. As I am quite new to matlab, I have no idea how to vectorize them and if its even possible.

What the inner loop does is, he brings the input data on specific identical x coordinates for the interpoaltion. The input data is a pressure distribution, but every pressure distribution has other x coordinates. For interpolation between the different pressure distributions, they need identical x_c coordinates.

So by finding the left and right neighbor of the x_c coordinate, I linearly interpolate between them, to find the pressure coefficient of the point, I need to know, so that for the interpolation at the and I have pressure distributions with identical x-coordinates. The outer loop is just doing this for all 280 x_c coordinates.

Thank you in advance. If you have any questions, let me know.

%% for all 280 points on the lower and upper airfoil
    for x_c=1:1:280
        x_c_inter=x_c0(x_c);                        % support points
        for a=1:1:3                                    % for all 3 points of the triangle of the delaunay triangulation 

            new_upper_airfoil_unique=new_upper_airfoil_loop(new_upper_airfoil_loop(:,1)==DataId_unique(a,:),2:3); % find the data of vertices of the triangle
            new_lower_airfoil_unique=new_lower_airfoil_loop(new_lower_airfoil_loop(:,1)==DataId_unique(a,:),2:3); 

            %% 145 points on the upper airfoil
            if x_c < 146
                FIND=new_upper_airfoil_unique(abs(new_upper_airfoil_unique(:,1)-x_c_inter)<0.0005,:);

                if(isempty(FIND))

                    leftNeighborIndex=nearestpoint(x_c_inter,new_upper_airfoil_unique(:,1),'previous'); % find left neighbor of the point
                    if isnan(leftNeighborIndex)                        
                        leftNeighborIndex=1;                                                             % if leftNeighbor is NaN, its because the prevoius neighbor is at index 0, so I use index=1 and extrapolate
                    end

                    rightNeighborIndex=leftNeighborIndex+1;

                    if rightNeighborIndex > size(new_upper_airfoil_unique,1)
                        rightNeighborIndex=leftNeighborIndex-1;                                           % if rightneighbor is bigger than the size of the vector, extrapolation
                    end
                    % interpolation of the pressure coefficient with the left and right neighbor, so it has the right x_c coordinate
                    average=interp1([new_upper_airfoil_unique(leftNeighborIndex,1),new_upper_airfoil_unique(rightNeighborIndex,1)], [new_upper_airfoil_unique(leftNeighborIndex,2),new_upper_airfoil_unique(rightNeighborIndex,2)], x_c_inter, 'linear','extrap');


                    % write it into a temporary array for the last interpolation
                    Temp(a,:)=[x_c_inter,average];

                else
                    % if a pint is found within the delta, just use this point, but change the x_c coordinate
                    FIND(1,1)=x_c_inter;
                    Temp(a,:)=FIND(1,:);

                end 
            else
                %% lower airfoil
                FIND=new_lower_airfoil_unique(abs(new_lower_airfoil_unique(:,1)-x_c_inter)<0.00008,:);
                if(isempty(FIND))

                    leftNeighborIndex=nearestpoint(x_c_inter,new_lower_airfoil_unique(:,1),'previous');
                    if isnan(leftNeighborIndex)
                        leftNeighborIndex=size(new_lower_airfoil_unique,1);
                    end
                    rightNeighborIndex=leftNeighborIndex-1;

                    average=interp1([new_lower_airfoil_unique(leftNeighborIndex,1),new_lower_airfoil_unique(rightNeighborIndex,1)], [new_lower_airfoil_unique(leftNeighborIndex,2),new_lower_airfoil_unique(rightNeighborIndex,2)], x_c_inter, 'linear','extrap');

                    Temp(a,:)=[x_c_inter,average];

                else
                    FIND(1,1)=x_c_inter;
                    Temp(a,:)=FIND(1,:);
                end  
            end
        
        
        end
        cp_inter=dot(bc',Temp(:,2)');  %% interpolation with barycentric coordinates


        
        DV(x_c,:)=[x_c_inter,cp_inter];

    end
end

[darova]

Can you please make a simple drawing or scheme of your idea?

[Gerrit]

I can try. I am doing an interpolation of pressure distributions of an airfoil. A pressure distribution is basically a 2D Graph with the coordinate X/C and the pressure coefficient cp.

The problem is, that for the interpolation, I need every pressure distribution of my input to have same x/c coordinates. Thats why i do a linear interpolation.

[darova]

Maybe you can attach an airfoil? Can you show the result you want to see?

[Mohammad Sami]

interp1 can handle vectors of x , v and inter_x, as long as your x is sorted. unless for some reason you are getting incorrect answer, you don't really need to use a for loop. interp1 will automatically take the two nearest point to interpolate.

[Image Analyst]

You have 3*280 = 840 iterations. I can do tens of millions of iterations in well under a second, so I don't think the for loop is what's slowing you down. Try to use some of the timing tools in MATLAB to see where the time is really spent.

[Gerrit]

Yeah, but I have 4 Interpolations, so 4*840 and that for 7500 cases, so 4*840*7500.... Thats why every millisecoud counts.

I did not know that interp1 is automatically using the nearestpoints. I will try to increase performance with this.