msh2inp.awk

Go to the documentation of this file.

#!/usr/bin/awk -We

# =======
# msh2inp
# =======
#
# :Version: 1.03
# :Date: 2009-04-14
# :Purpose: Convert gmsh msh files to UCD/inp format.
# :Language: awk, tested with gawk, mawk, nawk.
# :Author: Stefan Tibus
# :Copyright: (C) 2009 Stefan Tibus
# :License: GNU General Public License
#
# Changelog
# =========
#
# :2009-04-07 v1.00: First version.
# :2009-04-07 v1.01: Renumber used cells and nodes by default, may be
#                    disabled. *mawk* compatibility broken due to use of
#                    ``asorti``.
#                    :ToDo: Sort order is alphabetical, thus not only
#                           the numbering, but actual order is changed,
#                           even when writing out all elements.
# :2009-04-07 v1.02: Creating the mapping in the scanning loops makes
#                    ``asorti`` unnecessary, results in numerical order
#                    and is compatible with *mawk* again.
# :2009-04-14 v1.03: Added usage information.
#
#
# Background Information
# ======================
#
# INP (UCD) file format:
#
# - http://people.sc.fsu.edu/~burkardt/html/ucd_format.html
#
# gmsh msh file format:
#
# - http://www.geuz.org/gmsh/doc/texinfo/gmsh_10.html#SEC55
# - http://www.geuz.org/gmsh/doc/texinfo/gmsh-full.html#SEC55

function print_usage() {
  print ""
  print "msh2inp v1.03"
  print "© Stefan Tibus"
  print ""
  print "Convert gmsh msh files to UCD/inp format."
  print ""
  print "Usage:"
  print "msh2inp [options] <file> [<file> [...]]"
  print "cat <file> | msh2inp [options]"
  print ""
  print "Options:"
  print "--help        Show usage information."
  print "--debug       Provide some informational output on /dev/stderr."
  print "--renumber    Renumber nodes and elements in output (default)."
  print "--norenumber  Keep numbering as in input files."
  print ""
  print "Output filtering options:"
  print "(none)        Output all element types as in input files."
  print "--pt          Points."
  print "--line        Lines."
  print "--tri         Triangles."
  print "--quad        Quadrangles."
  print "--tet         Tetrahedra."
  print "--hex         Hexahedra."
  print "--prism       Prisms."
  print "--pyr         Pyramids."
  print ""
  print "Note: gmsh saves physical mesh elements only by default. Remember to assign\n      physical element numbers to all required element types."
  print ""
}

