Schmid Factor Code

1,125 views
Skip to first unread message

grandr...@gmail.com

unread,
Jul 12, 2016, 1:01:09 PM7/12/16
to mtex...@googlegroups.com
I wanted to create a code that would take an EBSD map, and when provided with a list of slip and twinning systems, and  the relative CRSS (critical resolved shear stress) to activate them determine which system is active and what the relative stress on that system is.  I have provided this below for anyone who is interested, set up for magnesium.  Data is assumed to be in the 'ebsd' variable.  This is undoubtedly not the most computationally efficient method  to use, but I wanted to be able to compare the systems against each other and check outputs at each stage of the process.  If anyone finds an error, please let me know.  The code to plot the figures is provided at the end with comments.
Jessica

EDIT: corrected Aug 11 2016 - error in twinning calculations.  There should be no use of antipodal anywhere in that section.  Error caused the opposite direction of force to give different schmid factors for twinning, when tension along +xvector and -xvector should give identical results.  Code below updated, some comments edited

%% General inputs
%Note that to simplify things, the schmid factor output will always be >0
%for slip and this code does not distinguish if the slip is in the negative
%direction.  For twinning, schmid factor output will always be -90 to +90,
%as it is directional.
%note that for magnesium a plane and the normal will not necessarily have
%the same index values (unlike cubic systems), so use the ,'hkl' suffix for
%plane normals and the 'uvw' suffix for directions.
%%define system crystallography

cs
=ebsd('Magnesium').CS;

%define tensile axis or for an arbitrary vector use r = vector3d(1,-1,0);
r
=xvector;

%define  CRSS of the basal, prismatic, pyramidal, extension twin so that the schmid
%factor can be scaled relatively between the systems
%based on Figure 11, Chapuis, A. & Driver, J. H. Temperature dependency of slip and twinning in plane strain compressed magnesium single crystals Acta Mater., 2011, 59, 1986-1994 in MPa. 


BasalFactor=4/4;

PrismaticFactor=18/4;

%type II pyramidal
PyramidalFactor=40/4;

%{10-12} extension twinning
ExtTwinFactor=11/4;


%% Basal slip calculations

%define slip plane and directions
% normal to the slip plane
m1
= Miller(0,0,0,1,cs,'hkl');
% slip direction in the slip plane
n1
= Miller(1,1,-2,0,cs,'uvw');
n2
= Miller(-2,1,1,0,cs,'uvw');
n3
= Miller(1,-2,1,0,cs,'uvw');

% calculate angle between force and plane normal for each point in the map
BasalPlaneNormal= ebsd('Magnesium').orientations*m1;
theta
=cos(angle(r, BasalPlaneNormal,'antipodal'));
%to check the angles between BasalPlaneNormal and r, use max(angle(r, BasalPlaneNormal,'antipodal')/degree)
%also check min, values should range from 0-90 degrees  We use antipodal
%symmetry to enforce this

