Coverage for fiqus/mesh_generators/MeshMultipole.py: 90%

1import os 

2import gmsh 

3from numpy import square as square 

4from numpy import sqrt as sqrt 

5 

6from fiqus.utils.Utils import GmshUtils 

7from fiqus.utils.Utils import FilesAndFolders as Util 

8from fiqus.data import DataFiQuS as dF 

9from fiqus.data import DataMultipole as dM 

10from fiqus.data.RegionsModelFiQuS import RegionsModel as rM 

11 

12 

13class Mesh: 

14 def __init__(self, data: dF.FDM() = None, sett: dF.FiQuSSettings() = None, mesh_folder: str = None, 

15 verbose: bool = False): 

16 """ 

17 Class to generate mesh 

18 :param data: FiQuS data model 

19 :param sett: settings data model 

20 :param verbose: If True more information is printed in python console. 

21 """ 

22 self.data: dF.FDM() = data 

23 self.set: dF.FiQuSSettings() = sett 

24 self.mesh_folder = mesh_folder 

25 self.verbose: bool = verbose 

26 

27 self.md = dM.MultipoleData() 

28 self.rm = rM() 

29 

30 self.gu = GmshUtils(self.mesh_folder, self.verbose) 

31 self.gu.initialize() 

32 self.occ = gmsh.model.occ 

33 self.mesh = gmsh.model.mesh 

34 

35 self.brep_iron_curves = {1: set(), 2: set(), 3: set(), 4: set()} 

36 self.mesh_parameters = dict.fromkeys(['SJ', 'SICN', 'SIGE', 'Gamma', 'nodes']) 

37 self.geom_folder = os.path.dirname(self.mesh_folder) 

38 self.model_file = f"{os.path.join(self.mesh_folder, self.data.general.magnet_name)}.msh" 

39 

40 self.colors = {'wedges': [86, 180, 233], # sky blue 

41 'half_turns_pos': [213, 94, 0], # vermilion 

42 'half_turns_neg': [255, 136, 42], # light vermilion 

43 'air': [240, 228, 66], # yellow 

44 'air_inf': [220, 208, 46], # dark yellow 

45 # yoke 

46 'BHiron1': [0, 114, 178], # blue 

47 'BHiron2': [0, 158, 115], # bluish green 

48 'BHiron4': [86, 180, 233], # sky blue 

49 # key 

50 'BHiron3': [220, 208, 46], # dark yellow 

51 # [230, 159, 0], # orange 

52 # collar 

53 'BHiron5': [204, 121, 167], # hopbush 

54 'BHiron6': [0, 114, 178], # blue 

55 'BHiron7': [204, 121, 167]} # reddish purple 

56 

57 def ending_step(self, gui: bool = False): 

58 if gui: 

59 self.gu.launch_interactive_GUI() 

60 else: 

61 gmsh.clear() 

62 gmsh.finalize() 

63 

64 def loadAuxiliaryFile(self): 

