function mfs_exp_tstress3d(fid, elres, t, title)

# usage: mfs_exp_tstress3d(fid, elres, t, title)

# Input  fid         File handle of output file
#        elres       Structure with element results
#        t(:)        Time steps
#        title       View name
#
# The function exports time-dependent stresses of 3-dimensional
# elements in Gmsh MSH file format.
#
# --------------------------------------------------------------------

# Define output formats

  f3e = " %-10.3e %-10.3e %-10.3e";
  fst = cstrcat("%-5d", f3e, f3e, f3e, "\n");

# Title modifications for shells

  suffix = {" (Upper Surface)", " (Lower Surface)"};

# Get number of elements with stresses

  nofelt  = 0;
  isshell = 0;

  for ntype = 1 : elres.noftyp
      if (elres.types(ntype).rtype > 1)
         nofelt += elres.types(ntype).nofelt;
         isshell = isshell || elres.types(ntype).rtype == 4;
      end
  end

  if (isshell)
     pass = 2;
  else
     pass = 1;
  end

# Loop over time steps

  for p = 1 : pass

      for l = 1 : elres.nofcol

#         Header

          fprintf(fid, "$ElementData\n");

#         Tags

          if (isshell)
             fprintf(fid, "1\n""%s""\n", [title, suffix{p}]);
          else
             fprintf(fid, "1\n""%s""\n", title);
          end
          fprintf(fid, "1\n%-10.3e\n", t(l));
          fprintf(fid, "3\n%-3d\n9\n%-5d\n", l - 1, nofelt);

          for ntype = 1 : elres.noftyp

              nelt = elres.types(ntype).nofelt;
              m1   = elres.types(ntype).ixelt1;
              m2   = m1 + nelt - 1;

              switch elres.types(ntype).rtype

                case 2  % rod
                  R = zeros(1, 9);
                  for m = m1 : m2
                      R(1) = elres.elem(m).results.sig(l);
                      fprintf(fid, fst, elres.elem(m).id, R);
                  end

                case 3  % plane stress
                  for m = m1 : m2
                      T = elres.elem(m).results.TE;
                      sig = squeeze(elres.elem(m).results.sig(:, l, :));
                      R = mfs_sig2dto3d(T, sig, elres.elem(m).nofpnt);
                      fprintf(fid, fst, elres.elem(m).id, R);
                  end

                case 4  % shell
                  for m = m1 : m2
                      T = elres.elem(m).results.TE;
                      sig = elres.elem(m).results.sig(:, l, p);
                      R = mfs_sig2dto3d(T, sig, elres.elem(m).nofpnt);
                      fprintf(fid, fst, elres.elem(m).id, R);
                  end

              end

          end   % Loop over element types

#         End

          fprintf(fid, "$EndElementData\n");

      end       % Loop over time steps

  end           % Loop upper/lower surface

end
