How to get a loop faster?

3 次查看（过去 30 天）
Carles Rafels Ybern
编辑： FirefoxMetzger 2016-8-23
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
  7 个评论
显示 5更早的评论
Carles Rafels Ybern
Carles Rafels Ybern 2016-8-15
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
KSSV 2016-8-16
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

FirefoxMetzger
FirefoxMetzger 2016-8-16
编辑：FirefoxMetzger 2016-8-23
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.

类别

在 Help Center 和 File Exchange 中查找有关 Data Type Identification 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by