65 self.md = Util.read_data_from_yaml(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.aux", dM.MultipoleData) 

66 

67 def updateAuxiliaryFile(self): 

68 md2 = Util.read_data_from_yaml(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.aux", dM.MultipoleData) 

69 md2.domains.physical_groups = self.md.domains.physical_groups 

70 Util.write_data_to_yaml(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.aux", md2.dict()) 

71 

72 def saveMeshFile(self): 

73 gmsh.write(self.model_file) 

74 

75 def saveRegionFile(self): 

76 Util.write_data_to_yaml(f"{os.path.join(self.mesh_folder, self.data.general.magnet_name)}.reg", self.rm.dict()) 

77 

78 def getIronCurvesTags(self): 

79 if self.data.magnet.geometry.with_iron_yoke: 

80 for quadrant, qq in self.md.geometries.iron.quadrants.items(): 

81 for area_name, area in qq.areas.items(): 

82 if area.surface: 

83 self.brep_iron_curves[quadrant] |= set(gmsh.model.getAdjacencies(2, area.surface)[1]) 

84 

85 def defineMesh(self): 

86 curves_list = [] 

87 for coil_nr, coil in self.md.geometries.coil.coils.items(): 

88 for pole_nr, pole in coil.poles.items(): 

89 for layer_nr, layer in pole.layers.items(): 

90 for winding_nr, winding in layer.windings.items(): 

91 for block_key, block in winding.blocks.items(): 

92 hts = block.half_turns 

93 for i in range(len(hts.areas) - 1): 

94 curves_list.append(self.occ.copy([(1, hts.lines[str(i + 1) + 'r'])])[0][1]) 

95 curves_list.append(self.occ.copy([(1, hts.lines['end'])])[0][1]) 

96 

97 self.occ.synchronize() 

98 if not self.data.magnet.mesh.default_mesh: 

99 curve_list_iron = [] 

100 if self.data.magnet.geometry.with_iron_yoke: 

101 # if self.from_brep: 

102 for quadrant, qq in self.brep_iron_curves.items(): 

103 for line in qq: 

104 curve_list_iron.append(self.occ.copy([(1, line)])[0][1]) 

105 # else: 

106 # for quadrant, qq in self.md.geometry.iron.quadrants.items(): 

107 # for line_name, line in qq.lines.items(): 

108 # curve_list_iron.append(self.occ.copy([(1, line)])[0][1]) 

109 

110 size_min = self.data.magnet.mesh.mesh_coil.SizeMin 

111 size_max = self.data.magnet.mesh.mesh_coil.SizeMax 

112 dist_min = self.data.magnet.mesh.mesh_coil.DistMin 

113 dist_max = self.data.magnet.mesh.mesh_coil.DistMax # iron max_radius 

114 

115 list_aux_type = "CurvesList" 

116 list_aux = curve_list_iron 

117 size_min_aux = self.data.magnet.mesh.mesh_iron.SizeMin 

118 size_max_aux = self.data.magnet.mesh.mesh_iron.SizeMax 

119 dist_min_aux = self.data.magnet.mesh.mesh_iron.DistMin 

120 dist_max_aux = self.data.magnet.mesh.mesh_iron.DistMax 

121 

122 else: 

123 min_height = 1. 

124 for name, cond in self.set.Model_Data_GS.conductors.items(): 

125 min_height = min(min_height, cond.cable.bare_cable_height_mean) 

126 

127 bore_centers = [self.occ.addPoint(coil.bore_center.x, coil.bore_center.y, 0.) 

128 for coil_nr, coil in self.md.geometries.coil.coils.items()] 

129 point = gmsh.model.getValue(0, self.md.geometries.coil.coils[1].poles[1].layers[1].windings[1].blocks[ 

130 1].half_turns.points['1i'], []) 

131 min_dist = sqrt(square(point[0] - self.md.geometries.coil.coils[1].bore_center.x) + 

132 square(point[1] - self.md.geometries.coil.coils[1].bore_center.y)) 

133 

134 size_min = min_height / 2 

135 size_max = min_height * 20 

136 dist_min = min_height * 2 

137 dist_max = self.md.geometries.iron.max_radius 

138 

139 list_aux_type = "PointsList" 

140 list_aux = bore_centers 

141 size_min_aux = min_height / 2 

142 size_max_aux = min_height * 20 

143 dist_min_aux = min_dist 

144 dist_max_aux = min_dist 

145 

146 self.occ.synchronize() 

147 

148 distance_coil = self.mesh.field.add("Distance") 

149 self.mesh.field.setNumbers(distance_coil, "CurvesList", curves_list) 

150 self.mesh.field.setNumber(distance_coil, "Sampling", 100) 

151 

152 threshold = self.mesh.field.add("Threshold") 

153 self.mesh.field.setNumber(threshold, "InField", distance_coil) 

154 self.mesh.field.setNumber(threshold, "SizeMin", size_min) 

155 self.mesh.field.setNumber(threshold, "SizeMax", size_max) 

156 self.mesh.field.setNumber(threshold, "DistMin", dist_min) 

157 self.mesh.field.setNumber(threshold, "DistMax", dist_max) 

158 

159 distance_aux = self.mesh.field.add("Distance") 

160 self.mesh.field.setNumbers(distance_aux, list_aux_type, list_aux) 

161 self.mesh.field.setNumber(distance_aux, "Sampling", 100) 

162 

163 threshold_aux = self.mesh.field.add("Threshold") 

164 self.mesh.field.setNumber(threshold_aux, "InField", distance_aux) 

165 self.mesh.field.setNumber(threshold_aux, "SizeMin", size_min_aux) 

166 self.mesh.field.setNumber(threshold_aux, "SizeMax", size_max_aux) 

167 self.mesh.field.setNumber(threshold_aux, "DistMin", dist_min_aux) 

168 self.mesh.field.setNumber(threshold_aux, "DistMax", dist_max_aux) 

169 

170 background = self.mesh.field.add("Min") 

171 self.mesh.field.setNumbers(background, "FieldsList", [threshold, threshold_aux] if ( 

172 self.data.magnet.geometry.with_iron_yoke and not self.data.magnet.mesh.default_mesh 

173 or self.data.magnet.mesh.default_mesh) else [threshold]) 

174 self.mesh.field.setAsBackgroundMesh(background) 

175 

176 def fragment(self): 

177 """ 

178 Fragment and group air domains 

179 """ 

180 holes = [] 

181 for group_name, surfaces in self.md.domains.groups_surfaces.items(): 

182 if group_name != 'air_inf': 

183 holes.extend([(2, s) for s in surfaces]) 

184 fragmented = self.occ.fragment([(2, self.md.geometries.air_inf.areas['inner'].surface)], holes)[1] 

185 self.md.domains.groups_surfaces['air'] = [] 

186 existing_domains = [e[0][1] for e in fragmented[1:]] 

187 for e in fragmented[0]: 

188 if e[1] not in existing_domains: 

189 self.md.domains.groups_surfaces['air'].append(e[1]) 

190 

191 self.occ.synchronize() 

192 

193 def createPhysicalGroups(self): 

194 """ 

195 Creates physical groups by grouping the mirrored entities according to the Roxie domains 

196 """ 

197 pg = self.md.domains.physical_groups 

198 for group_name, surfaces in self.md.domains.groups_surfaces.items(): 

199 pg.surfaces[group_name] = gmsh.model.addPhysicalGroup(2, surfaces) 

200 gmsh.model.setPhysicalName(2, pg.surfaces[group_name], group_name) 

201 if group_name[:5] == 'block': 

202 color = self.colors['half_turns_pos'] if group_name[-3:] == 'pos' else self.colors['half_turns_neg'] 

203 elif group_name[:5] == 'wedge': 

204 color = self.colors['wedges'] 

205 else: 

206 color = self.colors[group_name] 

207 gmsh.model.setColor([(2, i) for i in surfaces], color[0], color[1], color[2]) 

208 

209 pg.curves['air_inf'] = gmsh.model.addPhysicalGroup(1, [self.md.geometries.air_inf.lines['outer']]) 

210 gmsh.model.setPhysicalName(1, pg.curves['air_inf'], 'air_inf') 

211 

212 def assignRegionsTags(self): 

213 self.rm.air_far_field.vol.radius_out = float(gmsh.model.getValue(0, gmsh.model.getAdjacencies( 

214 1, self.md.geometries.air_inf.lines['outer'])[1][0], []).max()) 

215 self.rm.air_far_field.vol.radius_in = float(gmsh.model.getValue(0, gmsh.model.getAdjacencies( 

216 1, self.md.geometries.air_inf.lines['inner'])[1][0], []).max()) 

217 

218 self.rm.air.vol.name = "Air" 

219 self.rm.air.vol.number = self.md.domains.physical_groups.surfaces['air'] 

220 self.rm.air_far_field.vol.names = ["AirInf"] 

221 self.rm.air_far_field.vol.numbers = [self.md.domains.physical_groups.surfaces['air_inf']] 

222 self.rm.powered.vol.names = [] 

223 self.rm.powered.vol.numbers = [] 

224 self.rm.iron.vol.names = [] 

225 self.rm.iron.vol.numbers = [] 

226 self.rm.induced.vol.names = [] 

227 self.rm.induced.vol.numbers = [] 

228 for group_name, surface in self.md.domains.physical_groups.surfaces.items(): 

229 if group_name[:5] == 'block': 

230 self.rm.powered.vol.names.append(group_name) 

231 self.rm.powered.vol.numbers.append(surface) 

232 # self.fd.half_turns.cross_sections.append(self.conductor_areas[name]) 

233 elif group_name[:2] == 'BH': 

234 self.rm.iron.vol.names.append(group_name) 

235 self.rm.iron.vol.numbers.append(surface) 

236 elif group_name[:5] == 'wedge': 

237 self.rm.induced.vol.names.append(group_name) 

238 self.rm.induced.vol.numbers.append(surface) 

239 

240 self.rm.air_far_field.surf.name = "Surface_Inf" 

241 self.rm.air_far_field.surf.number = self.md.domains.physical_groups.curves['air_inf'] 

242 

243 def setMeshOptions(self): 

244 """ 

245 Meshes the generated domain 

246 """ 

247 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) 

248 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) 

249 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0) 

