Added basic strain functions for general geometry and any crystal system.

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@446 4c865b51-4357-4376-afb4-474e03ccb993

zUtil_Strain2/gtStrainGrainsGeneral.m

+function grain = gtStrainGrainsGeneral(grain, cryst, energy, Bmat, stdlim)
+% GTSTRAINGRAINSGENERAL grain = gtStrainGrainsGeneral(grain, cryst, energy,
+%                                                     Bmat, stdlim)
+%
+% Runs a strain calculation for all the grains in the cell input 'grain',
+% and puts the results in the field grain{i}.strain.
+% Needs the 'cryst' data as in the output of Indexter
+% (in file 4_grains/phase_##/index.mat) and the beam energy in keV.
+% It filters outliers based on their principal strains being away from the 
+% median value of all grains more than 'stdlim' times its std.
+% The fields 'Bmat' and 'stdlim' are optional. 
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set input
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+tic
+
+lambda = gtConvEnergyToWavelength(energy);
+
+if ~exist('Bmat','var') || isempty(Bmat)
+    Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+end
+
+if ~exist('stdlim','var') || isempty(stdlim)
+    stdlim = Inf;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Loop through grains
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+parfor ii = 1:length(grain)
+    
+    % Deformation tensor
+    DT = gtStrainSingleGrain(grain{ii}, cryst, lambda, Bmat);
+    grain{ii}.strain.DT = DT;
+    
+    % Strain tensor - symmetric part
+    grain{ii}.strain.ST = 0.5*(DT + DT');    
+    
+    % Rotation tensor - skew part
+    grain{ii}.strain.ROT = 0.5*(DT - DT');
+    
+    
+    % Principal strains - eigenvalues and eigenvectors        
+    [ivec, ival] = eig(grain{ii}.strain.ST);
+    
+    [ivalsort, ivalinds] = sort([ival(1,1), ival(2,2), ival(3,3)], 2, 'descend');
+    
+    % principal strain values ordered
+    grain{ii}.strain.pstrain = ivalsort;
+    
+    % principal strain directions in columns
+    grain{ii}.strain.pstraindir = ivec(:, ivalinds);
+    
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Filter outliers
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+med_ps = median(gtIndexAllGrainValues(grain,'strain','pstrain',1));
+std_ps = std(gtIndexAllGrainValues(grain,'strain','pstrain',1));
+
+
+parfor ii = 1:length(grain)
+    
+    if (grain{ii}.strain.pstrain(1) > med_ps + stdlim*std_ps) || ...
+       (grain{ii}.strain.pstrain(1) < med_ps - stdlim*std_ps)   
+        
+        grain{ii}.strain.DT(:)         = NaN;
+        grain{ii}.strain.ST(:)         = NaN;
+        grain{ii}.strain.ROT(:)        = NaN;
+        grain{ii}.strain.pstrain(:)    = NaN;
+        grain{ii}.strain.pstraindir(:) = NaN;
+        
+    end
+    
+end
+
+
+toc
+
+
+end % of function
\ No newline at end of file

zUtil_Strain2/gtStrainPlaneNormals.m

+function [pldef, drel] = gtStrainPlaneNormals(pl, G)
+% FUNCTION [pldef, drel] = gtStrainPlaneNormals(pl, G)
+%
+% Returns the new orientations of plane normals 'pl' after a deformation 
+% described by the deformation tensor 'G' (9 distinct components).
+% Also returns the elongations 'drel' due to the deformation in the 
+% directions of the plane normals. In a crystal this equals to the ratios 
+% of the d-spacings after and before the deformation.
+% The results are exact, no linear approximation is used.
+% The 'pl' and 'G' coordinates must be given in the same reference frame.
+%
+
+
+% Get an arbitrary vector e1 perpendicular to n:
+[~, maxcomp] = max(abs(pl), [], 1);
+
+
+ee = zeros(3,size(pl,2)*3);
+
+zero = zeros(1,size(pl,2));
+
+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);
+
+
+% 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],:);
+% 
+% ne1 = sqrt(sum(e1.*e1, 1));
+% e1n = e1./ne1([1 1 1],:);
+%  
+% 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,:)];
+% sum(ch.*pl, 1) < 0
+
+
+% e-s in deformed state:
+e1def = G*e1 + e1;
+e2def = G*e2 + e2;
+
+
+% 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,:)];
+
+     
+% normalise pldef
+normpldef = sqrt(sum(pldef.*pldef, 1));
+pldef     = pldef./normpldef([1 1 1],:);
+
+
+% check flip
+% if sum(pldef.*pl, 1) < 0
+% 	pldef = ...-pldef...;
+% end
+%sum(pldef.*pl, 1) >= 0
+
+
+
+% new relative d-spacing (exact calculation, non-linear in the G components)
+fvec = pl + G*pl;   % change along original plane normal
+
+% new normalised d-spacing
+drel = sum(pldef.*fvec, 1);
+
+
+
+end % of function
+

