voxelise solidworks data with full tensors

.gitignore

 *.pyc
+*.DS_Store
 *.csv
 *.gcode
 *.json
+*.stl

main.jl

@@ -460,7 +460,7 @@ end
 T0 = 215 # extrusion temp
 Tc = 25 # room temp
 Tcutoff = 100 # temperature above which strain isn't happening
-k = 0.005 #value unknown  - assumes cooling to 50C from 215 in 8 minutes
+k = 0.02 #value unknown
 Temp(t::Number) = Tc + (T0-Tc)*ℯ^(-k*t)
 # see shape of temp function
 #x = 1:100
@@ -481,6 +481,7 @@ function calc_cost(rollout::Vector{Int}, cdata::contourdata, vd::voxdata)
     considered_voxels = (1:length(vd.below))[(below.!=0) .& (.!isnan.(voxtimes))] # cannot calculate cost is no voxel underneath
     considered_voxels = considered_voxels[.!isnan.(voxtimes[below[considered_voxels]])]
     Δt = voxtimes[considered_voxels] - voxtimes[below[considered_voxels]]
+    ΔT = Tcutoff .- Temp.(Δt)
 
     # calculate stresses at each voxel
 

utils/__init__.py

+from .model import *

utils/model.py

+import numpy as np
+import trimesh as tm
+
+
+def get_model_height_arr(model_mesh, layer_height):
+    """
+    Get the z-height increments for the bounds of the model
+    Dependent on layer height
+    """
+    model_min_height, model_max_height = get_model_min_max_height(model_mesh)
+    model_height = get_model_height(model_mesh)
+    height_step_num = int(model_height / layer_height)
+    heights = np.round(np.linspace(model_min_height+layer_height, model_max_height-layer_height, height_step_num), 3)
+    return heights
+
+
+def get_model_height(model_mesh):
+    """
+    Returns the height of the model
+    """
+    model_min_height, model_max_height = get_model_min_max_height(model_mesh)
+    model_height = model_max_height - model_min_height
+    return model_height
+
+
+def get_model_min_max_height(model_mesh):
+    """
+    Get the minimum and maximum z values of the mesh
+    """
+    model_min_height = get_model_bounds(model_mesh)[4]
+    model_max_height = get_model_bounds(model_mesh)[5]
+    return model_min_height, model_max_height
+
+
+def get_model_xy_grid(model_mesh, resolution):
+    """
+    Get the xy grid used to generate the surfaces for slicing
+    """
+    bounds = get_model_bounds(model_mesh)
+    min_x, max_x = bounds[0], bounds[1]
+    min_y, max_y = bounds[2], bounds[3]
+    x = np.linspace(min_x-1e-3, max_x+1e-3, resolution)
+    y = np.linspace(min_y-1e-3, max_y+1e-3, resolution)
+    x, y = np.meshgrid(x, y)
+    return x, y
+
+
+def get_model_bounds(model_mesh):
+    """
+    Get the bounds of the given model (xmin, xmax, ymin, ymax, zmin, zmax)
+    """
+    return model_mesh.bounds
+
+
+def transform_to_point(model_mesh, point):
+    """
+    Generates the transform required to move from mesh centre to given point
+    """
+    mesh_centre = model_mesh.center
+    delta = np.array(point) - np.array(mesh_centre)
+    transform = tm.transformations.translation_matrix(delta)
+    model_mesh.transform(transform)
+
+
+def transform_to_bed_centre(model_mesh, bed_centre):
+    """
+    Generate the transform to move the model mesh to the centre of the bed
+    """
+    transform_to_point(model_mesh, bed_centre)
+
+
+def transform_to_z0(model_mesh):
+    """
+    Generate the transform required to put the bottom of the model on the build plate
+    """
+    model_centre = model_mesh.center
+    model_height = get_model_height(model_mesh)
+    new_loc = model_centre
+    new_loc[2] = model_height/2
+    transform_to_point(model_mesh, new_loc)
+
+
+def get_centroidal_axis_surf(model_mesh):
+    """
+    Slices a model at it's vertical centroid and returns the x-y grid defining the slice
+    """
+    cen_slice = model_mesh.slice(normal='z', origin=[0, 0, get_model_height(model_mesh)/2])
+    return cen_slice

voxelise.py

@@ -9,6 +9,7 @@ from scipy.spatial import KDTree
 from scipy.interpolate import LinearNDInterpolator
 import pandas as pd
 import json
