function [nodx, eltx] = mfs_exp_entities(fid, cmp, set2phys, typ2dim,
                                         eset)

# usage: [nodx, eltx] = mfs_exp_entities(fid, cmp, set2phys, typ2dim,
#                                        eset)
#
# Input  fid              File handle of output file
#        cmp              Structure with component
#        set2phys         Structure relating sets to physical names
#        typ2dim(noftyp)  Array with element dimensions
#        eset             Structure with modified element sets
#
# Output nodx             Structure with nodal points to export
#        eltx             Structure with elements to export
#
# This function builds the entities and writes them to the output
# file.
#
# ------------------------------------------------------------------------

# Check arguments

  if (nargin != 5 || nargout != 2)
     print_usage();
  endif

# Initialize

  nodx = struct(); eltx = struct();

  nofnod = cmp.nodes.nofnod;
  nofelt = cmp.elements.nofelt;
  ncoor  = cmp.nodes.ncoor;
  nnset  = numfields(cmp.nset);
  neset  = numfields(cmp.eset);

# Build array with information on sets the elements belong to

  len          = 3 + neset;
  etype        = cmp.elements.index(:, 3);
  eltinf       = zeros(nofelt, len);
  eltinf(:, 1) = 1 : nofelt;                         %  index
  eltinf(:, 2) = typ2dim(etype);                     %  dim.
  eltinf(:, 3) = [cmp.elements.types.gmshid](etype); %  elementType 

  setnames = fieldnames(eset);
  for n = 1 : neset
      eindx = eset.(setnames{n});
      eltinf(eindx, n + 3) = 1;
  endfor

  eltinf = sortrows(eltinf, 2 : len);

# Build element entities

  ndim    = zeros(1, 4);  % Points, Curves, Surfaces, Volumes
  nbl     = 0;            % counts element blocks
  nebl    = 1;            % counts elements in block
  ix1     = 1;            % Position of first element in block

  lentup = len - 1;
  acttup = eltinf(1, 2 : end);
  eltinf = [eltinf; -1(ones(1, len))];  % End marker
  blocks = struct();

  for l = 2 : nofelt + 1

      if (sum(eltinf(l, 2 : len) == acttup) == lentup)

         nebl++;

      else

         entityDim   = acttup(1);
         elementType = acttup(2);
         ixdim       = entityDim + 1;
         ndim(ixdim) += 1;
         m           = ndim(ixdim);

         blocks(++nbl) = struct("entityDim", entityDim,
                                "entityTag", m,
                                "elementType", elementType,
                                "numElementsInBlock", nebl,
                                "elementTag", eltinf(ix1 : l - 1, 1));

         if (len > 3)
            sets = find(acttup(3 : lentup));
            numPhys = length(sets);
            if (numPhys)
                physTags = [set2phys.eset.physTag](sets);
            else
                physTags = 0;
            endif
         else
            numPhys  = 0;
            physTags = 0;
         endif

         if (ixdim == 1)  % Add fake XYZ
            entities{ixdim}(m) = struct("Tag", m,
                                        "XYZ", [0, 0, 0],
                                        "numPhys", numPhys,
                                        "physTags", physTags);                         
         else
            entities{ixdim}(m) = struct("Tag", m,
                                        "numPhys", numPhys,
                                        "physTags", physTags);                         
         endif

         nebl = 1; ix1  = l;
         acttup = eltinf(l, 2 : len);

      endif

  endfor

# Build eltx structure (node indices are replaced later)

  for n = 1 : nbl
     eindx = blocks(n).elementTag;
     blocks(n).elementTag = cmp.elements.index(eindx, 1);
     ixe    = cmp.elements.index(eindx, 2);
     nodel  = [cmp.elements.elem(ixe).nodes];
     nebl   = blocks(n).numElementsInBlock;
     nelnod = length(nodel) / nebl;
     blocks(n).nodeTag = reshape(nodel, nelnod, nebl);
  endfor

  eltx = struct("numEntityBlocks", nbl, "numElements", nofelt,
                "minElementTag", cmp.elements.index(1, 1),
                "maxElementTag", cmp.elements.index(end, 1),
                "blocks", blocks);

# Assign nodal points to entities

  nodinf = zeros(nofnod, 3);
  nodinf(:, 1) = 1 : nofnod;  % Index

  for n = nbl : -1 : 1
      nodix = blocks(n).nodeTag;
      nodinf(nodix, 2) = blocks(n).entityDim;
      nodinf(nodix, 3) = blocks(n).entityTag;
  endfor