zUtil_Strain2/gtStrainSingleGrain.m

+function DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+% GTSTRAINSINGLEGRAIN  DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+% 
+% Computes the deformation tensor 'DT' for a single indexed grain
+% (crystal system is irrelevant). The 9 deformation matrix components are
+% computed in the same reference as the plane normals.
+% The deformation components are found by fitting all the interplanar 
+% angles and distances observed (no weighting of the equations). No linear
+% approximation is used.
+% If 'Bmat' is empty, it will be calculated inside the function. 
+%
+
+
+% All measured signed hkl-s in the grain 
+shkl = [];
+for ii = 1:length(grain.shklind)
+    shkl(:,ii) = cryst.shklfam{grain.hklind(ii)}(grain.shklind(ii),:)';
+end
+
+
+% Measured d-spacings
+dsp = 0.5*lambda./sind(grain.theta);
+
+% D-spacings in the undeformed reference crystal
+dsp_ref = cryst.dspacing(grain.hklind);
+
+% Required relative elongations to achieve unstrained state from deformed state
+elo_ref = dsp_ref./dsp;
+
+
+% Get B matrix for hkl -> Cartesian coordinate transforms
+if isempty(Bmat)
+    Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
+end
+
+% Cartesian coordinates of plane normals in the unstrained reference crystal
+% (output plane normals are normalized)
+pl_ref = gtCrystHKL2Cartesian(shkl, Bmat);
+
+% All combinations of plane normal dot products in the unstrained reference crystal
+dotpro_ref = abs(pl_ref'*pl_ref);
+dotpro_ref(dotpro_ref>1) = 1;  % correction for numerical errors
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Optimization
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Initial guess values for fitting the 9 components
+scguess    = [0 0 0 0 0 0 0 0 0];
+
+% Lower and upper bounds for the components solution
+boundslow  = -0.1*ones(1,9); 
+boundshigh =  0.1*ones(1,9);
+
+
+% Optimization parameters
+par.tolx     = 1e-6;  % precision for components (1e-6 is reasonable and fast)
+par.tolfun   = 1e-12; % optimization function target value (should be very low)
+par.typicalx = 1e-5*ones(1,9); % typical component values for calculating gradients
+
+
+
+% Find the "undeformation tensor" components 'sc' that would bring the
+% observed, deformed grain into the undeformed reference state.
+% Using Matlab's built-in non-linear least square solver. 
+scopt = lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
+        scguess, boundslow, boundshigh,...
+        optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
+        'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
+
+%[scopt, resnorm, residual, exitflag, optimout, lagrange, jacobian] = ...
+%    lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
+%    scguess, boundslow, boundshigh,...
+%   optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
+%    'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
+
+
+% "Undeformation tensor"
+DTundef = reshape(scopt,3,3);
+
+I = [1 0 0; 0 1 0; 0 0 1];
+
+% The real deformation tensor is a kind of inverse of DTundef
+DT = inv(I + DTundef) - I;
+
+
+end % of main function
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Sub-function
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Gives a vector of the deviations in interplanar angles and distances 
+% as a function of the deformation components 'sc'. These differences
+% should be eliminated by finding the right 'sc'. 
+
+
+function dev = sfDeviation(sc, pl0, elo_ref, dotpro_ref)
+
+    DT = reshape(sc,3,3);
+    
+    % 'pl0' is the as observed plane normals; 'pl' is 'pl0' after
+    % deformation 'DT'
+    [pl, elo] = gtStrainPlaneNormals(pl0, DT);
+    
+    
+    % Difference in elongation between computed and reference  
+    elodiff = elo - elo_ref;
+    
+    
+    % Difference in dot products between computed and reference
+    dotprodiff = abs(pl'*pl) - dotpro_ref;
+    
+    % Take upper triangle of matrix 
+    inds       = triu(true(size(dotprodiff)),1);
+    dotprodiff = dotprodiff(inds);
+    
+    
+    % All differences in one vector
+    dev = [elodiff(:); dotprodiff];
+
+end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+