function [dofs, rc] = mfs_new_dofs(msg, mxdofpnt, nodes, elements, constraints,
                                   lincon)

# usage: [dofs, rc] = mfs_new_dofs(msg, mxdofpnt, nodes, elements, constraints,
#                                  lincon)
#
# Input  msg          File handle of message file 
#        mxdofpnt     Maximum number of degrees of freedom at a nodal point
#        nodes        Structure with nodal point data
#        elements     Structure with element data
#        constraints  Structure with constraints
#        lincon       Structure array with linear constraints
#
# Output dofs         Structure with degree of freedom data
#        rc           Return code: 0 means no erros
#                                  1 means errors
#
# The function builds a structure with degree of freedom data.
#
# Fields of structure dofs:
#
#    mxdofpnt          Maximum number of degrees of freedom at a nodal
#                      point
#    ndofg             Number of global degrees of freedom
#    ndofl             Number of local degrees of freedom
#    ndofp             Number of prescribed degrees of freedom
#    ndofd             Number of dependent degrees of freedom
#    ndofa             Number of autonomous degrees of freedom
#    dofl(ndofl)       Indices of local degrees of freedom within global
#                      degrees of freedom
#    dofp(ndofp)       Indices of prescribed degrees of freedom within
#                      global degrees of freedom
#    dofd(ndofd)       Indices of dependent degrees of freedom within
#                      global degrees of freedom
#    dofa(ndofa)       Indices of autonomus degrees of freedom within
#                      global degrees of freedom
#    C(ndofg, ndofg)   Constraint matrix
#
# -----------------------------------------------------------------------------

# Initialize

  rc    = 0;
  dofs  = struct();
  ndofp = 0;
  ndofd = 0;
  ndofa = 0;

  prescribed_motion = 0;

# Check arguments

  if (nargin != 6 || nargout != 2)
     print_usage();
  end

  isconst = isstruct(constraints);
  islinc  = isstruct(lincon);

# Initialize global degrees of freedom

  ndofg = nodes.nofnod * mxdofpnt;
  dofg  = zeros(ndofg, 1);

# Determine dependent and autonomous degrees of freedom

  if (islinc)
     dofax  = dofg; dofdx = dofg;
     noflnc = length(lincon);
     for n = 1 : noflnc
         ix = (lincon(n).dofd(1, :) - 1) * mxdofpnt;
         ix += lincon(n).dofd(2, :);
         bad = find(dofdx(ix)); nbad = length(bad);
         if (nbad)
             rc = 1;
             for m = 1 : nbad
                 j = bad(m);
                 k = lincon(n).dofd(1, j);
                 fprintf(msg, "*E* Node %6d, dof %1d dependent ",
                         nodes.ids(k), lincon(n).dofd(2, j));
                 fprintf(msg, "in more than one linear constraint\n");
             endfor
         endif
         dofdx(ix) = 1;
         lincon(n).ixd = ix;
         ix = (lincon(n).dofa(1, :) - 1) * mxdofpnt;
         ix += lincon(n).dofa(2, :);
         dofax(ix) = 1;
         lincon(n).ixa = ix;
     endfor
     dofd = find(dofdx); ndofd = length(dofd);
     dofax(dofd) = 0;
     dofa = find(dofax); ndofa = length(dofa);
  endif

# Mark available degrees of freedom

  elem = elements.elem;
  for n = 1 : elements.nofelt
      ects = elem(n).ects;
      dofg(ects) = 1;
  endfor

  if (islinc)
     dofg(dofa) =  1;
     dofg(dofd) = -1;
  endif

# Process constraints

  if (isconst)

#    Prescribed motion

     if (isfield(constraints, "prescribed"));
        prescribed_motion = 1;
        if (mfs_checkfields1(msg, constraints.prescribed, {"id", "dofs"},
            "Prescribed displacements"));
           rc = 1;
        else
           for p = constraints.prescribed
               if (p.id)
                  id = p.id;
                  n = lookup(nodes.ids, id, "m");
                  if (n)
                      dgfs = p.dofs;
                      if (dgfs)
                         if (max(dgfs) > mxdofpnt || min(dgfs) < 1)
                            fprintf(msg, "*E* Illegal dof identifer used ");
                            fprintf(msg, "to prescribe motion ")
                            fprintf(msg, "at node %6d\n", id);
                            rc = 1;
                         else
                            for k = dgfs
                                l = (n - 1) * mxdofpnt + k;
                                if (dofg(l) == 1 || dofg(l) == 2)
                                   dofg(l) = 2;
                                elseif (dofg(l) == -1)
                                   rc = 1;
                                   fprintf(msg, "*E* Prescribed motion ");
                                   fprintf(msg, "defined for dependent dof %1d ", k);
                                   fprintf(msg, "at node %6d\n", id); 
                                else
                                   fprintf(msg, "*W* Prescribed motion ");
                                   fprintf(msg, "defined for undefined dof %1d ", k);
                                   fprintf(msg, "at node %6d\n", id); 
                                endif
                            endfor
                         endif
                      else
                         fprintf(msg, "*E* Prescribed dofs undefined ");
                         fprintf(msg, "at node %6d\n", id);
                         rc = 1;
                      endif
                  else
                     fprintf(msg, "*E* Prescribed motion defined");
                     fprintf(msg, " for undefined node %6d\n", id);
                     rc = 1;
                  endif
               else
                  fprintf(msg, "*E* Prescribed displacements: ");
                  fprintf(msg, "Node not defined\n");
                  rc = 1;
               endif
           endfor
        endif
     endif

  endif

  if (rc) return; end

# Identify local and prescribed degrees of freedom

  dofl  = find(dofg == 1);
  ndofl = length(dofl);

  if (prescribed_motion)
     dofp  = find(dofg == 2);
     ndofp = length(dofp);
  endif

# Build the constraint matrix

  if (islinc)

     C = sparse(ndofg, ndofg);
     dofix = dofg;        dofix(dofd) = 0;
     dofi  = find(dofix); ndofi       = length(dofi);
     C(dofi, dofi) = speye(ndofi);

     for n = 1 : noflnc
         ixd         = lincon(n).ixd;
         ixa         = lincon(n).ixa;
         C(ixd, ixa) = lincon(n).C;
     endfor

     Cdd = C(dofd, dofd);
     while (nnz(Cdd))
        C(dofd, dofi) += Cdd * C(dofd, dofi);
        Cdd = Cdd * Cdd;
     endwhile
     C(dofd, dofd) = Cdd;

  endif

# Build the degree of freedom structure

  dofs  = struct("mxdofpnt", mxdofpnt, "ndofg", ndofg,
                 "ndofl", ndofl, "ndofp", ndofp, "ndofd", ndofd,
                 "ndofa", ndofa, "dofl", dofl);
  if (prescribed_motion)
     dofs.dofp = dofp;
  endif
  if (islinc)
     dofs.dofd = dofd;
     dofs.dofa = dofa;
     dofs.C    = C;
  endif

endfunction
