function [elements, loads, constraints, nodesets, elsets, rc] = ...
         mfs_imp_elements(msg, fid, physnam, tdata)

# usage: [elements, loads, constraints, nodesets, rc] = ...
#        mfs_imp_elements(msg, fid, physnam, tdata)
#
# Input  msg          File handle of message file
#        fid          Input file handle
#        physnam      Structure with phyical group names
#        tdata        Structure with translation data
#
# Output elements     Structure array with element data
#        loads        Structure with load data
#                     (0 if there are no loads)
#        constraints  Structure with contraint data
#                     (0 if there are no constraints)
#        nodesets     Structure with node sets
#        elsets       Structure with element sets
#        rc           Return code: 0 = no errors, 1 = errors found
#
# The function reads the elements from from the input file (version 2.2
# file format). 
#
# Depending on the translation data, elements, loads or constraints are
# built from these data.
#
# ------------------------------------------------------------------------

  rc  = 1;

# Check arguments

  if (nargin != 4 || nargout != 6)
     print_usage();
  end

# Initialize

  nofelt = 0;  % Number of elements
  nofpnt = 0;  % Number of point loads
  nofdsp = 0;  % Number of prescribed displacements
  nofvel = 0;  % Number of prescribed velocities
  nofacc = 0;  % Number of prescribed accelerations
  nofprs = 0;  % Number of constraints of type prescribed

  dsp.k  = 0;  % Counts data entries for prescribed displacements
  vel.k  = 0;  % Counts data entries for prescribed velocities
  acc.k  = 0;  % Counts data entries for prescribed accelerations

  nodesets = struct();
  elsets   = struct();

# Header record

  line = fgetl(fid);
  if (! strcmp(line, "$Elements"))
     fprintf(msg, "*E* $Elements expected but not found\n");
     return;
  end

# Number of lines in $Element block

  line   = fgetl(fid);
  noflin = sscanf(line, "%d", 1);

# Element data

  rc = 0;
  l  = 1;

  [id, etype, groupid, nodes] = mfs_readelt(fid);

  while (l <= noflin)

      current_group = groupid;
      group_index   = lookup(physnam.ids, groupid);
      current_name  = physnam.names{group_index};
      group         = getfield(tdata, current_name); 

      switch group.type

      case "elements"

         if ((elist = iscell(group.name)))
            na = length(group.name);
            for n = 1 : na
                words = strsplit(group.name{n}, "=");
                gmshtyp(n) = sscanf(words{1}, "%d");
                mfstyp{n}  = strtrim(words{2});
            end
            [gmshtyp, I] = sort(gmshtyp);
            mfstyp = mfstyp(I);
         end

         isgeom = isfield(group, "geom");
         ismat  = isfield(group, "mat");

         while (groupid == current_group)
            elements(++nofelt).id  = id;
            elements(nofelt).nodes = nodes;
            eltgid(nofelt, 1) = id;
            eltgid(nofelt, 2) = group_index;
            if (elist)
               m = lookup(gmshtyp, etype, "m");
               if (m)
                  elements(nofelt).type = mfstyp{m};
               else
                  fprintf(msg, "*E* Physical group %s: No translation ", ...
                          current_name);
                  fprintf(msg, "defined for Gmsh element type %d\n", etype);
                  rc  = 1;
               end
            else
               elements(nofelt).type  = group.name;
            end
            if (isgeom)
               elements(nofelt).geom  = group.geom;
            end
            if (ismat)
               elements(nofelt).mat = group.mat;
            end
            if (++l <= noflin)
               [id, etype, groupid, nodes] = mfs_readelt(fid);
            else
               groupid = 0;
            end
         end

      case "loads"

         switch group.name

         case "point"

            while (groupid == current_group)
                if (etype != 15)
                   fprintf(msg, "*E* Physical group %s: Load name ", ...
                           current_name);
                   fprintf(msg, "=point= can only be applied to nodes\n")
                   rc = 1;
                else
                   point(++nofpnt).id = nodes(1);
                   point(nofpnt).data = group.data;
                   if (isfield(group, "lc"));
                      point(nofpnt).lc = group.lc;
                   end
                end
                if (++l <= noflin)
                   [id, etype, groupid, nodes] = mfs_readelt(fid);
                else
                   groupid = 0;
                end
            end

         case "disp"

            while (groupid == current_group)
                dsp.data{++dsp.k}  = group.data;
                for k = 1 : length(nodes)
                   dsp.ids(++nofdsp, 1) = nodes(k);
                   if (isfield(group, "lc"));
                      lc = group.lc;
                   else
                      lc = 1;
                   end
                   dsp.ids(nofdsp, 2) = lc;
                   dsp.ids(nofdsp, 3) = dsp.k;
                end
                if (++l <= noflin)
                   [id, etype, groupid, nodes] = mfs_readelt(fid);
                else
                   groupid = 0;
                end
            end

         case "velo"

            while (groupid == current_group)
                vel.data{++vel.k}  = group.data;
                for k = 1 : length(nodes)
                   vel.ids(++nofvel, 1) = nodes(k);
                   if (isfield(group, "lc"));
                      lc = group.lc;
                   else
                      lc = 1;
                   end
                   vel.ids(nofvel, 2) = lc;
                   vel.ids(nofvel, 3) = vel.k;
                end
                if (++l <= noflin)
                   [id, etype, groupid, nodes] = mfs_readelt(fid);
                else
                   groupid = 0;
                end
            end

         case "acce"

            while (groupid == current_group)
                acc.data{++acc.k}  = group.data;
                for k = 1 : length(nodes)
                   acc.ids(++nofacc, 1) = nodes(k);
                   if (isfield(group, "lc"));
                      lc = group.lc;
                   else
                      lc = 1;
                   end
                   acc.ids(nofacc, 2) = lc;
                   acc.ids(nofacc, 3) = acc.k;
                end
                if (++l <= noflin)
                   [id, etype, groupid, nodes] = mfs_readelt(fid);
                else
                   groupid = 0;
                end
            end

         end

      case "constraints"

         switch group.name

         case "prescribed"

            while (groupid == current_group)
                pdofs = mfs_packdofs(group.dofs);
                for k = 1 : length(nodes)
                   plist(++nofprs, 1) = nodes(k);
                   plist(nofprs, 2)   = pdofs;
                end
                if (++l <= noflin)
                   [id, etype, groupid, nodes] = mfs_readelt(fid);
                else
                   groupid = 0;
                end
            end

         end

      case "nodeset"

         while (groupid == current_group)
             if (etype != 15)
                fprintf(msg, "*E* Physical group %s: Nodeset ", ...
                        current_name);
                fprintf(msg, "does not contain points\n")
                rc = 1;
             else
                if (! isfield(nodesets, current_name))
                   nodesets.(current_name) = nodes(1);
                else
                   nodesets.(current_name) = ...
                   [nodesets.(current_name), nodes(1)];
                end
             end
             if (++l <= noflin)
                [id, etype, groupid, nodes] = mfs_readelt(fid);
             else
                groupid = 0;
             end
         end

      case "ignore"

         while (groupid == current_group)
             if (++l <= noflin)
                [id, etype, groupid, nodes] = mfs_readelt(fid);
             else
                groupid = 0;
             end
         end

      end

  end

