%  Calculate weights of shared boundary nodes between two domains 
function wght = WghtFDEqs(Domain)
  wght={};
  round5=@(z)eval(mat2str(z,5));
  Af=@(mn,MN)((1<mn(1))+(mn(1)<MN(1)))*((1<mn(2))+(mn(2)<MN(2)))/4;
  D=length(Domain);
  for i=1:D
    [xx{i}, yy{i}, rCw{i}]=Domain{i}{end-2:end};
    MN{i}=[length(yy{i}), length(xx{i})];
    dA{i}=abs(diff(xx{i}(1:2)))*abs(diff(yy{i}(1:2)));
    xx{i}=round5(xx{i});
    yy{i}=round5(yy{i});
  end
  nn=1;                   % unique node counter
  for i=1:D-1
    for j=i+1:D
      xint=intersect(xx{i},xx{j});
      yint=intersect(yy{i},yy{j});
      if isempty(xint) || isempty(yint)
        continue;
      end
      for x=xint
        for y=yint
          if nn>1 && ~isempty(find(all([x;y]==done)))
            continue;
          end
          dd=[i j];
          mn=[find(y==yy{i}), find(y==yy{j});
              find(x==xx{i}), find(x==xx{j})];
          ww=[dA{i}*Af(mn(:,1),MN{i})*rCw{i}, ...
              dA{j}*Af(mn(:,2),MN{j})*rCw{j}];
          for k=j+1:D    % check rest of regions
            xk=intersect(x,xx{k});
            yk=intersect(y,yy{k});
            if isempty(xk) || isempty(yk)
              continue;
            end
            dd(end+1)=k;
            mn(:,end+1)=[find(yk==yy{k}); find(xk==xx{k})];
            ww(end+1)=dA{k}*Af(mn(:,end),MN{k})*rCw{k};
          end
          ww=ww/sum(ww); % normalize weights
          done(:,nn)=[x;y]; % add to done list
          wght{nn}={dd;mn';ww};
          nn=nn+1;
        end
      end
    end
  end
end