+import utils as tu
 
 
 class Voxelizer:
@@ -56,11 +57,11 @@ class Voxelizer:
         print(gs-c)
 
         print("Interpolating pointcloud to grid stresses")
-        self.interpolator = LinearNDInterpolator(self.mesh.points, self.mesh['StressValues'])
+        #self.interpolator = LinearNDInterpolator(self.mesh.points, self.mesh['StressValues'])
         print(time.time()-gs) # up to 100s for mesh with 300k points
 
         print("Performing interpolations")
-        self.cstress = self.interpolator(self.centers)
+        #self.cstress = self.interpolator(self.centers)
 
 
     def check_interp1s(self, node_is):
@@ -90,16 +91,83 @@ class Voxelizer:
         return csv_name
 
 
+def voxelised_csv(mesh_file_name,
+            csv_fname,
+            voxels=None,
+            skip_lines=4,
+            voxel_dim=3,
+            layer_height=1
+        ):
+    """
+    takes in input from solidworks csv and voxelises.
+    """
+
+    s = time.time()
+
+    print("Reading mesh")
+    mesh = pv.read(mesh_file_name)
+    bed_centre = [110, 110, 0]
+    tu.transform_to_bed_centre(mesh, bed_centre)
+    tu.transform_to_z0(mesh)
+    #boundary = mesh.decimate_boundary()
+    v = time.time()
+    print(v-s)
+
+    print("Voxelizing mesh")
+    hwratio = voxel_dim/layer_height
+    if voxels is None:
+        m = mesh.copy()
+        m.points[:,2] *= hwratio
+        grid = pv.voxelize(m, density=voxel_dim, check_surface=False)
+        grid.points[:,2] *= 1/hwratio
+
+        print("Getting Voxel centers")
+        #points_to_add = []
+        temptree = KDTree(grid.points, leafsize=5) #kd tree for fast lookup of adjacent grid points
+        # TODO This logic is not necessarily accurate - voxel could still be empty
+        # but should be fine for our case
+        d1, ii  = temptree.query(grid.points + [voxel_dim, voxel_dim, layer_height])
+        d2, ii  = temptree.query(grid.points + [voxel_dim, 0, layer_height])
+        d3, ii  = temptree.query(grid.points + [0, voxel_dim, 0])
+        centers = grid.points[d1 + d2 + d3 < 1e-4] + [voxel_dim/2, voxel_dim/2, layer_height/2]
+        c = time.time()
+        print(c-s)
+        v2 = pd.DataFrame(centers, columns=['x','y','z'])
+        v2.to_csv('temp_vox_centers.csv', index=False)
+    else:
+        centers= pd.read_csv(voxels).to_numpy()
+
+    data = pd.read_csv(csv_fname, header=skip_lines)
+    locs = data.iloc[:,1:4].to_numpy()
+    locs_pv = pv.PolyData(locs)
+    tu.transform_to_bed_centre(locs_pv, bed_centre)
+    tu.transform_to_z0(locs_pv)
+    locs = locs_pv.points
+    ys = data.iloc[:,4:].to_numpy()
+
+    interpolator = LinearNDInterpolator(locs, ys)
+    stress = interpolator(centers)
+
+    # output values
+    out = pd.DataFrame(np.hstack((centers, stress)), columns=['x', 'y', 'z', 'Sx', 'Sy', 'Sz', 'Txy', 'Tyz', 'Txz'])
+    out = out.dropna()
+    return out
+
+
 def contour2dict(c):
     return {'pos':c.pos, 'time':c.time}
 
+
 if __name__ == '__main__':
-    obj = 'tensile'
+    obj = 'M1'
+    out = voxelised_csv(obj + '.stl', obj + '_raw.csv', voxels=obj+'_vox_centers.csv')
+    out.to_csv(obj + 'voxels.csv')
 
     #pl = Voxelizer(obj + '.vtk', 1, 1) #0.25)
-    #pl.export_voxels('tensile')
+    #pl.export_voxels(obj)
+
+    #contours = gc.decode_gcode(obj + '.gcode')
+    #outfile = open(obj + 'contours.json','w')
+    #json.dump([contour2dict(c) for c in contours], outfile)
+    #outfile.close()
 
-    contours = gc.decode_gcode(obj + '.gcode')
-    outfile = open(obj + 'contours.json','w')
-    json.dump([contour2dict(c) for c in contours], outfile)
-    outfile.close()