Modification in strain macros: only strain components are fitted, rotation removed.

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

zUtil_Strain2/gtStrainGrainsGeneral.m

@@ -34,16 +34,8 @@ end
 
 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');
-    
+    % Strain tensor
+    grain{ii}.strain.ST = gtStrainSingleGrain(grain{ii}, cryst, lambda, Bmat);
     
     % Principal strains - eigenvalues and eigenvectors        
     [ivec, ival] = eig(grain{ii}.strain.ST);
@@ -72,9 +64,7 @@ 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;
         

zUtil_Strain2/gtStrainSingleGrain.m

-function DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
-% GTSTRAINSINGLEGRAIN  DT = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+function ST = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
+% GTSTRAINSINGLEGRAIN  ST = 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
+% Computes the strain tensor 'ST' for a single indexed grain
+% (crystal system is irrelevant). The 6 strain 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
@@ -48,23 +48,25 @@ 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];
+% Initial guess values for fitting the 6 components
+scguess    = [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);
+boundslow  = -0.1*ones(1,6); 
+boundshigh =  0.1*ones(1,6);
 
 
 % 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
+par.typicalx = 1e-5*ones(1,6); % typical component values for calculating gradients
 
 
 
-% Find the "undeformation tensor" components 'sc' that would bring the
-% observed, deformed grain into the undeformed reference state.
+% Find the "unstrain tensor" components 'sc' that would bring the
+% observed, deformed grain into the undeformed reference state. Note that
+% rotations have no effect on the interplanar angles or distances,
+% therefore not needed for the fitting.
 % Using Matlab's built-in non-linear least square solver. 
 scopt = lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
         scguess, boundslow, boundshigh,...
@@ -78,13 +80,15 @@ scopt = lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
 %    'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
 
 
-% "Undeformation tensor"
-DTundef = reshape(scopt,3,3);
+% "Unstrain tensor"
+STundef = [scopt(1) scopt(6) scopt(5) ; ...
+           scopt(6) scopt(2) scopt(4) ; ...
+           scopt(5) scopt(4) scopt(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;
+% The real strain tensor is a kind of inverse of STundef
+ST = inv(I + STundef) - I;
 
 
 end % of main function
@@ -96,17 +100,19 @@ end % of main 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'. 
+% as a function of the 6 strain components 'sc'. These differences
+% should be eliminated by finding the right 'sc' in the fitting. 
 
 
 function dev = sfDeviation(sc, pl0, elo_ref, dotpro_ref)
 
-    DT = reshape(sc,3,3);
+    ST = [sc(1) sc(6) sc(5) ; ...
+          sc(6) sc(2) sc(4) ; ...
+          sc(5) sc(4) sc(3)];
     
     % 'pl0' is the as observed plane normals; 'pl' is 'pl0' after
-    % deformation 'DT'
-    [pl, elo] = gtStrainPlaneNormals(pl0, DT);
+    % deformation 'ST'
+    [pl, elo] = gtStrainPlaneNormals(pl0, ST);
     
     
     % Difference in elongation between computed and reference  
@@ -116,7 +122,8 @@ function dev = sfDeviation(sc, pl0, elo_ref, dotpro_ref)
     % Difference in dot products between computed and reference
     dotprodiff = abs(pl'*pl) - dotpro_ref;
     
-    % Take upper triangle of matrix 
+    % Take upper triangle of matrix (this avoids loops and seems to be 
+    % the fastest way)
     inds       = triu(true(size(dotprodiff)),1);
     dotprodiff = dotprodiff(inds);