function relem = mfs_reselt(nofelt, elem, disp, class, varargin)

# usage: relem = mfs_reselt(nofelt, elem, disp, class, varargin)
#
# Input  nofelt     Number of elements 
#        elem(:)    Structure array with element data
#        disp(:, :) Displacement matrix
#        class      Response class
#
# Output relem(:)   Structure array containing stress results
#
# with   relem.id                   Element identifier
#        relem.nofpnt               Number of result points
#        relem.coor(:, 2)           Coordinates of results points
#        relem.rtype                Result type: 3 = plane stress
#        relem.results.sig(3, :, :) Stresses: x, y, xy
#        relem.results.eps(3, :, :) Strains:  x, y, xy
#
# The function computes the results of q4r elements. The last index
# of sig and eps correponds to result points.
#
# The variable argument list contains additional arguments depending
# on the class:
#
#    statresp:    none
#    freqresp:    varargin{1}      Circular excitation frequency
#
# -------------------------------------------------------------------

  ib = [2, 3, 4, 1];     % Indices of edge end nodes

# Stress points

  ng = 4;
  [r, s, w] = mfs_gauss2d(ng);

# Loop over elements

  for ne = 1 : nofelt

      coor = elem(ne).coor;
      u    = disp(elem(ne).ects, :);
      E    = mfs_emat2d(elem(ne).mat);

      T    = mfs_rot2mid2d(coor);
      ue   = T * u;

      coor = [coor; 0.5 * (coor + coor(ib, :))];
      x    = coor(:, 1); 
      y    = coor(:, 2);

#   Loop over stress points

      for k = 1 : ng
          NE = mfs_NE_q4(r(k), s(k));
          BE = mfs_BE_q8(r(k), s(k), x, y);
          csp(k, :) = NE * elem(ne).coor;
          eps(:, :, k) = BE * ue;
          sig(:, :, k) = E * eps(:, :, k);
      end

      results   = struct("sig", sig, "eps", eps);
      relem(ne) = struct("id", elem(ne).id, "nofpnt", ng,
                          "coor", csp, "rtype", 3, "results", results);

  end

end
