function mfs_mergex(fnin1, fnin2, fnout)

# usage: mfs_mergex(fnin1, fnin2, fnout)
#
# Input  fnin1     Name of first input file
#        fnin2     Name of second input file
#        fnout     Name of output file
#
# The function merges two msh-files into one. The input files must both be
# Gmsh version 4.1 files. The output file will also be a Gmsh version 4.1
# file.
#
# The following datablocks are supported:
#
#   PhysicalNames
#   Entities
#   Nodes
#   Elements
#   NodeData
#
# -----------------------------------------------------------------------------

  numPhysTagPos = [5, 8, 8, 8]; % position of number of physical tags
  EntityFields  = {"Points", "Curves", "Surfaces", "Volumes"};

# Check arguments

  if (nargin != 3 || nargout != 0)
     print_usage();
  endif

# Read the first input file

  [fid, vers, rc] = mfs_mshfileopen(fnin1);
  if (rc) return; endif

  if (! strcmp(vers, "4.1"))
     printf("*E* mfs_merge: file %s is not version 4.1\n", fnin1);
     printf("               using format \"msh22\" might work\n");
     return
  endif

  printf("    mfs_merge: reading file %s\n", fnin1);
  data1 = mfs_merge_read(fid);

  fclose(fid);

# Read the second input file

  [fid, vers, rc] = mfs_mshfileopen(fnin2);
  if (rc) return; endif

  if (! strcmp(vers, "4.1"))
     printf("*E* mfs_merge: file %s is not version 4.1\n", fnin2);
     printf("               using format \"msh22\" might work\n");
     return
  endif

  printf("    mfs_merge: reading file %s\n", fnin2);
  data2 = mfs_merge_read(fid);

  fclose(fid);

# Merge Physical Names

  data = struct();

  numPhys1 = data1.PhysicalNames.numPhysicalNames;
  numPhys2 = data2.PhysicalNames.numPhysicalNames;
  pinc     = numPhys1;

  if (numPhys1)
     if (numPhys2)
        numPhys = numPhys1 + numPhys2;
        dim     = [data1.PhysicalNames.dimension, ...
                   data2.PhysicalNames.dimension];
        Tag     = [data1.PhysicalNames.physicalTag, ...
                   data2.PhysicalNames.physicalTag + pinc];
        name    = [data1.PhysicalNames.name, ...
                   data2.PhysicalNames.name];
        [~, ix] = sort(dim);
        Tag     = Tag(ix);
        name    = name(ix);
        data.PhysicalNames = struct("numPhysicalNames", numPhys,
                                    "dimension"       , dim,
                                    "physicalTag"     , Tag,
                                    "name"            , {name});
     else
        data.PhysicalNames = data1.PhysicalNames;
     endif
  else
     data.PhysicalNames = data2.PhysicalNames;
  endif

# Merge Entities

  numGeom1 = data1.Entities.numGeom;
  numGeom2 = data2.Entities.numGeom;
  numGeom  = numGeom1 + numGeom2;

  data.Entities = struct("numGeom", numGeom);

  for m = 1 : 4

      k = numPhysTagPos(m);

      if (numGeom1(m))

         if (numGeom2(m))

            Geom = data2.Entities.(EntityFields{m});
            for n = 1 : numGeom2(m)
                Geom{n}(1) += numGeom1(m);
                numPhys = Geom{n}(k);
                if (numPhys)
                   Geom{n}(k + 1 : k + numPhys) += pinc;
                endif
            endfor

            data.Entities.(EntityFields{m}) = ...
              [data1.Entities.(EntityFields{m}), Geom];

         else
            data.Entities.(EntityFields{m}) = ...
              data1.Entities.(EntityFields{m});
         endif

      elseif (numGeom2(m))

         Geom = data2.Entities.(EntityFields{m});
         for n = 1 : numGeom2(m)
             numPhys = Geom{n}(k);
             if (numPhys)
                Geom{n}(k + 1 : k + numPhys) += pinc;
             endif
         endfor

         data.Entities.(EntityFields{m}) = Geom;

      endif

  end

# Merge Nodes

  numBlk1 = data1.Nodes.numEntityBlocks;
  numBlk2 = data2.Nodes.numEntityBlocks;

  if (numBlk1 && numBlk2)

     ninc       = data1.Nodes.maxNodeTag;

     numBlk     = numBlk1 + numBlk2;
     numNodes   = data1.Nodes.numNodes + data2.Nodes.numNodes;
     minNodeTag = data1.Nodes.minNodeTag;
     maxNodeTag = data2.Nodes.maxNodeTag + ninc;

     blk2 = data2.Nodes.blocks;
     for m = 1 : numBlk2
         k = blk2(m).entityDim + 1;
         blk2(m).entityTag += numGeom1(k);
         blk2(m).nodeTag   += ninc;
     endfor

     blk = [data1.Nodes.blocks, blk2];

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

  elseif (numBlk1)

     printf("*E* mfsmerge: $Nodes only in file %s\n", fnin1); 
     return;

  elseif (numBlk2)

     printf("*E* mfsmerge: $Nodes only in file %s\n", fnin2); 
     return;

  endif

