From 85f4513b219a79c8e3a695f920de79c641fe42ae Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Mon, 14 Dec 2015 15:52:17 +0100
Subject: [PATCH] gtStrainPlaneNormals: optimized for deformation of single
plane normal
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Strain2/gtStrainPlaneNormals.m | 74 ++++++++++++++++------------
1 file changed, 42 insertions(+), 32 deletions(-)
diff --git a/zUtil_Strain2/gtStrainPlaneNormals.m b/zUtil_Strain2/gtStrainPlaneNormals.m
index 92f5e5fc..0f34ffb7 100644
--- a/zUtil_Strain2/gtStrainPlaneNormals.m
+++ b/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
--
GitLab