Asked by Tom
on 8 Aug 2019

I start with a 3 x 46 array where each 3 x1 column in the array is viewed as a position vector of a point, let's say the first vector in the array is p_1, the next one is p_2.

I also have a 3 x 100 array, which I am viewing as the position vectors of 100 points (where these points are always of a lower 'radius' than the points in the 3 x 46 array). I basically need to create a matrix which looks as follows, where each distance is the row vector obtained by vector subtraction of the two vectors:

[Distance from p_1 to first vector in the 3 x 100 array Distance from p_1 to second vector etc ]

[Distance from p_2 to first vector in the 3 x 100 array Distance from p_2 to second vector etc ]

[ Distance from p_3 ..... ]

and so on, until all the positions in the 3 x 46 array are covered, let me know if my intention is not clear.

Answer by Guillaume
on 8 Aug 2019

Edited by Guillaume
on 8 Aug 2019

Accepted Answer

This can be achieved easily without a loop:

%demo data

m = 46

n = 100

A = rand(3,m)

B = rand(3,n)

distance = sqrt(sum((permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])).^2, 3));

The result is a mxn matrix. It basically moves the (X, Y, Z) of each point into the 3rd dimension for both matrix, and move the 2nd dimension of the 1st matrix into the 1st. So points in A are along the rows, and points in B are along the column (and coordinates across the 3rd for both). We now have compatible array sizes, so can create a 3D matrix of (Xa-Xb, Ya-Yb, Za-Zb) for all the points. Square that, sum across (X, Y, Z) so it reduces to a 2D matrix, and take the square root for the distance.

edit: Note that if you want to store the vectors between the points as suggested by Jon instead of the distance (norm of the vectors), then it's simply:

vectors = permute(A, [2, 3, 1]) - permute(B, [3, 2, 1])

vectors(p1, p2, :) is the vector between point A(:, p1) and point B(:, p2).

Tom
on 10 Aug 2019

Another slight subtlety is that some of the functions depend on the normal vector of the point on the radial line, if the position vector of a point on the radial line is called p, this as defined as p divided by the modulus of p.

That means f{1} might not just be a function of rij (the distance from a point on the radial line to the inner points), but also a function of the position vector of the radial point which is being considered, I suppose this can also be added in to the arrayfun?

Tom
on 11 Aug 2019

@Guillame: I've afraid I can't get the above to work, I will try later but I'm running out of time and would be willing to just loops now. Basically I have these functions, where rij is the distance we discussed and normal is just the position vector for one of the points on the radial line divided by the modulus.

function [GV]=gradletSlip(rij,normal)

global Kn;