%calculate the orientation of the 3 slip directions
BasalSlip1=ebsd('Magnesium').orientations*n1;
BasalSlip2=ebsd('Magnesium').orientations*n2;
BasalSlip3=ebsd('Magnesium').orientations*n3;
%calculate the angle between the slip and the force
lmda1
=cos(angle(r,BasalSlip1,'antipodal'));
lmda2
=cos(angle(r,BasalSlip2,'antipodal'));
lmda3
=cos(angle(r,BasalSlip3,'antipodal'));
%to check the angles between BasalSlip1 and r, use max(angle(r, BasalSlip1,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,


% Now calculate the schmid factor
tau1
=theta.*lmda1;
tau2
=theta.*lmda2;
tau3
=theta.*lmda3;
%combine the variables and calculate the minimum value and the index where
%it occurs. max(BasalArray) to test.  The first 4 columns should be 0 to
%0.5, and the final should be an integer between 1 and 3
BasalArray=[tau1,tau2,tau3];
[val,I]=max(BasalArray,[],2);
BasalArray=[tau1,tau2,tau3,val,I];

%% Prismatic slip calculations

%define slip plane and directions
% normal to the slip plane
m4
= Miller(1,0,-1,0,cs,'hkl');
m5
= Miller(0,1,-1,0,cs,'hkl');
m6
= Miller(-1,1,0,0,cs,'hkl');
% slip direction in the slip plane
n4
= Miller(-1,2,-1,0,cs,'uvw');
n5
= Miller(2,-1,-1,0,cs,'uvw');
n6
= Miller(-1,-1,2,0,cs,'uvw');

%to check the angle use test=angle(m5,n5,'noSymmetry')/degree

% calculate Schmid factors for each point in the map
%angle between force and normal
PrismaticPlaneNormal4= ebsd('Magnesium').orientations*m4;
PrismaticPlaneNormal5= ebsd('Magnesium').orientations*m5;
PrismaticPlaneNormal6= ebsd('Magnesium').orientations*m6;
theta4
=cos(angle(r, PrismaticPlaneNormal4,'antipodal'));
theta5
=cos(angle(r, PrismaticPlaneNormal5,'antipodal'));
theta6
=cos(angle(r, PrismaticPlaneNormal6,'antipodal'));
%to check the angles between PrismaticPlaneNormal4 and r, use max(angle(r, PrismaticPlaneNormal4,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,


%calcualte the orientation of the 3 slip directions
PrismaticSlip4=ebsd('Magnesium').orientations*n4;
PrismaticSlip5=ebsd('Magnesium').orientations*n5;
PrismaticSlip6=ebsd('Magnesium').orientations*n6;
%calculate the angle between the slip and the force
lmda4
=cos(angle(r,PrismaticSlip4,'antipodal'));
lmda5
=cos(angle(r,PrismaticSlip5,'antipodal'));
lmda6
=cos(angle(r,PrismaticSlip6,'antipodal'));
%to check the angles between PrismaticPlaneNormal4 and r, use max(angle(r, PrismaticSlip4,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,



% Now calculate the schmid factor
tau4
=theta4.*lmda4;
tau5
=theta5.*lmda5;
tau6
=theta6.*lmda6;

%combine the variables and calculate the minimum value and the index where
%it occurs. The first 4 columns should be 0 to
%0.5, and the final should be an integer between 1 and 3
PrismaticArray=[tau4,tau5,tau6];
[val,I]=max(PrismaticArray,[],2);
PrismaticArray=[tau4,tau5,tau6,val,I];

%% Pyramidal slip calculations
%define slip plane and directions

% normal to the slip plane is
m7
= Miller(1,1,-2,2,cs,'hkl');  
m8
= Miller(1,-2,1,2,cs,'hkl');
m9
= Miller(-2,1,1,2,cs,'hkl');
m10
= Miller(-1,-1,2,2,cs,'hkl');
m11
= Miller(-1,2,-1,2,cs,'hkl');
m12
= Miller(2,-1,-1,2,cs,'hkl');

% slip direction in the slip plane
n7
= Miller(-1,-1,2,3,cs,'uvw');
n8
= Miller(-1,2,-1,3,cs,'uvw');
n9
= Miller(2,-1,-1,3,cs,'uvw');
n10
= Miller(1,1,-2,3,cs,'uvw');
n11
= Miller(1,-2,1,3,cs,'uvw');
n12
= Miller(-2,1,1,3,cs,'uvw');

% calculate Schmid factors for each point in the map
%angle between force and normal
PyrPlaneNormal7= ebsd('Magnesium').orientations*m7;
PyrPlaneNormal8= ebsd('Magnesium').orientations*m8;
PyrPlaneNormal9= ebsd('Magnesium').orientations*m9;
PyrPlaneNormal10= ebsd('Magnesium').orientations*m10;
PyrPlaneNormal11= ebsd('Magnesium').orientations*m11;
PyrPlaneNormal12= ebsd('Magnesium').orientations*m12;

theta7
=cos(angle(r, PyrPlaneNormal7,'antipodal'));
theta8
=cos(angle(r, PyrPlaneNormal8,'antipodal'));
theta9
=cos(angle(r, PyrPlaneNormal9,'antipodal'));
theta10
=cos(angle(r, PyrPlaneNormal10,'antipodal'));
theta11
=cos(angle(r, PyrPlaneNormal11,'antipodal'));
theta12
=cos(angle(r, PyrPlaneNormal12,'antipodal'));
%to check the angles between PyrPlaneNormal12 and r, use max(angle(r, PyrPlaneNormal12,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,



%calculate the orientation of the 3 slip directions
PyrSlip7=ebsd('Magnesium').orientations*n7;
PyrSlip8=ebsd('Magnesium').orientations*n8;
PyrSlip9=ebsd('Magnesium').orientations*n9;
PyrSlip10=ebsd('Magnesium').orientations*n10;
PyrSlip11=ebsd('Magnesium').orientations*n11;
PyrSlip12=ebsd('Magnesium').orientations*n12;



%calculate the angle between the slip and the force
lmda7
=cos(angle(r,PyrSlip7,'antipodal'));
lmda8
=cos(angle(r,PyrSlip8,'antipodal'));
lmda9
=cos(angle(r,PyrSlip9,'antipodal'));
lmda10
=cos(angle(r,PyrSlip10,'antipodal'));
lmda11
=cos(angle(r,PyrSlip11,'antipodal'));
lmda12
=cos(angle(r,PyrSlip12,'antipodal'));

%to check the angles between PyrSlip7 and r, use max(angle(r, PyrSlip7,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,


% Now calculate the schmid factor
tau7
=theta7.*lmda7;
tau8
=theta8.*lmda8;
tau9
=theta9.*lmda9;
tau10
=theta10.*lmda10;
tau11
=theta11.*lmda11;
tau12
=theta12.*lmda12;

%combine the variables and calculate the minimum value and the index where
%it occurs. The first 6 columns should be 0 to
%0.5, and the final should be an integer between 1 and 6
PyrArray=[tau7,tau8,tau9,tau10,tau11,tau12];
[val,I]=max(PyrArray,[],2);
PyrArray=[tau7,tau8,tau9,tau10,tau11,tau12,val,I];
%% %% Ext Twin calculations
%define slip plane and directions
% normal to the slip plane
m13
= Miller(-1,0,1,2,cs,'hkl');
m14
= Miller(1,0,-1,2,cs,'hkl');
m15
= Miller(1,-1,0,2,cs,'hkl');
m16
= Miller(-1,1,0,2,cs,'hkl');
m17
= Miller(0,1,-1,2,cs,'hkl');
m18
= Miller(0,-1,1,2,cs,'hkl');


% slip direction in the slip plane
n13
= Miller(1,0,-1,1,cs,'uvw');
n14
= Miller(-1,0,1,1,cs,'uvw');
n15
= Miller(-1,1,0,1,cs,'uvw');
n16
= Miller(1,-1,0,1,cs,'uvw');
n17
= Miller(0,-1,1,1,cs,'uvw');
n18
= Miller(0,1,-1,1,cs,'uvw');

% calculate Schmid factors for each point in the map
%angle between force and normal
ExtTwinNormal13= ebsd('Magnesium').orientations*m13;
ExtTwinNormal14= ebsd('Magnesium').orientations*m14;
ExtTwinNormal15= ebsd('Magnesium').orientations*m15;
ExtTwinNormal16= ebsd('Magnesium').orientations*m16;
ExtTwinNormal17= ebsd('Magnesium').orientations*m17;
ExtTwinNormal18= ebsd('Magnesium').orientations*m18;


%note that the 'antipodal' command is not used because twinning is
%directional.  this will give negative values of cos (range 1 to -1) when angle >90 degrees
%calculate the angle between the twin shear plane normal and the force
theta13=cos(angle(r, ExtTwinNormal13));
theta14=cos(angle(r, ExtTwinNormal14));
theta15=cos(angle(r, ExtTwinNormal15));
theta16=cos(angle(r, ExtTwinNormal16));
theta17=cos(angle(r, ExtTwinNormal17));
theta18=cos(angle(r, ExtTwinNormal18));

%to check the angles between ExtTwinNormal16 and r, use max(angle(r, ExtTwinNormal16)/degree)
%also check min, values should range from 0 to +180 degrees, because antipodal is not used




%calculate the orientation of the 6 slip directions
ExtTwinDir13=ebsd('Magnesium').orientations*n13;
ExtTwinDir14=ebsd('Magnesium').orientations*n14;
ExtTwinDir15=ebsd('Magnesium').orientations*n15;
ExtTwinDir16=ebsd('Magnesium').orientations*n16;
ExtTwinDir17=ebsd('Magnesium').orientations*n17;
ExtTwinDir18=ebsd('Magnesium').orientations*n18;


%note that the 'antipodal' command is not used because twinning is
%directional.  this will give negative values of cos (range 1 to -1) when angle >90 degrees
%calculate the angle between the slip and the force
lmda13
=cos(angle(r,ExtTwinDir13));
lmda14
=cos(angle(r,ExtTwinDir14));
lmda15
=cos(angle(r,ExtTwinDir15));
lmda16
=cos(angle(r,ExtTwinDir16));
lmda17
=cos(angle(r,ExtTwinDir17));
lmda18
=cos(angle(r,ExtTwinDir18));

%to check the angles between ExtTwinDir16 and r, use max(angle(r, ExtTwinDir16,'antipodal')/degree)
%also check min, values should range from 0-90 degrees,



% Now calculate the schmid factor
tau13
=theta13.*lmda13;
tau14
=theta14.*lmda14;
tau15
=theta15.*lmda15;
tau16
=theta16.*lmda16;
tau17
=theta17.*lmda17;
tau18
=theta18.*lmda18;

%combine the variables and calculate the minimum value and the index where
%it occurs.
ExtTwinArray=[tau13,tau14,tau15,tau16,tau17,tau18];
[val,I]=max(ExtTwinArray,[],2);
ExtTwinArray=[tau13,tau14,tau15,tau16,tau17,tau18,val,I];

%% Scale outputs, assemble



TotalArray=[tau1/BasalFactor,tau2/BasalFactor,tau3/BasalFactor,tau4/PrismaticFactor,tau5/PrismaticFactor,tau6/PrismaticFactor,tau7/PyramidalFactor,tau8/PyramidalFactor,tau9/PyramidalFactor,tau10/PyramidalFactor,tau11/PyramidalFactor,tau12/PyramidalFactor,tau13/ExtTwinFactor,tau14/ExtTwinFactor,tau15/ExtTwinFactor,tau16/ExtTwinFactor,tau17/ExtTwinFactor,tau18/ExtTwinFactor];
[val,I]=max(TotalArray,[],2);
TotalArray=[tau1/BasalFactor,tau2/BasalFactor,tau3/BasalFactor,tau4/PrismaticFactor,tau5/PrismaticFactor,tau6/PrismaticFactor,tau7/PyramidalFactor,tau8/PyramidalFactor,tau9/PyramidalFactor,tau10/PyramidalFactor,tau11/PyramidalFactor,tau12/PyramidalFactor,tau13/ExtTwinFactor,tau14/ExtTwinFactor,tau15/ExtTwinFactor,tau16/ExtTwinFactor,tau17/ExtTwinFactor,tau18/ExtTwinFactor,val,I];
%% to plot these values

%to see the basal schmid factors
%to see the basal schmid factors
figure;plot(ebsd('Magnesium'),BasalArray(:, 4))
%or scaled
figure;plot(ebsd('Magnesium'),BasalArray(:, 4)/BasalFactor)

%to see the prismatic schmid factors
figure;plot(ebsd('Magnesium'),PrismaticArray(:, 4))
%or scaled
figure;plot(ebsd('Magnesium'),PrismaticArray(:, 4)/PrismaticFactor)

%to see the pyramidal schmid factors
figure;plot(ebsd('Magnesium'),PyrArray(:, 7))
%or scaled
figure;plot(ebsd('Magnesium'),PyrArray(:, 7)/PyramidalFactor)

%to see the extension twin schmid factors-note negative values will not produce
%twinning undeer any level of stress
figure;plot(ebsd('Magnesium'),ExtTwinArray(:, 6))
top=caxis;
top=top(1,2);
setColorRange([0, top]) ;


%or scaled
figure;plot(ebsd('Magnesium'),ExtTwinArray(:, 6)/ExtTwinFactor)
top=caxis;
top=top(1,2)
setColorRange([0, top]) ;

%to see the highest schmid factors of all systems, scaled relative to each
%other
figure
;plot(ebsd('Magnesium'),TotalArray(:, 19))

%to map the active slip or twinning system
figure
;
%with colorcoding-Basal red, prismatic yellow, pyramidal green, ext twin
%Blue
colors
= zeros(size(TotalArray(:, 14)));  
colors
(TotalArray(:, 20) > 0 & TotalArray(:, 20) <= 3) = 1;     % red (1)
colors
(TotalArray(:, 20) > 3 & TotalArray(:, 20) <= 7) = 40;  % yellow (2)
colors
(TotalArray(:, 20) > 7 & TotalArray(:, 20) <= 12) = 65;      % green (3)
colors
(TotalArray(:, 20) > 12 & TotalArray(:, 20) <= 18) = 100;    %blue(4)
plot
(ebsd('Magnesium'),colors)
cmap
= [1 0 0; 1 1 0; 0 1 0; 0 0 1];  
colormap
(cmap);

%plot an IPF in X
figure
;oM = ipdfHSVOrientationMapping(ebsd('Magnesium')); oM.inversePoleFigureDirection =xvector;color = oM.orientation2color(ebsd('Magnesium').orientations); plot(ebsd('Magnesium'),color);
%plot a Pole figure- contour type
h
=[Miller(0,0,1,ebsd('Magnesium').CS), Miller(1,0,0,ebsd('Magnesium').CS)];figure;o=ebsd('Magnesium').orientations; cs= ebsd('Magnesium').CS;  odf = calcODF(ebsd('Magnesium').orientations,'halfwidth',5*degree); plotPDF(odf,h,'upper','projection','eangle','contourf','minmax');set(gcf,'name','ebsd','NumberTitle','off');

%plot a Pole figure- discrete type
h
=[Miller(0,0,1,ebsd('Magnesium').CS), Miller(1,0,0,ebsd('Magnesium').CS)];figure; plotPDF(ebsd('Magnesium').orientations,h,'upper','projection','eangle','points','all', 'MarkerSize',1);set(gcf,'name','ebsd','NumberTitle','off');

save
;



Ralf Hielscher

unread,
Jul 12, 2016, 5:55:08 PM7/12/16
to MTEX
Hi Jessica,

thank you for posting this. It helped me a lot. If you have a look at mtex 4.4.alpha3 you will observe that there is a new class - slipsSystem. This should allow you to write your code in a much more direct manner. Like this

Ralf.

PS: Would you like to check this?

%% General inputs
% Note that to simplify things, the schmid factor output will always be >0
% for slip and this code does not distinguish if the slip is in the
% negative direction.  For twinning, schmid factor output will always be
% -90 to +90, as it is directional. note that for magnesium a plane and the
% normal will not necessarily have the same index values (unlike cubic
% systems), so use the ,'hkl' suffix for plane normals and the 'uvw' suffix
% for directions. %define system crystallography


cs
=ebsd('Titanium').CS;



%define tensile axis or for an arbitrary vector use r = vector3d(1,-1,0);
r
=xvector;


%% define slip systems and twinnings
% define  CRSS of the basal, prismatic, pyramidal, extension twin so that
% the schmid factor can be scaled relatively between the systems based on
% Figure 11, Chapuis, A. & Driver, J. H. Temperature dependency of slip and
% twinning in plane strain compressed magnesium single crystals Acta
% Mater., 2011, 59, 1986-1994 in MPa.  BasalFactor is 4.  Either input
% directly

sSBasal
= slipSystem.basal(cs,4/4);
sSBasal
= sSBasal.symmetrise('antipodal')

sSPrismatic
= slipSystem.prismatic2A(cs,18/4);
sSPrismatic
= sSPrismatic.symmetrise('antipodal')

% type II pyramidal
sSPyramidal
= slipSystem.pyramidal2CA(cs,40/4);
sSPyramidal
= sSPyramidal.symmetrise('antipodal')

%{10-12} extension twinning
sSExtTwin
= slipSystem.tensilTwinT1(cs,11/4);
sSExtTwin
= sSExtTwin.symmetrise('antipodal')


%% Basal slip calculations


rCS
= ebsd('Titanium').orientations \ r;

sFBasal
= abs(sSBasal.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFBasal,[],2))

%% Prismatic slip calculations

sFPrismatic
= abs(sSPrismatic.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFPrismatic,[],2))


%% Pyramidal slip calculations

sFPyramidal
= abs(sSPyramidal.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFPyramidal,[],2))


%% %% Ext Twin calculations


sFExtTwin
= abs(sSPyramidal.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFExtTwin,[],2))


%% Scale outputs, assemble


sFTotal
= [sFBasal ./ sSBasal.CRSS, sFPrismatic ./ sSPrismatic.CRSS,...
  sFPyramidal
./ sSPyramidal.CRSS, sFExtTwin ./ sSExtTwin.CRSS];


[sFmax,sFid] = max(sFTotal,[],2);



%% to plot these values


%to see the highest schmid factors of all systems, scaled relative to each
%
other
figure
;plot(ebsd('Titanium'),sFmax)



%to map the active slip or twinning system
figure
;



plot
(ebsd(sFid > 0 & sFid <= 3),'facecolor','red','DisplayName','basal')
hold on
plot
(ebsd(sFid > 3 & sFid <= 6),'facecolor','yellow','DisplayName','prismatic')
plot
(ebsd(sFid > 6 & sFid <= 12),'facecolor','green','DisplayName','pyramidal')
plot
(ebsd(sFid > 12 & sFid <= 18),'facecolor','green','DisplayName','ext twin')
hold off


grandr...@gmail.com

unread,
Jul 13, 2016, 4:39:45 PM7/13/16
to MTEX

I downloaded 4.4.alpha3, and my original code works but can't get the new code to run at all - I get an error with the first line;

Error using slipSystem.basal
Too many input arguments.

I switched all instances of Titanium to Magnesium to fit my data so that may have caused an issue...
Jessica

Ralf Hielscher

unread,
Jul 16, 2016, 8:48:59 AM7/16/16
to MTEX
Hi Jessica,

you are right. I did some modifications to MTEX such that it work. You can now download MTEX 4.4.alpha.5. With this the following code


%% General inputs
% Note that to simplify things, the schmid factor output will always be >0
% for slip and this code does not distinguish if the slip is in the
% negative direction.  For twinning, schmid factor output will always be
% -90 to +90, as it is directional. note that for magnesium a plane and the
% normal will not necessarily have the same index values (unlike cubic
% systems), so use the ,'hkl' suffix for plane normals and the 'uvw' suffix
% for directions. %
define system crystallography
mtexdata ti
cs
=ebsd('Titanium').CS;



%define tensile axis or for an arbitrary vector use r = vector3d(1,-1,0);
r
=xvector;


%% define slip systems and twinnings
% define  CRSS of the basal, prismatic, pyramidal, extension twin so that
% the schmid factor can be scaled relatively between the systems based on
% Figure 11, Chapuis, A. & Driver, J. H. Temperature dependency of slip and
% twinning in plane strain compressed magnesium single crystals Acta
% Mater., 2011, 59, 1986-1994 in MPa.  BasalFactor is 4.  Either input
% directly

sSBasal
= slipSystem.basal(cs,4/4);
sSBasal
= sSBasal.symmetrise('antipodal')

sSPrismatic
= slipSystem.prismatic2A(cs,18/4);
sSPrismatic
= sSPrismatic.symmetrise('antipodal')

% type II pyramidal
sSPyramidal
= slipSystem.pyramidal2CA(cs,40/4);
sSPyramidal
= sSPyramidal.symmetrise('antipodal')

%{10-12} extension twinning
sSExtTwin
= slipSystem.tensilTwinT1(cs,11/4);
sSExtTwin
= sSExtTwin.symmetrise('antipodal')

% tensor direction in crystal coordinates
rCS
= ebsd('Titanium').orientations \ r;

%% Basal slip calculations

sFBasal
= abs(sSBasal.SchmidFactor(rCS));

plot
(ebsd('Titanium'),max(sFBasal,[],2))

%% Prismatic slip calculations

sFPrismatic
= abs(sSPrismatic.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFPrismatic,[],2))

%% Pyramidal slip calculations

sFPyramidal
= abs(sSPyramidal.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFPyramidal,[],2))

%% Ext Twin
calculations

sFExtTwin
= abs(sSExtTwin.SchmidFactor(rCS));
plot
(ebsd('Titanium'),max(sFExtTwin,[],2))

%% compare SchmidFactors scaled by CRSS

% combine all slip systems
sS
= [sSBasal;sSPrismatic;sSPyramidal;sSExtTwin];

% compute the relative Schmid factors scaled by CRSS
sFRelative
= abs(sS.SchmidFactor(rCS,'relative'));

%%

% compute the maximum relative Schmid factors
[sFmax,sFid] = max(sFRelative,[],2);

% plot the relatibe
figure
(1); plot(ebsd('Titanium'),sFmax)


%to map the active slip or
twinning system
figure
(2);

plot
(ebsd(sFid > 0 & sFid <= 3),'facecolor',[0.5 0.5 1],'DisplayName','basal')
hold on
plot
(ebsd(sFid > 3 & sFid <= 6),'facecolor',[0.5 1 0.5],'DisplayName','prismatic')
plot
(ebsd(sFid > 6 & sFid <= 12),'facecolor',[1 0.5 0.5],'DisplayName','pyramidal')
plot
(ebsd(sFid > 12 & sFid <= 18),'facecolor','yellow','DisplayName','ext twin')
hold off



should work (of course you have to replace the phase name to magnesium). This is active development. Hence, I would be interested in your opinion how to make things as smooth as possible.

Ralf.

grandr...@gmail.com

unread,
Jul 16, 2016, 10:44:21 PM7/16/16
to MTEX
I might be missing something, but the symmetry doesn't look right.  I'm assuming 'hkl' or 'hkil' designates the plane normal.  So for basal it looks good and

sSBasal = slipSystem.basal(cs,4/4);
>> sSBasal = sSBasal.symmetrise('antipodal')

results in

sSBasal = slipSystem (show methods, plot)
 mineral: Magnesium (6/mmm, X||a*, Y||b, Z||c)
 CRSS: 1
 size: 3 x 1
   U   V   T   W | H   K   I   L
   1   1  -2   0   0   0   0   1
   1  -2   1   0   0   0   0   1
  -2   1   1   0   0   0   0   1

However for prismatic, the code

sSPrismatic = slipSystem.prismatic2A(cs,18/4);
sSPrismatic = sSPrismatic.symmetrise('antipodal')
 
gives the following output
sSPrismatic = slipSystem (show methods, plot)
 mineral: Magnesium (6/mmm, X||a*, Y||b, Z||c)
 CRSS: 4.5
 size: 3 x 1
   U   V   T   W | H   K   I   L
   0   1  -1   0   2  -1  -1   0
  -1   1   0   0   1   1  -2   0
  -1   0   1   0   1  -2   1   0

Which I think is backwards.  The slip direction on the basal plane and the slip direction on the prismatic plane should be identical, but other than that it looks correct.

Pyramidal and Extension twinning look good, but keep in mind the definition of what is an 'extension twin' and what is a 'compression twin' will change with the c/a ratio.  You may want to use some sort of other terminology, for example 'extension twins' are the same in both Mg and Ti, but 'compression twins' for these two materials actually have different crystallography.  If you can write this software you obviously have an excellent understanding of this topic, but for an extra reference {Niewczas, M. Lattice correspondence during twinning in hexagonal close-packed crystals Acta Mater., 2010, 58, 5848-5857} is nice.


comparing the Schmid factor values calculated to my earlier code, they are identical for basal and pyramidal.

I think I found an issue with extension twin. The line
sFExtTwin = abs(sSExtTwin.SchmidFactor(rCS));
This doesn't seem correct, you don't want to take the absolute value, because in twinning the reverse system does not activate, so it should not be considered as an option.  Any negative Schmid factor value will not happen.


If you correct these, I'll compare the new results and let you know what I get.

Thank you for doing this work!  Without your software I would have no way to do this sort of calculation and the flexibility is extremely helpful.

Jessica

Ralf Hielscher

unread,
Jul 17, 2016, 1:23:04 AM7/17/16
to mtex...@googlegroups.com
Hi Jessica,

you can also define a slip system directly by

sSPrismatic  = slipSystem(Miller(2,-1,-1,0,cs,'uvtw'), Miller(0,1,-1,1,cs,'hkl'),CRSS);

I just predefined some common slip systems which might not be correct. You can check the predefined slip systems at the end of the file

mtex/geometry/@slipSystem/slipSystem.m

Thank you very much for your help. Since this is new functionality I would like to add this to the MTEX documentation and in particular I think about your example here.

All the best, Ralf.




********************************************************************
Ralf Hielscher                   Tel: +371-531-38556
Fakultät für Mathematik               +371-531-22200 (Sekr.)
Technische Universität Chemnitz  Fax: +371-531-22109
Reichenhainer Str. 39            E-mail: ralf.hi...@mathematik.tu-chemnitz.de
D-09126 Chemnitz                 http://www.tu-chemnitz.de/~rahi
********************************************************************

--
If you want to reduce the number of emails you get through this forum login to https://groups.google.com/forum/?fromgroups=#!forum/mtexmail, click "My membership" and select "Don't send me email updates". You can still get emails on selected topics by staring them.
---
You received this message because you are subscribed to the Google Groups "MTEX" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mtexmail+u...@googlegroups.com.
Visit this group at https://groups.google.com/group/mtexmail.
For more options, visit https://groups.google.com/d/optout.

Ivanka Mitrovic

unread,
Jul 29, 2016, 12:59:32 PM7/29/16
to MTEX

Hi all,

I am trying out new slipSystem options for calcite
(CS = crystalSymmetry('-3m1', [4.99 4.99 17.064], 'X||a', 'Y||b*', 'Z||c', 'mineral', 'Calcite', 'color', 'light blue')).

Considering pre-defined sS is not working out well for me,
I defined my own slip systems, for example:
% twin e <40-41>{-1018}
sStwinE = slipSystem(Miller(4,0,-4,1,CS,'UVTW'),Miller(-1,0,1,8,CS,'HKIL'));
sSTwinE = sStwinE.symmetrise('antipodal')

But I get a strange thing:
sSTwinE = slipSystem (show methods, plot)
 mineral: Calcite (-3m1, X||a, Y||b*, Z||c)
 CRSS: 1
 size: 3 x 1
         u         v         w       | h         k         l         u
   2.66667   1.33333  0.333333        -1         0         1         8
   1.33333  -1.33333 -0.333333        -1         1         0        -8
  -1.33333  -2.66667  0.333333         0         1        -1         8

For example, in case of r-slip, it seemed okay:
sShighR = slipSystem(Miller(-1,2,-1,0,CS,'UVTW'),Miller(1,0,-1,4,CS,'HKIL'));
sSHighR = sShighR.symmetrise('antipodal')
sSHighR = slipSystem (show methods, plot)
 mineral: Calcite (-3m1, X||a, Y||b*, Z||c)
 CRSS: 1
 size: 3 x 1
   u   v   w | h   k   l   u
   0   1   0   1   0  -1   4
  -1  -1   0   1  -1   0  -4
   1   0   0   0  -1   1   4

Any ideas what am I doing wrong?

Also, when calculating "basic" Schmid factor:
[tauMax,mActive,nActive,tau,ind] = calcShearStress(sigmaRot,m,n,'symmetrise') 
depending on how I write m and n, results are very different:
for example, f<2-201> and f<02-21> have drastically different results, 
in first case the max is 0.46, and on second it's 0.2. And this is with symmerise option.
What is the best way to approach this issue?

Thanks a lot!
Ivanka

Ralf Hielscher

unread,
Jul 29, 2016, 5:35:07 PM7/29/16
to MTEX
Hi Ivanka,

this is a bug in MTEX. Do

edit slipSystem/display

and then change in line 70 'trogonal' to 'trigonal'

This should fix your problem.

Thank you for reporting this.

Ralf.

Ivanka Mitrovic

unread,
Aug 1, 2016, 6:27:10 AM8/1/16
to MTEX
Works great! Thanks!

MTEXNewbie

unread,
Aug 31, 2018, 5:33:38 AM8/31/18
to MTEX
Hi Jessica,

Would it be possible to update the above code (Schmid factor for slip and twinning) for MTEX 5.x version that has improved built-in functions? The complexity might reduce a lot, the current script is difficult for beginner users.

Additionally, a revision (as per MTEX 5.x) of your MTEX manual would be very helpful to the community as well. :)

grandr...@gmail.com

unread,
Aug 31, 2018, 10:21:02 AM8/31/18
to MTEX
As far as I'm aware, my original code still works.  While it could be simplified by using the new functions, I feel that it still serves as a good 'first principles' method of computing schmid factor.  Each section spells out explicitly line-by-line what is done, and then the results are all assembled. This makes it a good starting point for users who are familiar with the crystallography but not MTEX.  If you want a shorter version, why not write it yourself?  It's a good exercise to stretch your skills.

For the Manual, if anyone finds any sections that are depreciated (i.e. no longer work) I'll update it.  As far as I'm aware, with version 5.x there are new functions available, but the old ones still work.  If I'm wrong about this, please let me know.

Currently I'm working on a significant upgrade to the app which will be very useful for new users: a scripting utility that generates an mtex script based on the plot actions taken in the GUI, so that's where I want to focus my efforts.

Regards, Jessica.
Reply all
Reply to author
Forward
0 new messages