Creating a surface from structured (x,y,z) topography data and meshing it

This post is a follow-up to a question in today’s webinar “Introduction to Python Scripting in Coreform Cubit” and is based upon this previous forum post.

The question from the webinar was:

How can I import a file containing structured x, y, z data (of a topographical surface) and create a CAD geoemtry from it – and mesh it?

Below is a Python script, that can be executed in Coreform Cubit, that demonstrates this. Note that I’m creating a “dat-file” that looks like:

0,-10.0,-10.0,0.020519746021903675
0,-10.0,-9.591836734693878,-0.006596306823113194
...
0,-10.0,10.0,0.007623109004567477
1,-9.591836734693878,-10.0,-0.021575477096699273
...
1,-9.591836734693878,10.0,0.030617838494543526
2,-9.183673469387756,-10.0,-0.004938534571641062
...
...
...
49,10.0,9.591836734693878,-0.03863820315232227
49,10.0,10.0,0.005197315777777938

This is a comma-separated values (CSV) file, with no header row, where each line encodes the data for a single vertex:

row_id, x, y, z

While you could define this file differently, it’s valuable that we have the row_id information to let us know how to connect these vertices together. The idea is represented in this image, where I’ve colored each vertex by the Cubit Color corresponding to each value of row_id:

This provides needed structure to identify how to create spline curves connecting each “line” of vertices, which we do in the command:

cubit.cmd( f"create curve spline location vertex {cubit.get_id_string( vertices_in_curve )} delete" )

and then we create the surface via “curve skinning” with the command:

cubit.cmd( f"create surface skin curve {cubit.get_id_string( curves_in_surface )}" )

Note that I’m utilizing several helpful Cubit-Python methods throughout the script:

  • cubit.cmd( <string> )
    • Executes the input string as a Cubit command
  • cubit.get_last_id( "<entity_type>" )
    • Gets the Cubit ID of the last-created entity_type, passed in as a string – e.g., "curve"
  • cubit.get_id_string( <list> )
    • Converts a list of integers into a Cubit-parseable string – e.g., [0, 1, 2] --> "0 1 2"

#!python
import numpy
import random
import pathlib

filename = pathlib.Path( "~/example_topography_data.csv" ).expanduser()

def main( N, mesh_size ):
  set_performance_mode( True )
  cubit.cmd( "reset" )
  create_dat_file( N, filename )
  data = read_dat_file( filename )
  create_surface( data )
  surface_id = cubit.get_last_id( "surface" )
  mesh_surface( surface_id, mesh_size )
  set_performance_mode( False )

def create_dat_file( N, dat_filename ):
  x = numpy.linspace( -10, 10, N )
  y = numpy.linspace( -10, 10, N )
  with open( dat_filename, "w+" ) as f:
    for row_id in range( 0, len( x ) ):
      for yi in y:
          z = 0.08 * ( random.random() - 0.5 )
          f.write( f"{row_id},{x[row_id]},{yi},{z}\n" )

def read_dat_file( dat_filename ):
  data = {}
  x, y, z = [ [], [], [] ]
  f = open( dat_filename, "r" )
  fLines = f.readlines()
  for fLine in fLines:
    fLine = fLine.strip().split( "," )
    row_id = int( fLine[0] )
    if ( row_id in data.keys() ) == False:
      data[row_id] = { "x": [], "y": [], "z": [] }
    data[row_id]["x"].append( float( fLine[1] ) )
    data[row_id]["y"].append( float( fLine[2] ) )
    data[row_id]["z"].append( float( fLine[3] ) )
  return data

def create_surface( data ):
  row_ids = list( data.keys() )
  curves_in_surface = []
  for row_id in row_ids:
    row_data = data[row_id]
    x, y, z = [ row_data["x"], row_data["y"], row_data["z"] ]
    vertices_in_curve = []
    for i in range( 0, len( x ) ):
      cubit.cmd( f"create vertex {x[i]} {y[i]} {z[i]}" )
      vertices_in_curve.append( cubit.get_last_id( "vertex" ) )
    cubit.cmd( f"create curve spline location vertex {cubit.get_id_string( vertices_in_curve )} delete" )
    curves_in_surface.append( cubit.get_last_id( "curve" ) )
  cubit.cmd( f"create surface skin curve {cubit.get_id_string( curves_in_surface )}" )
  cubit.cmd( f"delete curve {cubit.get_id_string( curves_in_surface )}" )

def mesh_surface( surface_id, mesh_size ):
  cubit.cmd( f"surface {surface_id} size {mesh_size}" )
  cubit.cmd( f"surface {surface_id} scheme map" )
  cubit.cmd( f"mesh surface {surface_id}" )

