anchor.diff:  END SUBROUTINE ElemM
   
   !------------------------------------------------------------------------------------------------------
  +!> Sets a list of DOF indices and the value these DOF should have
  +!! NOTE: need p%Reacts to have an updated first column that uses indices and not JointID
  +SUBROUTINE InitConstr(Init, p)
  +   TYPE(SD_InitType     ),INTENT(INOUT) :: Init
  +   TYPE(SD_ParameterType),INTENT(IN   ) :: p
  +   !
  +   INTEGER(IntKi) :: I,J
  +
  +   Init%BCs = 0
  +   DO I = 1, p%NReact
  +      DO J = 1, 6
  +         Init%BCs( (I-1)*6+J, 1) = (p%Reacts(I,1)-1)*6+J; ! DOF Index, looping through Joints in index order
  +         Init%BCs( (I-1)*6+J, 2) = p%Reacts(I, J+1);
  +      ENDDO
  +   ENDDO
  +END SUBROUTINE InitConstr
  +
   !> Apply constraint (Boundary conditions) on Mass and Stiffness matrices
   SUBROUTINE ApplyConstr(Init,p)
      TYPE(SD_InitType     ),INTENT(INOUT):: Init
  
positive.code:       ! deallocate temp matrices
        IF (ALLOCATED(TempPropsB))  DEALLOCATE(TempPropsB)
        IF (ALLOCATED(TempMembers)) DEALLOCATE(TempMembers)
        IF (ALLOCATED(TempReacts))  DEALLOCATE(TempReacts)
     END SUBROUTINE CleanUp_Discrt
  
  END SUBROUTINE SD_Discrt
  
negative.code: 
     ! --- Initialize Nodes
     Init%Nodes = 0   
     do I = 1,Init%NJoints
        Init%Nodes(I, 1) = I                ! JointID replaced by index I
        Init%Nodes(I, 2) = Init%Joints(I, 2)
        Init%Nodes(I, 3) = Init%Joints(I, 3)
        Init%Nodes(I, 4) = Init%Joints(I, 4)
     enddo
  
     ! --- Re-Initialize Reactions, pointing to index instead of JointID
     do I = 1, p%NReact
        JointID=p%Reacts(I,1)
        p%Reacts(I,1) = FINDLOCI(Init%Joints(:,1), JointID ) ! Replace JointID with Index
        if (p%Reacts(I,1)<=0) then
           CALL Fatal('Reaction joint table: line '//TRIM(Num2LStr(I))//' refers to JointID '//trim(Num2LStr(JointID))//' which is not in the joint list!')
           return
        endif
     enddo
  
     ! --- Re-Initialize interface joints, pointing to index instead of JointID
     Init%IntFc = 0
     do I = 1, Init%NInterf
        JointID=Init%Interf(I,1)
        Init%Interf(I,1) = FINDLOCI(Init%Joints(:,1), JointID )
        if (Init%Interf(I,1)<=0) then
           CALL Fatal('Interface joint table: line '//TRIM(Num2LStr(I))//' refers to JointID '//trim(Num2LStr(JointID))//' which is not in the joint list!')
           return
        endif
     enddo
  
     ! Change numbering in concentrated mass matrix
     do I = 1, Init%NCMass
        JointID = Init%CMass(I,1)
        Init%CMass(I,1) = FINDLOCI(Init%Joints(:,1), JointID )
        if (Init%Interf(I,1)<=0) then
           CALL Fatal('Concentrated mass table: line '//TRIM(Num2LStr(I))//' refers to JointID '//trim(Num2LStr(JointID))//' which is not in the joint list!')
           return
        endif
     enddo
  
  
     ! --- Initialize Elems, starting with each member as an element (we'll take NDiv into account later)
     p%Elems = 0
     ! --- Replacing "MemberID"  "JointID", and "PropSetID" by simple index in this tables
     DO iMem = 1, p%NMembers
        ! Column 1  : member index (instead of MemberID)
        p%Elems(iMem,     1)  = iMem
  