I am writing code to do finite volume method operations on an NI*NJ size grid, where all the values are stored in a 1D array, defined as in comment section 1. The double loop moves through each grid element and performs calculations for the grid element, which are various functions of the surrounding elements (to the east, north, west and south).
I have an example code block here below, and my challenge is to clean up the code to eliminate some of the repitition I see. I wonder if somehow employing classes that contain properties e, w, n, s
would be a beneficial approach to the problem, although I'm not well-versed in class constructions.
As can be seen, I have defined an anonymous function and then called it to calculate some value at each side (e,w,n,s). This pattern recurs numerous times throughout the entire code and it would be nice to have some way of combining such code blocks of quadruplets into one line of simple code.
arrayRe = initialize; % goes from 1:(NI*NJ)
arrayRn = initialize;
for i = 2:NI-1
for j = 2:NJ-1
% 1. coordinates (P is the center; all others are located around it (north, east, west, south))
iP = NI*(j-1) + i;
iE = iP + 1;
iW = iP - 1;
iN = iP + NI;
iS = iP - NI;
% 2. interpolation ratios
R.e = arrayRe(iP);
R.w = 1 - arrayRe(iW);
R.n = arrayRn(iP);
R.s = 1 - arrayRn(iS);
% interpolate a given value
interpolatevalue = @(array, R, i) ...
R*array(i) + (1-R)*array(iP);
% ax calc
ax_e = interpolatevalue(axArray, R.e, iE);
ax_w = interpolatevalue(axArray, R.w, iW);
ay_n = interpolatevalue(ayArray, R.n, iN);
ay_s = interpolatevalue(ayArray, R.s, iS);
% bx calc
bx_e = interpolatevalue(bxArray, R.e, iE);
bx_w = interpolatevalue(bxArray, R.w, iW);
by_n = interpolatevalue(byArray, R.n, iN);
by_s = interpolatevalue(byArray, R.s, iS);
end
end
Something I tried was to cut out the 4 lines of a particular block and put them into a function m file that returned a struct. For instance:
function a = axycalc(axArray, ayArray, R, iStruct, interpfun)
iE = iStruct.e ;
iW = iStruct.w ;
iS = iStruct.s ;
iN = iStruct.n ;
a.x_e = interpfun(axArray, R_e, iE);
a.x_w = interpfun(axArray, R_w, iW);
a.y_n = interpfun(ayArray, R_n, iN);
a.y_s = interpfun(ayArray, R_s, iS);
end
And then calling the function as follows:
a = axycalc(axArray, ayArray, R, struct('e',iE,'w',iW,'n',iN,'s',iS), @interpolatevalue)
I don't know whether this adds much in terms of readability, and creates a large amount of other m 10-line m files that I then need to keep track of. Is there a better way?