DBC_solve (K,L,force,in_data)

function [u] = DBC_solve (K,L,force,in_data)


dim    = size(K,1);
un = setdiff([1:dim],L);

un = setdiff(un,in_data.DBC.dofs);
DBC_dofs  = [in_data.DBC.dofs' ; un'];
DBC_displ = [in_data.DBC.displ' ; zeros(length(un),1)];

clear un
unknowns = setdiff([1:dim],DBC_dofs);

if  isfield(in_data,'slaves') 
    [T,slaves_yes] = MFC_ (K,L,force,in_data.slaves,in_data.master);
    K=T'*K*T;            
    force=T'*force';     
    [cc,idx] = setdiff(unknowns,slaves_yes);
    clear cc
    unknowns = unknowns(sort(idx));
    counter = 1;
    for i=unknowns unknowns(counter)=find(T(i,:)); counter=counter+1; end
    [cc,idx] = setdiff(DBC_dofs,slaves_yes);
    clear cc
    DBC_dofs=DBC_dofs(sort(idx)); 
    counter = 1;
    for i=DBC_dofs' DBC_dofs(counter)=find(T(i,:)); counter=counter+1; end
    m = length(DBC_dofs);
    n = length(unknowns);
    u = [ K(unknowns,unknowns)  spalloc(n,m,1); spalloc(m,n,1) speye(m,m) ]\...
        [force(unknowns)-K(unknowns,DBC_dofs)*DBC_displ; DBC_displ];
    u([unknowns' ; DBC_dofs]) = u;
    u = T*u;
else
    m = length(DBC_dofs);
    n = length(unknowns);

    u = [ K(unknowns,unknowns)  spalloc(n,m,1); spalloc(m,n,1) speye(m,m) ]\...
        [force(unknowns)'-K(unknowns,DBC_dofs)*DBC_displ; DBC_displ];
    u([unknowns' ; DBC_dofs]) = u;
end



function [T,slaves_yes] = MFC_ (K,L,force,slaves,master)

dim    = size(K,1);
T      = speye(dim);
not_slaves  = [1:dim];

for i=1:size(slaves,2)
    T(slaves(i).a,master(i).a)=1;
    not_slaves = setdiff(not_slaves,slaves(i).a);
end
slaves_yes = setdiff([1:dim],not_slaves);
T=T(:,not_slaves);