Coverage for src / sdynpy / demo / frame_wing.py: 0%
144 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 21:21 +0000
1"""Generates a simplified model of the frame/wing structure from the SEM/IMAC Substructuring Focus
2Group. It's not a great model because it uses beam finite elements, but it should be able to be
3used for example problems where the dynamics are somewhat similar to those in the real model.
4"""
6import numpy as np
8# from ..core.sdynpy_coordinate import (from_nodelist as _from_nodelist,
9# coordinate_array as _coordinate_array)
10# from ..fem.sdynpy_beam import beamkm as _beamkm, rect_beam_props as _rect_beam_props
11# from ..core.sdynpy_geometry import (Geometry as _Geometry,
12# node_array as _node_array,
13# traceline_array as _traceline_array,
14# coordinate_system_array as _coordinate_system_array)
15# from ..core.sdynpy_system import System as _System
16from sdynpy.core.sdynpy_coordinate import from_nodelist as _from_nodelist
17from sdynpy.fem.sdynpy_beam import beamkm as _beamkm, rect_beam_props as _rect_beam_props
18from sdynpy.core.sdynpy_geometry import (
19 Geometry as _Geometry,
20 node_array as _node_array,
21 traceline_array as _traceline_array,
22 coordinate_system_array as _coordinate_system_array,
23 element_array as _element_array,
24)
25from sdynpy.core.sdynpy_system import System as _System
28def _build_models():
29 # Define in inches, which is how the drawing is built, we will convert
30 # to meters later.
31 frame_thickness = 0.5 # inches
32 frame_length = 7.75 * 2 # inches
33 frame_width = 2.75 * 2 # inches
34 frame_bays = 4
35 elem_per_bay_length = 3
36 elem_per_bay_width = 4
37 wing_width = 22 # inches
38 wing_thickness = 0.25 * 1.3 # inches
39 wing_length = 4.385 # inches
40 wing_attachment_bay = 2
41 wing_attachment_offset = wing_length / 2 - 1.9376
42 E = 69e9 # [N/m^2],
43 nu = 0.33 # [-],
44 rho = 2830 # [kg/m^3]
45 IN2M = 0.0254
47 # Create the frame
48 node_x = np.linspace(-frame_length / 2, frame_length / 2, frame_bays * elem_per_bay_length + 1)
49 node_y = np.linspace(-frame_width / 2, frame_width / 2, elem_per_bay_width + 1)
51 node_indices = (
52 [(i, 0) for i in range(len(node_x))]
53 + [(i, len(node_y) - 1) for i in range(len(node_x))]
54 + [
55 (i * (elem_per_bay_length), j)
56 for i in range(len(node_x) // elem_per_bay_length + 1)
57 for j in range(len(node_y))
58 ]
59 )
60 node_indices = sorted(list(set(node_indices)))
62 node_ids = np.array([(i + 1) * 10 + j + 1 for i, j in node_indices])
64 # Here we convert to meters
65 frame_node_coordinates = (
66 np.array([[node_x[i], node_y[j], 0] for i, j in node_indices]) * IN2M
67 ) # Meters
69 # Create connectivity arrays
70 frame_connectivity = []
71 for i in range(len(node_x)):
72 for j in range(len(node_y)):
73 if not (i, j) in node_indices:
74 continue
75 idx0 = node_indices.index((i, j))
76 if (i + 1, j) in node_indices:
77 idx1 = node_indices.index((i + 1, j))
78 frame_connectivity.append((idx0, idx1))
79 if (i, j + 1) in node_indices:
80 idx1 = node_indices.index((i, j + 1))
81 frame_connectivity.append((idx0, idx1))
82 frame_connectivity = np.array(frame_connectivity)
84 # Create the beam models
85 frame_bend_direction_1 = np.array(
86 (
87 np.zeros(len(frame_connectivity)),
88 np.zeros(len(frame_connectivity)),
89 np.ones(len(frame_connectivity)),
90 )
91 ).T
93 width = frame_thickness * IN2M # Meters
94 height = frame_thickness * IN2M # Meters
95 frame_mat_props = _rect_beam_props(E, rho, nu, width, height, len(frame_connectivity))
96 K, M = _beamkm(
97 frame_node_coordinates, frame_connectivity, frame_bend_direction_1, **frame_mat_props
98 )
99 coordinates = _from_nodelist(node_ids, directions=[1, 2, 3, 4, 5, 6])
100 frame_system = _System(coordinates, M, K)
102 frame_geometry = _Geometry(
103 _node_array(node_ids, frame_node_coordinates),
104 _coordinate_system_array(1),
105 _traceline_array(
106 np.arange(frame_connectivity.shape[0]) + 1,
107 connectivity=[arr for arr in node_ids[frame_connectivity]],
108 ),
109 )
110 frame_modes = frame_system.eigensolution(num_modes=50)
111 frame_geometry.plot_shape(frame_modes)
113 # Create the wing
115 def structured_tri_mesh(
116 width, height, dx, dy=None, origin=(0.0, 0.0), diag="alternating", return_edges=False
117 ):
118 """
119 Structured triangular mesh on an axis-aligned rectangle using a rect grid + diagonal splits.
121 Returns
122 pts: (N,2) array
123 tri: (M,3) array of indices into pts
124 edges (optional): (E,2) array of unique undirected edges (i<j)
125 """
126 if dy is None:
127 dy = dx
128 x0, y0 = origin
130 nx = int(np.ceil(width / dx))
131 ny = int(np.ceil(height / dy))
133 xs = x0 + np.linspace(0.0, width, nx + 1)
134 ys = y0 + np.linspace(0.0, height, ny + 1)
135 X, Y = np.meshgrid(xs, ys, indexing="xy")
136 pts = np.column_stack([X.ravel(), Y.ravel()])
138 def vid(i, j): # i in [0..nx], j in [0..ny]
139 return j * (nx + 1) + i
141 tris = []
142 edges = []
143 for j in range(ny):
144 for i in range(nx):
145 v00 = vid(i, j)
146 v10 = vid(i + 1, j)
147 v01 = vid(i, j + 1)
148 v11 = vid(i + 1, j + 1)
150 edges.append([v00, v10])
151 edges.append([v00, v11])
152 edges.append([v00, v01])
153 edges.append([v10, v11])
154 edges.append([v10, v01])
155 edges.append([v01, v11])
157 if diag == "same" or diag == "x":
158 tris.append([v00, v10, v11])
159 tris.append([v00, v11, v01])
160 elif diag == "y":
161 tris.append([v00, v10, v01])
162 tris.append([v10, v11, v01])
163 elif diag == "alternating":
164 if (i + j) % 2 == 0:
165 tris.append([v00, v10, v11])
166 tris.append([v00, v11, v01])
167 else:
168 tris.append([v00, v10, v01])
169 tris.append([v10, v11, v01])
170 else:
171 raise ValueError("diag must be one of: same, alternating, x, y")
173 tri = np.asarray(tris, dtype=int)
175 if not return_edges:
176 return pts, tri
178 # Build all edges from triangles, canonicalize (i<j), then unique
179 all_edges = np.array(edges)
180 all_edges = np.sort(all_edges, axis=1) # undirected canonical form
181 edges = np.unique(all_edges, axis=0) # unique rows
183 return pts, tri, edges
185 wing_coordinates, wing_tris, wing_edges = structured_tri_mesh(
186 wing_length * IN2M,
187 wing_width * IN2M,
188 dy=frame_width / elem_per_bay_width * IN2M,
189 dx=frame_length / frame_bays / elem_per_bay_length * IN2M * 1.2,
190 origin=np.array(
191 (
192 -frame_length / 2
193 + frame_length / frame_bays * (wing_attachment_bay - 1)
194 - wing_attachment_offset,
195 -wing_width / 2,
196 )
197 )
198 * IN2M,
199 return_edges=True,
200 )
202 wing_coordinates = np.concatenate(
203 (
204 wing_coordinates,
205 (frame_thickness / 2 + wing_thickness / 2)
206 * IN2M
207 * np.ones((wing_coordinates.shape[0], 1)),
208 ),
209 axis=-1,
210 )
211 wing_ids = np.arange(wing_coordinates.shape[0]) + 1001
213 wing_geometry = _Geometry(
214 _node_array(wing_ids, wing_coordinates),
215 _coordinate_system_array(1),
216 # _traceline_array(np.arange(wing_edges.shape[0])+1,
217 # connectivity = [arr for arr in wing_ids[wing_edges]]),
218 element=_element_array(
219 np.arange(wing_tris.shape[0]) + 1,
220 type=41,
221 connectivity=[arr for arr in wing_ids[wing_tris]],
222 ),
223 )
225 # Find the nodes in the second bay to make them coincident
226 bay_length = frame_length / frame_bays
227 bay_xmin = (
228 -frame_length / 2
229 + (wing_attachment_bay - 1) * bay_length
230 - bay_length / (elem_per_bay_length * 2)
231 )
232 bay_xmax = (
233 -frame_length / 2
234 + (wing_attachment_bay) * bay_length
235 + bay_length / (elem_per_bay_length * 2)
236 )
237 # Get nodes in this range
238 nodes_in_bay = frame_geometry.node.id[
239 (frame_geometry.node.coordinate[:, 0] > bay_xmin * IN2M)
240 & (frame_geometry.node.coordinate[:, 0] < bay_xmax * IN2M)
241 ]
242 attachment_coordinates = frame_geometry.node(nodes_in_bay).coordinate
243 equivalent_wing_indices = np.array(
244 [
245 np.where(
246 wing_geometry.node.id
247 == wing_geometry.node_by_global_position(attachment_coordinate).id
248 )[0][0]
249 for attachment_coordinate in attachment_coordinates
250 ]
251 )
252 wing_geometry.node.coordinate[equivalent_wing_indices, :2] = attachment_coordinates[:, :2]
254 # Now build complete mass and stiffness matrices
255 wing_node_coordinates = wing_geometry.node.coordinate
256 wing_connectivity = wing_edges
257 wing_bend_direction_1 = np.array(
258 (
259 np.zeros(len(wing_connectivity)),
260 np.zeros(len(wing_connectivity)),
261 np.ones(len(wing_connectivity)),
262 )
263 ).T
264 width = wing_thickness * IN2M # Meters
265 height = wing_thickness * IN2M # Meters
266 wing_mat_props = _rect_beam_props(E, rho, nu, width, height, len(wing_connectivity))
267 K, M = _beamkm(
268 wing_node_coordinates, wing_connectivity, wing_bend_direction_1, **wing_mat_props
269 )
270 coordinates = _from_nodelist(wing_geometry.node.id, directions=[1, 2, 3, 4, 5, 6])
271 wing_system = _System(coordinates, M, K)
273 wing_modes = wing_system.eigensolution(num_modes=50)
274 wing_geometry.plot_shape(wing_modes)
276 # Now let's create a coupled system
277 equivalent_frame_indices = np.array(
278 [
279 np.where(
280 frame_geometry.node.id
281 == frame_geometry.node_by_global_position(attachment_coordinate).id
282 )[0][0]
283 for attachment_coordinate in attachment_coordinates
284 ]
285 )
287 connection_connectivity = []
288 for frame_index, wing_index in zip(equivalent_frame_indices, equivalent_wing_indices):
289 connection_connectivity.append([frame_index, wing_index + frame_node_coordinates.shape[0]])
291 connection_connectivity = np.array(connection_connectivity)
292 connection_bend_direction_1 = np.array(
293 (
294 np.ones(len(connection_connectivity)),
295 np.zeros(len(connection_connectivity)),
296 np.zeros(len(connection_connectivity)),
297 )
298 ).T
299 width = frame_thickness * IN2M # Meters
300 height = frame_thickness * IN2M # Meters
301 connection_mat_props = _rect_beam_props(E, rho, nu, width, height, len(connection_connectivity))
303 full_connectivity = np.concatenate(
304 (
305 frame_connectivity,
306 wing_connectivity + frame_node_coordinates.shape[0],
307 connection_connectivity,
308 )
309 )
310 full_coordinates = np.concatenate((frame_node_coordinates, wing_node_coordinates))
311 full_bend_1_direction = np.concatenate(
312 (frame_bend_direction_1, wing_bend_direction_1, connection_bend_direction_1)
313 )
314 full_mat_props = {}
315 for mat_prop in (frame_mat_props, wing_mat_props, connection_mat_props):
316 for key, value in mat_prop.items():
317 if key not in full_mat_props:
318 full_mat_props[key] = []
319 full_mat_props[key] = np.concatenate((full_mat_props[key], value))
320 K, M = _beamkm(full_coordinates, full_connectivity, full_bend_1_direction, **full_mat_props)
321 coordinates = np.concatenate(
322 (
323 _from_nodelist(frame_geometry.node.id, directions=[1, 2, 3, 4, 5, 6]),
324 _from_nodelist(wing_geometry.node.id, directions=[1, 2, 3, 4, 5, 6]),
325 )
326 )
327 full_system = _System(coordinates, M, K)
328 full_modes = full_system.eigensolution(num_modes=50)
329 full_geometry = frame_geometry + wing_geometry
331 full_geometry.plot_shape(full_modes)
333 return full_geometry, full_system, frame_geometry, frame_system, wing_geometry, wing_system
336geometry, system, frame_geometry, frame_system, wing_geometry, wing_system = _build_models()