Mesh Mess

FEM, CFD, CAE, and all that.

Prohibit the metric system!

Through Irreal I found this eloquent argument against metric system in favor of US system ( I don’t call it imperial since UK has officially adopted metric). 8 signatures have been gathered toward 25000 goal, astonishingly less than the ratio of nations endorsing the almighty US system (which is just US, essentially).

Advertisements

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)

The citation management in MS Word is from hell

I was going to say the MS Word by itself is from hell, at least for technical writings. Recently I had to convert my tex writing into word, everything was tolerable, including the conversion of the equations since Mathtype recognizes tex, until at the end when I was converting the reference and citations. There is a “manage source” option in the Word, but it only works for Endnote. For the record, that is bullshit: one product tries to force its user to adopt another equally pathetic commercial product? What an ally!

Then I decited to do it myself.

I regetted it in five minutes.

For any bibtex user, I can not even describe the frustration looking at the fancy menus that Word’s reference management “ribbon” (yeah, they really named it that way!) demonstrates. I don’t think any latex user who hasn’t lost his mind would achieve the task of converting 100+ references manually. For those who are familiar with this kind of problem, the bib2endnote by Trend Apted didn’t work: somehow the generated XML was denied by the Endnote parser in the Word.

So what did I do? I copied the texts of the bibliography and pasted them in the versatile Word, since I don’t need to edit them in the fugure anyway: it’s for whoever asks for the word version, no matter how pathetic it is.

Convert pdf to png

I’ve realized at a point of my life that to prepare a scientific writing in MS Word (if you have to), it takes longer in doing it completely in Word, than doing it using emacs+Auctex+latex then transfering the document to .doc file manually, since no pdf-word converters work perfectly for equations.¬† This usually happens when someone makes unrefusable¬† requests for .doc files, such as demanded by my advisor.

One thing I need to do during this process is to convert all the pdf plots generated by R and Gnuplot to bitmaps. This is when imagemagick comes in. Without pain of GUI, one can just do

for f in *.pdf; do convert -density 800x800 -quality 90 $f ${f/.pdf/.png};done

in bash.

Store a polynomial in non-symbolic language

Sometimes one wants to store a polynomial in a non-symbolic language, such as Fotran or C. A easy way to achieve this is to take advantage of the isomorphism between polynomial space P(n) and space R^{n+1}. Namely, we can build an array with entries as the coefficients of the polynomial to be stored. Furthermore, considering these coefficients are just the coordinates of the polynomial in the space spanned by the canonical basis 1,x,x^2,\dots, we can express the polynomial using any basis expansion.