Source code for pyNastran.converters.stl.stl

#pylint:  disable=C0111
import copy
from struct import unpack, Struct, pack
from collections import defaultdict
from typing import Optional

import numpy as np
import scipy
from scipy.spatial import KDTree

from cpylog import get_logger2, SimpleLogger
from pyNastran.utils import is_binary_file


[docs] def read_stl(stl_filename: str, remove_elements_with_bad_normals: bool=False, log: Optional[SimpleLogger]=None, debug: bool=False): """ Reads an STL file Parameters ---------- stl_filename : str the filename to read remove_elements_with_bad_normals : bool; default=False removes elements with NAN normal Returns ------- model : STL() the stl model """ model = STL(log=log, debug=debug) model.read_stl(stl_filename) if remove_elements_with_bad_normals: model.remove_elements_with_bad_normals() return model
[docs] class STL: #model_type = 'stl' #is_structured = False #is_outward_normals = True def __init__(self, log: Optional[SimpleLogger]=None, debug: bool=False): """ Initializes the STL object Parameters ---------- debug : bool/None; default=True used to set the logger if no logger is passed in True: logs debug/info/error messages False: logs info/error messages None: logs error messages log : logging module object / None if log is set, debug is ignored and uses the settings the logging object has """ self.log = get_logger2(log, debug=debug) self.nodes = None self.elements = None self.header = '' self.infilename = None
[docs] def write_stl(self, stl_out_filename: str, is_binary: bool=False, float_fmt: str='%6.12f', normalize_normal_vectors: bool=False, stop_on_failure: bool=True) -> None: """ Writes an STL file Parameters ---------- stl_out_filename : str the filename to write is_binary : bool; default=False should a binary file be written float_fmt : str; default='%6.12f' the format to use if an ASCII file is used normalize_normal_vectors : bool; default=False should the vectors be normalized """ self.log.info(f'---writing STL...{stl_out_filename!r}---') self._validate() solid_name = 'dummy_name' if is_binary: self.write_binary_stl(stl_out_filename, normalize_normal_vectors=normalize_normal_vectors, stop_on_failure=stop_on_failure) else: self.write_stl_ascii(stl_out_filename, solid_name, float_fmt=float_fmt, normalize_normal_vectors=normalize_normal_vectors, stop_on_failure=stop_on_failure)
[docs] def read_stl(self, stl_filename: str) -> None: """ Reads an STL file Parameters ---------- stl_filename : str the filename to read """ self.infilename = stl_filename self.log.info(f'---reading STL...{self.infilename}---') if is_binary_file(stl_filename): self.read_binary_stl(stl_filename) else: self.read_ascii_stl(stl_filename)
#self.log.info("nodes=%s nelements=%s" % (self.nodes, self.nelements)) #assert self.nodes > 0, 'nodes=%s' % self.nodes #assert self.nelements > 0, 'nelements=%s' % self.nelements
[docs] def write_binary_stl(self, stl_filename: str, normalize_normal_vectors: bool=False, stop_on_failure=True) -> None: """ Write an STL binary file Parameters ---------- stl_filename : str the filename to read normalize_normal_vectors : bool; default=False should the vectors be normalized """ self._validate() with open(stl_filename, 'wb') as infile: if hasattr(self, 'header'): self.header.ljust(80, '\0') header = '%-80s' % self.header[:80] else: header = '%-80s' % stl_filename[:80] infile.write(pack(b'80s', header.encode('ascii'))) #avector = [0., 0., 0.] #bvector = [0., 0., 0.] #cvector = [0., 0., 0.] nelements = self.elements.shape[0] infile.write(pack('i', nelements)) elements = self.elements p1 = self.nodes[elements[:, 0], :] p2 = self.nodes[elements[:, 1], :] p3 = self.nodes[elements[:, 2], :] avector = p2 - p1 bvector = p3 - p1 n = np.cross(avector, bvector) del avector, bvector if normalize_normal_vectors: normals_norm = np.linalg.norm(n, axis=1) inotnan = np.where(normals_norm != 0.)[0] inan = np.where(normals_norm == 0.)[0] n[inan, :] = np.nan n[inotnan, :] /= normals_norm[inotnan, np.newaxis] s = Struct('12fH') for eid, unused_element in enumerate(elements): data = s.pack(n[eid, 0], n[eid, 1], n[eid, 2], p1[eid, 0], p1[eid, 1], p1[eid, 2], p2[eid, 0], p2[eid, 1], p2[eid, 2], p3[eid, 0], p3[eid, 1], p3[eid, 2], 0) infile.write(data)
[docs] def read_binary_stl(self, stl_filename: str) -> None: """ Read an STL binary file Parameters ---------- stl_filename : str the filename to read """ with open(stl_filename, 'rb') as infile: data = infile.read() ndata = len(data) j = 0 while j < ndata: self.log.info(f' read_binary_stl: j={j} ndata={ndata}') self.header = data[j:j+80] nelements, = unpack('i', data[j+80:j+84]) j += 84 inode = 0 nodes_dict = {} assert nelements > 0, f'nelements={nelements}' elements = np.zeros((nelements, 3), 'int32') s = Struct('12fH') for ielement in range(nelements): (unused_nx, unused_ny, unused_nz, ax, ay, az, bx, by, bz, cx, cy, cz, unused_i) = s.unpack(data[j:j+50]) t1 = (ax, ay, az) t2 = (bx, by, bz) t3 = (cx, cy, cz) if t1 in nodes_dict: i1 = nodes_dict[t1] else: i1 = inode nodes_dict[t1] = inode inode += 1 if t2 in nodes_dict: i2 = nodes_dict[t2] else: i2 = inode nodes_dict[t2] = inode inode += 1 if t3 in nodes_dict: i3 = nodes_dict[t3] else: i3 = inode nodes_dict[t3] = inode inode += 1 elements[ielement] = [i1, i2, i3] j += 50 assert inode > 0, inode nnodes = inode + 1 # accounting for indexing self.elements = elements nodes = np.zeros((nnodes, 3), 'float64') for node, inode in nodes_dict.items(): nodes[inode] = node self.nodes = nodes
def _get_normals_data(self, elements): """ This is intended as a submethod to help handle the problem of bad normals """ nodes = self.nodes #self.log.debug("get_normals...elements.shape %s" % str(elements.shape)) p1 = nodes[elements[:, 0]] p2 = nodes[elements[:, 1]] p3 = nodes[elements[:, 2]] v12 = p2 - p1 v13 = p3 - p1 v123 = np.cross(v12, v13) normals_norm = np.linalg.norm(v123, axis=1) inan = np.where(normals_norm == 0.)[0] return v123, normals_norm, inan
[docs] def remove_elements_with_bad_normals(self): """removes dot and line elements""" elements = self.elements v123, normals_norm, inan = self._get_normals_data(elements) if len(inan): inotnan = np.where(normals_norm != 0)[0] self.elements = elements[inotnan, :] normals_norm = normals_norm[inotnan] v123 = v123[inotnan] self.log.info('removing %i elements with coincident nodes' % len(inan)) normals = v123 normals[:, 0] /= normals_norm normals[:, 1] /= normals_norm normals[:, 2] /= normals_norm return normals
[docs] def get_area(self, elements, stop_on_failure: bool=True): unused_v123, normals_norm, inan = self._get_normals_data(elements) if stop_on_failure: msg = 'Failed Elements: %s\n' % inan if len(inan) > 0: for inani in inan: msg += ' eid=%s nodes=%s\n' % (inani, elements[inani, :]) for ni in elements[inani]: msg += ' nid=%s node=%s\n' % (ni, self.nodes[ni, :]) raise RuntimeError(msg) return 0.5 * normals_norm
[docs] def get_normals(self, elements, stop_on_failure: bool=True): """ Parameters ---------- elements : (n, 3) ndarray ints the elements to get the normals for nodes : (n, ) ndarray; default=None -> all a subset of the nodes stop_on_failure : bool (default=True) True: crash if there are coincident points False: delete elements """ nodes = self.nodes v123, normals_norm, inan = self._get_normals_data(elements) if stop_on_failure: msg = 'Failed Elements: %s\n' % inan if len(inan) > 0: for ifail, inani in enumerate(inan): msg += ' eid=%s nodes=%s\n' % (inani, elements[inani, :]) for ni in elements[inani]: msg += ' nid=%s node=%s\n' % (ni, nodes[ni, :]) if ifail > 10: break msg += 'Failed Elements: %s; n=%s\n' % (inan, len(inan)) raise RuntimeError(msg) # we need to divide our (n,3) array in 3 steps normals = v123 # / normals_norm normals[:, 0] /= normals_norm normals[:, 1] /= normals_norm normals[:, 2] /= normals_norm else: inotnan = np.where(normals_norm != 0.)[0] #inan = np.where(normals_norm == 0.)[0] if len(inan): #normals_norm[inan] = np.array([1., 0., 0.]) normals_norm[inan] = 1. #normals_norm[inan, [1,2]] = 0. #elements = elements[inotnan, :] #normals_norm = normals_norm[inotnan] #v123 = v123[inotnan] # we need to divide our (n,3) array in 3 steps if 0: normals = v123 # / normals_norm normals[:, 0] /= normals_norm normals[:, 1] /= normals_norm normals[:, 2] /= normals_norm normals[inan, :] = -1.01 else: normals = v123 # / normals_norm normals[inotnan, 0] /= normals_norm[inotnan] normals[inotnan, 1] /= normals_norm[inotnan] normals[inotnan, 2] /= normals_norm[inotnan] return normals
[docs] def flip_normals(self, i=None) -> None: """ Flips the normals of the specified elements. Parameters ---------- i : (n, ) ndarray ints; default=None -> all the indicies to flip """ self.log.info("---flip_normals---") if i is None: elements = self.elements else: elements = self.elements[i, :] n0, n1, n2 = elements[:, 0], elements[:, 1], elements[:, 2] elements2 = elements.copy() elements2[:, 0] = n0 elements2[:, 1] = n2 elements2[:, 2] = n1 if i is None: self.elements = elements2 else: self.elements[i, :] = elements2 #[i, :]
[docs] def get_normals_at_nodes(self, normals=None, nid_to_eid=None): """ Calculates the normal vector of the nodes based on the average element normal. Parameters ---------- normals : (n, 3) ndarray floats The elemental normals nid_to_eid : dict[int] = [int, int, ... ] key = node_id value = list of element_ids Returns ------- normals_at_nodes : (nnodes, 3) ndarray ints the normals """ elements = self.elements nodes = self.nodes if normals is None: normals = self.get_normals(elements) if nid_to_eid is None: nid_to_eid = defaultdict(list) eid = 0 for (n1, n2, n3) in elements: nid_to_eid[n1].append(eid) nid_to_eid[n2].append(eid) nid_to_eid[n3].append(eid) eid += 1 del eid, n1, n2, n3 normals_at_nodes = np.zeros(nodes.shape, 'float64') eid = 0 for nid, elementsi in nid_to_eid.items(): pe = normals[elementsi] m = pe.mean(axis=0) normals_at_nodes[nid] = m / np.linalg.norm(m) eid += 1 return normals_at_nodes
[docs] def equivalence_nodes(self, tol: float=1e-5) -> None: """equivalences the nodes of the model and updates the elements""" nnodes = self.nodes.shape[0] # build the kdtree kdt = KDTree(self.nodes) # find the node ids of interest nids_new = np.unique(self.elements.ravel()) nids_new.sort() # check the closest 10 nodes for equality unused_deq, ieq = kdt.query(self.nodes[nids_new, :], k=10, distance_upper_bound=tol) # get the ids of the duplicate nodes slots = np.where(ieq[:, 1:] < nnodes) replacer = np.unique(ieq[slots]) # update the duplcated node id with it's partner # we'll pick the minimum ID for r in replacer: ip = np.where(ieq[r, :] < nnodes)[0] possible = ieq[r, ip] # node 11 can become node 10, but node 10 cannot become node 11 ip2 = np.where(r > possible)[0] if len(ip2): # replace the node ids possible2 = possible[ip2] r_new_nid = possible2.min() ireplace = np.where(self.elements == r) self.elements[ireplace] = r_new_nid
def _validate(self) -> None: assert len(self.nodes) > 0, 'No nodes were found in the model' assert len(self.elements) > 0, 'No nodes were found in the model'
[docs] def write_stl_ascii(self, out_filename: str, solid_name: str, float_fmt: str='%.6f', normalize_normal_vectors: bool=False, stop_on_failure: bool=True) -> None: """ Writes an STL in ASCII format solid solid_name facet normal -6.601157e-001 6.730213e-001 3.336009e-001 outer loop vertex 8.232952e-002 2.722531e-001 1.190414e+001 vertex 8.279775e-002 2.717848e-001 1.190598e+001 vertex 8.557653e-002 2.745033e-001 1.190598e+001 endloop endfacet end solid """ self.log.info("---write_stl_ascii...%r---" % out_filename) self._validate() noormal_format = ' facet normal %s %s %s\n' % (float_fmt, float_fmt, float_fmt) vertex_format = ' vertex %s %s %s\n' % (float_fmt, float_fmt, float_fmt) msg = 'solid %s\n' % solid_name normals = self.get_normals(self.elements, stop_on_failure=stop_on_failure) nodes = self.nodes elements = self.elements with open(out_filename, 'w') as out: out.write(msg) for element, normal in zip(elements, normals): try: p1, p2, p3 = nodes[element] except IndexError: print(element) raise #msg += ' facet normal -6.601157e-001 6.730213e-001 3.336009e-001\n' msg = noormal_format % tuple(normal) msg += ' outer loop\n' msg += vertex_format % tuple(p1) msg += vertex_format % tuple(p2) msg += vertex_format % tuple(p3) msg += ' endloop\n' msg += ' endfacet\n' #print(msg) out.write(msg) msg = 'endsolid\n' out.write(msg)
[docs] def read_ascii_stl(self, stl_filename: str) -> None: """ Reads an STL that's in ASCII format """ with open(stl_filename, 'r') as infile: line = infile.readline() inode = 0 ielement = 0 nodes_dict = {} elements = [] while line: if 'solid' in line[:6].lower(): line = infile.readline() # facet while 'facet' in line.strip()[:5].lower(): #facet normal -6.665299e-001 6.795624e-001 3.064844e-001 # outer loop # vertex 8.142845e-002 2.731541e-001 1.190024e+001 # vertex 8.186898e-002 2.727136e-001 1.190215e+001 # vertex 8.467505e-002 2.754588e-001 1.190215e+001 # endloop #endfacet infile.readline() # outer loop L1 = infile.readline().strip() L2 = infile.readline().strip() L3 = infile.readline().strip() v1 = L1.split()[1:] v2 = L2.split()[1:] v3 = L3.split()[1:] infile.readline() # endloop infile.readline() # endfacet t1 = tuple(v1) t2 = tuple(v2) t3 = tuple(v3) assert len(v1) == 3, '%r' % L1 assert len(v2) == 3, '%r' % L2 assert len(v3) == 3, '%r' % L3 if t1 in nodes_dict: i1 = nodes_dict[t1] else: i1 = inode nodes_dict[t1] = inode inode += 1 if t2 in nodes_dict: i2 = nodes_dict[t2] else: i2 = inode nodes_dict[t2] = inode inode += 1 if t3 in nodes_dict: i3 = nodes_dict[t3] else: i3 = inode nodes_dict[t3] = inode inode += 1 element = [i1, i2, i3] elements.append(element) ielement += 1 line = infile.readline() # facet #print "end of solid..." elif 'endsolid' in line.lower(): line = infile.readline() elif line.strip() == '': line = infile.readline() else: self.log.error(line) #line = f.readline() raise NotImplementedError(f'multiple solids are not supported; line={line!r}') #break assert inode > 0, inode nnodes = inode + 1 # accounting for indexing self.elements = np.array(elements, 'int32') nodes = np.zeros((nnodes, 3), 'float64') for node, inode in nodes_dict.items(): nodes[inode] = node self.nodes = nodes
[docs] def scale_nodes(self, xscale, yscale=None, zscale=None): """ Scales the model Parameters ---------- xscale : float the scaling factor for the x axis; also the default scaling factor yscale/zscale : float; default=xscale the scaling factors for the y/z axes """ if yscale is None: yscale = xscale if zscale is None: zscale = xscale self.nodes[:, 0] *= xscale self.nodes[:, 1] *= yscale self.nodes[:, 2] *= zscale
[docs] def shift_nodes(self, xshift, yshift, zshift): """Shifts the model""" self.nodes[:, 0] += xshift self.nodes[:, 1] += yshift self.nodes[:, 2] += zshift
[docs] def flip_axes(self, axes, scale): """ Swaps the axes Parameters ---------- axes : str 'xy', 'yz', 'xz' scale : float why is this here, but is not applied to all axes? """ if axes == 'xy': x = copy.deepcopy(self.nodes[:, 0]) y = copy.deepcopy(self.nodes[:, 1]) self.nodes[:, 0] = y * scale self.nodes[:, 1] = x * scale elif axes == 'yz': y = copy.deepcopy(self.nodes[:, 1]) z = copy.deepcopy(self.nodes[:, 2]) self.nodes[:, 1] = z * scale self.nodes[:, 2] = y * scale elif axes == 'xz': x = copy.deepcopy(self.nodes[:, 0]) z = copy.deepcopy(self.nodes[:, 2]) self.nodes[:, 0] = z * scale self.nodes[:, 2] = x * scale
[docs] def create_mirror_model(self, xyz, tol: float) -> None: """ Creates a mirror model. Parameters ---------- xyz : str {x, y, z} the direction of symmetry tol: float the tolerance for symmetry plane nodes .. note:: All elements on the symmetry plane will be removed """ assert xyz in ['x', 'y', 'z'], 'xyz=%r' % xyz assert tol >= 0.0, 'tol=%s' % tol nnodes = self.nodes.shape[0] if xyz == 'x': xyzi = 0 elif xyz == 'y': xyzi = 1 elif xyz == 'z': xyzi = 2 else: raise RuntimeError(xyz) # the nodes on the symmetry plane i = np.where(self.nodes[:, xyzi] < tol)[0] # smash the symmetry nodes to 0.0 self.nodes[i, xyzi] = 0. nodes_sym = copy.deepcopy(self.nodes) nodes_sym[:, xyzi] *= -1. # we're lazy and duplicating all the nodes # but will only write out a subset of them nodes = np.vstack([self.nodes, nodes_sym]) # create the symmetrical elements elements2 = [] elements3 = [] for element in self.elements: # the 3 "y" locations for the element epoints = nodes[element, xyzi] # [0] je = np.where(epoints <= tol)[0] if len(je) < 3: # not a symmetry element, so we save it elements2.append(element) # duplicate the node if it's not on the symmetry plane element3 = [elementi if elementi in i else (elementi + nnodes) for elementi in element] # the normal is now backwards, so we flip it element3.reverse() elements3.append(element3) self.nodes = nodes self.elements = np.array(elements2 + elements3, dtype='int32')
def _rotate_model(stl: STL) -> None: # pragma: no cover nodes = stl.nodes elements = stl.elements if 0: # rotate the model x, y, z = nodes[:, 0], nodes[:, 1], nodes[:, 2] #i = where(y > 0.0)[0] R = x**2 + y**2 theta = np.arctan2(y, x) iRz = np.where(R == 0)[0] theta[iRz] = 0.0 min_theta = min(theta) unused_dtheta = max(theta) - np.pi / 4 theta2 = theta + min_theta x2 = R * np.cos(theta2) y2 = R * np.sin(theta2) #print("x.shape", x.shape, y2.shape) nodes_rotated = np.transpose(np.vstack([x2, y2, z])) #print("nodes.shape", nodes_rotated.shape) #print(nodes_rotated) if 0: # project the volume (unused_nodes2, unused_elements2) = stl.project_mesh(nodes_rotated, elements) # write the model stl_geom_out = 'rotated.stl' stl.write_stl_ascii(stl_geom_out, 'sphere') if __name__ == '__main__': # pragma: no cover from pyNastran.converters.stl.stl_reshape import main main()