A water tank mesh generater in R

by Yi Zhang

I have been using LS-DYNA ICFD solver for several years. Basically, it’s a incompressible CFD solver addon over the platform of LS-DYNA, which has its original flow solver based on compressible flows. It’s nice to work with people in LSTC, the develop company of LS-DYNA, but their pre-post processor really doesn’t work well with me. For my own research, during the development of a numerical wave tank, most of the time what I need is just a tank mesh. In ICFD, mesh generation is done automatically within the flow domain, so the mesh information that is requried at the file input is just the meshed surface: the walls, the free surface, the bottom, and the other  domain boundaries. Since open lsprepost to build a mesh file from scratch is the last thing I’d like to do when every time certain mesh parameters are changed, such as in convergence test, I just write a script to do that work. The output plots for mesh example from this R script are like this:

The first plot shows all the faces of generated. Plot 2 shows the mesh for bottom, free surface and the “ceiling”, together with “side walls”. Plot 3 shows free surface mesh and two ends of the tank, where the wave are generated and dissipated.

R turns out to be handy for this job, for its way of handling data, especially the recycling and the mapping due to its embracing of functional paradigm. As the script below shows. The code may be fragile, still saves a lot of time from GUI.

##########################################################################
## Copyright (C) 2012-2013 Yi Zhang                                     ##
##                                                                      ##
## Author: Yi Zhang                                                     ##
## Created: 30 Nov 2012                                                 ##
##                                                                      ##
## Purpose: This file generates a surface mesh as used by ICFD          ##
## in LS-DYNA keyfile format.                                           ##
##                                                                      ##
## This file is free software; you can redistribute it and/or modify    ##
## it under the terms of the GNU General Public License as published by ##
## the Free Software Foundation; either version 2, or (at your option)  ##
## any later version.                                                   ##
##########################################################################

rm(list=ls())
meshfile <- "~/tmp/tank3d_mesh.k"

d <- c(10., 1., 0.5, 0.3) # dim: (x, y, z, z-air)
m <- c(75,14,10,10)           # nb. of elements in (x,y,z,z1)
n <- m + 1                # nb. of nodes
h <- d/m

## total nb. of elems. & nodes.
mtol <- 3*m[1]*m[2] + 2*(m[2]*(m[3]+m[4]) + m[1]*(m[3]+m[4]))
ntol <- 3*n[1]*n[2] + 2*(n[1]*(m[3]+m[4]-2)++(m[2]-1)*(m[3]+m[4]-2))

## beginning and the end id's of the each face
facegen <- function (members) {
  faces <- matrix(c(1, members[1]), nrow=1, ncol=2)
  for (i in 2:length(members)) {
    faces <- rbind(faces, c(tail(faces,1)[2]+1,
                            tail(faces,1)[2]+members[i]))}
  return(faces)}
edgegen <- function (i, k) {
  ## j = (1, 2, 3, 4)  (S, E, N, W)
  edge <- list(nface[i,1] : (nface[i,1]+m[k]),
               seq(nface[i,1]+m[k], nface[i,2], n[k]),
               (nface[i,2]-m[k]) : nface[i,2],
               seq(nface[i,1], nface[i,2], n[k]))
  return(edge) }
nbn   <- c(n[1]*n[2], n[1]*n[2], n[1]*n[2], n[1]*n[3],
           n[1]*n[3], n[1]*n[4], n[1]*n[4], n[2]*n[3],
           n[2]*n[3], n[2]*n[4], n[2]*n[4]) 
nbm   <- c(m[1]*m[2], m[1]*m[2], m[1]*m[2], m[1]*m[3],
           m[1]*m[3], m[1]*m[4], m[1]*m[4], m[2]*m[3],
           m[2]*m[3], m[2]*m[4], m[2]*m[4]) 
nface <- facegen(nbn)
mface <- facegen(nbm)
nedge <- list(edgegen(1, 1), edgegen(2, 1), edgegen(3, 1),
              edgegen(4, 1), edgegen(5, 1), edgegen(6, 1),
              edgegen(7, 1), edgegen(8, 2), edgegen(9, 2),
              edgegen(10, 2), edgegen(11, 2))

## nodes
x   <- rep(seq(0, d[1], h[1]), n[2])
y   <- as.vector(sapply(seq(-d[2]/2, d[2]/2, h[2]), function (x) rep(x, n[1])))
nod <- rbind(cbind(x, y, -d[3]),
             cbind(x, y, 0.   ),
             cbind(x, y,  d[4]))