BEGIN {
  # Simple option handling.
  inp_filter_output = 0
  DEBUG = 0
  inp_renumber = 1
  for (i = 1; i < ARGC; i++) {
    if (ARGV[i] == "--renumber") {
      inp_renumber = 1
    } else
    if (ARGV[i] == "--norenumber") {
      inp_renumber = 0
    } else
    if (ARGV[i] == "--debug") {
      DEBUG = 1
    } else
    if (ARGV[i] == "--line") {
      inp_filter_output = 1
      inp_output[1] = 1
    } else
    if (ARGV[i] == "--tri") {
      inp_filter_output = 1
      inp_output[2] = 1
    } else
    if (ARGV[i] == "--quad") {
      inp_filter_output = 1
      inp_output[3] = 1
    } else
    if (ARGV[i] == "--tet") {
      inp_filter_output = 1
      inp_output[4] = 1
    } else
    if (ARGV[i] == "--hex") {
      inp_filter_output = 1
      inp_output[5] = 1
    } else
    if (ARGV[i] == "--prism") {
      inp_filter_output = 1
      inp_output[6] = 1
    } else
    if (ARGV[i] == "--pyr") {
      inp_filter_output = 1
      inp_output[7] = 1
    } else
    if (ARGV[i] == "--pt") {
      inp_filter_output = 1
      inp_output[15] = 1
    } else
    if (ARGV[i] == "--help") {
      print_usage()
      error = -127
      exit error
    } else
    if (ARGV[i] ~ /^--./) {
      print "Unknown option \"" ARGV[i] "\"." > "/dev/stderr"
      print_usage()
      error = 1
      exit error
    } else
    # Stop option handling on first non-option.
    break
    # Delete handled option from argument list.
    delete ARGV[i]
  }

  msh_read_part = 0

  msh_part_meshformat = 1
  msh_part_nodes = 2
  msh_part_elements = 3
  msh_part_physicalnames = 4
#  msh_part_nodedata = 5
#  msh_part_elementdata = 6
#  msh_part_elementnodedata = 7
  msh_part_unknown = -1

  # Element type conversions (first order elements only).
  msh_elementtype_inp[1] = "line"
  msh_elementtype_nodes[1] = 2
  msh_elementtype_inp[2] = "tri"
  msh_elementtype_nodes[2] = 3
  msh_elementtype_inp[3] = "quad"
  msh_elementtype_nodes[3] = 4
  msh_elementtype_inp[4] = "tet"
  msh_elementtype_nodes[4] = 4
  msh_elementtype_inp[5] = "hex"
  msh_elementtype_nodes[5] = 8
  msh_elementtype_inp[6] = "prism"
  msh_elementtype_nodes[6] = 6
  msh_elementtype_inp[7] = "pyr"
  msh_elementtype_nodes[7] = 5
  msh_elementtype_inp[15] = "pt"
  msh_elementtype_nodes[15] = 1
  # Node order seems to be equivalent. (gmsh, GiD, VTK)
  # Unsure about msh_PRISM = VTK_WEDGE.
  msh_elementtype_str[1] = "lines"
  msh_elementtype_str[2] = "triangles"
  msh_elementtype_str[3] = "quadrangles"
  msh_elementtype_str[4] = "tetrahedra"
  msh_elementtype_str[5] = "hexahedra"
  msh_elementtype_str[6] = "prisms"
  msh_elementtype_str[7] = "pyramids"
  msh_elementtype_str[15] = "points"
}

(FNR == 1) {
  # Handle multiple files.
  if (NR > 1) {
    # New file, output current data.
    check_msh_finished()
    if (inp_filter_output) print_inp_filtered()
    else print_inp()
  }

  # Reset.
  msh_count_meshformat = 0
  msh_version = 0
  msh_type = 0
  msh_datasize = 0
  msh_count_nodes = 0
  delete nodes
  delete nodeX
  delete nodeY
  delete nodeZ
  msh_count_elements = 0
  delete elements
  delete elementType
 msh_count_type[1] = 0
 msh_count_type[2] = 0
 msh_count_type[3] = 0
 msh_count_type[4] = 0
 msh_count_type[5] = 0
 msh_count_type[6] = 0
 msh_count_type[7] = 0
 msh_count_type[15] = 0
  delete elementNodes
 delete msh_node_usedtype
  delete elementTags
  delete elementPhysical
  delete elementGeometry
  delete elementPartition
  delete elementNode
  msh_count_physicalnames = 0
  delete physicalnames
  delete physicalnameName
  msh_part_NR = 0

  # Filename handling.
  if (FILENAME == "-") {
    inpfile = "/dev/stdout"
  } else {
    inpfile = FILENAME
    if (inpfile ~ /\.msh$/) {
      sub(/\.msh$/, ".inp", inpfile)
    } else {
      inpfile = inpfile ".inp"
    }
    if (FILENAME == inpfile) {
      print "Filename " FILENAME " could not be handled." > "/dev/stderr"
      error = 1
      exit error
    }
  }
  if (DEBUG) {
    print "Reading from \"" FILENAME "\", output to \"" inpfile "\"." > "/dev/stderr"
  }
}