250 if not self.data.magnet.mesh.default_mesh: 

251 # gmsh.option.setNumber("Mesh.AngleToleranceFacetOverlap", self.data.one_lab.mesh.AngleToleranceFacetOverlap) 

252 # gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", self.data.one_lab.mesh.MeshSizeFromCurvature) 

253 # gmsh.option.setNumber("Mesh.MeshSizeMin", self.data.one_lab.mesh.MeshSizeMin) 

254 # gmsh.option.setNumber("Mesh.MeshSizeMax", self.data.one_lab.mesh.MeshSizeMax) 

255 gmsh.option.setNumber("Mesh.Algorithm", self.data.magnet.mesh.Algorithm) 

256 gmsh.option.setNumber("Mesh.Optimize", self.data.magnet.mesh.Optimize) 

257 gmsh.option.setNumber("Mesh.ElementOrder", self.data.magnet.mesh.ElementOrder) 

258 else: 

259 gmsh.option.setNumber("Mesh.Algorithm", 6) 

260 gmsh.option.setNumber("Mesh.Optimize", 1) 

261 gmsh.option.setNumber("Mesh.ElementOrder", 2) 

262 

263 def generateMesh(self): 

264 self.mesh.generate(2) 

265 self.mesh.removeDuplicateNodes() 

266 

267 def checkMeshQuality(self): 

268 tags = self.mesh.getElements(2)[1][0] 

269 

270 self.mesh_parameters['SJ'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSJ')) 

271 self.mesh_parameters['SICN'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSICN')) 

272 self.mesh_parameters['SIGE'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSIGE')) 

273 self.mesh_parameters['Gamma'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='gamma')) 

274 self.mesh_parameters['nodes'] = len(self.mesh.getNodes()[0]) 

275 

276 # gmsh.plugin.setNumber("AnalyseMeshQuality", "JacobianDeterminant", 1) 

277 # gmsh.plugin.setNumber("AnalyseMeshQuality", "CreateView", 100) 

278 # test = gmsh.plugin.run("AnalyseMeshQuality") 

279 # test2 = gmsh.view.getModelData(test, test) 

280 

281 # gmsh.logger.getLastError() 

282 # gmsh.logger.get()