pymeshmaker.py

# -*- coding: utf-8 -*-

""" Shapefile to .Geo file converter

Module to convert a shapefile '.shp' to a geometry file of gmsh '.geo'.
!!! IMPORTANT !!! It can handle inner polygons BUT it can't place boundaries in the inner polygons!

author: Ignace Pelckmans
            (University of Antwerp, Belgium)
"""

import os
import pickle
import time
import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from shapely.geometry import Point
import pandas as pd
import geopandas as gpd

matplotlib.use('Agg')

from shp2mpol import *
from pol2pol import *
from remove_vertices import *
from remove_duplicates import *
from reverse_polygon import *
from polplot import *
from np2bgfield import *

class pymesh:
    """
    Class which describes an a domain (geometry) and evenually a mesh (.msh) generated by gmsh

    author: Ignace Pelckmans
                        (University of Antwerp)
    """

    def __init__(self, name = ''):
        # attribute with the name of the mesh
        self.name = name

        # automatically set the directory to the default directory
        self.setDirectory()
        # automatically set the defaultb boundary and break lists (empty)
        self.setBoundariesAndBreaks()

    def setDirectory(self, dir = '/home/ignace/Documents/Ph D/TELEMAC/Playing_Around/First_Steps/Test_Cases/'):
        """
        Method to set a directory to save all files

        Args:
            dir: (Optional, defaults to .../Test_Cases/input/v1/) string with the directory to save all generated files into

        """

        self.dir = dir

    def loadShp(self, fn, proj = 32717, remove_vertices = True, remove_duplicates = True, reverse_polygon = False):
        """
        Method to load a shapefile (.shp) and tranform to shapely Polygon, A Numpy array with dimension n x 2 with the
        exterior coordinates and a list of Numpy arrays with dimensions m x 3 with all coordinates of the inner rings.

        Args:
            fn: (Required) string with path to shapefile
            proj: (Optional, defaults to 32717) epsg code of the desired projection (hint: use metric projections, for instane, avoid epsg:4326)
                   Default is set at 32717 which is the epsg code for WGS 84 / UTM 17S.
            remove_vertices: (Optional, defaults to True) True/False boolean to indicate whether excessive vertices should be removed
                             Excessive is defined as vertices where there is no change of direction.
            remove_duplicates: (Optional, defaults to True) True/False boolean to indicate whether duplica vertices should be removed
        """

        # load .shp and reproject
        pol, epsg = shp2Mpol(fn, return_coordinate_system = True)
        if epsg != proj:
            print('reprojecting...')
            pol = pol2Pol(pol, epsg, proj).geoms[0]
            print('reprojected!')

        if isinstance(pol, MultiPolygon):
            pol = pol.geoms[0]

        # remove excessive vertices
        if remove_vertices: pol = removeVertices(pol)
        # Remove duplicates (highly recommended!)
        if remove_duplicates: pol = pol.simplify(0)

        if reverse_polygon:
            pol = reversePolygon(pol)

        # convert to Numpy Array and rotate 90° do the shape becoms n x 2, remove z-coordinates
        exter = np.rot90(pol.exterior.xy)[:-1,0:2]
        inter = [np.rot90(i.xy)[:-1,0:2] for i in pol.interiors]

        # number of exterior vertices
        n = len(exter)

        # set attributes
        self.pol = pol
        self.exter = exter
        self.inter = inter
        self.n = n
        self.n_inner_rings = len(inter)
        self.proj = proj
    
    def loadGPKG(self, fn, proj = 0, remove_vertices = True, remove_duplicates = True, reverse_polygon = False):
        """
        Method to load a geopackage (.gpkg) and tranform to shapely Polygon, A Numpy array with dimension n x 2 with the
        exterior coordinates and a list of Numpy arrays with dimensions m x 3 with all coordinates of the inner rings.

        Args:
            fn: (Required) string with path to shapefile
            proj: (Optional, defaults to 32717) epsg code of the desired projection (hint: use metric projections, for instane, avoid epsg:4326)
                   Default is set at 32717 which is the epsg code for WGS 84 / UTM 17S.
            remove_vertices: (Optional, defaults to True) True/False boolean to indicate whether excessive vertices should be removed
                             Excessive is defined as vertices where there is no change of direction.
            remove_duplicates: (Optional, defaults to True) True/False boolean to indicate whether duplica vertices should be removed
        """

        # load .gpkg
        gdf = gpd.read_file(fn, engine='pyogrio')
        epsg = int(gdf.crs.srs.split(':')[-1])

        # if need, reproject
        if epsg != proj and proj != 0:
            gdf = gdf.to_crs(epsg = proj)
            print('reprojected!')

        roi = gdf.iloc[0]
        pol = roi.geometry

        if isinstance(pol, MultiPolygon):
            pol = pol.geoms[0]

        # remove excessive vertices
        if remove_vertices: pol = removeVertices(pol)
        # Remove duplicates (highly recommended!)
        if remove_duplicates: pol = pol.simplify(0)

        if reverse_polygon:
            pol = reversePolygon(pol)

        # convert to Numpy Array and rotate 90° do the shape becoms n x 2, remove z-coordinates
        exter = np.rot90(pol.exterior.xy)[:-1,0:2]
        inter = [np.rot90(i.xy)[:-1,0:2] for i in pol.interiors]

        # number of exterior vertices
        n = len(exter)

        # set attributes
        self.pol = pol
        self.exter = exter
        self.inter = inter
        self.n = n
        self.n_inner_rings = len(inter)
        self.proj = proj

    def loadPolygon(self, pol, proj = 32717, remove_duplicates = True):
        """
        Method to load a Shapely polygon.

        Args:
            pol: (Required) Shapely polygon
            remove_duplicates: (Optional, defaults to True) True/False boolean to indicate whether duplica vertices should be removed
        """

        if isinstance(pol, MultiPolygon):
            pol = pol.geoms[0]

        # Remove duplicates (highly recommended!)
        if remove_duplicates: pol = pol.simplify(0)

        # convert to Numpy Array and rotate 90° do the shape becoms n x 2, remove z-coordinates
        exter = np.rot90(pol.exterior.xy)[:-1,0:2]
        inter = [np.rot90(i.xy)[:-1,0:2] for i in pol.interiors]

        # number of exterior vertices
        n = len(exter)

        # set attributes
        self.pol = pol
        self.exter = exter
        self.inter = inter
        self.n = n
        self.n_inner_rings = len(inter)
        self.proj = proj

    def findBoundaries(self, boundaries = []):
        """
        Method to find boundaries based on x,y pairs
        :param boundaries: list of boundaries where each boundary follows the format: [[x1,y1],[x2,y2]]
        :return:
        """
        bs = []

        for boundary in boundaries:
            b0x, b0y = boundary[0]
            b1x, b1y = boundary[1]
            i0 = np.argmin((self.exter[:,0]-b0x)**2 + (self.exter[:,1]-b0y)**2)+1
            i1 = np.argmin((self.exter[:,0]-b1x)**2 + (self.exter[:,1]-b1y)**2)+1
            i0,i1 = np.sort([i0,i1])

            bs.append([i0,i1])
        
        boundaries = sorted(bs, key=lambda x: x[0])
        self.setBoundariesAndBreaks(boundaries = boundaries)

    def setBoundariesAndBreaks(self, boundaries = [], breaks = []):
        """
        Method to set the boundaries and breaks as attributes:

        Args:
            boundaries: (Optional, defaults to an empty list) List of pairs of node id's (starting & ending node) of the boundaries
                        Node id's can be looked up through the plotting method
            breaks: (Optional, defaults to an empty list) List of node id's to force on the geometry
        """

        self.boundaries = boundaries
        self.breaks = breaks

    def vert2shp(self, fn):
        """
        Method to create a shapefile containing points with their labels

        Args:
            fn: (Required) filename/directory path string to store the shapefile
        """

        # put all coordinates in one array
        coor = np.asarray(self.exter)
        for i in self.inter:
            coor = np.vstack((coor, np.asarray(i)))

        # create a CRS object
        source = osr.SpatialReference()
        source.ImportFromEPSG(self.proj)

        # Now convert it to a shapefile with OGR
        driver = ogr.GetDriverByName('Esri Shapefile')
        if fn[-4:] != '.shp': fn += '.shp'
        ds = driver.CreateDataSource(fn)
        layer = ds.CreateLayer('', source, geom_type = ogr.wkbPoint)
        # Add one attribute
        layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
        defn = layer.GetLayerDefn()

        for i in range(coor.shape[0]):
            # Create a new feature (attribute and geometry)
            feat = ogr.Feature(defn)
            feat.SetField('id', i+1)

            # Make a geometry, from Shapely object
            p = Point(coor[i,0], coor[i,1])
            geom = ogr.CreateGeometryFromWkb(p.wkb)
            geom.AssignSpatialReference(source)
            feat.SetGeometry(geom)

            layer.CreateFeature(feat)
            feat = geom = None  # destroy these

        # Save and close everything
        ds = layer = feat = geom = None

    def plotShp(self, figfile=None, plot_inner=True, plot_vert=True, plot_vertlabels=False,
                plot_each_x_labels=10, figure_size=(50, 50), font_size=10, set_bounds=False, axis=None):
        """
        Method to plot the original shapefile (after reprojecting)

        Args:
            figfile: (Optional) Path (string) to save the figure. Required if no axis is provided.
            plot_inner: (Optional) True/False to plot inner rings
            plot_vert: (Optional) True/False to show vertices
            plot_vertlabels: (Optional) True/False to show vertex labels
            plot_each_x_labels: (Optional) Int for label interval
            figure_size: (Optional) Tuple for figure size
            font_size: (Optional) Font size for vertex labels
            set_bounds: (Optional) False or list of bounds [min_lat, max_lat, min_lon, max_lon]
            axis: (Optional) Matplotlib axis to plot into; if None, creates a new figure
        """
        is_external_axis = axis is not None
        if not is_external_axis:
            fig, axis = plt.subplots(figsize=figure_size)
            axis.set_aspect('equal')

        # plot polygon
        if plot_inner:
            polPlot(self.exter, XY_inner=self.inter, show_vertices=plot_vert, empty=False, plot_on_axis=axis,
                    show_vertices_labels=plot_vertlabels, show_vertices_labels_interval=plot_each_x_labels,
                    vertices_color='maroon', font_size=font_size)
        else:
            polPlot(self.exter, show_vertices=plot_vert, empty=False, plot_on_axis=axis,
                    show_vertices_labels=plot_vertlabels, show_vertices_labels_interval=plot_each_x_labels,
                    vertices_color='maroon', font_size=font_size)

        # indicate boundaries
        for b in self.boundaries:
            X, Y = [], []
            for i in range(b[0], b[1] + 1):
                x, y = self.exter[i - 1]
                X.append(x)
                Y.append(y)
            axis.plot(X, Y, color='purple', label='boundaries')

        # indicate break points
        for b in self.breaks:
            x, y = self.exter[b - 1]
            axis.scatter(x, y, c='pink', label='breaks')

        if set_bounds and isinstance(set_bounds, (list, tuple)) and len(set_bounds) == 4:
            axis.set_xlim(set_bounds[2], set_bounds[3])
            axis.set_ylim(set_bounds[0], set_bounds[1])

        # Save figure only if we created one
        if not is_external_axis and figfile:
            fig.savefig(figfile)
            plt.close(fig)

    def writeGeoFile(self, filename, cell_size = 250, spline_max_len = False, inner_rings = True, include_inner_rings = False,
                    XY_field = ' ', background_file = '', multiplier = 1, minres = 25, maxres = 1e9,
                    zerovalues = 250, outside_mesh_size = 250, plot_arr = False,
                    buffer_interpolation = 50, algo = 'Frontal-Delaunay', line_type = 'Splines', no_mesh_def = False,
                    prevent_add_nodes_to_edge = False, create_background_file = True, add_physical_lines = True,
                    physical_line_labels = ['closed', 'TidalEntrance', 'ClosedInner'], inner_rings_splines = True,):
        """
        Method to write a .geo file

        Args:
            filename: (Required) Directory path (string) to store the .geo file
            cell_size: (Optional, defaults to 250) Int or string 'hetero' to indicate the type of meshing, if number is given the mesh will be homogeneous
            spline_max_len: (Optional, defaults to False) False/int to indicate the max length of a spline (if False, there is no limit)
            inner_rings: (Optional, defaults to True) True/False boolean to add the points and splines of the inner rings
            include_inner_rings: (Optional, defaults to False) True/False to include the inner rings in the domain.
                                !!! different combinations of the parameters inner_rings and include_inner_rings lead to different results:
                                - inner_rings = True and include_inner_rings = False:
                                                    Inner rings will be holes (or islands) in the final mesh
                                - inner_rings = True and include_inner_rings = True:
                                                    Inner rings won't be holes but the mesh will be forced on the edges of the inner ring so
                                                    that the mesh is completely the same + mesh in the holes as it would have been with the holes
                                - inner_rings = False and include_inner_rings = True/False:
                                                    Inner rings will be ignored.
            XY_field: (Required if cell_size = 'hetero') directory path string to a pickle which hold a list with the following elements:
                        - 2d numpy array representing a spatial raster with values for the mesh size
                        - coordinate pair with the coordinates of the left bottom corner of the 2d numpy array
                        - x resolution
                        - y resolution
            multiplier: (Optional, defaults 1) value to multiply the original array with
            minres: (Optional, defaults to 25) minimum mesh size
            zerovalues: (Optional, defaults to 250) the mesh size in areas with value 0
            outside_mesh_size: (Optional, defaults to 250) the mesh size outside the background field
            plot_arr: (Optional, defaults to False) False/file directory to store a plot of the original raster
            buffer_interpolation: (Optional, defaults to False) False/float to indicate the size of buffer around the channels to interpolate
            algo: (Optional, defaults to Frontal-Delaunay) String (Frontal-Delaunay or Delaunay) to set the 2D meshing algoritm
            line_type: (Optional, defaults to Splines) 'Splines' or 'Lines' to indicate the type lines used !!! only works if there are no boundaries, otherwise it's splines !!!
            no_mesh_def: (Optional, defaults to False) True/False boolean which if True, ignores all mesh size fields
            prevent_add_nodes_to_edge: (Optional, defaults to False) True/False boolean which if True, prevents adding nodes to the domain edge
            create_background_file: (Optional, defaults to True) True/False boolean which indicates whether to create a background file
            add_physical_lines: (Optional, defaults to True) True/False boolean which indicates whether to create physical lines
            physical_line_labels: (Optional, defaults to ['closed', 'TidalEntrance', 'ClosedInner']) dictionary of label names for the physical lines.
            inner_rings_splines: (Optional, defaults to True) True/False boolean to indicate wheter to use splines (True) or lines (False) to create the inner rings
        """
        # set attributes
        self.algo = algo
        if type(cell_size) == int: self.holescellsize = cell_size
        else: self.holescellsize = zerovalues
        self.inner = inner_rings


        # Open txt file and allow both writing and reading
        f = open(filename, "w+")
        self.geofile = filename

        # write header to indicate the right engine
        f.write("SetFactory(\"OpenCASCADE\");\n")
        f.write("Geometry.OCCParallel = 1;\n")
        f.write("General.NumThreads = 6;\n//+\n")

        # initialize counters for each features
        self.n_Point = 1
        self.n_Spline = 1
        self.n_LineLoop = 1
        self.n_Surface = 1

        # attribute if there are inner rings
        self.inner_rings = inner_rings

        # reorganize the boundary node id's, exterior node id's and break id's to put make sure the first boundary node has id 1
        self.reorganize()

        # write the point lines
        f = self.writeGeoFileAddPoints(f)

        # write the spline/line lines
        if line_type == 'Splines':
            f = self.writeGeoFileAddSplines(f, spline_max_len = spline_max_len, line_type = line_type)
        else:
            f = self.writeGeoFileAddLines(f)

        # write points and inner rings
        if inner_rings:
            f = self.writeGeoFileInnerRings(f, inner_rings_splines = inner_rings_splines)

        # write physical boundaries
        if add_physical_lines:
            f = self.writeGeoFilePhysicalLines(f, physical_line_labels = physical_line_labels)

        # A priori resolution
        f = self.writeGeoFileAPrioriResolution(f)

        # if indicated precent adding nodes to the domain edges
        if prevent_add_nodes_to_edge:
            f.write('Transfinite Line "*" = 2;\n')

        # Line loop and surface of exterior ring
        f = self.writeGeoFileLineLoopsSurface(f, inner_rings = inner_rings)

        print('background file: ', background_file)
        # set mesh cell size
        if not no_mesh_def:
            if type(cell_size) == str:
                if cell_size == 'hetero':
                    if len(background_file) > 0: create_background_file = False;
                    f = self.writeGeoFileHeterogeneousMesh(f, XY_field, multiplier = multiplier, minres = minres, maxres = maxres,
                        zerovalues = zerovalues, outside_mesh_size = outside_mesh_size, plot_arr =  plot_arr,
                        buffer_interpolation = buffer_interpolation, create_background_file = create_background_file, background_file = background_file)
            elif type(cell_size) == float or type(cell_size) == int:
                f = self.writeGeoFileHomogeneousMesh(f, mesh_size = cell_size)

        # set meshing algorithm
        f = self.writeGeoFileDefineAlgorithm(f, algorithm = algo)

        # close the .geo file
        f.close()

    def reorganize(self):
        """
        Method to reorganize XY list, break list and boundaries to put boundary nodes first

        There is a spline-break at the first vertex anyway, so if there is a boundary, place the first vertex of that boundary
        as the first vertex so the break will be at a boundary. To do so, we go over all vertices id's and store them in a list
        when we pass the vertex id which is the first vertex of the first boundary, we start storing the vertex id's in a new list.
        The first list will now be stiched after the second one.

        For example, if we have an exterior of 10 vertices and there is one boundary between 4 - 5:
        original exterior vertex id's list: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
        new exterior vertex id's list: [4, 5, 6, 7, 8, 10, 1, 2, 3]

        Note that the inner rings does not have to be reorganized since they cannot hold a boundary per definition.

        Args:
            -
        """
        # only proceed if there are boundaries
        if len(self.boundaries) > 0:

            # id of the first vertex of the first boundary
            b0 = self.boundaries[0][0]

            # initiate two lists
            exter_re1, exter_re2 = [], []

            # loop over all vertices in the exterior
            for i in range(len(self.exter)):
                # if i + i (vertex id's start at 1) is smaller than the first vertex of the first boundary ...
                if i + 1 < b0:
                    # ... append that vertex to the first list
                    exter_re1.append(self.exter[i])
                else:
                    # otherwise, append to the other list
                    exter_re2.append(self.exter[i])

            # rearange order by placing the second list first
            exter_re2.extend(exter_re1)
            # exter is now the updated exter list
            self.exter = exter_re2

            # this means that the id's are reorganized and so, also the boundaries and breaks should be reorganized
            boundaries_upd = []
            # reorganize all boundary id's
            #
            # example: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] with boundaries at 4 - 5 and 8 - 10
            # Here, the new boundaries will be: 1 - 2 and 5 - 7

            for b_l, b_r in self.boundaries:
                boundaries_upd.append([b_l - b0 + 1, b_r - b0+1])
            self.boundaries = boundaries_upd

            # similar for all break points
            # however, now break id's can be smaller than 0 and thus - b0 + 1 can result in updated breaks equal to 0 or 1
            # in that case, add the total number of exterior breaks
            #
            # example: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
            #              *     ^''^  *     ^'''''^
            # with boundaries at 4 - 5 and 8 - 10 and breaks at 2 and 6
            #
            # updated example: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
            #                   ^''^  *     ^'''''^     *
            # Here, the new boundaries will be: 1 - 2 and 5 - 7
            # and the new breaks will be: 3
            breaks_upd = []

            breaks = [b-b0+1 for b in self.breaks]
            self.breaks = [b+n if b < 1 else b for b in breaks ]

    def writeGeoFileAddPoints(self, f):
        """
        Method to add all points of the exterior points to the .geo file

        Args:
            f: (Required) open text file representing the .geo file
        """

        # add a header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Points\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')

        # each exterior vertex should be added as a Point in the format:
        # Point(1) = {x, y, 0}
        # Point(2) = {x, y, 0}
        # ...

        # loop over all points in the exterior
        for i in range(np.shape(self.exter)[0]):
                f.write("Point(%d) = {%f, %f, 0};\n" % (self.n_Point, self.exter[i][0], self.exter[i][1]))
                # keep track of the number of points
                self.n_Point += 1

        return f

    def writeGeoFileAddSplines(self, f, spline_max_len = False, line_type = 'Splines'):
        """
        Method to write the splines and lins of the exterior in the .geo file

        Args:
            f: (Required) open text file representing the .geo file
            spline_max_len: (Optional, defaults to False) False/int to indicate the max length of a spline (if False, there is no limit)
            line_type: (Optional, defaults to Splines) 'Splines' or 'Lines' to indicate the type lines used !!! only works if there are no boundaries, otherwise it's splines !!!
        """

        # add a header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Lines & Splines\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')

        # initate local variables
        boundaries = self.boundaries
        breaks = self.breaks
        n_Spline = self.n_Spline
        n_Point = self.n_Point

        # if there are boundaries
        if len(boundaries) > 0:
            # initialize list to store the id's of the boundaries (or 'physical lines')
            physical_lines = []

            # start with first boundary since the first vertices will always be boundary vertices in case of boundaries

            # first boundary
            bound1 = boundaries[0]
            # length of first boundary
            bound1diff = bound1[1]-bound1[0]

            # write a line object in the .geofile
            # format example:
            # Line(1) = {1, 2}
            # where the numbers inbetween curly brackets refer to .geo point id's
            f.write('Line(%d) = {1, ' % n_Spline)

            # add this id to the list of physical boundaries
            physical_lines.append(n_Spline)

            # add a line per segment of the first boundary
            # for instance if there is a boundary 1 - 3:
            # Line(1) = {1, 2}
            # Line(2) = {2, 3}
            for i in range(0,bound1diff-1):
                f.write('%d};\n' % (i+2))
                # keep track of number of lines and splines
                n_Spline += 1
                f.write('Line(%d) = {%d,' % (n_Spline, i+2))
                # add all these lines to the list of physical boundaries
                physical_lines.append(n_Spline)

            # finish the open Line definition
            f.write('%d};\n' % (bound1diff+1))
            n_Spline += 1
            # write the last line of the first boundary
            f.write('Spline(%d) = {%d, ' % (n_Spline, bound1diff+1))
            n_Spline += 1

            # variable to keep track if the loop is inside a boundary or not
            b = False
            # initiate a variable to count the number of segments in a spline
            spline_counter = 0
            # loop over all points except for the points in the first boundary
            for i in range(bound1diff+2,n_Point):
                # if the point is the starting vertex of a boundary, close whatever was open en start a new line
                # Example:
                #
                # before: Spline(n) = {1, 2, 3 .
                # newly added:
                # , 4};
                # Line(n + 1) = Line{4,
                #
                if list(np.asarray(boundaries)[:,0]).__contains__(i):
                    # start a new line
                    f.write("%d};\nLine(%d) = {%d, " % (i, n_Spline, i))
                    physical_lines.append(n_Spline)
                    n_Spline += 1
                    # the loop is inside a boundary
                    b = True

                    # if the boundary only exist out of one vertex
                    if list(np.asarray(boundaries)[:, 1]).__contains__(i+1):
                        b = False

                # if the vertex is part of a boundary
                elif b:
                    # close the upper line and start a new one
                    # Example:
                    #
                    # before:
                    # Line(n + 1) = Line{4,
                    # new:
                    # Line(n + 1) = Line{4, 5};
                    # Line(n + 2) = Line{5,
                    f.write('%d};\nLine(%d) = {%d,' % (i, n_Spline, i))
                    physical_lines.append(n_Spline)
                    n_Spline += 1
                    # close the boundary flag if this vertex is the second last vertex of a boundary
                    if list(np.asarray(boundaries)[:, 1]).__contains__(i+1):
                        b = False

                # if the vertex is the closing vertex of a boundary
                elif list(np.asarray(boundaries)[:, 1]).__contains__(i):
                    # close the upper line and start a new one
                    # Example:
                    #
                    # before:
                    # Line(n + 2) = Line{5,
                    # new:
                    # Line(n + 2) = Line{5, 6};
                    # Line(n + 3) = Line{6,
                    f.write("%d};\nSpline(%d) = {%d," % (i, n_Spline,i))
                    n_Spline += 1
                    # if this point is also the last point
                    # Example:
                    #
                    # before:
                    # Line(n + 3) = Line{6,
                    # new:
                    # Line(n + 3) = Line{6, 7};

                    if i == n_Point-1:
                        f.write("%d, 1};\n" % (n_Point-1))


                # if this point is also the last point (and no part of a boundary)
                # Example:
                #
                # before:
                # Line(n + 3) = Line{6,
                # new:
                # Line(n + 3) = Line{6, 7};
                elif i == n_Point-1:
                    f.write("%d, 1};\n" % (n_Point-1))

                # if the vertex is just a regular vertex and no part of a boundary
                else:
                    # count the number of vertices in a spline
                    spline_counter += 1
                    # if there is a max spline length indicated
                    if spline_max_len:
                        # if the number of vertices in the spline is lower than than the allowed number
                        if spline_counter <= spline_max_len:
                            # if the vertex is a break
                            if breaks.__contains__(i):
                                # close the upper spline and start a new one
                                f.write("%d};\nSpline(%d) = {%d, " % (i, n_Spline, i))
                                n_Spline += 1
                            else:
                                # if the vertex is not a break and the number of vertices is lower than the allowed number, just add the point to the open spline
                                f.write("%d, " % (i))
                        else:
                            # if the number of vertices in this spline would exceed the allowed number, start a new spline
                            f.write("%d};\nSpline(%d) = {%d, " % (i, n_Spline, i))
                            n_Spline += 1
                            spline_counter = 0
                    # if there is not a maximum number of vertices in a spline indicated, check if it is a break vertex
                    elif breaks.__contains__(i):
                            f.write("%d};\nSpline(%d) = {%d, " % (i, n_Spline, i))
                            n_Spline += 1
                    # if there is not a maximum number of vertices in a spline and if it is not a break vertex, just add it to the open spline
                    else:
                        f.write("%d, " % (i))

        # if there are no boundaries indicates
        else:
            if line_type == 'Splines':
                # create first spline
                f.write('Spline(1) = {')
                # add all points to this spline, except for breaks. In that case, start a new spline
                for i in range(1,n_Point):
                    # in case of breaker vertices
                    if breaks.__contains__(i):
                            f.write("%d};\nSpline(%d) = {%d, " % (i, n_Spline, i))
                            n_Spline += 1
                    f.write('%d,' % (i))
                f.write('%d};\n' % (1))
                n_Spline += 1
            elif line_type == 'Lines':
                # create first Line
                f.write('Line(1) = {1,')
                n_Spline += 1
                # add all points to this spline, except for breaks. In that case, start a new spline
                for i in range(2,n_Point-1):
                    f.write("%d};\nLine(%d) = {%d, " % (i, n_Spline, i))
                    n_Spline += 1
                f.write('%d};\n' % (n_Point-1))
                f.write('Line(%d) = {%d,1};\n' % (n_Spline, n_Point-1))
                n_Spline += 1


        # globalize local variables
        self.n_Spline = n_Spline
        if len(boundaries) > 0:
            self.physical_lines = physical_lines

        return f

    def writeGeoFileAddLines(self, f):
        """
        Method to write and lines of the exterior in the .geo file

        Args:
            f: (Required) open text file representing the .geo file
          """

        # add a header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Lines & Splines\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')


        physical_lines = []

        for i in range(1, self.n_Point-1):
            f.write('Line(%d) = {%d, %d};\n' % (i,i,i+1))

            for boundary in self.boundaries:
                if i >= boundary[0] and i+1 <= boundary[1]:
                    physical_lines.append(i)


        f.write('Line(%d) = {%d, 1};\n' % (i+1,i+1))

        # globalize local variables
        self.n_Spline = i+2
        if len(self.boundaries) > 0:
            self.physical_lines = physical_lines

        return f

    def writeGeoFileInnerRings(self, f, inner_rings_splines = True):
        """
        Method to include inner rings (holes) in the .geo file

        Args:
            f: (Required) open text file representing the .geo file
            inner_rings_splines: (Optional, defaults to True) True/False boolean to indicate to use splines for the inner ring, the alternative are lines
        """

        # initate local variables
        boundaries = self.boundaries
        breaks = self.breaks
        n_Spline = self.n_Spline
        n_Point = self.n_Point

        # check if there are inner rings
        if self.inner:
            # write a header
            f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                    '//+ Holes Points\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                    '++++++++++++++++++\n')

            # save all inner points in a list per ring
            Inner_Points = []
            # loop over all inter rings
            for I in self.inter:
                Inner_Points.append([])
                # loop over all vertices in an inner ring
                for i in range(0,len(I)):
                    # add as a point
                    f.write("Point(%d) = {%f, %f, 0};\n" % (n_Point, I[i][0], I[i][1]))
                    # save the point id's in a list of list per inner ring
                    Inner_Points[-1].append(n_Point)
                    n_Point += 1

            # start a new counter to include all inner rings
            n_Spline_inter = n_Spline

        if self.inner:
            # create header
            f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                    '//+ Holes Splines\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                    '++++++++++++++++++\n')

            # initiate a list to store all inner spline ids
            Inner_Splines_ids = []
            Inner_Splines_ids_perLL = []
            # create a spline per inner ring (no boundaries are allowed at this point)
            # create splines if indicated
            if inner_rings_splines:
                for I in Inner_Points:
                    # open a spline
                    f.write('Spline(%d) = {' % (n_Spline_inter))
                    # add the id to the list of inner splines
                    Inner_Splines_ids.append(n_Spline_inter)
                    Inner_Splines_ids_perLL.append([n_Spline_inter])
                    # count
                    n_Spline_inter += 1
                    # add all points to this spline
                    for i in I:
                        if not i in self.breaks or i == I[0]:
                            f.write('%d,' % (i))
                        else:
                            f.write('%d};\n' % (i))
                            f.write('Spline(%d) = {%d,' % (n_Spline_inter, i))
                            # count
                            Inner_Splines_ids.append(n_Spline_inter)
                            Inner_Splines_ids_perLL[-1].append(n_Spline_inter)
                            n_Spline_inter += 1
                    # end with the first point of the ring
                    f.write('%d};\n' % (I[0]))

            # otherwise, create rings of lines
            elif not inner_rings_splines:
                for I in Inner_Points:
                    # open a spline
                    f.write('Line(%d) = {%d,' % (n_Spline_inter, I[0]))
                    Inner_Splines_ids.append(n_Spline_inter)
                    Inner_Splines_ids_perLL.append([n_Spline_inter])
                    n_Spline_inter += 1
                    # add all points to this spline
                    for n in range(1,len(I)):
                        i = I[n]
                        f.write('%d};\nLine(%d) = {%d,' % (i, n_Spline_inter, i))
                        Inner_Splines_ids.append(n_Spline_inter)
                        Inner_Splines_ids_perLL[-1].append(n_Spline_inter)
                        n_Spline_inter += 1
                    # end with the first point of the ring
                    f.write('%d};\n' % (I[0]))

        # globalize local variables
        if self.inner:
            self.Inner_Points = Inner_Points
            self.Inner_Splines_ids = Inner_Splines_ids
            self.n_Spline_inter = n_Spline_inter
            self.Inner_Splines_ids_perLL = Inner_Splines_ids_perLL

        self.n_Spline = n_Spline
        self.n_Point = n_Point

        return f

    def writeGeoFilePhysicalLines(self, f, physical_line_labels = ['closed', 'TidalEntrance', 'ClosedInner']):
        """
        Method to set the physical properties of the boundaries

        Args:
            f: (Required) open text file representing the .geo file
            physical_line_labels: (Optional, defaults to {'closed', 'TidalEntrance', 'ClosedInner'}) dictionary of label names for the physical lines.
        """

        # set local variables based on global variables
        boundaries = self.boundaries
        n_Spline = self.n_Spline
        if self.inner:
            Inner_Splines_ids = self.Inner_Splines_ids
        physical_lines = self.physical_lines

        # add a header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Boundaries\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')
        # if there are no boundaries
        if len(boundaries)==0:
            f.write('Physical Line("%s") = {' % physical_line_labels[0])
            for i in range(n_Spline-2):
                f.write(' %d, ' % (i+1))
            f.write('%d};\n' % (n_Spline-1))
        # if there are boundaries
        else:
            # initialize a list to store all id's of what should be closed lines
            closed_lines = []
            # loop over all
            for i in range(1, n_Spline):
                # if it is not a physical line
                if not physical_lines.__contains__(i):
                    closed_lines.append(i)
        # add all boundary lines to the physical line
        if len(boundaries) > 0:
            j = 1
            f.write('Physical Line("%s") = {1' % physical_line_labels[j])
            for i in range(1,len(physical_lines)):
                if physical_lines[i] - physical_lines[i-1] == 1:
                    # add each pysical line
                    f.write(', %d' % (physical_lines[i]))
                else:
                    f.write('};\n')
                    j += 1
                    f.write('Physical Line("%s") = {%d' % (physical_line_labels[j], physical_lines[i]))


            # close the physical line object
            f.write('};\n')

            # non-boundary lines are 'closed' boundaries
            f.write('Physical Line("%s") = {' % physical_line_labels[0])
            for i in range(len(closed_lines)-1):
                f.write(' %d,' % (closed_lines[i]))
            f.write(' %d};\n' % (closed_lines[-1]))

        if self.inner_rings:
            # all inner holes are ClosedInner boundaries
            f.write('Physical Line("%s") = {' % physical_line_labels[-1])
            for i in range(len(Inner_Splines_ids)-1):
                f.write(' %d,' % (Inner_Splines_ids[i]))
            f.write(' %d};\n' % (Inner_Splines_ids[-1]))

        return f

    def writeGeoFileAPrioriResolution(self, f):
        """
        Method to set a a priori resolution

        Args:
            f: (Required) open text file representing the .geo file
            autocoherence: (Optional, defaults to False) True/False boolean to indicate autocoherence
        """

        # add header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ A Priori Resolution\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')
        # Accuracy of evaluation of the LC field for 1D mesh generation
        # Prohibit boundary mesh size interpolation (suggested if the cell size is indicated by a field)
        f.write('Mesh.LcIntegrationPrecision = 1e-3;\nMesh.CharacteristicLengthExtendFromBoundary = 0;\n')

        return f

    def writeGeoFileLineLoopsSurface(self, f, inner_rings = True):
        """
        Method to write the lines indicating the line loop and surface of the exterior ring
        Args:
            f: (Required) open text file representing the .geo file
            inner_rings: (Optional, defaults to True) True/False boolean to add the points and splines of the inner rings
        """

        # set global parameters as local
        n_LineLoop = self.n_LineLoop
        n_Surface = self.n_Surface
        n_Point = self.n_Point
        n_Spline = self.n_Spline
        if inner_rings:
            Inner_Splines_ids_perLL = self.Inner_Splines_ids_perLL

        # create a header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Create Line Loop and Surface\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')

        # create a line loop with all splines of the exterior
        # open the line loop
        f.write('Line Loop(%d) = {' % (n_LineLoop))
        n_LineLoop += 1
        # add all the splines
        for i in range(1,n_Spline-1):
            f.write('%d, '% (i))
        # close the line loop
        f.write('%d};\n' % (n_Spline-1))

        # check if there are inner rings
        if inner_rings:
            # create header
            f.write('//+ Inner loops\n')
            # counter
            t = 0
            # go over all inner rings ( Inner_Points is a list of list and so, its length is equal to the number of rings)
            for I in self.Inner_Points:
                # per ring, write a line loop, plane and physical surface
                f.write('Line Loop(%d) = {%d' % (n_LineLoop, Inner_Splines_ids_perLL[t][0]))
                if len(Inner_Splines_ids_perLL[t]) > 1:
                    for i in range(1, len(Inner_Splines_ids_perLL[t])):
                        f.write(', %d' % Inner_Splines_ids_perLL[t][i])
                f.write('};\n')

                # count!
                t += 1
                n_LineLoop += 1

        # create header
        f.write('//+ Plane Surface\n')
        # create a plane surface from the lineloop
        f.write('Plane Surface(%d) = {1' % (n_Surface))

        if inner_rings:
            for i in range(2,n_LineLoop):
                f.write(', %d' % i)
        f.write('};\n');

        # create a physical surface from the plane surface (no water can pass this surface)
        f.write('Physical Surface("inside") = {%d};\n' % (n_Surface))

        # count!
        n_LineLoop += 1
        n_Surface += 1

        # globalize local parameters
        self.n_LineLoop = n_LineLoop
        self.n_Surface = n_Surface
        self.n_Point = n_Point
        self.n_Spline = n_Spline

        return f

    def writeGeoFileHomogeneousMesh(self, f, mesh_size = 250):
        """
        Method to set a homogeneous cell size in the mesh

        Args:
            f: (Required) open text file representing the .geo file
            mesh_size: (Optional, defaults to 250) float representing the desired mesh size

        """

        # add header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Set Homogeneous Mesh\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')
        # create a field with the set homogeneous cell size everywhere
        f.write('//+\nField[1] = MathEval;\nField[1].F = "%f";\n' % (mesh_size))
        # set this field as the background field
        f.write('Background Field = 1;')

        return f

    def createBackgroundFile(self, XY_field, multiplier = 1, minres = 25, maxres = 1e9, zerovalues = 250, plot_arr = './test.png', buffer_interpolation = 50 ):
        """
        Method to create a background file

        Args:
            XY_field: (Required) directory path string or list of strings to a pickle which hold a list with the following elements:
                        - 2d numpy array representing a spatial raster with values for the mesh size (!!! non-channels cells should be zero)
                        - coordinate pair with the coordinates of the left bottom corner of the 2d numpy array
                        - x resolution
                        - y resolution
            multiplier: (Optional, defaults 1) value to multiply the original array with
            minres: (Optional, defaults to 25) minimum mesh size
            zerovalues: (Optional, defaults to 250) the mesh size in areas with value 0
            plot_arr: (Optional, defaults to False) False/file directory to store a plot of the original raster
            buffer_interpolation: (Optional, defaults to False) False/float to indicate the size of buffer around the channels to interpolate

        """
        print('Creating a background file...')
        # open the pickled file with the needed information for a heterogenuous background field
        print(XY_field)

        with open(XY_field, 'rb') as input:
            xy = pickle.load(input)
        # background file
        background_file = XY_field[:-3] + 'bg'
        # function to create a valid background field from a Numpy 2D array
        np2BackgroundField(xy[0], background_file , multiplier = multiplier, minres = minres, maxres = maxres,
            zerovalues = zerovalues, zeropoint= xy[1], xres = xy[2], yres = xy[3], plot_arr = plot_arr,
            buffer_interpolation = buffer_interpolation)

        print('Background file created!')
        print('-----------------------------------')
        return background_file

    def writeGeoFileHeterogeneousMesh(self, f, XY_field, multiplier = 1, minres = 25, maxres = 1e9,
        zerovalues = 250, outside_mesh_size = 250, plot_arr = './backgroundfield.png',
        buffer_interpolation = 50, create_background_file = True, background_file = ''):
        """
        Method to set a heterogenuous cell size in the mesh

        Args:
            f: (Required) open text file representing the .geo file
            XY_field: (Required)directory path string to a pickle which hold a list with the following elements:
                        - 2d numpy array representing a spatial raster with values for the mesh size
                        - coordinate pair with the coordinates of the left bottom corner of the 2d numpy array
                        - x resolution
                        - y resolution
            multiplier: (Optional, defaults 1) value to multiply the original array with
            minres: (Optional, defaults to 25) minimum mesh size
            zerovalues: (Optional, defaults to 250) the mesh size in areas with value 0
            plot_arr: (Optional, defaults to False) False/file directory to store a plot of the original raster
            buffer_interpolation: (Optional, defaults to False) False/float to indicate the size of buffer around the channels to interpolate
            create_background_file: (Optional, defaults to True) True/False boolean which indicates whether to create a background file

        """
        if create_background_file:
            if type(XY_field) != list:
                bf = [self.createBackgroundFile(XY_field, multiplier = multiplier, minres = minres, maxres = maxres, zerovalues = zerovalues, plot_arr = plot_arr, buffer_interpolation = buffer_interpolation)]
            else:
                bf = []
                for field in XY_field:
                    bf.append(self.createBackgroundFile(field, multiplier = multiplier, minres = minres,maxres = maxres, zerovalues = zerovalues, plot_arr = plot_arr, buffer_interpolation = buffer_interpolation))
        else:
            if len(background_file) == 0:
                bf = [XY_field.split('.')[0]+'.bg']
            else:
                if type(background_file) == list:
                    bf = background_file
                else: bf = [background_file]

        # add header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Set Heterogenuous Mesh\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')

        # complete the background directory if it's a local directory
        bf_correctpaths = []
        for background_file in bf:
            if background_file[0] == '.':
                background_file = os.getcwd() + background_file[1:]
            bf_correctpaths.append(background_file)

        # make sure that the required parameters are in a list format
        if type(outside_mesh_size) != list: outside_mesh_size = [outside_mesh_size]

        n = 0
        for background_file in bf_correctpaths:
            n += 1
            # add this background field as the structured background field
            f.write('Field[%d] = Structured;\nField[%d].FileName = "%s";\n' % (n, n, background_file))
            # set the cell size outside the background field
            f.write('Field[%d].OutsideValue = %f;\n' % (n, outside_mesh_size[n-1]))
            f.write('Field[%d].SetOutsideValue = %d;\n' % (n,n))
            # set the format of this background field and set is at the background file to apply the meshing on
            f.write('Field[%d].TextFormat = 1;\n' % (n))

        if n == 1:
            f.write('Background Field = 1;\n//+')
        else:
            n += 1
            f.write('Field[%d] = Min;\n' % n)
            f.write('Field[%d].FieldsList = {1' % n)
            for i in range(2,n):
                f.write(', %d' % i)
            f.write('};\n')
            f.write('Background Field = %d;\n//+' % n)

        return f

    def writeGeoFileDefineAlgorithm(self, f, algorithm = 'Frontal-Delaunay'):
        """
        Method to define the 2D meshing algorithm (https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options-list)

        Args:
            f: (Required) open text file representing the .geo file
            algoritm: (Optional, defaults to 'Frontal-Delaunay') string with the 2D algoritm:
                - Frontal-Delaunay
                - Delaunay

        """
        # add header
        f.write('//+\n//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n'
                '//+ Meshing Algorithm\n//++++++++++++++++++++++++++++++++++++++++++++++++++++'
                '++++++++++++++++++\n')

        if algorithm == 'Frontal-Delaunay': algo = 6
        elif algorithm == 'Delaunay': algo = 5
        else: algo = 6

        f.write("Mesh.Algorithm = %d;\n" % algo)

        return f

    def meshTheDomain(self, outputfile_mesh, geo_file = 'self', logfile = False, ncores = 6):
        """
        Method to mesh the .geo file

        Args:
            outputfile_mesh: (Required) directory path string to save the .mesh file
            geo_file: (Optional, defaults to the geofile of the object) directory path string of the .geo file

        """
        if geo_file == 'self':
            geo_file = self.geofile


        self.meshfile = outputfile_mesh

        if logfile:
            os.system('gmsh ' + geo_file + ' -nt %d -2 -o ' % ncores + outputfile_mesh + ' -format msh2 -log '+ logfile)
        else:
            # mesh the .geo file
            os.system('gmsh ' + geo_file + ' -nt %d -2 -o ' % ncores + outputfile_mesh + ' -format msh2')

    def getCoordinatesInnerHoleNodes(self, meshfile):
        """
        Method to retreive the coordinate pairs of all nodes which are located on the edge of an inner hole in a mesh (generated by gmsh)

        Args:
            meshfile: (Required) directory path string to the mesh file

        Returns:
            List of Numpy arrays of dimensions n x 2 with the coordinate pairs of all nodes which are located on the edge of an inner hole

        """

        # load mesh file
        mesh = open(meshfile)
        # convert to list of strings
        mesh_lines = mesh.readlines()

        # loop over all lines
        for i in range(len(mesh_lines)):
            line = mesh_lines[i]

            # save the line where the node sections starts and ends
            if line.startswith('$Nodes'):
                line_start_nodes = i + 2
            if line.startswith('$EndNodes'):
                line_end_nodes = i

            # save the line where the elements sections starts and ends
            if line.startswith('$Elements'):
                line_start_elements = i + 2
            if line.startswith('$EndElements'):
                line_end_elements = i

        # initialize two numpy arrays to store the codes on nodes and elements
        nodes = np.zeros([line_end_nodes - line_start_nodes, 4])
        elements = np.zeros([line_end_elements - line_start_elements, 7])

        # initiate counters
        t_nodes = 0
        t_elements = 0

        # again, loop over all lines
        for i in range(len(mesh_lines)):
            line = mesh_lines[i]

            # if it is a line within the nodes section
            if line_start_nodes <= i < line_end_nodes:
                # split the line into the different numbers and save as numpy array
                row = np.float64(line.split(" "))
                # save this line node codes in the numpy array
                nodes[t_nodes, :] = row
                t_nodes += 1

            if line_start_elements <= i < line_end_elements:
                # split the line into the different numbers and save as numpy array
                row = np.int64(line.split(" "))
                # only consider 1D line elements (these are located on the edge of the domain, and thus on the edge of inner holes)
                if int(row[1]) == 1:
                    # save this line node codes in the numpy array
                    elements[t_elements, :] = row
                    t_elements += 1

        # remove the zero-rows in the elements array. These rows would have been rows describing elements in the interior of the domain.
        # There is no need to save them and this would slow down the program, hence we didn't fill these rows and thus, they can be deleted.
        elements = elements[~np.all(elements == 0, axis=1)]
        # get all mesh partition id's
        mesh_parts = np.unique(elements[:,-3])
        # get the id's of the mesh partitions which belong to an inner hole
        inner_rings = mesh_parts[-self.n_inner_rings:]

        # initialize list to store the Numpy array of coordinate pairs per inner ring
        Inner_ring_nodes_xy = []

        # loop over all id's of the inner rings
        for id in inner_rings:
            # search for the rows in elements which are part of this inner ring
            r = np.where(elements[:,-3] == id)
            # pick the nodes of these elements
            nodes1 = np.int32(elements[r,-1])
            nodes2 = np.int32(elements[r,-2])
            ring_nodes = np.concatenate((nodes1, nodes2))
            ring_nodes = np.unique(ring_nodes)
            # look for all the coordinates of these nodes and append them to the list of coordinate pairs per inner ring
            ring_nodes_xy = nodes[np.isin(nodes[:,0], ring_nodes), -3:-1]
            Inner_ring_nodes_xy.append(ring_nodes_xy)

        return Inner_ring_nodes_xy

    def meshInnerHoles(self, outputfile_mesh, Inner_ring_nodes_xy):
        """
        Method to mesh all the inner holes and save the .msh files

        Args:
            Inner_ring_nodes_xy: (Required) List of n x 2 Numpy arrays with the coordinate pairs of the nodes on the edges of the inner holes
            original_mesh: (Required) Directory path string of the original mesh with holes
            outputfile_mesh: (Required) Directory path string to store the newly generated mesh with the holes included
        """
        # save the original exterior list
        exter_original = self.exter
        # path to original mesh
        original_mesh = self.meshfile
        fn = original_mesh[:-4]

        # loop over all inner rings
        for i in range(len(Inner_ring_nodes_xy)):
                xy = Inner_ring_nodes_xy[i]
                self.exter = xy
                # set the boundaries and breaks of each ring at 0
                self.setBoundariesAndBreaks(boundaries = [], breaks = [])
                # write a .geo file per inner ring, use the same field for cell size as the full domain, use lines instead of splines
                # and prevent the adding of new nodes to the edge, do not add new physical lines
                self.writeGeoFile(fn + '_hole_%d.geo' % (i+1), cell_size = 'hetero', inner_rings = False,
                                  algo = self.algo, line_type = 'Lines', prevent_add_nodes_to_edge = True,
                                  create_background_file = False, add_physical_lines = False)
                # mesh this geometry file
                self.meshTheDomain(fn + '_hole_%d.msh' % (i+1),fn + '_hole_%d.geo' % (i+1))

        # write a gsmh script to merge the meshes
        # !!! removing excessive physical lines messes up the mesh
        fm = open('delete_physicals.gmsh', 'w+')
        fm.write('Merge "%s";\n' % (original_mesh[:-4]+'_hole_1.msh'))
        fm.write('Merge "%s";\n' % (original_mesh[:-4]+'_hole_2.msh'))
        fm.write('Merge "%s";\n' % original_mesh)
        fm.write('Coherence Mesh;\n')
        fm.close()

        #os.system('gmsh delete_physicals.gmsh')
        os.system('gmsh delete_physicals.gmsh -2 -o ' + outputfile_mesh + ' -format msh2')

    def meshOneInnerHole(self, i, fn):
        """
        Method to mesh an inner hole

        Args:
            i: (Required) Integer representing which inner hole to mesh
            fn: (Required) Directory path string to store the .geo and .msh file (do not add extension)
        """

        exter = self.exter
        breaks = self.breaks
        boundaries = self.boundaries
        inner = self.inter[i]

        try: geofile = self.geofile
        except: pass

        self.exter = inner
        self.breaks = []
        self.boundaries = []

        self.writeGeoFile(fn+'.geo', cell_size = 50, inner_rings = False, add_physical_lines = False)
        self.meshTheDomain(fn+'.msh', geo_file = fn+'.geo')

        self.boundaries = boundaries
        self.breaks = breaks
        self.exter = exter
        try: self.geofile = geofile
        except: pass

    @staticmethod
    def nodes2CSV(fn, fn_mesh):
        """
        Method to translate all nodes in a gsmh file to a csv with two columns X & Y coordinates

        Args:
            fn: (Required) path string to save the csv file
            fn_mesh: (Required) path string to load the mesh from
        """
        with open(fn, 'w+') as fcsv:
            fcsv.write('X,Y\n')
            with open(fn_mesh, 'r') as fmesh:
                lines = fmesh.readlines()
                nodes_n = 0
                within_nodes = False
                for line in lines:
                    if line.startswith('$EndNodes'):
                        within_nodes = False

                    if within_nodes:
                        if nodes_n == 0:
                            nodes_n = int(line)
                            print('%d nodes' % nodes_n)
                        else:
                            line = line.split(' ')
                            x = float(line[1])
                            y = float(line[2])
                            fcsv.write('%f,%f\n' % (x,y))
                            line = line[0]

                    if line.startswith('$Nodes'):
                        within_nodes = True

    def calcEquilaterality(self, fn):
        """
        Method to calculate a csv with coordinate pairs and Equilaterality values for each triangle

        Args:
        - fn: (Required) Directory path string to save the csv
        """

        with open(self.meshfile, 'r') as f:

                lines = f.readlines()

                nodes = False
                els = False

                Nodes = []
                Triangles = []

                for line in lines:
                    if line.startswith('$Nodes'):
                        nodes = True

                    elif line.startswith('$EndNodes'):
                        nodes = False

                    elif nodes:
                        if len(line.split(' ')) > 1:
                            x = float(line.split(' ')[1])
                            y = float(line.split(' ')[2])
                            Nodes.append([x,y])

                    elif line.startswith('$Elements'):
                        els = True

                    elif line.startswith('$EndElements'):
                        els = False

                    elif els:
                        if len(line.split(' ')) > 1:
                            if '2 2 3 1' in line:
                                t1, t2, t3 = [int(i) for i in line.split(' ')[-3:]]
                                Triangles.append([t1-1, t2-1, t3-1])

                Nodes = np.asarray(Nodes)
                Equilateralities = []
                Triangles_mp = []

                for triangle in Triangles:
                    p1, p2, p3 = Nodes[triangle,:]
                    ptriangle = np.array([p1,p2,p3])
                    mp = np.mean(ptriangle, axis = 0)
                    Triangles_mp.append(mp)


                    d12 = ((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2)**0.5
                    d13 = ((p1[1]-p3[1])**2 + (p1[0]-p3[0])**2)**0.5
                    d23 = ((p2[1]-p3[1])**2 + (p2[0]-p3[0])**2)**0.5

                    h1 = np.arccos((d23**2 - d12**2 - d13**2)/(-2*d12*d13))
                    h2 = np.arccos((d13**2 - d12**2 - d23**2)/(-2*d12*d23))
                    h3 = np.arccos((d12**2 - d13**2 - d23**2)/(-2*d13*d23))
                    angles = np.array([h1, h2, h3])
                    equi = np.max(angles)/np.min(angles)
                    Equilateralities.append(equi)

                Equilateralities = np.asarray(Equilateralities)
                Triangles_mp = np.asarray(Triangles_mp)

                df = pd.DataFrame({'X': Triangles_mp[:,0], 'Y': Triangles_mp[:,1], 'Equi': Equilateralities})
                df.to_csv(fn, index = False)

    @staticmethod
    def checkForDuplicatesInGeoFile(geofile_fn):

        Points = []
        Point_IDs = []
        with open(geofile_fn) as f:
            lines = f.readlines()
            for line in lines:
                if line.startswith('Point'):
                    id = int(line.split(' = ')[0].replace('Point(', '')[:-1])
                    coor = line.split('=')[1].split(',')
                    x = float(coor[0].replace('{', ''));
                    y = float(coor[1])
                    Point_IDs.append(id)
                    Points.append([x, y])

        Point_IDs = np.asarray(Point_IDs)
        Points = np.asarray(Points)
        values, ids, counts = np.unique(Points, axis=0, return_counts=True, return_index=True)
        duplicates = values[counts > 1, :]
        duplicate_ids = Point_IDs[ids[counts > 1]]

        for d,id in enumerate(duplicate_ids):
            print('Duplicate detected for Point(%d):' % id)
            print('    X: %.2f' % duplicates[d, 0])
            print('    Y: %.2f' % duplicates[d, 1])
            print('-----------------------------')

        if len(duplicate_ids) == 0:
            print('No duplicates found.')

        return len(duplicate_ids)==0

    """
    @staticmethod
    def calcEquilaterality_static(mesh_fn, fn):
        
        Method to calculate a csv with coordinate pairs and Equilaterality values for each triangle

        Args:
        - fn: (Required) Directory path string to save the csv
        

        with open(mesh_fn, 'r') as f:

                lines = f.readlines()

                nodes = False
                els = False

                Nodes = []
                Triangles = []

                for line in lines:
                    if line.startswith('$Nodes'):
                        nodes = True

                    elif line.startswith('$EndNodes'):
                        nodes = False

                    elif nodes:
                        if len(line.split(' ')) > 1:
                            x = float(line.split(' ')[1])
                            y = float(line.split(' ')[2])
                            Nodes.append([x,y])

                    elif line.startswith('$Elements'):
                        els = True

                    elif line.startswith('$EndElements'):
                        els = False

                    elif els:
                        if len(line.split(' ')) > 1:
                            if line.split(' ')[1] ==  '2':
                                t1, t2, t3 = [int(i) for i in line.split(' ')[-3:]]
                                Triangles.append([t1-1, t2-1, t3-1])

                Nodes = np.asarray(Nodes)
                print(len(Triangles))
                Equilateralities = []
                Triangles_mp = []

                t = 0
                print('Calculating Triangle angles:')
                for triangle in Triangles:

                    print('   %.2f %%' % (t/len(Triangles)*100), end = '\r')
                    t += 1

                    p1, p2, p3 = Nodes[triangle,:]
                    ptriangle = np.array([p1,p2,p3])
                    mp = np.mean(ptriangle, axis = 0)
                    Triangles_mp.append(mp)

                    d12 = ((p1[1]-p2[1])**2 + (p1[0]-p2[0])**2)**0.5
                    d13 = ((p1[1]-p3[1])**2 + (p1[0]-p3[0])**2)**0.5
                    d23 = ((p2[1]-p3[1])**2 + (p2[0]-p3[0])**2)**0.5

                    h1 = np.arccos((d23**2 - d12**2 - d13**2)/(-2*d12*d13))
                    h2 = np.arccos((d13**2 - d12**2 - d23**2)/(-2*d12*d23))
                    h3 = np.arccos((d12**2 - d13**2 - d23**2)/(-2*d13*d23))
                    angles = np.array([h1, h2, h3])
                    equi = np.max(angles)/np.min(angles)
                    Equilateralities.append(equi)

                Equilateralities = np.asarray(Equilateralities)
                Triangles_mp = np.asarray(Triangles_mp)

                df = pd.DataFrame({'X': Triangles_mp[:,0], 'Y': Triangles_mp[:,1], 'Equi': Equilateralities})
                df.to_csv(fn, index = False)"""