function [numEquations, F_ext, freeDofs, fixedDofs]=directConstraintEnformcement(nDofPerNode,equationNumber_noConstraints, F_ext_nodesDofs,fixedBC_nodesDofs,constraints)

% Performs the constraint enforcement using the direct method
numEquations = equationNumber_noConstraints;

% Determine the total DOFs for loads and BCs
nDofTot = max(cellfun(@(x) max(x), numEquations)) + 1;
F_ext_DOFs = [];
F_ext_values = [];

for inode = 1:length(F_ext_nodesDofs)
    node = F_ext_nodesDofs{inode}{1};
    for iload = 2:length(F_ext_nodesDofs{inode})
        dof = F_ext_nodesDofs{inode}{iload}(1);
        value = F_ext_nodesDofs{inode}{iload}(2);
        dofTot = (node-1) * nDofPerNode + dof;
        F_ext_DOFs = [F_ext_DOFs, dofTot];
        F_ext_values = [F_ext_values, value];
    end
end

fixedDofs = [];
for inode = 1:length(fixedBC_nodesDofs)
    node = fixedBC_nodesDofs{inode}(1);
    for idof = 2:length(fixedBC_nodesDofs{inode})
        dof = fixedBC_nodesDofs{inode}(idof);
        dofTot = (node-1) * nDofPerNode + dof;
        fixedDofs = [fixedDofs, dofTot];
    end
end

if ~isempty(constraints)
    indexRemove = 0;
    previousMasterGlobalDofTot = []; % without the removal of the previous slave dofs
    previousMasterGlobalDof = []; % with the removal of the previous slave dofs
    for i = 1:length(constraints)
        line = constraints{i};
        masterNode = line(1);
        slaveNode = line(2);
        extraRemove = 0;
        
        for j = 3:length(line)  % Skip the first two elements
            masterGlobalDofTot = (masterNode-1) * nDofPerNode + (line(j));
            if ismember(masterGlobalDofTot, previousMasterGlobalDofTot)
                masterGlobalDOF = previousMasterGlobalDof(find(previousMasterGlobalDofTot == masterGlobalDofTot));
            else
                masterGlobalDOF = masterGlobalDofTot - indexRemove;
            end
                masterGlobalDOF=masterGlobalDOF(1);

            slaveGlobalDOF = (slaveNode-1) * nDofPerNode + (line(j)) - indexRemove;
            
            % Replace slave DOFs with master DOFs
            for eqIdx = 1:length(numEquations)
                eq = numEquations{eqIdx};
                for jj=1:length(eq)
                   if eq(jj)==slaveGlobalDOF
                      eq(jj)=masterGlobalDOF; 
                   end
                end
                numEquations{eqIdx}=eq;
            end
%             for iDof = 1:length(F_ext_DOFs)
%                 dof=F_ext_DOFs(iDof);
%                 if dof==slaveGlobalDOF
%                    F_ext_DOFs(iDof)=masterGlobalDOF; 
%                 end
%             end
            
            
            F_ext_DOFs(F_ext_DOFs == slaveGlobalDOF) = masterGlobalDOF;
            fixedDofs(fixedDofs == slaveGlobalDOF) = masterGlobalDOF;
            
            extraRemove = extraRemove + 1;
            previousMasterGlobalDofTot = [previousMasterGlobalDofTot, masterGlobalDofTot];
            previousMasterGlobalDof = [previousMasterGlobalDof, masterGlobalDOF];
        end
        
        if extraRemove > 0
            for eqIdx = 1:length(numEquations)
                eq = numEquations{eqIdx};
                for jj=1:length(eq)
                   if eq(jj)>slaveGlobalDOF
                      eq(jj)=eq(jj)-extraRemove; 
                   end
                end
                numEquations{eqIdx}=eq;
            end
            F_ext_DOFs(F_ext_DOFs > slaveGlobalDOF) = F_ext_DOFs(F_ext_DOFs > slaveGlobalDOF) - extraRemove;
            fixedDofs(fixedDofs > slaveGlobalDOF) = fixedDofs(fixedDofs > slaveGlobalDOF) - extraRemove;
        end
        indexRemove = indexRemove + extraRemove;
    end
end

% Remove duplicate indices in each line of numEquations
for i = 1:length(numEquations)
    % Find unique values and their counts
    [uniqueVals, ~, idx] = unique(numEquations{i}, 'stable');
    counts = histc(idx, 1:length(uniqueVals));

    % Keep only values that appear exactly once
    numEquations{i} = uniqueVals(counts == 1);
end

nDofTot = max(cellfun(@(x) max(x), numEquations));
F_ext = zeros(1, nDofTot);
for i = 1:length(F_ext_DOFs)
    F_ext(F_ext_DOFs(i)) = F_ext_values(i);
end

allDofs = 1:nDofTot;
fixedDofs = fixedDofs(:);
freeDofs = setdiff(allDofs, fixedDofs);

% Output
F_ext = F_ext';
freeDofs = freeDofs';
fixedDofs = fixedDofs';

end