def set_performance_mode( status=True ):
  if status == True:
    cubit.cmd( "undo off" )
    cubit.cmd( "info off" )
    cubit.cmd( "echo off" )
    cubit.cmd( "warning off" )
    cubit.cmd( "journal off" )
    cubit.cmd( "graphics off" )
  elif status == False:
    cubit.cmd( "undo on" )
    cubit.cmd( "info on" )
    cubit.cmd( "echo on" )
    cubit.cmd( "warning on" )
    cubit.cmd( "journal on" )
    cubit.cmd( "graphics on" )

main( N=50, mesh_size=0.2 )

@gvernon , this thread opened by you looks super relevant to me as I very often deal with complex topolography data. I wanted to reproduce your work for my own surface points but it was not successful. I can see two lines of code (cubit.cmd( f"create curve spline location vertex {cubit.get_id_string( vertices_in_curve )} delete" ) and cubit.cmd( f"create surface skin curve {cubit.get_id_string( curves_in_surface )}" )) are doing a lot of majic but I could not use them with my data.
I have attached the csv containing points that represent the surface. I very much appreciate if anyone have time to let me know how to make a surface for my complex topography. At the end, I want to make a brick ranging from min (x) to max (x), min (y) to max (y) and 1000 to -1000 (in z direction). I want to use the surface to webcut the brick. I could only make the vertices via the following:

#!python
import numpy as np

point_data = np.genfromtxt(r"C:\Users\po7517\Desktop\Ali\Modelling\My_results\cubit_mesh\surface_generator\data\surface.csv", delimiter=',')
for point in point_data:
 cubit.cmd(f"create vertex {point[1]} {point[2]} {point[3]}")

surface.csv (198.6 KB)

@ali_geo – it seems to work for me.

If I modify my code to remove the point-generation code (since you have your own point data file), and then use a coarser mesh size I get:

#!python
import numpy
import random
import pathlib

filename = pathlib.Path( "~/surface.csv" ).expanduser()

def main( mesh_size ):
  set_performance_mode( True )
  cubit.cmd( "reset" )
  data = read_dat_file( filename )
  create_surface( data )
  surface_id = cubit.get_last_id( "surface" )
  mesh_surface( surface_id, mesh_size )
  set_performance_mode( False )

def read_dat_file( dat_filename ):
  data = {}
  x, y, z = [ [], [], [] ]
  f = open( dat_filename, "r" )
  fLines = f.readlines()
  for fLine in fLines:
    fLine = fLine.strip().split( "," )
    row_id = int( fLine[0] )
    if ( row_id in data.keys() ) == False:
      data[row_id] = { "x": [], "y": [], "z": [] }
    data[row_id]["x"].append( float( fLine[1] ) )
    data[row_id]["y"].append( float( fLine[2] ) )
    data[row_id]["z"].append( float( fLine[3] ) )
  return data

def create_surface( data ):
  row_ids = list( data.keys() )
  curves_in_surface = []
  for row_id in row_ids:
    row_data = data[row_id]
    x, y, z = [ row_data["x"], row_data["y"], row_data["z"] ]
    vertices_in_curve = []
    for i in range( 0, len( x ) ):
      cubit.cmd( f"create vertex {x[i]} {y[i]} {z[i]}" )
      vertices_in_curve.append( cubit.get_last_id( "vertex" ) )
    cubit.cmd( f"create curve spline location vertex {cubit.get_id_string( vertices_in_curve )} delete" )
    curves_in_surface.append( cubit.get_last_id( "curve" ) )
  cubit.cmd( f"create surface skin curve {cubit.get_id_string( curves_in_surface )}" )
  cubit.cmd( f"delete curve {cubit.get_id_string( curves_in_surface )}" )

def mesh_surface( surface_id, mesh_size ):
  cubit.cmd( f"surface {surface_id} size {mesh_size}" )
  cubit.cmd( f"surface {surface_id} scheme map" )
  cubit.cmd( f"mesh surface {surface_id}" )

def set_performance_mode( status=True ):
  if status == True:
    cubit.cmd( "undo off" )
    cubit.cmd( "info off" )
    cubit.cmd( "echo off" )
    cubit.cmd( "warning off" )
    cubit.cmd( "journal off" )
    cubit.cmd( "graphics off" )
  elif status == False:
    cubit.cmd( "undo on" )
    cubit.cmd( "info on" )
    cubit.cmd( "echo on" )
    cubit.cmd( "warning on" )
    cubit.cmd( "journal on" )
    cubit.cmd( "graphics on" )

main( mesh_size=50)