Mesh Mess

FEM, CFD, CAE, and all that.

Tag: Mesh generation

A water tank mesh generater in R

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)

You get what you mesh for

In finite element modeling, in fact, in all the numerical modeling, the resolution of the results depends on the space where the problem is projected. As to PDE’s numerical solution, the projection is from space with infinite dimension to space with finite dimension. So no matter it’s finite element or finite volume or finite difference, the size of the mesh is what you can see at the finest scale. Beyond that, it’s either physical modeling (such as turbulence subscale modeling) or some other educated guess.

Here is an example.

Water surface as wave coming in

This plot shows a finite element solution of Navier-Stokes equation, in a water waves problem. The free surface is tracked using Level Set method. Taking closer look reveals the following not-so-smooth “water surface”:

Free surface, enlarged

The resolution of the free surface is restricted by the size of the element. For any mesh, examing close enough reveals similar wiggles. To put it another way, in the numerical simulation, free surface is not really a surface, it’s more like a layer with finite thickness. If I plot the centroid of the elements that reside on the free surface, it looks like this:

Elements on the free surface

This is what a non-polished numerical result should look like. Smoothing such as with cubic splines sweepts everything under the table, including the message: “you get what you mesh for”.

3D alpha shape test

Example of isolating tetrahedron elements within the boundary of a cube ant its inscribed sphere. The tetrahedra element edge length is around 0.5 in general. It is indicated that the higher resolution of boundary recognition is provided by alpha=0.5, close to the element edge length.

Three types of elements with alpha=0.5

Three types of elements with alpha=0.5

alpha=0.5

alpha=0.5

alpha=1

alpha=1

alpha=2

alpha=2

2D Alpha shape test

An example of  implementing alpha shape method in boundary recognition, using Mathematica.

Boundary recognition

Based on the idea of alpha shape, one can devise a boundary recognition scheme by eliminating the elements whose \alpha  circumsphere contains no nodes which are within the boundary. This means, first one must have those nodes identified, then those nodes and the coefficient \alpha (defines the extended version of circumsphere) are used for boundary recognition. Here is a example of identifying the boundary of those nodes lying inside a cube but outside its inscribed sphere.

mesh_cube_inscb_para1mesh_cube_inscb_para2mesh_cube_inscb_para3

The first four plots show original meshed cube with inside view, the next six plots are when elements inside the sphere are eliminated. The rest is what we have when we change \alpha, i.e. the coefficient determining whether one tetrahedron should be eliminated because its extended circumsphere, whose radius depends on \alpha, contains nodes within the boundary. Here \alpha is reduces, resulting in more elements erased from the picture.