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

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 

14import fiqus.data.DataFiQuSPancake3D as dm 

15from fiqus.utils.Utils import GmshUtils 

16from fiqus.parsers.ParserGetDPTimeTable import ParserGetDPTimeTable 

17from fiqus.parsers.ParserGetDPOnSection import ParserGetDPOnSection 

18 

19logger = logging.getLogger(__name__) 

20 

21 

22class TimeSeriesData: 

23 def __init__( 

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

25 ): 

26 """ """ 

27 self.t = timeValues 

28 self.values = values 

29 self.quantityName = quantityName 

30 self.units = units 

31 self.title = title 

32 self.fileName = fileName 

33 

34 def plot(self, dir): 

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

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

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

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

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

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

41 if self.title: 

42 ax.set_title(self.title) 

43 

44 # Save the plot: 

45 if self.fileName is None: 

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

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

48 else: 

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

50 

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

52 plt.close("all") 

53 

54 

55class MagneticField: 

56 def __init__( 

57 self, 

58 pointsForMagnitudes, 

59 magnitudes, 

60 pointsForVectors, 

61 vectors, 

62 units, 

63 title=None, 

64 interpolationMethod="linear", 

65 colormap="jet", 

66 planeNormal=None, 

67 ): 

68 self.pointsForMagnitudes = pointsForMagnitudes 

69 self.magnitudes = magnitudes 

70 self.pointsForVectors = pointsForVectors 

71 self.vectors = vectors 

72 self.units = units 

73 self.title = title 

74 self.interpolationMethod = interpolationMethod 

75 self.colormap = colormap 

76 self.planeNormal = planeNormal 

77 

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

79 start_time = timeit.default_timer() 

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

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

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

83 

84 points = self.pointsForMagnitudes 

85 magnitudes = self.magnitudes 

86 

87 xArray = points[:, 0] 

88 yArray = points[:, 1] 

89 

90 cropFactorVertical = 0.75 

91 cropFactorHorizontal = 0.66 

92 num_of_grid_points_in_y = 90 

93 

94 # Put parsed values on a grid: 

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

96 xLinSpace = np.linspace( 

97 xArray.min() * cropFactorHorizontal, 

98 xArray.max() * cropFactorHorizontal, 

99 round(num_of_grid_points_in_y * aspectRatio), 

100 ) 

101 yLinSpace = np.linspace( 

102 yArray.min() * cropFactorVertical, 

103 yArray.max() * cropFactorVertical, 

104 num_of_grid_points_in_y, 

105 ) 

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

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

108 

109 # Then interpolate the grid values: 

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

111 grid_values = grid_z.flatten() 

112 xLinSpace = np.linspace( 

113 xArray.min() * cropFactorHorizontal, 

114 xArray.max() * cropFactorHorizontal, 

115 round(num_of_grid_points_in_y * aspectRatio) * 6, 

116 ) 

117 yLinSpace = np.linspace( 

118 yArray.min() * cropFactorVertical, 

119 yArray.max() * cropFactorVertical, 

120 num_of_grid_points_in_y * 6, 

121 ) 

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

123 grid_z = griddata( 

124 grid_points, 

125 grid_values, 

126 (grid_x, grid_y), 

127 method=self.interpolationMethod, 

128 ) 

129 

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

131 

132 # cs = ax.contourf( 

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

134 # ) 

135 cs = ax.imshow( 

136 grid_z, 

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

138 extent=[ 

139 xArray.min() * cropFactorHorizontal * 1000, 

140 xArray.max() * cropFactorHorizontal * 1000, 

141 yArray.min() * cropFactorVertical * 1000, 

142 yArray.max() * cropFactorVertical * 1000, 

143 ], 

144 origin="lower", 

145 aspect="auto", 

146 ) 

147 ax.set_aspect("equal") 

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

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

150 

151 fig.colorbar( 

152 cs, 

153 ax=ax, 

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

155 location="bottom", 

156 ) 

157 

158 # # Create a Rectangle patch for each coil: 

159 # 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') 

160 # 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') 

161 # # Add the patch to the Axes 

162 # ax.add_patch(rect1) 

163 # ax.add_patch(rect2) 

164 

165 if streamlines: 

166 points = self.pointsForVectors 

167 vectors = self.vectors 

168 uValues = vectors[:, 0] 

169 vValues = vectors[:, 1] 

170 grid_u = griddata( 

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

172 ) 

173 grid_v = griddata( 

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

175 ) 

176 

177 ax.streamplot( 

178 grid_x * 1000, 

179 grid_y * 1000, 

180 grid_u, 

181 grid_v, 

182 color="white", 

183 linewidth=0.8, 

184 density=1.0, 

185 ) 

186 

187 if self.title: 

188 ax.set_title(self.title) 

189 

190 # Save the plot: 

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

192 plotFileName = os.path.join( 

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

194 ) 

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

