Skip to content
Snippets Groups Projects
Commit c01f7d5e authored by Peter Reischig's avatar Peter Reischig Committed by Nicola Vigano
Browse files

Now input is loaded if not specified. Updated help.

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@595 4c865b51-4357-4376-afb4-474e03ccb993
parent a3fd6608
No related branches found
No related tags found
No related merge requests found
......@@ -4,14 +4,43 @@ function [grain, copl, out] = gtStrainGrainsGeneral(grain, cryst, ...
%
% grain = gtStrainGrainsGeneral(grain, cryst, energy, Bmat, errfilt, stdlim)
%
% Runs a strain calculation for all the grains in the cell input 'grain',
% and puts the results in the field grain{i}.strain.
% Runs a strain calculation for all the grains in the cell array 'grain{i}',
% and puts the results in the field 'grain{i}.strain'.
%
% It calls 'gtStrainSingleGrain' to do the strain calculations (see that
% help for more details).
%
% Grains with all coplanar plane normals or with outlier strain values will
% have their strain components set to NaN.
%
% It filters outliers based on their smallest or largest principal strain
% component being away from the median value of all grains more than
% 'stdlim' times the corresponding standard deviation.
%
% OPTIONAL INPUT
% grain - multpile grain data in cell array as output from Indexter
% cryst - crystallography data as in 4_grains/phase_##/index.mat file
% or parameters file
% energy - energy of the X-ray beam in keV
% Bmat - matrix to transform HKL indices into Cartesian coordinates
% errfilt - error filter for pairs (0 < errfilt <= 1); this best fraction
% of Friedel pairs will be used for error fitting;
% e.g. 0.75 meaning the best 75% based on their path distance
% from the grain center of mass; (default is 0.75)
% stdlim - limits in std of all grains' smallest and largest
% principal strain components to reject outlier grains
% (default is Inf)
%
% Either 'grain.hklsp' or 'cryst.shklfam' field is required.
%
% OUTPUT
% grain.strain - new field added containing strain data for each grain
% copl - indices of grains having plane normals that are all
% close to being coplanar
% out - grain indices that have their principal strain components
% identified as outliers
%
%
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......@@ -20,21 +49,53 @@ function [grain, copl, out] = gtStrainGrainsGeneral(grain, cryst, ...
tic
lambda = gtConvEnergyToWavelength(energy);
if ~exist('grain','var') || isempty(grain)
disp('Loading grain info from file 4_grains/phase_01/index.mat.')
index = load('4_grains/phase_01/index.mat');
grain = index.grain;
end
if ~exist('cryst','var') || isempty(cryst)
disp('Loading crystallographic info from file 4_grains/phase_01/index.mat.')
index = load('4_grains/phase_01/index.mat');
cryst = index.cryst;
end
if ~exist('energy','var') || isempty(energy)
disp('Loading beam energy from file parameters.mat.')
parameters = load('parameters.mat');
parameters = parameters.parameters;
energy = parameters.acq.energy;
end
if ~exist('Bmat','var') || isempty(Bmat)
Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
end
if ~exist('errfilt','var') || isempty(errfilt)
errfilt = 0.75;
end
if ~exist('stdlim','var') || isempty(stdlim)
stdlim = Inf;
end
lambda = gtConvEnergyToWavelength(energy);
fprintf('Using the parameters:\n')
fprintf(' Beam energy (keV): %g\n', energy)
fprintf(' Beam wavelength (A): %g\n', lambda)
fprintf(' Error filter: %g%%\n', errfilt*100)
fprintf(' Strain outliers (std): %g\n', stdlim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Strain fitting
%% Run strain fitting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Running strain fitting...')
parfor ii = 1:length(grain)
% Strain fitting (returns 'strain' structure)
......
function strain = gtStrainSingleGrain(grain, cryst, lambda, Bmat, errfilt)
% GTSTRAINSINGLEGRAIN Computes the strain tensor for a single grain.
%
% ST = gtStrainSingleGrain(grain, cryst, lambda, Bmat)
% strain = gtStrainSingleGrain(grain, cryst, lambda, Bmat, errfilt)
%
% Computes the strain tensor 'ST' for a single indexed grain
% Computes the average elastic strain tensor 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).
% computed in the same reference as the plane normals. Small deformations
% are assumed: the strain tensor is the approximated as the symmetric part
% of the deformation gradient tensor. No further approximation is used.
%
% Small deformations are assumed so that the strain tensor is considered to
% be the symmetric part of the deformation tensor, but no further
% approximation is used.
% The best fraction of Friedel pairs (determined by 'errfilt')
% based on their path distance from the center of mass of the grain is
% used for the fitting. Additional pairs are added if needed to avoid all
% plane normals being coplanar. If all plane normals are coplanar,
% there is no strain sensitivity out of that plane, and the results
% are expected to be poor.
%
% It uses the function 'lsqnonlin' from Matlab's Optimization Toolbox.
% The strain components are found by fitting the interplanar angles and
% interplanar distances from the pairs against an undeformed reference
% crystal. Currently there is no weighting of the contribution from the
% different pairs.
%
% The function 'lsqnonlin' from Matlab's Optimization Toolbox is used for
% the fitting.
%
%
% INPUT
% grain - a single grain data as from Indexter
% cryst - crystallography data as in parameters file
% lambda - X-ray beam wavelength
% Bmat - matrix to transform HKL indices into Cartesian coordinates
% (if empty, it will be calculated inside the function)
% grain - a single grain data as from Indexter
% cryst - crystallography data as in index.mat file or parameters file
% lambda - wavelength of the X-ray beam
% Bmat - matrix to transform HKL indices into Cartesian coordinates
% (if empty, it will be calculated inside the function)
% errfilt - error filter for pairs (0 < errfilt <= 1); this best fraction
% of Friedel pairs will be used for error fitting;
% e.g. 0.75 meaning the best 75% based on their path distance
% from the grain center of mass
%
% Either 'grain.hklsp' or 'cryst.shklfam' field is required.
%
% OUTPUT
% ST - strain tensor (approximated as the symmetric part of the
% deformation tensor)
% strain.
% strainT - symmetric small strain tensor
% Eps - strain components in 'strainT' as a vector
% [ST(1,1) ST(2,2) ST(3,3) ST(2,3) ST(1,3) ST(1,2)]
% pstrain - the three principal strain values in descending order
% pstraindir - the three principal directions as column vectors
% corresponding to 'pstrain'
% coplanar - true if the all the plane normals in the grain are close
% to being in the same plane (poor results expected)
% plfitted - true for the plane normals which were used for the
% fitting
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Filter out the worst pairs
%% Parameters to use
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Minimum number of Friedel pairs to be used:
minn = 5;
% Lower limit for the determinant value that estimates how coplanar the
% fitted plane normals are. Should be a value significantly larger than 0.
detlim = 1;
% Initial guess values for fitting the 6 strain components
scguess = [0 0 0 0 0 0];
% Lower and upper bounds for the strain components solution
boundslow = -0.1*ones(1,6);
boundshigh = 0.1*ones(1,6);
% For Matlab's 'lsqnonlin' optimization routine
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,6); % typical component values for calculating gradients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Filter out the worst pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plane indices used for strain fitting
plind = false(size(grain.hklind));
......@@ -47,7 +93,7 @@ if (~exist('errfilt','var') || isempty(errfilt) || (errfilt >= 1) || (errfilt <=
% Check that plane normals are not all in one plane (determinant is
% significantly larger than 0 or some tolerance)
detpl = det(pl(:,1:maxind)*pl(:,1:maxind)');
detpl = det(grain.pl*grain.pl');
if (detpl >= detlim)
coplanar = false;
else
......@@ -134,24 +180,9 @@ dotpro_ref(dotpro_ref>1) = 1; % correction for numerical errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optimization
%% Run optimization routine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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,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,6); % typical component values for calculating gradients
% 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,
......@@ -163,9 +194,9 @@ scopt = lsqnonlin( @(sc) sfDeviation(sc, pl, elo_ref, dotpro_ref),...
'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
%[scopt, resnorm, residual, exitflag, optimout, lagrange, jacobian] = ...
% lsqnonlin( @(sc) sfDeviation(sc, grain.pl, elo_ref, dotpro_ref),...
% lsqnonlin( @(sc) sfDeviation(sc, pl, elo_ref, dotpro_ref),...
% scguess, boundslow, boundshigh,...
% optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
% optimset('TolX',par.tolx, 'TolFun',par.tolfun,...
% 'TypicalX',par.typicalx, 'FunValCheck','off', 'Display','off'));
......@@ -210,7 +241,7 @@ end % of main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Sub-function
%% Sub-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gives a vector of the deviations in interplanar angles and distances
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment