Coverage for fiqus/parsers/ParserGetDPOnSection.py: 9%

1import numpy as np 

2 

3import re 

4import math 

5from typing import Literal 

6 

7 

8class ParserGetDPOnSection: 

9 """ 

10 This class parses GetDP's TimeTable format output files. 

11 """ 

12 

13 def __init__(self, filePath, data_type: Literal["scalar", "vector"], depth): 

14 self.time_values = [] 

15 self.data_type = data_type 

16 self.depth = depth 

17 

18 if self.depth not in [0, 1]: 

19 raise NotImplementedError("Only depth = 0 and depth = 1 is implemented.") 

20 

21 # Check GMSH documentation for object types: 

22 if self.depth == 0: 

23 if self.data_type == "scalar": 

24 lineName = "SP" # scalar point 

25 elif self.data_type == "vector": 

26 lineName = "VP" # vector point 

27 elif self.depth == 1: 

28 if self.data_type == "scalar": 

29 lineName = "ST" 

30 elif self.data_type == "vector": 

31 lineName = "VT" 

32 

33 # Parse data: 

34 with open(filePath) as file: 

35 data = file.read() 

36 

37 time_values_line = re.search(r"TIME\{(.*)\}", data)[0] 

38 time_values = re.findall(r"TIME\{(.*)\}", time_values_line) 

39 self.time_values = [ 

40 float(time_value) for time_value in time_values[0].split(",") 

41 ] 

42 data = data.replace(time_values_line, "") 

43 

44 points = re.findall(lineName + r"\((.*)\){.*\..*}", data) 

45 points = [point.split(",") for point in points] 

46 self.points = np.array(points, dtype=float) 

47 

48 if self.depth == 1: 

49 length = np.shape(self.points)[1] 

50 step = int(length / 3) 

51 self.points = ( 

52 self.points[:, 0:step] 

53 + self.points[:, step : 2 * step] 

54 + self.points[:, 2 * step : 3 * step] 

55 ) / 3 

56 

57 values = re.findall(lineName + r"\(.*\){(.*\..*)}", data) 

58 values = [value.split(",") for value in values] 

59 self.values = np.array(values, dtype=float) 

60 

61 if self.depth == 1: 

62 length = np.shape(self.values)[1] 

63 step = int(length / 3) 

64 self.values = ( 

65 self.values[:, 0:step] 

66 + self.values[:, step : 2 * step] 

67 + self.values[:, 2 * step : 3 * step] 

68 ) / 3 

69 

70 # Somehow, even with depth 0, we get duplicate magnitudes for different points. 

71 # We need to remove them for better plotting: 

