%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-01-10 11:15:50
 % @LastEditTime: 2021-04-04 16:37:42

 % optBodyF = 1 : Null body force 
 % optBodyF = 2 : Effective body force 
 % optBodyF = 3 : Pointwise body force 
%--------------------------
function [Ue] = microFEA_QBC_2mat(D0, D1, ke0, ke1, iK, jK, edofMat, dx,dy, ...
                optBodyF, ufixed, freedofs, fixeddofs, xPhys, cx, cy)

    testNum = size(ufixed, 2);
    e0 = eye(testNum);  
    F  = zeros(size(ufixed)); 
    
    gaussB = shapeDerivatives(dx/2, dy/2,2);
    jac = dx*dy/4;    

    [nely, nelx] = size(xPhys);
    sK = zeros(64, nelx*nely);
    
    D = 0;
    for i = 1:nelx
        for j = 1:nely       
            ele = j + (i-1)*nely;   
            if xPhys(ele)==1
                De = D1;
                sK(:, ele) = ke1(:);
            else
                De = D0;
                sK(:, ele) = ke0(:);
            end
                        
            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

    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; axis equal;
%         end
end