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.

[edeyejedi]

Just wondering if there are any updates progress on this, or thoughts?
I consistently find myself in a loop of mesh generation in OM2D and testing in DFM.
The latest iteration of this had me visiting this thread again as my patience is running out!
I have previously used the "orthogonalize grid" operation in RGFGrid to just get a mesh working, but it feels a little corrupt as I have no idea what it is doing to depth values once it starts moving nodes around. My assumption is nothing, but I bring the output back in, re-interpolate the depths to be sure, then save off the files again.
It would be so cool (for me at least) to have a clear workflow from OM2D that appeases DFM.

[NadaSulaiman]

Hi. Unfortunately no progress on my end.
I have came to the conclusion that for DFLOWFM, the unstructured grid type (essentially a curvilinear grid tied together by triangles for refinement) is more stable and favorable than triangular grid. This paper discuss it in depth. https://www.researchgate.net/publication/225759861_Efficient_scheme_for_the_shallow_water_equations_on_unstructured_grids_with_application_to_the_Continental_Shelf

[krober10nd]

Hi, never had the motivation to finish this as new things always come through. It's unfortunate the scheme is so sensitive to the mesh. This is almost a fatal flaw IMO.

I found this repo recently though with this code that calculates orthogonality. Someone could translate it to matlab.

https://github.com/eigemx/neatmesh/blob/0633e3744b37b310b116189b62d707cf7d52a4c7/src/neatmesh/_analyzer.py#L119

Someone (myself included) needs to write a mesh non-orthongality calculator. The optimization approach would basically involve flipping edges and moving nodes in patches that don't meet the quality threshold.

[edeyejedi]

Thanks for the responses.
@NadaSulaiman that's understood. The curvilinear grids take a lot to put together as well, and do not always provide resolution where it is needed - which is one of the things that makes OM2D so good.
@krober10nd yeah I would love to help and dive in to this and suggest some fixes but unfortunately I am at the applied end, and my coding is not likely good enough :)
Re "flipping edges". I saw this in RGFGrid as an option "Flip Edges". Although "Flip Lines" in the manual:
Flip Lines
Minimise the number of edges connected to a node. The optimal number of edges to a node
is six.
Nodes that are connected to more than, say, six other nodes, are typically enclosed by cells
of highly non-uniform shape and wildly varying sizes. This motivates to improve the mesh
connectivity by selecting Operation→Flip Lines.

This is from the D-Flow manual:
Orthogonality is very important for accuracy: advised orthogonality values
for your grid are around 0.01, preferrably lower. The current treshold is already very
high at a value of 0.5. Use RGFGRID to improve your grid orthogonality (and smoothness).

I am guessing it would have to be a pretty big task to have an orthogonality check, constraint and adjustment step in the mesh generation.

[krober10nd]

Actually I don't think it would be too much work we just need to dedicate time. I also work on the applied side. I would be willing to set up a small project to work on this.

[edeyejedi]

Ok, sounds good. I spent some time yesterday using @NadaSulaiman's code, and some other bits, to inspect orthogonality values. As @NadaSulaiman mentioned, getting a match between RGFGRID and OM2D outputs is not straight forward and even then I still didn't get an exact match. Also, I spent ages manually adjusting a grid to get all orthog values below 0.5 and it still would not let me run the simulation. A function to established a rgfgrid-style_net.nc structure/variables with out manually saving off in RGFGrid may be a starting point. Happy to help where i can.

[edeyejedi]

Also found this that may provide some ideas: https://github.com/Deltares/MeshKernelPy/tree/main

[krober10nd]

One can compute mesh orthogonality as the dot product between the normal of each triangle's edge and the vector connecting each centroid to each adjacent triangle. This means that each triangle will have multiple orthongality values with the minimum value being perhaps recorded for each element.

A value of 1 for orthogonality for an element indicates all vectors described above are orientated in the same direction which is ideal. I assume then DelftFM measures non-orthogonality as they enforce a minimum of 0.5

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import matplotlib.colors as mcolors

# Function to compute the centroid of a triangle
def compute_centroid(vertices, element):
    x_coords = vertices[element, 0]
    y_coords = vertices[element, 1]
    centroid = np.array([np.mean(x_coords), np.mean(y_coords)])
    return centroid

