function [lincon, mxdofpnt, rc] = ...
mfs_new_lincon(msg, subtype, constraints, nodes)

# usage: [lincon, mxdofpnt, rc] = ...
#        mfs_new_lincon(msg, subtype, constraints, nodes)
#
# Input  msg          File handle of message file
#        subtype      Model subtype
#        constraints  Structure with contraints 
#        nodes        Structure with nodal point data
#
# Output lincon       Structure array with contraint matrices
#        mxdofpnt     Maximum number of degrees of freedom per nodal
#                     point
#        rc           Return code: 0 means no errors
#                                  1 means errors
#
# The function processes the linear constraints. It builds the constraint
# matrices and determines mxdofpnt based on the linear constraints. 
# mxdofpnt is later adjusted base on the elements.
#
# Fields of structure array lincon:
#
#    dofd(2, :)   Dependent degrees of freedom:
#                   internal nodal point identifier, dof identifier
#    dofa(2, :)   Autonomous degrees of freedom:
#                   internal nodal point identifier, dof identifier
#    C(:, :)      Constraint matrix of one constraint
#                 (rows = dep. dofs, cols = autonomous dofs)
#
# -----------------------------------------------------------------------------

  ctypes = {"prescribed"};

# Check arguments

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

# Initialize

  lincon   = 0;
  mxdofpnt = 0;
  rc       = 0;

# Get list of linear constraints implemented

  fname = strcat(mfilename(), ".m");
  bpath = file_in_loadpath(fname);
  len   = length(bpath) - length(fname);
  lcdir = strcat(bpath(1 : len), "lincon/", subtype, "/");

  lcfound = readdir(lcdir);
  lctypes = lcfound(3 : end);

# Check fields of structure constraints

  legal_fields   = [ctypes; lctypes];
  defined_fields = fieldnames(constraints);
  defleg         = ismember(defined_fields, legal_fields);

  illeg = find(! defleg);

  if ((lnilleg = length(illeg)))
     rc = 1;
     for n = 1 : lnilleg
         fprintf(msg, "*E* Unknown constraint %s found\n",
                 defined_fields{illeg(n)});
     end
  end

  if (rc) return; end

# Check if there are linear constraints

  deflin = ismember(defined_fields, lctypes);

  if (! sum(deflin)) return; end

# Process linear constraints

  lincon = struct();
  nlc    = 0;

  for n = 1 : length(deflin)
      if (deflin(n))
         lcpath = [lcdir, defined_fields{n}];
         lcinp  = getfield(constraints, defined_fields{n});
         addpath(lcpath);
         for l = 1 : length(lcinp)
             [lincon(++nlc), mxdp, rc1] = ...
                mfs_cs(msg, l, lcinp(l), nodes);
             rc = rc || rc1;
             mxdofpnt = max(mxdofpnt, mxdp);
         end
         rmpath(lcpath);
      end
  end

end