# Merge Elements

  numBlk1 = data1.Elements.numEntityBlocks;
  numBlk2 = data2.Elements.numEntityBlocks;

  if (numBlk1 && numBlk2)

     einc          = data1.Elements.maxElementTag;

     numBlk        = numBlk1 + numBlk2;
     numElements   = data1.Elements.numElements + ...
                     data2.Elements.numElements;
     minElementTag = data1.Elements.minElementTag;
     maxElementTag = data2.Elements.maxElementTag + einc;

     blk2 = data2.Elements.blocks;
     for m = 1 : numBlk2
         k = blk2(m).entityDim + 1;
         blk2(m).entityTag  += numGeom1(k);
         blk2(m).elementTag += einc;
         blk2(m).nodeTag    += ninc;
     endfor

     blk = [data1.Elements.blocks, blk2];

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

  elseif (numBlk1)

     printf("*E* mfs_merge: $Elements only in file %s\n", fnin1); 
     return;

  elseif (numBlk2)

     printf("*E* mfs_merge: $Elements only in file %s\n", fnin2); 
     return;

  endif

# Merge NodeData

  numData1 = data1.NodeData.numNodeData;
  numData2 = data2.NodeData.numNodeData;

  if (numData1 && numData2)

     if (numData1 != numData2)
        printf("*E* mfs_merge: inconsistent number of $NodeData blocks\n");
        printf("               file %s contains %d blocks\n", ...
               fnin1, numData1);
        printf("               file %s contains %d blocks\n", ...
               fnin2, numData2);
        return;
     endif

     numData = numData1;
     nodinc  = max(data1.NodeData.blocks(1).values(1, :), [], 2);

     for n = 1 : numData
         blck1 = data1.NodeData.blocks(n);
         blck2 = data2.NodeData.blocks(n);
         if ((numComp = blck1.numComp) != blck2.numComp)
            printf("*E* mfs_merge: inconsistent number components found");
            printf(" (%d / %2)\n", blck1.numComp, blck2.numComp);
            printf("               Data block titles = \n");
            printf("               %s\n", blck1.title);
            printf("               %s\n", blck2.title);
            return
         endif
         if (blck1.tstepno != blck2.tstepno)
            printf("*E* mfs_merge: inconsistent time step numbers found");
            printf(" (%d / %2)\n", blck1.tstepno, blck2.tstepno);
            printf("               Data block titles = \n");
            printf("               %s\n", blck1.title);
            printf("               %s\n", blck2.title);
            return
         endif
         blck2.values(1, :) += nodinc;
         numNodes = blck1.numNodes + blck2.numNodes;
         blocks(n) = struct("title"   , blck1.title,
                            "tstep"   , blck1.tstep,
                            "tstepno" , blck1.tstepno,
                            "numComp" , blck1.numComp,
                            "numNodes", numNodes,
                            "values"  , [blck1.values, blck2.values]);
     endfor

     data.NodeData = struct("numNodeData", numData,
                            "blocks"     , blocks);

  elseif (numData1)

     printf("*E* mfs_merge: $NodeData only in file %s\n", fnin1); 
     return;

  elseif (numData2)

     printf("*E* mfs_merge: $NodeData only in file %s\n", fnin2); 
     return;

  endif

# Write the output file

  if ((fout = fopen(fnout, "wt")) == -1)
     printf("*E* mfs_merge: could not open file %s\n", fnout);
     return; 
  end

  printf("    mfs_merge: writing file %s\n", fnout);

  fprintf(fout, "$MeshFormat\n4.1 0 %d8\n$EndMeshFormat\n", sizeof(1));

  % Physical Names

  if ((numPhys = data.PhysicalNames.numPhysicalNames))
     fprintf(fout, "$PhysicalNames\n");
     fprintf(fout, "%d\n", numPhys);
     dim  = data.PhysicalNames.dimension;
     Tag  = data.PhysicalNames.physicalTag;
     name = data.PhysicalNames.name;
     for n = 1 : numPhys
         fprintf(fout, "%d %d %s\n", dim(n), Tag(n), name{n});
     end
     fprintf(fout, "$EndPhysicalNames\n");
  end

  % Entities

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

  if ((ngeom = numGeom(1)))
     Geom = data.Entities.(EntityFields{1});
     k    = numPhysTagPos(1);
     for n = 1 : ngeom
        fprintf(fid, "%d %f %f %f %d", Geom{n}(1 : k));
        if ((numPhys = Geom{n}(k)))
           fprintf(fid, " %d", Geom{n}(k + 1 : k + numPhys));
        endif
        fprintf(fid, "\n");
     endfor
  endif

  for m = 2 : 4
     if ((ngeom = numGeom(m)))
        Geom = data.Entities.(EntityFields{m});
        k    = numPhysTagPos(m);
        for n = 1 : ngeom
           fprintf(fid, "%d %f %f %f %f %f %f %d", Geom{n}(1 : k));
           if ((numPhys = Geom{n}(k)))
              fprintf(fid, " %d", Geom{n}(k + 1 : k + numPhys));
           endif
           fprintf(fid, " 0\n");
        endfor
     endif
  endfor

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

  % Nodes

  if (isfield(data, "Nodes"))
     mfs_exp_nodes(fid, data.Nodes);
  endif

  % Elements

  if (isfield(data, "Elements"))
     mfs_exp_elts(fid, data.Elements);
  endif

  % NodeData

  if (isfield(data, "NodeData"))
     mfs_exp_nodedata(fid, data.NodeData);
  endif

  fclose(fout);

end