x   <- rep(seq(0, d[1], h[1]), n[3])
z   <- as.vector(sapply(seq(-d[3], 0, h[3]), function (x) rep(x, n[1])))
nod <- rbind(nod,
             cbind(x, -d[2]/2, z),
             cbind(x,  d[2]/2, z))

x   <- rep(seq(0, d[1], h[1]), n[4])
z   <- as.vector(sapply(seq(0, d[4], h[4]), function (x) rep(x, n[1])))
nod <- rbind(nod,
             cbind(x, -d[2]/2, z),
             cbind(x,  d[2]/2, z))

y   <- rep(seq(-d[2]/2, d[2]/2, h[2]), n[3])
z   <- as.vector(sapply(seq(-d[3], 0, h[3]), function (x) rep(x, n[2])))
nod <- rbind(nod,
             cbind(0.,   y, z),
             cbind(d[1], y, z))

y   <- rep(seq(-d[2]/2, d[2]/2, h[2]), n[4])
z   <- as.vector(sapply(seq(0, d[4], h[4]), function (x) rep(x, n[2])))
nod <- rbind(nod,
             cbind(0.,   y, z),
             cbind(d[1], y, z))

nod <- cbind(1:(dim(nod)[1]), nod)

## elements
f <- function (i, j) {
  p <- nface[i,1] : nface[i,2]
  a <- p[(p-nface[i,1]+1) %% n[j] != 0];
  b <- p[(p-nface[i,1]+1) %% n[j] != 1]
  return(cbind(a[a <  nface[i,2] - n[j]],
               b[b   nface[i,1] + n[j]],
               a[a >  nface[i,1] + n[j] - 1]))
}
ele <- cbind(1:mtol,
             unlist(sapply(1:11, function (i) rep(i,
                                                  mface[i,2]-mface[i,1]+1))), 
             rbind(f(1, 1), f(2, 1), f(3, 1), f(4, 1),
                   f(5, 1), f(6, 1), f(7, 1),
                   f(8, 2), f(9, 2), f(10, 2), f(11, 2)))

## Eliminate redundant nodes (at corners/edges)
e <- ele[,3:6]
nr <- function (ma, a, b) {
  for(i in 1:length(a)) ma[ma == a[i]] <- b[i]
  return(ma) }
a <- rbind(c(11, 2, 7, 2), c(11, 4, 6, 2),
           c(10, 2, 7, 4), c(10, 4, 6, 4),
           c(9 , 2, 5, 2), c(9 , 4, 4, 2),
           c(8 , 2, 5, 4), c(8 , 4, 4, 4),
           c(11, 1, 2, 2), c(11, 3, 3, 2),
           c(10, 1, 2, 4), c(10, 3, 3, 4),
           c(9 , 1, 1, 2), c(9 , 3, 2, 2),
           c(8 , 1, 1, 4), c(8 , 3, 2, 4),
           c(7 , 1, 2, 3), c(7 , 3, 3, 3),
           c(6 , 1, 2, 1), c(6 , 3, 3, 1),
           c(5 , 1, 1, 3), c(5 , 3, 2, 3),
           c(4 , 1, 1, 1), c(4 , 3, 2, 1))
for(i in 1:dim(a)[1]){
  f1 <- a[i,1]; s1 <- a[i,2]; f2 <- a[i,3]; s2 <- a[i,4]
  e <- nr(e, unlist(nedge[[f1]][s1]), unlist(nedge[[f2]][s2]))}

nid <- sort(unique(as.vector(e)))
nod <- nod[nid,]
nod[,1] <- 1:ntol
e <- nr(e, nid, 1:ntol)
ele[,3:6] <- e

## Write file in LS-DYNA keyfile format
write("*KEYWORD", file=meshfile) 
write("*MESH_SURFACE_ELEMENT", file=meshfile, append=TRUE) 
write.table(cbind(format(ele[,1], width=8,justify="right"),
                  format(ele[,2:6],
                         width=7,justify="right")),
            file=meshfile, append=TRUE, quote=FALSE,
            row.names=FALSE, col.names=FALSE)  
write("*MESH_SURFACE_NODE", file=meshfile, append=TRUE) 
write.table(cbind(format(nod[,1],
                         width=8, digits=1, justify="right"),
                  format(nod[,2:4], width=15, digits=7,
                         justify="right")),
            file=meshfile, append=TRUE, quote=FALSE,
            row.names=FALSE, col.names=FALSE)  
write("*END", file=meshfile, append=TRUE)
Advertisements