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

[gvernon]

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 )