利用MATLAB进行abaqus结果文件后处理：IPF map及曲面截面‌

# ---------------------------------------------------------------------
# The following scripts are to plot the inverse pole figure contour for
#     the Euler angles obtained from CP-FEM simulation result
# ---------------------------------------------------------------------
# *Step-1
# Extract the Euler angle from odb file of CP-FEM, this script should be run in the abaqus
jobName = 'E:\MyPapers\TiSpherInden\CPFEMSimu\Job-1'
odb = session.openOdb(name=jobName + '.odb')
EAFile = open('EAofGP.csv', 'w')
EA_phi1 = odb.steps['Step-1'].frames[-1].fieldOutputs['SDV130'].values
EA_phi = odb.steps['Step-1'].frames[-1].fieldOutputs['SDV131'].values
EA_phi2 = odb.steps['Step-1'].frames[-1].fieldOutputs['SDV132'].values
for i in range(len(EA_phi1) + 1):
    if i == 0:
        EAFile.write('%2s, %2s, %5s, %5s, %5s\n'
                  % ('Element', 'GPoint', 'phi1', 'PHI', 'phi2'))
    else:
        EAFile.write('%6i, %6i, %8.4F, %8.4F, %8.4F\n'
                  %(EA_phi1[i - 1].elementLabel, EA_phi1[i - 1].integrationPoint,
                   EA_phi1[i - 1].data, EA_phi[i - 1].data, EA_phi2[i - 1].data))
else:
    EAFile.close()
# Extract the coordinates of integration Pt from odb file of CP-FEM,
# this script should be run in the abaqus
# It should be noted that abaqus does not
# output the coordinate of int Pt by default and there is NO GUI access
# to set the output of coordinate of int Pt. To output the coordinate of int Pt,
# .inp file should be modified:
# <*Element Output, directions=YES
# LE, PE, PEEQ, PEMAG, S, SDV,COORD>
# i.e. add "COORD" under "*Element Output"
COORDintPtFile = open('COORDofGP.csv', 'w')
COORDIntPt = odb.steps['Step-1'].frames[-1].fieldOutputs['COORD'].values
for i in range(len(COORDIntPt) + 1):
    if i == 0:
        COORDintPtFile.write('%2s, %2s, %5s, %5s, %5s\n'
                        % ('Element', 'GPoint', 'coordX', 'coordY', 'coordZ'))
    else:
        COORDintPtFile.write(
            '%6i, %6i, %8.4F, %8.4F, %8.4F\n' %
            (COORDIntPt[i - 1].elementLabel, COORDIntPt[i - 1].integrationPoint,
             COORDIntPt[i - 1].data[0], COORDIntPt[i - 1].data[1],
             COORDIntPt[i - 1].data[2]))
else:
    COORDintPtFile.close()

%% 根据晶体塑性模拟所得欧拉角画IPF map
% -----------------------------------------------------------------------
% -- the coordinates and euler angles of each gauss point store in the
%        COORDofGP.csv and EAofGP.csv
%        file COORDofGP.csv : Element  GPoint  coordX  coordY  coordZ
%        file EAofGP.csv : Element  GPoint  phi1  PHI  phi2
% -- COORDofGP.csv and EAofGP.csv are extract from the .odb file generated by
%        crystal plastic FE simulation by a python script storing at
%         E:\MyPapers\MyPythonFuncs\EulerAnglePostProc.py
% -- this script needs mtex
% -----------------------------------------------------------------------

clear all
% --- read the coordinates and euler angle of gauss points from the file :
%       COORDofGP.csv and EAofGP.csv
COORDofGP = readtable("E:\MyPapers\TiSpherInden\CPFEMSimu\COORDofGP.csv");
EAofGP = readtable("E:\MyPapers\TiSpherInden\CPFEMSimu\EAofGP.csv");

% --- transform the readed table to matrix
EAofGP = table2array(EAofGP);
COORDofGP = table2array(COORDofGP);

% --- define the crystalSymmetry of pure Titanium
cs = crystalSymmetry('6/mmm',[3,3,4.761],'x||a','mineral','Titanium (Alpha)');

