function postWriter3D()

% create GiD.post.msh file (only needed if post-processing only)
nInfo = load('nodeInfo.dat');
nNode = length(nInfo);
mID = fopen('test.post.msh','w');
fprintf(mID,'MESH "solidMesh" dimension 3 ElemType Hexahedra Nnode 8\n');
fprintf(mID,'# color 0 255 0\n');
fprintf(mID,'Coordinates\n');
for k = 1:nNode
    fprintf(mID,'%i %12.6g %12.6g %12.6g\n',nInfo(k,:));
end
fprintf(mID,'End Coordinates\n');
eInfo = load('elementInfo.dat');
nElem = length(eInfo);
fprintf(mID,'Elements\n');
for k = 1:nElem
    fprintf(mID,'%i %i %i %i %i %i %i %i %i\n',eInfo(k,1:9));
end
fprintf(mID,'End Elements\n');

bInfo = load('NodesInfoB.dat');
bNode = length(bInfo);
fprintf(mID,'MESH "beamMesh" dimension 3 ElemType Linear Nnode 2\n');
fprintf(mID,'# color 0 0 255\n');
fprintf(mID,'Coordinates\n');
for k = 1:bNode
    fprintf(mID,'%i %12.6g %12.6g %12.6g\n',bInfo(k,:));
end
fprintf(mID,'End Coordinates\n');
bElem = load('PileInfo.dat');
nBeam = length(bElem);
fprintf(mID,'Elements\n');
for k = 1:nBeam
    fprintf(mID,'%i %i %i\n',bElem(k,:));
end
fprintf(mID,'End Elements\n');
fclose(mID);

% create GiD.post.msh file (contains all of the post-process results)
fID = fopen('d2w1z3t4.post.res','w');
fprintf(fID,'GiD Post Results File 1.0\n');
fprintf(2,'* %s\n','Creating flavia.res from FEA. THIS MAY TAKE A FEW MINUTES ...');

% LOAD AND PROCESS DISPLACEMENT DATA
disp = load('displacement.out');
time = disp(:,1);
disp(:,1) = [];
[nStep,nDisp] = size(disp);
bDisp = load('beamDisp.out');
bDisp(:,1) = [];
for i = 1:nStep
    fprintf(fID,'ResultGroup "Loading Analysis" %12.6g OnNodes\n',time(i));
    fprintf(fID,'ResultDescription "Nodal Displacements" Vector:3\n');
    fprintf(fID,'ComponentNames "X-Displacement" "Y-Displacement" "Z-Displacement"\n');
    fprintf(fID,'Values\n');
    % solid nodes
    for k = 1:nNode
        fprintf(fID,'%i %12.6g %12.6g %12.6g\n',nInfo(k,1),disp(i,(3*k-2):(3*k)));
    end
    % beam nodes
    for k = 1:bNode
        fprintf(fID,'%i %12.6g %12.6g %12.6g\n',bInfo(k,1),bDisp(i,(6*k-5):(6*k-3)));
    end
    fprintf(fID,'End Values\n');
end
fprintf(2,'* %s\n','Done with displacements ...')

% LOAD AND PROCESS SOLID ELEMENT RESULTS
stress = load('stress.out');
stress(:,1) = [];
strain = load('strain.out');
strain(:,1) = [];

fprintf(fID,'GaussPoints "elemGP" ElemType Hexahedra\n');
fprintf(fID,'Number Of Gauss Points: 1\n');
fprintf(fID,'Natural Coordinates: Internal\n');
fprintf(fID,'End GaussPoints\n');
for i = 1:nStep
    fprintf(fID,'ResultGroup "Loading Analysis" %12.6g OnGaussPoints "elemGP"\n',time(i));
    fprintf(fID,'ResultDescription "Gauss Point Stress" Matrix\n');
    fprintf(fID,'ComponentNames "sxx" "syy" "szz" "sxy" "syz" "sxz"\n');
    fprintf(fID,'ResultDescription "Gauss Point Yield Ratio" Scalar\n');
    fprintf(fID,'ResultDescription "Gauss Point Mean Stress" Scalar\n');
    fprintf(fID,'ResultDescription "Gauss Point Norm(Eta)" Scalar\n');
    fprintf(fID,'ResultDescription "Gauss Point Strain" Matrix\n');
    fprintf(fID,'ComponentNames "exx" "eyy" "ezz" "gxy" "gyz" "gzx"\n');
    fprintf(fID,'Values\n');
    for k = 1:nElem
        fprintf(fID,'%i %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g',eInfo(k,1),stress(i,(7*k-6):(7*k)));
        sig = stress(i,(7*k-6):(7*k-1));
        p = sum(sig(1:3))/3;
        normS = ((sig(1)-p)^2 + (sig(2)-p)^2 + (sig(3)-p)^2 + 2*sum(sig(4:6).*sig(4:6)))^0.5;
        fprintf(fID,' %12.6g %12.6g',p,normS);
        fprintf(fID,' %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g\n',strain(i,(6*k-5):(6*k)));
    end
    fprintf(fID,'End Values\n');
end
fprintf(2,'* %s\n','Done with stress and strain ...')

% LOAD AND PROCESS PILE SHEAR FORCE AND BENDING MOMENT
beam = load('globalForces.out');
beam(:,1) = [];
bTally = max(bInfo(1,:));
nPtr = zeros(bTally,1);
for k = 1:bNode
   nPtr(bInfo(k,1)) = k;
end
Px = zeros(bNode,1);
Py = Px; Pz = Px; Mx = Px; My = Px; Mz = Px;
bCount = Px;
for j = 1:nStep
    for k = 1:nBeam
        temp = beam(j,(12*k-11):(12*k));
        I = nPtr(bElem(k,2));
        J = nPtr(bElem(k,3));

        Px(I) = Px(I) - temp(1);
        Py(I) = Py(I) - temp(2);
        Pz(I) = Pz(I) - temp(3);
        Mx(I) = Mx(I) - temp(4);
        My(I) = My(I) - temp(5);
        Mz(I) = Mz(I) - temp(6);
        
        Px(J) = Px(J) + temp(7);
        Py(J) = Py(J) + temp(8);
        Pz(J) = Pz(J) + temp(9);
        Mx(J) = Mx(J) + temp(10);
        My(J) = My(J) + temp(11);
        Mz(J) = Mz(J) + temp(12);

        bCount(I) = bCount(I) + 1;
        bCount(J) = bCount(J) + 1;
    end
    fprintf(fID,'ResultGroup "Loading Analysis" %12.6g OnNodes\n',time(j));
    fprintf(fID,'ResultDescription "Pile Internal Forces" Vector\n');
    fprintf(fID,'ComponentNames "P_x" "P_y" "P_z"\n');
    fprintf(fID,'ResultDescription "Pile Internal Moments" Vector\n');
    fprintf(fID,'ComponentNames "M_x" "M_y" "M_z"\n');
    fprintf(fID,'Values\n');
    for k = 1:bNode
        mWrite = [Px(k)/bCount(k),...
                  Py(k)/bCount(k),...
                  Pz(k)/bCount(k),...
                  Mx(k)/bCount(k),...
                  My(k)/bCount(k),...
                  Mz(k)/bCount(k)];
        fprintf(fID,'%i %12.6g %12.6g %12.6g %12.6g %12.6g %12.6g\n',bInfo(k,1),mWrite);
    end
    fprintf(fID,'End Values\n');
end
fprintf(2,'* %s\n','Done with pile shear and moments ...')

fclose(fID);

return
