diff --git a/zUtil_Crack/GtCrackAnalysis.m b/zUtil_Crack/GtCrackAnalysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..a6121c16b2a215b360ec3f0be222580a93793eea
--- /dev/null
+++ b/zUtil_Crack/GtCrackAnalysis.m
@@ -0,0 +1,750 @@
+classdef GtCrackAnalysis < handle
+% Complete procedure
+% - create a new crack analysis:
+% >> study = GtCrackAnalysis();
+%
+% - load the mesh of the crack
+% >> study.loadCrackMesh('crack_mesh.vtk');
+%
+% - set the scale factor applied when segmenting the crack (compared to DCT
+% voxel size)
+% >> study.setScaleFactor(factor);
+%
+% - load the binary volume of the crack
+% >> study.loadCrackVolume('crack_binary_volume.tif');
+%
+% - load the grain IDs from the DCT volume
+% >> study.loadGrainsVolume('dct.edf');
+%
+% - match crack mesh faces with grains ID
+% >> study.matchCrackFacesWithGrains();
+%
+% - load the Rodigues Vectors of each grain
+% >> study.loadRVectors('r_vectors.mat');
+%
+% - set the crystal system of each phase
+% >> study.setPhasesCrystalSystem({'cubic', 'hexagonal'});
+
+% - play with the functions
+% >> study.computeNormalSample();
+% >> study.computeNormalCrystal();
+% >> study.computeNormalCrystalSST();
+% >> study.computeSchmidFactors();
+
+properties
+ name; % Name of the study
+ %N; % Number of cycles
+
+ Nphases; % Number of phases
+ p_crystalSystem; % Crystal system of each phase
+ p_symm; % Symmetry operators of each phase
+ p_grainIDs; % List of grain IDs for each phase
+ p_faces; % List of grain IDs for each phase
+ p_slipSystemsFamilies;% List of slip system families of each phase
+
+ binCrack; % Binary volume of the crack
+ binCrackVoxelSize; % Voxel size of crack binary volume (µm)
+ binBBox; % BoundingBox of the binary volume of the crack
+
+ dctVol; % DCT volume containing the grain IDs
+ dctVoxelSize; % Voxel size of DCT volume
+ r_vectors; % Rodrigues vectors
+ scaleFactor; % Scale factor between DCT volume and crack volume
+
+ mesh; % Mesh for figures
+ Nvertices; % Number of vertices
+ Nfaces; % Number of faces
+ f_grainID; % Grain ID of each face
+ f_GB % Faces indices corresponding to Grain Boundaries
+
+ Ngrains; % Number of grains
+ grainIDs; % List of grain IDs
+ maxGrainID; % Maximal grain ID
+ g_phaseID; % Phase related to grain ID
+ g_orientation; % Orientation matrix g for each grain ID
+ g_schmidFactors; % Schmid factors for each grains
+
+ norm_samp; % Normal vector for each face in sample CS
+ norm_cryst; % Normal vector for each face in crystal CS
+ %norm_cryst_sst; % Normal vector for each face in crystal SST
+
+ facecol_phys; % Face orientation color in sample CS
+ facecol_cryst; % Face orientation color in crystal CS
+ facecol_cryst_sst; % Face orientation color in crystal SST
+
+ grainIDColorMap; % Color map for grain ID/label
+ grainIPFColorMap; % Color map for grain inverse pole figure
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PUBLIC METHODS
+methods (Access = public)
+ function obj = GtCrackAnalysis()
+ % GTCRACKANALYSIS/GTCRACKANALYSIS Contructor
+ obj.name = '';
+
+ obj.Nphases = 0;
+ obj.p_crystalSystem = {};
+ obj.p_symm = {};
+ obj.p_grainIDs = {};
+ obj.p_faces = {};
+ obj.p_slipSystemsFamilies = {};
+
+ obj.binCrack = [];
+ obj.binCrackVoxelSize = 0;
+ obj.binBBox = [];
+
+ obj.dctVol = [];
+ obj.dctVoxelSize = 0;
+ obj.r_vectors = [];
+ obj.scaleFactor = 0;
+
+ obj.mesh = struct;
+ obj.Nvertices = 0;
+ obj.Nfaces = 0;
+ obj.f_grainID = [];
+ obj.f_GB = [];
+
+ obj.Ngrains = 0;
+ obj.grainIDs = [];
+ obj.maxGrainID = 0;
+ obj.g_phaseID = [];
+ obj.g_orientation = {};
+ obj.g_schmidFactors = [];
+
+ obj.norm_samp = [];
+ obj.norm_cryst = [];
+ %obj.norm_cryst_sst = [];
+
+ obj.facecol_phys = [];
+ obj.facecol_cryst = [];
+ obj.facecol_cryst_sst = [];
+
+ obj.grainIDColorMap = [];
+ obj.grainIPFColorMap = [];
+ end
+
+ function delete(obj)
+ % GTCRACKANALYSIS/DELETE Destructor
+ clear obj;
+ end
+
+ function loadCrackMesh(obj, filename)
+ % GTCRACKANALYSIS/LOADCRACKMESH Load the mesh of the crack.
+ if ~exist('filename', 'var') || isempty(filename)
+ [f, p] = uigetfile('*', 'Select the file containing the mesh ' ...
+ of the crack');
+ filename = fullfile(f, p);
+ end
+ [fpath, fname, fext] = fileparts(filename);
+ if strcmp(fext, '.vtk')
+ disp([' Loading crack mesh from ' filename ' ... ']);
+ obj.mesh = gtVTKMeshReader(filename);
+ obj.Nvertices = length(obj.mesh.vertices);
+ obj.Nfaces = length(obj.mesh.faces);
+ disp([' number of vertices : ' num2str(obj.Nvertices) ]);
+ disp([' number of faces : ' num2str(obj.Nfaces) ]);
+ disp(' ');
+ else
+ gtError(['Cannot load mesh from ' fext ' file!']);
+ end
+ end
+
+ function setScaleFactor(obj, factor)
+ % GTCRACKANALYSIS/SETUPSCALEFACTOR Set the scale factor applied to the
+ % initial volume.
+ if obj.binCrackVoxelSize == 0
+ obj.binCrackVoxelSize = inputwdefaultnumeric('Please enter ' ...
+ 'the voxel size of binary crack volume.', '0.7');
+ end
+ if obj.dctVoxelSize == 0
+ obj.dctVoxelSize = inputwdefaultnumeric('Please enter the ' ...
+ 'voxel size of the DCT volume.', '1.4');
+ end
+ obj.scaleFactor = binCrackVoxelSize/dctVoxelSize;
+ %obj.scaleFactor = factor;
+ end
+
+ function loadCrackVolume(obj, filename, voxelSize)
+ % GTCRACKANALYSIS/LOADCRACKVolume Load the binary volume of the crack.
+ if ~exist('voxelSize', 'var') || isempty(voxelSize)
+ voxelSize = 0.7;
+ disp(['Warning: crack volume voxel size not given, setting it ' ...
+ 'to default: ' num2str(voxelSize)]);
+ end
+ obj.binCrackVoxelSize = voxelSize;
+
+ if ~exist('filename', 'var') || isempty(filename)
+ [f, p] = uigetfile('*', 'Select the file containing the binary ' ...
+ 'volume of the crack');
+ filename = fullfile(f, p);
+ end
+ [fpath, fname, fext] = fileparts(filename);
+ disp([' Loading crack binary volume from ' filename ' ... ']);
+ switch fext
+ case '.edf'
+ bin = logical(edf_read(filename));
+ case {'.tif', '.tiff'}
+ bin = logical(gtTIFVolReader(filename));
+ case {'.raw', '.vol'}
+ bin = logical(gtHSTVolReader(filename));
+ otherwise
+ gtError('Unsupported type of volume data!');
+ end
+ obj.binCrack = bin;
+ disp([' size of crack binary volume: ' num2str(size(bin))]);
+ disp(' ');
+
+ bboxFile = fullfile(fpath, [fname '.bbox']);
+ disp([' Loading crack binary volume from ' bboxFile ' ... ']);
+ if exist(bboxFile, 'file')
+ fid = fopen(bboxFile, 'r');
+ obj.binBBox = fscanf(fid, '[%d, %d, %d] [%d, %d, %d]')';
+ fclose(fid);
+ end
+ disp([' BoundingBox of crack volume: ' num2str(obj.binBBox)]);
+ disp(' ');
+ end
+
+ function loadGrainsVolume(obj, volume, voxelSize)
+ % GTCRACKANALYSIS/LOADGRAINSVOLUME Load the grains volume containing the
+ % grains ID.
+ if ~exist('voxelSize', 'var') || isempty(voxelSize)
+ voxelSize = 1.4;
+ disp(['Warning: DCT volume voxel size not given, setting it to ' ...
+ 'default: ' num2str(voxelSize)]);
+ end
+ obj.dctVoxelSize = voxelSize;
+
+ %if obj.scaleFactor == 0
+ % obj.setScaleFactor();
+ %end
+ if isnumeric(volume) && ndims(volume == 3)
+ dct = volume;
+ else
+ if ~exist('volume', 'var') || isempty(volume)
+ [f, p] = uigetfile('*', 'Select the file containing the ' ...
+ 'binary volume of the crack');
+ volume = fullfile(f, p);
+ end
+ [fpath, fname, fext] = fileparts(volume);
+ disp([' Loading grains volume from ' volume ' ... ']);
+ switch fext
+ case '.edf'
+ dct = uint16(edf_read(volume));
+ case '.tif'
+ dct = uint16(gtTIFVolReader(volume));
+ case {'.raw', '.vol'}
+ dct = uint16(gtHSTVolReader(volume));
+ otherwise
+ gtError('Unsupported type of volume data!');
+ end
+ end
+
+ if isempty(obj.Nphases)
+ disp(' Phases not set, please set it!');
+ obj.setPhasesCrystalSystem();
+ end
+
+ obj.grainIDs = unique(dct)';
+ obj.grainIDs = obj.grainIDs(obj.grainIDs>0);
+ obj.maxGrainID = max(obj.grainIDs);
+ obj.Ngrains = length(obj.grainIDs);
+ obj.dctVol = dct;
+ disp([' Grains volume size : ' num2str(size(obj.dctVol))]);
+ disp([' Number of grains : ' num2str(obj.Ngrains) ]);
+ disp([' Number of grain IDs : ' num2str(obj.maxGrainID) ]);
+ disp(' ');
+
+ % The multi-phase material is not completely managed right now
+ obj.g_phaseID = ones(1, obj.maxGrainID);
+ end
+
+ function setPhaseGrainID(obj)
+ % GTCRACKANALYSIS/SETPHASEGRAINID Set grain IDs of each phases
+ if isempty(obj.dctVol)
+ gtCrackAnalysis.loadGrainsVolume();
+ end
+ if isempty(obj.Nphases)
+ disp(' Phases not set, please set it!');
+ obj.setPhasesCrystalSystem();
+ end
+ obj.p_grainIDs = cell(1, obj.Nphases);
+ for iphase=1:obj.Nphases
+ obj.p_grainIDs{iphase} = uint16(obj.grainIDs);
+ end
+ end
+
+ function matchCrackFacesWithGrains(obj)
+ % GTCRACKANALYSIS/MATCHCRACKFACESWITHGRAINS Match the crack mesh faces
+ % with the grains.
+ if obj.Nvertices == 0 || obj.Nfaces == 0
+ gtCrackAnalysis.loadCrackMesh();
+ end
+ if isempty(obj.dctVol)
+ gtCrackAnalysis.loadGrainsVolume();
+ end
+
+ volSize = size(obj.dctVol);
+ obj.f_grainID = zeros(obj.Nfaces, 1, 'uint16');
+ disp([' Grain volume size is: ' num2str(volSize)]);
+
+ p1 = round(obj.mesh.vertices(obj.mesh.faces(:, 1), :) / obj.dctVoxelSize);
+ p2 = round(obj.mesh.vertices(obj.mesh.faces(:, 2), :) / obj.dctVoxelSize);
+ p3 = round(obj.mesh.vertices(obj.mesh.faces(:, 3), :) / obj.dctVoxelSize);
+
+ minIndex = ones(obj.Nfaces, 3);
+ maxIndex = volSize(ones(1, obj.Nfaces), :);
+ ok = between(p1, minIndex, maxIndex) & ...
+ between(p2, minIndex, maxIndex) & ...
+ between(p3, minIndex, maxIndex);
+
+ if any(~ok)
+ disp(['Problem with ' num2str(sum(~ok)) ' vertices:']);
+ end
+
+ for iface=1:length(obj.mesh.faces)
+ if ok(iface)
+ igrain1 = obj.dctVol(p1(iface, 1), p1(iface, 2), p1(iface, 3));
+ igrain2 = obj.dctVol(p2(iface, 1), p2(iface, 2), p2(iface, 3));
+ igrain3 = obj.dctVol(p3(iface, 1), p3(iface, 2), p3(iface, 3));
+ if igrain1 == igrain2 && igrain1 == igrain3
+ obj.f_grainID(iface) = igrain1;
+ else
+ % GB labeled as maxGrainID+1
+ obj.f_grainID(iface) = obj.maxGrainID + 1;
+ end
+ end
+ end
+ obj.f_GB = find(obj.f_grainID == obj.maxGrainID+1);
+ %obj.f_GB = obj.f_grainID == obj.maxGrainID+1;
+ disp(' ');
+ end
+
+ function loadRVectors(obj, filename)
+ % GTCRACKANALYSIS/LOADRVECTORS Load the list of R vectors and compute
+ % orientation matrices.
+ if ~exist('filename', 'var') || isempty(filename)
+ if exist('r_vectors.mat', 'file')
+ filename = fullfile(pwd, 'r_vectors.mat');
+ else
+ [f, p] = uigetfile('*.mat', 'Select the Matlab file ' ...
+ 'containing the R vectors');
+ filename = fullfile(f, p);
+ end
+ end
+ [fpath, fname, fext] = fileparts(filename);
+ if isempty(fext)
+ filename = fullfile(fpath, fname, '.mat');
+ end
+ tmp = load(filename, 'r_vectors');
+ obj.r_vectors = tmp.r_vectors(:, 2:4);
+
+ % Might be better to check that grain IDs match with array index...
+ %if length(obj.r_vectors) ~= obj.maxGrainID
+ %end
+
+ % Compute orientation matrices g
+ obj.g_orientation = gtMathsRod2RotMat(obj.r_vectors);
+ end
+
+ function setPhasesCrystalSystem(obj, crystal_systems)
+ % GTCRACKANALYSIS/SETPHASESCRYSTALSYSTEM Set the crystal system of each
+ % phase.
+ if ~exist('crystal_systems', 'var') || isempty(crystal_systems)
+ N = inputwdefaultnumeric('Number of phases', '1');
+ crystal_systems = cell(1, N);
+ for iphase=1:N
+ crystal_systems{iphase} = inputwdefault([...
+ 'Crystal systems for phase ' num2str(iphase) ''], 'cubic');
+ end
+ elseif ~iscell(crystal_systems) && ischar(crystal_systems)
+ crystal_systems = {crystal_systems};
+ end
+ obj.Nphases = length(crystal_systems);
+
+ for iphase=1:obj.Nphases
+ % Get crystal system
+ try
+ validatestring(crystal_systems{iphase}, {'cubic', ...
+ 'hexagonal', 'trigonal'});
+ catch Mexc
+ throw(Mexc);
+ end
+ syst = crystal_systems{iphase};
+
+ % Get symmetry operators
+ tmp = gtCrystGetSymmetryOperators(syst);
+ if isfield(tmp, 'g3')
+ symm{iphase} = {tmp.g3};
+ else
+ symm{iphase} = {tmp.g};
+ end
+
+ % Store phase crystal system and symmetry operators
+ obj.p_crystalSystem{iphase} = syst;
+ obj.p_symm{iphase} = symm{iphase};
+ end
+ end
+
+ function computeNormalSampleColor(obj)
+ % GTCRACKANALYSIS/COMPUTENORMALSAMPLE Compute the normal of all faces in
+ % the sample CS.
+ if obj.Nvertices == 0 || obj.Nfaces == 0
+ gtCrackAnalysis.loadCrackMesh();
+ end
+ %obj.norm_samp = zeros(obj.Nfaces, 3);
+ p1 = obj.mesh.vertices(obj.mesh.faces(:, 1), :);
+ p2 = obj.mesh.vertices(obj.mesh.faces(:, 2), :);
+ p3 = obj.mesh.vertices(obj.mesh.faces(:, 3), :);
+ %v1 = p1 - p2;
+ %v2 = p3 - p2;
+ %n = cross( v1 , v2 );
+ %n = cross(p3 - p2, p1 - p2);
+ n = cross(p1 - p2, p3 - p2);
+ n_norm = sqrt(sum(n.^2, 2));
+ obj.norm_samp = n./n_norm(:, [1 1 1]);
+ norm_zero = find(n_norm < 1.e-5);
+ obj.norm_samp(norm_zero, :) = repmat([1 0 0], size(norm_zero, 1) 1]);
+
+ % Get the vector orientation color in sample
+ rgb = gtVectorOrientationColor(obj.norm_samp);
+
+ % Label GBs in black
+ obj.facecol_phys = obj.colorGBFaces(rgb);
+ end
+
+ function computeNormalCrystalColor(obj)
+ % GTCRACKANALYSIS/COMPUTENORMALCRYST Compute the normal of all faces in
+ % the crystal CS.
+ if isempty(obj.norm_samp)
+ obj.computeNormalSampleColor();
+ end
+ if isempty(obj.r_vectors)
+ obj.loadRVectors();
+ end
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+
+ obj.norm_cryst = zeros(obj.Nfaces, 3);
+ for iphase=1:obj.Nphases
+ for igrain=obj.p_grainIDs{iphase}
+ g_facesID = find(obj.f_grainID == igrain);
+ if g_facesID
+ g = obj.g_orientation{igrain};
+ obj.norm_cryst(g_facesID, :) = ...
+ gtVectorLab2Cryst(obj.norm_samp(g_facesID, :), g);
+ end
+ end
+ end
+
+ %%%%% TEST
+ % Turn all normals so that they lay in upper half of sphere
+ %lower_sphere = obj.norm_cryst(:,3) < 0;
+ %obj.norm_cryst(lower_sphere,:) = -obj.norm_cryst(lower_sphere,:);
+ %%%%% TEST
+
+ % Get the vector orientation color in crystal
+ rgb = gtVectorOrientationColor(obj.norm_cryst);
+
+ % Label GBs in black
+ obj.facecol_cryst = obj.colorGBFaces(rgb);
+ end
+
+ function computeNormalCrystalSSTColor(obj)
+ % GTCRACKANALYSIS/COMPUTENORMALCRYSTSST Compute the normal of all faces
+ % in the crystal CS SST.
+ if isempty(obj.norm_cryst)
+ obj.computeNormalCrystalColor();
+ end
+
+ obj.facecol_cryst_sst = zeros(obj.Nfaces, 3);
+ for iphase=1:obj.Nphases
+ for igrain=obj.p_grainIDs{iphase}
+ g_facesID = find(obj.f_grainID == igrain);
+ if length(g_facesID)
+ obj.facecol_cryst_sst(g_facesID, :) = gtCrystVector2SST(...
+ obj.norm_cryst(g_facesID, :), ...
+ obj.p_crystalSystem{iphase}, obj.p_symm{iphase}, false);
+ %obj.norm_cryst(g_facesID, :), ...
+ %obj.p_crystalSystem{iphase}, obj.p_symm{iphase});
+ end
+ end
+ end
+ end
+
+ function computeSchmidFactors(obj, load)
+ % GTCRACKANALYSIS/COMPUTESCHMIDFACTORS Compute Schmid factors depending on
+ % the given load direction.
+ if isempty(obj.r_vectors)
+ obj.loadRVectors();
+ end
+ if ~exist('load', 'var') || isempty(load)
+ load = [0 0 1];
+ end
+ load = load / norm(load);
+ for iphase=1:obj.Nphases
+ grainIDs = obj.p_grainIDs{iphase};
+ loadCryst = gtVectorLab2Cryst(load, {obj.g_orientation{grainIDs}});
+ slipFamilies = obj.p_slipSystemsFamilies{iphase};
+ slipSystems = gtGetSlipSystems(slipFamilies);
+ for ifam=1:length(slipFamilies)
+ syst = getfield(slipSystems, slipFamilies{ifam});
+ schmid = zeros(1, length(syst));
+ for igrain=1:obj.maxGrainID
+ for isyst=1:length(syst)
+ theta = dot(F, syst(isyst).n');
+ kappa = dot(F, syst(isyst).l');
+ schmid(isyst) = abs(theta * kappa);
+ end
+ end
+ setfield(obj.g_schmidFactors(iphase), slipFamilies{ifam}, schmid);
+ end
+ end
+ end
+
+ function [hf, ha, hp] = drawCrackEdgesHeight(obj)
+ % GTCRACKANALYSIS/DRAWCRACKFACESHEIGHT
+ disp(' Drawing faces height map ...');
+ disp(' ');
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+ heightColor = squeeze(obj.mesh.vertices(:, 3));
+ [hf, ha, hp] = obj.drawCrackEdgesColor(heightColor, 'Name', ...
+ 'Height of the crack vertices');
+ end
+
+ function [hf, ha, hp] = drawCrackGrainLabel(obj)
+ % GTCRACKANALYSIS/DRAWCRACKGRAINLABEL Draw the crack faces grain ID/label
+ disp(' Drawing grain label map ...');
+ disp(' ');
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+ if isempty(obj.grainIDColorMap)
+ obj.createGrainIDColorMap();
+ end
+ [hf, ha, hp] = obj.drawCrackFacesColor(obj.f_grainID, false, ...
+ 'Name', 'Grain labels', 'Colormap', obj.grainIDColorMap);
+ end
+
+ function [hf, ha, hp] = drawCrackGrainIPF(obj)
+ % GTCRACKANALYSIS/DRAWCRACKGRAINIPF Draw the crack faces grain inverse
+ % pole figure
+ disp(' Drawing grain inverse pole figure map ...');
+ disp(' ');
+ if isempty(obj.grainIPFColorMap)
+ obj.createGrainIPFColorMap();
+ end
+ [hf, ha] = obj.drawCrackFacesColor(obj.f_grainID, false, 'Name', ...
+ 'Grain inverse pole figure', 'Colormap', obj.grainIPFColorMap);
+ end
+
+ function [hf, ha, hp] = drawCrackFacesHeight(obj)
+ % GTCRACKANALYSIS/DRAWCRACKFACESHEIGHT
+ disp(' Drawing faces height map ...');
+ disp(' ');
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+ f_height = squeeze(mean(reshape(obj.mesh.vertices(...
+ obj.mesh.faces(:, :), 3), [obj.Nfaces 3 1]), 2));
+ Ncolor = 50;
+ colorGB = (1-1/Ncolor) * min(f_height) - max(f_height)/Ncolor;
+ heightColor = obj.colorGBFaces(f_height, colorGB);
+ cmap = [0 0 0 ; jet(Ncolor)];
+ [hf, ha, hp] = obj.drawCrackFacesColor(heightColor, false, ...
+ 'ColorMap', cmap, ...
+ 'Name', 'Height of the crack faces');
+ set(hp, 'CDataMapping', 'scaled');
+ end
+
+ function [hf, ha, hp] = drawCrackFacesOrientationSample(obj)
+ % GTCRACKANALYSIS/DRAWCRACKFACESORIENTATIONSAMPLE Draw the crack faces
+ % with the sample orientation color code.
+ disp(' Drawing faces orientation in sample CS ...');
+ disp(' ');
+ if isempty(obj.facecol_phys) || isempty(obj.norm_samp)
+ obj.computeNormalSampleColor();
+ end
+ [hf, ha, hp] = obj.drawCrackFacesColor(obj.facecol_phys, false, ...
+ 'Name', 'Orientation of the crack faces in the sample');
+ end
+
+ function [hf, ha, hp] = drawCrackFacesOrientationCrystal(obj)
+ % GTCRACKANALYSIS/DRAWCRACKFACESORIENTATIONCRYSTAL Draw the crack faces
+ % with the crystal orientation color code.
+ disp(' Drawing faces orientation in crystal CS ...');
+ disp(' ');
+ if isempty(obj.facecol_cryst) || isempty(obj.norm_cryst)
+ obj.computeNormalCrystalColor();
+ end
+ [hf, ha, hp] = obj.drawCrackFacesColor(obj.facecol_cryst, false, ...
+ 'Name', 'Orientation of the crack faces in the crystal');
+ end
+
+ function [hf, ha, hp] = drawCrackFacesOrientationCrystalSST(obj, sat)
+ % GTCRACKANALYSIS/DRAWCRACKFACESORIENTATIONCRYSTALSST Draw the crack
+ % faces with the crystal SST orientation color code.
+ if ~exist('sat', 'var') || isempty(sat)
+ saturation = true;
+ end
+ disp(' Drawing faces orientation in crystal SST ...');
+ disp(' ');
+ if isempty(obj.facecol_cryst_sst) || isempty(obj.norm_cryst) || ...
+ isempty(obj.norm_cryst)
+ obj.computeNormalCrystalSSTColor();
+ end
+ [hf, ha, hp] = obj.drawCrackFacesColor(obj.facecol_cryst_sst, sat, ...
+ 'Name', 'Orientation of the crack faces in the crystal SST');
+ end
+
+end % end of public methods
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PROTECTED METHODS
+methods (Access = protected)
+ function x = colorGBFaces(obj, x, color)
+ % GTCRACKANALYSIS/COLORGBFACES Change the color of GB faces to black or
+ % white.
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+ Ncomp = size(x, 2);
+ if ~exist('color', 'var') || isempty(color)
+ color = zeros(1, Ncomp);
+ elseif ~size(color) == [1 Ncomp]
+ gtError('GtCrackAnalysis:colorGBFaces', ...
+ 'Wrong RGB color vector size!');
+ end
+ if Ncomp == 1
+ x(obj.f_GB) = color(ones(1, size(obj.f_GB, 1)));
+ elseif Ncomp == 3
+ x(obj.f_GB, :) = color(ones(1, size(obj.f_GB, 1)), :);
+ end
+ end
+
+ function createGrainIDColorMap(obj)
+ % GTCRACKANALYSIS/CREATEGRAINIDCOLORMAP Creates a colormap for grain
+ % labels.
+ cmap = gtRandCmap(obj.maxGrainID);
+ cmap(1, :) = [ 1 1 1 ];
+ obj.grainIDColorMap = [ cmap ; 0 0 0 ];
+ end
+
+ function createGrainIPFColorMap(obj, LD)
+ % GTCRACKANALYSIS/CREATEGRAINIPFCOLORMAP Creates a colormap for grain
+ % inverse pole figure for LD direction.
+ if isempty(obj.f_grainID)
+ obj.matchCrackFacesWithGrains();
+ end
+ if isempty(obj.r_vectors)
+ obj.loadRVectors();
+ end
+ if ~exist('LD', 'var') || isempty(LD)
+ LD = [ 0 0 1 ];
+ end
+
+ obj.setPhaseGrainID();
+ map = zeros(obj.maxGrainID, 3);
+ for iphase=1:obj.Nphases
+ grainIDs = obj.p_grainIDs{iphase};
+ LDc = gtVectorLab2Cryst(LD, {obj.g_orientation{grainIDs}});
+ map(grainIDs, :) = gtCrystVector2SST(LDc, ...
+ obj.p_crystalSystem{iphase}, obj.p_symm{iphase});
+ end
+
+ obj.grainIPFColorMap = [ 1 1 1 ; map ; 0 0 0 ];
+ end
+
+ function [hf, ha, hp] = drawCrackEdgesColor(obj, edgeColor, varargin)
+ % GTCRACKANALYSIS/DRAWCRACKORIENTATION Draw the crack faces orientation
+ % with given color.
+ par.Color = 'w';
+ par.ColorMap = [];
+ par.Name = [];
+ par = parse_pv_pairs(par, varargin);
+
+ % Create figure
+ hf = figure();
+ ha = axes();
+
+ % Create a patch from the crack mesh
+ hp = patch(obj.mesh);
+ set(hp, 'FaceVertexCData', edgeColor, ...
+ 'FaceColor', 'interp', ...
+ 'CDataMapping', 'scaled', ...
+ 'EdgeColor', 'interp');
+ %'LineWidth', 5 ...)
+
+ minColor = min(edgeColor);
+ maxColor = max(edgeColor);
+
+ shading(ha, 'flat');
+ axis equal;
+ axis vis3d;
+ set(ha, 'XTickLabel', '', 'YTickLabel', '', 'ZTickLabel', '', ...
+ 'Xcolor', 'w', 'Ycolor', 'w', 'Zcolor', 'w', ...
+ 'GridLineStyle', 'none', ...
+ 'CLim', [minColor maxColor]);
+ if exist('title', 'var') && ~isempty(title)
+ set(hf, 'Name', title);
+ end
+
+ parFields = fieldnames(par);
+ for i=1:numel(parFields)
+ val = par.(parFields{i});
+ if ~isempty(val)
+ set(hf, parFields{i}, val);
+ end
+ end
+ end
+
+ function [hf, ha, hp] = drawCrackFacesColor(obj, faceColor, sat, varargin)
+ % GTCRACKANALYSIS/DRAWCRACKORIENTATION Draw the crack faces orientation
+ % with given color.
+ % Set and get figure parameters
+ par.Color = 'w';
+ par.ColorMap = [];
+ par.Name = [];
+ par = parse_pv_pairs(par, varargin);
+
+ if sat
+ ok = any(faceColor ~= 0, 2);
+ saturated = max(faceColor(ok, :), [], 2);
+ faceColor(ok, :) = faceColor(ok, :)./saturated(:, [1 1 1]);
+ end
+
+ % Create figure
+ hf = figure();
+ ha = axes();
+
+ % Create a patch from the crack mesh
+ hp = patch(obj.mesh);
+ set(hp, 'FaceVertexCData', faceColor);
+
+ shading(ha, 'flat');
+ axis equal;
+ axis vis3d;
+ set(ha, 'XTickLabel', '', 'YTickLabel', '', 'ZTickLabel', '', ...
+ 'Xcolor', 'w', 'Ycolor', 'w', 'Zcolor', 'w', ...
+ 'GridLineStyle', 'none');
+ if exist('title', 'var') && ~isempty(title)
+ set(hf, 'Name', title);
+ end
+
+ parFields = fieldnames(par);
+ for i=1:numel(parFields)
+ val = par.(parFields{i});
+ if ~isempty(val)
+ set(hf, parFields{i}, val);
+ end
+ end
+ shg;
+ end
+
+end % end of protected methods
+
+end % end of class
diff --git a/zUtil_Cryst/gtGetSlipSystems.m b/zUtil_Cryst/gtGetSlipSystems.m
new file mode 100644
index 0000000000000000000000000000000000000000..e51aad2754c8acb29e8df60fa190407bad12d8d0
--- /dev/null
+++ b/zUtil_Cryst/gtGetSlipSystems.m
@@ -0,0 +1,107 @@
+function slipSystems = gtGetSlipSystems(slipFamilies)
+% Based on Z-set slip systems definition
+
+slipSystems = struct;
+
+for ifam:1:length(slipFamilies)
+ switch slipFamilies{ifam}
+ case {'fcc', 'FCC', 'octahedral'}
+
+ fcc(1).n = [ 1 1 1]; fcc(1).l = [-1 0 1]; % Bd / B4
+ fcc(2).n = [ 1 1 1]; fcc(2).l = [ 0 -1 1]; % Ba / B2
+ fcc(3).n = [ 1 1 1]; fcc(3).l = [-1 1 0]; % Bc / B5
+
+ fcc(4).n = [ 1 -1 1]; fcc(4).l = [-1 0 1]; % Db / D4
+ fcc(5).n = [ 1 -1 1]; fcc(5).l = [ 0 1 1]; % Dc / D1
+ fcc(6).n = [ 1 -1 1]; fcc(6).l = [ 1 1 0]; % Da / D6
+
+ fcc(7).n = [-1 1 1]; fcc(7).l = [ 0 -1 1]; % Ab / A2
+ fcc(8).n = [-1 1 1]; fcc(8).l = [ 1 1 0]; % Ad / A6
+ fcc(9).n = [-1 1 1]; fcc(9).l = [ 1 0 1]; % Ac / A3
+
+ fcc(10).n = [ 1 1 -1]; fcc(10).l = [-1 1 0]; % Cb / C5
+ fcc(11).n = [ 1 1 -1]; fcc(11).l = [ 1 0 1]; % Ca / C3
+ fcc(12).n = [ 1 1 -1]; fcc(12).l = [ 0 1 1]; % Cd / C1
+
+ slipSystems.fcc = fcc;
+
+ case {'bcc', 'BCC'}
+
+ bcc(1).n = [-1 0 1]; bcc(1).l = [ 1 1 1]; % Bd / B4
+ bcc(2).n = [ 0 -1 1]; bcc(2).l = [ 1 1 1]; % Ba / B2
+ bcc(3).n = [-1 1 0]; bcc(3).l = [ 1 1 1]; % Bc / B5
+
+ bcc(4).n = [-1 0 1]; bcc(4).l = [ 1 -1 1]; % Db / D4
+ bcc(5).n = [ 0 1 1]; bcc(5).l = [ 1 -1 1]; % Dc / D1
+ bcc(6).n = [ 1 1 0]; bcc(6).l = [ 1 -1 1]; % Da / D6
+
+ bcc(7).n = [ 0 -1 1]; bcc(7).l = [-1 1 1]; % Ab / A2
+ bcc(8).n = [ 1 1 0]; bcc(8).l = [-1 1 1]; % Ad / A6
+ bcc(9).n = [ 1 0 1]; bcc(9).l = [-1 1 1]; % Ac / A3
+
+ bcc(10).n = [-1 1 0]; bcc(10).l = [ 1 1 -1]; % Cb / C5
+ bcc(11).n = [ 1 0 1]; bcc(11).l = [ 1 1 -1]; % Ca / C3
+ bcc(12).n = [ 0 1 1]; bcc(12).l = [ 1 1 -1]; % Cd / C1
+
+ slipSystems.bcc = bcc;
+
+ case {'bcc112', 'BCC112'}
+
+ bcc112(1).n = [ 1 1 2]; bcc112(1).l = [ 1 1 -1];
+ bcc112(2).n = [-1 1 2]; bcc112(2).l = [ 1 -1 1];
+ bcc112(3).n = [ 1 -1 2]; bcc112(3).l = [ -1 1 1];
+
+ bcc112(4).n = [ 1 1 -2]; bcc112(4).l = [ 1 1 1];
+ bcc112(5).n = [ 1 2 1]; bcc112(5).l = [ 1 -1 1];
+ bcc112(6).n = [-1 2 1]; bcc112(6).l = [ 1 1 -1];
+
+ bcc112(7).n = [ 1 -2 1]; bcc112(7).l = [ 1 1 1];
+ bcc112(8).n = [ 1 2 -1]; bcc112(8).l = [ -1 1 1];
+ bcc112(9).n = [ 2 1 1]; bcc112(9).l = [ -1 1 1];
+
+ bcc112(10).n = [-2 1 1]; bcc112(10).l = [ 1 1 1];
+ bcc112(11).n = [ 2 -1 1]; bcc112(11).l = [ 1 1 -1];
+ bcc112(12).n = [ 2 1 -1]; bcc112(12).l = [ 1 -1 1];
+
+ slipSystems.bcc112 = bcc112;
+
+ case {'bcc123', 'BCC123'}
+
+ bcc123(1).n = [ 1 2 3]; bcc123(1).l = [ 1 1 -1];
+ bcc123(2).n = [-1 2 3]; bcc123(2).l = [ 1 -1 1];
+ bcc123(3).n = [ 1 -2 3]; bcc123(3).l = [-1 1 1];
+
+ bcc123(4).n = [ 1 2 -3]; bcc123(4).l = [ 1 1 1];
+ bcc123(5).n = [ 3 1 2]; bcc123(5).l = [-1 1 1];
+ bcc123(6).n = [-3 1 2]; bcc123(6).l = [ 1 1 1];
+
+ bcc123(7).n = [ 3 -1 2]; bcc123(7).l = [ 1 1 -1];
+ bcc123(8).n = [ 3 1 -2]; bcc123(8).l = [ 1 -1 1];
+ bcc123(9).n = [ 2 3 1]; bcc123(9).l = [ 1 -1 1];
+
+ bcc123(10).n = [-2 3 1]; bcc123(10).l = [ 1 1 -1];
+ bcc123(11).n = [ 2 -3 1]; bcc123(11).l = [ 1 1 1];
+ bcc123(12).n = [ 2 3 -1]; bcc123(12).l = [-1 1 1];
+
+ bcc123(13).n = [ 1 3 2]; bcc123(13).l = [ 1 -1 1];
+ bcc123(14).n = [-1 3 2]; bcc123(14).l = [ 1 1 -1];
+ bcc123(15).n = [ 1 -3 2]; bcc123(15).l = [ 1 1 1];
+
+ bcc123(16).n = [ 1 3 -2]; bcc123(16).l = [-1 1 1];
+ bcc123(17).n = [ 2 1 3]; bcc123(17).l = [ 1 1 -1];
+ bcc123(18).n = [-2 1 3]; bcc123(18).l = [ 1 -1 1];
+
+ bcc123(19).n = [ 2 -1 3]; bcc123(19).l = [-1 1 1];
+ bcc123(20).n = [ 2 1 -3]; bcc123(20).l = [ 1 1 1];
+ bcc123(21).n = [ 3 2 1]; bcc123(21).l = [-1 1 1];
+
+ bcc123(22).n = [-3 2 1]; bcc123(22).l = [ 1 1 1];
+ bcc123(23).n = [ 3 -2 1]; bcc123(23).l = [ 1 1 -1];
+ bcc123(24).n = [ 3 2 -1]; bcc123(24).l = [ 1 -1 1];
+
+ slipSystems.bcc123 = bcc123;
+
+ otherwise
+ gtError('gtGetSlipSystems:wrong_input', ...
+ 'Unsupported slip systems family!');
+end