docs: 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