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 element results
#
# with
#
#        relem.id                  Element identifier
#        relem.nofpnt              Number of result points
#        relem.coor(2)             Coordinates of result point
#        relem.rtype               Result type: 1 = beam
#        relem.results.rslt(3, :)  N, Qy, Mz
#
# This function computes the results of 2-dimensional b2 elements.
#
# The variable argument list contains additional arguments depending
# on the class:
#
#    statresp:    none
#    freqresp:    varargin{1}      Circular excitation frequency
#
# -------------------------------------------------------------------

# Loop over elements

  for nel = 1 : nofelt

     d = elem(nel).coor(2, :) - elem(nel).coor(1, :);
     L = sqrt(d(1)^2 + d(2)^2);
     c = d(1) / L;
     s = d(2) / L;

     T = eye(6, 6, "double");
     T(1 : 2, 1 : 2) = [c, s; -s, c];
     T(4 : 5, 4 : 5) = [c, s; -s, c];

     uE = T * disp(elem(nel).ects, :);

     coor = mean(elem(nel).coor);

     L6 = 6 * L;

     kN = elem(nel).mat.E * elem(nel).geom.A / L;
     kM = elem(nel).mat.E * elem(nel).geom.I / L;
     kQ = kM / L^2;

     R(1, :) = kN * (uE(4, :) - uE(1, :));                       % N
     R(2, :) = -kQ * ([12, L6, -12, L6] * uE([2; 3; 5; 6], :));  % Qy
     R(3, :) = kM * (uE(6, :) - uE(3, :));                       % Mz

     results    = struct("rslt", R);
     relem(nel) = struct("id", elem(nel).id, "nofpnt", 1, 
                         "coor", coor, "rtype", 1, "results", results);

  end

end
