Coverage for src/meshpy/mesh_creation_functions/beam_fibers_in_rectangle.py: 95%
55 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 generate multiple parallel fibers within a
23rectangle.
25This can for example be used to create fiber reinforced composite
26plates.
27"""
29import numpy as _np
31from meshpy.core.geometry_set import GeometryName as _GeometryName
32from meshpy.core.geometry_set import GeometrySet as _GeometrySet
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 check_node_by_coordinate as _check_node_by_coordinate
39def _intersect_line_with_rectangle(
40 length, width, start_line, direction_line, fail_if_no_intersection=True
41):
42 """Calculate the intersection points between a line and a rectangle.
44 Args
45 ----
46 length: scalar
47 Rectangle length in x direction.
48 width: scalar
49 Rectangle width in y direction.
50 start_line: 2d-list
51 Point on the line.
52 direction_line: 2d-list
53 Direction of the line.
54 fail_if_no_intersection: bool
55 If this is true and no intersections are found, an error will be
56 thrown.
58 Return
59 ----
60 (start, end, projection_found)
61 start: 2D vector
62 Start point of intersected line.
63 end: 2D vector
64 End point of intersected line.
65 projection_found: bool
66 True if intersection is valid.
67 """
69 # Convert the input values to np.arrays.
70 start_line = _np.asarray(start_line)
71 direction_line = _np.asarray(direction_line)
73 # Set definition for the boundary lines of the rectangle. The director is
74 # chosen in a way, that the values [0, 1] for the line parameters alpha are
75 # valid.
76 # boundary_lines = [..., [start, dir], ...]
77 boundary_lines = [
78 [[0, 0], [length, 0]],
79 [[0, 0], [0, width]],
80 [[0, width], [length, 0]],
81 [[length, 0], [0, width]],
82 ]
83 # Convert to numpy arrays.
84 boundary_lines = [
85 [_np.array(item) for item in boundary] for boundary in boundary_lines
86 ]
88 # Loop over the boundaries.
89 alpha_list = []
90 for start_boundary, direction_boundary in boundary_lines:
91 # Set up the linear system to solve the intersection problem.
92 A = _np.transpose(_np.array([direction_line, -direction_boundary]))
94 # Check if the system is solvable.
95 if _np.abs(_np.linalg.det(A)) > 1e-10:
96 alpha = _np.linalg.solve(A, start_boundary - start_line)
97 if 0 <= alpha[1] and alpha[1] <= 1:
98 alpha_list.append(alpha[0])
100 # Check that intersections were found.
101 if len(alpha_list) < 2:
102 if fail_if_no_intersection:
103 raise ValueError("No intersections found!")
104 return (None, None, False)
106 # Return the start and end point on the line.
107 return (
108 start_line + min(alpha_list) * direction_line,
109 start_line + max(alpha_list) * direction_line,
110 True,
111 )
114def create_fibers_in_rectangle(
115 mesh,
116 beam_class,
117 material,
118 length,
119 width,
120 angle,
121 fiber_normal_distance,
122 fiber_element_length,
123 *,
124 reference_point=None,
125 fiber_element_length_min=None,
126):
127 """Create multiple fibers in a rectangle.
129 Args
130 ----
131 mesh: Mesh
132 Mesh that the fibers will be added to.
133 beam_class: Beam
134 Class that will be used to create the beam elements.
135 material: Material
136 Material for the beam.
137 length: float
138 Length of the rectangle in x direction (starting at x=0)
139 width: float
140 Width of the rectangle in y direction (starting at y=0)
141 angle: float
142 Angle of the fibers in degree.
143 fiber_normal_distance: float
144 Normal distance between the parallel fibers.
145 fiber_element_length: float
146 Length of a single beam element. In general it will not be possible to
147 exactly achieve this length.
148 reference_point: [float, float]
149 Specify a single point inside the rectangle that one of the fibers will pass through.
150 Per default the reference point is in the middle of the rectangle.
151 fiber_element_length_min: float
152 Minimum fiber length. If a fiber is shorter than this value, it will not be created.
153 The default value is half of fiber_element_length.
154 """
156 if reference_point is None:
157 reference_point = 0.5 * _np.array([length, width])
158 else:
159 if (
160 reference_point[0] < 0.0
161 or reference_point[0] > length
162 or reference_point[1] < 0.0
163 or reference_point[1] > width
164 ):
165 raise ValueError("The reference point has to lie within the rectangle")
167 if fiber_element_length_min is None:
168 fiber_element_length_min = 0.5 * fiber_element_length
169 elif fiber_element_length_min < 0.0:
170 raise ValueError("fiber_element_length_min must be positive!")
172 # Get the fiber angle in rad.
173 fiber_angle = angle * _np.pi / 180.0
174 sin = _np.sin(fiber_angle)
175 cos = _np.cos(fiber_angle)
177 # Direction and normal vector of the fibers.
178 fiber_direction = _np.array([cos, sin])
179 fiber_normal = _np.array([-sin, cos])
181 # Get an upper bound of the number of fibers in this layer.
182 diagonal = _np.sqrt(length**2 + width**2)
183 fiber_n_max = int(_np.ceil(diagonal / fiber_normal_distance)) + 1
185 # Go in both directions from the start point.
186 for direction_sign, n_start in [[-1, 1], [1, 0]]:
187 # Create a fiber as long as an intersection is found.
188 for i_fiber in range(n_start, fiber_n_max):
189 # Get the start and end point of the line.
190 start, end, projection_found = _intersect_line_with_rectangle(
191 length,
192 width,
193 reference_point
194 + fiber_normal * i_fiber * fiber_normal_distance * direction_sign,
195 fiber_direction,
196 fail_if_no_intersection=False,
197 )
199 if projection_found:
200 # Calculate the length of the line.
201 fiber_length = _np.linalg.norm(end - start)
203 # Create the beams if the length is not smaller than the fiber
204 # distance.
205 if fiber_length >= fiber_element_length_min:
206 # Calculate the number of elements in this fiber.
207 fiber_nel = int(_np.round(fiber_length / fiber_element_length))
208 fiber_nel = _np.max([fiber_nel, 1])
209 _create_beam_mesh_line(
210 mesh,
211 beam_class,
212 material,
213 _np.append(start, 0.0),
214 _np.append(end, 0.0),
215 n_el=fiber_nel,
216 )
217 else:
218 # The current search position is already outside of the rectangle, no need to continue.
219 break
221 return_set = _GeometryName()
222 return_set["north"] = _GeometrySet(
223 mesh.get_nodes_by_function(_check_node_by_coordinate, 1, width),
224 )
225 return_set["east"] = _GeometrySet(
226 mesh.get_nodes_by_function(_check_node_by_coordinate, 0, length),
227 )
228 return_set["south"] = _GeometrySet(
229 mesh.get_nodes_by_function(_check_node_by_coordinate, 1, 0)
230 )
231 return_set["west"] = _GeometrySet(
232 mesh.get_nodes_by_function(_check_node_by_coordinate, 0, 0)
233 )
234 return_set["all"] = _GeometrySet(mesh.elements)
235 return return_set