# Detect whicht part of the file we are about to read.
(msh_read_part == 0) {
  if ($0 == "$MeshFormat") {
    msh_read_part = msh_part_meshformat
  } else
  if ($0 == "$Nodes") {
    msh_read_part = msh_part_nodes
  } else
  if ($0 == "$Elements") {
    msh_read_part = msh_part_elements
  } else
  if ($0 == "$PhysicalNames") {
    msh_read_part = msh_part_physicalnames
  } else
#  if ($0 == "$NodeData") {
#    msh_read_part = msh_part_nodedata
#  } else
#  if ($0 == "$ElementData") {
#    msh_read_part = msh_part_elementdata
#  } else
#  if ($0 == "$ElementNodeData") {
#    msh_read_part = msh_part_elementnodedata
#  } else
  if ($0 ~ /^\$/) {
    msh_read_part = msh_part_unknown
  }

  msh_part_end = $0
  sub(/^\$/, "$End", msh_part_end)

  msh_part_NR = 0
  next
}

# Read header.
# version-number file-type data-size
(msh_read_part == msh_part_meshformat) {
  if ($0 == msh_part_end) {
    msh_read_part = 0
    next
  }
  if (msh_count_meshformat == 0) {
    msh_version = $1
    msh_type = $2
    msh_datasize = $3
    msh_count_meshformat++
    if (msh_version != 2) {
      print "Unknown file version on line " FNR "." > "/dev/stderr"
      error = 2
      exit error
    }
    if (msh_type != 0) {
      print "Unknown file type on line " FNR "." > "/dev/stderr"
      error = 3
      exit error
    }
  } else {
    print "Multiple format specifications found on line " FNR "." > "/dev/stderr"
    error = 4
    exit error
  }
  next
}

# Read nodes.
# number-of-nodes
# node-number x-coord y-coord z-coord
# ...
(msh_read_part == msh_part_nodes) {
  if ($0 == msh_part_end) {
    if (msh_nodes_to_read > 0) {
      print "Node missing on line " FNR "." > "/dev/stderr"
      error = 5
      exit error
    }
    msh_read_part = 0
    next
  }
  msh_part_NR++
  if (msh_part_NR == 1) {
    msh_nodes_to_read = $1
    if (DEBUG) {
      print "Reading " msh_nodes_to_read " nodes." > "/dev/stderr"
    }
    next
  }
  if (msh_nodes_to_read < 1) {
    print msh_part_end " missing in \"" FILENAME "\" on line " FNR "." > "/dev/stderr"
    error = 6
    exit error
  }
  msh_count_nodes++
  nodes[msh_count_nodes] = $1
  nodeX[$1] = $2
  nodeY[$1] = $3
  nodeZ[$1] = $4
  msh_nodes_to_read--
  next
}

# Read elements.
# number-of-elements
# elm-number elm-type number-of-tags < tag > ... node-number-list
# ...
(msh_read_part == msh_part_elements) {
  if ($0 == msh_part_end) {
    if (msh_elements_to_read > 0) {
      print "Element missing on line " FNR "." > "/dev/stderr"
      error = 7
      exit error
    }
    msh_read_part = 0
    next
  }
  msh_part_NR++
  if (msh_part_NR == 1) {
    msh_elements_to_read = $1
    if (DEBUG) {
      print "Reading " msh_elements_to_read " elements." > "/dev/stderr"
    }
    next
  }
  if (msh_elements_to_read < 1) {
    print msh_part_end " missing on line " FNR "." > "/dev/stderr"
    error = 8
    exit error
  }
  msh_count_elements++
  elements[msh_count_elements] = $1
  elementType[$1] = $2
  msh_count_type[$2]++
  elementNodes[$1] = msh_elementtype_nodes[$2]
  if (!elementNodes[$1]) {
    print "Unknown element type on line " FNR "." > "/dev/stderr"
    error = 9
    exit error
  }
  elementTags[$1] = $3
  if ($3 != 3) {
    print "Can handle 3 default tags only (line " FNR ")." > "/dev/stderr"
    error = 10
    exit error
  }
  elementPhysical[$1] = $4
  elementGeometry[$1] = $5
  elementPartition[$1] = $6
  j = 6
  for (i = 3; i <= $3; i++) {
    j = i + 3
  }
  for (i = 1; i <= elementNodes[$1]; i++) {
    k = i + j
    elementNode[$1,i] = $k
    msh_node_usedtype[$k,$2]++
  }
  msh_elements_to_read--
  next
}

