comparing columns between differently sized matrices

Hi, I am writing a finite element analysis code in MATLAB using tetrahedral elements. In the formulation of the problem, I have a 3xM surface connectivity matrix containing global node numbers and I have a 4xN volume connectivity matrix where M is the number of surfaces and N is the number of elements (M≠N). Each column in the surface connectivity matrix represents the surface number and each column in the volume connectivity matrix represents the element number. I need to compare each column of the 3xM surface connectivity matrix to each column of the 4xN volume connectivity matrix to determine which element number corresponds to each surface. Example matrices are given below.

       surface connectivity matrix (3xM)       
           1     4     4     1
           3     1     7     4
           2     2     2     7
       Volume connectivity matrix (4xN)
           1     1     1     1     3     4
           4     8     7     7     2     7  
           7     3     8     9     8     2
           2     2     2     8    10    10

As can be seen in this example surface #1 corresponds to element #2, surface #2 corresponds to element #1, surface #3 corresponds to elements #1 and #6 (this surface is shared between two elements) etc. To complicate matters, the global node ordering for the surfaces is different than the ordering for the element numbers, but the ordering doesn’t matter. I can perform the comparison using a for loop (looping M times) and searching the entire 4xN matrix during each loop iteration, but since M and N are >10^5, the computation time is unreasonable. Is there a way to perform this comparison operation in a vectorized way to avoid loops?

Answer 1

[Fangjun Jiang]

SCM=[1     4     4     1
     3     1     2     4
     2     2     7     7];
VCM=[1     1     1     1     3     4
     4     8     7     7     2     7  
     7     3     8     9     8     2
     2     2     2     8    10    10];
M=size(SCM,2);
N=size(VCM,2);

CORRESP=false(M,N);
for k=1:M
    for j=1:N
        CORRESP(k,j)=all(ismember(SCM(:,k),VCM(:,j)));
    end
end

CORRESP =

       0     1     0     0     0     0
       1     0     0     0     0     0
       1     0     0     0     0     1
       1     0     0     0     0     0

It shows surface #1 corresponds to element #2, surface #2 corresponds to element #1, surface #3 corresponds to elements #1 and #6, surface #4 corresponds to elements #1.

[Joel]

Yes thank you I edited the question to fix the typo.

When I use a loop I do in fact use ismember. My solution only requires one loop over the M surfaces. However, since M is large the loop is very slow. Is there a way to do this without a loop?

Thank you.

[Fangjun Jiang]

Can you show your code so we understand how the loop is done?

Answer 2

[Joel]

My main question above is very general. A solution to it would be ideal. However my current solution is slightly less general. In general each face can be connected to at most 2 elements (as I've dicussed above). In this solution I search the mesh in such a way as to eliminate one element (e.g. by applying a boundary condition). I don't think the details are too relvant as it will only complicate things. So in the below code, each face is connected to only one element (SCM = surface connectivity matrix and VCM = volume connectivity matrix).

     for i=1:size(SCM,2) % loop over M (number of columns in SCM)
         b=double(ismember(VCM,SCM(:,i))); % search the VCM
         elem(i)=find(sum(b)==3); % find element number 
     end

Again my end goal is simply to compare two columns of different sizes within matrices of different sizes without using a loop.

Thank you.

[Fangjun Jiang]

Not using a for-loop doesn't necessarily mean fast. There are many factor can impact the efficiency. In your code, double() and sum() seem to be unnecessary. You can use all() for that. I hope you pre-allocate elem.

Again, show your actual code so others can point out any place that can be improved. There are many examples on this forum that the best approach is application dependent.

[Joel]

I've pre-allocated elem. And you are correct that double() wasn't necessary, but all() will not work since matrix b is a 4xN matrix and each column contains at least one zero. Thus I use sum() to find the column that equals 3 (i.e. all three global node numbers exist).

Although I agree the above code may not be optimized perfectly, the slow down is not in the double() or sum(); it is in the use of a loop.

I still think that showing my entire code (which is long) is unnecessary and will only complicate matters. The above code snipit describes exactly what I am trying to do (but using a loop) and thus is application specific. The only other detail for my current specific application is that each face will correspond to only one element (not two which is general).

I do appreciate your time and your help thus far.

[Fangjun Jiang]

So your algorithm is to compare the whole VCM matrix(4xN) with each column of SCM matrix(3xM) at a time, the result b is also a 4xN matrix. You then check if any of the b column sums to 3.

When N is large (10^5), I think this might not be efficient. I think you agree that my double-loop algorithm works too. Why don't you time both to see the difference. I'll probably do it too when I have spare time.

[Jan]

SUM operates on LOGICAL arrays, so there is no need to cast to DOUBLE in "b=double(ismember( ..."

[Joel]

Using a quite coarse mesh (N=16844 and M=2090) Fangjun Jiang's code took 903.06 s and my version took 2.05 s.

Using my above loop technique, much of the time is actually taken up by the ismember function. Again all these point to the necessity of avoiding the ismember function within a loop.

I would appreciate any ideas in this regard. Thank you.

[Fangjun Jiang]

Now I thought it more, I think your approach is probably most efficient. I can't think of any other ways. Don't be obsessed with "no-loop". The loop in new versions of MATLAB is very efficient. There are functions like cellfun(), arrayfun(), bxsfun() etc that can be used to eliminate the loop but I don't think there will be a good use for your case.

If you are still interested in improving the efficiency of your code, many other experts can probably point out other areas to improve if you post your code.