Specifying mesh orthogonality threshold

[WPringle]

Delft FM has a strict grid orthogonality requirement to simulate.

Equation for orthogonality error: ....

Would like to be able to specify the maximum allowable orthogonality error for mesh generation and cleaning so that it can pass the Delft FM pre-simulation tests.

[NadaSulaiman]

• The orthogonality is defined as the cosine of the angle φ between a flowlink (a line segment connecting two cell circumcenters) and a netlink (an edge of a cell, connecting the corners of a cell). [Ref: D-FLOW Manual, p 167]

• Orthogonality threshold in DELFT FM cannot be bypassed. Current threshold is at 0.5. Advised orthogonality value for the grid should be lower than 0.01 for accuracy in solution. [Ref: D-FLOW Manual, p 129]

• Examples of perfect and poor orthogonality from the manual are below:

Note: For the example below, I have used the mesh generated by (Example_1_NZ.m) for the ease of reproducibility.

• Mesh was converted from .14 (ADCIRC) to .nc (D-FLOW FM Format).
• Colored dots represent orthogonality as calculated by the model:

Close ups of problematic elements (red), and the value of orthogonality displayed on the edge.

• The grid generation tool inside DELFT FM (RGFGRID) do offer an automatic function where the nodes are moved around to satisfy the minimum value, but it rarely gets the grid below the 0.5 threshold, and often distort the grid along the process, making it have worse properties, hence, the user is left to manually tweak the location of the problematic nodes to pass the orthogonality check, which is extremely time consuming.

• The model also will not run a grid which has circumcenters close to each other. I couldn’t find the specific value which the model calculates the distance based on, but the workaround suggested by the model is to use the “merge circumcenters” function, which leads to undesirable element shape in my opinion. [REF: RGFGRID Manual, p 63]

• Colored dots represent elements with close circumcenters

• Two pictures below are close ups of the elements before and after using the “merge circumcenters” tool.

[krober10nd]

Thank you Nada. We will need to implement this formula inside the msh class to test mesh generation and improvement strategies.

I would try reducing the gradation the g parameter to 0.15 in the mesh size function. This will improve the orthogonality constraint making it more amenable for simulation.

[CHLNDDEV]

I am attempting to compute the orthogonality of the New Zealand msh obj (acquired by running the example 1), however I am unsure which edge of the triangle is the "netlink" in your description.

Is the netlink the edge that intersects the flowlink? If so then my code below calculates the ortho. but it's almost near zero indicating a good mesh, but perhaps my calculation is incorrect in some way.

Note this script can be read with all the functionality inside MATLAB and the utils folder of OceanMesh2D on the general path.

@NadaSulaiman Any ideas what may be wrong here?

% Attempt at calculate ortho. of mesh cells
clear all; clc; clearvars; close all;

load nz

TR = triangulation(m.t,m.p);

C = circumcenter(TR);

ee = edges(TR);

orthogonality =  NaN(length(m.t),1);

for i = 1 : length(ee)
    
    startID = ee(i,1);
    endID   = ee(i,2);
    
    ID = edgeAttachments(TR,startID,endID);
    
    if length(ID{1}) > 1
        
        % flowlink (a line segment connecting two cell circumcenters)
        flowlink = [C(ID{1}(2),1),C(ID{1}(2),2), 0] - [C(ID{1}(1),1),C(ID{1}(1),2), 0];
        
        % and a netlink (an edge of a cell, connecting the corners of a cell)
        netlink = [m.p(endID,1),m.p(endID,2), 0] - [m.p(startID,1),m.p(startID,2), 0];
        
        % angle between two flowlink and netlink
        x = C(ID{1}(:),1);
        y = C(ID{1}(:),2);
        m1 = diff(y)./diff(x); 
        
        x = m.p(ee(i,:),1);
        y = m.p(ee(i,:),2);
        m2 = diff(y)./diff(x); 
        
        theta=atan(m1) - atan(m2);
        
        %theta=atan( ( m1 - m2 )./(1+(m1*m2)));
        rad2deg( theta )
        % orthogonality
        orthogonality(i)=abs(cos(theta));
        %
        disp(['The orthogonality is ',num2str(orthogonality(i))]);
        
        figure;
        triplot(m.t(ID{1},:),m.p(:,1),m.p(:,2));
        hold on; plot(m.p(ee(i,:),1),m.p(ee(i,:),2),'r-');
        hold on; plot(C(ID{1},1),C(ID{1},2),'g-s');
        title(['The orthogonality is ',num2str(orthogonality(i))])
        axis equal
        pause;
        close all;
        
        
    end