# Build entities from nodal point sets
# (Only sets with only one nodal point are considered.)

  setnames = fieldnames(cmp.nset);
  m        = ndim(1);
  coor     = zeros(1, 3);

  for n = 1 : nnset

      nodix = cmp.nset.(setnames{n});

      if (length(nodix) == 1)

         physTag = set2phys.nset(n).physTag;

         if (nodinf(nodix, 3) >= 0)
            coor(1 : ncoor)  = cmp.nodes.coor(nodix, :);
            entities{1}(++m) = struct("Tag", m,
                                      "XYZ", coor,
                                      "numPhys", 1,
                                      "physTags", physTag);
            nodinf(nodix, 2) = 0;
            nodinf(nodix, 3) = -m;
            ndim(1) += 1;
         else
            mm = -nodinf(nodix, 3);
            numPhys  = entities{1}(mm).numPhys + 1;
            physTags = [entities{1}(mm).physTags, PhysTag];
            entities{1}(mm).numPhys  = numPhys;
            entities{1}(mm).physTags = physTags;
         endif

      endif

  endfor

  nodinf(:, 3) = abs(nodinf(:, 3));

# Build nodx structure (node indices are replaced later)

  nodinf = sortrows(nodinf, 2 : 3);

  nbl    = 0;            % counts element blocks
  ix1    = 1;
  nnbl   = 1;
  acttup = nodinf(1, 2 : 3);
  nodinf = [nodinf; [-1, -1, -1]];
  blocks = struct();

  for l = 2 : nofnod + 1

      if (sum(nodinf(l, 2 : 3) == acttup) == 2)

         nnbl++;

      else

         entityDim = acttup(1);
         entityTag = acttup(2);

         blocks(++nbl) = struct("entityDim", entityDim,
                                "entityTag", entityTag,
                                "numNodesInBlock", nnbl,
                                "nodeTag", nodinf(ix1 : l - 1, 1));

         nnbl = 1; ix1  = l;
         acttup = nodinf(l, 2 : 3);

      endif

  endfor

  for n = 1 : nbl
      nodix = blocks(n).nodeTag;
      coor  = cmp.nodes.coor(nodix, :);
      if (ncoor == 2)
         coor(:, 3) = 0;
      endif
      blocks(n).xyz = coor';
  end

  nodx = struct("numEntityBlocks", nbl, "numNodes", nofnod,
                "minNodeTag", cmp.nodes.ids(1),
                "maxNodeTag", cmp.nodes.ids(end),
                "blocks", blocks);

# Add coordinate bounds to entities

  minXYZ = zeros(1, 3); maxXYZ = zeros(1, 3);
  minXYZ(1 : ncoor) = cmp.nodes.mincor;
  maxXYZ(1 : ncoor) = cmp.nodes.maxcor;

  for l = 2 : 4
      for m = 1 : ndim(l)
          entities{l}(m).minXYZ = minXYZ;
          entities{l}(m).maxXYZ = maxXYZ;
      endfor
  endfor

# Replace nodal point indices by identifiers

  for n = 1 : length(nodx.blocks)
      nodix = nodx.blocks(n).nodeTag;
      nodx.blocks(n).nodeTag = cmp.nodes.ids(nodix);
  endfor

  for n = 1 : length(eltx.blocks)
      nodix = eltx.blocks(n).nodeTag;
      if (columns(nodix) > 1)
         eltx.blocks(n).nodeTag = cmp.nodes.ids(nodix);
      else
         eltx.blocks(n).nodeTag = cmp.nodes.ids(nodix)';
      endif
  endfor

# Output entities

  fprintf(fid, "$Entities\n");
  fprintf(fid, "%d %d %d %d\n", ndim);

  for m = 1 : ndim(1)
      enty = entities{1}(m);
      fprintf(fid, "%d %f %f %f %d",
              enty.Tag, enty.XYZ, enty.numPhys);
      if (enty.numPhys)
         fprintf(fid, " %d", enty.physTags);
      endif
      fprintf(fid, "\n");
      
  endfor

  for l = 2 : 4
      for m = 1 : ndim(l)
          enty = entities{l}(m);
          fprintf(fid, "%d %f %f %f %f %f %f %d",
                  enty.Tag, enty.minXYZ, enty.maxXYZ, enty.numPhys);
          if (enty.numPhys)
             fprintf(fid, " %d", enty.physTags);
          endif
          fprintf(fid, " 0\n");
      endfor
  endfor

  fprintf(fid, "$EndEntities\n");

end
