Skip to content

Instantly share code, notes, and snippets.

@GregVernon
Created November 12, 2023 06:39
Show Gist options
  • Save GregVernon/e063f9430d4f9288089a8428594f90a1 to your computer and use it in GitHub Desktop.
Save GregVernon/e063f9430d4f9288089a8428594f90a1 to your computer and use it in GitHub Desktop.
import numpy
import sys
from matplotlib import pyplot as plt
sys.path.append( r"C:\Program Files\Coreform Cubit 2023.11\bin")
import cubit
def main( angles ):
fea_results, coarse_results = do_study( angles )
plot_results( angles, fea_results, coarse_results )
def plot_results( angles, fea_results, coarse_results ):
# Tris vs Angle
fig, ax = plt.subplots()
ax.plot( fea_results["num_tris"], angles, "-o", label="regular" )
ax.plot( coarse_results["num_tris"], angles, "-o", label="coarse" )
ax.set_xlabel( "Number of Triangles" )
ax.set_ylabel( "Requested Angle" )
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid( visible=True, which="both", axis="both" )
ax.legend()
plt.show()
# Tris vs Tri Surf Err
fig, ax = plt.subplots()
ax.plot( fea_results["num_tris"], fea_results["tri_surf_dist_err"], "-o", label="regular" )
ax.plot( coarse_results["num_tris"], coarse_results["tri_surf_dist_err"], "-o", label="coarse" )
ax.set_xlabel( "Number of Triangles" )
ax.set_ylabel( "Maximum Triangle-Surface Distance" )
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid( visible=True, which="both", axis="both" )
ax.legend()
plt.show()
# Tris vs Tri Vol Err
fig, ax = plt.subplots()
ax.plot( fea_results["num_tris"], fea_results["tri_vol_err"], "-o", label="regular" )
ax.plot( coarse_results["num_tris"], fea_results["tri_vol_err"], "-o", label="coarse" )
ax.set_xlabel( "Number of Triangles" )
ax.set_ylabel( "Triangulated Volume Error" )
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid( visible=True, which="both", axis="both" )
ax.legend()
plt.show()
def do_study( angles ):
num_tris = []
tri_surf_dist_err = []
tri_vol_err = []
for i in range( 0, len( angles ) ):
num_tris.append( make_cylinder_mesh( "fea", angles[i] ) )
tri_surf_dist_err.append( compute_tri_surf_dist_err() )
tri_vol_err.append( compute_tri_vol_err() )
fea_results = { "num_tris": num_tris, "tri_surf_dist_err": tri_surf_dist_err, "tri_vol_err": tri_vol_err }
num_tris = []
tri_surf_dist_err = []
tri_vol_err = []
for i in range( 0, len( angles ) ):
num_tris.append( make_cylinder_mesh( "coarse", angles[i] ) )
tri_surf_dist_err.append( compute_tri_surf_dist_err() )
tri_vol_err.append( compute_tri_vol_err() )
coarse_results = { "num_tris": num_tris, "tri_surf_dist_err": tri_surf_dist_err, "tri_vol_err": tri_vol_err }
return fea_results, coarse_results
def make_cylinder_mesh( method, angle ):
cubit.cmd( "reset" )
cubit.cmd( "cylinder radius 1.0 height 4" )
approx_size = 10 * ( 2 * 1.0 * numpy.tan( numpy.deg2rad( angle ) / 2.0 ) )
cubit.cmd( f"surface all size {approx_size}" )
match method:
case "fea":
cubit.cmd( f"surface all scheme trimesh geometry approximation angle {angle}" )
cubit.cmd( "set trimesher geometry sizing on" )
cubit.cmd( "set trimesher coarse off" )
case "coarse":
cubit.cmd( "surface all scheme trimesh" )
cubit.cmd( f"set trimesher coarse on ratio 100 angle {angle}" )
cubit.cmd( "mesh surface all" )
num_tris = len( cubit.get_entities( "tri" ) )
return num_tris
def compute_tri_surf_dist_err():
T = cubit.get_entities( "tri" )
S = cubit.get_entities( "surface" )
tri_surf_dist = numpy.zeros( len( T ), dtype="double" )
for sid in S:
surface = cubit.surface( sid )
T = cubit.parse_cubit_list( "tri", f"in surface {sid}" )
for tid in T:
tri_center = cubit.get_center_point( "tri", tid )
surf_loc = surface.closest_point_trimmed( tri_center )
tri_surf_dist[tid-1] = numpy.sqrt( (tri_center[0] - surf_loc[0])**2 + (tri_center[1] - surf_loc[1])**2 + (tri_center[2] - surf_loc[2])**2 )
max_tri_surf_dist = max( tri_surf_dist )
return max_tri_surf_dist
def compute_tri_vol_err():
cubit.cmd( "vol all scheme tetmesh" )
cubit.cmd( "mesh volume all" )
V = cubit.get_entities( "volume" )
cad_vol = numpy.zeros( len( V ), dtype="double" )
mesh_vol_err = numpy.zeros( len( V ), dtype="double" )
for v in range( 0, len( V ) ):
vid = V[v]
cad_vol[v] = cubit.volume( vid ).volume()
cubit.cmd( "volume all scheme tetmesh" )
cubit.cmd( "mesh volume all" )
mesh_vol = cubit.get_meshed_volume_or_area( "volume", (vid,) )
mesh_vol_err[v] = abs( mesh_vol - cad_vol ) / cad_vol
return mesh_vol_err
if __name__ == "__main__":
dev_angles = numpy.logspace( -1, 1, 11 )
main( dev_angles )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment