function sig = mfs_beamstress3d(fid, cmp, cdata, ixelt, coor, ldc)

# usage: sig = mfs_beamstress3d(fid, cmp, cdata, ixelt, coor, ldc)
#
# Input  fid             File handle of output file
#        cmp             Structure with component data
#        cdata           Structure with data of class processed
#        ixelt(:)        Indices of elements in elem
#        coor(:, 2)      Coordinates of stress points
#        ldc(:)          Loadcase list
#
# Output sig(:, :, :)    Bending stresses (optional):
#                         1st index: elements
#                         2nd index: stresspoints
#                         3rd index: loadcases
#
# The function computes the normal stresses of 3-dim. beams.
#
# ---------------------------------------------------------------------

  if (nargin != 6 || nargout != 1)
     print_usage();
  endif

# Initialize

  nofelt = length(ixelt);
  nofldc = length(ldc);
  np     = rows(coor);
  sig    = zeros(nofelt, np, nofldc);
  fmt    = "  %10.3e  %10.3e  %10.3e\n";
  elem   = cdata.elres.elem;

# Compute results

  for l = 1 : nofldc

      lc = ldc(l);
      fprintf(fid, "Beam stresses of loadcase %d\n", lc);

      fprintf(fid, "\n  element  point      y           z");
      fprintf(fid, "         sig\n");
      fprintf(fid, "  -----------------------------------");
      fprintf(fid, "---------------\n");

      for k = 1 : nofelt

          ix = ixelt(k);

          if (elem(ix).rtype == 1)

             fprintf(fid, "   %6d  ", elem(ix).id);

             geom = cmp.elements.elem(ix).geom;

             NQM  = elem(ix).results.rslt(:, lc);
             sigN = NQM(1) / geom.A;
             My   = NQM(5);
             Mz   = NQM(6);

             if (geom.Iyz)
                D   = 1 - geom.Iyz^2 / (geom.Iy * geom.Iz);
                Mys = (My - Mz * geom.Iyz / geom.Iz) / D;
                Mzs = (Mz - My * geom.Iyz / geom.Iy) / D;
                My  = Mys;
                Mz  = Mzs;
             endif

             MyIy = My / geom.Iy;
             MzIz = Mz / geom.Iz;

             for n = 1 : np
                 sigp = sigN + MyIy * coor(n, 2)- MzIz * coor(n, 1);
                 if (n > 1)
                    fprintf(fid, "           ");
                 endif
                 fprintf(fid, "  %2.0d", n);
                 fprintf(fid, fmt, coor(n, :), sigp);
                 sig(k, n, l) = sigp;
             endfor

          else
             fprintf(fid, "is not a beam element\n");
          endif
              
      endfor

      fprintf(fid, "\n");

  endfor

endfunction