# Function to compute the normal vector of an edge
def compute_normal(vertices, node1, node2):
    edge_vector = vertices[node2] - vertices[node1]
    normal_vector = np.array([-edge_vector[1], edge_vector[0]])
    return normal_vector / np.linalg.norm(normal_vector)

# Function to compute orthogonality for each edge of a triangle
def check_orthogonality(vertices, elements):
    orthogonality_results = []
    
    for element in elements:
        centroid = compute_centroid(vertices, element)
        
        for i in range(3):
            node1 = element[i]
            node2 = element[(i + 1) % 3]
            edge_centroid = (vertices[node1] + vertices[node2]) / 2
            normal_vector = compute_normal(vertices, node1, node2)
            centroid_vector = centroid - edge_centroid
            centroid_vector /= np.linalg.norm(centroid_vector)
            
            orthogonality = abs(np.dot(normal_vector, centroid_vector))  # Absolute value to ensure positive range
            orthogonality_results.append((node1, node2, orthogonality, edge_centroid, normal_vector, centroid_vector))
    
    return orthogonality_results

# Function to plot the mesh and highlight orthogonality
def plot_mesh(vertices, elements, orthogonality_results):
    fig, ax = plt.subplots()
    
    # Plot the triangular mesh
    triangles = tri.Triangulation(vertices[:, 0], vertices[:, 1], elements)
    ax.triplot(triangles, 'go-', lw=0.5)
    
    norm = mcolors.Normalize(vmin=0, vmax=0.5)
    sm = plt.cm.ScalarMappable(cmap='coolwarm', norm=norm)
    
    # Highlight edges based on orthogonality
    for result in orthogonality_results:
        node1, node2, orthogonality, edge_centroid, normal_vector, centroid_vector = result
        x = [vertices[node1, 0], vertices[node2, 0]]
        y = [vertices[node1, 1], vertices[node2, 1]]
        color = plt.cm.coolwarm(norm(orthogonality))  # Color from colormap based on orthogonality
        ax.plot(x, y, color=color, lw=2)
        ax.quiver(edge_centroid[0], edge_centroid[1], normal_vector[0], normal_vector[1], 
                  color='blue', scale=5, scale_units='xy')
        ax.quiver(edge_centroid[0], edge_centroid[1], centroid_vector[0], centroid_vector[1], 
                  color='green', scale=5, scale_units='xy')
    
    ax.set_aspect('equal')
    plt.colorbar(sm, ax=ax, label='Orthogonality (0 to 0.5)')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.title('Mesh Orthogonality with Vectors')
    plt.show()

# Example usage
vertices = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.5, np.sqrt(3)/2],
    [1.5, np.sqrt(3)/2],
    [0.5, 1.5]
])

elements = np.array([
    [0, 1, 2],  # Ideal orthogonality
    [1, 3, 2],  # Ideal orthogonality
    [0, 2, 4],  # Non-orthogonal element
    [2, 4, 3]   # Poor orthogonality
])

orthogonality_results = check_orthogonality(vertices, elements)
plot_mesh(vertices, elements, orthogonality_results)

[alexander0042]

Not directly related to OceanMesh, but I've worked out a solution to this using the Deltares MeshKernal/ UGrid tools, which contain routines for orthogonalization. So the process is:

1. Use the ever amazing OceanMesh2D to create an ADCIRC grid
2. Read that into Python with adcircpy
3. Create a UGrid mesh
4. Run an initial, whole domain orthogonalization (mesh2d_compute_orthogonalization)
5. Identify edges with poor orthogonality (mesh2d_get_orthogonality)
6. Loop through each edge:

• Wiggle nodes around
• Run a local orthogonalization over each problematic edge

8. Save back as UGrid for Delft3D-FM

I'll post my terribly Python below in case it's helpful to someone immediately, and hopefully post a repo with a bit more detail on this method. I'm happy to submit a PR here, but since it's Python and relies on a ton of (windows only) external libraries, not sure it's a good fit for OceanMesh.

#%% Script to try to orthoardinalize the an OcenaMesh2D Grid
# Alexander Rey, July 5, 2024

