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

# usage: relem = mfs_reselt(nofelt, elin, disp, class, varargin)
#
# Input  nofelt      Number of elements
#        elem(:)     Structure array with element data
#        disp(:, :)  Displacmenent matrix
#        class       Response class
#
# Output relem(:)    Structure array containing element results
#
# with   relem.id                   Element identifier
#        relem.nofpnt               Number of result points
#        relem.coor(:, 2)           Coordinates of result 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 functions computes stresses and strains of 2-dimensional 
# q4 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
#
# --------------------------------------------------------------------

# Stress points

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

# Loop over elements

  for ne = 1 : nofelt

      x  = elem(ne).coor(:, 1);
      y  = elem(ne).coor(:, 2);
      ue = disp(elem(ne).ects, :);
      E  = mfs_emat2d(elem(ne).mat);

#   Loop over stress points

      for k = 1 : ng
          NE = mfs_NE_q4(r(k), s(k));
          BE = mfs_BE_q4(r(k), s(k), x, y);
          coor(k, :) = NE * [x, y];
          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", coor, "rtype", 3, "results", results);

  end

end
