import numpy as np

def directConstraintEnformcement(nDofPerNode,equationNumber_noConstraints, F_ext_nodesDofs,fixedBC_nodesDofs,constraints=None):
    #Performs the constraint enforcement using the direct method
    numEquations = equationNumber_noConstraints

    #Determine the total DOFs for loads and BCs
    nDofTot=max(max(sublist) for sublist in numEquations)+1
    F_ext_DOFs=[]
    F_ext_values=[]
    for inode in range(len(F_ext_nodesDofs)):
        node=F_ext_nodesDofs[inode][0]
        for iload in range(1,len(F_ext_nodesDofs[inode])):
            dof=F_ext_nodesDofs[inode][iload][0]
            value=F_ext_nodesDofs[inode][iload][1]
            dofTot=node*nDofPerNode+dof
            F_ext_DOFs.append(dofTot)
            F_ext_values.append(value)
    fixedDofs=[]
    for inode in range(len(fixedBC_nodesDofs)):
        node=fixedBC_nodesDofs[inode][0]
        for idof in range(1,len(fixedBC_nodesDofs[inode])):
            dof=fixedBC_nodesDofs[inode][idof]
            dofTot=node*nDofPerNode+dof
            fixedDofs.append(dofTot)


    if constraints is not None and len(constraints)>0:
        indexRemove = 0
        previousMasterGlobalDofTot=[] #without the removal of the previous slave dofs
        previousMasterGlobalDof=[] #with the removal of the previous slave dofs
        for line in constraints:
            masterNode = line[0]
            slaveNode = line[1]
            extraRemove = 0
            for i in range(2, len(line)):  # Skip the first two elements
                masterGlobalDofTot=masterNode * nDofPerNode + line[i]
                if masterGlobalDofTot in previousMasterGlobalDofTot:
                    masterGlobalDOF = previousMasterGlobalDof[previousMasterGlobalDofTot.index(masterGlobalDofTot)]
                else:
                    masterGlobalDOF = masterGlobalDofTot - indexRemove
               
                slaveGlobalDOF = slaveNode * nDofPerNode + line[i]  - indexRemove
                # Replace slave DOFs with master DOFs 
                for eq in numEquations:
                    for j in range(len(eq)):
                        if eq[j] == slaveGlobalDOF:
                            eq[j] = masterGlobalDOF
                for dof in F_ext_DOFs:
                    if dof == slaveGlobalDOF:
                        F_ext_DOFs[F_ext_DOFs.index(dof)] = masterGlobalDOF
                for dof in fixedDofs:
                    if dof == slaveGlobalDOF:
                        fixedDofs[fixedDofs.index(dof)] = masterGlobalDOF
                extraRemove+=1
                previousMasterGlobalDofTot.append(masterGlobalDofTot)
                previousMasterGlobalDof.append(masterGlobalDOF)
            if extraRemove >0:
                for eq in numEquations:
                    for j in range(len(eq)):
                        if eq[j] > slaveGlobalDOF:
                            eq[j] -= extraRemove   
                for dof in F_ext_DOFs:
                    if dof > slaveGlobalDOF:
                        F_ext_DOFs[F_ext_DOFs.index(dof)] -= extraRemove
                for dof in fixedDofs:
                    if dof > slaveGlobalDOF:
                        fixedDofs[fixedDofs.index(dof)] -= extraRemove         
            indexRemove += extraRemove

    # Remove duplicate indices in each line of numEquations
    for i in range(len(numEquations)):
        count = {}
        for index in numEquations[i]:
            if index in count:
                count[index] += 1
            else:
                count[index] = 1
        numEquations[i] = [index for index in numEquations[i] if count[index] == 1]

    nDofTot=max(max(sublist) for sublist in numEquations)+1
    F_ext=np.zeros(nDofTot)
    for i in range(len(F_ext_DOFs)):
        F_ext[F_ext_DOFs[i]]=F_ext_values[i]
    allDofs=np.arange(nDofTot)
    fixedDofs=np.array(fixedDofs)
    freeDofs=np.setdiff1d(allDofs, fixedDofs)

    return numEquations, F_ext, freeDofs, fixedDofs