fixing for K of length 1

@msh/private/GridData.m

@@ -135,7 +135,7 @@
     disp(['Bounding maximum depth to ',num2str(maxdepth), ' meters.']) ;  
 end
 
-if strcmp(type,'slope')
+if strcmp(type,'slope') 
     disp(['Computing cell-averaged slope using the ' slope_calc ' method'])
     if isempty(obj.b) || all(obj.b == 0) 
         error(['You must have the bathymetry on the mesh before ' ...
@@ -145,6 +145,11 @@
                  'affect the computation of the slopes'])
     end
 end
+if length(K) == 1
+    if strcmp(type,'slope') || strcmp(type,'all')
+        error('Cannot compute slopes on bathymetry with K of length 1')      
+    end
+end
 
 if strcmp(NaNs,'fill') || strcmp(NaNs,'fillinside')
     disp('Fill in NaNs using nearest neighbour interpolation.')
@@ -266,19 +271,22 @@
 
 % delta is max and min bounds of the connecting element centers (or if only
 % one connecting element then do difference with the vertex itself
-
-% kjr 10/17/2018 enlarge the cell-averaging stencil by a factor N 
-pcon_max = max(squeeze(max(pcon,[],1)),2*obj.p(K,:)-squeeze(min(pcon,[],1)));
-pcon_min = min(squeeze(min(pcon,[],1)),2*obj.p(K,:)-squeeze(max(pcon,[],1)));
+if length(K) == 1
+    pcon_max = max(max(squeeze(pcon)),2*obj.p(K,:)-max(squeeze(pcon)));
+    pcon_min = min(min(squeeze(pcon)),2*obj.p(K,:)-min(squeeze(pcon)));
+else
+    pcon_max = max(squeeze(max(pcon,[],1)),2*obj.p(K,:)-squeeze(min(pcon,[],1)));
+    pcon_min = min(squeeze(min(pcon,[],1)),2*obj.p(K,:)-squeeze(max(pcon,[],1)));
+end
 delta = pcon_max - pcon_min;
 
+% kjr 10/17/2018 enlarge the cell-averaging stencil by a factor N 
 pcon_max = N*pcon_max+(1-N)*obj.p(K,:) ;
 pcon_min = N*pcon_min+(1-N)*obj.p(K,:) ;
 
-
 %% Now read the bathy (only what we need)
 if ~isa(geodata,'geodata')
-    BufferL = max(delta); 
+    BufferL = max(delta,[],1); 
     lon_min = min(obj.p(K,1)) - BufferL(1);
     lon_max = max(obj.p(K,1)) + BufferL(1);
     lat_min = min(obj.p(K,2)) - BufferL(2);
@@ -339,7 +347,6 @@
 end
 
 %% Calculate N - number of DEM grids to average - for each node
-tic
 if strcmp(interp,'CA')
     [~,IDXX,IDXY] = FindLinearIdx(obj.p(K,1),obj.p(K,2),DEM_X,DEM_Y);
     [~,IDXR,IDXT] = FindLinearIdx(pcon_max(:,1),pcon_max(:,2),DEM_X,DEM_Y);
@@ -437,7 +444,6 @@
         b = F(obj.p(K,1),obj.p(K,2));
     end
 end
-toc
 %% Put in the global msh class
 if strcmp(type,'depth') || strcmp(type,'all')
     if isempty(obj.b)

utilities/FindLinearIdx.m

@@ -9,6 +9,11 @@
 nx = size(lon,2);
 np = numel(x);
 
+if ny == 1 && ny == 1
+   IX = 1; IX1 = 1; IX2 = 1;
+   return
+end
+
 % make sure entry points are column vectors
 x = x(:);
 y = y(:);