From 244f5d110da2e0687aa0b367e40dcd6d9a96453b Mon Sep 17 00:00:00 2001 From: Yuxiang Date: Wed, 20 Dec 2023 02:13:44 +0100 Subject: [PATCH] save updates --- .vscode/settings.json | 2 +- .../src/geometries/composite_geometry.py | 268 +++++++++++------- doublezigzag_150x150x11x150.stl | Bin 0 -> 6284 bytes honeycomb_150x150x10x150.stl | Bin 0 -> 6284 bytes honeycomb_700x150x10x150.stl | Bin 0 -> 13884 bytes .../try_new_thought/trail_3.py | 145 ++++++++-- .../try_new_thought/trail_double_zigzag.py | 167 +++++++++++ 7 files changed, 451 insertions(+), 131 deletions(-) create mode 100644 doublezigzag_150x150x11x150.stl create mode 100644 honeycomb_150x150x10x150.stl create mode 100644 honeycomb_700x150x10x150.stl create mode 100644 some_thoughts_20230822_new/try_new_thought/trail_double_zigzag.py diff --git a/.vscode/settings.json b/.vscode/settings.json index 56e1baf..b7648e4 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,6 +1,6 @@ { "[python]": { - "editor.defaultFormatter": "ms-python.python" + "editor.defaultFormatter": "ms-python.black-formatter" }, "python.formatting.provider": "none", "docify.commentService.docstringFormat": "Doxygen", diff --git a/amworkflow/src/geometries/composite_geometry.py b/amworkflow/src/geometries/composite_geometry.py index 9f9b9cc..2829442 100644 --- a/amworkflow/src/geometries/composite_geometry.py +++ b/amworkflow/src/geometries/composite_geometry.py @@ -44,10 +44,9 @@ from amworkflow.src.utils.meter import timer -def polygon_maker(side_num: int, - side_len: float, - rotate: float = None, - bound: bool = False): +def polygon_maker( + side_num: int, side_len: float, rotate: float = None, bound: bool = False +): """ @brief Creates a regular polygon. The polygon is oriented counterclockwise around the origin. If bound is True the polygon will be only the boundary (TopoDS_Wire) the polygon. @param side_num Number of sides of the polygon. @@ -58,8 +57,8 @@ def polygon_maker(side_num: int, """ vertices = [] sides = [] - r = side_len / (2 * np.sin(np.pi/side_num)) - ang = np.linspace(0,2*np.pi, side_num + 1) + r = side_len / (2 * np.sin(np.pi / side_num)) + ang = np.linspace(0, 2 * np.pi, side_num + 1) if rotate != None: ang += rotate for i in range(side_num + 1): @@ -69,26 +68,29 @@ def polygon_maker(side_num: int, t_pnt = gp_Pnt(x, y, 0) vertices.append(t_pnt) if i != 0: - if i == side_num : + if i == side_num: t_edge = create_edge(pnt1=vertices[0], pnt2=vertices[-1]) sides.append(t_edge) wire = create_wire(wire, t_edge) elif i == 1: - t_edge = create_edge(pnt1=vertices[i-1], pnt2=t_pnt) + t_edge = create_edge(pnt1=vertices[i - 1], pnt2=t_pnt) sides.append(t_edge) wire = create_wire(t_edge) else: - t_edge = create_edge(pnt1=vertices[i-1], pnt2=t_pnt) + t_edge = create_edge(pnt1=vertices[i - 1], pnt2=t_pnt) sides.append(t_edge) wire = create_wire(wire, t_edge) - + face = create_face(wire) if bound: return wire else: return reverse(face) -def hexagon_multiplier(side_num: int, side_len: float, iter_num: int, wall: float, center: gp_Pnt = None) -> TopoDS_Face: + +def hexagon_multiplier( + side_num: int, side_len: float, iter_num: int, wall: float, center: gp_Pnt = None +) -> TopoDS_Face: """ @brief Creates a hexagon with multiplier. This is an iterative approach to the topological sorting algorithm. @param side_num Number of sides in the hexagon. @@ -101,18 +103,20 @@ def hexagon_multiplier(side_num: int, side_len: float, iter_num: int, wall: floa def multiplier_unit(face_odd, face_even, unit_len): """ - @brief Fuse a unit of mass to an odd or even side number. - @param face_odd face with odd side number ( numpy array ) - @param face_even face with even side number ( numpy array ) - @param unit_len length of side ( numpy array ) - @return one unit of hexagon multiplication. + @brief Fuse a unit of mass to an odd or even side number. + @param face_odd face with odd side number ( numpy array ) + @param face_even face with even side number ( numpy array ) + @param unit_len length of side ( numpy array ) + @return one unit of hexagon multiplication. """ if side_num % 2 == 0: face = face_even - ang = np.linspace(0,2*np.pi, side_num + 1) + ang = np.linspace(0, 2 * np.pi, side_num + 1) else: face = face_odd - ang = np.linspace(np.pi / side_num,2*np.pi + np.pi / side_num, side_num + 1) + ang = np.linspace( + np.pi / side_num, 2 * np.pi + np.pi / side_num, side_num + 1 + ) cnt = get_face_center_of_mass(face) len = unit_len * 0.5 / np.tan(np.pi / side_num) * 2 face_collect = [face] @@ -125,16 +129,19 @@ def multiplier_unit(face_odd, face_even, unit_len): face_collect.append(rot_face) fuse = fuser(fuse, rot_face) return fuse + # This function will generate a hollow carver for each iteration. for i in range(iter_num): # This function computes the fuse of the carver. if i == 0: - face_even = hollow_carver(polygon_maker(side_num, side_len, rotate=-np.pi / side_num), wall) - face_odd = hollow_carver(polygon_maker(side_num, side_len),wall) + face_even = hollow_carver( + polygon_maker(side_num, side_len, rotate=-np.pi / side_num), wall + ) + face_odd = hollow_carver(polygon_maker(side_num, side_len), wall) fuse = multiplier_unit(face_odd, face_even, unit_len=side_len) else: face_even = fuse - # face_odd = rotate_face(fuse,angle = -np.pi / side_num) + # face_odd = rotate_face(fuse,angle = -np.pi / side_num) face_odd = fuse fuse = multiplier_unit(face_odd, face_even, unit_len=3 * i * side_len) # translate center to use list of coordinates @@ -142,22 +149,23 @@ def multiplier_unit(face_odd, face_even, unit_len): translate(fuse, list(center.Coord())) return fuse -def isoceles_triangle_maker(bbox_len:float, bbox_wid: float, thickness: float = None): + +def isoceles_triangle_maker(bbox_len: float, bbox_wid: float, thickness: float = None): """ - @brief (Having problem with wall thickness now.) Create isoceles triangulation. This is a function to create isoceles triangulation of a bounding box and its widest corner - @param bbox_len length of bounding box of the triangle - @param bbox_wid width of bounding box of the triangle - @param thickness thickness of the wall of the triangle - @return a hollowed triangle face. + @brief (Having problem with wall thickness now.) Create isoceles triangulation. This is a function to create isoceles triangulation of a bounding box and its widest corner + @param bbox_len length of bounding box of the triangle + @param bbox_wid width of bounding box of the triangle + @param thickness thickness of the wall of the triangle + @return a hollowed triangle face. """ h = bbox_wid l = bbox_len t = thickness if thickness != None else 0 - ary = np.array([h,l]) - r1 = np.linalg.norm(ary*0.5, 2) - ang1 = np.arccos((h*0.5) / r1) + ary = np.array([h, l]) + r1 = np.linalg.norm(ary * 0.5, 2) + ang1 = np.arccos((h * 0.5) / r1) ang2 = np.pi - ang1 - r = [h*0.5, r1, r1] + r = [h * 0.5, r1, r1] thet = np.arctan(l * 0.5 / h) thet1 = np.arctan(h / l) t1 = t / np.sin(thet) @@ -168,33 +176,35 @@ def isoceles_triangle_maker(bbox_len:float, bbox_wid: float, thickness: float = r2 = np.linalg.norm(0.5 * ary1, 2) ri = [0.5 * h1, r2, r2] ang = [0, ang2, -ang2] + def worker(r): """ - @brief Creates a triangle in a polarized coordinate system. - @param r The radius of each vertex. - @return a triangle + @brief Creates a triangle in a polarized coordinate system. + @param r The radius of each vertex. + @return a triangle """ pt_lst = [] # Add a point on the plane. for i in range(3): x = r[i] * np.sin(ang[i]) y = r[i] * np.cos(ang[i]) - pnt = gp_Pnt(x,y,0) + pnt = gp_Pnt(x, y, 0) pt_lst.append(pnt) # Create a wire for each point in the list of points. for i in range(3): # Create a wire for the i th point of the point. if i == 0: - edge = create_edge(pt_lst[i], pt_lst[i+1]) + edge = create_edge(pt_lst[i], pt_lst[i + 1]) wire = create_wire(edge) elif i != 2: - edge = create_edge(pt_lst[i], pt_lst[i+1]) + edge = create_edge(pt_lst[i], pt_lst[i + 1]) wire = create_wire(wire, edge) else: edge = create_edge(pt_lst[i], pt_lst[0]) wire = create_wire(wire, edge) face = create_face(wire) return reverse(face) + outer_face = worker(r) # The face to be cuttered. if t != 0: @@ -203,19 +213,23 @@ def worker(r): return new_face else: return outer_face - -def create_sym_hexagon1_infill(total_len: float, total_wid:float, height:float, th: float) : + + +def create_sym_hexagon1_infill( + total_len: float, total_wid: float, height: float, th: float +): """ - @brief Create an infill pattern using symmetrical hexagon with defined len, height and numbers. - @param total_len total length of the bounding box. - @param total_wid total wid of the bounding box. - @param height height of the prism. This is the same as height of the hexagon. - @param th thickness of the wall of the hexagon. - @return + @brief Create an infill pattern using symmetrical hexagon with defined len, height and numbers. + @param total_len total length of the bounding box. + @param total_wid total wid of the bounding box. + @param height height of the prism. This is the same as height of the hexagon. + @param th thickness of the wall of the hexagon. + @return """ - p0 = [0,th * 0.5] + p0 = [0, th * 0.5] p1 = [] + def find_intersect(lines: np.ndarray) -> np.ndarray: parallel = False coplanarity = False @@ -226,18 +240,18 @@ def find_intersect(lines: np.ndarray) -> np.ndarray: L2 = (pt4 - pt3) / np.linalg.norm(pt4 - pt3) V1 = pt4 - pt1 D1 = np.linalg.norm(V1) - if np.isclose(np.dot(L1, L2),0) or np.isclose(np.dot(L1, L2),-1): + if np.isclose(np.dot(L1, L2), 0) or np.isclose(np.dot(L1, L2), -1): parallel = True print("Two lines are parallel.") - return np.full((3,1), np.nan) + return np.full((3, 1), np.nan) indicate = np.linalg.det(np.array([V1, L1, L2])) if np.abs(indicate) < 1e-8: coplanarity = True else: print("lines are not in the same plane.") - return np.full((1,3), np.nan) + return np.full((1, 3), np.nan) if coplanarity and not parallel: - if np.isclose(D1,0): + if np.isclose(D1, 0): return pt1 else: pt5_pt4 = np.linalg.norm(np.cross(V1, L1)) @@ -251,18 +265,22 @@ def find_intersect(lines: np.ndarray) -> np.ndarray: o = L1 * pt1_o + pt1 return o + def index2array(ind: list, array: np.ndarray): real_array = [] - for i,v in enumerate(ind): + for i, v in enumerate(ind): item = [] for j, vv in enumerate(v): item.append(array[vv]) real_array.append(item) return real_array -def polygon_interpolater(plg: np.ndarray, step_len: float = None, num: int = None, isclose: bool = True): + +def polygon_interpolater( + plg: np.ndarray, step_len: float = None, num: int = None, isclose: bool = True +): def deter_dum(line: np.ndarray): - ratio = step_len / np.linalg.norm(line[0]-line[1]) + ratio = step_len / np.linalg.norm(line[0] - line[1]) if ratio > 0.75: num = 0 elif (ratio > 0.4) and (ratio <= 0.75): @@ -275,9 +293,10 @@ def deter_dum(line: np.ndarray): num = 4 elif (ratio > 0.14) and (ratio <= 0.19): num = 5 - elif ratio <=0.14: + elif ratio <= 0.14: num = 7 return num + new_plg = plg pos = 0 n = 1 @@ -286,18 +305,19 @@ def deter_dum(line: np.ndarray): for i, pt in enumerate(plg): if i == len(plg) - n: break - line = np.array([pt, plg[i+1]]) + line = np.array([pt, plg[i + 1]]) if num is not None: p_num = num else: p_num = deter_dum(line) insert_p = linear_interpolate(line, p_num) - new_plg = np.concatenate((new_plg[:pos+1],insert_p, new_plg[pos+1:])) - pos +=p_num+1 + new_plg = np.concatenate((new_plg[: pos + 1], insert_p, new_plg[pos + 1 :])) + pos += p_num + 1 return new_plg + # @timer -class CreateWallByPointsUpdate(): +class CreateWallByPointsUpdate: def __init__(self, coords: list, th: float, height: float, is_close: bool = True): self.coords = Pnt(coords).coords self.height = height @@ -314,6 +334,7 @@ def __init__(self, coords: list, th: float, height: float, is_close: bool = True self.ths = [] self.lft_coords = [] self.rgt_coords = [] + self.volume = 0 self.side_coords: list self.create_sides() self.pnts = Segments(self.side_coords) @@ -322,22 +343,24 @@ def __init__(self, coords: list, th: float, height: float, is_close: bool = True self.loop_generator = nx.simple_cycles(self.G) self.check_pnt_in_wall() self.postprocessing() - + def create_sides(self): if self.R is not None: self.coords = polygon_interpolater(self.coords, self.interpolate) self.coords = bender(self.coords, self.R) self.coords = [i for i in self.coords] self.th *= 0.5 - for i,p in enumerate(self.coords): + for i, p in enumerate(self.coords): if i != len(self.coords) - 1: - self.central_segments.append([self.coords[i], self.coords[i+1]]) - a1 = self.coords[i+1] - self.coords[i] + self.central_segments.append([self.coords[i], self.coords[i + 1]]) + a1 = self.coords[i + 1] - self.coords[i] if i == 0: if self.is_close: dr = angular_bisector(self.coords[-1] - p, a1) # ang = angle_of_two_arrays(dir_vecs[i-1],dr) - ang2 = angle_of_two_arrays(laterality_indicator(p - self.coords[-1], True), dr) + ang2 = angle_of_two_arrays( + laterality_indicator(p - self.coords[-1], True), dr + ) ang_th = ang2 if ang2 > np.pi / 2: dr *= -1 @@ -347,8 +370,10 @@ def create_sides(self): dr = laterality_indicator(a1, True) nth = self.th else: - dr = angular_bisector(-self.vecs[i-1], a1) - ang2 = angle_of_two_arrays(laterality_indicator(self.vecs[i-1], True), dr) + dr = angular_bisector(-self.vecs[i - 1], a1) + ang2 = angle_of_two_arrays( + laterality_indicator(self.vecs[i - 1], True), dr + ) ang_th = ang2 if ang2 > np.pi / 2: dr *= -1 @@ -358,12 +383,14 @@ def create_sides(self): if self.is_close: self.central_segments.append([self.coords[i], self.coords[0]]) a1 = self.coords[0] - self.coords[i] - dr = angular_bisector(-self.vecs[i-1], a1) - ang2 = angle_of_two_arrays(laterality_indicator(self.vecs[i-1], True), dr) + dr = angular_bisector(-self.vecs[i - 1], a1) + ang2 = angle_of_two_arrays( + laterality_indicator(self.vecs[i - 1], True), dr + ) ang_th = ang2 if ang2 > np.pi / 2: dr *= -1 - ang_th = np.pi - ang2 + ang_th = np.pi - ang2 nth = np.abs(self.th / np.cos(ang_th)) else: dr = laterality_indicator(a1, True) @@ -382,48 +409,52 @@ def create_sides(self): else: self.rgt_coords = self.rgt_coords[::-1] self.side_coords = self.lft_coords + self.rgt_coords + [self.lft_coords[0]] - + def check_pnt_in_wall(self): for pnt, coord in self.pnts.pts_index.items(): for vec in self.central_segments: - lmbda, dist = shortest_distance_point_line(vec,coord) + lmbda, dist = shortest_distance_point_line(vec, coord) if dist < 0.9 * self.th: print(f"pnt:{pnt},dist:{dist},lmbda:{lmbda}, vec:{vec}") - self.in_wall_pts_list.update({pnt:True}) + self.in_wall_pts_list.update({pnt: True}) break - - def visualize(self, display_polygon: bool = True,display_central_path:bool = False,all_polygons:bool = False): + def visualize( + self, + display_polygon: bool = True, + display_central_path: bool = False, + all_polygons: bool = False, + ): # Extract the x and y coordinates and IDs a = self.pnts.pts_index x = [coord[0] for coord in a.values()] y = [coord[1] for coord in a.values()] ids = list(a.keys()) # Get the point IDs # Create a scatter plot in 2D - plt.subplot(1,2,1) + plt.subplot(1, 2, 1) # plt.figure() plt.scatter(x, y) # Annotate points with IDs for i, (xi, yi) in enumerate(zip(x, y)): - plt.annotate(f'{ids[i]}', (xi, yi), fontsize=12, ha='right') - + plt.annotate(f"{ids[i]}", (xi, yi), fontsize=12, ha="right") + if display_polygon: if all_polygons: display_loops = nx.simple_cycles(self.G) else: display_loops = self.result_loops - for lp in display_loops: + for lp in display_loops: coords = [self.pnts.pts_index[i] for i in lp] x = [point[0] for point in coords] y = [point[1] for point in coords] - plt.plot(x + [x[0]], y + [y[0]], linestyle='-', marker='o') + plt.plot(x + [x[0]], y + [y[0]], linestyle="-", marker="o") if display_central_path: talist = np.array(self.coords).T x1 = talist[0] y1 = talist[1] - plt.plot(x1, y1, 'bo-', label='central path', color = "b") - + plt.plot(x1, y1, "bo-", label="central path", color="b") + # Create segments by connecting consecutive points # a_subtitute = np.array(self.side_coords) @@ -431,32 +462,46 @@ def visualize(self, display_polygon: bool = True,display_central_path:bool = Fal # x2 = toutput1[0] # y2 = toutput1[1] # plt.plot(x2, y2, 'ro-', label='outer line') - + # Set labels and title - plt.xlabel('X-axis') - plt.ylabel('Y-axis') - plt.title('Points With Polygons detected') + plt.xlabel("X-axis") + plt.ylabel("Y-axis") + plt.title("Points With Polygons detected") - plt.subplot(1,2,2) + plt.subplot(1, 2, 2) # layout = nx.spring_layout(self.G) layout = nx.circular_layout(self.G) # Draw the nodes and edges - nx.draw(self.G, pos=layout, with_labels=True, node_color='skyblue', font_size=10, node_size=300) + nx.draw( + self.G, + pos=layout, + with_labels=True, + node_color="skyblue", + font_size=10, + node_size=300, + ) plt.title("Multi-Digraph") - plt.tight_layout() + plt.tight_layout() # Show the plot plt.show() - + def get_loops(self): return [i for i in self.all_loops if len(i) > 2] def visualize_graph(self): layout = nx.spring_layout(self.G) # Draw the nodes and edges - nx.draw(self.G, pos=layout, with_labels=True, node_color='skyblue', font_size=10, node_size=500) + nx.draw( + self.G, + pos=layout, + with_labels=True, + node_color="skyblue", + font_size=10, + node_size=500, + ) plt.title("NetworkX Graph Visualization") plt.show() - + def Shape(self): loop_r = self.rank_result_loops() print(loop_r) @@ -471,14 +516,14 @@ def Shape(self): poly_r = cutter3D(poly_r, poly_c) self.poly = poly_r if not np.isclose(self.height, 0): - wall_compound = create_prism(self.poly, [0,0,self.height]) + wall_compound = create_prism(self.poly, [0, 0, self.height]) faces = topo_explorer(wall_compound, "face") wall_shell = sewer(faces) self.wall = solid_maker(wall_shell) return self.wall else: return self.poly - + def postprocessing(self): correct_loop_count = 0 for lp in self.loop_generator: @@ -487,28 +532,32 @@ def postprocessing(self): in_wall_pt_count = 0 virtual_vector_count = 0 if len(lp) > 2: - for i,pt in enumerate(lp): + for i, pt in enumerate(lp): if i == 0: if pt in self.in_wall_pts_list: in_wall_pt_count += 1 - if (lp[-1],pt) in self.pnts.virtual_vector: + if (lp[-1], pt) in self.pnts.virtual_vector: virtual_vector_count += 1 else: if pt in self.in_wall_pts_list: in_wall_pt_count += 1 - if (lp[i-1],pt) in self.pnts.virtual_vector: + if (lp[i - 1], pt) in self.pnts.virtual_vector: virtual_vector_count += 1 - if (in_wall_pt_count > 0) or (virtual_vector_count > 1) or ((in_wall_pt_count == 0) and (virtual_vector_count > 0)): - # if (in_wall_pt_count > 0): + if ( + (in_wall_pt_count > 0) + or (virtual_vector_count > 1) + or ((in_wall_pt_count == 0) and (virtual_vector_count > 0)) + ): + # if (in_wall_pt_count > 0): visible_loop = False - break + break else: - real_loop = False + real_loop = False if real_loop and visible_loop: self.result_loops.append(lp) for pt in lp: if pt not in self.in_loop_pts_list: - self.in_loop_pts_list.update({pt:[correct_loop_count]}) + self.in_loop_pts_list.update({pt: [correct_loop_count]}) else: self.in_loop_pts_list[pt].append(correct_loop_count) correct_loop_count += 1 @@ -532,18 +581,25 @@ def postprocessing(self): filtered_lp = np.where(loop_counter > 1)[0].tolist() print("filtered:", filtered_lp) print("vote:", loop_counter) - self.result_loops = [v for i,v in enumerate(self.result_loops) if i not in filtered_lp] + self.result_loops = [ + v for i, v in enumerate(self.result_loops) if i not in filtered_lp + ] print("result:", self.result_loops) def rank_result_loops(self): areas = np.zeros(len(self.result_loops)) - for i,lp in enumerate(self.result_loops): + for i, lp in enumerate(self.result_loops): lp_coord = [self.pnts.pts_index[i] for i in lp] area = p_get_face_area(lp_coord) areas[i] = area rank = np.argsort(areas).tolist() - self.result_loops = sorted(self.result_loops, key=lambda x: rank.index(self.result_loops.index(x)),reverse=True) + self.volume = (2 * np.max(areas) - np.sum(areas)) * self.height * 1e-6 + self.result_loops = sorted( + self.result_loops, + key=lambda x: rank.index(self.result_loops.index(x)), + reverse=True, + ) return self.result_loops - - def occ_pnt(self,coord) -> gp_Pnt: - return gp_Pnt(*coord) \ No newline at end of file + + def occ_pnt(self, coord) -> gp_Pnt: + return gp_Pnt(*coord) diff --git a/doublezigzag_150x150x11x150.stl b/doublezigzag_150x150x11x150.stl new file mode 100644 index 0000000000000000000000000000000000000000..0a8377302e64b322c6316742d2c3b94bdf928f89 GIT binary patch literal 6284 zcma)AJ&zkz6dekBDt8Ds3tiX}DtmX51|d;oJA_c!wSo|eS6PIJl;Vj7DFrkrP*NZf zL6Mar(c0@O6@MVVfQE*ag8O~Xn|YSUQr6cq_uPB#=bOphvk&KYzPVgK{AzVE|LoEH zqs!HUd4F=ZKe>HpezrRQ{K5Lm_5DZlk1tkV-&)hn^Ltn4_byhq)(`K0(v3fl+C_9A%!fpX+(kyVfgpiqF?y7LN;Lsp zh;jYl>Gh85x{7;r2_TC2=&qOaEjg(WUy3PKOv3s8{4V!_Sy$@ zy;XCHW9fGgeY6*0c^6hSL#b>dSBMyAx|Ln?F40;XViy^5jQ%Mb%cy9fKSaAWe@g|I z6uJ07fBSt{t^`qQWe!I^wr~QlzZ=!aUO0sC{gM542%#}yWm|~Y=S#|lRZZfM`1T_( z)9S;~iHZ4W*yT9=PGS+vOQXqry*v2T4Kb{SJ%%yS5`IY5}Uh$JPk4H z(-I@xK1*Y}xUX2AhF~q(-5;h zEivNnZ{MWnDa|nM*OWPOfA$$7@=lxpQvbXTtgIgSgq%#eF69A&>d`_3tDSa9HG}$m zE==&LIA!yZQ`0SgKg;^jv(|PW#8Evb()9tN zv=Om)4PkvH^+>+4YsFt$VjqVQaYj{fa5NN)Zh#x%t6OcSvcCyu9EzyU?RY;UEktXo zVZ2P8RiBb#*y!UQ3I{B=Bg1nP0p5{^hb5Zn1oePZrt(s{E_GpLID(qMJP-#4vf!~ z<*O`kBjb|==CDTC29VxG1X`*QRs&&Wc11?E`w`aRf~d02kfJz>Jl+lJvk(1=ZXhg+ zBUZ?^M&RzlXwkRw(Z60npM!h~L-Yu-)ptwoxT^k3ee)i5qV$#j|7UhdjV3RJ{?SKR zZEGL?RE)b_l;w{j_$}30VNCU=W=IR6TA20p&_HYl40EE2;%XmN)#44gohyw4nN( HTi5*$s{5HE literal 0 HcmV?d00001 diff --git a/honeycomb_150x150x10x150.stl b/honeycomb_150x150x10x150.stl new file mode 100644 index 0000000000000000000000000000000000000000..c7d3a4f22de96da702cbbfdacd90685afa5e74bb GIT binary patch literal 6284 zcmb7|zmHT^6vtm|XtXdfTTfvyu~;!Cbh5KE(ves|L(~Fy2P_CfkQhlcx#2IcA`zM? zj1(kV7)#yRrLd-;v!dpY;QYMjd+xrs?=i*ho%hb?{J!Vj-K*EG%r5_UcYW{v>Tve$ zgV~pNS9fO1{j1CUFD}ontq#7sv%bB)^UnvnSHuM0$cEneOH{qMT&lD9?__6mt4 zRo)}oV?v_aKBqU`Q+|dXCVK|(uNnvXE^DMKP+B5Mf$d6&jpzpYnxe%_2g7^O`-c@k zw297y=;G+<)7ImmJ=eEh^Ub_w5pAGt1B@si)a+ty9{>7^a^}6@2q;;^b|c2DO@VgeIvLpwtVtY*>#NG7b?PtgF>IO;r~Q-0J#sM8q) zc&L#(9gJreZPZ+CiLx_sn1D`99Xd#$RPxY>qBRG-b)K`t&=if+tChk%&+qh8C!BTz z%#Y!jh!Y2Q5$<`{Z<7I{#eI$)QwdtMyH+Ng)V%Gha`fW!E1xKfIpMYM4`18KIYO|M zYvT3AcmLTXkL-kYU1Sk*73+p5s>?j7gxlcfyuSbBmk+hx`CJN5I70MloY+~$>BBb8 zzy5je=BMxc%bLR9Z0~Ie(wky>_VNFJ)Li8=V;=R`*qXf>CrGMlVxsEjrAz>?Gwa1E z%Ij$aOQBbkt0T`9c6fL~rju9W1l6OSC2a59lRtm+nyYp7J*22MZ%KF@9=Y|d-d zZiv!PR2w-C$rbvJbgUtcGYFDOkyEm*Z#kmP+;6Wo)fO5jKj*uT#xX{)lxyO(O=_Jb zIxST%-d^f@8sRqFE1RMkdDT_tpzgtXpJ#fUAgQW}37fl;FPpETH=^ZiVdUdPo$a>I zHQ)S`vYbvxxzC+w4g4;_Ud4y&f+synA&o>nw_Dn$4-rMFvX{CVJj|VlHWHy3)i{6- zd5Ii)6!m6rgDHD4sq1t;j$A<;24w>G#IOyj?Uw#iWHo?{4#$B6N?{fv4B3l*Z2tzT!P0)hT&pHCRYSp*EOx0HSx5NGCGJ7Yx+{b|-q3 zPXi?KUgIq2?IqC+6#gN71-lca(6-9zPs#3N!`?@WSd?--VO5y!6(<3blMoF|BZ7AAx!Mq*YWgx51~W)lXom^27Fi)5xv5J?X+A!J~{ zLqu?)vkW0Wpcfg0G$JCdOnBYBFzCwQDsKJt2LuEmi`1`E^{sQy?Ry6@edgWssc%)C zAGgkVovY7X*}44M&4XL7?BCdV;m*#pH}_xO*<84~xp3+7&U5?MUwZlArw1?I+4ycQ55QFeYv{L>amS_43bwuzX2`>bv%8+RF&;RV^t| zGzJc2Y?plZ&pnPgg6z7b@M zOJFb5&P2=~`W*8!8F?bG7v{}H*sDRFyiT*hHMBiR?EU%yAHkGS^e|~4ff-G7C5q(g z&tLQrL;@{UfPea!i|>|Q8INvug)sGr$SosPq_74Os9_$erAmm#y+3?PD^%CVlMquk z4$;_q_?}Z&=NnZf6CUL=K}ijUb_1rftH-0M%h{~iQHlm4WykYErGp`K--Bajgej^j^29HbA`T7 zM}F%oZ+otgKucM}ddZB#EOg^2GOu(5M+91GG@kgmj~aT|^-4AF|AGizwMg6pOI;1g z=pN9^2-A?5kj@(Bul9m0bR(BIkU$OdO2;2WXiwBJ7Y($FGypVrhvCPZU9pv;mCsMG zH@CY)hYqf6|Mli&A3@Z!c1Poz|2;>9u3Dsz{`w92ys2c0%&R}{(i2^^NT8+gi3#f^ zzN1%2J0TiKV6Q#Uc7hu89(??L8VijBJ(-RJ3G8+BlS8u9L@pX;SH>JY?0O|}nB5V9 zmU0|X*A7GS3ig;Y2=EWI z>|(F!IE)ZmNT8)qV6W)OG5uFhrn_}CmFUg(F;RO>&5|*pyPaxid_AxHoJn5Q`ty)I z?C0vZo>FmaEDp)~oJO^a840H5%DmEgN90)<(8PZ`rtF`yH7nt>o42IlW??`5vxto z8@_t!{KGm6>FhMVH8{yRWQ}=ZwFyg1od4vL^EFp-9$0N~)S}%GC5T~&Uexw8bsIU> z>d{8MmO#5fpq>%iIv#jl3Qsy-%@C_iSYnN9AK&F%$(W-jcw&k6GGes}OGo3M>CF~PgG5<;rl=Ptd&bdY3h5?)tn6 ze`%I@FKjO(Xq|hlZN$3nS9oWIH$~B2M$j(mw{1*Z`}ni>>V6frh5_%aqP>iu_lMuf zfmiabwI?(iXs>lou-y=)Pq>E7D;y8J|BUu(ljA@lo{%SSFCL#;y~RT9Hwu&l>Xt^-5mZ`-3Ad7)P|dk}(Y=tj3`yv{JR|lPtvN>J>l7I#GR( zUS4;398OrB*b+4k{+u=FUF&PmDjD&_>XN3tjPPp7#=~a2&JwM>#IDs0l7(oq#ymk9 zZpjGVuSCNt7O|UR@;kYPD1E|fAIlZJYdKe<5&rTJ#l3j6S!15?8qPIjPZW)?H4Nxs zv{hq{pkD4FM84&@g7=2ertJD1X-mYDOvmTo)kp*HETg@Qpz+ZA>zG${%uA@#E(*RI zajLT&*I+zB!ea`8cQ}b7{H0Y3_xRDaT+I+P=K7t5#X;{u-`nj8siP;@wp>-B^auz^!AnnVu)AM)e9SnXpCsxt_Kprf4%gPteSF51DxI-KX!@e1$EJq(+PFWrUvx zZCgZ^Ydoc(Ub2)l z;>40qLLMIbRzC4=?Qp(nO5sn0}1SfQg6BgRVEQopc z3JL6GmyEH zUc<@ym?rFr%!W=JvI??BocsMGBs69u6LFiJy+Us@f#*!5rwOd1kU&eKQB;ZM6;xi> z%QQOc18Yvq8)hWOkqFQmX&B~Xt{H_hIqGGNT;0XZKW0IQq8#6{L_LVTFe6#RYED}@ zp@EjN2A)jtyPIq&RLUApCP0qe#@!vwCxeMPxp@-3* zH4@?98)P~YW=qDr<5g!3;#mf1j4x|Y#d=J?WpT>yKHOI@rZWhL+$vJL^Nes{e-z_7V>-dKIn`u7QMBuQ*QDS$@6~6~32{2B7uIZ*jcVn#mQ;SA5%@p#d`=>4w0Q5Y$1I Qk<-oy0)f>M?smiQIT2MeegFUf literal 0 HcmV?d00001 diff --git a/some_thoughts_20230822_new/try_new_thought/trail_3.py b/some_thoughts_20230822_new/try_new_thought/trail_3.py index 5cd1093..a38890d 100644 --- a/some_thoughts_20230822_new/try_new_thought/trail_3.py +++ b/some_thoughts_20230822_new/try_new_thought/trail_3.py @@ -1,7 +1,9 @@ import numpy as np +from scipy.optimize import fsolve from amworkflow.api import amWorkflow as aw from amworkflow.src.geometries.composite_geometry import CreateWallByPointsUpdate +from amworkflow.src.utils.writer import stl_writer # th = 8 # l = 10 @@ -45,10 +47,15 @@ # pmfo = np.vstack((pmf, pout_nd)) - def honeycomb_infill( - overall_length: float, line_width: float, honeycomb_num: int = 1 -) -> np.ndarray: + overall_length: float = None, + overall_width: float = None, + line_width: float = None, + honeycomb_num: int = 1, + angle: float = 1.1468, + regular: bool = False, + side_len: float = None, +): """Create honeycomb geometry. Args: @@ -61,7 +68,23 @@ def honeycomb_infill( """ - def half_honeycomb(origin: np.ndarray, side_length: float) -> np.ndarray: + def cal_len_wid(variable): + x, y = variable + x_rad = np.radians( + x + ) # Convert degrees to radians for trigonometric calculations + # eq1 = 2 * (np.cos(x_rad) + 1) * y + 1.5 * line_width - overall_length + eq1 = ( + (2 * np.cos(x_rad) + 2) * y * honeycomb_num + + 1.5 * line_width + - overall_length + ) + eq2 = 3 * line_width + 2 * y * np.sin(x_rad) - overall_width + return [eq1, eq2] + + def half_honeycomb( + origin: np.ndarray, side_length1: float, side_length2, angle: float + ) -> np.ndarray: """Create half of a honeycomb geometry. Args: @@ -76,26 +99,85 @@ def half_honeycomb(origin: np.ndarray, side_length: float) -> np.ndarray: """ points = np.zeros((5, 2)) points[0] = origin - points[1] = origin + np.array([0.5 * side_length, np.sqrt(3) * side_length / 2]) - points[2] = points[1] + np.array([side_length, 0]) - points[3] = points[2] + np.array([0.5 * side_length, -np.sqrt(3) * side_length / 2]) - points[4] = points[3] + np.array([0.5 * side_length, 0]) + # points[1] = origin + np.array([0.5 * side_length, np.sqrt(3) * side_length / 2]) + # points[2] = origin + np.array([length, 0]) + # points[3] = origin + np.array([0.5 * length, -np.sqrt(3) * side_length / 2]) + # points[4] = origin + np.array([0.5 * length, 0]) + points[1] = origin + np.array( + [side_length1 * np.cos(angle), side_length1 * np.sin(angle)] + ) + points[2] = points[1] + np.array([side_length2, 0]) + points[3] = points[2] + np.array( + [side_length1 * np.cos(angle), -side_length1 * np.sin(angle)] + ) + points[4] = points[3] + np.array([side_length2 * 0.5, 0]) return points - length = (overall_length - (3 * line_width)) / (3 * honeycomb_num) - start_point = np.array([0, line_width * 0.5]) - offset = np.array([1.5 * line_width, 0]) - overall_width = 3 * line_width + np.sqrt(3) * length - point_num = 16 + (honeycomb_num - 1) * 10 - half_points = np.zeros((int(point_num / 2) - 2, 2)) + if not regular: + if overall_width is not None: + overall_length -= line_width + overall_width -= line_width + initial_guesses = [ + (x, y) for x in range(0, 89, 10) for y in range(1, overall_length, 10) + ] + # Set for storing unique solutions + updated_solutions = set() + # Loop through each initial guess to find all possible solutions + for guess in initial_guesses: + sol = fsolve(cal_len_wid, guess) + # Round the solutions to 2 decimal places to avoid minor variations + rounded_sol = (round(sol[0], 10), round(sol[1], 10)) + + # Check if the solution is within the specified ranges and not already found + if 0 <= rounded_sol[0] <= 90 and 0 <= rounded_sol[1] <= 150: + updated_solutions.add(rounded_sol) + if updated_solutions: + ideal_sol = min(updated_solutions, key=lambda x: np.abs(x[0] - 45)) + else: + raise ValueError("No solution found") + length = ideal_sol[1] + angle = np.radians(ideal_sol[0]) + else: + length = ( + (overall_length - 1.5 * line_width) + / (2 * np.cos(angle) + 2) + / honeycomb_num + ) + overall_width = 3 * line_width + 2 * np.sin(angle) * length + start_point = np.array([0, line_width * 0.5]) + offset = np.array([line_width * 0.5 + length * 0.5, 0]) + else: + length = side_len + start_point = np.array([0, line_width * 0.5]) + offset = np.array([line_width * 0.5 + length * 0.5, 0]) + regular_length = ( + 1.5 * line_width + (2 * np.cos(angle) + 2) * side_len * honeycomb_num + ) + overall_length = regular_length + overall_width = 3 * line_width + 2 * np.sin(angle) * length + + # def cal_angle(x): + # return np.cos(x) + 1 - np.sin(x) - 0.75 * line_width / length + + # angle = fsolve(cal_angle, 0)[0] + + # unit_width = 4 * line_width + np.sqrt(3) * length + point_num = 12 + (honeycomb_num - 1) * 10 + half_points = np.zeros((int(point_num / 2), 2)) half_points[0] = start_point for i in range(honeycomb_num): if i == 0: start = start_point + offset - honeycomb_unit = half_honeycomb(start, length) + honeycomb_unit = half_honeycomb(start, length, length, angle) else: - start = start_point + offset + np.array([3 * length * i, 0]) - honeycomb_unit = half_honeycomb(start, length) + start = ( + start_point + + offset + + np.array([(2 * np.cos(angle) + 2) * length * i, 0]) + ) + honeycomb_unit = half_honeycomb(start, length, length, angle) + # if i == honeycomb_num - 1: + # honeycomb_unit[-1][0] += line_width * 0.5 half_points[i * 5 + 1 : i * 5 + 6] = honeycomb_unit another_half = half_points.copy() another_half[:, 1] = -another_half[:, 1] @@ -113,11 +195,26 @@ def half_honeycomb(origin: np.ndarray, side_length: float) -> np.ndarray: return points -ppt = honeycomb_infill(150, 8, 1) - - -wall = CreateWallByPointsUpdate(ppt, 8, 2) +# ppt = honeycomb_infill(150, 8, 3) +# ppt = honeycomb_infill( +# regular=True, side_len=63.26, angle=np.deg2rad(84.79), honeycomb_num=2, line_width=8 +# ) +ppt = honeycomb_infill( + overall_length=700, overall_width=150, line_width=10, honeycomb_num=3 +) +# 60 degree: 1.04719 +# 150x150x150x10: volume (L): 1.520399999948474 +# 700x150x10x150: volume (L): 5.014000002962625 +wall = CreateWallByPointsUpdate(ppt, 10, 150) +print(ppt) wall.visualize(all_polygons=False, display_central_path=True) +wall_shape = wall.Shape() +# stl_writer( +# wall_shape, +# "honeycomb_700x150x10x150", +# store_dir="/Users/yuxianghe/Documents/BAM/amworkflow_restructure", +# ) +print("volume (L):", wall.volume) # lft_coords = wall.lft_coords # rgt_coords = wall.rgt_coords # pieces = [] @@ -138,9 +235,9 @@ def half_honeycomb(origin: np.ndarray, side_length: float) -> np.ndarray: # def find_contours(image): # contours, _ = cv.findContours(image, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) - + # return contours - + # poly = wall.Shape() # wall.visualize(all_polygons=False, display_central_path=False) # aw.tool.write_stl(poly, "sucess_new_scheme",store_dir="/home/yhe/Documents/new_am2/amworkflow/some_thoughts_20230822_new/try_new_thought") @@ -154,4 +251,4 @@ def half_honeycomb(origin: np.ndarray, side_length: float) -> np.ndarray: # cv.drawContours(contour_image, contours, -1, (0, 255, 0), 2) # cv.imshow("Contours", contour_image) # cv.waitKey(0) -# cv.destroyAllWindows() \ No newline at end of file +# cv.destroyAllWindows() diff --git a/some_thoughts_20230822_new/try_new_thought/trail_double_zigzag.py b/some_thoughts_20230822_new/try_new_thought/trail_double_zigzag.py new file mode 100644 index 0000000..2ebbe7f --- /dev/null +++ b/some_thoughts_20230822_new/try_new_thought/trail_double_zigzag.py @@ -0,0 +1,167 @@ +import numpy as np +from scipy.optimize import fsolve + +from amworkflow.src.geometries.composite_geometry import CreateWallByPointsUpdate +from amworkflow.src.utils.writer import stl_writer + + +def zigzag_infill( + overall_length: float = None, + overall_width: float = None, + line_width: float = None, + zigzag_num: int = 1, + angle: float = 1.1468, + regular: bool = False, + side_len: float = None, +): + def cal_len_wid(variable): + x, y = variable + x_rad = np.radians( + x + ) # Convert degrees to radians for trigonometric calculations + # eq1 = 2 * (np.cos(x_rad) + 1) * y + 1.5 * line_width - overall_length + eq1 = ( + ((2 * np.cos(x_rad)) * y + line_width * np.sin(x_rad) * 2) * zigzag_num + + 2 * line_width + - overall_length + ) + eq2 = 3 * line_width + 2 * y * np.sin(x_rad) - overall_width + return [eq1, eq2] + + def half_zigzag( + origin: np.ndarray, side_length1: float, side_length2, angle: float + ) -> np.ndarray: + """Create half of a honeycomb geometry. + + Args: + origin: Origin of the half honeycomb. + length: Length of the honeycomb infill. + width: Width of the honeycomb infill. + line_width: Width of the honeycomb lines. + + Returns: + points: list of points defining half of the honeycomb geometry. + + """ + points = np.zeros((5, 2)) + points[0] = origin + # points[1] = origin + np.array([0.5 * side_length, np.sqrt(3) * side_length / 2]) + # points[2] = origin + np.array([length, 0]) + # points[3] = origin + np.array([0.5 * length, -np.sqrt(3) * side_length / 2]) + # points[4] = origin + np.array([0.5 * length, 0]) + points[1] = origin + np.array( + [side_length1 * np.cos(angle), side_length1 * np.sin(angle)] + ) + points[2] = points[1] + np.array([side_length2, 0]) + points[3] = points[2] + np.array( + [side_length1 * np.cos(angle), -side_length1 * np.sin(angle)] + ) + points[4] = points[3] + np.array([side_length2 * 0.5, 0]) + return points + + if not regular: + if overall_width is not None: + overall_length -= line_width + overall_width -= line_width + initial_guesses = [ + (x, y) for x in range(0, 89, 10) for y in range(1, overall_length, 10) + ] + # Set for storing unique solutions + updated_solutions = set() + # Loop through each initial guess to find all possible solutions + for guess in initial_guesses: + sol = fsolve(cal_len_wid, guess) + # Round the solutions to 2 decimal places to avoid minor variations + rounded_sol = (round(sol[0], 16), round(sol[1], 16)) + + # Check if the solution is within the specified ranges and not already found + if 0 <= rounded_sol[0] <= 90 and 0 <= rounded_sol[1] <= 150: + updated_solutions.add(rounded_sol) + if updated_solutions: + ideal_sol = min(updated_solutions, key=lambda x: np.abs(x[0] - 45)) + else: + raise ValueError("No solution found") + length = ideal_sol[1] + angle = np.radians(ideal_sol[0]) + else: + length = ( + (overall_length - 1.5 * line_width) / zigzag_num + ) - 2 * line_width * np.sin(angle) / (2 * np.cos(angle)) + overall_width = 3 * line_width + 2 * np.sin(angle) * length + start_point = np.array([0, line_width * 0.5]) + offset = np.array([line_width * 1 + line_width * np.sin(angle) * 0.5, 0]) + else: + length = side_len + start_point = np.array([0, line_width * 0.5]) + offset = np.array([line_width * 0.5 + length * 0.5, 0]) + regular_length = ( + 1.5 * line_width + (2 * np.cos(angle) + 2) * side_len * zigzag_num + ) + overall_length = regular_length + overall_width = 3 * line_width + 2 * np.sin(angle) * length + + # def cal_angle(x): + # return np.cos(x) + 1 - np.sin(x) - 0.75 * line_width / length + + # angle = fsolve(cal_angle, 0)[0] + + # unit_width = 4 * line_width + np.sqrt(3) * length + point_num = 12 + (zigzag_num - 1) * 10 + half_points = np.zeros((int(point_num / 2), 2)) + half_points[0] = start_point + for i in range(zigzag_num): + if i == 0: + start = start_point + offset + honeycomb_unit = half_zigzag( + start, length, line_width * np.sin(angle), angle + ) + else: + start = ( + start_point + + offset + + np.array( + [ + (2 * np.cos(angle) * length + 2 * line_width * np.sin(angle)) + * i, + 0, + ] + ) + ) + honeycomb_unit = half_zigzag( + start, length, line_width * np.sin(angle), angle + ) + # if i == zigzag_num - 1: + # honeycomb_unit[-1][0] += line_width * 0.5 + half_points[i * 5 + 1 : i * 5 + 6] = honeycomb_unit + another_half = half_points.copy() + another_half[:, 1] = -another_half[:, 1] + another_half = np.flipud(another_half) + points = np.concatenate((half_points, another_half), axis=0) + outer_points = np.array( + [ + [0, -overall_width * 0.5], + [overall_length, -overall_width * 0.5], + [overall_length, overall_width * 0.5], + [0, overall_width * 0.5], + ] + ) + points = np.concatenate((points, outer_points), axis=0) + return points + + +th = 11 + +ppt = zigzag_infill(overall_length=150, overall_width=150, line_width=th, zigzag_num=1) +# 60 degree: 1.04719 +# 150x150x150x10: volume(L): 1.3706792917581947 +# 150x150x150x11: volume(L): 1.4895852945607906 +wall = CreateWallByPointsUpdate(ppt, th, 150) +print(ppt) +wall.visualize(all_polygons=False, display_central_path=True) +wall_shape = wall.Shape() +stl_writer( + wall_shape, + "doublezigzag_150x150x11x150", + store_dir="/Users/yuxianghe/Documents/BAM/amworkflow_restructure", +) +print(wall.volume)