Source code for pyNastran.converters.stl.dev.stl_mesh

"""
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.479.8237&rep=rep1&type=pdf
http://math.stackexchange.com/questions/4322/check-whether-a-point-is-within-a-3d-triangle
https://en.wikipedia.org/wiki/Barycentric_coordinate_system
http://en.wikipedia.org/wiki/Line-plane_intersection
http://math.stackexchange.com/questions/544946/determine-if-projection-of-3d-point-onto-plane-is-within-a-triangle
"""
import numpy as np

from pyNastran.converters.stl.stl import read_stl
from scipy.spatial import cKDTree
import scipy.interpolate

[docs]def projected_barycentric_coord(p, q, u, v): r""" points p, q vector u, v 3 *v / \ <----p q*----*u 1 2 u = p2 - p1 v = p3 - p1 """ n = np.cross(u, v) one_over_4_area_squared = 1.0 / (n @ n) w = p - q b[2] = (np.cross(u, w) @ n) * one_over_4_area_squared b[1] = (np.cross(w, v) @ n) * one_over_4_area_squared b[0] = 1.0 - b[1] - b[2] return b
[docs]def project_points_onto_stl(stl, points): """ Parameters ---------- nodes : (n, 3) ndarray floats The nodes on the surface. elements : (n, 3) ndarray ints The elements on the surface. """ nodes = stl.nodes elements = stl.elements if not hasattr(stl, 'centroids'): n1 = elements[:, 0] n2 = elements[:, 1] n3 = elements[:, 2] p1 = nodes[n1, :] p2 = nodes[n2, :] p3 = nodes[n3, :] centroids = (p1 + p2 + p3) / 3. stl.centroids = centroids tree = cKDTree(centroids, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) stl.tree = tree #tree = scipy.spatial.KDTree(data, leafsize=10) #tree.query_ball_point(x, r, p=2., eps=0) #tree.query_ball_tree(other, r, p=2., eps=0) #tree.query_pairs(r, p=2., eps=0) #tree.sparse_distance_matrix(other, max_distance, p=2.) tree = stl.tree #d : array of floats #The distances to the nearest neighbors. #If x has shape tuple+(self.m,), then d has shape tuple+(k,). #Missing neighbors are indicated with infinite distances. #i : ndarray of ints #The locations of the neighbors in self.data. #If `x` has shape tuple+(self.m,), then `i` has shape tuple+(k,). #Missing neighbors are indicated with self.n. dist, i = tree.query(points, k=1, eps=0, p=2, distance_upper_bound=np.inf, n_jobs=1) # distance from centroid to point, such that we get the element id directly print('dist =', dist) print('i =', i) n1 = elements[i, 0] n2 = elements[i, 1] n3 = elements[i, 2] p1 = nodes[n1, :] p2 = nodes[n2, :] p3 = nodes[n3, :] u = p2 - p1 v = p3 - p1 #w = points_rotated - p1 n = np.cross(u, v) #n2 = 1 / n**2 #gamma_a = (np.cross(u, w) @ n) / n2 #gamma_b = (np.cross(u, w) @ n) / n2 try: nmag = np.linalg.norm(n, axis=1) except ValueError: print('n.shape =', n.shape) raise #area = nmag / 2. assert nmag.size == n1.size, 'nmag.size=%s n1.size=%s' % (nmag.size, n1.size) print('n1 =', n1) print('n2 =', n2) print('n3 =', n3) p = points a = p1 b = p2 c = p3 # http://math.stackexchange.com/questions/544946/determine-if-projection-of-3d-point-onto-plane-is-within-a-triangle pc = p - c alpha = np.linalg.norm(np.cross(p - b, pc)) / nmag beta = np.linalg.norm(np.cross(pc, p - a)) / nmag gamma = 1 - alpha - beta #print('alpha =', alpha) #print('beta =', beta) #print('gamma =', gamma) #print('a*p =', alpha[:, np.newaxis] * p1) p_prime = alpha[:, np.newaxis] * p1 + beta[:, np.newaxis] * p2 + gamma[:, np.newaxis] * p3 #print('p_prime =\n', p_prime) #tree.query_ball_point(x, r, p=2., eps=0) #tree.query_ball_tree(other, r, p=2., eps=0) #tree.query_pairs(r, p=2., eps=0) #tree.sparse_distance_matrix(other, max_distance, p=2.) return p_prime
[docs]def project_line_onto_stl(stl, pa, pb, npoints=11): """top down projection""" normal = np.array([0., 0., -1.], dtype='float32') #max_z = nodes[:, 2].max() #min_z = nodes[:, 2].min() # TODO: rotate if want a new normal #dz = max_z - min_z #dzi = dz / 20. #points_rotated = points #out_points = project_points_onto_stl(stl, points) # TODO: rotate if want a new normal p = np.linspace(0., 1., num=npoints, endpoint=True) p21 = pb - pa ratio = p21 / np.linalg.norm(p21) print('p =', p) print('ratio =', ratio) points = pa + p[:, np.newaxis] * ratio print('points =', points) out_points = project_points_onto_stl(stl, points) return out_points
[docs]def project_curve_onto_stl(stl, points, npoints=11): """top down projection""" normal = np.array([0., 0., -1.], dtype='float32') #max_z = nodes[:, 2].max() #min_z = nodes[:, 2].min() # TODO: rotate if want a new normal #dz = max_z - min_z #dzi = dz / 20. #points_rotated = points #out_points = project_points_onto_stl(stl, points) # TODO: rotate if want a new normal # create interpolation curve from points p2 = points[1:, :] p1 = points[:-1, :] dx = np.linalg.norm(p2 - p1, axis=1) assert dx.size == p1.shape[0] t = dx.sum() pa = points[0, :] dp = points - pa dx2 = np.linalg.norm(dp, axis=1) t = dx2.sum() # http://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.interpolate.interp1d.html func = scipy.interpolate.interp1d(t, dx2, kind='cubic', axis=-1, copy=True, bounds_error=None, fill_value=np.nan, assume_sorted=False) # cubic spline p = np.linspace(0., t, num=npoints, endpoint=True) t2 = func(p) dx = func(t2) + pa #p21 = pb - pa #ratio = p21 / np.linalg.norm(p21) #print('p =', p) #print('ratio =', ratio) points = pa + dx print('points =', points) out_points = project_points_onto_stl(stl, points) return out_points
[docs]def main(): # pragma: no cover import os import pyNastran PKG_PATH = pyNastran.__path__[0] stl_filename = os.path.join(PKG_PATH, 'converters', 'stl', 'sphere.stl') stl = read_stl(stl_filename) #XYZ Global = (2.0035907914418716, 1.3287668328026303, 2.873731014735773) #NodeID = 142; xyz=(1.88823, 1.5, 2.94889) #lineNo=2110 annotate_cell_picker() #XYZ Global = (1.9419959964242275, 1.141259948469464, 2.869267723165781) #NodeID = 141; xyz=(1.93018, 1.02165, 2.85504) #lineNo=2110 annotate_cell_picker() #XYZ Global = (2.1320656653448338, 1.4367816967143772, 2.83778333777658) #NodeID = 137; xyz=(2.25, 1.5, 2.79904) # nids = [142, 137, 141] # 2.0035907914418716, 1.3287668328026303, 2.873731014735773 points = np.array([ [2.0035907914418716, 1.3287668328026303, 2.873731014735773], [2.25, 1.5, 2.79904], [2.25, 1.5, 2.79903], ], dtype='float32') pa = points[0, :] pb = points[1, :] out_points = project_points_onto_stl(stl, points) out_points2 = project_line_onto_stl(stl, pa, pb, npoints=11)
#out_points3 = project_curve_onto_stl(stl, points, npoints=11)
[docs]def build(): from pyNastran.bdf.bdf import BDF model = BDF(debug=False) xyz1 = [0., 0., 0.] xyz2 = [0., 1., 0.] xyz3 = [1., 1., 0.] xyz4 = [1., 0., 0.] model.add_grid(1, xyz=xyz1) model.add_grid(2, xyz=xyz2) model.add_grid(3, xyz=xyz3) model.add_grid(4, xyz=xyz4) model.add_cquad4(eid=1, pid=1000, nids=[1, 2, 3, 4]) model.add_pshell(pid=100, mid1=1000, t=0.1) model.add_mat1(mid=1000, E=1e7, G=None, nu=0.3)
if __name__ == '__main__': # pragma: no cover build() main()