R/gmsh.R

In telemac: R Interface to the TELEMAC Model Suite

#' Read Gmsh file
#'
#' Reads the mesh file generated by Gmsh (only format version 2 supported).
#'
#' @param fname \code{character}, name of a Gmsh file.
#'
#' @return A \code{list} with the following elements:
#' \describe{
#'   \item{nelem}{Number of mesh elements (triangles)}
#'   \item{npoin}{Number of mesh points}
#'   \item{ikle}{A \code{nelem x 3} integer matrix of point indices (referring to \code{x} and \code{y}) defining the mesh elements}
#'   \item{ipobo}{An integer vector of length \code{npoin} defining the mesh boundaries (inner boundaries are zero, outer boundaries numbered)}
#'   \item{x}{A vector of length \code{npoin} giving the x coordinates of the mesh points}
#'   \item{y}{A vector of length \code{npoin} giving the y coordinates of the mesh points}
#' }
#'
#' @note So far only mesh file format version 2 is supported. Besides, only elements
#' 'line' (considered as boundary) and 'triangle' (treated as actual mesh elements)
#' can be handled.
#'
#' @export
read_msh <- function(fname) {
  # read data
  dat <- readLines(fname)

  # check file format version
  gmshver <- tonum(dat[grep("$MeshFormat", dat, fixed = T)+1])[1]
  if (floor(gmshver) != 2)
    stop("Only MSH file format version 2 is supported!", call. = F)

  # get node points
  inodes <- grep("$Nodes", dat, fixed = T)
  npoin <- as.numeric(dat[inodes + 1])
  nodes <- sapply(dat[seq(inodes + 2, inodes + npoin + 1)], tonum, USE.NAMES = F)

  # analyse elements
  ielem <- grep("$Elements", dat, fixed = T)
  nelem_t <- as.numeric(dat[ielem + 1])
  jelem <- seq(ielem + 2, ielem + nelem_t + 1)
  types <- sapply(dat[jelem], function(x) tonum(x)[2], USE.NAMES = F)
  if (any(types > 2))
    stop("Can only handle element types 'line' and 'triangle' in MSH file!", call. = F)

  # get boundary nodes (ipobo)
  nbnd <- length(which(types == 1))
  ipobo <- sapply(dat[jelem[types == 1]], function(x) as.integer(tonum(x)[1]), USE.NAMES = F)

  # get triangle definitions (ikle)
  nelem <- length(which(types == 2))
  ikle <- sapply(dat[jelem[types == 2]], function(x) as.integer(tail(tonum(x), 3)), USE.NAMES = F)

  # output
  return(list(nelem = as.integer(nelem),
              npoin = as.integer(npoin),
              ikle = t(ikle),
              ipobo = ipobo,
              x = nodes[2,],
              y = nodes[3,]))
}