function [cs, mxdofpnt, rc] = mfs_cs(msg, nc, lcinp, nodes)

# usage: [cs, mxdofpnt, rc] = mfs_cs(msg, nc, lcinp, nodes)
#
# Input  msg      File handle of message file
#        nc       Number of constraint
#        lcinp    Structure with definition of linear constraint
#        nodes    Structure with nodal point data 
#
# Output cs       Structure with constraint matrix
#        mxdofpnt Maximum number of degrees of freedom per nodal
#                 point
#        rc       Return code: 0 means no errors
#                              1 means errors
#
# The function computes the constraint matrix of a linear constraint
# of type "connect".
#
# Fields of structure lcinp:
#   dofs(:)     Degree of freedom identifiers
#   noda(:)     External nodal point identifiers of autonomous nodes
#   nodd(:)     External nodal point identifiers of dependent nodes
#
# Fields of structure cs: 
#   dofa(2, :)  Autonomous degrees of freedom:
#                 internal nodal point identifier, dof identifier
#   dofd(2, :)  Dependent degrees of freedom:
#                 internal nodal point identifier, dof identifier
#    C(:, :)   Constraint matrix of one constraint
#              (rows = dep. dofs, cols = autonomous dofs)
#
# -------------------------------------------------------------------

  legal_fields = {"dofs", "noda", "nodd"};
  nlegal = length(legal_fields);

# Initialize

  mxdofpnt = 0;
  cs       = struct("dofa", [], "dofd", [], "C", []);
  rc       = 0;

# Check for unknown fields

  fields = fieldnames(lcinp);
  defleg = ismember(fields, legal_fields);
  if (sum(defleg) != nlegal)
     rc = 1;
     for n = 1 : nlegal
         if (! defleg(n))
            fprintf(msg, "*E* constraints.connect(%3.0d): ", nc);
            fprintf(msg, "unknown field %s\n", fields{n});
         end
     end
  end

# Check degree of freedom identifiers

  if (! isfield(lcinp, "dofs") || isempty(lcinp.dofs))
     rc = 1;
     fprintf(msg, "*E* constraints.connect(%3d): ", nc);
     fprintf(msg, "dofs undefined\n");
  elseif ((mxdofpnt = max(lcinp.dofs)) > 6 || min(lcinp.dofs) < 1)
     fprintf(msg, "*E* constraints.connect(%3.0d): ", nc);
     fprintf(msg, "illegal dof identifier found\n");
  end

# Check identifiers of autonomous nodes

  if (! isfield(lcinp, "noda") || isempty(lcinp.noda))
     rc = 1;
     fprintf(msg, "*E* constraints.connect(%3d): ", nc);
     fprintf(msg, "noda undefined\n");
  else
     nnoda = length(lcinp.noda);
     ixa   = lookup(nodes.ids, lcinp.noda, "m");
     bad   = find(ixa == 0);
     nbad  = length(bad);
     if (nbad)
        rc = 1;
        for n = 1 : nbad
            fprintf(msg, "*E* constraints.connect(%3d).noda(%2d): ",
                    nc, bad(n));
            fprintf(msg, "node %6.0d does not exist\n", 
                    lcinp.noda(bad(n)));
        end
     end
  end

# Check identifiers of dependent nodes

  if (! isfield(lcinp, "nodd") || isempty(lcinp.nodd))
     rc = 1;
     fprintf(msg, "*E* constraints.connect(%3d): ", nc);
     fprintf(msg, "nodd undefined\n");
  else
     nnodd = length(lcinp.nodd);
     ixd   = lookup(nodes.ids, lcinp.nodd, "m");
     bad   = find(ixd == 0);
     nbad  = length(bad);
     if (nbad)
        rc = 1;
        for n = 1 : nbad
            fprintf(msg, "*E* constraints.connect(%3d).nodd(%2d): ",
                    nc, bad(n));
            fprintf(msg, "node %6.0d does not exist\n", 
                    lcinp.nodd(bad(n)));
        end
     end
  end

  if (rc) return; end

  if (nnoda != nnodd)
     rc = 1;
     fprintf(msg, "*E* constraints.connect(%3d): ", nc);
     fprintf(msg, "number of nodes does not match\n");
  end

# Check if nodal points are both autonomous and dependent

  bad  = intersect(lcinp.noda, lcinp.nodd);
  nbad = length(bad);
  if (nbad)
     rc = 1;
     for n = 1 : nbad
         fprintf(msg, "*E* constraints.connect(%3d): ", nc);
         fprintf(msg, "node %6.0d in both noda and nodd\n",
                 bad(n));
     end
  end

  if (rc) return; end

# Build structure with constaint matrix

  dofs  = sort(lcinp.dofs);
  ndofs = length(lcinp.dofs);
  ndoft = ndofs * nnoda;

  dofd = zeros(2, ndoft); dofa = zeros(2, ndoft);

  j = 0;
  for n = 1 : nnoda
      for m = 1 : ndofs
          dofa(1, ++j) = ixa(n);
          dofa(2, j)   = dofs(m);         
          dofd(1, j)   = ixd(n);
          dofd(2, j)   = dofs(m);         
      end
  end

  cs = struct("dofa", dofa, "dofd", dofd, "C", speye(ndoft));

end
