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

# usage: [elements, loads, constraints, nodesets, rc] = ...
#        mfs_imp_elements(msg, fid, physnam, entities, tdata)
#
# Input  msg          File handle of message file
#        fid          Input file handle
#        physnam      Structure with phyical group names
#        entities     Cell array with entity tables
#        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 4.1
# file format). 
#
# Depending on the translation data, elements, loads or constraints are
# built from these data.
#
# Cell array of entity tables contains the entity tables. The entity
# tables are structure arrays containing the fields
#     table(:).entity      entity tag
#     table(:).numphys     number of physical tags
#     table(:).phystag(:)  physical tags
# The tables are sorted according to ascending entity tags.
#
# ------------------------------------------------------------------------

  rc  = 1;

# Check arguments

  if (nargin != 5 || 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 = struct(); vel = struct(); acc = struct();

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

  loads = 0; constraints = 0;

# Header record

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

# Number of blocks

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

# Process blocks

  rc = 0;

  for nB = 1 : numBlocks

      line  = fgetl(fid);
      words = strsplit(line);
      entityDim = sscanf(words{1}, "%d");
      entityTag = sscanf(words{2}, "%d");
      elemType  = sscanf(words{3}, "%d");

      % Get Elements

      ne = sscanf(words{4}, "%d");
      [ids, nofnod, nodes] = mfs_readelt(fid, ne);

      % Get physical data

      entityTable = entities{entityDim + 1};
      index       = lookup([entityTable.entity], entityTag);
      numphys     = entityTable(index).numphys;
      phystag     = entityTable(index).phystag;

      for np = 1 : numphys

          ixtag = lookup([physnam.ids], phystag(np));
          pname = physnam.names{ixtag};
          group = tdata.(pname);

          switch group.type

          case "elements"

             if (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
                m = find(gmshtyp == elemType);
                if (m)
                   etype = mfstyp{m};
                else
                   etype = 0;
                   fprintf(msg, "*E* Physical group %s: No translation ",
                           pname);
                   fprintf(msg, "defined for Gmsh element type %d\n", elemType);
                   rc  = 1;
                end
             else
                etype = group.name;
             end

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

             for n = 1 : ne 
                 elements(++nofelt).id  = ids(n);
                 elements(nofelt).nodes = nodes(n, :);
                 elements(nofelt).type  = etype;
                 if (isgeom)
                    elements(nofelt).geom = group.geom;
                 end
                 if (ismat)
                    elements(nofelt).mat = group.mat;
                 end
                 eltgid(nofelt, 1) = ids(n);
                 eltgid(nofelt, 2) = ixtag;
             end

          case "loads"

             switch group.name

             case "point"

                if (elemType != 15)
                   fprintf(msg, "*E* Physical group %s: Load name ", pname);
                   fprintf(msg, "=point= can only be applied to nodes\n")
                   rc = 1;
                else
                   for n = 1 : ne 
                       point(++nofpnt).id = nodes(n, 1);
                       point(nofpnt).data = group.data;
                       if (isfield(group, "lc"))
                          point(nofpnt).lc = group.lc;
                       end
                   end
                end

             case "disp"
                [dsp, nofdsp] = mfs_readmotion(nofdsp, dsp, group, nodes);

             case "velo"
                [vel, nofvel] = mfs_readmotion(nofvel, vel, group, nodes);

             case "acce"
                [acc, nofacc] = mfs_readmotion(nofacc, acc, group, nodes);

             end

          case "constraints"

             switch group.name

             case "prescribed"

                pdofs = mfs_packdofs(group.dofs);

                for n = 1 : ne 
                    for k = 1 : nofnod
                        plist(++nofprs, 1) = nodes(n, k);
                        plist(nofprs, 2)   = pdofs;
                    end
                end

             end

         case "nodeset"

            newnodes = reshape(nodes, 1, numel(nodes));
            if (isfield(nodesets, pname))
               nodesets.(pname) = [nodesets.(pname), newnodes];
            else
               nodesets.(pname) = newnodes;
            end
            nodesets.(pname) = unique(nodesets.(pname)); 

         case "ignore"

         end

      end  % Loop over physical properties

  end % Block loop

# Process prescribed displacements

  if (nofdsp)
     disp = mfs_buildmotion(nofdsp, dsp);
  end

# Process prescribed velocities

  if (nofvel)
     velo = mfs_buildmotion(nofvel, vel);
  end

# Process prescribed accelerations

  if (nofacc)
     acce = mfs_buildmotion(nofacc, acc);
  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
  end

# Build the constraints structure

  if (nofprs)
     constraints = struct("prescribed", prs);
  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
