How to get a loop faster?

Hi, I'll try to explain briefly my problem...

I need to create the coordinates and connectivities files for getting a calculation by finite element method. In other words, I have to create the FEM mesh. For that, I need to create the connectivities that relates to the coordinates file but these last can not have duplications.

I start from a 3D regular hexahedron model with Nk, Nj, Ni-dimensions. Total of cells = Nj*Nk*Ni. So, I have got three matrices (one for each direction X-Y-Z) of three-dimensional (2*Ni,2*Nj,2*Nk). I.e, gives 1 cell (I,J,K) you can get the eight values of X/Y/Z that form 1 cell by using these three matrices.

To create the connectivities and coordinates of my model, I have to check if each coordinate (X,Y,Z) has already been created by a previous cell and then, gives that same id-node. And then, I don't create repeated coordinates.

I've been looking for a heap information about how can I get to make faster the loops (Preallocating variables, Parallel Computing but I think in that case I cannot use,...).

My last case executed was made for 250000 cells and the Elapsed time was 900seg. (aprox. 15 min). The problem is when I need to run a greater example, for example, more than 4000000 cells, time is untenable. I need help, please!!

I show you my script:

A) Some definitions:

• Ni=model.dim(1); Nj=model.dim(2); Nk=model.dim(3);
• colCoo: vector that gives the order of how to pick the coordinates of each cell.
• ivox: each cell.
• model.fX & model.fY & model.fZ: are the coordinate matrices of X/Y/Z values.
• Each model.variable is preallocating for more speed.

B) Script:

