Coverage for src/meshpy/mesh_creation_functions/beam_stent.py: 98%
108 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 stent according to Auricchio 2012."""
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_arc_segment_via_rotation as _create_beam_mesh_arc_segment_via_rotation,
32)
33from meshpy.mesh_creation_functions.beam_basic_geometry import (
34 create_beam_mesh_line as _create_beam_mesh_line,
35)
36from meshpy.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. This cell is on the x-y plane.
55 Args
56 ----
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 the total cell.
63 height: float
64 Height of the total cell.
65 fac_bottom: the ratio of the bottom's width to the cell's width
66 fac_neck: the ratio of the neck's width to the cell's width
67 fac_radius: the ratio of the S1's radius to the cell's width
68 alpha: radiant
69 The angle between the lines and horizontal line
70 n_el: int
71 Number of elements per beam line.
72 S1, S2, S3: bool
73 This check weather the curve S1, S2 or S3 will be created.
74 If the cell is on bottom of the stent flat S1 and S2 won't
75 be created. If the cell is on top of the flat S1 and S3
76 won't be created
78 ( these variables are described in a file )
80 Return
81 ----
82 mesh: Mesh
83 A mesh with this structure
84 """
86 mesh = _Mesh()
88 def add_line(pointa, pointb, n_el_line):
89 """Shortcut to add line."""
90 return _create_beam_mesh_line(
91 mesh, beam_class, material, pointa, pointb, n_el=n_el_line
92 )
94 def add_segment(center, axis_rotation, radius, angle, n_el_segment):
95 """Shortcut to add arc segment."""
96 return _create_beam_mesh_arc_segment_via_rotation(
97 mesh,
98 beam_class,
99 material,
100 center,
101 axis_rotation,
102 radius,
103 angle,
104 n_el=n_el_segment,
105 )
107 neck_width = width * fac_neck
108 bottom_width = width * fac_bottom
109 top_width = 2 * neck_width - bottom_width
110 radius = width * fac_radius
112 if S1:
113 neck_point = _np.array([-neck_width, height * 0.5, 0])
114 d = (height * 0.5 / _np.tan(alpha) + bottom_width - neck_width) / _np.sin(alpha)
115 CM = _np.array([-_np.sin(alpha), -_np.cos(alpha), 0]) * (d - radius)
116 MO = _np.array([_np.cos(alpha), -_np.sin(alpha), 0]) * _np.sqrt(
117 radius**2 - (d - radius) ** 2
118 )
119 S1_angle = _np.pi / 2 + _np.arcsin((d - radius) / radius)
120 S1_center1 = CM + MO + neck_point
121 S1_axis_rotation1 = _Rotation([0, 0, 1], 2 * _np.pi - S1_angle - alpha)
122 add_segment(S1_center1, S1_axis_rotation1, radius, S1_angle, n_el)
123 add_line([-bottom_width, 0, 0], mesh.nodes[-1].coordinates, 2 * n_el)
125 S1_center2 = 2 * neck_point - S1_center1
126 S1_axis_rotation2 = _Rotation([0, 0, 1], _np.pi - alpha - S1_angle)
128 add_segment(S1_center2, S1_axis_rotation2, radius, S1_angle, n_el)
129 add_line(
130 mesh.nodes[-1].coordinates, 2 * neck_point - [-bottom_width, 0, 0], 2 * n_el
131 )
133 if S3:
134 S3_radius = (width - bottom_width + height * 0.5 / _np.tan(alpha)) * _np.tan(
135 alpha / 2
136 )
137 S3_center = [-width, height * 0.5, 0] + S3_radius * _np.array([0, 1, 0])
138 S3_axis_rotation = _Rotation()
139 S3_angle = _np.pi - alpha
140 add_segment(S3_center, S3_axis_rotation, S3_radius, S3_angle, n_el)
141 add_line(mesh.nodes[-1].coordinates, [-bottom_width, height, 0], 2 * n_el)
143 if S2:
144 S2_radius = (height * 0.5 / _np.tan(alpha) + top_width) * _np.tan(alpha * 0.5)
145 S2_center = [0, height * 0.5, 0] - S2_radius * _np.array([0, 1, 0])
146 S2_angle = _np.pi - alpha
147 S2_axis_rotation = _Rotation([0, 0, 1], _np.pi)
148 add_segment(S2_center, S2_axis_rotation, S2_radius, S2_angle, 2 * n_el)
149 add_line(mesh.nodes[-1].coordinates, [-top_width, 0, 0], 2 * n_el)
151 mesh.translate([width, -height / 2, 0])
152 return mesh
155def create_stent_column(
156 beam_class, material, width, height, n_height, n_el=1, **kwargs
157):
158 """Create a column of completed cells. A completed cell consists of one
159 cell, that is created with the create cell function and it's reflection.
161 Args
162 ----
163 beam_class: Beam
164 Class that will be used to create the beam elements.
165 material: Material
166 Material for the beam.
167 width: float
168 Width of the total cell.
169 height: float
170 Height of the total cell.
171 n_height: int
172 The number of cells in the column
173 n_el: int
174 Number of elements per beam line.
175 ( these variables are described in a file )
177 Return
178 ----
179 mesh: Mesh
180 A mesh with this structure.
181 """
183 mesh_column = _Mesh()
184 for i in range(n_height):
185 S1 = True
186 S2 = True
187 S3 = True
188 if i == 0:
189 S1 = False
190 S2 = False
191 if i == 1:
192 S2 = False
193 if i == n_height - 1:
194 S1 = False
195 S3 = False
197 if i == n_height - 2:
198 S3 = False
199 unit_cell = create_stent_cell(
200 beam_class,
201 material,
202 width,
203 height,
204 S1=S1,
205 S2=S2,
206 S3=S3,
207 n_el=n_el,
208 **kwargs,
209 )
210 unit_cell.translate([0, i * height, 0])
211 mesh_column.add(unit_cell)
213 column_copy = mesh_column.copy()
214 column_copy.reflect(normal_vector=[1, 0, 0], origin=[width, 0, 0])
215 mesh_column.add(column_copy)
217 return mesh_column
220def create_beam_mesh_stent_flat(
221 beam_class, material, width_flat, height_flat, n_height, n_column, n_el=1, **kwargs
222):
223 """Create a flat stent structure on the x-y plane.
225 Args
226 ----
227 beam_class: Beam
228 Class that will be used to create the beam elements.
229 material: Material
230 Material for the beam.
231 width_flat: float
232 The width of the flat structure.
233 height_flat: float
234 The height of the flat structure.
235 n_height: int
236 The number of cells in y direction.
237 n_column: int
238 The number of columns in x direction.
239 n_el: int
240 Number of elements per beam line.
242 Return
243 ----
244 mesh: Mesh
245 A mesh with this structure
246 """
248 mesh_flat = _Mesh()
249 width = width_flat / n_column / 2
250 height = height_flat / n_height
251 column_mesh = create_stent_column(
252 beam_class, material, width, height, n_height, n_el=n_el, **kwargs
253 )
254 for i in range(n_column):
255 column_copy = column_mesh.copy()
256 column_copy.translate([2 * width * i, 0, 0])
257 mesh_flat.add(column_copy)
259 for i in range(n_column // 2):
260 for j in range(n_height - 1):
261 _create_beam_mesh_line(
262 mesh_flat,
263 beam_class,
264 material,
265 [4 * i * width, j * height, 0],
266 [4 * i * width, (j + 1) * height, 0],
267 n_el=2 * n_el,
268 )
269 _create_beam_mesh_line(
270 mesh_flat,
271 beam_class,
272 material,
273 [(4 * i + 2) * width, 0, 0],
274 [(4 * i + 2) * width, height, 0],
275 n_el=2 * n_el,
276 )
277 return mesh_flat
280def create_beam_mesh_stent(
281 mesh,
282 beam_class,
283 material,
284 length,
285 diameter,
286 n_axis,
287 n_circumference,
288 add_sets=False,
289 **kwargs,
290):
291 """Create a stent structure around cylinder, The cylinder axis will be the
292 z-axis.
294 Args
295 ----
296 mesh: Mesh
297 Mesh that the stent will be added to.
298 beam_class: Beam
299 Class that will be used to create the beam elements.
300 material: Material
301 Material for the beam.
302 length: float
303 The length of this stent.
304 diameter: float
305 The diameter of the stent's cross section.
306 n_axis: int
307 Number of cells in axial-direction.
308 n_circumference: int
309 Number of cells around the diameter.
310 ( these variables are described in a file )
311 add_sets: bool
312 If this is true the sets are added to the mesh and then displayed
313 n eventual VTK output, even if they are not used for a boundary
314 condition or coupling.
316 Return
317 ----
318 return_set: GeometryName
319 Set with nodes on the top, bottom boundaries. Those
320 sets only contains end nodes of lines, not the middle ones.
321 The set 'all' contains all nodes.
322 """
324 # Only allow even number of columns.
325 if n_circumference % 2 == 1:
326 raise ValueError("has to be even even number!")
328 # Set the Parameter for other functions
329 height_flat = length
330 width_flat = _np.pi * diameter
331 n_height = n_axis
332 n_column = n_circumference
334 mesh_stent = create_beam_mesh_stent_flat(
335 beam_class, material, width_flat, height_flat, n_height, n_column, **kwargs
336 )
337 mesh_stent.rotate(_Rotation([1, 0, 0], _np.pi / 2))
338 mesh_stent.rotate(_Rotation([0, 0, 1], _np.pi / 2))
339 mesh_stent.translate([diameter / 2, 0, 0])
340 mesh_stent.wrap_around_cylinder()
342 # List of nodes from the stent that are candidates for connections.
343 stent_nodes = [node for node in mesh_stent.nodes if node.is_end_node]
345 # Add connections for the nodes with same positions.
346 mesh_stent.couple_nodes(nodes=stent_nodes)
348 # Get min and max nodes of the honeycomb.
349 min_max_nodes = _get_min_max_nodes(stent_nodes)
351 # Return the geometry set.
352 return_set = _GeometryName()
353 return_set["top"] = min_max_nodes["z_max"]
354 return_set["bottom"] = min_max_nodes["z_min"]
355 return_set["all"] = _GeometrySet(mesh_stent.elements)
357 mesh.add_mesh(mesh_stent)
359 if add_sets:
360 mesh.add(return_set)
361 return return_set