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(:, :)  Displacmenent matrix
#        class       Response class
#
# Output relem(:)   Structure containing element results
#
# with   relem.id                    Element identifier
#        relem.nofpnt                Number of result points
#        relem.coor(3)               Coordinates of result point
#        relem.rtype                 Result type: 4 = shell results
#        relem.results.T33(3, 3)     Transformation matrix to element
#                                    coordinate system
#        relem.results.sig(3, :, 2)  Stresses: x, y, xy
#        relem.results.eps(3, :, 2)  Strains:  x, y, xy
#        relem.results.rslt(8, :)    Resultants: Nx, Ny, Nxy,
#                                    Mx, My, Mxy, Qx, Qy
#
# The functions computes stresses, strains and stress resultants of
# s3 elements.
#
# The variable argument list contains additional arguments depending
# on the class:
#
#    statresp:    none
#    freqresp:    varargin{1}      Circular excitation frequency
#
# --------------------------------------------------------------------

  mdofs = [1, 2, 7, 8, 13, 14];             % Membrane
  bdofs = [4, 5, 10, 11, 16, 17];           % Bending
  sdofs = [3, 4, 5, 9, 10, 11, 15, 16, 17]; % Shear
  kS    = 5 / 6;  % Shear correction factor

# Initialize

  ng = 1;
  r  = 1 / 3; s = r;
  NE = [r, r, r];

# Loop over elements

  for ne = 1 : nofelt

      geom = elem(ne).geom;
      h    = geom.t; h2 = h^2 / 6;

      [coore, TE, T33] = mfs_esys(elem(ne).coor, geom);
      ue               = TE * disp(elem(ne).ects, :);
      E                = mfs_emat2d(elem(ne).mat);
      Ghs              = kS * h * E(3, 3);

      x = coore(:, 1); y = coore(:, 2);

      um = ue(mdofs, :); ub = ue(bdofs, :); us = ue(sdofs, :);

      [BM, BB, BS0, BS1, JE] = mfs_BE(x, y);

      Ghs = mfs_Gscale(Ghs, JE, elem(ne).geom);

#     Membrane 

      coor = NE * elem(ne).coor;
      eps(:, :, 1) = BM * um;
      sig(:, :, 1) = E * eps(:, :, 1);
      eps(:, :, 2) = eps(:, :, 1);
      sig(:, :, 2) = sig(:, :, 1);
      rslt(1 : 3, :) = h * sig(:, :, 1);

#     Bending

      epsb = 0.5 * h * BB * ub;
      sigb = E * epsb;
      eps(:, :, 1) += epsb;
      eps(:, :, 2 ) -= epsb;
      sig(:, :, 1) += sigb;
      sig(:, :, 2) -= sigb;
      rslt(4 : 6, :) = h2 * sigb;

#     Transverse shear force

      BS = 2 * (BS0 + [s; r] .* BS1);
      gamma = BS * us;
      rslt(7 : 8, :) = Ghs * gamma;

      results   = struct("TE", T33, "sig", sig, "eps", eps, "rslt", rslt);
      relem(ne) = struct("id", elem(ne).id, "nofpnt", ng,
                         "coor", coor, "rtype", 4, "results", results);

  end

end