# Read physical names.
# number-of-names
# physical-number "physical-name"
# ...
(msh_read_part == msh_part_physicalnames) {
  if ($0 == msh_part_end) {
    if (msh_physicalnames_to_read > 0) {
      print "Physical name missing on line " FNR "." > "/dev/stderr"
      error = 11
      exit error
    }
    msh_read_part = 0
    next
  }
  msh_part_NR++
  if (msh_part_NR == 1) {
    msh_physicalnames_to_read = $1
    print "Reading " msh_physicalnames_to_read " physical names." > "/dev/stderr"
    next
  }
  if (msh_physicalnames_to_read < 1) {
    print msh_part_end " missing on line " FNR "." > "/dev/stderr"
    error = 12
    exit error
  }
  msh_count_physicalnames++
  physicalnames[msh_count_physicalnames] = $1
  physicalnameName[$1] = $2
  msh_physicalnames_to_read--
  next
}

# Skip unknown part.
(msh_read_part == msh_part_unknown) {
  if ($0 == msh_part_end) {
    msh_read_part = 0
    next
  }
  next
}

# Check if all parts have been read.
function check_msh_finished(        i, s) {
  if (msh_read_part != 0) {
    print msh_part_end " missing on line " FNR "." > "/dev/stderr"
    error = 13
    exit error
  }
  if (DEBUG) {
    print "Nodes:          " msh_count_nodes > "/dev/stderr"
    print "Elements:       " msh_count_elements > "/dev/stderr"
    print "Physical Names: " msh_count_physicalnames > "/dev/stderr"
    print "Points:         " msh_count_type[15] > "/dev/stderr"
    print "Lines:          " msh_count_type[1] > "/dev/stderr"
    print "Triangles:      " msh_count_type[2] > "/dev/stderr"
    print "Quadrangles:    " msh_count_type[3] > "/dev/stderr"
    print "Tetrahedra:     " msh_count_type[4] > "/dev/stderr"
    print "Hexahedra:      " msh_count_type[5] > "/dev/stderr"
    print "Prisms:         " msh_count_type[6] > "/dev/stderr"
    print "Pyramids:       " msh_count_type[7] > "/dev/stderr"
  }
  s = 0
  for (i in msh_count_type) {
    s += msh_count_type[i]
  }
  if (s != msh_count_elements) {
    print "Sum of element counts not equal to total element count in \"" FILENAME "\"." > "/dev/stderr"
    error = 14
    exit error
  }
}

# Print inp output of current data.
function print_inp(        i,inode,icell,j) {
  if (DEBUG) {
    print "No output filtering." > "/dev/stderr"
    print "Output " msh_count_nodes " nodes and " msh_count_elements " cells." > "/dev/stderr"
  }
  # Header: nnodes, ncells, nnodedata, ncelldata, nmodeldata
  print msh_count_nodes, msh_count_elements, 0, 0, 0 > inpfile
  # Nodes: n, x, y, z
  for (i = 1; i <= msh_count_nodes; i++) {
    inode = nodes[i]
    print inode, nodeX[inode], nodeY[inode], nodeZ[inode] > inpfile
  }
  # Cells: n, material, type, nodes
  for (i = 1; i <= msh_count_elements; i++) {
    icell = elements[i]
    oldORS = ORS
    ORS = ""
    print icell, elementPhysical[icell], msh_elementtype_inp[elementType[icell]] > inpfile
    for (j = 1; j <= elementNodes[icell]; j++) {
      print "", elementNode[icell,j] > inpfile
    }
    ORS = oldORS
    print "" > inpfile
  }
  if (DEBUG) print "" > "/dev/stderr"
}

