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 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: 1 = beam results
#        relem.results.rslt(6, :) N, Qy, Qz, Mx, My, Mz
#
# This function computes the results of 3-dimensional b2 elements.
#
# The variable argument list contains additional arguments depending
# on the class:
#
#    statresp:    none
#    freqresp:    varargin{1}      Circular excitation frequency
#    transresp:   varargin{1}      Time steps
#
# -------------------------------------------------------------------

  for nel = 1 : nofelt

# Geometry

     if (isfield(elem(nel).geom, "Q"))
        v = elem(nel).geom.Q - elem(nel).coor(1, :);
     else
        v = elem(nel).geom.v;
     end
 
     d = elem(nel).coor(2, :) - elem(nel).coor(1, :);
     L = norm(d);

     ex = d' / L;
     ez = v' - (v * ex) * ex;
     ez = ez / norm(ez);
     ey = cross(ez, ex);

# Transformation matrix from structural to element system

     T33 = [ex'; ey'; ez'];
     T   = zeros(12, 12, "double");
     T( 1 :  3,  1 : 3) = T33;
     T( 4 :  6,  4 : 6) = T33;
     T( 7 :  9,  7 : 9) = T33;
     T(10 : 12, 10 :12) = T33;

     if (isfield(elem(nel).geom, "CE"))
        TE  = eye(12, 12);
        TE(2,  4) = -elem(nel).geom.CE(2);
        TE(3,  4) =  elem(nel).geom.CE(1);
        TE(8, 10) =  TE(2, 4);
        TE(9, 10) =  TE(3, 4);
        T = TE * T;
     end

     if (isfield(elem(nel).geom, "P"))
        TA = eye(12, 12);
        TA(1, 5) = -elem(nel).geom.P(2);
        TA(1, 6) =  elem(nel).geom.P(1);
        TA(2, 4) = -TA(1, 5);
        TA(3, 4) = -TA(1, 6);
        TA(7, 11) = -elem(nel).geom.P(2);
        TA(7, 12) =  elem(nel).geom.P(1);
        TA(8, 10) = -TA(7, 11);
        TA(9, 10) = -TA(7, 12);
        T = TA * T;
     end

# Displacements in element system

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

# Material data

     E = elem(nel).mat.E;
     G = E / (2 * (1 + elem(nel).mat.ny));

# Normal force N => R(1)

     R(1, :) = E * elem(nel).geom.A * (uE(7, :) - uE(1, :)) / L;

# Shear forces Qy and Qz => R(2 : 3)

     d3v  = 12 * (uE(2, :) - uE(8, :)) / L^3 + 6 * (uE(6, :) + uE(12, :)) / L^2;
     d3w  = 12 * (uE(3, :) - uE(9, :)) / L^3 - 6 * (uE(5, :) + uE(11, :)) / L^2;

     R(2, :) =  E * (elem(nel).geom.Iyz * d3w - elem(nel).geom.Iz * d3v);
     R(3, :) = -E * (elem(nel).geom.Iy * d3w - elem(nel).geom.Iyz * d3v);

# Torsional moment Mx => R(4)

     R(4, :) = G * elem(nel).geom.IT * (uE(10, :) - uE(4, :)) / L;

# Bending moments My and Mz => R(5 : 6)

     EL    = E / L;
     Edphi = EL * (uE(11, :) - uE(5, :));
     Edpsi = EL * (uE(12, :) - uE(6, :));

     R(5, :) = elem(nel).geom.Iy * Edphi + elem(nel).geom.Iyz * Edpsi;
     R(6, :) = elem(nel).geom.Iyz * Edphi + elem(nel).geom.Iz * Edpsi;

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

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

  end

end
