Coverage for fiqus/post_processors/PostProcessPancake3D.py: 16%

1import os 

2import logging 

3import timeit 

4 

5import gmsh 

6import matplotlib.pyplot as plt 

7import matplotlib.patches as patches 

8import numpy as np 

9from scipy.interpolate import griddata 

10import matplotlib 

11import matplotlib.pyplot as plt 

12 

13from fiqus.mains.MainPancake3D import Base 

14from fiqus.utils.Utils import GmshUtils 

15from fiqus.parsers.ParserGetDPTimeTable import ParserGetDPTimeTable 

16from fiqus.parsers.ParserGetDPOnSection import ParserGetDPOnSection 

17 

18logger = logging.getLogger(__name__) 

19 

20 

21class TimeSeriesData: 

22 def __init__( 

23 self, timeValues, values, quantityName, units, title=None, fileName=None 

24 ): 

25 """ """ 

26 self.t = timeValues 

27 self.values = values 

28 self.quantityName = quantityName 

29 self.units = units 

30 self.title = title 

31 self.fileName = fileName 

32 

33 def plot(self, dir): 

34 plt.rcParams.update({"font.size": 10}) 

35 plt.style.use("_mpl-gallery") 

36 fig, ax = plt.subplots(figsize=(5.95, 3)) 

37 ax.plot(self.t, self.values, marker="o") 

38 ax.set_xlabel(f"Time (s)") 

39 ax.set_ylabel(f"{self.quantityName} ({self.units})") 

40 if self.title: 

41 ax.set_title(self.title) 

42 

43 # Save the plot: 

44 if self.fileName is None: 

45 quantityName = self.quantityName.replace(" ", "_") 

46 filePath = os.path.join(dir, f"{quantityName}_TimeSeriesPlot.pdf") 

47 else: 

48 filePath = os.path.join(dir, f"{self.fileName}_TimeSeriesPlot.pdf") 

49 

50 plt.savefig(filePath, bbox_inches="tight", pad_inches=0) 

51 plt.close("all") 

52 

53 

54class MagneticField: 

55 def __init__( 

56 self, 

57 pointsForMagnitudes, 

58 magnitudes, 

59 pointsForVectors, 

60 vectors, 

61 units, 

62 title=None, 

63 interpolationMethod="linear", 

64 colormap="jet", 

65 planeNormal=None, 

66 ): 

67 self.pointsForMagnitudes = pointsForMagnitudes 

68 self.magnitudes = magnitudes 

69 self.pointsForVectors = pointsForVectors 

70 self.vectors = vectors 

71 self.units = units 

72 self.title = title 

73 self.interpolationMethod = interpolationMethod 

74 self.colormap = colormap 

75 self.planeNormal = planeNormal 

76 

77 def plot(self, dir, streamlines=True): 

78 start_time = timeit.default_timer() 

79 logger.info(f"Plotting {self.title} has been started.") 

80 plt.rcParams.update({"font.size": 10}) 

81 plt.style.use("_mpl-gallery") 

82 

83 points = self.pointsForMagnitudes 

84 magnitudes = self.magnitudes 

85 

86 xArray = points[:, 0] 

87 yArray = points[:, 1] 

88 

89 cropFactorVertical = 0.75 

90 cropFactorHorizontal = 0.66 

91 num_of_grid_points_in_y = 90 

92 

93 # Put parsed values on a grid: 

94 aspectRatio = (xArray.max() - xArray.min()) / (yArray.max() - yArray.min()) 

95 xLinSpace = np.linspace( 

96 xArray.min() * cropFactorHorizontal, 

97 xArray.max() * cropFactorHorizontal, 

98 round(num_of_grid_points_in_y * aspectRatio), 

99 ) 

100 yLinSpace = np.linspace( 

101 yArray.min() * cropFactorVertical, 

102 yArray.max() * cropFactorVertical, 

103 num_of_grid_points_in_y, 

104 ) 

105 grid_x, grid_y = np.meshgrid(xLinSpace, yLinSpace) 

106 grid_z = griddata(points, magnitudes, (grid_x, grid_y), method="linear") 

107 

108 # Then interpolate the grid values: 

109 grid_points = np.array([grid_x.flatten(), grid_y.flatten()]).T 

110 grid_values = grid_z.flatten() 

111 xLinSpace = np.linspace( 

112 xArray.min() * cropFactorHorizontal, 

113 xArray.max() * cropFactorHorizontal, 

114 round(num_of_grid_points_in_y * aspectRatio) * 6, 

115 ) 

116 yLinSpace = np.linspace( 

117 yArray.min() * cropFactorVertical, 

118 yArray.max() * cropFactorVertical, 

119 num_of_grid_points_in_y * 6, 

120 ) 

121 grid_x, grid_y = np.meshgrid(xLinSpace, yLinSpace) 

