function [uglobal, deformed] = recover_full_u(lx, nelx, nely, nelz, microx, edofMat_ma, edofMat_mi, optDesign, ...
    vF, regulated_matrix, deform_nodes_full, dofid, U_bezier)
    
    globaly = nely * microx;
    globalx = nelx * microx;
    globalz = nelz * microx;
    [~, ~, edofMat_full] = prepare_assemble(globalx, globaly, globalz);
    uglobal = zeros(3 * (globalx + 1) * (globaly + 1) * (globalz + 1), 1);

    for i = 1:nelx
        for j = 1:nely
            for k = 1:nelz
                ele  = (k-1)*nelx*nely + (i-1)*nely + j;
                edof = edofMat_ma(ele, :);
                ue = U_bezier(edof);
                he = regulated_matrix{ele} * ue; % he order: [4-borders, inner]
                umicro = he(dofid);

                for ii = 1:microx
                    for jj = 1:microx
                        for kk = 1:microx
                            % which micro element in full domain
                            fulli = (i-1)*microx + ii;
                            fullj = (j-1)*microx + jj;
                            fullk = (k-1)*microx + kk;
                            global_ele = (fullk-1)*globalx*globaly + (fulli-1)*globaly + fullj;
                            global_edof = edofMat_full(global_ele, :);

                            % which micro element in this macro-ele
                            local_ele = (kk-1)*microx^2 + (ii-1)*microx + jj;
                            edof = edofMat_mi(local_ele, :);
                            uglobal(global_edof) = umicro(edof);
                        end
                    end
                end
            end
        end        
    end
    ref_nodes_full = mesh3D(nelx, nely, nelz, microx, 0, lx, 1);
    ufull = (deform_nodes_full - ref_nodes_full)';
    ufull = ufull(:);

    [freedofs, ~, ~, ~] = design_domain(1, optDesign, nelx, nely, nelz, microx, vF, 0, ref_nodes_full, lx, 1);
    avgerr = sum((ufull(freedofs) - uglobal(freedofs)).^2) ./ sum(ufull(freedofs).^2)

    deformed = reshape(uglobal, 3, [])' + ref_nodes_full;   
end