% --- define the ipfKey, the TSL key is the color used by common EBSD tests
%     reference : https://mtex-toolbox.github.io/EBSDAdvancedMaps.html
ipfKey = ipfTSLKey(cs);
plot(ipfKey)
figure
for i = 1 : length(EAofGP)
    phi1 = EAofGP(i,3);
    Phi = EAofGP(i,4);
    phi2 = EAofGP(i,5);
    % -- define the orientation by Euler angle
    ori = orientation.byEuler(phi1*degree,Phi*degree,phi2*degree,cs);

    % -- calculate the rbg color according to orientation using ipfKey
    rgb = ipfKey.orientation2color(ori);

    % -- plot IPF with colors
    plotIPDF(ori,vector3d.Z,'MarkerFaceColor',rgb,'MarkerSize',20)
    hold on
end

parfor i = 1 : length(EAofGP)
    phi1 = EAofGP(i,3);
    Phi = EAofGP(i,4);
    phi2 = EAofGP(i,5);
    ori = orientation.byEuler(phi1*degree,Phi*degree,phi2*degree,cs);
    rgb = ipfKey.orientation2color(ori)
    colors(i,:) = rgb;
end

% --- 在 COORDofGP 的最大与最小值范围内生成一个空间点阵
x = linspace(min(COORDofGP(:,3)),max(COORDofGP(:,3)),6);
y = linspace(min(COORDofGP(:,4)),max(COORDofGP(:,4)),6);
z = linspace(min(COORDofGP(:,5)),max(COORDofGP(:,5)),6);

% --- 去掉 x, y, z 的最大与最小值，因为在边界处对 color 插值会得到 NaN
x(1) = [];x(end) = [];
y(1) = [];y(end) = [];
z(1) = [];z(end) = [];

[X, Y, Z] = meshgrid(x,y,z);

% --- 通过插值获得 X, Y, Z 散点处的 colors 值
Cred =  ...
  griddata(COORDofGP(:,3),COORDofGP(:,4),COORDofGP(:,5),colors(:,3),X,Y,Z);
Cblue  ...
= griddata(COORDofGP(:,3),COORDofGP(:,4),COORDofGP(:,5),colors(:,4),X,Y,Z);
Cgreen =  ...
  griddata(COORDofGP(:,3),COORDofGP(:,4),COORDofGP(:,5),colors(:,5),X,Y,Z);

% --- section plane
xsurf_x = linspace(x(2),x(end-1),4);
ysurf_y = linspace(y(2),y(end-1),4);
[xsurf,ysurf] = meshgrid(xsurf_x, ysurf_y);

% --- zsurf = xsurf 表明定义的 section plane is Z = X 平面，这里也可以定义
%     section plane 为曲面.
%     可以参考web(fullfile(docroot, 'matlab/ref/slice.html'))中的马鞍面
zsurf = xsurf;


surfaceObjRed = slice(X,Y,Z,Cred,xsurf,ysurf,zsurf);
% --- surfaceObjRed 是 slice 生成的一个 surface 对象。由于 slice 命令中每个
%     网格面片额颜色是由四个顶点处的数值决定的，不能直接设置每个网格面片的
%     rbg颜色，所以这里先用 slice 以不同颜色画出截面，然后提取截面上的顶点
%     坐标和颜色。

Xdata = surfaceObjRed.XData;
Ydata = surfaceObjRed.YData;
Zdata = surfaceObjRed.ZData;

Cred2 = surfaceObjRed.CData;

surfaceObjblue = slice(X,Y,Z,Cblue,xsurf,ysurf,zsurf);
Cblue2 = surfaceObjblue.CData;
surfaceObjgreen = slice(X,Y,Z,Cgreen,xsurf,ysurf,zsurf);
Cgreen2 = surfaceObjgreen.CData;

% --- 最后用 patch 逐个画出面片，每个面片都可以根据刚刚所得的顶点颜色值指定
%     一个颜色。(当网格面片尺寸较小，四个顶点处的颜色值差异不大，这里采用
%     左上角 1 号顶点的 rgb 作为该处面片的颜色)
for i = 1 : length(Xdata(:,1))-1
    for j = 1 : length(Xdata(:,1))-1
        v = [Xdata(j,i)     Ydata(j,i)     Zdata(j,i);...
             Xdata(j+1,i)   Ydata(j+1,i)   Zdata(j+1,i); ...
             Xdata(j+1,i+1) Ydata(j+1,i+1) Zdata(j+1,i+1); ...
             Xdata(j,i+1)   Ydata(j,i+1)   Zdata(j,i+1)];
        f = [1 2 3 4];
        figure(2)
        patch('Faces',f,'Vertices',v,'FaceColor', ...
                                [Cred2(j,i),Cblue2(j,i),Cgreen2(j,i)])
    end
end
view(3)