122 grid_z = griddata( 

123 grid_points, 

124 grid_values, 

125 (grid_x, grid_y), 

126 method=self.interpolationMethod, 

127 ) 

128 

129 fig, ax = plt.subplots(figsize=(5.95, 4.165)) 

130 

131 # cs = ax.contourf( 

132 # grid_x, grid_y, grid_z, cmap=matplotlib.colormaps[self.colormap], levels=np.linspace(0, 0.8, 16) 

133 # ) 

134 cs = ax.imshow( 

135 grid_z, 

136 cmap=matplotlib.colormaps[self.colormap], 

137 extent=[ 

138 xArray.min() * cropFactorHorizontal * 1000, 

139 xArray.max() * cropFactorHorizontal * 1000, 

140 yArray.min() * cropFactorVertical * 1000, 

141 yArray.max() * cropFactorVertical * 1000, 

142 ], 

143 origin="lower", 

144 aspect="auto", 

145 ) 

146 ax.set_aspect("equal") 

147 ax.set_xlabel("x (mm)") 

148 ax.set_ylabel("z (mm)") 

149 

150 fig.colorbar( 

151 cs, 

152 ax=ax, 

153 label=f"Magnitude of the Magnetic Field ({self.units})", 

154 location="bottom", 

155 ) 

156 

157 # # Create a Rectangle patch for each coil: 

158 # rect1 = patches.Rectangle((5.0e-3*1000, -4.042e-3/2*1000), 120.0e-6*60*1000, 4.042e-3*1000, linewidth=2, edgecolor='k', facecolor='none') 

159 # rect2 = patches.Rectangle((-5.0e-3*1000, -4.042e-3/2*1000), -120.0e-6*60*1000, 4.042e-3*1000, linewidth=2, edgecolor='k', facecolor='none') 

160 # # Add the patch to the Axes 

161 # ax.add_patch(rect1) 

162 # ax.add_patch(rect2) 

163 

164 if streamlines: 

165 points = self.pointsForVectors 

166 vectors = self.vectors 

167 uValues = vectors[:, 0] 

168 vValues = vectors[:, 1] 

169 grid_u = griddata( 

170 points, uValues, (grid_x, grid_y), method=self.interpolationMethod 

171 ) 

172 grid_v = griddata( 

173 points, vValues, (grid_x, grid_y), method=self.interpolationMethod 

174 ) 

175 

176 ax.streamplot( 

177 grid_x * 1000, 

178 grid_y * 1000, 

179 grid_u, 

180 grid_v, 

181 color="white", 

182 linewidth=0.8, 

183 density=1.0, 

184 ) 

185 

186 if self.title: 

187 ax.set_title(self.title) 

188 

189 # Save the plot: 

190 time = self.title.split(" ")[-2] 

191 plotFileName = os.path.join( 

192 dir, f"Magnetic_Field_On_A_Cut_Plane_at_{time}_s.pdf" 

193 ) 

194 plt.savefig(plotFileName, bbox_inches="tight", pad_inches=0) 

195 plt.close("all") 

196 end_time = timeit.default_timer() 

197 time_span = end_time - start_time 

198 logger.info(f"Plotting {self.title} has been finished in {time_span:.2f} s.") 

199 

200 

201class Postprocess(Base): 

202 def __init__( 

203 self, 

204 fdm, 

205 geom_folder, 

206 mesh_folder, 

207 solution_folder, 

208 ) -> None: 

209 super().__init__(fdm, geom_folder, mesh_folder, solution_folder) 

210 

211 def plotTimeSeriesPlots(self): 

212 start_time = timeit.default_timer() 

213 logger.info("Plotting time series plots has been started.") 

214 

215 for timeSeriesPlot in self.pp.timeSeriesPlots: 

216 filePath = os.path.join( 

217 self.solution_folder, timeSeriesPlot.fileName + "-TimeTableFormat.csv" 

218 ) 

219 

220 if hasattr(timeSeriesPlot, 'position'): 

221 if hasattr(timeSeriesPlot.position, 'turnNumber'): 

222 title = ( 

223 f"{timeSeriesPlot.quantityProperName} at turn" 

224 f" {timeSeriesPlot.position.turnNumber}" 

225 ) 

226 else: 

227 title = ( 

228 f"{timeSeriesPlot.quantityProperName} at" 

229 f" ({timeSeriesPlot.position.x}, {timeSeriesPlot.position.y}," 

230 f" {timeSeriesPlot.position.z})" 

231 ) 

232 else: 

233 title = f"{timeSeriesPlot.quantityProperName}" 

234 

235 # Parse data: 

236 parser = ParserGetDPTimeTable(filePath) 

237 timeValues = parser.time_values 

238 scalarvalues = parser.get_equivalent_scalar_values() 

239 

240 if parser.data_type == "vector": 

241 title = title + " (Magnitudes of 3D Vectors)" 

