function [ixdofs, elements, eltids, rc] = ...
         mfs_prepelts(elsel, types, indx, elem, eset)

# usage: [ixdofs, elements, eltids, rc] = ...
#        mfs_prepelts(eltids, types, indx, elem, eset)
#
# Input  elsel        List of element identifiers or name of an element set
#        types(:)     Structure array with element type data
#        indx(:, 3)   Element index
#        elem(:)      Structure array with element data
#        eset         Structure with element sets
#
# Output ixdofs(:)    Indices of degrees of freedom referenced by elements
#                     in list
#        elements     Structure with data of elements in list
#        eltids(:)    External element identifiers of selected elements
#        rc           Return code : 0 = no errors
#                                   1 = errors
#
# The function prepares the data needed to compute element results for
# selected elements.
#
# -------------------------------------------------------------------------

# Initialize

  ixdofs   = [];
  elements = struct(); 
  eltids   = [];
  rc       = 0;

# Get element indices

  if (ischar(elsel))

  % from element set

     if (isfield(eset, elsel))
        ixe    = eset.(elsel);
        eltids = indx(ixe, 1);
        ixe    = sort(ixe);
     else
        rc = 1;
        printf("*E* Element set %s does not exist\n", elsel);
     end

  elseif (iscell(elsel))

  % from cell aray of element sets

     n = length(elsel); ixe = []; 
     for n = 1 : n 
        if(isfield(eset, elsel{n}))
           ixe = [ixe, eset.(elsel{n})];
        else
           rc = 1;
           printf("*E* element set %s does not exist\n", elsel{n});
        end
     end 
     eltids = indx(ixe, 1);
     ixe    = sort(ixe);

  else

  % from list of element identifiers

     eltids = elsel;
     ixe    = lookup(indx(:, 1), eltids, "m");
     ixb    = find(ixe == 0);

     if (! isempty(ixb))
        rc = 1;
        for k = 1 : length(ixb)
            badid = eltids(ixb(k));
            printf("*E* Element %6.0 does not exist\n", badid);
        end
     end

  end

  if (rc) return; end

# Get degrees of freedom referenced by elements

  dofs   = [];
  nofelt = length(ixe);

  for n = 1 : nofelt
      dofs = [dofs, elem(indx(ixe(n), 2)).ects];
  end

  ixdofs = unique(dofs);

# Build structure elements containing only elements in list

  elements.nofelt = nofelt;

  % Preliminary structure types

  ixtypes = unique(indx(ixe, 3));
  elements.noftyp = length(ixtypes);
  elements.types  = types(ixtypes);

  % Structure index

  elements.index(:, 1) = indx(ixe, 1);
  elements.index(:, 3) = lookup(ixtypes, indx(ixe, 3));
  [ixe, ixpos] = sort(ixe);
  elements.index(:, 2) = ixpos;

  % Adjust types.nofelt and types.ixelt1

  for n = 1 : elements.noftyp
      ix = find(elements.index(:, 3) == n);
      elements.types(n).nofelt = length(ix);
      elements.types(n).ixelt1 = min(ix);
  end

  % Structure elem

  elements.elem = elem(indx(ixe, 2));
  for n = 1 : nofelt
      elements.elem(n).ects = lookup(ixdofs, elements.elem(n).ects, "m");
  end

end
