def read_inp(path):
  """
  Reads Abaqus inp file
  """
  
  def lineInfo(line):
    out =  {"type": "data"}
    if line[0] == "*":
      if line[1] == "*": 
        out["type"] = "comment"
        out["text"] = line[2:]
      else:
        out["type"] = "command"
        words = line[1:].split(",")
        out["value"] = words[0].strip()
        out["options"] = {}
        for word in words[1:]:
          key, value =  [s.strip() for s in word.split("=")]
          out["options"][key] = value
    return out
  
  def elementMapper(inpeltype):
    if inpeltype == "t3d2": return "Line2"
    if inpeltype[:3] in ["cps", "cpe", "cax"]:
      if inpeltype[3] == "3": return "tri3"
      if inpeltype[3] == "4": return "quad4"
    if inpeltype[:3] in ["c3d"]:
      if inpeltype[3] == "4": return "tetra4"
      if inpeltype[3] == "5": return "pyra5"
      if inpeltype[3] == "6": return "prism6"
      if inpeltype[3] == "8": return "hexa8"
    
  nlabels      = []
  coords       = []
  nsets        = {}
  elabels      = []
  etypes       = []
  connectivity = []
  esets        = {}
  surfaces     = {}
  
  # File preprocessing
  lines = np.array([l.strip().lower() for l in open(path).readlines()])
  lines = [line for line in  lines if len(line) != 0]
  # Data processing
  env, setlabel = None, None
  for line in lines: 
    d = lineInfo(line)
    if d["type"] == "command": 
      env = d["value"]
      # Nodes
      if env == "node":
        opt = d["options"]
        currentset = None
        if "nset" in opt.keys(): 
          currentset = opt["nset"]
          nsets[currentset] = []
           
      # Elements
      if env == "element":
        opt = d["options"]
        eltype = elementMapper(opt["type"])
        currentset = None
        if "elset" in opt.keys(): 
          currentset = opt["elset"]
          esets[currentset] = []
          
      # Nsets
      if env == "nset":
        opt = d["options"]      
        currentset = opt["nset"]
        nsets[currentset] = []
        
      # Elsets     
      if env == "elset":
        opt = d["options"]      
        currentset = opt["elset"]
        esets[currentset] = []
      
      # Surfaces
      if env == "surface":
        opt = d["options"]
        currentsurface = opt["name"]
        if opt["type"] == "element":
          surfaces[currentsurface] = []  
                    
             
    if d["type"] == "data": 
      words = line.strip().split(",")
      if env == "node":
        label = int(words[0]) 
        nlabels.append(label) 
        coords.append(
            np.array([np.float64(w) for w in words[1:4]])
                     )
        if currentset != None: nsets[currentset].append(label)
            
              
      if env == "element": 
        label  = int(words[0])
        elabels.append(label)
        connectivity.append(
            np.array( [np.int32(w) for w in words[1:] if len(w) != 0 ])
                           )
        etypes.append(eltype)                   
        if currentset != None: esets[currentset].append(label)
      
      if env == "nset":
        nsets[currentset] += [int(w) for w in words if len(w) != 0]   
        
      if env == "elset":
        esets[currentset] += [int(w) for w in words if len(w) != 0]  
      
      if env == "surface":
        if opt["type"] == "element":
          surfaces[currentsurface].append([w.strip() for w in words])
  
  surfaces2 = {}        
  for tag, surface in surfaces.items():
    surfaces2[tag] = []
    for sdata in surface:
      labels = esets[sdata[0]]
      face = int(sdata[1].split("s")[1].strip())-1
      for label in labels:
        surfaces2[tag].append((label, face))             
  
  return Mesh(nlabels = nlabels,
              coords  = coords,
              nsets   = nsets,
              elabels = elabels,
              etypes  = etypes,
              connectivity = connectivity,
              esets = esets,)
Reads Abaqus inp file