compare_facet_techniques.py

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 )