# Import
from meshkernel import (
    GeometryList,
    MeshKernel,
    OrthogonalizationParameters,
    ProjectionType,
    Mesh2d,
)
from ugrid import UGrid

import matplotlib.pyplot as plt
import numpy as np

#%% Read in the .14 mesh
from adcircpy import AdcircMesh
input_file = r"C:\Users\reya\Documents\Projects\Microplastics\Beaufort\Beaufort_TestO5_agg_01Lim.14"
mesh = AdcircMesh.open(input_file, crs=4326)

#%% Create a mesh2d kernel
mesh2d_mk = Mesh2d()
mesh2d_mk.node_x = mesh.x.to_numpy(np.float64)
mesh2d_mk.node_y = mesh.y.to_numpy(np.float64)
mesh2d_mk.face_nodes = mesh.triangles.to_numpy(np.int32).reshape(-1)
mesh2d_mk.edge_nodes = mesh.triangulation.edges.astype(np.int32).reshape(-1)

mk = MeshKernel(ProjectionType.SPHERICAL)
mk.mesh2d_set(mesh2d_mk)
mesh2d_output = mk.mesh2d_get()

#%% Get orthogonalization parameters
orthogonality = mk.mesh2d_get_orthogonality().values

selecting_polygon = GeometryList(
    np.empty(0, dtype=np.double), np.empty(0, dtype=np.double)
)
land_boundaries = GeometryList(
    np.empty(0, dtype=np.double), np.empty(0, dtype=np.double)
)

selecting_polygon.num_coordinates = 0
land_boundaries.num_coordinates = 0

orthoNodes = np.where(mk.mesh2d_get_orthogonality().values > 0.3)[0]
mk.mesh2d_delete_small_flow_edges_and_small_triangles(0.01,0.1)

