dream2abahex.m

function dream2abahex(voxFileName, inpFileName)
    % *dream2abahex* generates ABAQUS input file 
    % containing C3D8 (brick) mesh and a text file with 
    % grain orientations of a microstructure synthesized 
    % in Dream.3D by reading Los Alamos FFT file
    % 
    % % Syntax
    % dream2abahex(voxFileName, inpFileName)
    %
    % % Input
    % voxFileName - full path to Los Alamos FFT file written by Dream.3D
    % inpFileName - desired path to Abaqus input file 
    % % Example:
    % dream2abahex('dp_64x64x64.vox', 'dp_64x64x64-aba.inp')
    % 
    % % Result
    % ABAQUS input file containing C3D8 mesh with 
    % - element sets with individual phases
    % - sections with individual phases
    % - element sets with grains 
    % - node sets of faces for easier assignment of BCs
    % and orientations (tex) file containing
    % Euler angle sets for each element in the mesh
    %
    % Read more at http://latmarat.net/blog/scripts/dream2abahex/
   
    % --------------------------
    % written by
    % Marat I. Latypov while at POSTECH
    % latmarat@postech.edu
    % September 2014
    % --------------------------
    
    tic  
    %% Digest data from Los Alamos FFT file generated in Dream.3D

    % open FFT file
    fid = fopen(voxFileName,'rt');
    rawData = textscan(fid, '%f %f %f %f %f %f %d %d','delimiter',' ');
    fclose(fid);

    % load euler angles, coordinates, grain and phase IDs
    euler   = cell2mat(rawData(1:3));
    xyz     = cell2mat(rawData(4:6));
    grains  = cell2mat(rawData(7));
    phases  = cell2mat(rawData(8));
    [~,order] = sortrows(xyz,[1,3,2]);
    
    % reorder grain IDs, phase IDs, and Euler angles according to abaqus convention
    grains = grains(order);
    phases = phases(order);
    euler = euler(order,:);
     
    % get the number of vox along x, y, z
    xVox = size(unique(xyz(:,1)),1);
    yVox = size(unique(xyz(:,2)),1);
    zVox = size(unique(xyz(:,3)),1);

    % get step size and boundaries for the mesh
    step   = zeros(1,3);
    boxmin = zeros(1,3);
    boxmax = zeros(1,3);
    for ii = 1:3
        step(ii) = min(diff(unique(xyz(:,ii))));
        boxmin(ii) = xyz(1,ii)-step(ii)/2;
        boxmax(ii) = xyz(end,ii)+step(ii)/2;
    end

    %% Generate 3D mesh
    % generate nodes 
    [x,y,z] = meshgrid(boxmin(1):step(1):boxmax(1),boxmin(2):step(2):boxmax(2),boxmin(3):step(3):boxmax(3));
    numNodes = numel(x);
    coord = [reshape(x,numNodes,1), reshape(y,numNodes,1), reshape(z,numNodes,1)];
    nodes = [(1:numNodes)', sortrows(coord,[1,3,2])];

    % allocate array for elements
    elem = zeros(size(xyz,1),9);
    count = 1;

    % start loop over voxel dimensions
    for ix = 1:xVox
        for iz = 1:zVox
            for iy = 1:yVox

                % get element label
                elem(count,1) = count;

                % nodes on the plane with lower x
                elem(count,2) = iy + (iz-1)*(yVox+1) + (ix-1)*(yVox+1)*(zVox+1);
                elem(count,3) = elem(count,2) + 1;
                elem(count,4) = elem(count,3) + yVox + 1;
                elem(count,5) = elem(count,2) + yVox + 1;

                % nodes on the plane with higher x
                elem(count,6) = iy + (iz-1)*(yVox+1) + ix*(yVox+1)*(zVox+1);
                elem(count,7) = elem(count,6) + 1;
                elem(count,8) = elem(count,7) + yVox + 1;
                elem(count,9) = elem(count,6) + yVox + 1;

                count = count+1;
            end
        end
    end

    %% Write inp file
    % open inp file and write keywords 
    inpFile = fopen(inpFileName,'wt');
    fprintf(inpFile,'** Generated by: dream2abahex.m\n');
    fprintf(inpFile,'**PARTS\n**\n');
    fprintf(inpFile,'*Part, name=DREAM\n');

    % write nodes
    fprintf(inpFile,'*NODE, NSET=AllNodes\n');
    fprintf(inpFile,'%d,\t%e,\t%e, \t%e\n',nodes');

    % write elements
    fprintf(inpFile,'*Element, type=C3D8, ELSET=AllElements\n');
    fprintf(inpFile,'%d,\t%d,\t%d,\t%d,\t%d,\t%d,\t%d,\t%d,\t%d\n',elem');

    % create element sets containing grains
    for ii = 1:numel(unique(grains))
        fprintf(inpFile,'\n*Elset, elset=Grain-%d\n',ii);
        fprintf(inpFile,'%d, %d, %d, %d, %d, %d, %d, %d, %d\n',elem(grains==ii)');
    end

	% create element sets containing phases
    uniPhases = unique(phases);
    for ii = 1:numel(unique(phases))
        fprintf(inpFile,'\n*Elset, elset=Phase-%d\n',ii);
        fprintf(inpFile,'%d, %d, %d, %d, %d, %d, %d, %d, %d\n',elem(phases==uniPhases(ii))');
    end
    
    % create sections for phases
    for ii = 1:numel(uniPhases)
        fprintf(inpFile,'\n**Section: Section_Phase-%d\n*Solid Section, elset=Phase-%d, material=Phase-%d\n',ii,ii,ii);
        fprintf(inpFile,'%d, %d, %d, %d, %d, %d, %d, %d, %d\n',elem(phases==uniPhases(ii))');
    end
    fprintf(inpFile,'\n');

    % create node sets containing surface nodes for BCs
    for ii = 1:3
        fprintf(inpFile,'\n**\n*Nset, nset=NODES-%d\n',ii);
        fprintf(inpFile,'%d, %d, %d, %d, %d, %d, %d, %d, %d\n',nodes(nodes(:,ii+1)==boxmin(ii))');
        fprintf(inpFile,'\n**\n*Nset, nset=NODES+%d\n',ii);
        fprintf(inpFile,'%d, %d, %d, %d, %d, %d, %d, %d, %d\n',nodes(nodes(:,ii+1)==boxmax(ii))');
    end

    % write a closing keyword
    fprintf(inpFile,'\n**\n*End Part\n');

    % close the file
    fclose(inpFile);  
    
    
    %% Write aeuler file    
    [inpFileDir,fileName,fileExt] = fileparts(inpFileName);
    if strcmp(inpFileDir,'') == 0
        eulFileName = [inpFileDir '\' fileName '.tex'];
    else
        eulFileName = [fileName '.tex'];
    end 
    eulFile = fopen(eulFileName,'wt');
    
    fprintf(eulFile,'# Euler angle sets (phi1, Phi, phi2) for C3D8 mesh in %s\n', [fileName fileExt]);
    fprintf(eulFile,'# Euler angle sets are sorted according to element labels\n');
    fprintf(eulFile,'# 1 element in mesh - 1 set of Euler angles. Generated by dream2abahex.m\n');
    fprintf(eulFile,'B\t%d\n', size(euler,1));
    fprintf(eulFile,'%13.4f%13.4f%13.4f\n', euler');
    fclose(inpFile);
    
    fprintf('***dream2abahex.m completed***\n');
    fprintf('ABAQUS mesh is written to %s\n', inpFileName);
    fprintf('Orientations are written to %s\n', eulFileName);
    
    toc
end