Parse a GMSH .msh file and return a dictionary containing the data
    neccessary to create CRTomo grids

def parse_gmsh(filename, boundary_file):
    """
    Parse a GMSH .msh file and return a dictionary containing the data
    neccessary to create CRTomo grids
    """
    mesh = {}

    fid = open(filename, 'r')
    line = fid.readline()
    while(line):
        if(line.startswith('$MeshFormat')):
            pass
        elif(line.startswith('$Nodes')):
            nodes = []
            line = fid.readline()
            nr_nodes = np.fromstring(line, dtype=int, count=1, sep=r'\n')
            nr_nodes
            while(line):
                line = fid.readline()
                if(line.startswith('$EndNodes')):
                    break
                node = np.fromstring(line, dtype=float, sep=' ')
                nodes.append(node)
            mesh['nodes'] = nodes
        elif(line.startswith('$Elements')):
            """
            Create a dictionary with the element types as keys. E.g.:
            elements['15'] provides all elements of type 15 (Points)
            """
            elements = {}
            line = fid.readline()
            nr_elements = np.fromstring(line, dtype=int, count=1, sep=r'\n')
            nr_elements
            while(line):
                line = fid.readline()
                if(line.startswith('$EndElements')):
                    break
                element = np.fromstring(line, dtype=int, sep=' ')
                # el_nr = element[0]
                el_type = element[1]
                el_nr_tags = element[2]
                # el_tags = element[3:3 + el_nr_tags]
                el_nodes = element[3 + el_nr_tags:]

                # now decide where to put it
                key = str(el_type)
                if(key in elements.keys()):
                    elements[key].append(el_nodes)
                else:
                    elements[key] = []
                    elements[key].append(el_nodes)

            mesh['elements'] = elements
        line = fid.readline()

    fid.close()

    # if boundary_file is != None, then sort the lines (element type 1)
    # according to the element types
    boundaries = {}

    if(boundary_file is not None):
        # load the original boundary lines
        # it is possible that GMSH added additional nodes on these lines, and
        # that is why we need to find all mesh lines that lie on these original
        # lines.
        bids = np.loadtxt(boundary_file)

        for btype in ('12', '11'):
            # select all original boundaries with this type
            a = np.where(bids[:, 4] == int(btype))[0]
            boundaries[btype] = []
            # for each of those lines, find all lines of the mesh that belong
            # here
            for orig_line in bids[a, :]:
                # print('Find all lines lying on the line: ')
                found_one_line = False
                # print(orig_line)
                # construct line equation

                # x1 == x2 ?
                # split into coordinates
                ox1 = orig_line[0]
                ox2 = orig_line[2]
                oy1 = orig_line[1]
                oy2 = orig_line[3]

                if(orig_line[0] == orig_line[2]):
                    # special case: we only need to find all lines with x1 ==
                    # x2 == x1_orig and y_min >= y_orig_min and y_max <=
                    # <_orig_max
                    for line in elements['1']:
                        if(btype == '11'):
                            if(line[0] == 48 and line[1] == 150):
                                pass
                                # print('Find all lines lying on the line: ')
                                # print('This is the line')

                        # it doesn't matter any more to be able to assign x ->
                        # y values. Thus we can sort the y values and just
                        # check
                        # if the new line lies in between the original one
                        oy1, oy2 = np.sort([orig_line[1], orig_line[3]])
                        x1, x2 = np.sort(
                            [
                                mesh['nodes'][line[0] - 1][1],
                                mesh['nodes'][line[1] - 1][1]
                            ]
                        )
                        y1, y2 = np.sort(
                            [
                                mesh['nodes'][line[0] - 1][2],
                                mesh['nodes'][line[1] - 1][2]
                            ]
                        )

                        if np.isclose(x1, x2) and np.isclose(x2, ox1):
                            if(y1 >= oy1 and y2 <= oy2):
                                found_one_line = True
                                boundaries[btype].append(line)

                else:
                    # print('checking with full line equation')
                    # no vertical line
                    # we need the full check using the line equation
                    slope = (orig_line[1] - orig_line[3]) / (
                        orig_line[0] - orig_line[2])
                    intersect = orig_line[1] - (slope * orig_line[0])
                    # print('Slope', slope, ' Intercept ', intersect)
                    for line in elements['1']:
                        x1 = mesh['nodes'][line[0] - 1][1]
                        y1 = mesh['nodes'][line[0] - 1][2]
                        x2 = mesh['nodes'][line[1] - 1][1]
                        y2 = mesh['nodes'][line[1] - 1][2]

                        # print(x1, x2, y1, y1)
                        check = False
                        # check if x coordinates of the test line fit in the
                        # original line
                        if(ox1 < ox2):
                            if(x1 < x2):
                                if((np.isclose(x1, ox1) or x1 > ox1) and
                                   (np.isclose(x2, ox2) or x2 < ox2)):
                                    check = True
                            else:
                                if((np.isclose(x2, ox1) or x2 >= ox1) and
                                   (np.isclose(x1, ox2) or x1 <= ox2)):
                                    check = True
                        else:
                            if(x1 < x2):
                                if((np.isclose(x1, ox2) or x1 >= ox2) and
                                   (np.isclose(x2, ox1) or x2 <= ox1)):
                                    check = True
                            else:
                                if((np.isclose(x2, ox2) or x2 >= ox2) and
                                   (np.isclose(x1, ox1) or x1 <= ox1)):
                                    check = True

                        # print('boundary check:', check)
                        if(check):
                            # the line lies within the x-range of the orig line
                            ytest1 = slope * x1 + intersect
                            ytest2 = slope * x2 + intersect
                            if(np.around(ytest1 - y1, 5) == 0 and
                               np.around(ytest2 - y2, 5) == 0):
                                boundaries[btype].append(line)
                                # found = True
                                found_one_line = True
                                # print('found it new', line)
                # add a weak check: we need to find at least one line in the
                # mesh corresponding to this boundary line:
                if not found_one_line:
                    raise Exception('no mesh line found for this boundary')

            print('Total number of boundaries of this type:',
                  len(boundaries[btype]))
    mesh['boundaries'] = boundaries
    return mesh