function [elements, mxdofpnt, rc] = ...
mfs_new_elements(msg, input_elem, nodes, subtype, mxdofpnt)

# usage: [elements, mxdofpnt, rc] = ...
#        mfs_new_elements(msg, input_elem, nodes, subtype, mxdofpnt)
#
# Input  msg         File handle of message file 
#        input_elem  Structure array with input data of elements
#        nodes       Structure with nodal point data
#        subtype     Subtype
#        mxdofpnt    Maximum number of degrees of freedom at a nodal point
#
# Output elements    Structure with element data
#        mxdofpnt    Maximum number of degrees of freedom at a nodal point
#        rc          Return code: 0 means no erros
#                                 1 means errors
#
# The function builds a structure with all element data from the
# input data.
#
# Input data:
#
#   input_elem(:).id         Element identifier
#   input_elem(:).type       Element type
#   input_elem(:).nodes(:)   Identifiers of element nodal points
#   input_elem(:).mat        Structure with element material data
#   input_elem(:).geom       Structure with element geometry data
#
#   nodes                    see mfs_new_nodes
#
# Output data:
#
#   elements.nofelt          Number of elements
#   elements.noftyp          Number of element types
#   elements.types(:)        Element type data
#   elements.index(:, 3)     Element index: identifier, index, type
#   elements.elem(:)         Element data
#
#   with
#
#   types.name               Type name
#   types.nelnod             Number of element nodal points
#   types.ndofnod            Number of degrees of freedom per nodal point
#   types.mat                indicates if material data are needed
#   types.geom               indicates if geometrical data are needed
#   types.gmshid             Geometry type in Gmsh
#   types.rtype              Result type
#   types.axes               Axes type
#   types.nofelt             Number of elements of this type
#   types.ixelt1             Index of first element of this type
#                            in elem(:)
#
#   elem.id                  Element identifier
#   elem.nodes(:)            Identifiers of element nodal points
#   elem.coor(:, ndim)       Coordinates of element nodal points
#   elem.ects(:)             Element connectivity table
#   elem.mat                 Structure with element material data
#   elem.geom                Structure with element geometry data 
#
# The element connectivity table contains the indices of the global degrees 
# of freedom the element connects to.
#
# -----------------------------------------------------------------------------

  elem_fields = {"id", "type", "nodes"};

# Check arguments

  if (nargin != 5 || nargout != 3)
     print_usage();
  end

# Initialize

  rc       = 0;
  elements = 0;

# Build path to element routines

  fname = [mfilename(), ".m"];
  epath = file_in_loadpath(fname);
  len   = length(epath) - length(fname);
  epath = strcat(epath(1 : len), "elem/", subtype, "/");

# Check input data

  if ((rc = mfs_checkfields1(msg, input_elem, elem_fields, "Elements")))
     return;
  end 

# Get element identifiers

  nofelt = length(input_elem);
  eltids = [input_elem.id];

  if (length(eltids) < nofelt)
     rc = 1;
     for n = 1 : nofelt
         if (isempty(input_elem(n).id))
            fprintf(msg, "*E* Identifier of element number %5d undefined\n", n);
         end
     end
     return;
 end

# Get element types and sort them

  element_types = cell(1, nofelt);
  for n = 1 : nofelt
      if (ischar(input_elem(n).type))
         element_types{n} = input_elem(n).type;
      else
         rc = 1;
         fprintf(msg, "*E* Element %6d: type not defined\n", eltids(n));
      end
  end

  if (rc) return; end

  [s, idxtype] = sort(element_types); 

# Get information on element types

  current_type = s{1};
  noftyp       = 1;
  ixelt1       = 1;
  nofelt_type  = 1;

  for ix = 2 : nofelt
      acttype = input_elem(idxtype(ix)).type;
      if (strcmp(current_type, acttype))
         nofelt_type++;
      else
         [typedata, rct] = mfs_eldata(msg, current_type, epath);
         if (rct)
            rc = 1;
         else
            typd = setfield(typedata, "nofelt", nofelt_type);
            typd = setfield(typd, "ixelt1", ixelt1);
            types(noftyp) = typd;
            mxdofpnt = max(mxdofpnt, typedata.ndofnod);
         end
         current_type = acttype;
         noftyp++;
         ixelt1 = ix;
         nofelt_type = 1;
      end
      ix++;
  end

  [typedata, rct] = mfs_eldata(msg, current_type, epath);
  if (rct)
     rc = 1;
  else
     typd = setfield(typedata, "nofelt", nofelt_type);
     typd = setfield(typd, "ixelt1", ixelt1);
     types(noftyp) = typd;
     mxdofpnt = max(mxdofpnt, typedata.ndofnod);
  end

  if (rc) return; end

