diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py
index 237fd99b..c7e6a4ac 100644
--- a/python/ctem/ctem/ctem_viewer_3d.py
+++ b/python/ctem/ctem/ctem_viewer_3d.py
@@ -1,48 +1,16 @@
 import numpy as np
-import open3d as o3d
-import trimesh
-
-class surf:
+import open3d
+class polysurface:
     """A model of a self-gravitating small body"""
     def __init__(self):
         self.polyfile = ""
-        self.mesh = trimesh.Trimesh()
         #used for Open3d module
-        self.o3d_mesh = o3d.geometry.TriangleMesh()
-        self.nface = 0
-        self.nvert = 0
-        return
-    
-    def io_read_ply(self,filename):
-        """Reads in PLY file into Open3d triangular mesh object and computes normals"""
-        self.polyfile = filename
-        self.mesh = trimesh.load(f"{self.polyfile}")
-        trimesh.repair.fix_normals(self.mesh)
-        self.mesh.vertex_colors = np.full((self.nvert,4),1.0)
-        self.nvert = len(self.mesh.vertices)
-        self.nface = len(self.mesh.faces)
-        self.convert_trimesh_to_o3d()
+        self.mesh = open3d.geometry.TriangleMesh()
         return
     
-    def render(self):
-        """Renders output"""
-        self.convert_trimesh_to_o3d()
-        o3d.visualization.draw_geometries([self.o3d_mesh])
-        return
-    
-    def convert_trimesh_to_o3d(self):
-        """Update Open3D mesh with new vertices and faces, and recompute all normals"""
-        self.o3d_mesh.vertices  = o3d.utility.Vector3dVector(self.mesh.vertices)
-        self.o3d_mesh.triangles = o3d.utility.Vector3iVector(self.mesh.faces)
-        self.o3d_mesh.compute_vertex_normals()
-        self.o3d_mesh.compute_triangle_normals()
-        self.o3d_mesh.vertex_colors = o3d.utility.Vector3dVector(self.mesh.vertex_colors[:,0:3])
-        return
-
 if __name__ == '__main__':
     import os
     import ctem_io_readers
-    import open3d
 
     #Create and initialize data dictionaries for parameters and options from CTEM.in
     notset = '-NOTSET-'
@@ -84,34 +52,32 @@ def convert_trimesh_to_o3d(self):
     dem_file = parameters['workingdir'] + 'surface_dem.dat'
     surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize'])
     # Build mesh grid
-    gridsize = parameters['gridsize']
+    s = parameters['gridsize']
     pix = parameters['pix']
-    ygrid, xgrid = np.mgrid[-gridsize / 2:gridsize / 2,
-                             -gridsize / 2:gridsize / 2] * pix
+    ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix
 
     xvals = xgrid.flatten()
     yvals = ygrid.flatten()
     zvals = surface_dem.flatten()
-    dem_vec = np.array((xvals, yvals, zvals*10)).T
-    pcd = open3d.geometry.PointCloud()
-    pcd.points = open3d.utility.Vector3dVector(dem_vec)
-    pcd.estimate_normals()
-
-    # estimate radius for rolling ball
-    distances = pcd.compute_nearest_neighbor_distance()
-    avg_dist = np.mean(distances)
-    radius = 1.5 * avg_dist
-
-    mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(
-        pcd,
-        o3d.utility.DoubleVector([radius, radius * 2]))
-
-    # create the triangular mesh with the vertices and faces from open3d
-    tri_mesh = trimesh.Trimesh(np.asarray(mesh.vertices), np.asarray(mesh.triangles),
-                               vertex_normals=np.asarray(mesh.vertex_normals))
-
-    trimesh.convex.is_convex(tri_mesh)
-
-    o3d.visualization.draw_geometries([trimesh])
-    
+    verts = np.array((xvals, yvals, zvals)).T
+    # pcd = open3d.geometry.PointCloud()
+    # pcd.points = open3d.utility.Vector3dVector(verts)
+    # pcd.estimate_normals()
+    # pcd.normalize_normals()
+    # o3d.visualization.draw_geometries([pcd], point_show_normal=True)
+    # tetra_mesh, tetra_pt_map = open3d.geometry.TetraMesh.create_from_point_cloud(pcd)
+    # o3d.visualization.draw_geometries([tetra_mesh])
+    faces = np.empty([2*(s-1)**2, 3], dtype=np.int64)
+    for j in range(s - 1):
+        for i in range(s - 1):
+            idx = (s-1)*j+i
+            #faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ]
+            faces[idx,:] = [ j*s+i, j*s+i+1,  (j+1)*s+i ]
+            idx += (s-1)**2
+            faces[idx,:] = [j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ]
+    surf = polysurface()
+    surf.mesh.vertices  = open3d.utility.Vector3dVector(verts)
+    surf.mesh.triangles = open3d.utility.Vector3iVector(faces)
+    surf.mesh.compute_vertex_normals()
+    open3d.visualization.draw_geometries([surf.mesh])
 
diff --git a/python/ctem/tests/viewer3d/polytest.py b/python/ctem/tests/viewer3d/polytest.py
index 451f8629..13cee521 100644
--- a/python/ctem/tests/viewer3d/polytest.py
+++ b/python/ctem/tests/viewer3d/polytest.py
@@ -3,6 +3,6 @@
 
 surf = ctem_viewer_3d.surf()
 surf.io_read_ply("216kleopatra-mesh.ply")
-v = np.asarray(surf.o3d_mesh)
-t = np.asarray(surf.o3d_mesh)
+v = np.asarray(surf.mesh)
+t = np.asarray(surf.mesh)
 surf.render()
\ No newline at end of file