%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-06-16 14:52:29
 % @LastEditTime: 2021-06-16 15:48:00

 % compute NH of \Omega_{\alpha}
 % input: interpolation matrix R:[ndofs, v], which is ordered as col by col
 % output: surface NH on this fine mesh, [nnodes_mi, 2, 2, nnodes_ma]
%--------------------------
function NH = NH_alpha(microx, edofMat_mi, nnodes, Nh, R)

    NH = zeros((microx+1)^2, 2, 2, nnodes);
    for ele_mi = 1:microx^2
        edof = edofMat_mi(ele_mi, :);
        
        NH_tmp = zeros(4, 2, 2, nnodes);
        for gp = 1:4
            tmp = Nh{gp} * R(edof, :);      % [2, v+p]
            tmp1 = reshape(tmp', 2, [], 2); % [2, nnodes, 2]
            NH_tmp(gp, :, :, :) = permute(tmp1, [3,1,2]); % [2,2,nnodes]
        end
                    
        NH(edof(2:2:end)/2, :, :, :) = NH_tmp([1,3,4,2], :, :, :);
    end
end