196 plt.close("all") 

197 end_time = timeit.default_timer() 

198 time_span = end_time - start_time 

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

200 

201 

202class Postprocess(Base): 

203 def __init__( 

204 self, 

205 fdm, 

206 geom_folder, 

207 mesh_folder, 

208 solution_folder, 

209 ) -> None: 

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

211 

212 def plotTimeSeriesPlots(self): 

213 start_time = timeit.default_timer() 

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

215 

216 for timeSeriesPlot in self.pp.timeSeriesPlots: 

217 filePath = os.path.join( 

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

219 ) 

220 # check if this is Pancake3DPositionInCoordinates or not: 

221 if isinstance( 

222 timeSeriesPlot, dm.Pancake3DPostprocessTimeSeriesPlotPositionRequired 

223 ): 

224 if isinstance( 

225 timeSeriesPlot.position, dm.Pancake3DPositionInCoordinates 

226 ): 

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 = ( 

234 f"{timeSeriesPlot.quantityProperName} at turn" 

235 f" {timeSeriesPlot.position.turnNumber}" 

236 ) 

237 else: 

238 title = f"{timeSeriesPlot.quantityProperName}" 

239 

240 # Parse data: 

241 parser = ParserGetDPTimeTable(filePath) 

242 timeValues = parser.time_values 

243 scalarvalues = parser.get_equivalent_scalar_values() 

244 

245 if parser.data_type == "vector": 

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

247 

248 elif parser.data_type == "tensor": 

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

250 

251 data = TimeSeriesData( 

252 timeValues=timeValues, 

253 values=scalarvalues, 

254 quantityName=timeSeriesPlot.quantityProperName, 

255 units=timeSeriesPlot.units, 

256 title=title, 

257 ) 

258 data.plot(dir=self.solution_folder) 

259 

260 end_time = timeit.default_timer() 

261 time_span = end_time - start_time 

262 logger.info( 

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

264 ) 

265 

266 if self.python_postprocess_gui: 

267 gmsh.initialize() 

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

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

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

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

272 posFiles = [ 

273 fileName 

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

275 if fileName.endswith(".pos") 

276 ] 

277 for posFile in posFiles: 

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

279 

280 self.gu = GmshUtils(self.solution_folder) 

281 self.gu.launch_interactive_GUI() 

282 gmsh.finalize() 

283 

284 def plotMagneticFieldOnCutPlane(self): 

285 start_time = timeit.default_timer() 

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

287 

288 magnitudeFilePath = os.path.join( 

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

290 ) 

291 vectorFilePath = os.path.join( 

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

293 ) 

294 

295 magnitudeParser = ParserGetDPOnSection( 

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

297 ) 

298 magnitudeParser.project_values_on_a_plane( 

299 plane_normal=self.pp.magneticFieldOnCutPlane.planeNormal, 

300 plane_x_axis_unit_vector=self.pp.magneticFieldOnCutPlane.planeXAxisUnitVector, 

301 ) 

302 

303 if self.pp.magneticFieldOnCutPlane.streamLines: 

304 vectorParser = ParserGetDPOnSection( 

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

306 ) 

307 vectorParser.project_values_on_a_plane( 

308 plane_normal=self.pp.magneticFieldOnCutPlane.planeNormal, 

309 plane_x_axis_unit_vector=self.pp.magneticFieldOnCutPlane.planeXAxisUnitVector, 

310 ) 

311 

312 end_time = timeit.default_timer() 

313 time_span = end_time - start_time 

314 logger.info( 

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

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

317 ) 

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

319 time = magnitudeParser.time_values[timeStep] 

320 if ( 

321 self.pp.magneticFieldOnCutPlane.timesToBePlotted is None 

322 or time in self.pp.magneticFieldOnCutPlane.timesToBePlotted 

323 ): 

324 magnitudePoints = magnitudeParser.points 

325 magnitudes = magnitudeParser.get_values_at_time_step(timeStep) 

326 vectorPoints = vectorParser.points 

327 vectors = vectorParser.get_values_at_time_step(timeStep) 

328 

329 magneticFieldOnCutPlane = MagneticField( 

330 pointsForMagnitudes=magnitudePoints, 

331 magnitudes=magnitudes, 

332 pointsForVectors=vectorPoints, 

333 vectors=vectors, 

334 units="T", 

335 title=( 

336 "Magnetic Field On Plane" 

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

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

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

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

341 ), 

342 interpolationMethod=self.pp.magneticFieldOnCutPlane.interpolationMethod, 

343 colormap=self.pp.magneticFieldOnCutPlane.colormap, 

344 planeNormal=self.pp.magneticFieldOnCutPlane.planeNormal, 

345 ) 

346 magneticFieldOnCutPlane.plot( 

347 dir=self.solution_folder, 

348 streamlines=self.pp.magneticFieldOnCutPlane.streamLines, 

349 )