while len(orthoNodes)>0:
    mk.mesh2d_compute_orthogonalization(
        project_to_land_boundary_option=0,
        orthogonalization_parameters=OrthogonalizationParameters(outer_iterations=2,
                                                                 boundary_iterations=10,
                                                                 inner_iterations=10,
                                                                 orthogonalization_to_smoothing_factor=0.95,
                                                                 orthogonalization_to_smoothing_factor_at_boundary=0.95,
                                                                 areal_to_angle_smoothing_factor=0.5),
        selecting_polygon=selecting_polygon,
        land_boundaries=land_boundaries,
    )

    mk.mesh2d_delete_small_flow_edges_and_small_triangles(0.01,
                                                          0.1)

    orthoVals  = mk.mesh2d_get_orthogonality().values
    orthoNodes = np.where( orthoVals> 0.3)[0]

    # Find the node with the highest orthogonality value
    orthoNodeMax = np.argmax(orthoVals)

    mesh2d_output = mk.mesh2d_get()

    # mk.mesh2d_delete_edge(mesh2d_output.edge_x[orthoNodeMax], mesh2d_output.edge_y[orthoNodeMax])
    # mk.mesh2d_delete_hanging_edges()
    print(orthoNodeMax)
    print(len(orthoNodes))
    print(orthoVals.max())

    # Make a polyon offset from the worst node

    # Plot the problem area
    fig, ax = plt.subplots()
    mesh2d_output.plot_edges(ax, color="black")
    ax.scatter(mesh2d_output.edge_x[orthoNodeMax], mesh2d_output.edge_y[orthoNodeMax], color="red")


    # Get ID of node with the highest orthogonality value
    nodeID = mesh2d_output.edge_nodes[orthoNodeMax*2]
    nodeIDb = mesh2d_output.edge_nodes[orthoNodeMax*2 + 1]

    # Move node a random amount
    edgeLength = np.sqrt((mesh2d_output.node_x[nodeID] - mesh2d_output.node_x[nodeIDb])**2 +
                            (mesh2d_output.node_y[nodeID] - mesh2d_output.node_y[nodeIDb])**2)
    delta_x = mesh2d_output.node_x[nodeID] + np.random.uniform(-edgeLength*0.1, edgeLength*0.1)
    delta_y = mesh2d_output.node_y[nodeID] + np.random.uniform(-edgeLength*0.1, edgeLength*0.1)

    # mk.mesh2d_move_node(delta_x, delta_y, nodeID)
    mesh2dMoved = Mesh2d()
    mesh2dMoved.node_x = mesh2d_output.node_x
    mesh2dMoved.node_y = mesh2d_output.node_y
    mesh2dMoved.face_nodes = mesh2d_output.face_nodes
    mesh2dMoved.edge_nodes = mesh2d_output.edge_nodes
    mesh2dMoved.face_edges = mesh2d_output.face_edges
    mesh2dMoved.node_x[nodeID] = delta_x
    mesh2dMoved.node_y[nodeID] = delta_y
    mk = MeshKernel(ProjectionType.SPHERICAL)
    mk.mesh2d_set(mesh2dMoved)


    selecting_polygon = GeometryList(
        x_coordinates=[mesh2d_output.edge_x[orthoNodeMax] - edgeLength*5,
                          mesh2d_output.edge_x[orthoNodeMax] + edgeLength*5,
                          mesh2d_output.edge_x[orthoNodeMax] + edgeLength*5,
                          mesh2d_output.edge_x[orthoNodeMax] - edgeLength*5,
                          mesh2d_output.edge_x[orthoNodeMax] - edgeLength*5],
        y_coordinates=[mesh2d_output.edge_y[orthoNodeMax] - edgeLength*5,
                            mesh2d_output.edge_y[orthoNodeMax] - edgeLength*5,
                            mesh2d_output.edge_y[orthoNodeMax] + edgeLength*5,
                            mesh2d_output.edge_y[orthoNodeMax] + edgeLength*5,
                            mesh2d_output.edge_y[orthoNodeMax] - edgeLength*5])


    ax.plot(selecting_polygon.x_coordinates, selecting_polygon.y_coordinates, color="red")

    # Plot the old and new node
    ax.scatter(mesh2d_output.node_x[nodeID], mesh2d_output.node_y[nodeID], color="blue")
    ax.scatter(mesh2d_output.node_x[nodeIDb], mesh2d_output.node_y[nodeIDb], color="blue")
    ax.scatter(delta_x, delta_y, color="green")


    # Zoom in
    ax.set_xlim([mesh2d_output.edge_x[orthoNodeMax] - edgeLength*10, mesh2d_output.edge_x[orthoNodeMax] + edgeLength*10])
    ax.set_ylim([mesh2d_output.edge_y[orthoNodeMax] - edgeLength*10, mesh2d_output.edge_y[orthoNodeMax] + edgeLength*10])
    plt.show()


#%% Save
mesh2d_ugrid = UGrid.from_meshkernel_mesh2d_to_ugrid_mesh2d(
    mesh2d=mk.mesh2d_get(), name="mesh2d", is_spherical=True
)

attribute_dict = {
    "name": "Unknown projected",
    "epsg": np.array([4326], dtype=int),
    "grid_mapping_name": "Unknown projected",
    "longitude_of_prime_meridian": np.array([0.0], dtype=float),
    "semi_major_axis": np.array([6378137.0], dtype=float),
    "semi_minor_axis": np.array([6356752.314245], dtype=float),
    "inverse_flattening": np.array([6356752.314245], dtype=float),
    "EPSG_code": "EPSG:4326",
    "value": "value is equal to EPSG code",
}
out_file = r"C:\Users\reya\Documents\Projects\Microplastics\Beaufort\PyOrthoB.nc"

with UGrid(out_file, "w+") as ug:
    # 1. Define a new mesh2d
    topology_id = ug.mesh2d_define(mesh2d_ugrid)
    # 3. Put a new mesh2d
    ug.mesh2d_put(topology_id, mesh2d_ugrid)
    # 3. Add crs to file
    ug.variable_int_with_attributes_define("wgs84", attribute_dict)
    # 4. Add conventions (global attributes)
    conventions = {
        "institution": "Deltares",
        "references": "Unknown",
        "source": "Unknown Unknown. Model: Unknown",
        "history": "Created on 2017-11-27T18:05:09+0100, Unknown",
        "Conventions": "CF-1.6 UGRID-1.0/Deltares-0.8",
    }
    ug.attribute_global_define(conventions)