function [D, penalties] = addClosuresToDiffOp(diffOp, D, BCs)
    penalties = cell(size(BCs));
    for i = 1:size(BCs,1)
        for j = 1:size(BCs,2)
            [closure, penalties{i,j}] = diffOp.boundary_condition(BCs{i,j}.boundary, BCs{i,j}.type);
            if ~isa(closure, 'function_handle')
                closure = @(v) closure*v;
            end
            D = @(v) D(v) + closure(v);
        end
    end
end