Coverage for src/meshpy/mesh_creation_functions/beam_honeycomb.py: 96%
85 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-28 04:21 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2025 MeshPy Authors
4#
5# Permission is hereby granted, free of charge, to any person obtaining a copy
6# of this software and associated documentation files (the "Software"), to deal
7# in the Software without restriction, including without limitation the rights
8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9# copies of the Software, and to permit persons to whom the Software is
10# furnished to do so, subject to the following conditions:
11#
12# The above copyright notice and this permission notice shall be included in
13# all copies or substantial portions of the Software.
14#
15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21# THE SOFTWARE.
22"""This file has functions to create a honeycomb beam mesh."""
24import numpy as _np
26from meshpy.core.geometry_set import GeometryName as _GeometryName
27from meshpy.core.geometry_set import GeometrySet as _GeometrySet
28from meshpy.core.mesh import Mesh as _Mesh
29from meshpy.core.rotation import Rotation as _Rotation
30from meshpy.mesh_creation_functions.beam_basic_geometry import (
31 create_beam_mesh_line as _create_beam_mesh_line,
32)
33from meshpy.utils.nodes import get_min_max_nodes as _get_min_max_nodes
36def create_beam_mesh_honeycomb_flat(
37 mesh,
38 beam_class,
39 material,
40 width,
41 n_width,
42 n_height,
43 *,
44 n_el=1,
45 closed_width=True,
46 closed_height=True,
47 create_couplings=True,
48 add_sets=False,
49):
50 """Add a flat honeycomb structure. The structure will be created in the x-y
51 plane.
53 Args
54 ----
55 mesh: Mesh
56 Mesh that the honeycomb will be added to.
57 beam_class: Beam
58 Class that will be used to create the beam elements.
59 material: Material
60 Material for the beam.
61 width: float
62 Width of one honeycomb.
63 n_width: int
64 Number of honeycombs in x-direction.
65 n_height: int
66 Number of honeycombs in y-direction.
67 n_el: int
68 Number of elements per beam line.
69 closed_width: bool
70 If the last honeycombs in x-direction will be closed.
71 closed_height: bool
72 If the last vertical lines in y-direction will be created.
73 create_couplings: bool
74 If the nodes will be connected in this function.
75 add_sets: bool
76 If this is true the sets are added to the mesh and then displayed
77 in eventual VTK output, even if they are not used for a boundary
78 condition or coupling.
80 Return
81 ----
82 return_set: GeometryName
83 Set with nodes on the north, south, east and west boundaries. Those
84 sets only contains end nodes of lines, not the middle ones. The set
85 'all' contains all nodes.
86 """
88 def add_line(pointa, pointb):
89 """Shortcut to add line."""
90 return _create_beam_mesh_line(
91 mesh_honeycomb, beam_class, material, pointa, pointb, n_el=n_el
92 )
94 # Geometrical shortcuts.
95 sin30 = _np.sin(_np.pi / 6)
96 cos30 = _np.sin(2 * _np.pi / 6)
97 a = width * 0.5 / cos30
98 nx = _np.array([1.0, 0.0, 0.0])
99 ny = _np.array([0.0, 1.0, 0.0])
100 zig_zag_x = nx * width * 0.5
101 zig_zag_y = ny * a * sin30 * 0.5
103 # Create the honeycomb structure
104 mesh_honeycomb = _Mesh()
105 origin = _np.array([0, a * 0.5 * sin30, 0])
106 for i_height in range(n_height + 1):
107 # Start point for this zig-zag line.
108 base_row = origin + (2 * zig_zag_y + a * ny) * i_height
110 # If the first node is up or down of base_row.
111 if i_height % 2 == 0:
112 direction = 1
113 else:
114 direction = -1
116 for i_width in range(n_width + 1):
117 base_zig_zag = base_row + direction * zig_zag_y + width * i_width * nx
119 # Do not add a zig-zag line on the last run (that one is only
120 # for the remaining vertical lines).
121 if i_width < n_width:
122 add_line(
123 base_zig_zag, base_zig_zag + zig_zag_x - 2 * direction * zig_zag_y
124 )
125 add_line(
126 base_zig_zag + zig_zag_x - 2 * direction * zig_zag_y,
127 base_zig_zag + nx * width,
128 )
130 # Check where the vertical lines start.
131 if i_height % 2 == 0:
132 base_vert = base_zig_zag
133 else:
134 base_vert = base_zig_zag + zig_zag_x - 2 * direction * zig_zag_y
136 # Only add vertical lines at the end if closed_width.
137 if (i_width < n_width) or (direction == 1 and closed_width):
138 # Check if the vertical lines at the top should be added.
139 if not (i_height == n_height) or (not closed_height):
140 add_line(base_vert, base_vert + ny * a)
142 # List of nodes from the honeycomb that are candidates for connections.
143 honeycomb_nodes = [node for node in mesh_honeycomb.nodes if node.is_end_node]
145 # Add connections for the nodes with same positions.
146 if create_couplings:
147 mesh_honeycomb.couple_nodes(nodes=honeycomb_nodes)
149 # Get min and max nodes of the honeycomb.
150 min_max_nodes = _get_min_max_nodes(honeycomb_nodes)
152 # Return the geometry set.
153 return_set = _GeometryName()
154 return_set["north"] = min_max_nodes["y_max"]
155 return_set["east"] = min_max_nodes["x_max"]
156 return_set["south"] = min_max_nodes["y_min"]
157 return_set["west"] = min_max_nodes["x_min"]
158 return_set["all"] = _GeometrySet(mesh_honeycomb.elements)
160 mesh.add(mesh_honeycomb)
161 if add_sets:
162 mesh.add(return_set)
163 return return_set
166def create_beam_mesh_honeycomb(
167 mesh,
168 beam_class,
169 material,
170 diameter,
171 n_circumference,
172 n_axis,
173 *,
174 n_el=1,
175 closed_top=True,
176 vertical=True,
177 add_sets=False,
178):
179 """Wrap a honeycomb structure around a cylinder. The cylinder axis will be
180 the z-axis.
182 Args
183 ----
184 mesh: Mesh
185 Mesh that the honeycomb will be added to.
186 beam_class: Beam
187 Class that will be used to create the beam elements.
188 material: Material
189 Material for the beam.
190 diameter: float
191 Diameter of the cylinder.
192 n_circumference: int
193 Number of honeycombs around the diameter. If vertical is False this
194 has to be an odd number.
195 n_axis: int
196 Number of honeycombs in axial-direction.
197 n_el: int
198 Number of elements per beam line.
199 closed_top: bool
200 If the last honeycombs in axial-direction will be closed.
201 vertical: bool
202 If there are vertical lines in the honeycomb or horizontal.
203 add_sets: bool
204 If this is true the sets are added to the mesh and then displayed
205 in eventual VTK output, even if they are not used for a boundary
206 condition or coupling.
208 Return
209 ----
210 return_set: GeometryName
211 Set with nodes on the top and bottom boundaries. Those sets only
212 contains end nodes of lines, not the middle ones. The set "all"
213 contains all nodes.
214 """
216 # Calculate the input values for the flat honeycomb mesh.
217 if vertical:
218 width = diameter * _np.pi / n_circumference
219 closed_width = False
220 closed_height = closed_top
221 rotation = _Rotation([0, 0, 1], _np.pi / 2) * _Rotation([1, 0, 0], _np.pi / 2)
222 n_height = n_axis
223 n_width = n_circumference
224 else:
225 if not n_circumference % 2 == 0:
226 raise ValueError(
227 "There has to be an even number of elements along the diameter in horizontal mode. "
228 "Given: {}!".format(n_circumference)
229 )
230 H = diameter * _np.pi / n_circumference
231 r = H / (1 + _np.sin(_np.pi / 6))
232 width = 2 * r * _np.cos(_np.pi / 6)
233 closed_width = closed_top
234 closed_height = False
235 rotation = _Rotation([0, 1, 0], -0.5 * _np.pi)
236 n_height = n_circumference - 1
237 n_width = n_axis
239 # Create the flat mesh, do not create couplings, as they will be added
240 # later in this function, where also the diameter nodes will be
241 # connected.
242 mesh_temp = _Mesh()
243 honeycomb_sets = create_beam_mesh_honeycomb_flat(
244 mesh_temp,
245 beam_class,
246 material,
247 width,
248 n_width,
249 n_height,
250 n_el=n_el,
251 closed_width=closed_width,
252 closed_height=closed_height,
253 create_couplings=False,
254 )
256 # Move the mesh to the correct position.
257 mesh_temp.rotate(rotation)
258 mesh_temp.translate([diameter / 2, 0, 0])
259 mesh_temp.wrap_around_cylinder()
261 # Add connections for the nodes with same positions.
262 honeycomb_nodes = [node for node in mesh_temp.nodes if node.is_end_node]
263 mesh_temp.couple_nodes(nodes=honeycomb_nodes)
265 # Return the geometry set'
266 return_set = _GeometryName()
267 return_set["all"] = honeycomb_sets["all"]
268 if vertical:
269 return_set["top"] = honeycomb_sets["north"]
270 return_set["bottom"] = honeycomb_sets["south"]
271 else:
272 return_set["top"] = honeycomb_sets["east"]
273 return_set["bottom"] = honeycomb_sets["west"]
274 if add_sets:
275 mesh_temp.add(return_set)
277 # Add to this mesh
278 mesh.add_mesh(mesh_temp)
280 return return_set