Coverage for src / sdynpy / signal_processing / sdynpy_integration.py: 14%

1# -*- coding: utf-8 -*- 

2""" 

3Functions for integrating equations of motion in time. 

4""" 

5""" 

6Copyright 2022 National Technology & Engineering Solutions of Sandia, 

7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S. 

8Government retains certain rights in this software. 

9 

10This program is free software: you can redistribute it and/or modify 

11it under the terms of the GNU General Public License as published by 

12the Free Software Foundation, either version 3 of the License, or 

13(at your option) any later version. 

14 

15This program is distributed in the hope that it will be useful, 

16but WITHOUT ANY WARRANTY; without even the implied warranty of 

17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

18GNU General Public License for more details. 

19 

20You should have received a copy of the GNU General Public License 

21along with this program. If not, see <https://www.gnu.org/licenses/>. 

22""" 

23 

24import numpy as np 

25import scipy.signal as signal 

26 

27 

28def MCK_to_StateSpace(M, C, K): 

29 

30 ndofs = M.shape[0] 

31 

32 # A = [[ 0, I], 

33 # [M^-1*K,M^-1*C]] 

34 

35 A_state = np.block([[np.zeros((ndofs, ndofs)), np.eye(ndofs)], 

36 [-np.linalg.solve(M, K), -np.linalg.solve(M, C)]]) 

37 

38 # B = [[ 0, M^-1]] 

39 

40 B_state = np.block([[np.zeros((ndofs, ndofs))], [np.linalg.inv(M)]]) 

41 

42 # C = [[ I, 0], # Displacements 

43 # [ 0, I], # Velocities 

44 # [M^-1*K,M^-1*C], # Accelerations 

45 # [ 0, 0]] # Forces 

46 

47 C_state = np.block([[np.eye(ndofs), np.zeros((ndofs, ndofs))], 

48 [np.zeros((ndofs, ndofs)), np.eye(ndofs)], 

49 [-np.linalg.solve(M, K), -np.linalg.solve(M, C)], 

50 [np.zeros((ndofs, ndofs)), np.zeros((ndofs, ndofs))]]) 

51 

52 # D = [[ 0], # Displacements 

53 # [ 0], # Velocities 

54 # [ M^-1], # Accelerations 

55 # [ I]] # Forces 

56 

57 D_state = np.block([[np.zeros((ndofs, ndofs))], 

58 [np.zeros((ndofs, ndofs))], 

59 [np.linalg.inv(M)], 

60 [np.eye(ndofs)]]) 

61 

62 return A_state, B_state, C_state, D_state 

63 

64 

65def MCK_incomplete_to_StateSpace(M, C, K, massless_dofs=[]): 

66 nomass = len(massless_dofs) 

67 mass = M.shape[0] - nomass 

68 mass_dofs = np.array([v for v in range(M.shape[0]) if v not in massless_dofs]) 

69 massless_dofs = np.array(massless_dofs) 

70 

71 Kmm = K[mass_dofs[:, np.newaxis], mass_dofs] 

72 Kmn = K[mass_dofs[:, np.newaxis], massless_dofs] 

73 Knm = K[massless_dofs[:, np.newaxis], mass_dofs] 

74 Knn = K[massless_dofs[:, np.newaxis], massless_dofs] 

75 

76 Cmm = C[mass_dofs[:, np.newaxis], mass_dofs] 

77 Cnm = C[massless_dofs[:, np.newaxis], mass_dofs] 

78 Cnn = C[massless_dofs[:, np.newaxis], massless_dofs] 

79 

80 Mmm = M[mass_dofs[:, np.newaxis], mass_dofs] 

81 

82 A = np.block([ 

83 [np.zeros((mass, mass)), np.eye(mass), np.zeros((mass, nomass))], 

84 [-np.linalg.solve(Mmm, Kmm), -np.linalg.solve(Mmm, Cmm), -np.linalg.solve(Mmm, Kmn)], 

85 [-np.linalg.solve(Cnn, Knm), -np.linalg.solve(Cnn, Cnm), -np.linalg.solve(Cnn, Knn)] 

86 ]) 

87 

88 B = np.block([ 

89 [np.zeros((mass, mass)), np.zeros((mass, nomass))], 

90 [np.linalg.inv(Mmm), np.zeros((mass, nomass))], 

91 [np.zeros((nomass, mass)), np.linalg.inv(Cnn)] 

92 ]) 

93 

