function mfs_exp_mesh_axes(fid, cmp)

# usage: mfs_exp_mesh_axes(fid, cmp)
#
# Input  fid     File handle of output file
#        cmp     Structure with component
#
# This function exports element axes to the Gmsh MSH file.
#
# ------------------------------------------------------------------------

# Check arguments

  if (nargin != 2)
     print_usage();
  end

# Check model subtype

  if (! strcmp(cmp.subtype, "3d"))
     printf("*W* mfs_export: element axes supported only for ");
     printf("3d solids\n");
     return;
  end

# Compute element axes

  ids = []; axes = [];

  for ntype = 1 : cmp.elements.noftyp

      m1 = cmp.elements.types(ntype).ixelt1;
      m2 = m1 + cmp.elements.types(ntype).nofelt - 1;

      switch cmp.elements.types(ntype).axes

      case  2 % beam
         [ids2, axes2] = mfs_axes_b2(cmp.elements.elem(m1 : m2));
         ids  = [ids; ids2];
         axes = [axes; axes2];

      case 3  % triangle
         [ids3, axes3] = mfs_axes_t3(cmp.elements.elem(m1 : m2));
         ids  = [ids; ids3];
         axes = [axes; axes3];

      case 4  % quadrangle
         [ids4, axes4] = mfs_axes_q4(cmp.elements.elem(m1 : m2));
         ids  = [ids; ids4];
         axes = [axes; axes4];

      end

  end     % Loop over element types

# Sort element axes according to ascending element identifiers

  [ids, ix] = sort(ids);
  axes = axes(ix, :);

# Write axes to the output file

  nel = length(ids);

  titles = {"Element x-Axis", "Element y-Axis", "Element z-Axis"};
  comps  = {1 : 3, 4 : 6, 7 : 9};

  for k = 1 : 3

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

      fprintf(fid, "1\n""%s""\n", titles{k});
      fprintf(fid, "1\n0.0000\n");
      fprintf(fid, "3\n0\n3\n%-5d\n", nel);

      c = comps{k};

      for n = 1 : nel
          fprintf(fid, "%-5d %-12.5f %-12.5f %-12.5f\n", ...
                  ids(n), axes(n, c));
      end

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

  end
  
end