for K=1:model.dim(3)

      for J=1:model.dim(2)
          for I=1:model.dim(1)
              model.indxconec(ivox,1)=ivox;
              colCoo=[iorder(1:3:end-2,ivox),iorder(2:3:end-1,ivox),iorder(3:3:end,ivox)];
              val=nan(8,3); %preallocating for speed
              for t=1:8
                  val(t,1)=model.fX(colCoo(t,1),colCoo(t,2),colCoo(t,3));
                  val(t,2)=model.fY(colCoo(t,1),colCoo(t,2),colCoo(t,3));
                  val(t,3)=model.fZ(colCoo(t,1),colCoo(t,2),colCoo(t,3));
              end
              indx=zeros(1,8);
              loc= find(any(all(bsxfun(@eq,reshape(val.',1,3,[]),model.coo),2),3));
              [~, index]=ismember(model.coo(loc,1:3),val,'rows');
              indx(1,index(1:1:end))=loc;
              [Posindx]=find(indx==0);
              model.coo(idcoo+1:idcoo+size(Posindx,2),1:3)=val(Posindx,:); clear val
              model.indxcoo(idcoo+1:idcoo+size(Posindx,2),1)=(idcoo+1:idcoo+size(Posindx,2));
              idcoo=idcoo+size(Posindx,2);% lenght of coo put before
              for e=1:8
                  if indx(1,e)==0
                      indconc=indconc+1;
                      indx(e)=indconc;
                  end
              end
              model.conec(ivox,1:8)=indx(1,1:8);
          end %end I-loop
      end %end J-loop
  end %end K-loop

[per isakson]

Did you run Profile to Improve Performance?

[Carles Rafels Ybern]

hi isakson,

yes and thanks for your interest. I ran "Profile to improve Performance" and the most of execution time is spent on lines:

loc= find(any(all(bsxfun(@eq,reshape(val.',1,3,[]),model.coo),2),3));

[~, index]=ismember(model.coo(loc,1:3),val,'rows');

but I don't know how can I reduce time into these internal function. Do you know if exist anyway to do the same faster? Especially, "ismember" function is slower when model.coo is greater. Is there any function that makes the same but faster?

Thanks, Carles R.

[Chad Greene]

In my experience, ismember tends to slow things down mostly due to sorting. It sorts everything you give it, every time you call it. Any way you can make sure everything is sorted before the loop, then do some sort of logical check?

[per isakson]

Didn't the Code Analyzer suggest that you skip find and let loc be a logical index? It should help a tiny bit.

The expression, bsxfun(@eq,reshape(val.',1,3,[]),model.coo), makes me wonder whether the values of val are whole numbers?

[per isakson]

Maybe this thread contains something of interest to you speedup of ismember

[Carles Rafels Ybern]

Thanks isakson for your help. Answering your question: no, are rational numbers since, val are the 8 coordinates (x,y,z-values) of each hexahedron cell. But, loc are whole numbers that give the position of repeated coordinate.

Moreover, I saw that link and it was where I found out that expression but thanks again.

Lastly, I reformulated all the script basing on other things and I've gotten to create the conectivity file faster.

[KSSV]

Dear Carles

Generating FEM mesh needs a meshing methodology to be followed. In the respective methodology, they discuss the nodal connectivities as well. In there, you need not to bother whether a node is repeating or not. I can say, generating 3D mesh is a bit tedious job. May I know what methodology you are following? What references you have?

Answer 1

[FirefoxMetzger]

Here is my result, with a lot more commenting. I couldn't test if I made a mistake, because ... no MWE. Also note that for some reason ivox is never generated from (K,J,I) - I assumed it would be, to make sense of the code

    % iorder stores hexahedron related stuff size(iorder) == [3*8 K*I*J]
    % iorder(:,idx) <-- a single hexahedron
    % Structure of a single hexahedron:
    % Point1 - Point2 - ...
    % XYZ XYZ XYZ packed

    % ivox : The current hexahedron

    for K=1:model.dim(3)
        for J=1:model.dim(2)
            for I=1:model.dim(1)
                % set an element of an usused variable to a funny value
                model.indxconec(ivox,1) = ivox; % what does this do !?

                % reshape the hexahedron coordinates
                hexIdxCoordinates=[iorder(1:3:end-2,ivox) iorder(2:3:end-1,ivox) iorder(3:3:end,ivox)];

                % look up real world coordinates of hexahedrons 
                WorldCoordinateIdx = sub2ind(size(model.fX),hexIdxCoordinates(:,1),hexIdxCoordinates(:,2),hexIdxCoordinates(:,3));
                hexWorldCoordinates = [model.fX(WorldCoordinateIdx) model.fY(WorldCoordinateIdx) model.fZ(WorldCoordinateIdx)];

                % % check if any of the hexahedron's points was previously 
                % % visited and find the corresponding indices
                indx=zeros(1,8);
                % compare each point to all points in model.coordinates
                tmp = reshape(hexWorldCoordinates.',1,3,[]);
                compAllPoints = bsxfun(@eq,tmp,model.coordinates);
                % find where the viseted points are stored in model.coordinates 
                idxPointInModel = find(any(all(compAllPoints,2),3));
                % find the corespoinding positions in the hexahedron's world
                % coordinates
                [~, idxPointInHex] = ismember(...
                    model.coordinates(idxPointInModel,:),...
                    hexWorldCoordinates,'rows');
                indx(1,idxPointInHex) = idxPointInModel;
                % find all points that have NOT been entered in
                % model.coordinates yet
                idxUnusedPoints=find(indx==0);

                % write the elements currently not in the Coordinates list into
                % the coordinates list
                newPointsIdx = corrdinatesEnd+(1:numel(idxUnusedPoints));
                model.coordinates(newPointsIdx,:) = hexWorldCoordinates(idxUnusedPoints,:);

                % add some funny numbers to an unused array
                model.indxcoo(newPointsIdx,1) = newPointsIdx;

                % save where a hexahedron is connected ?
                indx(indx == 0) = newPointsIdx;
                model.conec(ivox,1:8) = indx; 

                % increase the pointer to the end of the model.coordinate array
                corrdinatesEnd = corrdinatesEnd+numel(idxUnusedPoints);% lenght of coo put before
            end %end I-loop
        end %end J-loop
    end %end K-loop

Currently runtime is O([num Hexahedrons]^2) due to the bsxfun. It also has a pretty big constant factor due to ismember(...) and because we operate on the 8 points not on hexahedrons.

The first thing that can be done is reducing the constant factor by refactoring. In particular we can remove the bsxfun, all the find calls and replace the one ismember call with

    [~, idxInPointCloud] = ismember(hexCoordinates,model.pointcloud(1:coordinatesEnd,:),'rows');

This returns a vector of the storage location of coordinates that were already processed. If they haven't been visited, the vector contains 0, see ismember.

Secondly remove calls to model.variable. It doesn't matter if model is a struct or class (class is worse though). It makes your code slower, not faster! Also be carefull with growing arrays in combination with structs.

    a = struct;
    for idx = 1:1e6
        a.variable(idx) = 42;
    end

will not show the array growth warning, but that is NOT because there is no growing array. It's because the code analyzer can't detect it. Internally a.variable still changes size during each iteration.

This results in:

    % iorder stores hexahedron coordinates size(iorder) == [3*8 K*I*J]
    % iorder(:,idx) <-- a single hexahedron
    % Structure of a single hexahedron:
    % Point1 - Point2 - ...
    % XYZ XYZ XYZ packed
    % Note: The coordinates are actually indices. They still need to be mapped
    % into the real world coordinates

    % ivox : The current poisition in the loop

    % preallocate variables _before_ any loop - for speed :)
    numHexahedrons = INSERT_CORRECT_NUMBER_HERE;

    dim = model.dim;
    fX = model.fX;
    fY = model.fY;
    fZ = model.fZ;
    pointcloud = nan(8*numHexahedrons,3);
    coordinatesEnd = 0;
    connectionArray = nan(numHexahedrons,8);

    for K=1:dim(3)
        for J=1:dim(2)
            for I=1:dim(1)
                % get ivox from K, J, I
                % TODO - FILL ME WITH THE RIGHT CODE

                % set an element of an usused array to some funny value
                model.indxconec(ivox,1) = ivox; % what does this do !?

                % reshape the hexahedron coordinates
                hexIdxCoordinates=[iorder(1:3:end-2,ivox) iorder(2:3:end-1,ivox) iorder(3:3:end,ivox)];

                % look up coordinates of the hexahedron
                CoordinateIdx = sub2ind(size(fX),hexIdxCoordinates(:,1),hexIdxCoordinates(:,2),hexIdxCoordinates(:,3));
                hexCoordinates = [fX(CoordinateIdx) fY(CoordinateIdx) fZ(CoordinateIdx)];

                % find points that occur in hex and list of all coordinates
                [~, idxInPointCloud] = ismember(hexCoordinates,pointcloud(1:coordinatesEnd,:),'rows');
                % idxInPointCloud == 0 --> Point not in Pointcloud

                % add points not in pointcloud
                numNewPoints = sum(idxInPointCloud == 0);
                newPointsIdx = corrdinatesEnd + (1:numNewPoints);
                pointcloud(newPointsIdx,:) = hexCoordinates(idxInPointCloud == 0,:);

                % update the index with the just added points
                % and store the connectivity
                idxInPointCloud(idxInPointCloud == 0) = newPointsIdx;
                connectionArray(ivox,:) = idxInPointCloud;

                % update end pointer of corrdinate list
                corrdinatesEnd = corrdinatesEnd + numNewPoints;
            end %end I-loop
        end %end J-loop
    end %end K-loop

    % add some funny numbers to an unused array
    model.indxcoo(:,1) = 1:coordinatesEnd;

    % clean up unused pointcloud values
    pointcloud(coordinatesEnd+(1:end),:) = [];

    % save data back to the (slow) model struct
    model.connec = connectionArray;
    model.pointcloud = pointcloud;

You can get even faster by writing your own function searching the pointcloud, which is currently done by ismember(). However I think this should already be faster then the previous version and this answer is way to long already :-).

Hope I could help.