% FEM initial
function [KE0,KE1,iK,jK,edofMat,nodes_in_ele] = FEAinit(nelx,nely, optKE,dx,dy,D0,D1)

    if nargin > 6
        % two mats
        KE0 = intKE(dx/2,dy/2,D0,optKE);
        KE1 = intKE(dx/2,dy/2,D1,optKE);
    else
        % one mat
        KE0 = intKE(dx/2,dy/2,D0,optKE);
        KE1 = 1;
    end        
    
    nodenrs = reshape(1:(1 + nelx) * (1 + nely), 1 + nely, 1 + nelx);
    edofVec = reshape(2 * nodenrs(1:end - 1, 1:end - 1) + 1, nelx * nely, 1);
    edofMat = repmat(edofVec,1,8)+repmat([0 1 2*nely+[2 3 0 1] -2 -1],nelx*nely,1);
    nodes_in_ele = edofMat(:, 2:2:8) ./ 2; % 4 nodes of a element

    iK = reshape(kron(edofMat, ones(8, 1))', 64 * nelx * nely, 1);
    jK = reshape(kron(edofMat, ones(1, 8))', 64 * nelx * nely, 1); 
end