Merge pull request #22 from CHLNDDEV/wpringle-dpolypatch

Wpringle dpolypatch

@geodata/geodata.m

@@ -682,7 +682,7 @@ function plot(obj,type,projection)
             end
             bufx = 0.2*(obj.bbox(1,2) - obj.bbox(1,1));
             bufy = 0.2*(obj.bbox(2,2) - obj.bbox(2,1));
-            if startsWith(projection,'ste')
+            if ~isempty(regexp(projection,'ste'))
                 m_proj(projection,'lat',min(obj.bbox(2,:)),...
                     'long',mean(obj.bbox(1,:)),'radius',...
                     min(179.9,1.20*max(diff(obj.bbox(2,:)))));

@geodata/private/Read_shapefile.m

@@ -44,17 +44,19 @@
             end
             if minus 
                bboxt(:,1) = bboxt(:,1) - 360; 
-            end
-            % The shaperead is much faster if it is available
+            end          
             % Read the structure
             try 
+                % The shaperead is faster if it is available
                 S = shaperead(fname{1},'BoundingBox',bboxt);
                 % Get rid of unwanted components;
                 D = struct2cell(S);
                 S = cell2struct(D(3:4,:)',{'X','Y'},2);
                 disp('Read shapefile with shaperead')
                 sr = 1;
             catch
+                % If only m_shaperead is available or if some error occured
+                % with shaperead (e.g., 3D shapefile, m_shaperead may work)
                 disp('Reading shapefile with m_shaperead')
                 % This uses m_map (slower but free)
                 S = m_shaperead(fname{1},reshape(bboxt',4,1));
@@ -165,11 +167,14 @@
 new_islandb_type = []; new_mainb_type = [];
 for i = 1 : size(tmpC,1)
     if sr
+        % using shaperead
         points = tmpC{i,1}(1:end-1,:) ;
         In     = tmpInC{i,1}(1:end-1) ;
     else
+        % using m_shaperead
         points = tmpC{i,1}(1:end,:) ;
         if size(points,2) == 3
+            % if 3-D shapefile
             height = points(:,3); 
             points = points(:,1:2); 
             type   = tmpC{i,2};

@meshgen/private/dpoly.m

@@ -83,19 +83,20 @@
     % in_outer are the points inside the outer meshing domain.
     %
     % "in" is defined as the intersection of the in = in_boubox && in_outer
-    %
-    % Note: the boubox is pre-pended when forming outer.
-    firstNaN = find(isnan(outer(:,1)),1,'first') ;
-    in_boubox = inpoly(p(inside,:),outer(1:firstNaN-1,:)) ;
 
-    % Check if the points are inside the meshing domain (omega)
-    edges = Get_poly_edges( outer );
     if sum(inside)~=0
+        % Check if the points are inside the meshing domain (omega)
+        % Note: the boubox is pre-pended when forming outer.
+        firstNaN = find(isnan(outer(:,1)),1,'first') ;
+        in_boubox = inpoly(p(inside,:),outer(1:firstNaN-1,:)) ;
+        % now check if inside the full boundary
+        edges = Get_poly_edges( outer );
         in_outer    = inpoly(p(inside,:),outer,edges);
         if inpoly_flip
             in_outer = ~in_outer;
         end
     else
+        in_boubox   = d_l*0;
         in_outer    = d_l*0;
     end
 
@@ -106,4 +107,5 @@
 
     d(inside) = d_l;
 end
+% EOF
 end

@msh/msh.m

@@ -1002,9 +1002,7 @@ function write(obj,fname,type)
                     varargino{end+1} = 0;
                 end
                 obj = clean(obj,varargino(:));
-            end
-
-            if opt.nscreen
+            elseif opt.nscreen
                 disp(['number of nodes is ' num2str(length(obj.p))])
                 disp(['mean quality is ' num2str(mq_m)])
                 disp(['min quality is ' num2str(mq_l)])

utilities/Calc_f13_inpoly.m

@@ -8,6 +8,7 @@
 %               Attributes currently supported:
 %              'Cf' ('quadratic_friction_coefficient_at_sea_floor')
 %              'EV' ('average_horizontal_eddy_viscosity_in_sea_water_wrt_depth')
+%              'Adv' ('advection_state')
 %               etc. 
 %            3) A cell-arry of polygons in which you would like to alter 
 %               the attribute.
@@ -36,6 +37,9 @@
 elseif strcmpi(attribute,'surface_submergence_state')
     attrname = 'surface_submergence_state';
     default_val = 0;
+elseif strcmpi(attribute,'Adv')
+    attrname = 'advection_state';
+    default_val = -1e4;
 else
     error(['Attribute ' attribute ' not currently supported'])
 end