Added checker for counterclockwise

- When making global mesh it is possible that elements are counterclockwise in the mesh projection (stereographic) but will not be counterclockwise in a cylindrical projection oft used for simulation. Made a routine to check for this add change such elements.
CHLNDDEV · Feb 4, 2019 · fe78650 · fe78650
1 parent e0c466b
commit fe78650
Show file tree

Hide file tree

Showing 4 changed files with 321 additions and 371 deletions.
diff --git a/@meshgen/meshgen.m b/@meshgen/meshgen.m
@@ -659,14 +659,23 @@
             fprintf(1,' ------------------------------------------------------->\n') ;
 
             %% Doing the final cleaning and fixing to the mesh...
-            [p,t] = fixmesh(p,t);
-            % Put the mesh class into the grd part of meshgen
-            obj.grd.p = p; obj.grd.t = t;
             % Clean up the mesh if specified
             if obj.cleanup 
+                % Put the mesh class into the grd part of meshgen
+                obj.grd.p = p; obj.grd.t = t;
                 obj = clean(obj);
+            else
+                % Fix mesh on the projected space
+                [p1(:,1),p1(:,2)] = m_ll2xy(p(:,1),p(:,2)); 
+                [p,t1] = fixmesh(p1,t);
+                % Put the mesh class into the grd part of meshgen
+                obj.grd.p = p; obj.grd.t = t1;
             end
 
+            % Check element order, important for the global meshes crossing
+            % -180/180 boundary
+            obj.grd = CheckElementOrder(obj.grd);
+
             if obj.plot_on
                 figure; plot(obj.qual,'linewi',2);
                 hold on
@@ -741,8 +750,10 @@
                [obj.pfix(:,1),obj.pfix(:,2)] = ...
                                      m_ll2xy(obj.pfix(:,1),obj.pfix(:,2)); 
             end
+            % transform coordinates to projected space and "fix"
             [obj.grd.p(:,1),obj.grd.p(:,2)] = ...
                                    m_ll2xy(obj.grd.p(:,1),obj.grd.p(:,2)); 
+            [obj.grd.p,obj.grd.t] = fixmesh(obj.grd.p,obj.grd.t);                   
             if db
                 % Begin by just deleting poor mesh boundary elements
                 tq = gettrimeshquan(obj.grd.p,obj.grd.t);
@@ -807,7 +818,8 @@
             [obj.grd.p(:,1),obj.grd.p(:,2)] = ...
                                    m_xy2ll(obj.grd.p(:,1),obj.grd.p(:,2));
             if ~isempty(obj.pfix)
-             [obj.pfix(:,1),obj.pfix(:,2)]=m_xy2ll(obj.pfix(:,1),obj.pfix(:,2)); 
+                [obj.pfix(:,1),obj.pfix(:,2)] = ...
+                                     m_xy2ll(obj.pfix(:,1),obj.pfix(:,2)); 
             end
         end
 

diff --git a/@msh/msh.m b/@msh/msh.m
@@ -1295,6 +1295,47 @@ function write(obj,fname,type)
             warning('on','all')
         end
 
+        function obj = CheckElementOrder(obj,proj)
+            if nargin == 1
+               proj = 0; 
+            end
+            vx = obj.p(:,1); vy = obj.p(:,2);
+            xt = [vx(obj.t(:,1)) vx(obj.t(:,2)) ...
+                  vx(obj.t(:,3)) vx(obj.t(:,1))];
+            yt = [vy(obj.t(:,1)) vy(obj.t(:,2)) ...
+                  vy(obj.t(:,3)) vy(obj.t(:,1))];
+            dxt = diff(xt,[],2);
+            dyt = diff(yt,[],2);
+            for ii = 1:3
+                iii = ii - 1; if iii == 0; iii = 3; end
+                I = abs(dxt(:,ii)) > 180 & abs(dxt(:,iii)) > 180;
+                xt(I,ii) = xt(I,ii) + 360*sign(dxt(I,ii)); 
+            end
+            xt(:,end) = xt(:,1);
+            % Get new diff
+            dxt = diff(xt,[],2);
+            Area = dxt(:,3).*-dyt(:,2) + dxt(:,2).*dyt(:,3);
+            if ~isempty(find(Area < 0, 1))
+                disp([num2str(sum(Area < 0)) ...
+                 ' elements found that are not in counterclockwise order'])
+                if proj
+                    global MAP_PROJECTION MAP_VAR_LIST MAP_COORDS
+                    % kjr 2018,10,17; Set up projected space imported from msh class
+                    MAP_PROJECTION = obj.proj ;
+                    MAP_VAR_LIST   = obj.mapvar ;
+                    MAP_COORDS     = obj.coord ;
+                    [pt(:,1),pt(:,2)] = m_ll2xy(obj.p(:,1),obj.p(:,2));
+                else
+                    pt = obj.p;
+                end
+                [~,tt] = fixmesh(pt,obj.t);
+                obj.t(Area < 0,:) = tt(Area < 0,:);
+                obj = CheckElementOrder(obj,~proj);
+            else
+                disp('All elements now in counterclockwise order')
+            end
+        end
+
         function [out1,barlen,bars] = CalcCFL(obj,dt)
             g      = 9.81;        % gravity
             [bars,barlen] = GetBarLengths(obj);

diff --git a/utilities/m_tricontour.m b/utilities/m_tricontour.m
@@ -1,4 +1,4 @@
-function [C,h] = m_tricontour(tri,p,z,N)
+function H = m_tricontour(tri,p,z,N,C)
 %  M_CONTOURF Adds filled contours to a map
 %    M_CONTOUR(LONG,LAT,DATA,...)
 %
@@ -12,6 +12,6 @@
 end
 
 [X,Y] = m_ll2xy(p(:,1),p(:,2),'clip','on');
-hold on; [C,h] = tricontour([X,Y],tri,z,N);
+hold on; H = tricontour(tri,X,Y,z,N,C);
 
 end