# Process element data

  nelt  = 1;

  for m = 1 : noftyp

      nelnod      = types(m).nelnod;
      nofelt_type = types(m).nofelt;
      ixelt1      = types(m).ixelt1;
      needgeom    = types(m).geom;
      needmat     = types(m).mat;

      if (needgeom && ! isfield(input_elem, "geom"))
         fprintf(msg, "*E* Elements: Field \"geom\" undefined\n");
         rc = 1;
      end

      if (needmat && ! isfield(input_elem, "mat"))
         fprintf(msg, "*E* Elements: Field \"mat\" undefined\n");
         rc = 1;
      end

      if (rc) return; end;

      ixelt2 = ixelt1 + nofelt_type - 1;
      ixelts = idxtype(ixelt1 : ixelt2);
      idelt  = eltids(ixelts);
      [idelt, ixs] = sort(idelt);

      for n = 1 : nofelt_type
          ix = ixelt1 + ixs(n) - 1;
          ix = idxtype(ix);
          id = idelt(n);
          nn = length(input_elem(ix).nodes);
          if (nn == 0)
             fprintf(msg, "*E* Nodes of element %6d not definied\n", id);
             rc = 1; 
          elseif (length(input_elem(ix).nodes) != nelnod)
             fprintf(msg, "*E* Element %6d has wrong number of nodes\n", id);
             rc = 1; 
          end
          if (needgeom && ! isstruct(input_elem(ix).geom))
             fprintf(msg, "*E* Element %6d: \"geom\" undefined or ", id);
             fprintf(msg, "not a structure\n");
             rc = 1;
          end
          if (needmat && ! isstruct(input_elem(ix).mat))
             fprintf(msg, "*E* Element %6d: \"mat\" undefined or ", id);
             fprintf(msg, "not a structure\n");
             rc = 1;
          end
          if (rc == 0)
             elem(nelt).id    = id;
             elem(nelt).nodes = input_elem(ix).nodes;
             if (needgeom)
                elem(nelt).geom  = input_elem(ix).geom;
             end
             if (needmat)
                elem(nelt).mat   = input_elem(ix).mat;
             end
             for k = 1 : nelnod - 1
                 for l = k + 1 : nelnod
                     if(elem(nelt).nodes(k) == elem(nelt).nodes(l))
                        fprintf(msg, "*E* Element %6d has duplicate node %6d\n",
                                elem(nelt).id, elem(nelt).nodes(k));
                        rc = 1;
                     end
                 end
             end
          end
          nelt++;
      end

  end

  if (rc) return; end

# Process element nodal points: get coordinates and connectivity

  coor = nodes.coor;

  for m = 1 : noftyp

      nofelt_type = types(m).nofelt;
      ix          = types(m).ixelt1;
      nelnod      = types(m).nelnod;
      ndofnod     = types(m).ndofnod;
      neldof      = nelnod * ndofnod;

      for n = 1 : nofelt_type

          nodno = lookup(nodes.ids, elem(ix).nodes, "m");
          ixb   = find(nodno == 0);

          if (isempty(ixb))
             elem(ix).nodes = nodno;
             elem(ix).coor  = coor(nodno, :);
             ixdof1 = (nodno' - 1) * mxdofpnt;
             ixdof2 = ixdof1 + ndofnod;
             ixdof1 += 1;
             ect = linspace(ixdof1, ixdof2, ndofnod)';
             elem(ix).ects = reshape(ect, neldof, 1)';
          else
             rc = 1;
             for k = ixb 
                 fprintf(msg, "*E* Element %6d references undefined node %6d\n",
                        elem(ix).id, elem(ix).nodes(k));

             end
          end

          ix++;

      end

  end

# Build element index

  index = zeros(nofelt, 3);
  index(:, 1 : 2) = [[elem.id]; [1 : nofelt]]';

  n1    = 1;

  for m = 1 : noftyp
      nofelt_type       = types(m).nofelt;
      n2                = n1 + nofelt_type - 1;
      index(n1 : n2, 3) = m;
      n1                = n2 + 1;
  end

  index = sortrows(index, 1);

# Check on duplicate element identifiers

  ixb = find(index(1 : nofelt-1, 1) == index(2 : nofelt, 1));
  if (! isempty(ixb))
     rc = 1;
     for n = ixb
         fprintf(msg, "*E* Duplicate element identifier %6d found\n",
                 index(n, 1));
     end
  end

  if (rc) return; end

# Build the element structure

  elements = struct("nofelt", nofelt, "noftyp", noftyp,
                    "types", types, "elem", elem, "index", index);

end
