gtStrainPlaneNormals: optimized for deformation of single plane normal

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Strain2/gtStrainPlaneNormals.m

@@ -21,73 +21,81 @@ function [pldef, drel] = gtStrainPlaneNormals(pl, defT)
 %   drel  - relative elongations along the (deformed) plane normals
 %
 
+num_pls = size(pl, 2);
+num_def = size(defT, 3);
 
 % Get an arbitrary vector e1 perpendicular to n:
 [~, maxcomp] = max(abs(pl), [], 1);
 
 
-ee = zeros(3,size(pl,2)*3);
+ee = zeros(3, num_pls * 3);
 
-zero = zeros(1,size(pl,2));
+zero = zeros(1, num_pls);
 
-ee(:,1:3:end) = [ pl(2,:); -pl(1,:);     zero];
-ee(:,2:3:end) = [    zero;  pl(3,:); -pl(2,:)];
-ee(:,3:3:end) = [-pl(3,:);     zero;  pl(1,:)];
+ee(:, 1:3:end) = [ pl(2, :); -pl(1, :);     zero];
+ee(:, 2:3:end) = [    zero;  pl(3, :); -pl(2, :)];
+ee(:, 3:3:end) = [-pl(3, :);     zero;  pl(1, :)];
 
-ind2 = (0:size(pl,2)-1)*3 + maxcomp;
-e1   = ee(:,ind2);
+ind2 = (0:num_pls-1) * 3 + maxcomp;
+e1   = ee(:, ind2);
 
 
-% Cross product of e1 and pl: e2=cross(e1,pl)
-e2 = [pl(2,:).*e1(3,:) - pl(3,:).*e1(2,:);
-      pl(3,:).*e1(1,:) - pl(1,:).*e1(3,:);
-      pl(1,:).*e1(2,:) - pl(2,:).*e1(1,:)];
+% Cross product of e1 and pl: e2=cross(e1, pl)
+e2 = [pl(2, :) .* e1(3, :) - pl(3, :) .* e1(2, :);
+      pl(3, :) .* e1(1, :) - pl(1, :) .* e1(3, :);
+      pl(1, :) .* e1(2, :) - pl(2, :) .* e1(1, :)];
 
   
 % Check   
 % ne2 = sqrt(sum(e2.*e2, 1));
-% e2n = e2./ne2([1 1 1],:);
+% e2n = e2./ne2([1 1 1], :);
 % 
 % ne1 = sqrt(sum(e1.*e1, 1));
-% e1n = e1./ne1([1 1 1],:);
+% e1n = e1./ne1([1 1 1], :);
 %  
-% cross(e1n,e2n) - pl
+% cross(e1n, e2n) - pl
 
 
 % check flip e2 if needed to get them right handed:
-% ch = [e1(2,:).*e2(3,:) - e1(3,:).*e2(2,:);
-%       e1(3,:).*e2(1,:) - e1(1,:).*e2(3,:);
-%       e1(1,:).*e2(2,:) - e1(2,:).*e2(1,:)];
+% ch = [e1(2, :).*e2(3, :) - e1(3, :).*e2(2, :);
+%       e1(3, :).*e2(1, :) - e1(1, :).*e2(3, :);
+%       e1(1, :).*e2(2, :) - e1(2, :).*e2(1, :)];
 % sum(ch.*pl, 1) < 0
 
 
 % e1 and e2 vectors in deformed state:
-if (size(defT,3) == 1)
+if (num_def == 1)
     e1def = defT*e1 + e1;
     e2def = defT*e2 + e2;
+elseif (num_pls == 1)
+    defT_t = permute(defT, [2 1 3]);
+    defT_t = reshape(defT_t, 3, [])';
+
+    e1def = reshape(defT_t * e1, 3, []) + e1(:, ones(num_def, 1));
+    e2def = reshape(defT_t * e2, 3, []) + e2(:, ones(num_def, 1));
 else
     % Expand vectors for multiplication element-wise
-    e1r   = reshape(e1,1,[]);
-    defTe = reshape(defT,3,[],1).*[e1r; e1r; e1r];
-    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
+    e1r   = reshape(e1, 1, []);
+    defTe = reshape(defT, 3, [], 1).*[e1r; e1r; e1r];
+    defTe = defTe(:, 1:3:end) + defTe(:, 2:3:end) + defTe(:, 3:3:end);
     e1def = defTe + e1;
 
-    e2r   = reshape(e2,1,[]);
-    defTe = reshape(defT,3,[],1).*[e2r; e2r; e2r];
-    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
+    e2r   = reshape(e2, 1, []);
+    defTe = reshape(defT, 3, [], 1).*[e2r; e2r; e2r];
+    defTe = defTe(:, 1:3:end) + defTe(:, 2:3:end) + defTe(:, 3:3:end);
     e2def = defTe + e2;
 end
 
 
 % new deformed pl
-pldef = [e1def(2,:).*e2def(3,:) - e1def(3,:).*e2def(2,:);
-         e1def(3,:).*e2def(1,:) - e1def(1,:).*e2def(3,:);
-         e1def(1,:).*e2def(2,:) - e1def(2,:).*e2def(1,:)];
+pldef = [e1def(2, :).*e2def(3, :) - e1def(3, :).*e2def(2, :);
+         e1def(3, :).*e2def(1, :) - e1def(1, :).*e2def(3, :);
+         e1def(1, :).*e2def(2, :) - e1def(2, :).*e2def(1, :)];
 
      
 % normalise pldef
 normpldef = sqrt(sum(pldef.*pldef, 1));
-pldef     = pldef./normpldef([1 1 1],:);
+pldef     = pldef./normpldef([1 1 1], :);
 
 
 % check flip
@@ -100,12 +108,14 @@ pldef     = pldef./normpldef([1 1 1],:);
 % New length along original plane normal
 % (exact calculation, non-linear in the defT components)
 
-if (size(defT,3) == 1)
+if (size(defT, 3) == 1)
     fvec = defT*pl + pl;
+elseif (num_pls == 1)
+    fvec = reshape(defT_t * pl, 3, []) + pl(:, ones(num_def, 1));
 else
-    plr   = reshape(pl,1,[]);
-    defTe = reshape(defT,3,[],1).*[plr; plr; plr];
-    defTe = defTe(:,1:3:end) + defTe(:,2:3:end) + defTe(:,3:3:end);
+    plr   = reshape(pl, 1, []);
+    defTe = reshape(defT, 3, [], 1).*[plr; plr; plr];
+    defTe = defTe(:, 1:3:end) + defTe(:, 2:3:end) + defTe(:, 3:3:end);
     
     fvec = defTe + pl;   
 end