72 _, unique_indices = np.unique( 

73 self.values[:, np.shape(self.values)[1] // 2], return_index=True 

74 ) 

75 self.values = self.values[unique_indices, :] 

76 self.points = self.points[unique_indices, :] 

77 

78 def get_values_at_time(self, time): 

79 """ 

80 Returns the values at the specified time. 

81 """ 

82 if self.points.shape[1] != 2: 

83 raise ValueError( 

84 "Use project_values_on_a_plane() before calling get_values_at_time()!" 

85 ) 

86 

87 for index, time_value in enumerate(self.time_values): 

88 if math.isclose(float(time_value), float(time), abs_tol=1e-10): 

89 return self.values[:, index] 

90 

91 def get_values_at_time_step(self, time_step): 

92 """ 

93 Returns the values at the specified time step. 

94 """ 

95 if self.points.shape[1] != 2: 

96 raise ValueError( 

97 "Use project_values_on_a_plane() before calling get_values_at_time()!" 

98 ) 

99 

100 if self.data_type == "scalar": 

101 return self.values[:, time_step] 

102 elif self.data_type == "vector": 

103 return self.values[:, time_step * 2 : time_step * 2 + 2] 

104 

105 def project_values_on_a_plane(self, plane_normal, plane_x_axis_unit_vector): 

106 """ """ 

107 

108 class unitVector: 

109 def __init__(self, u, v, w) -> None: 

110 length = math.sqrt(u**2 + v**2 + w**2) 

111 self.u = u / length 

112 self.v = v / length 

113 self.w = w / length 

114 

115 def rotate(self, theta, withRespectTo): 

116 # Rotate with respect to the withRespectTo vector by theta degrees: 

117 # https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle 

118 a = withRespectTo.u 

119 b = withRespectTo.v 

120 c = withRespectTo.w 

121 

122 rotationMatrix = np.array( 

123 [ 

124 [ 

125 math.cos(theta) + a**2 * (1 - math.cos(theta)), 

126 a * b * (1 - math.cos(theta)) - c * math.sin(theta), 

127 a * c * (1 - math.cos(theta)) + b * math.sin(theta), 

128 ], 

129 [ 

130 b * a * (1 - math.cos(theta)) + c * math.sin(theta), 

131 math.cos(theta) + b**2 * (1 - math.cos(theta)), 

132 b * c * (1 - math.cos(theta)) - a * math.sin(theta), 

133 ], 

134 [ 

135 c * a * (1 - math.cos(theta)) - b * math.sin(theta), 

136 c * b * (1 - math.cos(theta)) + a * math.sin(theta), 

137 math.cos(theta) + c**2 * (1 - math.cos(theta)), 

138 ], 

139 ] 

140 ) 

141 vector = np.array([[self.u], [self.v], [self.w]]) 

142 rotatedVector = rotationMatrix @ vector 

143 return unitVector( 

144 rotatedVector[0][0], 

145 rotatedVector[1][0], 

146 rotatedVector[2][0], 

147 ) 

148 

149 def __pow__(self, otherUnitVector): 

150 # Cross product: 

151 u = self.v * otherUnitVector.w - self.w * otherUnitVector.v 

152 v = self.w * otherUnitVector.u - self.u * otherUnitVector.w 

153 w = self.u * otherUnitVector.v - self.v * otherUnitVector.u 

154 return unitVector(u, v, w) 

155 

156 def __mul__(self, otherUnitVector) -> float: 

157 # Dot product: 

158 return ( 

159 self.u * otherUnitVector.u 

160 + self.v * otherUnitVector.v 

161 + self.w * otherUnitVector.w 

162 ) 

163 

164 if len(plane_normal) != 3: 

165 raise ValueError( 

166 "planeNormal for magneticFieldOnCutPlane must be a list of" 

167 " three numbers!" 

168 ) 

169 

170 if len(plane_x_axis_unit_vector) != 3: 

171 raise ValueError( 

172 "planeXAxis for magneticFieldOnCutPlane must be a list of" 

173 " three numbers!" 

174 ) 

175 

176 plane_normal = unitVector(plane_normal[0], plane_normal[1], plane_normal[2]) 

177 plane_x_axis = unitVector( 

178 plane_x_axis_unit_vector[0], 

179 plane_x_axis_unit_vector[1], 

180 plane_x_axis_unit_vector[2], 

181 ) 

182 

183 # Rotate perperndicular vector with respect to the plane's normal vector 

184 # by 90 degrees to find the second perpendicular vector: 

185 plane_y_axis = plane_x_axis.rotate(math.pi / 2, plane_normal) 

186 

187 # Build the transformation matrix to change from the global coordinate 

188 # system to the plane's coordinate system: 

189 transformationMatrix = np.array( 

190 [ 

191 [plane_x_axis.u, plane_x_axis.v, plane_x_axis.w], 

192 [plane_y_axis.u, plane_y_axis.v, plane_y_axis.w], 

193 [plane_normal.u, plane_normal.v, plane_normal.w], 

194 ] 

195 ) 

196 points = self.points.transpose() 

197 new_points = transformationMatrix @ points 

198 new_points = new_points.transpose() 

199 self.points = new_points[:, 0:2] 

200 

201 if self.data_type == "vector": 

202 reshaped_values = self.values 

203 reshaped_values = reshaped_values.reshape( 

204 (len(self.values) * len(self.time_values), 3) 

205 ) 

206 reshaped_values = reshaped_values.transpose() 

207 new_values = transformationMatrix @ reshaped_values 

208 new_values = new_values.transpose() 

209 new_values = new_values[:, 0:2] 

210 self.values = new_values.reshape( 

211 np.shape(self.values)[0], int(np.shape(self.values)[1] / 3 * 2) 

212 )