94 C = np.block([ 

95 [np.eye(mass), np.zeros((mass, mass)), np.zeros((mass, nomass))], # Massive displacements 

96 [np.zeros((mass, mass)), np.eye(mass), np.zeros((mass, nomass))], # Massive velocities 

97 [-np.linalg.solve(Mmm, Kmm), -np.linalg.solve(Mmm, Cmm), - \ 

98 np.linalg.solve(Mmm, Kmn)], # Massive Accelerations 

99 [np.zeros((mass, mass)), np.zeros((mass, mass)), 

100 np.zeros((mass, nomass))], # Massive Inputs 

101 [np.zeros((nomass, mass)), np.zeros((nomass, mass)), 

102 np.eye(nomass)], # Nonmassive displacements 

103 [-np.linalg.solve(Cnn, Knm), -np.linalg.solve(Cnn, Cnm), - \ 

104 np.linalg.solve(Cnn, Knn)], # Nonmassive velocities 

105 [np.zeros((nomass, mass)), np.zeros((nomass, mass)), 

106 np.zeros((nomass, nomass))], # Nonmassive inputs 

107 ]) 

108 

109 D = np.block([ 

110 [np.zeros((mass, mass)), np.zeros((mass, nomass))], # Massive displacements 

111 [np.zeros((mass, mass)), np.zeros((mass, nomass))], # Massive velocities 

112 [np.linalg.inv(Mmm), np.zeros((mass, nomass))], # Massive Accelerations 

113 [np.eye(mass), np.zeros((mass, nomass))], # Massive Inputs 

114 [np.zeros((nomass, mass)), np.zeros((nomass, nomass))], # Nonmassive displacements 

115 [np.zeros((nomass, mass)), np.linalg.inv(Cnn)], # Nonmassive velocities 

116 [np.zeros((nomass, mass)), np.eye(nomass)], # Nonmassive inputs 

117 ]) 

118 

119 doflist = {} 

120 doflist['state'] = [(dof, modifier) for modifier in ['', 'd'] 

121 for dof in mass_dofs] + [(dof, '') for dof in massless_dofs] 

122 doflist['input'] = [(dof, 'in') for dof in mass_dofs] + [(dof, 'in') for dof in massless_dofs] 

123 doflist['output'] = [(dof, modifier) for modifier in ['', 'd', 'dd', 'in'] for dof in mass_dofs] + [ 

124 (dof, modifier) for modifier in ['', 'd', 'in'] for dof in massless_dofs] 

125 

126 return A, B, C, D, doflist 

127 

128 

129def integrate_MCK(M, C, K, times, forces, x0=None): 

130 A, B, C, D = MCK_to_StateSpace(M, C, K) 

131 linear_system = signal.StateSpace(A, B, C, D) 

132 

133 times_out, sys_out, x_out = signal.lsim(linear_system, forces, times, x0) 

134 

135 sys_disps = sys_out[:, 0 * M.shape[0]:1 * M.shape[0]] 

136 sys_vels = sys_out[:, 1 * M.shape[0]:2 * M.shape[0]] 

137 sys_accels = sys_out[:, 2 * M.shape[0]:3 * M.shape[0]] 

138 sys_forces = sys_out[:, 3 * M.shape[0]:4 * M.shape[0]] 

139 

140 return sys_disps, sys_vels, sys_accels, sys_forces 

141 

142 

143def modal_parameters_to_MCK(natural_frequencies, damping_ratios): 

144 nmodes = natural_frequencies.size 

145 M = np.eye(nmodes) 

146 K = np.diag((natural_frequencies * 2 * np.pi)**2) 

147 C = np.diag((2 * 2 * np.pi * natural_frequencies * damping_ratios)) 

148 return M, C, K 

149 

150 

151def integrate_modes(natural_frequencies, damping_ratios, modal_forces, times, return_accel=False): 

152 M, C, K = modal_parameters_to_MCK(natural_frequencies, damping_ratios) 

153 q = integrate_MCK(M, C, K, times, modal_forces) 

154 if return_accel: 

155 return q[2] 

156 else: 

157 return q[0] 

158 

159 

160def frequency_domain_differentiation(abscissa, signal, order=1, axis=-1, fmin=-1, fmax=np.inf): 

161 dt = np.mean(np.diff(abscissa)) 

162 n = signal.shape[axis] 

163 freqs = np.fft.fftfreq(n, dt) 

164 frequency_indices = np.logical_and(np.abs(freqs) > fmin, np.abs(freqs) < fmax) 

165 fft = np.fft.fft(signal, axis=axis) 

166 # Construct tuples for multidimensional analysis 

167 slice_array = tuple([frequency_indices if i == axis or (i - axis) == 

168 fft.ndim else np.newaxis for i in range(fft.ndim)]) 

169 freq_slice_array = tuple([frequency_indices if i == axis or ( 

170 i - axis) == fft.ndim else slice(None, None, None) for i in range(fft.ndim)]) 

171 not_freq_slice_array = tuple([np.logical_not(frequency_indices) if i == axis or ( 

172 i - axis) == fft.ndim else slice(None, None, None) for i in range(fft.ndim)]) 

173 factor = (1j * 2 * np.pi * freqs[slice_array])**order 

174 # Perform integration/differentiation 

175 fft[freq_slice_array] *= factor 

176 fft[not_freq_slice_array] = 0 + 0j 

177 return np.real(np.fft.ifft(fft, axis=axis))