Coverage for src/sdynpy/core/sdynpy_system.py: 39%
710 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-18 21:32 +0000
1# -*- coding: utf-8 -*-
2"""
3Defines a system of matrices representing a structure.
5The System consists of mass, stiffness, and (if necessary) damping matrices.
6The System also contains a CoordinateArray defining the degrees of freedom of
7the System, as well as a transformation that takes the System from its internal
8state degrees of freedom to physical degrees of freedom.
9"""
10"""
11Copyright 2022 National Technology & Engineering Solutions of Sandia,
12LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
13Government retains certain rights in this software.
15This program is free software: you can redistribute it and/or modify
16it under the terms of the GNU General Public License as published by
17the Free Software Foundation, either version 3 of the License, or
18(at your option) any later version.
20This program is distributed in the hope that it will be useful,
21but WITHOUT ANY WARRANTY; without even the implied warranty of
22MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23GNU General Public License for more details.
25You should have received a copy of the GNU General Public License
26along with this program. If not, see <https://www.gnu.org/licenses/>.
27"""
29import numpy as np
30from .sdynpy_coordinate import CoordinateArray, from_nodelist, outer_product, coordinate_array
31from ..fem.sdynpy_beam import beamkm, rect_beam_props
32from ..fem.sdynpy_shaker import Shaker4DoF
33from ..signal_processing import frf as spfrf
34from ..signal_processing import generator
35from scipy.linalg import eigh, block_diag, null_space, eig
36from scipy.signal import lsim, StateSpace, resample, butter, filtfilt
37import copy
38import netCDF4 as nc4
39import matplotlib.pyplot as plt
40import warnings
41from pathlib import Path
44class System:
45 """Matrix Equations representing a Structural Dynamics System"""
47 def __init__(self, coordinate: CoordinateArray, mass, stiffness, damping=None,
48 transformation=None, enforce_symmetry = True):
49 """
50 Create a system representation including mass, stiffness, damping, and
51 transformation matrices.
53 Parameters
54 ----------
55 coordinate : CoordinateArray
56 Physical degrees of freedom in the system.
57 mass : np.ndarray
58 2D array consisting of the mass matrix of the system
59 stiffness : np.ndarray
60 2D array consisting of the stiffness matrix of the system
61 damping : np.ndarray, optional
62 2D array consisting of the damping matrix of the system. If not
63 specified, the damping will be zero.
64 transformation : np.ndarray, optional
65 A transformation between internal "state" degrees of freedom and
66 the physical degrees of freedom defined in `coordinate`. The
67 default transformation is the identity matrix.
68 enforce_symmetry : bool, optional
69 If True, raise a ValueError if the matrices provided are not
70 symmetric.
72 Raises
73 ------
74 ValueError
75 If inputs are improperly sized
77 Returns
78 -------
79 None.
81 """
82 mass = np.atleast_2d(np.array(mass))
83 stiffness = np.atleast_2d(np.array(stiffness))
84 if not mass.shape == stiffness.shape:
85 raise ValueError('Mass and Stiffness matrices must be the same shape')
86 if not mass.ndim == 2 or (mass.shape[0] != mass.shape[1]):
87 raise ValueError('Mass should be a 2D, square array')
88 if not stiffness.ndim == 2 or (stiffness.shape[0] != stiffness.shape[1]):
89 raise ValueError('Stiffness should be a 2D, square array')
90 if damping is None:
91 damping = np.zeros(stiffness.shape)
92 else:
93 damping = np.atleast_2d(np.array(damping))
94 if not damping.shape == stiffness.shape:
95 raise ValueError('Damping and Stiffness matrices must be the same shape')
96 if not damping.ndim == 2 or (damping.shape[0] != damping.shape[1]):
97 raise ValueError('Damping should be a 2D, square array')
98 if transformation is None:
99 transformation = np.eye(mass.shape[0])
100 else:
101 transformation = np.atleast_2d(np.array(transformation))
102 if not transformation.ndim == 2:
103 raise ValueError('transformation must be 2D')
104 if not transformation.shape[-1] == mass.shape[0]:
105 raise ValueError(
106 'transformation must have number of columns equal to the number of rows in the mass matrix')
107 coordinate = np.atleast_1d(coordinate)
108 if not isinstance(coordinate, CoordinateArray):
109 raise ValueError('coordinate must be a CoordinateArray object')
110 if not coordinate.ndim == 1 or coordinate.shape[0] != transformation.shape[0]:
111 raise ValueError(
112 'coordinate must be 1D and have the same size as transformation.shape[0] or mass.shape[0] if no transformation is specified')
113 # Check symmetry
114 if not np.allclose(mass, mass.T) and enforce_symmetry:
115 raise ValueError('mass matrix must be symmetric')
116 if not np.allclose(stiffness, stiffness.T, atol=1e-6*stiffness.max()) and enforce_symmetry:
117 raise ValueError('stiffness matrix must be symmetric')
118 if not np.allclose(damping, damping.T) and enforce_symmetry:
119 raise ValueError('damping matrix must be symmetric')
121 self._coordinate = coordinate
122 self._mass = mass
123 self._stiffness = stiffness
124 self._damping = damping
125 self._transformation = transformation
126 self._enforce_symmetry = enforce_symmetry
128 def __repr__(self):
129 return 'System with {:} DoFs ({:} internal DoFs)'.format(self.ndof_transformed, self.ndof)
131 def spy(self, subplots_kwargs={'figsize': (10, 3)}, spy_kwargs={}):
132 """
133 Plot the structure of the system's matrices
135 Parameters
136 ----------
137 subplots_kwargs : dict, optional
138 Default arguments passed to `matplotlib.pyplot`'s `subplots` function.
139 The default is {'figsize':(10,3)}.
140 spy_kwargs : dict, optional
141 Default arguments passed to `matplotlib.pyplot`'s 'spy' function.
142 The default is {}.
144 Returns
145 -------
146 ax : Axes
147 Axes on which the subplots are defined.
149 """
150 fig, ax = plt.subplots(1, 3, squeeze=True, **subplots_kwargs)
151 ax[0].spy(abs(self.transformation), **spy_kwargs)
152 ax[1].spy(abs(self.mass) + abs(self.stiffness) + abs(self.damping), **spy_kwargs)
153 ax[2].spy(abs(self.transformation.T), **spy_kwargs)
154 ax[0].set_title('Output Transformation')
155 ax[1].set_title('Internal State Matrices')
156 ax[2].set_title('Input Transformation')
157 fig.tight_layout()
158 # Trying to figure out how to scale the subfigures identically...
159 # plt.pause(0.01)
160 # # Now adjust the sizes of the plots
161 # sizes = []
162 # for a in ax:
163 # bbox = a.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
164 # width, height = bbox.width, bbox.height
165 # sizes.append((width, height))
166 # scales = [sizes[0][0],sizes[1][1],sizes[2][1]]
167 # plt.pause(0.01)
168 # scales = [scale/min(scales) for scale in scales]
169 # for a,scale in zip(ax,scales):
170 # position = a.get_position()
171 # center_x = (position.xmax + position.xmin)/2
172 # center_y = (position.ymax + position.ymin)/2
173 # new_left = center_x - (center_x - position.xmin)/2
174 # new_bottom = center_y - (center_y - position.ymin)/2
175 # a.set_position([new_left,new_bottom,position.width/scale,position.height/scale])
176 return ax
178 @property
179 def transformation(self):
180 """Get or set the transformation matrix"""
181 return self._transformation
183 @transformation.setter
184 def transformation(self, value):
185 """Set the transformation matrix"""
186 if not value.ndim == 2:
187 raise ValueError('transformation must be 2D')
188 if not value.shape[-1] == self.mass.shape[0]:
189 raise ValueError(
190 'transformation must have number of columns equal to the number of rows in the mass matrix')
191 self._transformation = value
193 @property
194 def coordinate(self):
195 """Get or set the degrees of freedom in the system"""
196 return self._coordinate
198 @coordinate.setter
199 def coordinate(self, value):
200 """Set the degrees of freedom in the system"""
201 if not isinstance(value, CoordinateArray):
202 raise ValueError('coordinate must be a CoordinateArray object')
203 if not value.ndim == 1 or value.shape[0] != self.mass.shape[0]:
204 raise ValueError('coordinate must be 1D and have the same size as mass.shape[0]')
205 self._coordinate = value
207 @property
208 def mass(self):
209 """Get or set the mass matrix of the system"""
210 return self._mass
212 @mass.setter
213 def mass(self, value):
214 """Set the mass matrix of the system"""
215 if not value.shape == self.stiffness.shape:
216 raise ValueError('Mass and Stiffness matrices must be the same shape')
217 if not value.ndim == 2 or (value.shape[0] != value.shape[1]):
218 raise ValueError('Mass should be a 2D, square array')
219 if not np.allclose(value, value.T) and self._enforce_symmetry:
220 raise ValueError('mass matrix must be symmetric')
221 self._mass = value
223 @property
224 def stiffness(self):
225 """Get or set the stffness matrix of the system"""
226 return self._stiffness
228 @stiffness.setter
229 def stiffness(self, value):
230 """Set the stiffness matrix of the system"""
231 if not value.shape == self.mass.shape:
232 raise ValueError('Mass and Stiffness matrices must be the same shape')
233 if not value.ndim == 2 or (value.shape[0] != value.shape[1]):
234 raise ValueError('Stiffness should be a 2D, square array')
235 if not np.allclose(value, value.T) and self._enforce_symmetry:
236 raise ValueError('stiffness matrix must be symmetric')
237 self._stiffness = value
239 @property
240 def damping(self):
241 """Get or set the damping matrix of the system"""
242 return self._damping
244 @damping.setter
245 def damping(self, value):
246 """Set the damping matrix of the system"""
247 if not value.shape == self.stiffness.shape:
248 raise ValueError('damping and Stiffness matrices must be the same shape')
249 if not value.ndim == 2 or (value.shape[0] != value.shape[1]):
250 raise ValueError('damping should be a 2D, square array')
251 if not np.allclose(value, value.T) and self._enforce_symmetry:
252 raise ValueError('damping matrix must be symmetric')
253 self._damping = value
255 M = mass
256 K = stiffness
257 C = damping
259 @property
260 def ndof(self):
261 """Get the number of internal degrees of freedom of the system"""
262 return self.mass.shape[0]
264 @property
265 def ndof_transformed(self):
266 """Get the number of physical degrees of freedom of the system"""
267 return self.transformation.shape[0]
269 @property
270 def massless_dofs(self):
271 """Gets the massless degrees of freedom in the system"""
272 return np.all(self.mass == 0, axis=0)
274 @property
275 def massive_dofs(self):
276 return ~self.massless_dofs
278 def to_state_space(self, response_coordinates = None, excitation_coordinates = None,
279 output_excitation_signals = True, extra_outputs = None):
280 """
281 Creates state space matrices A, B, C, and D from a SDynPy System object
283 Parameters
284 ----------
285 response_coordinates : dict, optional
286 A dictionary with keys 0, 1, and 2 representing the displacement
287 derivative and CoordinateArray values representing the degrees of
288 freedom desired at that degree of freedom. For example, a
289 dictionary {0:sdpy.coordinate_array([101,102],['X+','X+']),
290 2:sdpy.coordinate_array([203,204],['X+','Y+'])} would give
291 displacement results at 101X+ and 102X+ and acceleration results at
292 203X+ and 204Y+. The value of None can also be supplied to get all
293 degrees of freedom at that derivative. For example, {2:None} would
294 get acceleration results at all degrees of freedom in the System.
295 If not specified, the output will be given for displacement,
296 velocity, and acceleration, for all degrees of freedom. Rows of
297 the output vector will be ordered first by derivative order then
298 by degree of freedom order, so all displacements come before all
299 velocities, and those before all accelerations.
300 excitation_coordinates : CoordinateArray, optional
301 The degrees of freedom to use as inputs in the state and output
302 equations. The default is to provide inputs at all degrees of freedom.
303 output_excitation_signals : bool, optional
304 If True, the inputs to the system will be passed through and output
305 as part of the output matrices C and D. The default is True. These
306 will be in the last rows of C and D, one for each of the excitation
307 degrees of freedom.
308 extra_outputs : dict
309 A dictionary containing data to allow additional outputs to be
310 obtained from the state space formulation. The dictionary must
311 have keys 'C' and 'D'. These must contain NumPy ndarray objects that
312 will be concatenated with the C and D matrices. Both matrices
313 should have number of rows equal to the number of extra outputs
314 requested. The number of columns of C should be equal to the
315 number of states (2x number of massive degrees of freedom + 1x
316 number of massless degrees of freedom). States are ordered
317 displacements of massive
318 degrees of freedom, velocities of massive degrees of freedom,
319 displacements of massless degrees of freedom. Within each group,
320 they will be ordered in the same order as the system matrices (not
321 the transformation matrix). The number of columns of the D matrix
322 should be equivalent to the number of excitation coordinates
323 provided.
325 Returns
326 -------
327 A : ndarray
328 The state matrix. States are ordered displacements of massive
329 degrees of freedom, velocities of massive degrees of freedom,
330 displacements of massless degrees of freedom. Within each group,
331 they will be ordered in the same order as the system matrices (not
332 the transformation matrix).
333 B : ndarray
334 The input matrix. Inputs are in the same order as the
335 excitation_coordinate argument. States are in the order as
336 described for the A return value.
337 C : ndarray
338 The output matrix. States are in the order as described for the A
339 return value. Outputs are ordered as follows: displacement
340 responses, velocity responses, acceleration responses,
341 excitation signals, extra outputs. Within each group,
342 the coordinates are ordered the same way as the provided
343 response_coordinate and excitation_coordinate arguments.
344 D : ndarray
345 The feedforward matrix. Input order is described in the return
346 value for B. Output order is described in the return value for C.
347 """
348 if response_coordinates is None:
349 response_coordinates = {0:self.coordinate, 1:self.coordinate, 2:self.coordinate}
350 phi_response = {}
351 for derivative in [0,1,2]:
352 if derivative in response_coordinates:
353 if isinstance(response_coordinates[derivative], str) and response_coordinates[derivative] == 'state':
354 # If you want the state variables back, the transformation is the
355 # identity matrix.
356 phi_response[derivative] = np.eye(self.transformation.shape[-1])
357 elif response_coordinates[derivative] is None:
358 phi_response[derivative] = self.transformation
359 else:
360 phi_response[derivative] = self.transformation_matrix_at_coordinates(response_coordinates[derivative])
361 else:
362 phi_response[derivative] = self.transformation_matrix_at_coordinates(CoordinateArray(0))
363 if excitation_coordinates is None:
364 phi_input = self.transformation
365 else:
366 phi_input = self.transformation_matrix_at_coordinates(excitation_coordinates)
367 tdofs_input = phi_input.shape[0]
369 massless_dofs = np.all(self.mass == 0, axis=0)
370 ix_11 = np.ix_(~massless_dofs,~massless_dofs)
371 ix_12 = np.ix_(~massless_dofs,massless_dofs)
372 ix_22 = np.ix_(massless_dofs,massless_dofs)
373 ix_21 = np.ix_(massless_dofs,~massless_dofs)
374 n2 = np.sum(massless_dofs)
375 n1 = np.sum(~massless_dofs)
376 M_11 = self.mass[ix_11]
377 K_11 = self.stiffness[ix_11]
378 K_22 = self.stiffness[ix_22]
379 K_12 = self.stiffness[ix_12]
380 K_21 = self.stiffness[ix_21]
381 C_11 = self.damping[ix_11]
382 C_22 = self.damping[ix_22]
383 C_12 = self.damping[ix_12]
384 C_21 = self.damping[ix_21]
385 Z_11 = np.zeros((n1,n1))
386 Z_12 = np.zeros((n1,n2))
387 Z_1i = np.zeros((n1,tdofs_input))
388 Z_2i = np.zeros((n2,tdofs_input))
389 I_1 = np.eye(n1)
390 I_i = np.eye(tdofs_input)
391 Phi_i1 = phi_input[:,~massless_dofs]
392 Phi_i2 = phi_input[:,massless_dofs]
394 tdofs_response = {}
395 Phi_r1 = {}
396 Phi_r2 = {}
397 Z_r1 = {}
398 Z_r2 = {}
399 Z_ri = {}
400 for derivative in [0,1,2]:
401 tdofs_response[derivative] = phi_response[derivative].shape[0]
402 Phi_r1[derivative] = phi_response[derivative][:,~massless_dofs]
403 Phi_r2[derivative] = phi_response[derivative][:,massless_dofs]
404 Z_r1[derivative] = np.zeros((tdofs_response[derivative],n1))
405 Z_r2[derivative] = np.zeros((tdofs_response[derivative],n2))
406 Z_ri[derivative] = np.zeros((tdofs_response[derivative],tdofs_input))
408 if np.any(Phi_r2[2] != 0):
409 warnings.warn('It looks like you are asking for the second order derivative of a degree of freedom that only has a first-order derivative defined. Acceleration of massless degrees of freedom cannot be computed and will be set to zero.')
411 # Degree of freedom ordering = [x1,v1,x2] (Disp 1, Vel 1, Disp 2)
412 # Input ordering [V] (Voltages)
414 # # x1, v1, x2
415 A = np.block([[Z_11,I_1, Z_12],
416 [np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,K_21) - np.linalg.solve(M_11,K_11),
417 -np.linalg.solve(M_11,C_11) + np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,C_21),
418 np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,K_22) - np.linalg.solve(M_11,K_12)],
419 [-np.linalg.solve(C_22,K_21), -np.linalg.solve(C_22,C_21), -np.linalg.solve(C_22,K_22)]])
421 B = np.block([[Z_1i],
422 [np.linalg.solve(M_11,Phi_i1.T)+np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,Phi_i2.T)],
423 [np.linalg.solve(C_22,Phi_i2.T)]])
425 # Contributions from the massive degrees of freedom
426 C1 = np.block([
427 # For displacements, just the states x1
428 [Phi_r1[0], Z_r1[0], Z_r2[0]],
429 # For velocities, just the states v1
430 [Z_r1[1], Phi_r1[1], Z_r2[1]],
431 # For acceleration, we want the row of A corresponding to v1dot
432 [Phi_r1[2] @ np.block([[np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,K_21) - np.linalg.solve(M_11,K_11),
433 -np.linalg.solve(M_11,C_11) + np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,C_21),
434 np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,K_22) - np.linalg.solve(M_11,K_12)]])],
435 # For the inputs, we pass through the inputs
436 [Z_1i.T,Z_1i.T,Z_2i.T]])
438 # Contributions from the massless degrees of freedom
439 C2 = np.block([
440 # For displacements, just the states x2
441 [Z_r1[0], Z_r1[0], Phi_r2[0]],
442 # For the velocities, we want the row of A corresponding to x2dot
443 [Phi_r2[1] @ np.block([[-np.linalg.solve(C_22,K_21), -np.linalg.solve(C_22,C_21), -np.linalg.solve(C_22,K_22)]])],
444 # For the accelerations, we will have zeros
445 [Z_r1[2], Z_r1[2], Z_r2[2]],
446 # For the inputs, we will pass through the inputs
447 [Z_1i.T,Z_1i.T,Z_2i.T]])
449 C = C1+C2
451 D = np.block([
452 # For the displacements, we will have zeros
453 [Z_ri[0]],
454 # For the velocities, we want the row of B corresponding to x2dot
455 [Phi_r2[1] @ np.linalg.solve(C_22,Phi_i2.T)],
456 # For the accelerations, we want the row of B corresponding to v1dot
457 [Phi_r1[2] @ (np.linalg.solve(M_11,Phi_i1.T)+np.linalg.solve(M_11,C_12)@np.linalg.solve(C_22,Phi_i2.T))],
458 # For the inputs, we will pass through the inputs
459 [I_i]
460 ])
462 if not output_excitation_signals:
463 C = C[:-tdofs_input]
464 D = D[:-tdofs_input]
466 if extra_outputs is not None:
467 C = np.concatenate((C,extra_outputs['C']))
468 D = np.concatenate((D,extra_outputs['D']))
470 return A,B,C,D
472 def time_integrate(self, forces, responses=None,
473 initial_state=None,
474 integration_oversample=1, extra_outputs = None):
475 """
476 Integrate a system to produce responses to an excitation
478 Parameters
479 ----------
480 forces : TimeHistoryArray
481 The forces applied to the system, which should be a
482 `TimeHistoryArray`. The time step and excitation degrees of freedom
483 will be taken from the `TimeHistoryArray`
484 responses : dict, optional
485 A dictionary with keys 0, 1, and 2 representing the displacement
486 derivative and CoordinateArray values representing the degrees of
487 freedom desired at that degree of freedom. For example, a
488 dictionary {0:sdpy.coordinate_array([101,102],['X+','X+']),
489 2:sdpy.coordinate_array([203,204],['X+','Y+'])} would give
490 displacement results at 101X+ and 102X+ and acceleration results at
491 203X+ and 204Y+. The value of None can also be supplied to get all
492 degrees of freedom at that derivative. For example, {2:None} would
493 get acceleration results at all degrees of freedom in the System.
494 If not specified, the output will be given for displacement,
495 velocity, and acceleration, for all degrees of freedom. Rows of
496 the output vector will be ordered first by derivative order then
497 by degree of freedom order, so all displacements come before all
498 velocities, and those before all accelerations.
499 initial_state : np.ndarray, optional
500 The initial conditions of the integration. The default is zero
501 displacement and zero velocity.
502 integration_oversample : int
503 The amount of oversampling that will be applied to the force by
504 zero-padding the fft. It is generally better to create the forces
505 such that they are oversampled than to use this argument.
506 extra_outputs : dict
507 A dictionary containing data to allow additional outputs to be
508 obtained from the state space formulation. The dictionary must
509 have keys 'C', 'D', and 'coordinate'.
510 Values for 'C' and 'D' must contain NumPy ndarray objects that
511 will be concatenated with the C and D matrices. See the
512 documentation for `to_state_space` for more information on these
513 keys. Additionally, the values of the 'coordinate' key must
514 contain a CoordinateArray object with the same size as the number
515 of rows of C and D. These will be assigned as the coordinates of
516 the extra signals in the outputted TimeHistoryArray.
518 Returns
519 -------
520 response_array : TimeHistoryArray
521 The responses of the system to the forces applied
523 """
524 from .sdynpy_data import data_array, FunctionTypes
525 dt = forces.abscissa_spacing
526 references = forces.coordinate[..., 0]
527 forces = forces.ordinate
528 references = np.atleast_1d(references)
529 A, B, C, D = self.to_state_space(responses, references, False,
530 extra_outputs=extra_outputs)
531 forces = np.atleast_2d(forces)
532 times = np.arange(forces.shape[-1]) * dt
533 linear_system = StateSpace(A, B, C, D)
534 if integration_oversample != 1:
535 forces, times = resample(forces, len(times) * integration_oversample, times, axis=-1)
536 times_out, time_response, x_out = lsim(linear_system, forces.T, times, initial_state)
537 if time_response.ndim == 1:
538 time_response = time_response[:, np.newaxis]
539 response_coordinates = []
540 if responses is None:
541 responses = {0:self.coordinate, 1:self.coordinate, 2:self.coordinate}
542 for derivative in [0,1,2]:
543 if derivative in responses:
544 if isinstance(responses[derivative], str) and responses[derivative] == 'state':
545 # If you want the state variables back, the transformation is the
546 # identity matrix.
547 response_coordinates.append(coordinate_array(np.arange(self.ndof)+1,0))
548 elif responses[derivative] is None:
549 response_coordinates.append(self.coordinate)
550 else:
551 response_coordinates.append(responses[derivative])
552 if extra_outputs is not None:
553 response_coordinates.append(extra_outputs['coordinate'])
554 response_coordinates = np.concatenate(response_coordinates)
555 response_array = data_array(FunctionTypes.TIME_RESPONSE,
556 times, time_response.T,
557 response_coordinates[:, np.newaxis])
558 if integration_oversample != 1:
559 response_array = response_array.extract_elements(
560 slice(None, None, integration_oversample))
561 return response_array
563 def eigensolution(self, num_modes=None, maximum_frequency=None, complex_modes=False, return_shape=True):
564 """
565 Computes the eigensolution of the system
567 Parameters
568 ----------
569 num_modes : int, optional
570 The number of modes of the system to compute. The default is to
571 compute all the modes.
572 maximum_frequency : float, optional
573 The maximum frequency to which modes will be computed.
574 The default is to compute all the modes.
575 complex_modes : bool, optional
576 Whether or not complex modes are computed. The default is False.
577 return_shape : bool, optional
578 Specifies whether or not to return a `ShapeArray` (True) or a reduced
579 `System` (False). The default is True.
581 Raises
582 ------
583 NotImplementedError
584 Raised if complex modes are specified.
586 Returns
587 -------
588 System or ShapeArray
589 If `return_shape` is True, the a ShapeArray will be returned. If
590 `return_shape` is False, a reduced system will be returned.
592 """
593 if complex_modes is False:
594 if num_modes is not None:
595 num_modes = [0, int(num_modes) - 1]
596 if maximum_frequency is not None:
597 maximum_frequency = (2 * np.pi * maximum_frequency)**2
598 maximum_frequency = [-maximum_frequency, maximum_frequency] # Convert to eigenvalue
599 lam, phi = eigh(self.K, self.M, subset_by_index=num_modes,
600 subset_by_value=maximum_frequency)
601 # Mass normalize the mode shapes
602 lam[lam < 0] = 0
603 freq = np.sqrt(lam) / (2 * np.pi)
604 normalized_mass = np.diag(phi.T @ self.M @ phi)
605 phi /= np.sqrt(normalized_mass)
606 # Ignore divide by zero if the frequency is zero
607 with np.errstate(divide='ignore', invalid='ignore'):
608 damping = np.diag(phi.T @ self.C @ phi) / (2 * (2 * np.pi * freq))
609 damping[np.isnan(damping)] = 0.0
610 damping[np.isinf(damping)] = 0.0
611 # Add in the transformation to get back to physical dofs
612 phi = self.transformation @ phi
613 if return_shape:
614 from .sdynpy_shape import shape_array
615 return shape_array(self.coordinate, phi.T, freq, damping)
616 else:
617 return System(self.coordinate, np.eye(freq.size), np.diag((2 * np.pi * freq)**2), np.diag(2 * (2 * np.pi * freq) * damping), phi)
618 else:
619 if self.ndof > 1000:
620 warnings.warn('The complex mode implementation currently computes all eigenvalues and eigenvectors, which may take a long time for large systems.')
621 # For convenience, assign a zeros matrix
622 Z = np.zeros(self.M.shape)
623 A = np.block([[ Z, self.M],
624 [self.M, self.C]])
625 B = np.block([[-self.M, Z],
626 [ Z, self.K]])
627 lam, E = eig(-B,A)
628 # Sort the eigenvalues such that they are increasing in frequency
629 isort = np.argsort(np.abs(lam))
630 lam = lam[isort]
631 E = E[:,isort]
632 # Eigenvalues will be in complex conjugate pairs. Let's only keep
633 # the ones that are greater than zero.
634 keep = lam.imag > 0
635 lam = lam[keep]
636 E = E[:,keep]
637 if A.shape[0]//2 != lam.size:
638 warnings.warn('The complex mode implementation currently does not do well with rigid body modes (0 Hz frequencies). Compute them from a real-modes solution then transform to complex.')
639 # Cull values we don't want
640 if num_modes is not None:
641 lam = lam[:num_modes]
642 E = E[:,:num_modes]
643 if maximum_frequency is not None:
644 keep = np.abs(lam) <= maximum_frequency*2*np.pi
645 lam = lam[keep]
646 E = E[:,keep]
647 # Mass normalize the mode shapes
648 E = E/np.sqrt(np.einsum('ji,jk,ki->i',E,A,E))
649 # Find repeated eigenvalues where the eigenvectors are not orthogonal
650 # TODO: Might have to do some orthogonalization for repeated
651 # eigenvalues
652 # A_modal = np.einsum('ji,jk,kl->il',E,A,E)
653 # Extract just the displacement partition
654 psi = E[E.shape[0]//2:,:]
655 # Add in a transformation to get to physical dofs
656 psi_t = self.transformation @ psi
657 frequency = np.abs(lam)/(2*np.pi)
658 damping = -np.real(lam)/np.abs(lam)
659 if return_shape:
660 from .sdynpy_shape import shape_array
661 return shape_array(self.coordinate, psi.T, frequency, damping)
662 else:
663 warnings.warn('Complex Modes will in general not diagonalize the system M, C, and K, matrices.')
664 return System(self.coordinate,
665 psi.T@self.M@psi,
666 psi.T@self.K@psi.T,
667 psi.T@self.C@psi.T,
668 psi_t)
670 def transformation_shapes(self, shape_indices=None):
671 from .sdynpy_shape import shape_array
672 if shape_indices is None:
673 shape_indices = slice(None)
674 shape_matrix = self.transformation[:, shape_indices]
675 return shape_array(self.coordinate, shape_matrix.T,
676 frequency=0, damping=0)
678 def remove_transformation(self):
679 return System(coordinate_array(np.arange(self.ndof)+1,0),
680 self.mass.copy(),self.stiffness.copy(),self.damping.copy())
682 def frequency_response(self, frequencies, responses=None, references=None,
683 displacement_derivative=0):
684 """
685 Computes frequency response functions at the specified frequency lines.
687 Parameters
688 ----------
689 frequencies : ndarray
690 A 1D array of frequencies.
691 responses : CoordinateArray, optional
692 A set of coordinates to compute responses. The default is to
693 create responses at all coordinates.
694 references : CoordinateArray, optional
695 A set of coordinates to use as inputs. The default is to use all
696 coordinates as inputs.
697 displacement_derivative : int, optional
698 The number of derivatives to apply to the response. The default is
699 0, which corresponds to displacement. 1 would be the first
700 derivative, velocity, and 2 would be the second derivative,
701 acceleration.
703 Returns
704 -------
705 frf : TransferFunctionArray
706 A TransferFunctionArray containing the frequency response function
707 for the system at the specified input and output degrees of freedom.
709 """
710 H = spfrf.sysmat2frf(frequencies, self.M, self.C, self.K)
711 H = (1j * (2 * np.pi * frequencies[:, np.newaxis, np.newaxis]))**displacement_derivative * H
712 # Apply transformations
713 if responses is None:
714 output_transform = self.transformation
715 output_coordinates = self.coordinate
716 else:
717 responses = np.atleast_1d(responses)
718 output_transform = self.transformation_matrix_at_coordinates(responses)
719 output_coordinates = responses
720 if references is None:
721 input_transform = self.transformation
722 input_coordinates = self.coordinate
723 else:
724 references = np.atleast_1d(references)
725 input_transform = self.transformation_matrix_at_coordinates(references)
726 input_coordinates = references
727 H = output_transform @ H @ input_transform.T
728 # Put it into a transfer function array
729 from .sdynpy_data import data_array, FunctionTypes
730 frf = data_array(FunctionTypes.FREQUENCY_RESPONSE_FUNCTION, frequencies,
731 np.moveaxis(H, 0, -1),
732 outer_product(output_coordinates, input_coordinates))
733 return frf
735 def assign_modal_damping(self, damping_ratios):
736 """
737 Assigns a damping matrix to the system that results in equivalent
738 modal damping
740 Parameters
741 ----------
742 damping_ratios : ndarray
743 An array of damping values to assign to the system
745 Returns
746 -------
747 None.
749 """
750 damping_ratios = np.array(damping_ratios)
751 if damping_ratios.ndim == 1:
752 shapes = self.eigensolution(num_modes=damping_ratios.size)
753 else:
754 shapes = self.eigensolution()
755 shapes.damping = damping_ratios
756 # Compute the damping matrix
757 modal_system = shapes.system()
758 shape_pinv = np.linalg.pinv(modal_system.transformation.T)
759 full_damping_matrix = shape_pinv@modal_system.damping@shape_pinv.T
760 self.damping[:] = full_damping_matrix
762 def save(self, filename):
763 """
764 Saves the system to a file
766 Parameters
767 ----------
768 filename : str
769 Name of the file in which the system will be saved.
771 Returns
772 -------
773 None.
775 """
776 np.savez(filename, mass=self.mass, stiffness=self.stiffness, damping=self.damping,
777 transformation=self.transformation, coordinate=self.coordinate.view(np.ndarray),
778 enforce_symmetry = self._enforce_symmetry)
780 @classmethod
781 def load(cls, filename):
782 """
783 Load a system from a file
785 Parameters
786 ----------
787 filename : str
788 Name of the file from which the system will be loaded.
790 Returns
791 -------
792 System
793 A system consisting of the mass, stiffness, damping, and transformation
794 in the file
796 """
797 if isinstance(filename, Path):
798 filename = str(filename)
799 data = np.load(filename)
800 return cls(data['coordinate'].view(CoordinateArray), data['mass'],
801 data['stiffness'], data['damping'], data['transformation'],
802 True if 'enforce_symmetry' not in data else data['enforce_symmetry'])
804 def __neg__(self):
805 new_system = copy.deepcopy(self)
806 new_system.mass *= -1
807 new_system.stiffness *= -1
808 new_system.damping *= -1
809 return new_system
811 def __mul__(self, value):
812 new_system = copy.deepcopy(self)
813 new_system.mass *= value
814 new_system.stiffness *= value
815 new_system.damping *= value
816 return new_system
818 def __rmul__(self, value):
819 new_system = copy.deepcopy(self)
820 new_system.mass *= value
821 new_system.stiffness *= value
822 new_system.damping *= value
823 return new_system
825 @classmethod
826 def concatenate(cls, systems, coordinate_node_offset=0):
827 """
828 Combine multiple systems together
830 Parameters
831 ----------
832 systems : iterable of System objects
833 Iterable of Systems that will be concatenated. Matrices will be
834 assembled in block diagonal format
835 coordinate_node_offset : int, optional
836 Offset applied to the coordinates so the nodes do not overlap.
837 The default is 0.
839 Returns
840 -------
841 System
842 A system consisting of the combintation of the provided systems.
844 """
845 coordinates = [system.coordinate.copy() for system in systems]
846 if coordinate_node_offset != 0:
847 for i in range(len(coordinates)):
848 coordinates[i].node += coordinate_node_offset * (i + 1)
849 all_coordinates = np.concatenate(coordinates)
850 return cls(all_coordinates,
851 block_diag(*[system.mass for system in systems]),
852 block_diag(*[system.stiffness for system in systems]),
853 block_diag(*[system.damping for system in systems]),
854 block_diag(*[system.transformation for system in systems]),
855 enforce_symmetry = True if all([s._enforce_symmetry for s in systems]) else False)
857 @classmethod
858 def substructure_by_position(cls, systems, geometries, distance_threshold=1e-8, rcond=None):
859 """
860 Applies constraints to systems by constraining colocated nodes together
862 Parameters
863 ----------
864 systems : iterable of System objects
865 A set of systems that will be combined and constrained together
866 geometries : iterable of Geometry objects
867 A set of geometries that will be combined together
868 distance_threshold : float, optional
869 The distance between nodes that are considered colocated.
870 The default is 1e-8.
871 rcond : float, optional
872 Condition number to use in the nullspace calculation on
873 the constraint matrix. The default is None.
875 Returns
876 -------
877 combined_system : System
878 System consisting of constraining the input systems together.
879 combined_geometry : Geometry
880 Combined geometry of the new system
882 """
883 from .sdynpy_geometry import Geometry
884 combined_geometry, node_offset = Geometry.overlay_geometries(
885 geometries, return_node_id_offset=True)
886 combined_system = cls.concatenate(systems, node_offset)
887 global_coords = combined_geometry.global_node_coordinate()
888 node_distances = np.linalg.norm(
889 global_coords[:, np.newaxis, :] - global_coords[np.newaxis, :, :], axis=-1)
890 # Find locations where the value is less than the tolerances, except for on the centerline
891 node_pairs = [[combined_geometry.node.id[index] for index in pair]
892 for pair in zip(*np.where(node_distances < distance_threshold)) if pair[0] < pair[1]]
893 # Find matching DoF pairs
894 constraint_matrix = []
895 for node1, node2 in node_pairs:
896 # Find the dofs associated with each node
897 system_1_dof_indices = np.where(combined_system.coordinate.node == node1)[0]
898 system_1_dofs = combined_system.coordinate[system_1_dof_indices]
899 system_1_transformation = combined_system.transformation[system_1_dof_indices]
900 global_deflections_1 = combined_geometry.global_deflection(system_1_dofs)
901 system_2_dof_indices = np.where(combined_system.coordinate.node == node2)[0]
902 system_2_dofs = combined_system.coordinate[system_2_dof_indices]
903 system_2_transformation = combined_system.transformation[system_2_dof_indices]
904 global_deflections_2 = combined_geometry.global_deflection(system_2_dofs)
905 # Split between translations and rotations
906 translation_map_1 = np.where((abs(system_1_dofs.direction) <= 3)
907 & (abs(system_1_dofs.direction) > 0))
908 rotation_map_1 = np.where(abs(system_1_dofs.direction) > 3)
909 translation_map_2 = np.where((abs(system_2_dofs.direction) <= 3)
910 & (abs(system_2_dofs.direction) > 0))
911 rotation_map_2 = np.where(abs(system_2_dofs.direction) > 3)
912 neutral_map_1 = np.where(abs(system_1_dofs.direction) == 0)
913 neutral_map_2 = np.where(abs(system_2_dofs.direction) == 0)
914 # Do translations and rotations separately
915 for map_1, map_2 in [[translation_map_1, translation_map_2],
916 [rotation_map_1, rotation_map_2],
917 [neutral_map_1, neutral_map_2]]:
918 deflections_1 = global_deflections_1[map_1].T
919 deflections_2 = global_deflections_2[map_2].T
920 transform_1 = system_1_transformation[map_1]
921 transform_2 = system_2_transformation[map_2]
922 full_constraint = deflections_1 @ transform_1 - deflections_2 @ transform_2
923 constraint_matrix.append(full_constraint)
924 constraint_matrix = np.concatenate(constraint_matrix, axis=0)
925 return combined_system.constrain(constraint_matrix, rcond), combined_geometry
927 def constrain(self, constraint_matrix, rcond=None):
928 """
929 Apply a constraint matrix to the system
931 Parameters
932 ----------
933 constraint_matrix : np.ndarray
934 A matrix of constraints to apply to the structure (B matrix in
935 substructuring literature)
936 rcond : float, optional
937 Condition tolerance for computing the nullspace. The default is None.
939 Returns
940 -------
941 System
942 Constrained system.
944 """
945 substructuring_transform_matrix = null_space(constraint_matrix, rcond)
946 new_mass = substructuring_transform_matrix.T @ self.mass @ substructuring_transform_matrix
947 new_stiffness = substructuring_transform_matrix.T @ self.stiffness @ substructuring_transform_matrix
948 new_damping = substructuring_transform_matrix.T @ self.damping @ substructuring_transform_matrix
949 new_transform = self.transformation @ substructuring_transform_matrix
950 return System(self.coordinate, new_mass, new_stiffness, new_damping, new_transform,
951 enforce_symmetry=self._enforce_symmetry)
953 def transformation_matrix_at_coordinates(self, coordinates):
954 """
955 Return the transformation matrix at the specified coordinates
957 Parameters
958 ----------
959 coordinates : CoordinateArray
960 coordinates at which the transformation matrix will be computed.
962 Raises
963 ------
964 ValueError
965 Raised if duplicate coordinates are requested, or if coordinates that
966 do not exist in the system are requested.
968 Returns
969 -------
970 return_value : np.ndarray
971 Portion of the transformation matrix corresponding to the
972 coordinates input to the function.
974 """
975 consistent_arrays, shape_indices, request_indices = np.intersect1d(
976 abs(self.coordinate), abs(coordinates), assume_unique=False, return_indices=True)
977 # Make sure that all of the keys are actually in the consistent array matrix
978 if consistent_arrays.size != coordinates.size:
979 extra_keys = np.setdiff1d(abs(coordinates), abs(self.coordinate))
980 if extra_keys.size == 0:
981 raise ValueError(
982 'Duplicate coordinate values requested. Please ensure coordinate indices are unique.')
983 raise ValueError(
984 'Not all indices in requested coordinate array exist in the system\n{:}'.format(str(extra_keys)))
985 # Handle sign flipping
986 multiplications = coordinates.flatten()[request_indices].sign(
987 ) * self.coordinate[shape_indices].sign()
988 return_value = self.transformation[shape_indices] * multiplications[:, np.newaxis]
989 # Invert the indices to return the dofs in the correct order as specified in keys
990 inverse_indices = np.zeros(request_indices.shape, dtype=int)
991 inverse_indices[request_indices] = np.arange(len(request_indices))
992 return_value = return_value[inverse_indices]
993 return return_value
995 def substructure_by_coordinate(self, dof_pairs, rcond=None,
996 return_constrained_system=True):
997 """
998 Constrain the system by connecting the specified degree of freedom pairs
1000 Parameters
1001 ----------
1002 dof_pairs : iterable of CoordinateArray
1003 Pairs of coordinates to be connected. None can be passed instead of
1004 a second degree of freedom to constrain to ground
1005 rcond : float, optional
1006 Condition threshold to use for the nullspace calculation on the
1007 constraint matrix. The default is None.
1008 return_constrained_system : bool, optional
1009 If true, apply the constraint matrix and return the constrained
1010 system, otherwise simply return the constraint matrix. The default
1011 is True.
1013 Returns
1014 -------
1015 np.ndarray or System
1016 Returns a System object with the constraints applied if
1017 `return_constrained_system` is True, otherwise just return the
1018 constraint matrix.
1020 """
1021 constraint_matrix = []
1022 for constraint_dof_0, constraint_dof_1 in dof_pairs:
1023 constraint = self.transformation_matrix_at_coordinates(constraint_dof_0)
1024 if constraint_dof_1 is not None:
1025 constraint -= self.transformation_matrix_at_coordinates(constraint_dof_1)
1026 constraint_matrix.append(constraint)
1027 constraint_matrix = np.concatenate(constraint_matrix, axis=0)
1028 if return_constrained_system:
1029 return self.constrain(constraint_matrix, rcond)
1030 else:
1031 return constraint_matrix
1033 def substructure_by_shape(self, constraint_shapes, connection_dofs_0,
1034 connection_dofs_1=None, rcond=None,
1035 return_constrained_system=True):
1036 """
1037 Constrain the system using a set of shapes in a least-squares sense.
1039 Parameters
1040 ----------
1041 constraint_shapes : ShapeArray
1042 An array of shapes to use as the basis for the constraints
1043 connection_dofs_0 : CoordinateArray
1044 Array of coordinates to use in the constraints
1045 connection_dofs_1 : CoordinateArray, optional
1046 Array of coordinates to constrain to the coordinates in
1047 `connection_dofs_0`. If not specified, the `connection_dofs_0`
1048 degrees of freedom will be constrained to ground.
1049 rcond : float, optional
1050 Condition threshold on the nullspace calculation. The default is None.
1051 return_constrained_system : bool, optional
1052 If true, apply the constraint matrix and return the constrained
1053 system, otherwise simply return the constraint matrix. The default
1054 is True.
1056 Returns
1057 -------
1058 np.ndarray or System
1059 Returns a System object with the constraints applied if
1060 `return_constrained_system` is True, otherwise just return the
1061 constraint matrix.
1063 """
1064 shape_matrix_0 = constraint_shapes[connection_dofs_0].T
1065 transform_matrix_0 = self.transformation_matrix_at_coordinates(connection_dofs_0)
1066 constraint_matrix = np.linalg.lstsq(shape_matrix_0, transform_matrix_0)[0]
1067 if connection_dofs_1 is not None:
1068 shape_matrix_1 = constraint_shapes[connection_dofs_1].T
1069 transform_matrix_1 = self.transformation_matrix_at_coordinates(connection_dofs_1)
1070 constraint_matrix -= np.linalg.lstsq(shape_matrix_1, transform_matrix_1)[0]
1071 if return_constrained_system:
1072 return self.constrain(constraint_matrix, rcond)
1073 else:
1074 return constraint_matrix
1076 def substructure_shakers(self,shakers,coordinates):
1077 """
1078 Uses substructuring to attach 4DoF Shaker Models to the Substructure
1080 Parameters
1081 ----------
1082 shakers : array of Shaker4DoF
1083 4-DoF electromechanical models of a shaker that are added to the
1084 System
1085 coordinates : CoordinateArray
1086 The degrees of freedom at which the shakers should be attached.
1088 Returns
1089 -------
1090 System
1091 A system with shakers attached to it. The shaker degrees of freedom
1092 will be assigned based on the next power of 10 bigger than the
1093 current maximum node number. Degrees of freedom end in 1 for
1094 armature force/displacement, 2 for body force/displacement, 3 for
1095 force gauge force/displacement, and 4 for drive voltage/current
1096 """
1097 constraint_pairs = []
1098 shaker_systems = []
1100 max_node = self.coordinate.node.max()
1101 shaker_node_offset = int(10**np.ceil(np.log10(max_node)))
1102 if shaker_node_offset == max_node:
1103 shaker_node_offset *= 10
1104 if shaker_node_offset < 100:
1105 shaker_node_offset = 100
1107 for i,(direction,shaker) in enumerate(zip(coordinates,shakers)):
1108 M,C,K = shaker.MCK()
1109 coords = coordinate_array(shaker_node_offset + i*10 + np.arange(4)+1,[1,1,1,0])
1110 shaker_systems.append(System(coords,M,K,C,enforce_symmetry = False))
1111 constraint_dof = coordinate_array(shaker_node_offset + i*10 + 3,1)
1112 constraint_pairs.append([direction,constraint_dof])
1114 constraint_pairs = np.array(constraint_pairs).view(CoordinateArray)
1116 combined_system = System.concatenate([self]+shaker_systems)
1118 constrained_system = combined_system.substructure_by_coordinate(constraint_pairs)
1120 return constrained_system
1122 def copy(self):
1123 """
1124 Returns a copy of the system object
1126 Returns
1127 -------
1128 System
1129 A copy of the system object
1131 """
1132 return System(self.coordinate.copy(), self.mass.copy(), self.stiffness.copy(),
1133 self.damping.copy(), self.transformation.copy(),
1134 self._enforce_symmetry)
1136 def set_proportional_damping(self, mass_fraction, stiffness_fraction):
1137 """
1138 Sets the damping matrix to a proportion of the mass and stiffness matrices.
1140 The damping matrix will be set to `mass_fraction*self.mass +
1141 stiffness_fraction*self.stiffness`
1143 Parameters
1144 ----------
1145 mass_fraction : float
1146 Fraction of the mass matrix
1147 stiffness_fraction : TYPE
1148 Fraction of the stiffness matrix
1150 Returns
1151 -------
1152 None.
1154 """
1155 self.damping = self.mass * mass_fraction + self.stiffness * stiffness_fraction
1157 @classmethod
1158 def beam(cls, length, width, height, num_nodes, E=None, rho=None, nu=None, material=None):
1159 """
1160 Create a beam mass and stiffness matrix
1162 Parameters
1163 ----------
1164 length : float
1165 Lenghth of the beam
1166 width : float
1167 Width of the beam
1168 height : float
1169 Height of the beam
1170 num_nodes : int
1171 Number of nodes in the beam.
1172 E : float, optional
1173 Young's modulus of the beam. If not specified, a `material` must be
1174 specified instead
1175 rho : float, optional
1176 Density of the beam. If not specified, a `material` must be
1177 specified instead
1178 nu : float, optional
1179 Poisson's ratio of the beam. If not specified, a `material` must be
1180 specified instead
1181 material : str, optional
1182 A specific material can be specified instead of `E`, `rho`, and `nu`.
1183 Should be a string 'steel' or 'aluminum'. If not specified, then
1184 options `E`, `rho`, and `nu` must be specified instead.
1186 Raises
1187 ------
1188 ValueError
1189 If improper materials are defined.
1191 Returns
1192 -------
1193 system : System
1194 A system object consisting of the beam mass and stiffness matrices.
1195 geometry : Geometry
1196 A Geometry consisting of the beam geometry.
1198 """
1199 from .sdynpy_geometry import Geometry, node_array, traceline_array, coordinate_system_array
1200 node_positions = np.array((np.linspace(0, length, num_nodes),
1201 np.zeros(num_nodes),
1202 np.zeros(num_nodes))).T
1203 node_connectivity = np.array((np.arange(num_nodes - 1), np.arange(1, num_nodes))).T
1204 bend_direction_1 = np.array((np.zeros(num_nodes - 1),
1205 np.zeros(num_nodes - 1),
1206 np.ones(num_nodes - 1))).T
1207 if material is None:
1208 if E is None or rho is None or nu is None:
1209 raise ValueError('Must specify material or E, nu, and rho')
1210 elif material.lower() == 'steel':
1211 E = 200e9 # [N/m^2],
1212 nu = 0.25 # [-],
1213 rho = 7850 # [kg/m^3]
1214 elif material.lower() == 'aluminum':
1215 E = 69e9 # [N/m^2],
1216 nu = 0.33 # [-],
1217 rho = 2830 # [kg/m^3]
1218 else:
1219 raise ValueError('Unknown Material {:}'.format(material))
1220 mat_props = rect_beam_props(E, rho, nu, width, height, num_nodes - 1)
1221 K, M = beamkm(node_positions, node_connectivity, bend_direction_1, **mat_props)
1222 coordinates = from_nodelist(np.arange(num_nodes) + 1, directions=[1, 2, 3, 4, 5, 6])
1223 system = cls(coordinates, M, K)
1224 nodelist = node_array(np.arange(num_nodes) + 1, node_positions)
1225 tracelines = traceline_array(connectivity=np.arange(num_nodes) + 1)
1226 coordinate_systems = coordinate_system_array()
1227 geometry = Geometry(nodelist, coordinate_systems, tracelines)
1228 return system, geometry
1230 def get_indices_by_coordinate(self, coordinates, ignore_sign=False):
1231 """
1232 Gets the indices in the transformation matrix corresponding coordinates
1234 Parameters
1235 ----------
1236 coordinates : CoordinateArray
1237 Coordinates to extract transformation indices
1238 ignore_sign : bool, optional
1239 Specify whether or not to ignore signs on the coordinates. If True,
1240 then '101X+' would match '101X+' or '101X-'. The default is False.
1242 Raises
1243 ------
1244 ValueError
1245 Raised if duplicate coordinates or coordinates not in the system
1246 are requested
1248 Returns
1249 -------
1250 np.ndarray
1251 Array of indices.
1253 """
1254 if ignore_sign:
1255 consistent_arrays, shape_indices, request_indices = np.intersect1d(
1256 abs(self.coordinate), abs(coordinates), assume_unique=False, return_indices=True)
1257 else:
1258 consistent_arrays, shape_indices, request_indices = np.intersect1d(
1259 self.coordinate, coordinates, assume_unique=False, return_indices=True)
1260 # Make sure that all of the keys are actually in the consistent array matrix
1261 if consistent_arrays.size != coordinates.size:
1262 extra_keys = np.setdiff1d(abs(coordinates), abs(self.coordinate))
1263 if extra_keys.size == 0:
1264 raise ValueError(
1265 'Duplicate coordinate values requested. Please ensure coordinate indices are unique.')
1266 raise ValueError(
1267 'Not all indices in requested coordinate array exist in the system\n{:}'.format(str(extra_keys)))
1268 # Handle sign flipping
1269 return_value = shape_indices
1270 # Invert the indices to return the dofs in the correct order as specified in keys
1271 inverse_indices = np.zeros(request_indices.shape, dtype=int)
1272 inverse_indices[request_indices] = np.arange(len(request_indices))
1273 return return_value[inverse_indices]
1275 def reduce(self, reduction_transformation):
1276 """
1277 Apply the specified reduction to the model
1279 Parameters
1280 ----------
1281 reduction_transformation : np.ndarray
1282 Matrix to use in the reduction
1284 Returns
1285 -------
1286 System
1287 Reduced system.
1289 """
1290 mass = reduction_transformation.T @ self.mass @ reduction_transformation
1291 stiffness = reduction_transformation.T @ self.stiffness @ reduction_transformation
1292 damping = reduction_transformation.T @ self.damping @ reduction_transformation
1293 transformation = self.transformation @ reduction_transformation
1294 # Force symmetry
1295 mass = (mass + mass.T) / 2
1296 stiffness = (stiffness + stiffness.T) / 2
1297 damping = (damping + damping.T) / 2
1298 return System(self.coordinate, mass, stiffness, damping, transformation, self._enforce_symmetry)
1300 def reduce_guyan(self, coordinates):
1301 """
1302 Perform Guyan reduction on the system
1304 Parameters
1305 ----------
1306 coordinates : CoordinateArray
1307 A list of coordinates to keep in the reduced system.
1309 Raises
1310 ------
1311 ValueError
1312 Raised the transformation matrix is not identity matrix.
1314 Returns
1315 -------
1316 System
1317 Reduced system.
1319 """
1320 if isinstance(coordinates, CoordinateArray):
1321 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1322 raise ValueError(
1323 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1324 keep_dofs = self.get_indices_by_coordinate(coordinates)
1325 else:
1326 keep_dofs = np.array(coordinates)
1327 discard_dofs = np.array([i for i in range(self.ndof) if i not in keep_dofs])
1328 I_a = np.eye(keep_dofs.size)
1329 K_dd = self.stiffness[discard_dofs[:, np.newaxis],
1330 discard_dofs]
1331 K_da = self.stiffness[discard_dofs[:, np.newaxis],
1332 keep_dofs]
1333 T_guyan = np.concatenate((I_a, -np.linalg.solve(K_dd, K_da)), axis=0)
1334 T_guyan[np.concatenate((keep_dofs, discard_dofs)), :] = T_guyan.copy()
1335 return self.reduce(T_guyan)
1337 def reduce_dynamic(self, coordinates, frequency):
1338 """
1339 Perform Dynamic condensation
1341 Parameters
1342 ----------
1343 coordinates : CoordinateArray
1344 A list of coordinates to keep in the reduced system.
1345 frequency : float
1346 The frequency to preserve in the dynamic reduction.
1348 Raises
1349 ------
1350 ValueError
1351 Raised if the transformation is not identity matrix.
1353 Returns
1354 -------
1355 System
1356 Reduced system.
1358 """
1359 if isinstance(coordinates, CoordinateArray):
1360 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1361 raise ValueError(
1362 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1363 keep_dofs = self.get_indices_by_coordinate(coordinates)
1364 else:
1365 keep_dofs = np.array(coordinates)
1366 discard_dofs = np.array([i for i in range(self.ndof) if i not in keep_dofs])
1367 I_a = np.eye(keep_dofs.size)
1368 D = self.stiffness - (2 * np.pi * frequency)**2 * self.mass
1369 D_dd = D[discard_dofs[:, np.newaxis],
1370 discard_dofs]
1371 D_da = D[discard_dofs[:, np.newaxis],
1372 keep_dofs]
1373 T_dynamic = np.concatenate((I_a, -np.linalg.solve(D_dd, D_da)), axis=0)
1374 T_dynamic[np.concatenate((keep_dofs, discard_dofs)), :] = T_dynamic.copy()
1375 return self.reduce(T_dynamic)
1377 def reduce_craig_bampton(self, connection_degrees_of_freedom: CoordinateArray,
1378 num_fixed_base_modes: int,
1379 return_shape_matrix: bool = False):
1380 """
1381 Computes a craig-bampton substructure model for the system
1383 Parameters
1384 ----------
1385 connection_degrees_of_freedom : CoordinateArray
1386 Degrees of freedom to keep at the interface.
1387 num_fixed_base_modes : int
1388 Number of fixed-base modes to use in the reduction
1389 return_shape_matrix : bool, optional
1390 If true, return a set of shapes that represents the transformation
1391 in addition to the reduced system. The default is False.
1393 Raises
1394 ------
1395 ValueError
1396 Raised if coordinate arrays are specified when there is already a
1397 transformation.
1399 Returns
1400 -------
1401 System
1402 Reduced system in craig-bampton form
1403 ShapeArray
1404 Shapes representing the craig-bampton transformation
1406 """
1407 # Construct craig bampton transformation
1408 if isinstance(connection_degrees_of_freedom, CoordinateArray):
1409 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1410 raise ValueError(
1411 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1412 connection_indices = self.get_indices_by_coordinate(connection_degrees_of_freedom)
1413 else:
1414 connection_indices = np.array(connection_degrees_of_freedom)
1415 other_indices = np.array([i for i in range(self.ndof) if i not in connection_indices])
1417 # Extract portions of the mass and stiffness matrices
1418 K_ii = self.K[other_indices[:, np.newaxis], other_indices]
1419 M_ii = self.M[other_indices[:, np.newaxis], other_indices]
1420 K_ib = self.K[other_indices[:, np.newaxis], connection_indices]
1422 # Compute fixed interface modes
1423 lam, Phi_ii = eigh(K_ii, M_ii, subset_by_index=[0, int(num_fixed_base_modes) - 1])
1424 # Normalize the mode shapes
1425 lam[lam < 0] = 0
1426 normalized_mass = np.diag(Phi_ii.T @ M_ii @ Phi_ii)
1427 Phi_ii /= np.sqrt(normalized_mass)
1428 Z_bi = np.zeros((connection_indices.size, num_fixed_base_modes))
1429 # Compute constraint modes
1430 Psi_ib = -np.linalg.solve(K_ii, K_ib)
1431 I_bb = np.eye(connection_indices.size)
1432 T_cb = np.block([[Phi_ii, Psi_ib],
1433 [Z_bi, I_bb]])
1434 T_cb[np.concatenate((other_indices, connection_indices)), :] = T_cb.copy()
1435 if return_shape_matrix:
1436 from .sdynpy_shape import shape_array
1437 freq = np.sqrt(lam) / (2 * np.pi)
1438 all_freqs = np.concatenate((freq, np.zeros(connection_indices.size)))
1439 shapes = shape_array(self.coordinate, T_cb.T, all_freqs, comment1=['Fixed Base Mode {:}'.format(
1440 i + 1) for i in range(num_fixed_base_modes)] + ['Constraint Mode {:}'.format(str(dof)) for dof in connection_degrees_of_freedom])
1441 return self.reduce(T_cb), shapes
1442 else:
1443 return self.reduce(T_cb)
1445 @classmethod
1446 def from_exodus_superelement(cls, superelement_nc4, transformation_exodus_file=None,
1447 x_disp='DispX', y_disp='DispY', z_disp='DispZ',
1448 x_rot=None, y_rot=None, z_rot=None,
1449 reduce_to_external_surfaces=False):
1450 """
1451 Creates a system from a superelement from Sierra/SD
1453 Parameters
1454 ----------
1455 superelement_nc4 : netCDF4.Dataset or string
1456 Dataset from which the superelement data will be loaded
1457 transformation_exodus_file : Exodus, ExodusInMemory, or str, optional
1458 Exodus data containing the transformation between the reduced
1459 superelement state and the physical space. If not specified, no
1460 transformation will be created.
1461 x_disp : str, optional
1462 Variable name to read for x-displacements in the transformation
1463 Exodus file. The default is 'DispX'.
1464 y_disp : str, optional
1465 Variable name to read for y-displacements in the transformation
1466 Exodus file. The default is 'DispY'.
1467 z_disp : str, optional
1468 Variable name to read for z-displacements in the transformation
1469 Exodus file. The default is 'DispZ'.
1470 x_rot : str, optional
1471 Variable name to read for x-rotations in the transformation
1472 Exodus file. The default is to not read rotations.
1473 y_rot : str, optional
1474 Variable name to read for y-rotations in the transformation
1475 Exodus file. The default is to not read rotations.
1476 z_rot : str, optional
1477 Variable name to read for z-rotations in the transformation
1478 Exodus file. The default is to not read rotations.
1479 reduce_to_external_surfaces : bool, optional
1480 If True, exodus results will be reduced to external surfaces
1482 Raises
1483 ------
1484 ValueError
1485 raised if bad data types are passed to the arguments.
1487 Returns
1488 -------
1489 system : System
1490 System containing the superelement representation
1491 geometry : Geometry
1492 Geometry that can be used to plot the system
1493 boundary_dofs : CoordinateArray
1494 Degrees of freedom that can be used to constrain the test article.
1496 """
1497 from ..fem.sdynpy_exodus import Exodus, ExodusInMemory, reduce_exodus_to_surfaces, read_sierra_matlab_matrix_file, read_sierra_matlab_map_file
1498 from .sdynpy_geometry import node_array, coordinate_system_array, Geometry
1499 if isinstance(superelement_nc4, str):
1500 ds = nc4.Dataset(superelement_nc4)
1501 elif isinstance(superelement_nc4, nc4.Dataset):
1502 ds = superelement_nc4
1503 else:
1504 raise ValueError('superelement_nc4 must be a string or a netCDF4 Dataset')
1505 cbmap = ds['cbmap'][:].data.copy()
1506 Kr = ds['Kr'][:].data.copy()
1507 Mr = ds['Mr'][:].data.copy()
1508 Cr = ds['Cr'][:].data.copy()
1509 num_constraint_modes = ds.dimensions['NumConstraints'].size
1510 num_fixed_base_modes = ds.dimensions['NumEig'].size
1511 boundary_dofs = coordinate_array(*cbmap[cbmap[:, 0] > 0].T)
1512 if transformation_exodus_file is None:
1513 transformation = None
1514 coordinate_nodes = cbmap[:, 0]
1515 coordinate_dirs = cbmap[:, 1]
1516 coordinate_nodes[coordinate_nodes == 0] = np.arange(num_fixed_base_modes) + 1
1517 coordinates = coordinate_array(coordinate_nodes, coordinate_dirs)
1518 cs_array = coordinate_system_array()
1519 n_array = node_array(ds['node_num_map'][:], np.array(
1520 [ds['coord{:}'.format(d)][:] for d in 'xyz']).T)
1521 geometry = Geometry(node=n_array, coordinate_system=cs_array)
1522 else:
1523 if isinstance(transformation_exodus_file, str):
1524 exo = Exodus(transformation_exodus_file)
1525 elif isinstance(transformation_exodus_file, Exodus):
1526 exo = transformation_exodus_file
1527 elif isinstance(transformation_exodus_file, ExodusInMemory):
1528 exo = transformation_exodus_file
1529 else:
1530 raise ValueError('transformation_exodus_file must be a string or a sdpy.Exodus')
1531 if reduce_to_external_surfaces:
1532 exo = reduce_exodus_to_surfaces(exo, variables_to_transform=[var for var in [x_disp, y_disp, z_disp, x_rot, y_rot, z_rot] if var is not None])
1533 from .sdynpy_shape import ShapeArray
1534 shapes = ShapeArray.from_exodus(exo, x_disp, y_disp, z_disp, x_rot, y_rot, z_rot)
1535 transformation = shapes.shape_matrix.T
1536 coordinates = shapes[0].coordinate
1537 geometry = Geometry.from_exodus(exo)
1538 system = cls(coordinates, Mr, Kr, Cr, transformation)
1539 return system, geometry, boundary_dofs
1541 @classmethod
1542 def from_sierra_sd_mfile_output(cls, maa_file, kaa_file, gid_file, aset_map_file):
1543 """
1544 Generates a System object from the MFILE, MAA, and KAA outputs from
1545 Sierra/SD
1547 Parameters
1548 ----------
1549 maa_file : str
1550 Path to the file containing the MAA matrix.
1551 kaa_file : str
1552 Path to the file containing the KAA matrix.
1553 gid_file : str
1554 Path to the file contiaining the global node id map
1555 aset_map_file : str
1556 Path to the A-set Map.
1558 Returns
1559 -------
1560 system : System
1561 A system object with the specified matrices and degrees of freedom.
1563 """
1564 M = read_sierra_matlab_matrix_file(maa_file)
1565 K = read_sierra_matlab_matrix_file(kaa_file)
1566 gid = read_sierra_matlab_map_file(gid_file)
1567 asetmap = read_sierra_matlab_map_file(aset_map_file).reshape(-1,9)
1568 valid_dofs = asetmap != -1
1569 # Check if any of the acoustic, temperature, or voltage are defined
1570 if np.any(valid_dofs[:,-3:]):
1571 warnings.warn('Sierra system contains acoustic, temperature, or voltage degrees of freedom. These will be ignored by SDynPy.')
1572 valid_dofs = valid_dofs[:,:-3]
1573 asetmap = asetmap[:,:-3]
1574 coordinates = coordinate_array(gid[:,np.newaxis],[1,2,3,4,5,6])[valid_dofs]
1575 dof_ordering = asetmap[valid_dofs]-1
1576 M = M[dof_ordering[:,np.newaxis],dof_ordering]
1577 K = K[dof_ordering[:,np.newaxis],dof_ordering]
1578 system = cls(coordinates,M,K)
1579 return system
1581 def simulate_test(
1582 self, # The system itself
1583 bandwidth,
1584 frame_length,
1585 num_averages,
1586 excitation,
1587 references,
1588 responses=None, # All Responses
1589 excitation_level=1.0,
1590 excitation_noise_level=0.0,
1591 response_noise_level=0.0,
1592 steady_state_time=0.0,
1593 excitation_min_frequency=None,
1594 excitation_max_frequency=None,
1595 signal_fraction=0.5,
1596 extra_time_between_frames=0.0,
1597 integration_oversample=10,
1598 antialias_filter_cutoff_factor=3,
1599 antialias_filter_order=4,
1600 multihammer_impact_spacing_factor = 4,
1601 extra_outputs = None,
1602 **generator_kwargs
1603 ):
1604 available_excitations = ['pseudorandom', 'random',
1605 'burst random', 'chirp', 'hammer', 'multi-hammer', 'sine']
1606 if not excitation.lower() in available_excitations:
1607 raise ValueError('Excitation must be one of {:}'.format(available_excitations))
1608 # Create the input signal
1609 references = np.atleast_1d(references)
1610 num_signals = references.size
1611 sample_rate = bandwidth * 2 * integration_oversample
1612 dt = 1 / sample_rate
1613 frame_time = dt * frame_length * integration_oversample
1614 # Create the signals
1615 if excitation.lower() == 'pseudorandom':
1616 if num_signals > 1:
1617 print('Warning: Pseudorandom generally not recommended for multi-reference excitation.')
1618 kwargs = {'fft_lines': frame_length // 2,
1619 'f_nyq': bandwidth,
1620 'signal_rms': excitation_level,
1621 'min_freq': excitation_min_frequency,
1622 'max_freq': excitation_max_frequency,
1623 'integration_oversample': integration_oversample,
1624 'averages': num_averages + int(np.ceil(steady_state_time / frame_time))}
1625 kwargs.update(generator_kwargs)
1626 signals = np.array([
1627 generator.pseudorandom(**kwargs)[1]
1628 for i in range(num_signals)
1629 ])
1630 elif excitation.lower() == 'random':
1631 kwargs = {'shape': (num_signals,),
1632 'n_samples': frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate),
1633 'rms': excitation_level,
1634 'dt': dt,
1635 'low_frequency_cutoff': excitation_min_frequency,
1636 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
1637 kwargs.update(generator_kwargs)
1638 signals = generator.random(**kwargs)
1639 elif excitation.lower() == 'burst random':
1640 kwargs = {'shape': (num_signals,),
1641 'n_samples': frame_length * integration_oversample,
1642 'on_fraction': signal_fraction,
1643 'delay_fraction': 0,
1644 'rms': excitation_level,
1645 'dt': dt,
1646 'low_frequency_cutoff': excitation_min_frequency,
1647 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
1648 kwargs.update(generator_kwargs)
1649 signal_list = [generator.burst_random(**kwargs) for i in range(num_averages)]
1650 full_list = []
1651 for i, signal in enumerate(signal_list):
1652 full_list.append(
1653 np.zeros((num_signals, int(extra_time_between_frames * sample_rate))))
1654 full_list.append(signal)
1655 signals = np.concatenate(full_list, axis=-1)
1656 elif excitation.lower() == 'chirp':
1657 if num_signals > 1:
1658 print('Warning: Chirp generally not recommended for multi-reference excitation.')
1659 kwargs = {'frequency_min': 0 if excitation_min_frequency is None else excitation_min_frequency,
1660 'frequency_max': bandwidth if excitation_max_frequency is None else excitation_max_frequency,
1661 'signal_length': frame_time,
1662 'dt': dt}
1663 kwargs.update(generator_kwargs)
1664 signals = np.array([
1665 generator.chirp(**kwargs)
1666 for i in range(num_signals)
1667 ]) * excitation_level
1668 signals = np.tile(signals, [1, num_averages +
1669 int(np.ceil(steady_state_time / frame_time))])
1670 elif excitation.lower() == 'hammer':
1671 if num_signals > 1:
1672 print(
1673 'Warning: Hammer impact generally not recommended for multi-reference excitation, consider multi-hammer instead')
1674 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
1675 signal_length = int(frame_length * integration_oversample * num_averages + (num_averages + 1)
1676 * extra_time_between_frames * sample_rate + 2 * pulse_width * sample_rate)
1677 pulse_times = np.arange(num_averages)[
1678 :, np.newaxis] * (frame_time + extra_time_between_frames) + pulse_width + extra_time_between_frames
1679 kwargs = {'signal_length': signal_length,
1680 'pulse_time': pulse_times,
1681 'pulse_width': pulse_width,
1682 'pulse_peak': excitation_level,
1683 'dt': dt,
1684 'sine_exponent': 2}
1685 kwargs.update(generator_kwargs)
1686 signals = generator.pulse(**kwargs)
1687 signals = np.tile(signals, [num_signals, 1])
1688 elif excitation.lower() == 'multi-hammer':
1689 signal_length = frame_length * integration_oversample
1690 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
1691 signals = []
1692 for i in range(num_signals):
1693 signals.append([])
1694 for j in range(num_averages):
1695 pulse_times = []
1696 last_pulse = 0
1697 while last_pulse < frame_time * signal_fraction:
1698 next_pulse = last_pulse + pulse_width * (np.random.rand() * multihammer_impact_spacing_factor + 1)
1699 pulse_times.append(next_pulse)
1700 last_pulse = next_pulse
1701 pulse_times = np.array(pulse_times)
1702 pulse_times = pulse_times[pulse_times <
1703 frame_time * signal_fraction, np.newaxis]
1704 kwargs = {'signal_length': signal_length,
1705 'pulse_time': pulse_times,
1706 'pulse_width': pulse_width,
1707 'pulse_peak': excitation_level,
1708 'dt': dt,
1709 'sine_exponent': 2}
1710 kwargs.update(generator_kwargs)
1711 signal = generator.pulse(**kwargs)
1712 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
1713 signals[-1].append(signal)
1714 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
1715 signals[-1] = np.concatenate(signals[-1], axis=-1)
1716 signals = np.array(signals)
1717 elif excitation.lower() == 'sine':
1718 if num_signals > 1:
1719 print(
1720 'Warning: Sine signal generally not recommended for multi-reference excitation')
1721 frequencies = excitation_max_frequency if excitation_min_frequency is None else excitation_min_frequency
1722 num_samples = frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate)
1723 kwargs = {'frequencies': frequencies,
1724 'dt': dt,
1725 'num_samples': num_samples,
1726 'amplitudes': excitation_level}
1727 kwargs.update(generator_kwargs)
1728 signals = np.tile(generator.sine(**kwargs), (num_signals, 1))
1729 # Set up the integration
1730 from .sdynpy_data import time_history_array
1731 abscissa = np.arange(signals.shape[-1])*dt
1732 references = time_history_array(abscissa,signals,references[:,np.newaxis]).copy()
1733 responses = self.time_integrate(
1734 references, responses, extra_outputs = extra_outputs)
1735 # Now add noise
1736 responses.ordinate += response_noise_level * np.random.randn(*responses.ordinate.shape)
1737 references.ordinate += excitation_noise_level * np.random.randn(*references.ordinate.shape)
1738 # Filter with antialiasing filters, divide filter order by 2 because of filtfilt
1739 if antialias_filter_order > 0:
1740 lowpass_b, lowpass_a = butter(antialias_filter_order // 2,
1741 antialias_filter_cutoff_factor * bandwidth, fs=sample_rate)
1742 responses.ordinate = filtfilt(lowpass_b, lowpass_a, responses.ordinate)
1743 references.ordinate = filtfilt(lowpass_b, lowpass_a, references.ordinate)
1744 if integration_oversample > 1:
1745 responses = responses.downsample(integration_oversample)
1746 references = references.downsample(integration_oversample)
1747 responses = responses.extract_elements_by_abscissa(steady_state_time, np.inf)
1748 references = references.extract_elements_by_abscissa(steady_state_time, np.inf)
1749 return responses, references
1751 # def reduce_serep(self,shapes_full,shapes_reduced):
1752 # if isinstance(shapes_full,)
1755substructure_by_position = System.substructure_by_position