end

egfix_mid = (m.p(ee(:,1),:) + ...
    m.p(ee(:,2),:))/2;

figure; fastscatter(egfix_mid(:,1),egfix_mid(:,2),orthogonality)

[krober10nd]

I think the correct angle calculation is

 theta=atan(m1-m2);

where m1 and m2 are the slopes of the flowlink and netlink, respectively.

[krober10nd]

hm further inspecting this code a little bit more tonight, still coming up with low orthogonality (all < 0.5) values for a Delaunay triangulation of a random set of points.

% calculate ortho. of mesh cells
clc; clearvars; close all;

p = rand(100,2); 
t = delaunay(p); 
TR = triangulation(t,p); 

% from a triangulation generated via OM2d
%load nz
%p = m.p; 
%t = m.t;
%TR = triangulation(m.t,m.p);

C = circumcenter(TR);

ee = edges(TR);

orthogonality =  NaN(length(ee),1);

figure; triplot(TR); 

for i = 1 : length(ee)
    
    startID = ee(i,1);
    endID   = ee(i,2);
        
    ID = edgeAttachments(TR,startID,endID);
    
    if length(ID{1}) > 1        
        % angle between two flowlink and netlink
        x1 = C(ID{1},1);
        y1 = C(ID{1},2);
        
        x2 = p(ee(i,:),1);
        y2 = p(ee(i,:),2);
        
        % angle between flowlink and netlink
        theta = mod(atan2d(diff(y1),diff(x1)),360) - mod(atan2d(diff(y2),diff(x2)),360);
        % orthogonality
        orthogonality(i)=cosd(theta);
        
        if orthogonality(i) > 0.50
            
            figure;
            triplot(t(ID{1},:),p(:,1),p(:,2));
            hold on; plot(p(ee(i,:),1),p(ee(i,:),2),'r-');
            hold on; plot(C(ID{1},1),C(ID{1},2),'g-s');
            title(['The orthogonality is ',num2str(orthogonality(i))])
            axis equal
            pause;
            close all;
        end
        
        
    end
end

disp(['Maximum orthogonality is ',num2str(max(orthogonality))])
figure; histogram(orthogonality)

[NadaSulaiman]

Hi all, many thanks for the efforts.

I have adjusted the code to obtain identical orthogonality values as it was displayed inside the DELFT FM grid generation module (RGFGRID).

Summary of the steps:

• ADCIRC grid (.14) was converted to DELFT FM grid (_net.nc) using adcirc2dflowfm.m script. Note: to run this code, you need to download some dependencies, such as readNet.m and writeNet from this portal.

• The grid was loaded into RGFGRID (version: 6.07.00.72783,) and exported again as _net.nc file. Re-saving the mesh directly from RGFGRID seem to export more information inside the .nc file variables.

• Attached is the converted DELFT mesh (delftexport_net.nc). Please load into the directory to run the code below.
  delftexport_net.zip

% calculate ortho. of mesh cells
clc; clearvars; close all;
 
%p = rand(100,2); 
 p(:,1)=ncread('delftexport_net.nc','mesh2d_node_x')';
 p(:,2)=ncread('delftexport_net.nc','mesh2d_node_y')';
 
%t = delaunay(p); 
t_hold=ncread('delftexport_net.nc','mesh2d_face_nodes')';
t=t_hold(:,1:3);
 
TR = triangulation(t,p); 
 
% from a triangulation generated via OM2d
%load nz
%p = m.p; 
%t = m.t;
%TR = triangulation(m.t,m.p);
 
%C = circumcenter(TR);
C(:,1)=ncread('delftexport_net.nc','mesh2d_face_x')';
C(:,2)=ncread('delftexport_net.nc','mesh2d_face_y')';
 