# Process prescribed displacements

  if (nofdsp)
     dsp.ids = sortrows(dsp.ids);
     nd = 0;
     id = dsp.ids(1, 1);
     lc = dsp.ids(1, 2);
     dt = dsp.data{dsp.ids(1, 3)};
     for k = 2 : nofdsp
        if (dsp.ids(k, 1) != id || dsp.ids(k, 2) != lc)
           disp(++nd).id = id;
           disp(nd).data = dt;
           disp(nd).lc   = lc;
           id = dsp.ids(k, 1);
           lc = dsp.ids(k, 2);
           dt = dsp.data{dsp.ids(k, 3)};
        end
     end
     disp(++nd).id = id;
     disp(nd).data = dt;
     disp(nd).lc   = lc;
  end

# Process prescribed velocities

  if (nofvel)
     vel.ids = sortrows(vel.ids);
     nd = 0;
     id = vel.ids(1, 1);
     lc = vel.ids(1, 2);
     dt = vel.data{vel.ids(1, 3)};
     for k = 2 : nofvel
        if (vel.ids(k, 1) != id || vel.ids(k, 2) != lc)
           velo(++nd).id = id;
           velo(nd).data = dt;
           velo(nd).lc   = lc;
           id = vel.ids(k, 1);
           lc = vel.ids(k, 2);
           dt = vel.data{vel.ids(k, 3)};
        end
     end
     velo(++nd).id = id;
     velo(nd).data = dt;
     velo(nd).lc   = lc;
  end

# Process prescribed accelerations

  if (nofacc)
     acc.ids = sortrows(acc.ids);
     nd = 0;
     id = acc.ids(1, 1);
     lc = acc.ids(1, 2);
     dt = acc.data{acc.ids(1, 3)};
     for k = 2 : nofacc
        if (acc.ids(k, 1) != id || acc.ids(k, 2) != lc)
           acce(++nd).id = id;
           acce(nd).data = dt;
           acce(nd).lc   = lc;
           id = acc.ids(k, 1);
           lc = acc.ids(k, 2);
           dt = acc.data{acc.ids(k, 3)};
        end
     end
     acce(++nd).id = id;
     acce(nd).data = dt;
     acce(nd).lc   = lc;
  end

# Process constraints

  if (nofprs)
     plist = sortrows(plist);
     np    = 0;
     id    = plist(1, 1);
     pdofs = plist(1, 2);
     for k = 2 : nofprs
         if (plist(k, 1) == id)
            pdofs = bitor(pdofs, plist(k, 2));
         else
            prs(++np).id = id;
            prs(np).dofs = mfs_unpackdofs(pdofs);
            id    = plist(k, 1);
            pdofs = plist(k, 2);
         end
     end
     prs(++np).id = id;
     prs(np).dofs = mfs_unpackdofs(pdofs);
  end

# End record

  line = fgetl(fid);
  if (! strcmp(line, "$EndElements"))
     fprintf(msg, "*E* $EndElements expected but not found\n");
     rc = 1;
     return;
  end

# Build the loads structure

  if (nofpnt || nofdsp || nofvel || nofacc)
     loads = struct();
     if (nofpnt)
        loads = setfield(loads, "point", point);
     end
     if (nofdsp)
        loads = setfield(loads, "disp", disp);
     end
     if (nofvel)
        loads = setfield(loads, "velo", velo);
     end
     if (nofacc)
        loads = setfield(loads, "acce", acce);
     end
  else
     loads = 0;
  end

# Build the constraints structure

  if (nofprs)
     constraints = struct();
     constraints = setfield(constraints, "prescribed", prs);
  else
     constraints = 0;
  end

# Build the element sets

  [s, ix] = sort(eltgid(:, 2)); 
  eltgids = eltgid(ix, :);
  eltgids(nofelt + 1, :) = [-1, -1];  % End marker

  n = 1;
  while (n < nofelt)
      n1  = n;
      gid = eltgids(n1, 2);
      while (gid == eltgids(++n, 2)) end
      n2 = n - 1;
      setnam = physnam.names{gid};
      elsets = setfield(elsets, setnam, eltgids(n1 : n2, 1)');
  end

end
