function phi1 = forcelip(K,phi,m,n) n1 = n-1; ksig = (K>0); phiks = pl(phi-ksig'); K1 = getK1(K,m,n); Kf = getK1(phiks'.*K,m,n); phi1 = zeros(1,m^n1); iip = find(K1>0); phi1(iip) = Kf(iip)./K1(iip); return function phib = termwisebd(K,phi,m,n) n1 = n-1; global Ainl Binl Ain1 Bin1 GG iif m n mylinopt % assume K>0 K1 = getK1(K,m,n); Kf = getK1(phi'.*K,m,n); phib = zeros(m^n1,1); if sum(K1)>1e-6 phib = min(Kf ./ (K1 + ~K1),n1); end %phib = ceil(phib); dims = repmat([m],[1 n1]); for val = n1:-1:1 %xind = find(phib==val); xind = find(ceil(phib)==val); for i=xind(:)' v = phib(i); x = ind2coord(i,dims); Y = allrho1(x,m); Yind = coord2ind(Y,dims); badY = find(v-phib(Yind) > 1); phib(badY) = 0*phib(badY) + v - 1; end end phib = phib'; return