Coverage for src / sdynpy / core / sdynpy_system.py: 40%
707 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +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_exodus import Exodus, ExodusInMemory, reduce_exodus_to_surfaces, read_sierra_matlab_matrix_file, read_sierra_matlab_map_file
33from ..fem.sdynpy_shaker import Shaker4DoF
34from ..signal_processing import frf as spfrf
35from ..signal_processing import generator
36from scipy.linalg import eigh, block_diag, null_space, eig
37from scipy.signal import lsim, StateSpace, resample, butter, filtfilt
38import copy
39import netCDF4 as nc4
40import matplotlib.pyplot as plt
41import warnings
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 data = np.load(filename)
798 return cls(data['coordinate'].view(CoordinateArray), data['mass'],
799 data['stiffness'], data['damping'], data['transformation'],
800 True if 'enforce_symmetry' not in data else data['enforce_symmetry'])
802 def __neg__(self):
803 new_system = copy.deepcopy(self)
804 new_system.mass *= -1
805 new_system.stiffness *= -1
806 new_system.damping *= -1
807 return new_system
809 def __mul__(self, value):
810 new_system = copy.deepcopy(self)
811 new_system.mass *= value
812 new_system.stiffness *= value
813 new_system.damping *= value
814 return new_system
816 def __rmul__(self, value):
817 new_system = copy.deepcopy(self)
818 new_system.mass *= value
819 new_system.stiffness *= value
820 new_system.damping *= value
821 return new_system
823 @classmethod
824 def concatenate(cls, systems, coordinate_node_offset=0):
825 """
826 Combine multiple systems together
828 Parameters
829 ----------
830 systems : iterable of System objects
831 Iterable of Systems that will be concatenated. Matrices will be
832 assembled in block diagonal format
833 coordinate_node_offset : int, optional
834 Offset applied to the coordinates so the nodes do not overlap.
835 The default is 0.
837 Returns
838 -------
839 System
840 A system consisting of the combintation of the provided systems.
842 """
843 coordinates = [system.coordinate.copy() for system in systems]
844 if coordinate_node_offset != 0:
845 for i in range(len(coordinates)):
846 coordinates[i].node += coordinate_node_offset * (i + 1)
847 all_coordinates = np.concatenate(coordinates)
848 return cls(all_coordinates,
849 block_diag(*[system.mass for system in systems]),
850 block_diag(*[system.stiffness for system in systems]),
851 block_diag(*[system.damping for system in systems]),
852 block_diag(*[system.transformation for system in systems]),
853 enforce_symmetry = True if all([s._enforce_symmetry for s in systems]) else False)
855 @classmethod
856 def substructure_by_position(cls, systems, geometries, distance_threshold=1e-8, rcond=None):
857 """
858 Applies constraints to systems by constraining colocated nodes together
860 Parameters
861 ----------
862 systems : iterable of System objects
863 A set of systems that will be combined and constrained together
864 geometries : iterable of Geometry objects
865 A set of geometries that will be combined together
866 distance_threshold : float, optional
867 The distance between nodes that are considered colocated.
868 The default is 1e-8.
869 rcond : float, optional
870 Condition number to use in the nullspace calculation on
871 the constraint matrix. The default is None.
873 Returns
874 -------
875 combined_system : System
876 System consisting of constraining the input systems together.
877 combined_geometry : Geometry
878 Combined geometry of the new system
880 """
881 from .sdynpy_geometry import Geometry
882 combined_geometry, node_offset = Geometry.overlay_geometries(
883 geometries, return_node_id_offset=True)
884 combined_system = cls.concatenate(systems, node_offset)
885 global_coords = combined_geometry.global_node_coordinate()
886 node_distances = np.linalg.norm(
887 global_coords[:, np.newaxis, :] - global_coords[np.newaxis, :, :], axis=-1)
888 # Find locations where the value is less than the tolerances, except for on the centerline
889 node_pairs = [[combined_geometry.node.id[index] for index in pair]
890 for pair in zip(*np.where(node_distances < distance_threshold)) if pair[0] < pair[1]]
891 # Find matching DoF pairs
892 constraint_matrix = []
893 for node1, node2 in node_pairs:
894 # Find the dofs associated with each node
895 system_1_dof_indices = np.where(combined_system.coordinate.node == node1)[0]
896 system_1_dofs = combined_system.coordinate[system_1_dof_indices]
897 system_1_transformation = combined_system.transformation[system_1_dof_indices]
898 global_deflections_1 = combined_geometry.global_deflection(system_1_dofs)
899 system_2_dof_indices = np.where(combined_system.coordinate.node == node2)[0]
900 system_2_dofs = combined_system.coordinate[system_2_dof_indices]
901 system_2_transformation = combined_system.transformation[system_2_dof_indices]
902 global_deflections_2 = combined_geometry.global_deflection(system_2_dofs)
903 # Split between translations and rotations
904 translation_map_1 = np.where((abs(system_1_dofs.direction) <= 3)
905 & (abs(system_1_dofs.direction) > 0))
906 rotation_map_1 = np.where(abs(system_1_dofs.direction) > 3)
907 translation_map_2 = np.where((abs(system_2_dofs.direction) <= 3)
908 & (abs(system_2_dofs.direction) > 0))
909 rotation_map_2 = np.where(abs(system_2_dofs.direction) > 3)
910 neutral_map_1 = np.where(abs(system_1_dofs.direction) == 0)
911 neutral_map_2 = np.where(abs(system_2_dofs.direction) == 0)
912 # Do translations and rotations separately
913 for map_1, map_2 in [[translation_map_1, translation_map_2],
914 [rotation_map_1, rotation_map_2],
915 [neutral_map_1, neutral_map_2]]:
916 deflections_1 = global_deflections_1[map_1].T
917 deflections_2 = global_deflections_2[map_2].T
918 transform_1 = system_1_transformation[map_1]
919 transform_2 = system_2_transformation[map_2]
920 full_constraint = deflections_1 @ transform_1 - deflections_2 @ transform_2
921 constraint_matrix.append(full_constraint)
922 constraint_matrix = np.concatenate(constraint_matrix, axis=0)
923 return combined_system.constrain(constraint_matrix, rcond), combined_geometry
925 def constrain(self, constraint_matrix, rcond=None):
926 """
927 Apply a constraint matrix to the system
929 Parameters
930 ----------
931 constraint_matrix : np.ndarray
932 A matrix of constraints to apply to the structure (B matrix in
933 substructuring literature)
934 rcond : float, optional
935 Condition tolerance for computing the nullspace. The default is None.
937 Returns
938 -------
939 System
940 Constrained system.
942 """
943 substructuring_transform_matrix = null_space(constraint_matrix, rcond)
944 new_mass = substructuring_transform_matrix.T @ self.mass @ substructuring_transform_matrix
945 new_stiffness = substructuring_transform_matrix.T @ self.stiffness @ substructuring_transform_matrix
946 new_damping = substructuring_transform_matrix.T @ self.damping @ substructuring_transform_matrix
947 new_transform = self.transformation @ substructuring_transform_matrix
948 return System(self.coordinate, new_mass, new_stiffness, new_damping, new_transform,
949 enforce_symmetry=self._enforce_symmetry)
951 def transformation_matrix_at_coordinates(self, coordinates):
952 """
953 Return the transformation matrix at the specified coordinates
955 Parameters
956 ----------
957 coordinates : CoordinateArray
958 coordinates at which the transformation matrix will be computed.
960 Raises
961 ------
962 ValueError
963 Raised if duplicate coordinates are requested, or if coordinates that
964 do not exist in the system are requested.
966 Returns
967 -------
968 return_value : np.ndarray
969 Portion of the transformation matrix corresponding to the
970 coordinates input to the function.
972 """
973 consistent_arrays, shape_indices, request_indices = np.intersect1d(
974 abs(self.coordinate), abs(coordinates), assume_unique=False, return_indices=True)
975 # Make sure that all of the keys are actually in the consistent array matrix
976 if consistent_arrays.size != coordinates.size:
977 extra_keys = np.setdiff1d(abs(coordinates), abs(self.coordinate))
978 if extra_keys.size == 0:
979 raise ValueError(
980 'Duplicate coordinate values requested. Please ensure coordinate indices are unique.')
981 raise ValueError(
982 'Not all indices in requested coordinate array exist in the system\n{:}'.format(str(extra_keys)))
983 # Handle sign flipping
984 multiplications = coordinates.flatten()[request_indices].sign(
985 ) * self.coordinate[shape_indices].sign()
986 return_value = self.transformation[shape_indices] * multiplications[:, np.newaxis]
987 # Invert the indices to return the dofs in the correct order as specified in keys
988 inverse_indices = np.zeros(request_indices.shape, dtype=int)
989 inverse_indices[request_indices] = np.arange(len(request_indices))
990 return_value = return_value[inverse_indices]
991 return return_value
993 def substructure_by_coordinate(self, dof_pairs, rcond=None,
994 return_constrained_system=True):
995 """
996 Constrain the system by connecting the specified degree of freedom pairs
998 Parameters
999 ----------
1000 dof_pairs : iterable of CoordinateArray
1001 Pairs of coordinates to be connected. None can be passed instead of
1002 a second degree of freedom to constrain to ground
1003 rcond : float, optional
1004 Condition threshold to use for the nullspace calculation on the
1005 constraint matrix. The default is None.
1006 return_constrained_system : bool, optional
1007 If true, apply the constraint matrix and return the constrained
1008 system, otherwise simply return the constraint matrix. The default
1009 is True.
1011 Returns
1012 -------
1013 np.ndarray or System
1014 Returns a System object with the constraints applied if
1015 `return_constrained_system` is True, otherwise just return the
1016 constraint matrix.
1018 """
1019 constraint_matrix = []
1020 for constraint_dof_0, constraint_dof_1 in dof_pairs:
1021 constraint = self.transformation_matrix_at_coordinates(constraint_dof_0)
1022 if constraint_dof_1 is not None:
1023 constraint -= self.transformation_matrix_at_coordinates(constraint_dof_1)
1024 constraint_matrix.append(constraint)
1025 constraint_matrix = np.concatenate(constraint_matrix, axis=0)
1026 if return_constrained_system:
1027 return self.constrain(constraint_matrix, rcond)
1028 else:
1029 return constraint_matrix
1031 def substructure_by_shape(self, constraint_shapes, connection_dofs_0,
1032 connection_dofs_1=None, rcond=None,
1033 return_constrained_system=True):
1034 """
1035 Constrain the system using a set of shapes in a least-squares sense.
1037 Parameters
1038 ----------
1039 constraint_shapes : ShapeArray
1040 An array of shapes to use as the basis for the constraints
1041 connection_dofs_0 : CoordinateArray
1042 Array of coordinates to use in the constraints
1043 connection_dofs_1 : CoordinateArray, optional
1044 Array of coordinates to constrain to the coordinates in
1045 `connection_dofs_0`. If not specified, the `connection_dofs_0`
1046 degrees of freedom will be constrained to ground.
1047 rcond : float, optional
1048 Condition threshold on the nullspace calculation. The default is None.
1049 return_constrained_system : bool, optional
1050 If true, apply the constraint matrix and return the constrained
1051 system, otherwise simply return the constraint matrix. The default
1052 is True.
1054 Returns
1055 -------
1056 np.ndarray or System
1057 Returns a System object with the constraints applied if
1058 `return_constrained_system` is True, otherwise just return the
1059 constraint matrix.
1061 """
1062 shape_matrix_0 = constraint_shapes[connection_dofs_0].T
1063 transform_matrix_0 = self.transformation_matrix_at_coordinates(connection_dofs_0)
1064 constraint_matrix = np.linalg.lstsq(shape_matrix_0, transform_matrix_0)[0]
1065 if connection_dofs_1 is not None:
1066 shape_matrix_1 = constraint_shapes[connection_dofs_1].T
1067 transform_matrix_1 = self.transformation_matrix_at_coordinates(connection_dofs_1)
1068 constraint_matrix -= np.linalg.lstsq(shape_matrix_1, transform_matrix_1)[0]
1069 if return_constrained_system:
1070 return self.constrain(constraint_matrix, rcond)
1071 else:
1072 return constraint_matrix
1074 def substructure_shakers(self,shakers,coordinates):
1075 """
1076 Uses substructuring to attach 4DoF Shaker Models to the Substructure
1078 Parameters
1079 ----------
1080 shakers : array of Shaker4DoF
1081 4-DoF electromechanical models of a shaker that are added to the
1082 System
1083 coordinates : CoordinateArray
1084 The degrees of freedom at which the shakers should be attached.
1086 Returns
1087 -------
1088 System
1089 A system with shakers attached to it. The shaker degrees of freedom
1090 will be assigned based on the next power of 10 bigger than the
1091 current maximum node number. Degrees of freedom end in 1 for
1092 armature force/displacement, 2 for body force/displacement, 3 for
1093 force gauge force/displacement, and 4 for drive voltage/current
1094 """
1095 constraint_pairs = []
1096 shaker_systems = []
1098 max_node = self.coordinate.node.max()
1099 shaker_node_offset = int(10**np.ceil(np.log10(max_node)))
1100 if shaker_node_offset == max_node:
1101 shaker_node_offset *= 10
1102 if shaker_node_offset < 100:
1103 shaker_node_offset = 100
1105 for i,(direction,shaker) in enumerate(zip(coordinates,shakers)):
1106 M,C,K = shaker.MCK()
1107 coords = coordinate_array(shaker_node_offset + i*10 + np.arange(4)+1,[1,1,1,0])
1108 shaker_systems.append(System(coords,M,K,C,enforce_symmetry = False))
1109 constraint_dof = coordinate_array(shaker_node_offset + i*10 + 3,1)
1110 constraint_pairs.append([direction,constraint_dof])
1112 constraint_pairs = np.array(constraint_pairs).view(CoordinateArray)
1114 combined_system = System.concatenate([self]+shaker_systems)
1116 constrained_system = combined_system.substructure_by_coordinate(constraint_pairs)
1118 return constrained_system
1120 def copy(self):
1121 """
1122 Returns a copy of the system object
1124 Returns
1125 -------
1126 System
1127 A copy of the system object
1129 """
1130 return System(self.coordinate.copy(), self.mass.copy(), self.stiffness.copy(),
1131 self.damping.copy(), self.transformation.copy(),
1132 self._enforce_symmetry)
1134 def set_proportional_damping(self, mass_fraction, stiffness_fraction):
1135 """
1136 Sets the damping matrix to a proportion of the mass and stiffness matrices.
1138 The damping matrix will be set to `mass_fraction*self.mass +
1139 stiffness_fraction*self.stiffness`
1141 Parameters
1142 ----------
1143 mass_fraction : float
1144 Fraction of the mass matrix
1145 stiffness_fraction : TYPE
1146 Fraction of the stiffness matrix
1148 Returns
1149 -------
1150 None.
1152 """
1153 self.damping = self.mass * mass_fraction + self.stiffness * stiffness_fraction
1155 @classmethod
1156 def beam(cls, length, width, height, num_nodes, E=None, rho=None, nu=None, material=None):
1157 """
1158 Create a beam mass and stiffness matrix
1160 Parameters
1161 ----------
1162 length : float
1163 Lenghth of the beam
1164 width : float
1165 Width of the beam
1166 height : float
1167 Height of the beam
1168 num_nodes : int
1169 Number of nodes in the beam.
1170 E : float, optional
1171 Young's modulus of the beam. If not specified, a `material` must be
1172 specified instead
1173 rho : float, optional
1174 Density of the beam. If not specified, a `material` must be
1175 specified instead
1176 nu : float, optional
1177 Poisson's ratio of the beam. If not specified, a `material` must be
1178 specified instead
1179 material : str, optional
1180 A specific material can be specified instead of `E`, `rho`, and `nu`.
1181 Should be a string 'steel' or 'aluminum'. If not specified, then
1182 options `E`, `rho`, and `nu` must be specified instead.
1184 Raises
1185 ------
1186 ValueError
1187 If improper materials are defined.
1189 Returns
1190 -------
1191 system : System
1192 A system object consisting of the beam mass and stiffness matrices.
1193 geometry : Geometry
1194 A Geometry consisting of the beam geometry.
1196 """
1197 from .sdynpy_geometry import Geometry, node_array, traceline_array, coordinate_system_array
1198 node_positions = np.array((np.linspace(0, length, num_nodes),
1199 np.zeros(num_nodes),
1200 np.zeros(num_nodes))).T
1201 node_connectivity = np.array((np.arange(num_nodes - 1), np.arange(1, num_nodes))).T
1202 bend_direction_1 = np.array((np.zeros(num_nodes - 1),
1203 np.zeros(num_nodes - 1),
1204 np.ones(num_nodes - 1))).T
1205 if material is None:
1206 if E is None or rho is None or nu is None:
1207 raise ValueError('Must specify material or E, nu, and rho')
1208 elif material.lower() == 'steel':
1209 E = 200e9 # [N/m^2],
1210 nu = 0.25 # [-],
1211 rho = 7850 # [kg/m^3]
1212 elif material.lower() == 'aluminum':
1213 E = 69e9 # [N/m^2],
1214 nu = 0.33 # [-],
1215 rho = 2830 # [kg/m^3]
1216 else:
1217 raise ValueError('Unknown Material {:}'.format(material))
1218 mat_props = rect_beam_props(E, rho, nu, width, height, num_nodes - 1)
1219 K, M = beamkm(node_positions, node_connectivity, bend_direction_1, **mat_props)
1220 coordinates = from_nodelist(np.arange(num_nodes) + 1, directions=[1, 2, 3, 4, 5, 6])
1221 system = cls(coordinates, M, K)
1222 nodelist = node_array(np.arange(num_nodes) + 1, node_positions)
1223 tracelines = traceline_array(connectivity=np.arange(num_nodes) + 1)
1224 coordinate_systems = coordinate_system_array()
1225 geometry = Geometry(nodelist, coordinate_systems, tracelines)
1226 return system, geometry
1228 def get_indices_by_coordinate(self, coordinates, ignore_sign=False):
1229 """
1230 Gets the indices in the transformation matrix corresponding coordinates
1232 Parameters
1233 ----------
1234 coordinates : CoordinateArray
1235 Coordinates to extract transformation indices
1236 ignore_sign : bool, optional
1237 Specify whether or not to ignore signs on the coordinates. If True,
1238 then '101X+' would match '101X+' or '101X-'. The default is False.
1240 Raises
1241 ------
1242 ValueError
1243 Raised if duplicate coordinates or coordinates not in the system
1244 are requested
1246 Returns
1247 -------
1248 np.ndarray
1249 Array of indices.
1251 """
1252 if ignore_sign:
1253 consistent_arrays, shape_indices, request_indices = np.intersect1d(
1254 abs(self.coordinate), abs(coordinates), assume_unique=False, return_indices=True)
1255 else:
1256 consistent_arrays, shape_indices, request_indices = np.intersect1d(
1257 self.coordinate, coordinates, assume_unique=False, return_indices=True)
1258 # Make sure that all of the keys are actually in the consistent array matrix
1259 if consistent_arrays.size != coordinates.size:
1260 extra_keys = np.setdiff1d(abs(coordinates), abs(self.coordinate))
1261 if extra_keys.size == 0:
1262 raise ValueError(
1263 'Duplicate coordinate values requested. Please ensure coordinate indices are unique.')
1264 raise ValueError(
1265 'Not all indices in requested coordinate array exist in the system\n{:}'.format(str(extra_keys)))
1266 # Handle sign flipping
1267 return_value = shape_indices
1268 # Invert the indices to return the dofs in the correct order as specified in keys
1269 inverse_indices = np.zeros(request_indices.shape, dtype=int)
1270 inverse_indices[request_indices] = np.arange(len(request_indices))
1271 return return_value[inverse_indices]
1273 def reduce(self, reduction_transformation):
1274 """
1275 Apply the specified reduction to the model
1277 Parameters
1278 ----------
1279 reduction_transformation : np.ndarray
1280 Matrix to use in the reduction
1282 Returns
1283 -------
1284 System
1285 Reduced system.
1287 """
1288 mass = reduction_transformation.T @ self.mass @ reduction_transformation
1289 stiffness = reduction_transformation.T @ self.stiffness @ reduction_transformation
1290 damping = reduction_transformation.T @ self.damping @ reduction_transformation
1291 transformation = self.transformation @ reduction_transformation
1292 # Force symmetry
1293 mass = (mass + mass.T) / 2
1294 stiffness = (stiffness + stiffness.T) / 2
1295 damping = (damping + damping.T) / 2
1296 return System(self.coordinate, mass, stiffness, damping, transformation, self._enforce_symmetry)
1298 def reduce_guyan(self, coordinates):
1299 """
1300 Perform Guyan reduction on the system
1302 Parameters
1303 ----------
1304 coordinates : CoordinateArray
1305 A list of coordinates to keep in the reduced system.
1307 Raises
1308 ------
1309 ValueError
1310 Raised the transformation matrix is not identity matrix.
1312 Returns
1313 -------
1314 System
1315 Reduced system.
1317 """
1318 if isinstance(coordinates, CoordinateArray):
1319 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1320 raise ValueError(
1321 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1322 keep_dofs = self.get_indices_by_coordinate(coordinates)
1323 else:
1324 keep_dofs = np.array(coordinates)
1325 discard_dofs = np.array([i for i in range(self.ndof) if i not in keep_dofs])
1326 I_a = np.eye(keep_dofs.size)
1327 K_dd = self.stiffness[discard_dofs[:, np.newaxis],
1328 discard_dofs]
1329 K_da = self.stiffness[discard_dofs[:, np.newaxis],
1330 keep_dofs]
1331 T_guyan = np.concatenate((I_a, -np.linalg.solve(K_dd, K_da)), axis=0)
1332 T_guyan[np.concatenate((keep_dofs, discard_dofs)), :] = T_guyan.copy()
1333 return self.reduce(T_guyan)
1335 def reduce_dynamic(self, coordinates, frequency):
1336 """
1337 Perform Dynamic condensation
1339 Parameters
1340 ----------
1341 coordinates : CoordinateArray
1342 A list of coordinates to keep in the reduced system.
1343 frequency : float
1344 The frequency to preserve in the dynamic reduction.
1346 Raises
1347 ------
1348 ValueError
1349 Raised if the transformation is not identity matrix.
1351 Returns
1352 -------
1353 System
1354 Reduced system.
1356 """
1357 if isinstance(coordinates, CoordinateArray):
1358 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1359 raise ValueError(
1360 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1361 keep_dofs = self.get_indices_by_coordinate(coordinates)
1362 else:
1363 keep_dofs = np.array(coordinates)
1364 discard_dofs = np.array([i for i in range(self.ndof) if i not in keep_dofs])
1365 I_a = np.eye(keep_dofs.size)
1366 D = self.stiffness - (2 * np.pi * frequency)**2 * self.mass
1367 D_dd = D[discard_dofs[:, np.newaxis],
1368 discard_dofs]
1369 D_da = D[discard_dofs[:, np.newaxis],
1370 keep_dofs]
1371 T_dynamic = np.concatenate((I_a, -np.linalg.solve(D_dd, D_da)), axis=0)
1372 T_dynamic[np.concatenate((keep_dofs, discard_dofs)), :] = T_dynamic.copy()
1373 return self.reduce(T_dynamic)
1375 def reduce_craig_bampton(self, connection_degrees_of_freedom: CoordinateArray,
1376 num_fixed_base_modes: int,
1377 return_shape_matrix: bool = False):
1378 """
1379 Computes a craig-bampton substructure model for the system
1381 Parameters
1382 ----------
1383 connection_degrees_of_freedom : CoordinateArray
1384 Degrees of freedom to keep at the interface.
1385 num_fixed_base_modes : int
1386 Number of fixed-base modes to use in the reduction
1387 return_shape_matrix : bool, optional
1388 If true, return a set of shapes that represents the transformation
1389 in addition to the reduced system. The default is False.
1391 Raises
1392 ------
1393 ValueError
1394 Raised if coordinate arrays are specified when there is already a
1395 transformation.
1397 Returns
1398 -------
1399 System
1400 Reduced system in craig-bampton form
1401 ShapeArray
1402 Shapes representing the craig-bampton transformation
1404 """
1405 # Construct craig bampton transformation
1406 if isinstance(connection_degrees_of_freedom, CoordinateArray):
1407 if not np.allclose(self.transformation, np.eye(*self.transformation.shape)):
1408 raise ValueError(
1409 'Coordinates can only be specified with a CoordinateArray if the transformation is identity')
1410 connection_indices = self.get_indices_by_coordinate(connection_degrees_of_freedom)
1411 else:
1412 connection_indices = np.array(connection_degrees_of_freedom)
1413 other_indices = np.array([i for i in range(self.ndof) if i not in connection_indices])
1415 # Extract portions of the mass and stiffness matrices
1416 K_ii = self.K[other_indices[:, np.newaxis], other_indices]
1417 M_ii = self.M[other_indices[:, np.newaxis], other_indices]
1418 K_ib = self.K[other_indices[:, np.newaxis], connection_indices]
1420 # Compute fixed interface modes
1421 lam, Phi_ii = eigh(K_ii, M_ii, subset_by_index=[0, int(num_fixed_base_modes) - 1])
1422 # Normalize the mode shapes
1423 lam[lam < 0] = 0
1424 normalized_mass = np.diag(Phi_ii.T @ M_ii @ Phi_ii)
1425 Phi_ii /= np.sqrt(normalized_mass)
1426 Z_bi = np.zeros((connection_indices.size, num_fixed_base_modes))
1427 # Compute constraint modes
1428 Psi_ib = -np.linalg.solve(K_ii, K_ib)
1429 I_bb = np.eye(connection_indices.size)
1430 T_cb = np.block([[Phi_ii, Psi_ib],
1431 [Z_bi, I_bb]])
1432 T_cb[np.concatenate((other_indices, connection_indices)), :] = T_cb.copy()
1433 if return_shape_matrix:
1434 from .sdynpy_shape import shape_array
1435 freq = np.sqrt(lam) / (2 * np.pi)
1436 all_freqs = np.concatenate((freq, np.zeros(connection_indices.size)))
1437 shapes = shape_array(self.coordinate, T_cb.T, all_freqs, comment1=['Fixed Base Mode {:}'.format(
1438 i + 1) for i in range(num_fixed_base_modes)] + ['Constraint Mode {:}'.format(str(dof)) for dof in connection_degrees_of_freedom])
1439 return self.reduce(T_cb), shapes
1440 else:
1441 return self.reduce(T_cb)
1443 @classmethod
1444 def from_exodus_superelement(cls, superelement_nc4, transformation_exodus_file=None,
1445 x_disp='DispX', y_disp='DispY', z_disp='DispZ',
1446 x_rot=None, y_rot=None, z_rot=None,
1447 reduce_to_external_surfaces=False):
1448 """
1449 Creates a system from a superelement from Sierra/SD
1451 Parameters
1452 ----------
1453 superelement_nc4 : netCDF4.Dataset or string
1454 Dataset from which the superelement data will be loaded
1455 transformation_exodus_file : Exodus, ExodusInMemory, or str, optional
1456 Exodus data containing the transformation between the reduced
1457 superelement state and the physical space. If not specified, no
1458 transformation will be created.
1459 x_disp : str, optional
1460 Variable name to read for x-displacements in the transformation
1461 Exodus file. The default is 'DispX'.
1462 y_disp : str, optional
1463 Variable name to read for y-displacements in the transformation
1464 Exodus file. The default is 'DispY'.
1465 z_disp : str, optional
1466 Variable name to read for z-displacements in the transformation
1467 Exodus file. The default is 'DispZ'.
1468 x_rot : str, optional
1469 Variable name to read for x-rotations in the transformation
1470 Exodus file. The default is to not read rotations.
1471 y_rot : str, optional
1472 Variable name to read for y-rotations in the transformation
1473 Exodus file. The default is to not read rotations.
1474 z_rot : str, optional
1475 Variable name to read for z-rotations in the transformation
1476 Exodus file. The default is to not read rotations.
1477 reduce_to_external_surfaces : bool, optional
1478 If True, exodus results will be reduced to external surfaces
1480 Raises
1481 ------
1482 ValueError
1483 raised if bad data types are passed to the arguments.
1485 Returns
1486 -------
1487 system : System
1488 System containing the superelement representation
1489 geometry : Geometry
1490 Geometry that can be used to plot the system
1491 boundary_dofs : CoordinateArray
1492 Degrees of freedom that can be used to constrain the test article.
1494 """
1495 from .sdynpy_geometry import node_array, coordinate_system_array, Geometry
1496 if isinstance(superelement_nc4, str):
1497 ds = nc4.Dataset(superelement_nc4)
1498 elif isinstance(superelement_nc4, nc4.Dataset):
1499 ds = superelement_nc4
1500 else:
1501 raise ValueError('superelement_nc4 must be a string or a netCDF4 Dataset')
1502 cbmap = ds['cbmap'][:].data.copy()
1503 Kr = ds['Kr'][:].data.copy()
1504 Mr = ds['Mr'][:].data.copy()
1505 Cr = ds['Cr'][:].data.copy()
1506 num_constraint_modes = ds.dimensions['NumConstraints'].size
1507 num_fixed_base_modes = ds.dimensions['NumEig'].size
1508 boundary_dofs = coordinate_array(*cbmap[cbmap[:, 0] > 0].T)
1509 if transformation_exodus_file is None:
1510 transformation = None
1511 coordinate_nodes = cbmap[:, 0]
1512 coordinate_dirs = cbmap[:, 1]
1513 coordinate_nodes[coordinate_nodes == 0] = np.arange(num_fixed_base_modes) + 1
1514 coordinates = coordinate_array(coordinate_nodes, coordinate_dirs)
1515 cs_array = coordinate_system_array()
1516 n_array = node_array(ds['node_num_map'][:], np.array(
1517 [ds['coord{:}'.format(d)][:] for d in 'xyz']).T)
1518 geometry = Geometry(node=n_array, coordinate_system=cs_array)
1519 else:
1520 if isinstance(transformation_exodus_file, str):
1521 exo = Exodus(transformation_exodus_file)
1522 elif isinstance(transformation_exodus_file, Exodus):
1523 exo = transformation_exodus_file
1524 elif isinstance(transformation_exodus_file, ExodusInMemory):
1525 exo = transformation_exodus_file
1526 else:
1527 raise ValueError('transformation_exodus_file must be a string or a sdpy.Exodus')
1528 if reduce_to_external_surfaces:
1529 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])
1530 from .sdynpy_shape import ShapeArray
1531 shapes = ShapeArray.from_exodus(exo, x_disp, y_disp, z_disp, x_rot, y_rot, z_rot)
1532 transformation = shapes.shape_matrix.T
1533 coordinates = shapes[0].coordinate
1534 geometry = Geometry.from_exodus(exo)
1535 system = cls(coordinates, Mr, Kr, Cr, transformation)
1536 return system, geometry, boundary_dofs
1538 @classmethod
1539 def from_sierra_sd_mfile_output(cls, maa_file, kaa_file, gid_file, aset_map_file):
1540 """
1541 Generates a System object from the MFILE, MAA, and KAA outputs from
1542 Sierra/SD
1544 Parameters
1545 ----------
1546 maa_file : str
1547 Path to the file containing the MAA matrix.
1548 kaa_file : str
1549 Path to the file containing the KAA matrix.
1550 gid_file : str
1551 Path to the file contiaining the global node id map
1552 aset_map_file : str
1553 Path to the A-set Map.
1555 Returns
1556 -------
1557 system : System
1558 A system object with the specified matrices and degrees of freedom.
1560 """
1561 M = read_sierra_matlab_matrix_file(maa_file)
1562 K = read_sierra_matlab_matrix_file(kaa_file)
1563 gid = read_sierra_matlab_map_file(gid_file)
1564 asetmap = read_sierra_matlab_map_file(aset_map_file).reshape(-1,9)
1565 valid_dofs = asetmap != -1
1566 # Check if any of the acoustic, temperature, or voltage are defined
1567 if np.any(valid_dofs[:,-3:]):
1568 warnings.warn('Sierra system contains acoustic, temperature, or voltage degrees of freedom. These will be ignored by SDynPy.')
1569 valid_dofs = valid_dofs[:,:-3]
1570 asetmap = asetmap[:,:-3]
1571 coordinates = coordinate_array(gid[:,np.newaxis],[1,2,3,4,5,6])[valid_dofs]
1572 dof_ordering = asetmap[valid_dofs]-1
1573 M = M[dof_ordering[:,np.newaxis],dof_ordering]
1574 K = K[dof_ordering[:,np.newaxis],dof_ordering]
1575 system = cls(coordinates,M,K)
1576 return system
1578 def simulate_test(
1579 self, # The system itself
1580 bandwidth,
1581 frame_length,
1582 num_averages,
1583 excitation,
1584 references,
1585 responses=None, # All Responses
1586 excitation_level=1.0,
1587 excitation_noise_level=0.0,
1588 response_noise_level=0.0,
1589 steady_state_time=0.0,
1590 excitation_min_frequency=None,
1591 excitation_max_frequency=None,
1592 signal_fraction=0.5,
1593 extra_time_between_frames=0.0,
1594 integration_oversample=10,
1595 antialias_filter_cutoff_factor=3,
1596 antialias_filter_order=4,
1597 multihammer_impact_spacing_factor = 4,
1598 extra_outputs = None,
1599 **generator_kwargs
1600 ):
1601 available_excitations = ['pseudorandom', 'random',
1602 'burst random', 'chirp', 'hammer', 'multi-hammer', 'sine']
1603 if not excitation.lower() in available_excitations:
1604 raise ValueError('Excitation must be one of {:}'.format(available_excitations))
1605 # Create the input signal
1606 references = np.atleast_1d(references)
1607 num_signals = references.size
1608 sample_rate = bandwidth * 2 * integration_oversample
1609 dt = 1 / sample_rate
1610 frame_time = dt * frame_length * integration_oversample
1611 # Create the signals
1612 if excitation.lower() == 'pseudorandom':
1613 if num_signals > 1:
1614 print('Warning: Pseudorandom generally not recommended for multi-reference excitation.')
1615 kwargs = {'fft_lines': frame_length // 2,
1616 'f_nyq': bandwidth,
1617 'signal_rms': excitation_level,
1618 'min_freq': excitation_min_frequency,
1619 'max_freq': excitation_max_frequency,
1620 'integration_oversample': integration_oversample,
1621 'averages': num_averages + int(np.ceil(steady_state_time / frame_time))}
1622 kwargs.update(generator_kwargs)
1623 signals = np.array([
1624 generator.pseudorandom(**kwargs)[1]
1625 for i in range(num_signals)
1626 ])
1627 elif excitation.lower() == 'random':
1628 kwargs = {'shape': (num_signals,),
1629 'n_samples': frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate),
1630 'rms': excitation_level,
1631 'dt': dt,
1632 'low_frequency_cutoff': excitation_min_frequency,
1633 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
1634 kwargs.update(generator_kwargs)
1635 signals = generator.random(**kwargs)
1636 elif excitation.lower() == 'burst random':
1637 kwargs = {'shape': (num_signals,),
1638 'n_samples': frame_length * integration_oversample,
1639 'on_fraction': signal_fraction,
1640 'delay_fraction': 0,
1641 'rms': excitation_level,
1642 'dt': dt,
1643 'low_frequency_cutoff': excitation_min_frequency,
1644 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
1645 kwargs.update(generator_kwargs)
1646 signal_list = [generator.burst_random(**kwargs) for i in range(num_averages)]
1647 full_list = []
1648 for i, signal in enumerate(signal_list):
1649 full_list.append(
1650 np.zeros((num_signals, int(extra_time_between_frames * sample_rate))))
1651 full_list.append(signal)
1652 signals = np.concatenate(full_list, axis=-1)
1653 elif excitation.lower() == 'chirp':
1654 if num_signals > 1:
1655 print('Warning: Chirp generally not recommended for multi-reference excitation.')
1656 kwargs = {'frequency_min': 0 if excitation_min_frequency is None else excitation_min_frequency,
1657 'frequency_max': bandwidth if excitation_max_frequency is None else excitation_max_frequency,
1658 'signal_length': frame_time,
1659 'dt': dt}
1660 kwargs.update(generator_kwargs)
1661 signals = np.array([
1662 generator.chirp(**kwargs)
1663 for i in range(num_signals)
1664 ]) * excitation_level
1665 signals = np.tile(signals, [1, num_averages +
1666 int(np.ceil(steady_state_time / frame_time))])
1667 elif excitation.lower() == 'hammer':
1668 if num_signals > 1:
1669 print(
1670 'Warning: Hammer impact generally not recommended for multi-reference excitation, consider multi-hammer instead')
1671 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
1672 signal_length = int(frame_length * integration_oversample * num_averages + (num_averages + 1)
1673 * extra_time_between_frames * sample_rate + 2 * pulse_width * sample_rate)
1674 pulse_times = np.arange(num_averages)[
1675 :, np.newaxis] * (frame_time + extra_time_between_frames) + pulse_width + extra_time_between_frames
1676 kwargs = {'signal_length': signal_length,
1677 'pulse_time': pulse_times,
1678 'pulse_width': pulse_width,
1679 'pulse_peak': excitation_level,
1680 'dt': dt,
1681 'sine_exponent': 2}
1682 kwargs.update(generator_kwargs)
1683 signals = generator.pulse(**kwargs)
1684 signals = np.tile(signals, [num_signals, 1])
1685 elif excitation.lower() == 'multi-hammer':
1686 signal_length = frame_length * integration_oversample
1687 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
1688 signals = []
1689 for i in range(num_signals):
1690 signals.append([])
1691 for j in range(num_averages):
1692 pulse_times = []
1693 last_pulse = 0
1694 while last_pulse < frame_time * signal_fraction:
1695 next_pulse = last_pulse + pulse_width * (np.random.rand() * multihammer_impact_spacing_factor + 1)
1696 pulse_times.append(next_pulse)
1697 last_pulse = next_pulse
1698 pulse_times = np.array(pulse_times)
1699 pulse_times = pulse_times[pulse_times <
1700 frame_time * signal_fraction, np.newaxis]
1701 kwargs = {'signal_length': signal_length,
1702 'pulse_time': pulse_times,
1703 'pulse_width': pulse_width,
1704 'pulse_peak': excitation_level,
1705 'dt': dt,
1706 'sine_exponent': 2}
1707 kwargs.update(generator_kwargs)
1708 signal = generator.pulse(**kwargs)
1709 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
1710 signals[-1].append(signal)
1711 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
1712 signals[-1] = np.concatenate(signals[-1], axis=-1)
1713 signals = np.array(signals)
1714 elif excitation.lower() == 'sine':
1715 if num_signals > 1:
1716 print(
1717 'Warning: Sine signal generally not recommended for multi-reference excitation')
1718 frequencies = excitation_max_frequency if excitation_min_frequency is None else excitation_min_frequency
1719 num_samples = frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate)
1720 kwargs = {'frequencies': frequencies,
1721 'dt': dt,
1722 'num_samples': num_samples,
1723 'amplitudes': excitation_level}
1724 kwargs.update(generator_kwargs)
1725 signals = np.tile(generator.sine(**kwargs), (num_signals, 1))
1726 # Set up the integration
1727 from .sdynpy_data import time_history_array
1728 abscissa = np.arange(signals.shape[-1])*dt
1729 references = time_history_array(abscissa,signals,references[:,np.newaxis]).copy()
1730 responses = self.time_integrate(
1731 references, responses, extra_outputs = extra_outputs)
1732 # Now add noise
1733 responses.ordinate += response_noise_level * np.random.randn(*responses.ordinate.shape)
1734 references.ordinate += excitation_noise_level * np.random.randn(*references.ordinate.shape)
1735 # Filter with antialiasing filters, divide filter order by 2 because of filtfilt
1736 if antialias_filter_order > 0:
1737 lowpass_b, lowpass_a = butter(antialias_filter_order // 2,
1738 antialias_filter_cutoff_factor * bandwidth, fs=sample_rate)
1739 responses.ordinate = filtfilt(lowpass_b, lowpass_a, responses.ordinate)
1740 references.ordinate = filtfilt(lowpass_b, lowpass_a, references.ordinate)
1741 if integration_oversample > 1:
1742 responses = responses.downsample(integration_oversample)
1743 references = references.downsample(integration_oversample)
1744 responses = responses.extract_elements_by_abscissa(steady_state_time, np.inf)
1745 references = references.extract_elements_by_abscissa(steady_state_time, np.inf)
1746 return responses, references
1748 # def reduce_serep(self,shapes_full,shapes_reduced):
1749 # if isinstance(shapes_full,)
1752substructure_by_position = System.substructure_by_position