Coverage for src/beamme/mesh_creation_functions/applications/beam_stent.py: 99%
106 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-28 15:20 +0000
1# The MIT License (MIT)
2#
3# Copyright (c) 2018-2026 BeamMe 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 stent according to Auricchio 2012."""
24import numpy as _np
26from beamme.core.geometry_set import GeometryName as _GeometryName
27from beamme.core.geometry_set import GeometrySet as _GeometrySet
28from beamme.core.mesh import Mesh as _Mesh
29from beamme.core.rotation import Rotation as _Rotation
30from beamme.mesh_creation_functions.beam_arc import (
31 create_beam_mesh_arc_segment_via_rotation as _create_beam_mesh_arc_segment_via_rotation,
32)
33from beamme.mesh_creation_functions.beam_line import (
34 create_beam_mesh_line as _create_beam_mesh_line,
35)
36from beamme.utils.nodes import get_min_max_nodes as _get_min_max_nodes
39def create_stent_cell(
40 beam_class,
41 material,
42 width,
43 height,
44 fac_bottom=0.6,
45 fac_neck=0.55,
46 fac_radius=0.36,
47 alpha=0.46 * _np.pi,
48 S1=True,
49 S2=True,
50 S3=True,
51 n_el=1,
52):
53 """Create a cell of the stent.
55 This cell is on the x-y plane.
57 Args
58 ----
59 beam_class: Beam
60 Class that will be used to create the beam elements.
61 material: Material
62 Material for the beam.
63 width: float
64 Width of the total cell.
65 height: float
66 Height of the total cell.
67 fac_bottom: the ratio of the bottom's width to the cell's width
68 fac_neck: the ratio of the neck's width to the cell's width
69 fac_radius: the ratio of the S1's radius to the cell's width
70 alpha: radiant
71 The angle between the lines and horizontal line
72 n_el: int
73 Number of elements per beam line.
74 S1, S2, S3: bool
75 This check weather the curve S1, S2 or S3 will be created.
76 If the cell is on bottom of the stent flat S1 and S2 won't
77 be created. If the cell is on top of the flat S1 and S3
78 won't be created
80 ( these variables are described in a file )
82 Return
83 ----
84 mesh: Mesh
85 A mesh with this structure
86 """
87 mesh = _Mesh()
89 def add_line(pointa, pointb, n_el_line):
90 """Shortcut to add line."""
91 return _create_beam_mesh_line(
92 mesh, beam_class, material, pointa, pointb, n_el=n_el_line
93 )
95 def add_segment(center, axis_rotation, radius, angle, n_el_segment):
96 """Shortcut to add arc segment."""
97 return _create_beam_mesh_arc_segment_via_rotation(
98 mesh,
99 beam_class,
100 material,
101 center,
102 axis_rotation,
103 radius,
104 angle,
105 n_el=n_el_segment,
106 )
108 neck_width = width * fac_neck
109 bottom_width = width * fac_bottom
110 top_width = 2 * neck_width - bottom_width
111 radius = width * fac_radius
113 if S1:
114 neck_point = _np.array([-neck_width, height * 0.5, 0])
115 d = (height * 0.5 / _np.tan(alpha) + bottom_width - neck_width) / _np.sin(alpha)
116 CM = _np.array([-_np.sin(alpha), -_np.cos(alpha), 0]) * (d - radius)
117 MO = _np.array([_np.cos(alpha), -_np.sin(alpha), 0]) * _np.sqrt(
118 radius**2 - (d - radius) ** 2
119 )
120 S1_angle = _np.pi / 2 + _np.arcsin((d - radius) / radius)
121 S1_center1 = CM + MO + neck_point
122 S1_axis_rotation1 = _Rotation([0, 0, 1], 2 * _np.pi - S1_angle - alpha)
123 add_segment(S1_center1, S1_axis_rotation1, radius, S1_angle, n_el)
124 add_line([-bottom_width, 0, 0], mesh.nodes[-1].coordinates, 2 * n_el)
126 S1_center2 = 2 * neck_point - S1_center1
127 S1_axis_rotation2 = _Rotation([0, 0, 1], _np.pi - alpha - S1_angle)
129 add_segment(S1_center2, S1_axis_rotation2, radius, S1_angle, n_el)
130 add_line(
131 mesh.nodes[-1].coordinates, 2 * neck_point - [-bottom_width, 0, 0], 2 * n_el
132 )
134 if S3:
135 S3_radius = (width - bottom_width + height * 0.5 / _np.tan(alpha)) * _np.tan(
136 alpha / 2
137 )
138 S3_center = [-width, height * 0.5, 0] + S3_radius * _np.array([0, 1, 0])
139 S3_axis_rotation = _Rotation()
140 S3_angle = _np.pi - alpha
141 add_segment(S3_center, S3_axis_rotation, S3_radius, S3_angle, n_el)
142 add_line(mesh.nodes[-1].coordinates, [-bottom_width, height, 0], 2 * n_el)
144 if S2:
145 S2_radius = (height * 0.5 / _np.tan(alpha) + top_width) * _np.tan(alpha * 0.5)
146 S2_center = [0, height * 0.5, 0] - S2_radius * _np.array([0, 1, 0])
147 S2_angle = _np.pi - alpha
148 S2_axis_rotation = _Rotation([0, 0, 1], _np.pi)
149 add_segment(S2_center, S2_axis_rotation, S2_radius, S2_angle, 2 * n_el)
150 add_line(mesh.nodes[-1].coordinates, [-top_width, 0, 0], 2 * n_el)
152 mesh.translate([width, -height / 2, 0])
153 return mesh
156def create_stent_column(
157 beam_class, material, width, height, n_height, n_el=1, **kwargs
158):
159 """Create a column of completed cells.
161 A completed cell consists of one cell, that is created with the create cell
162 function and it's reflection.
164 Args
165 ----
166 beam_class: Beam
167 Class that will be used to create the beam elements.
168 material: Material
169 Material for the beam.
170 width: float
171 Width of the total cell.
172 height: float
173 Height of the total cell.
174 n_height: int
175 The number of cells in the column
176 n_el: int
177 Number of elements per beam line.
178 ( these variables are described in a file )
180 Return
181 ----
182 mesh: Mesh
183 A mesh with this structure.
184 """
185 mesh_column = _Mesh()
186 for i in range(n_height):
187 S1 = True
188 S2 = True
189 S3 = True
190 if i == 0:
191 S1 = False
192 S2 = False
193 if i == 1:
194 S2 = False
195 if i == n_height - 1:
196 S1 = False
197 S3 = False
199 if i == n_height - 2:
200 S3 = False
201 unit_cell = create_stent_cell(
202 beam_class,
203 material,
204 width,
205 height,
206 S1=S1,
207 S2=S2,
208 S3=S3,
209 n_el=n_el,
210 **kwargs,
211 )
212 unit_cell.translate([0, i * height, 0])
213 mesh_column.add(unit_cell)
215 column_copy = mesh_column.copy()
216 column_copy.reflect(normal_vector=[1, 0, 0], origin=[width, 0, 0])
217 mesh_column.add(column_copy)
219 return mesh_column
222def create_beam_mesh_stent_flat(
223 beam_class, material, width_flat, height_flat, n_height, n_column, n_el=1, **kwargs
224):
225 """Create a flat stent structure on the x-y plane.
227 Args
228 ----
229 beam_class: Beam
230 Class that will be used to create the beam elements.
231 material: Material
232 Material for the beam.
233 width_flat: float
234 The width of the flat structure.
235 height_flat: float
236 The height of the flat structure.
237 n_height: int
238 The number of cells in y direction.
239 n_column: int
240 The number of columns in x direction.
241 n_el: int
242 Number of elements per beam line.
244 Return
245 ----
246 mesh: Mesh
247 A mesh with this structure
248 """
249 mesh_flat = _Mesh()
250 width = width_flat / n_column / 2
251 height = height_flat / n_height
252 column_mesh = create_stent_column(
253 beam_class, material, width, height, n_height, n_el=n_el, **kwargs
254 )
255 for i in range(n_column):
256 column_copy = column_mesh.copy()
257 column_copy.translate([2 * width * i, 0, 0])
258 mesh_flat.add(column_copy)
260 for i in range(n_column // 2):
261 for j in range(n_height - 1):
262 _create_beam_mesh_line(
263 mesh_flat,
264 beam_class,
265 material,
266 [4 * i * width, j * height, 0],
267 [4 * i * width, (j + 1) * height, 0],
268 n_el=2 * n_el,
269 )
270 _create_beam_mesh_line(
271 mesh_flat,
272 beam_class,
273 material,
274 [(4 * i + 2) * width, 0, 0],
275 [(4 * i + 2) * width, height, 0],
276 n_el=2 * n_el,
277 )
278 return mesh_flat
281def create_beam_mesh_stent(
282 mesh,
283 beam_class,
284 material,
285 length,
286 diameter,
287 n_axis,
288 n_circumference,
289 **kwargs,
290):
291 """Create a stent structure around cylinder, The cylinder axis will be the z-axis.
293 Args
294 ----
295 mesh: Mesh
296 Mesh that the stent will be added to.
297 beam_class: Beam
298 Class that will be used to create the beam elements.
299 material: Material
300 Material for the beam.
301 length: float
302 The length of this stent.
303 diameter: float
304 The diameter of the stent's cross section.
305 n_axis: int
306 Number of cells in axial-direction.
307 n_circumference: int
308 Number of cells around the diameter.
309 ( these variables are described in a file )
311 Return
312 ----
313 return_set: GeometryName
314 Set with nodes on the top, bottom boundaries. Those
315 sets only contains end nodes of lines, not the middle ones.
316 The set 'all' contains all nodes.
317 """
318 # Only allow even number of columns.
319 if n_circumference % 2 == 1:
320 raise ValueError("has to be even even number!")
322 # Set the Parameter for other functions
323 height_flat = length
324 width_flat = _np.pi * diameter
325 n_height = n_axis
326 n_column = n_circumference
328 mesh_stent = create_beam_mesh_stent_flat(
329 beam_class, material, width_flat, height_flat, n_height, n_column, **kwargs
330 )
331 mesh_stent.rotate(_Rotation([1, 0, 0], _np.pi / 2))
332 mesh_stent.rotate(_Rotation([0, 0, 1], _np.pi / 2))
333 mesh_stent.translate([diameter / 2, 0, 0])
334 mesh_stent.wrap_around_cylinder()
336 # List of nodes from the stent that are candidates for connections.
337 stent_nodes = [node for node in mesh_stent.nodes if node.is_end_node]
339 # Add connections for the nodes with same positions.
340 mesh_stent.couple_nodes(nodes=stent_nodes)
342 # Get min and max nodes of the honeycomb.
343 min_max_nodes = _get_min_max_nodes(stent_nodes)
345 # Return the geometry set.
346 return_set = _GeometryName()
347 return_set["top"] = min_max_nodes["z_max"]
348 return_set["bottom"] = min_max_nodes["z_min"]
349 return_set["all"] = _GeometrySet(mesh_stent.elements)
351 mesh.add_mesh(mesh_stent)
353 return return_set