%ee = edges(TR);
ee=double(ncread('delftexport_net.nc','mesh2d_edge_nodes')');
 
%orthogonality =  NaN(length(ee),1);
 
figure; triplot(TR); 
 
%Midpoint coordinates - for ploting
midpoint(:,1)=double(ncread('delftexport_net.nc','mesh2d_edge_x')');
midpoint(:,2)=double(ncread('delftexport_net.nc','mesh2d_edge_y')');
 
ID=ncread('delftexport_net.nc','mesh2d_edge_faces')';
indices = find(ID(:,2)==0);
ID(indices,:) = [];
ee(indices,:) = [];
midpoint(indices,:) = [];
 
 
orthogonality =  NaN(length(ID),1);
 
%for i = 1 : length(ee)
for i = 1 : length(ID)    
%     startID = ee(i,1);
%     endID   = ee(i,2);
%         
%     ID = edgeAttachments(TR,startID,endID);
    
%     if length(ID{1}) > 1        
 
        % angle between two flowlink and netlink
%         x1 = C(ID{1},1);
%         y1 = C(ID{1},2);
        x1 = C(ID(i,:),1);
        y1 = C(ID(i,:),2);
        
        x2 = p(ee(i,:),1);
        y2 = p(ee(i,:),2);
        
        % angle between flowlink and netlink
        theta = mod(atan2d(diff(y1),diff(x1)),360) - mod(atan2d(diff(y2),diff(x2)),360);
        % orthogonality
       % orthogonality(i)=cosd(theta);
        orthogonality(i)=abs(cosd(theta));
 
        
        if orthogonality(i) > 0.50
            
%             Uncomment to plot 
 
%             figure;
%             %triplot(t(ID{1},:),p(:,1),p(:,2));
%             triplot(t(ID(i,:),:),p(:,1),p(:,2));
%             hold on; plot(p(ee(i,:),1),p(ee(i,:),2),'r-');
%             %hold on; plot(C(ID{1},1),C(ID{1},2),'g-s');
%             hold on; plot(C(ID(i,:),1),C(ID(i,:),2),'g-s');
%             title(['The orthogonality is ',num2str(orthogonality(i))])
%             axis equal
%             pause;
%             close all;
 
 
 
        %end
        
        
    end
end
 
disp(['Maximum orthogonality is ',num2str(max(orthogonality))])
figure; histogram(orthogonality)
 
 
%Plot grid with orthogonality values
figure; triplot(TR); 
hold on
textscatter(midpoint(:,1),midpoint(:,2),string(round(orthogonality,3)))

In conclusion, the variables t, ee, and C which were imported directly from the DELFT FM grid seem to be different in values and/or dimensions from those computed in MATLAB. Very strange.

Here are some side-by-side close-ups to compare between the orthogonality obtained from RGFGRID and the modified MATLAB code. I believe they are identical.

[NadaSulaiman]

I am attempting to compute the orthogonality of the New Zealand msh obj (acquired by running the example 1), however I am unsure which edge of the triangle is the "netlink" in your description.

Is the netlink the edge that intersects the flowlink? If so then my code below calculates the ortho. but it's almost near zero indicating a good mesh, but perhaps my calculation is incorrect in some way.

Note this script can be read with all the functionality inside MATLAB and the utils folder of OceanMesh2D on the general path.

@NadaSulaiman Any ideas what may be wrong here?

% Attempt at calculate ortho. of mesh cells
clear all; clc; clearvars; close all;

load nz

TR = triangulation(m.t,m.p);

C = circumcenter(TR);

ee = edges(TR);

orthogonality =  NaN(length(m.t),1);

for i = 1 : length(ee)
    
    startID = ee(i,1);
    endID   = ee(i,2);
    
    ID = edgeAttachments(TR,startID,endID);
    
    if length(ID{1}) > 1
        
        % flowlink (a line segment connecting two cell circumcenters)
        flowlink = [C(ID{1}(2),1),C(ID{1}(2),2), 0] - [C(ID{1}(1),1),C(ID{1}(1),2), 0];
        
        % and a netlink (an edge of a cell, connecting the corners of a cell)
        netlink = [m.p(endID,1),m.p(endID,2), 0] - [m.p(startID,1),m.p(startID,2), 0];
        
        % angle between two flowlink and netlink
        x = C(ID{1}(:),1);
        y = C(ID{1}(:),2);
        m1 = diff(y)./diff(x); 
        
        x = m.p(ee(i,:),1);
        y = m.p(ee(i,:),2);
        m2 = diff(y)./diff(x); 
        
        theta=atan(m1) - atan(m2);
        
        %theta=atan( ( m1 - m2 )./(1+(m1*m2)));
        rad2deg( theta )
        % orthogonality
        orthogonality(i)=abs(cos(theta));
        %
        disp(['The orthogonality is ',num2str(orthogonality(i))]);
        
        figure;
        triplot(m.t(ID{1},:),m.p(:,1),m.p(:,2));
        hold on; plot(m.p(ee(i,:),1),m.p(ee(i,:),2),'r-');
        hold on; plot(C(ID{1},1),C(ID{1},2),'g-s');
        title(['The orthogonality is ',num2str(orthogonality(i))])
        axis equal
        pause;
        close all;
        
        
    end
end

egfix_mid = (m.p(ee(:,1),:) + ...
    m.p(ee(:,2),:))/2;

figure; fastscatter(egfix_mid(:,1),egfix_mid(:,2),orthogonality)

Apologies for the unclear description. Perhaps this figure from Bomers et al. (2019) better illustrates the idea.

[NadaSulaiman]

Alternatively, instead of saving the mesh as .14. and then converting it to _net.nc, we could write it directly as DELFT mesh by adding those three lines to the end of Example_1_NZ.m script. (writeNet.m needs to be downloaded first)

TR = triangulation(m.t,m.p); 
ee = edges(TR);
writeNet('NZ_net.nc',m.p(:,1),m.p(:,2),ee')

[NadaSulaiman]

Since

C = circumcenter(TR)

does not equal

C(:,1)=ncread('delftexport_net.nc','mesh2d_face_x')';
C(:,2)=ncread('delftexport_net.nc','mesh2d_face_y')';

Could anyone have a quick look at the definition of how the circumcenter of a triangle is computed from the D-Flow Flexible Mesh Technical Reference Manual (p.17) and confirm if it's the same method used for the MATLAB circumcenter function? I couldn't tell.

[krober10nd]

Hey Nada, I'll take a deeper look but there should be only one definition of the circumcenter. Perhaps the grid is projected into a map projection by Delft RF grid?

Could you share the netcdf files necessary to run the snippets of code you posted?

The element table may also be orientated in a different way by RFgrid

[NadaSulaiman]

Hi Keith. The grid is uploaded in the zip file here [delftexport_net.zip]. You should be able to run the code with it.

[krober10nd]

Indeed @NadaSulaiman the circumcenter calculation in MATLAB using the triangulation object produces different values for the circumcenters than what the utility you used above does. Interestingly some of the circumcenters are the same and some are not.
Blue asterisks are calculated in MATLAB and the red squares are what is found in the file "delftexport_net.zip". MATLAB's circumcenters are somehow always orthogonal by construction which seems incorrect to me.

[krober10nd]

I think DFM is placing the circumcenter on the boundary of the element if it's outside of the element while the MATLAB one is not.

[NadaSulaiman]

You are correct. I think it places it to the edge closest to the location of the circumcenter outside the element.
Suppose we adjust the code above to get the location of the modified circumcenter, could it then be added as a criteria when generating a mesh in OceanMesh2D?

[krober10nd]

@NadaSulaiman Yes, we can do that.

One naïve way would be:

1. Check if the circumcenter is inside the triangle (use https://www.mathworks.com/help/matlab/ref/triangulation.pointlocation.html)
2. If it isn't, find the closest edge to the circumcenter.
3. Assign the circumcenter to the nearest edge midpoint.
4. Calculate the orthogonality with the circumcenters.

It seems that if the circumcenter falls outside the cell, then DFM sets the triangle's circumcenter to the edge's midpoint nearest to the circumcenter.