242 

243 elif parser.data_type == "tensor": 

244 title = title + " (Von Misses Equivalents of Rank 2 Tensors)" 

245 

246 data = TimeSeriesData( 

247 timeValues=timeValues, 

248 values=scalarvalues, 

249 quantityName=timeSeriesPlot.quantityProperName, 

250 units=timeSeriesPlot.units, 

251 title=title, 

252 ) 

253 data.plot(dir=self.solution_folder) 

254 

255 end_time = timeit.default_timer() 

256 time_span = end_time - start_time 

257 logger.info( 

258 f"Plotting time series plots has been finished in {time_span:.2f} s." 

259 ) 

260 

261 if self.python_postprocess_gui: 

262 gmsh.initialize() 

263 gmsh.option.setNumber("Geometry.Volumes", 0) 

264 gmsh.option.setNumber("Geometry.Surfaces", 0) 

265 gmsh.option.setNumber("Geometry.Curves", 0) 

266 gmsh.option.setNumber("Geometry.Points", 0) 

267 posFiles = [ 

268 fileName 

269 for fileName in os.listdir(self.solution_folder) 

270 if fileName.endswith(".pos") 

271 ] 

272 for posFile in posFiles: 

273 gmsh.open(os.path.join(self.solution_folder, posFile)) 

274 

275 self.gu = GmshUtils(self.solution_folder) 

276 self.gu.launch_interactive_GUI() 

277 gmsh.finalize() 

278 

279 def plotMagneticFieldOnCutPlane(self): 

280 start_time = timeit.default_timer() 

281 logger.info("Parsing magnetic field on the cut plane has been started.") 

282 

283 magnitudeFilePath = os.path.join( 

284 self.solution_folder, "magneticFieldOnCutPlaneMagnitude-DefaultFormat.pos" 

285 ) 

286 vectorFilePath = os.path.join( 

287 self.solution_folder, "magneticFieldOnCutPlaneVector-DefaultFormat.pos" 

288 ) 

289 

290 magnitudeParser = ParserGetDPOnSection( 

291 magnitudeFilePath, data_type="scalar", depth=1 

292 ) 

293 magnitudeParser.project_values_on_a_plane( 

294 plane_normal=self.pp.magneticFieldOnCutPlane.planeNormal, 

295 plane_x_axis_unit_vector=self.pp.magneticFieldOnCutPlane.planeXAxisUnitVector, 

296 ) 

297 

298 if self.pp.magneticFieldOnCutPlane.streamLines: 

299 vectorParser = ParserGetDPOnSection( 

300 vectorFilePath, data_type="vector", depth=1 

301 ) 

302 vectorParser.project_values_on_a_plane( 

303 plane_normal=self.pp.magneticFieldOnCutPlane.planeNormal, 

304 plane_x_axis_unit_vector=self.pp.magneticFieldOnCutPlane.planeXAxisUnitVector, 

305 ) 

306 

307 end_time = timeit.default_timer() 

308 time_span = end_time - start_time 

309 logger.info( 

310 "Parsing magnetic field on the cut plane has been finished in" 

311 f" {time_span:.2f} s." 

312 ) 

313 for timeStep in range(2, len(magnitudeParser.time_values)): 

314 time = magnitudeParser.time_values[timeStep] 

315 if ( 

316 self.pp.magneticFieldOnCutPlane.timesToBePlotted is None 

317 or time in self.pp.magneticFieldOnCutPlane.timesToBePlotted 

318 ): 

319 magnitudePoints = magnitudeParser.points 

320 magnitudes = magnitudeParser.get_values_at_time_step(timeStep) 

321 vectorPoints = vectorParser.points 

322 vectors = vectorParser.get_values_at_time_step(timeStep) 

323 

324 magneticFieldOnCutPlane = MagneticField( 

325 pointsForMagnitudes=magnitudePoints, 

326 magnitudes=magnitudes, 

327 pointsForVectors=vectorPoints, 

328 vectors=vectors, 

329 units="T", 

330 title=( 

331 "Magnetic Field On Plane" 

332 f" {self.pp.magneticFieldOnCutPlane.planeNormal[0]}x + " 

333 f" {self.pp.magneticFieldOnCutPlane.planeNormal[1]}y + " 

334 f" {self.pp.magneticFieldOnCutPlane.planeNormal[2]}z = 0 at" 

335 f" {time:.2f} s" 

336 ), 

337 interpolationMethod=self.pp.magneticFieldOnCutPlane.interpolationMethod, 

338 colormap=self.pp.magneticFieldOnCutPlane.colormap, 

339 planeNormal=self.pp.magneticFieldOnCutPlane.planeNormal, 

340 ) 

341 magneticFieldOnCutPlane.plot( 

342 dir=self.solution_folder, 

343 streamlines=self.pp.magneticFieldOnCutPlane.streamLines, 

344 )