function data = mfs_merge_read(fid)

# usage: data = mfs_merge_read(fid)
#
# Input  fname  File handle of file to read
#
# Output data   Structure with contents of input file
#
# The function reads a file in Gmsh msh format 4.1 The following
# datablocks are supported:
#
#   PhysicalNames
#   Entities
#   Nodes
#   Elements
#   NodeData
#
# Fields of structure data:
#
#   PhysicalNames.numPhysicalNames
#                .dimension(:)
#                .physicalTag(:)
#                .name{:}
#   Entities.numGeom(4)
#           .Points{:}
#           .Curves{:}
#           .Surfaces{:}
#           .Volumes{:}
#   Nodes.numEntityBlocks
#        .numNodes
#        .minNodeTag
#        .maxNodeTag
#        .blocks(:).entityDim
#                  .entityTag
#                  .numNodesInBlock
#                  .nodeTag(:)
#                  .xyz(:, 3)
#   Elements.numEntityBlocks
#           .numElements
#           .minElementTag
#           .maxElementTag
#           .blocks(:).entityDim
#                     .entityTag
#                     .elementType
#                     .numElementsInBlock
#                     .elementTag(:)
#                     .nodeTag(:, :)
#   NodeData.numNodeData
#           .blocks(:).title
#                     .tstep
#                     .tstepno
#                     .numComp
#                     .numNodes
#                     .values(numComp + 1, numNodes)
#
# --------------------------------------------------------------------

  EntityFields = {"Points", "Curves", "Surfaces", "Volumes"};

# Check arguments

  if (nargin != 1 || nargout != 1)
     print_usage();
  endif

# Initialize

  PhysicalNames = struct("numPhysicalNames", 0);
  Entities      = struct("numGeom", [0, 0, 0, 0]);
  Nodes         = struct("numEntityBlocks", 0);
  Elements      = struct("numEntityBlocks", 0);
  NodeData      = struct("numNodeData", 0); 
  data          = struct("PhysicalNames", PhysicalNames,
                         "Entities"     , Entities,
                         "Nodes"        , Nodes,
                         "Elements"     , Elements,
                         "NodeData"     , NodeData);

# Process the contents

  while ((line = fgetl(fid)) != -1)

     key = sscanf(line, "%s");

     switch key

     case "$PhysicalNames"

        printf("               found datablock $PhyscialNames\n");
        numphys = fscanf(fid, "%d\n", 1);
        for n = 1 : numphys
            [dim(n), tag(n), name{n}] = ...
            fscanf(fid, "%d %d %s\n", "C");
        endfor
        data.PhysicalNames = struct("numPhysicalNames", numphys,
                                    "dimension", dim,
                                    "physicalTag", tag,
                                    "name", {name});

        line = fgetl(fid);

     case "$Entities"

        printf("               found datablock $Entities\n");
        numGeom     = fscanf(fid, "%d %d %d %d\n", 4);
        data.Entities = struct("numGeom", numGeom);
        for n = 1 : 4
            if (numGeom(n))
               enty = cell(1, numGeom(n));
               for m = 1 : numGeom(n)
                   line  = fgetl(fid);
                   enty{m} = sscanf(line, "%f");
               endfor
               data.Entities.(EntityFields{n}) = enty;
            endif
        endfor

        line = fgetl(fid);

     case "$Nodes"

        printf("               found datablock $Nodes\n");

        [numBlocks, numNodes, minNodeTag, maxNodeTag] = ...
           fscanf(fid, "%d %d %d %d\n", "C");

        blocks = struct();
        nn = 0;
        for n = 1 : numBlocks
            [entityDim, entityTag, dummy, numNodesInBlock] = ...
               fscanf(fid, "%d %d %d %d\n", "C");
            if (numNodesInBlock)
               nodeTag = fscanf(fid, "%d", numNodesInBlock);
               coor    = fscanf(fid, "%f\n", [3, numNodesInBlock]);
               blocks(++nn) = struct("entityDim", entityDim,
                                     "entityTag", entityTag,
                                     "numNodesInBlock", numNodesInBlock,
                                     "nodeTag", nodeTag,
                                     "xyz", coor);
            endif
        endfor
        numBlocks = nn;

        data.Nodes = struct("numEntityBlocks", numBlocks,
                            "numNodes"       , numNodes,
                            "minNodeTag"     , minNodeTag,
                            "maxNodeTag"     , maxNodeTag,
                            "blocks"         , blocks);

        fskipl(fid);

     case "$Elements"

        printf("               found datablock $Elements\n");

        [numBlocks, numElements, minElementTag, maxElementTag] = ...
           fscanf(fid, "%d %d %d %d\n", "C");

        blocks = struct();
        for n = 1 : numBlocks
            [entityDim, entityTag, elemType, numElementsInBlock] = ...
              fscanf(fid, "%d %d %d %d\n", "C");
            [elementTag, ~, nodeTag] = ...
              mfs_readelt(fid, numElementsInBlock);
            blocks(n) = struct("entityDim"         , entityDim,
                               "entityTag"         , entityTag,
                               "elementType"       , elemType,
                               "numElementsInBlock", numElementsInBlock,
                               "elementTag"        , elementTag,
                               "nodeTag"           , nodeTag');
        endfor

        data.Elements = struct("numEntityBlocks", numBlocks,
                               "numElements"    , numElements,
                               "minElementTag"  , minElementTag,
                               "maxElementTag"  , maxElementTag,
                               "blocks"         , blocks);

        fskipl(fid);

     case "$NodeData"

        printf("               found datablock $NodeData\n");
        n = ++data.NodeData.numNodeData;
        fskipl(fid);
        title = fgetl(fid);
        fskipl(fid, 1);
        tstep = fscanf(fid, "%f\n", 1);
        fskipl(fid, 1);
        tstepno  = fscanf(fid,"%d\n", 1);
        numComp  = fscanf(fid,"%d\n", 1);
        numNodes = fscanf(fid,"%d\n", 1);
        values   = fscanf(fid, "%f\n", [numComp + 1, numNodes]);
        data.NodeData.blocks(n) = struct("title"   , title,
                                         "tstep"   , tstep,
                                         "tstepno" , tstepno,
                                         "numComp" , numComp,
                                         "numNodes", numNodes,
                                         "values"  , values);

        fskipl(fid);

     otherwise

        mfs_skip(fid);

     endswitch

  endwhile

endfunction
