function stress = mfs_getstress3d(elem, ixelt, cols)

# usage: stress = mfs_getstress3d(elem, ixelt, cols)
#
# Input  elem(:, :)  Structure array with element results
#        ixelt(:)    List of element indices
#        cols(:)     Column indices
#
# Output stress{:}   Cell array with element stresses
#
# The function returns a cell array of structures containing the
# stresses of 3d-elements. The fields of the structure depend on the
# result type.
#
# Column indices correspond to loadcases or excitation frequencies,
# depending on the class of results.
#
# --------------------------------------------------------------------

  nofelt = length(ixelt);
  stress = cell(nofelt, 1);

  for n = 1 : nofelt

      ixe  = ixelt(n);
      res  = elem(ixe);
      rtype = res.rtype;

      if (rtype == 0)
         printf("*E* Requested results not implemented for element %5.0d\n",
                res.id);
         continue;
      end

#     Fields common to all result types

      s = struct("id", res.id, "nofpnt", res.nofpnt,
                 "nofcol", length(cols), "coor", res.coor,
                  "rtype", rtype);

#     Fields depending on result type

      switch rtype

        case 1 % beam

           % stresses computed by mfs_beamstress

        case 2 % rod

           s.sig = res.results.sig(cols);

        case 3 % membrane

           s.TE = res.results.TE;
           sig  = res.results.sig;
           if (res.nofpnt == 1)
              s.sigx  = sig(1, cols);
              s.sigy  = sig(2, cols);
              s.tauxy = sig(3, cols);
           else
              s.sigx  = squeeze(sig(1, cols, :)).';
              s.sigy  = squeeze(sig(2, cols, :)).';
              s.tauxy = squeeze(sig(3, cols, :)).';
           end

        case 4 % shell

           s.TE     = res.results.TE;
           sig      = res.results.sig;
           s.sigxu  = sig(1, cols, 1);
           s.sigyu  = sig(2, cols, 1);
           s.tauxyu = sig(3, cols, 1);
           s.sigxl  = sig(1, cols, 2);
           s.sigyl  = sig(2, cols, 2);
           s.tauxyl = sig(3, cols, 2);

      end

      stress{n} = s;

  end

end
