Generating an OpenFOAM mesh from a point cloud (x,y,z)

I was recently provided a file containing x,y,z points defining the bed elevation extracted from a flume study. I needed a way to efficiently convert this to an OpenFOAM mesh file while having control over the mesh resolution. I am not that familiar with the snappyHexMesh program (yet) but one option would be to convert the point cloud to an STL file. The STL file can then be used to generate a mesh using snappyHexMesh. Note that there is a nice program in Matlab that will perform the conversion from a point cloud to STL for you. Incidentally there is a great tutorial on using snappyHexMesh that can be found here.

Since I am familar with the blockMesh utility and GIS, I decided to try another approach. Using ArcGIS I converted the point cloud first to a TIN and then a DEM. A DEM is essentially a regular matrix of elevations and which makes it very easy to work with as opposed to a randomly spaced point cloud. I then exported the DEM to an ASCII DEM. Finally I wrote a PYTHON script to read the ASCII DEM and write out the blockMeshDict file. The script can be found below. Currently the script requires to user to specify the number of rows and columns as well as the cell size however these could also be read from the ASCII DEM header. The script is nice in that it is easy to make modifications to extract subsets. For example I have a version used to extract a 2D slice of the bed.

There is one issue with the code to be aware of. The code treats each centroid of the DEM as a vertex in the blockMesh which is not quite correct. This may or may not be an issue depending on your problem, but can be corrected. If I get around to correcting this I will post the code.

Here is an image of the mesh I was able to generate for my point cloud:

'''
Created on May 10, 2013

@author: Patrick Grover
'''
from numpy import *

# note that matrix starts at 1,1 and goes up to numrows,numcols
def getvertex(row, col, numrows):
    index = numrows*(col-1) + row
    return index

def getBlock(row, col, numrows,numcols, blockSettings):
    retval =  "\thex ("
    retval = retval + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1) + " "    
    
    numvertices = numrows*numcols
    
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + ") "
    retval = retval + blockSettings     
    return retval

def getTopWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + ") "
    return retval

def getBottomWall(row, col, numrows,numcols):
    retval = "(" + str(getvertex(row,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows)-1 ) + " "    
    retval = retval + str(getvertex(row,col+1,numrows)-1 )   + ") "
    return retval

def getInlet(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col,numrows)+ numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1)    
    retval = retval + ") "
    return retval

def getOutlet(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1)      
    retval = retval + ") "
    return retval

def getLeftWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "  
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1)      
    retval = retval + ") "
    return retval

def getRightWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col,numrows)+ numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1)
    retval = retval + ") "
    return retval


#--------------------------------------------------------------------------------
# BlockMesh Parameters

# Manually set these. Could read these from the ASCII header.
depth = 2.0
numrows = 67
numcols = 1000
cellsize = 0.03   

# Meshing details - describes the mesh
blockSettings = "(2 1 20) simpleGrading (1 1 10)\n"

# Input ASCII DEM

asciiFilePath = '/home/geeko/working/GIS/dunes_3cm_67x1000.asc'
# Path to the blockMeshDict

blockMeshPath = "/home/geeko/working/BlockMesh/TestMesh2/constant/polyMesh/blockMeshDict"

# Read in ASCII file
data = genfromtxt(asciiFilePath,delimiter=' ',skiprows=6,unpack=True, filling_values=-9999)

print data

fo= open(blockMeshPath,'w')
fo.write("FoamFile\n{version\t2.0;\n\tformat\tascii;\n\tclass\tdictionary;\n\tobject\tblockMeshDict;\n}\n\n")

fo.write("convertToMeters 1.0;\n\n")


# Write out vertices
fo.write("vertices\n")
fo.write("(\n")

for y in range(0,numcols):
    for x in range(0,numrows):
        #print "x= " + str(x) + "y= " + str(y)
        if data[y,x]==-9999:
            print "!!!!!!!!!!!!!!!! Error!!!!!!!!!!!!!!!!!!!!!"
        #fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(data[numrows-y-1,x]) + ")\n")
        fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(data[y-1,x]) + ")\n")
for y in range(0,numcols):
    for x in range(0,numrows):
        fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(depth) + ")\n")
    
fo.write(");\n")



# Write out blocks
fo.write("blocks\n")
fo.write("(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write(getBlock(row, col, numrows,numcols, blockSettings))
fo.write(");\n")
fo.close

#edges
fo.write("edges\n")
fo.write("(\n")
fo.write(");\n");

# Boundaries
fo.write("boundary\n")
fo.write("(\n")


# Top Wall
fo.write("\ttopWall\n")
fo.write("\t{\n")
fo.write("\ttype symmetryPlane;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write("\t\t" + getTopWall(row, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# getBottomWall
fo.write("\tbottomWall\n")
fo.write("\t{\n")
fo.write("\ttype wall;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write("\t\t" + getBottomWall(row, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# Inlet
fo.write("\tinlet\n")
fo.write("\t{\n")
fo.write("\ttype patch;\n")
#fo.write("\ttype cyclic;\n")
#fo.write("\tneighbourPatch outlet;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range(1,numrows):
    fo.write("\t\t" + getInlet(row, 1, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# Outlet
fo.write("\toutlet\n")
fo.write("\t{\n")
fo.write("\ttype patch;\n")
#fo.write("\ttype cyclic;\n")
#fo.write("\tneighbourPatch inlet;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range(1,numrows):
    fo.write("\t\t" + getOutlet(row, numcols, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")



# Front and Back
fo.write("\tfrontAndBack\n")
fo.write("\t{\n")
fo.write("\ttype empty;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
fo.write("\t// Left side\n")
for col in range (1,numcols):
    fo.write("\t\t" + getLeftWall(1, col, numrows,numcols) + "\n")
fo.write("\t// Right side\n")
for col in range (1,numcols):
    fo.write("\t\t" + getRightWall(numrows, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

fo.write(");\n")


#mergePatchPairs
fo.write("mergePatchPairs\n")
fo.write("(\n")
fo.write(");\n");


print "Done!"