# Print filtered inp output of current data.
function print_inp_filtered(        oldORS,t,ncells,i,icell,cell_used,cell_map,nnodes,inode,node_used,node_map,node_remap,j) {
  if (DEBUG) {
    oldORS = ORS
    ORS = ""
    print "Output filtering enabled: " > "/dev/stderr"
    for (t in inp_output) {
      if (inp_output[t]) print msh_elementtype_str[t] " " > "/dev/stderr"
    }
    ORS = oldORS
    print "" > "/dev/stderr"
  }

  # Find filtered cells and create mapping.
  ncells = 0
  delete cell_used
  delete cell_map
  for (i = 1; i <= msh_count_elements; i++) {
    icell = elements[i]
    for (t in inp_output) {
      if (inp_output[t]) {
        if (elementType[icell] == t) {
          cell_used[icell]++
          if (cell_used[icell] == 1) {
            ncells++
            cell_map[ncells] = icell
          }
        }
      }
    }
  }
  if (DEBUG) {
    print ncells " cells found." > "/dev/stderr"
  }

  # Find filtered nodes and create mapping.
  nnodes = 0
  delete node_used
  delete node_map
  for (i = 1; i <= msh_count_nodes; i++) {
    inode = nodes[i]
    for (t in inp_output) {
      if (inp_output[t]) {
        if (msh_node_usedtype[inode,t]) {
          node_used[inode]++
          if (node_used[inode] == 1) {
            nnodes++
            node_map[nnodes] = inode
          }
        }
      }
    }
  }
  if (DEBUG) {
    print nnodes " nodes found." > "/dev/stderr"
    print "Output " nnodes " nodes and " ncells " cells." > "/dev/stderr"
    if (inp_renumber) print "Renumbering used nodes and cells." > "/dev/stderr"
  }

  delete node_remap
  if (inp_renumber) {
    # Header: nnodes, ncells, nnodedata, ncelldata, nmodeldata
    print nnodes, ncells, 0, 0, 0 > inpfile
    # Nodes: n, x, y, z
    for (i = 1; i <= nnodes; i++) {
      inode = node_map[i]
      node_remap[inode] = i
      print i, nodeX[inode], nodeY[inode], nodeZ[inode] > inpfile
    }
    # Cells: n, material, type, nodes
    for (i = 1; i <= ncells; i++) {
      icell = cell_map[i]
      oldORS = ORS
      ORS = ""
      print i, elementPhysical[icell], msh_elementtype_inp[elementType[icell]] > inpfile
      for (j = 1; j <= elementNodes[icell]; j++) {
        print "", node_remap[elementNode[icell,j]] > inpfile
      }
      ORS = oldORS
      print "" > inpfile
    }
  } else {
    # Header: nnodes, ncells, nnodedata, ncelldata, nmodeldata
    print nnodes, ncells, 0, 0, 0 > inpfile
    # Nodes: n, x, y, z
    for (i = 1; i <= msh_count_nodes; i++) {
      inode = nodes[i]
      if (node_used[inode]) {
        print inode, nodeX[inode], nodeY[inode], nodeZ[inode] > inpfile
      }
    }
    # Cells: n, material, type, nodes
    for (i = 1; i <= msh_count_elements; i++) {
      icell = elements[i]
      if (cell_used[icell]) {
        oldORS = ORS
        ORS = ""
        print icell, elementPhysical[icell], msh_elementtype_inp[elementType[icell]] > inpfile
        for (j = 1; j <= elementNodes[icell]; j++) {
          print "", elementNode[icell,j] > inpfile
        }
        ORS = oldORS
        print "" > inpfile
      }
    }
  }
  if (DEBUG) print "" > "/dev/stderr"
}

END {
  if (error == -127) exit 0
  if (error) exit error

  check_msh_finished()
  # Print inp.
  if (inpfile) {
    if (inp_filter_output) print_inp_filtered()
    else print_inp()
  }
}

# vim:set sts=2 sw=2: