% optBodyF = 1 : Null body force 
% optBodyF = 2 : Effective body force 
% optBodyF = 3 : Pointwise body force 
function [Ue] = microFEA_QBC(De, ke, iK, jK, edofMat, dx,dy, ...
        optBodyF, ufixed, freedofs, fixeddofs, xPhys, cx,cy, e, emin)    
    
    testNum = size(ufixed, 2);
    e0 = eye(testNum);  
    F = zeros(size(ufixed)); 
    gaussB = shapeDerivatives(dx,dy);
    jac = dx*dy/4;    
    
    [nely, nelx] = size(xPhys);
      
    D = 0;
    for i = 1:nelx
        for j = 1:nely       
            ele = j + (i-1)*nely;     
            De = De * xPhys(j,i);
                        
            if optBodyF==2
                D = D + De;

            elseif optBodyF==3
                D = De;

                fe = zeros(8,testNum);
                for gp = 1:4
                    for k = 4:testNum
                        eta = [e0(4, k), e0(8, k); 
                               e0(5, k), e0(7, k); 
                               e0(9, k), e0(6, k)] * [cx(j,i);cy(j,i)]; 
                        fe(:,k) = fe(:,k) + gaussB{gp}'*D*eta*jac;      
                    end
                end
                edof = edofMat(ele,:);
                F(edof,1:testNum) = F(edof,1:testNum) + fe;
            end
        end
    end
    
    if optBodyF==2
        D = D/nelx/nely;
        for i = 1:nelx
            for j = 1:nely       
                ele = j + (i-1)*nely;

                fe = zeros(8,testNum);
                for gp = 1:4
                    for k = 4:testNum
                        eta = [e0(4, k), e0(8, k); 
                               e0(5, k), e0(7, k); 
                               e0(9, k), e0(6, k)] * [cx(j,i);cy(j,i)];
                        fe(:,k) = fe(:,k) + gaussB{gp}'*D*eta*jac;      
                    end
                end

                edof = edofMat(ele,:);
                F(edof,1:testNum) = F(edof,1:testNum) + fe;
            end
        end
    end

    if nargin>14
        sK = ke(:) * (emin + xPhys(:)'*(e-emin)); 
    else
        sK = ke(:) * xPhys(:)'; % xPhys = sum(H(De)^2)/4, thus not need penalty
    end
    K = sparse(iK, jK, sK(:)); 
    K = (K + K') / 2;        

    fr = -K(freedofs, fixeddofs)*ufixed(fixeddofs,:) + F(freedofs,:);
    
    Ue = zeros(size(ufixed));
    Ue(fixeddofs, :) = ufixed(fixeddofs,:);
    Ue(freedofs, :)  = K(freedofs, freedofs) \ fr;     
    
%     % test
%     for j = 1:testNum
%         figure;       
%         coorx = reshape(Ue(1:2:end,j),nely+1,nelx+1) + cx;
%         coory = reshape(Ue(2:2:end,j),nely+1,nelx+1) + cy;
%         scatter(coorx(:),coory(:)); hold on;
%     end
end