GV=((1/(8*pi*Kn)*(eye(3)/norm(rij)+(rij'*rij)/((norm(rij))^3)))+sqrt(pi/2)*((-1/(4*pi))*(rij'*normal* ...

(eye(3)/norm(rij)^3 -3*(rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal)))+(3*Kn/(8*pi)*(eye(3)/norm(rij)^3-3* ...

(rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal))/5+(3*Kn/(pi*20)*(-eye(3)/norm(rij)^3+3*(rij'*rij)/norm(rij)^5)));

end

function [FV]=thermletSlip(rij,normal)

global Kn;

FV=((-Kn/(5*pi)*normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5)*(eye(3)-normal'*normal))') ...

*sqrt(pi/2)+((1/(4*pi*norm(rij)^3)*rij*(eye(3)-normal'*normal))')/5;

end

function [GT]=gradletJump(rij,normal)

global Kn;

GT=((3*Kn/(8*pi)*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5)*normal')*sqrt(pi/8)+(-1/(4*pi)*rij'* ...

((normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5))*normal'))/4);

end

function [FT]=thermletJump(rij,normal)

global Kn;

FT=((1/(15*pi*Kn*norm(rij)))+(normal*rij'/(4*pi*norm(rij)^3))*sqrt(pi/8)+(-Kn/(5*pi) ...

*((normal*(eye(3)/norm(rij)^3-3*(rij'*rij)/norm(rij)^5))*normal'))/4);

end

which evaluate to give matrices, columns in the order you gave. I basically just want to evaluate these in the way I explained above with loops if necessary to get a 400 x 138 matrix. I am trying to create the matrix with this code where p is 3 x 46 and sing is 3 x 100 (where totalNodes = 100), but I think something is wrong with the indices and I am getting a 400 x 400 matrix. Bear in mind some functions have 'normal' or 'in' which is just the vector for the position divided by the modulus.

for II=1:46

for IJ=1:totalNodes

in=(p(1:3,II)/(norm(p(1:3,II))))';

r=(p(1:3,II)-sing(1:3,IJ))';

A((1+(II-1)*3):(3+(II-1)*3),(1+(IJ-1)*3):(3+(IJ-1)*3))=gradletSlip(r,in)';

end

for IJ=3*totalNodes+1:4*totalNodes

in=(p(1:3,II)/(norm(p(1:3,II))))';

r=(p(1:3,II)-sing(1:3,IJ-3*totalNodes))';

A((1+(II-1)*3):(3+(II-1)*3),IJ)=thermletSlip(r,in);

end

end

for II=1:46

for IJ=1:totalNodes

in=(p(1:3,II)/(norm(p(1:3,II))))';

r=(p(1:3,II)-sing(1:3,IJ))';

A(3*totalNodes+II,(1+(IJ-1)*3):(3+(IJ-1)*3))=gradletJump(r,in);

end

for IJ=3*totalNodes+1:4*totalNodes

in=(p(1:3,II)/(norm(p(1:3,II))))';

r=(p(1:3,II)-sing(1:3,IJ-3*totalNodes))';

A(3*totalNodes+II,IJ)=thermletJump(r,in);

end

end

Guillaume
on 12 Aug 2019

arrayfun is just a for loop in disguise. It just saves you on having to compute the bounds yourself (more on that later). With either method, if speed is not an issue, I'd create a cell array rather than index into the destination. The code is a lot easier to read this way.

So, I'd do:

matrixbuilder = @(rij, normal) [gradletSlip(rij, normal), thermletSlip(rij, normal);

gradletJump(rij, normal), thermletJump(rij, normal)];

[colp, colsing] = ndgrid(1:size(colp, 2), 1:size(colsing, 2))

distance = arrayfun(@(colp, colsing) matrixbuilder((p(:, colp) - sing(:, colsing))', p(:, colp)' / norm(p(:, colp))), colp, colsing, 'UniformOutput', false)

distmat = cell2mat(distance);

Sign in to comment.

Answer by Jon
on 8 Aug 2019

Edited by Jon
on 8 Aug 2019

So from your description you actually want to find and store the vectors that go from each point in the first matix to each point in the second matrix. Here is a way you could do that

m = 46

n = 100

A = rand(3,m)

B = rand(3,n)

% preallocate 3d matrix to hold results

C = zeros(m,3,n);

% loop through columns of A calculating "distance" vectors and storing

% them

for k = 1:m

C(k,:,:) = B - A(:,k);

end

% to look, for example, at "distance vectors" from the 3rd column of A to all of the

% columns of B you can use

disp(squeeze(C(3,:,:)))

Note, the squeeze function gets rid of extraneous single dimension in multidimensional arrays, which would otherwise make it difficult to see the result

If instead you just want the scalar Euclidean distance from each point in A to each point in B, you could do the following

D = zeros(m,n);

for k = 1:m

D(k,:) = vecnorm(B - A(:,k));

end

If you wanted both you could put both calculations in the same loop.

It may be possible to vectorize this further and eliminate the loop, but it is not immediately obvious to me how to do this. Maybe someone else has a suggestion.

Jon
on 8 Aug 2019

I don't know what your application is so I can't give you any advice on whether you want to define the "distance" vectors as P-S or S-P. I don't think it changes the main logic that I outlined. If I'm understanding correctly it just changes some algebraic signs.

You could probably work out how to permute the 4d array that I have obtained above to get the 3n x 3m array that you want. Otherwise, here is how you can store it this way in the first place. You should check the details, to make sure it is right, but this should give you enough of an idea to go on.

m = 4

n = 5

A = rand(3,m)

B = rand(3,n)

% with 3d result for each pairing

% preallocate 2d array to hold results

J = zeros(n*3,m*3);

for i = 1:m

r = B - A(:,i);

rnorm = vecnorm(r);

for j = 1:n

startRow = (j-1)*3 + 1; % starting row where this 3x3 result is stored

startCol = (i-1)*3 + 1; % starting column where this 3x3 result is stored

J(startRow:startRow+2,startCol:startCol+2) = ...

1/(8*pi*0.7*sqrt(2/pi))*(eye(3)/rnorm(j)+(r(:,j)'*r(:,j))/rnorm(j)^3);

end

end

Tom
on 9 Aug 2019

Guillaume
on 9 Aug 2019

Working in N-D is no more harder than working in 2-D. And in your case it would make it much easier to find the matrix related to a pair of point. You can either have the indices of the points as the first 2 indices, so:

distmat = squeeze(J(p1, p2, :, :)); %3x3 distance matrix between A(:, p1) and B(:, p2)

or as the last 2 indices:

distmat = J(:, :, p1, p2); %3x3 distance matrix between A(:, p1) and B(:, p2)

Storing the same in a 2D array, it's a lot less clear:

distmat = J((p1-1)*3 + (1:3), (p2-1)*3 + (1:3));

Another option, as I suggested is to store the matrices in a 2D cell array, in which case:

distmat = J{p1, p2};

The disadvantage of cell arrays is the memory and speed overhead. For an easy way to do this, see comment to my answer.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.