% 2D discontinuous coefficients
% Ndsr                     - number of the discontinuity subregions
% Vdsr(k), k=1,...Ndsr     - value of the discontinuity in the sunregions
%  dsr(k,1:4), k=1,...Ndsr - the coordinates of the discontinuity subregions,
%              given as [x1,x2,y1,y2]
%        
function val=disco_rule(i,j,indicator,dsr,Vdsr,Ndsr,hstep)

      x = Xdir(i,hstep);
      y = Ydir(j,hstep);
      if indicator==1,
         for k=1:Ndsr
           if (((dsr(k,1)<x)&(x<=dsr(k,2)))  & ...
               ((dsr(k,3)<y)&(y<=dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      if indicator==2, 
         for k=1:Ndsr
          if (((dsr(k,1)<x)&(x<=dsr(k,2)))  & ...
              ((dsr(k,3)<=y)&(y<dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      if indicator==3,
         for k=1:Ndsr
          if (((dsr(k,1)<x)&(x<=dsr(k,2)))  & ...
              ((dsr(k,3)<=y)&(y<dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      if indicator==4, 
         for k=1:Ndsr
          if (((dsr(k,1)<=x)&(x<dsr(k,2)))  & ...
              ((dsr(k,3)<y)&(y<=dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      if indicator==5, 
         for k=1:Ndsr
          if (((dsr(k,1)<=x)&(x<dsr(k,2)))  & ...
              ((dsr(k,3)<y)&(y<=dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      if indicator==6, 
         for k=1:Ndsr
          if (((dsr(k,1)<=x)&(x<dsr(k,2)))  & ...
              ((dsr(k,3)<=y)&(y<dsr(k,4)))) 
                      val = Vdsr(k);
                      return
           end
         end
      end
      val = 1.0d0;
%