Coverage for fiqus/mesh_generators/MeshConductorAC_Strand_RutherfordCopy.py: 11%

1import os 

2import pickle 

3import logging 

4 

5import gmsh 

6import numpy as np 

7 

8from fiqus.data import RegionsModelFiQuS 

9from fiqus.utils.Utils import GmshUtils, FilesAndFolders 

10from fiqus.data.RegionsModelFiQuS import RegionsModel 

11 

12from fiqus.geom_generators.GeometryConductorAC_Strand import TwistedStrand, Line 

13 

14occ = gmsh.model.occ 

15 

16logger = logging.getLogger('FiQuS') 

17 

18class Mesh: 

19 def __init__(self, fdm, verbose=True): 

20 """ 

21 A base-class used to manage the mesh for the CAC Strand and Rutherford models. 

22 

23 :ivar fdm: The fiqus data model for input parameters. 

24 :vartype fdm: dict 

25 :ivar cacdm: The magnet section from the fdm. 

26 :vartype cacdm: object 

27 :ivar mesh_folder: The path to the folder where the mesh files are stored. 

28 :vartype mesh_folder: str 

29 :ivar magnet_name: The name of the magnet model. 

30 :vartype magnet_name: str 

31 :ivar mesh_file: The path to the .msh file for the mesh. 

32 :vartype mesh_file: str 

33 :ivar regions_file: The path to the .regions file for the mesh. 

34 :vartype regions_file: str 

35 :ivar verbose: If True, the class will print additional information during execution. 

36 :vartype verbose: bool 

37 :ivar gu: An instance of the GmshUtils class for managing the gmsh utility. 

38 :vartype gu: GmshUtils 

39 """ 

40 self.fdm = fdm 

41 self.cacdm = fdm.magnet 

42 self.mesh_folder = os.path.join(os.getcwd()) 

43 self.magnet_name = fdm.general.magnet_name 

44 self.mesh_file = os.path.join(self.mesh_folder, f"{self.magnet_name}.msh") 

45 self.regions_file = os.path.join(self.mesh_folder, f"{self.magnet_name}.regions") 

46 self.verbose = verbose 

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

48 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh) 

49 

50 gmsh.option.setNumber("General.Terminal", verbose) 

51 

52 def save_mesh(self, gui: bool = False): 

53 """ 

54 Saves the mesh to a .msh file. If gui is True, the mesh is also loaded in the gmsh GUI. 

55 """ 

56 gmsh.write(self.mesh_file) 

57 if gui: 

58 self.gu.launch_interactive_GUI() 

59 else: 

60 if gmsh.isInitialized(): 

61 gmsh.clear() 

62 gmsh.finalize() 

63 

64 def load_mesh(self, gui : bool = False): 

65 """ 

66 Loads a previously generated mesh. 

67 """ 

68 gmsh.clear() 

69 gmsh.open(self.mesh_file) 

70 

71 if gui: 

72 self.gu.launch_interactive_GUI() 

73 

74class StrandMesh(Mesh): 

75 """ 

76 A subclass of Mesh that handles mesh generation for the twisted strand cross-section geometries. 

77 

78 :ivar geometry_class: An instance of the TwistedStrand class that stores all the information about the model structure. Initially set to None and should be loaded or assigned before mesh generation methods are called. 

79 :vartype geometry_class: TwistedStrand 

80 

81 :param fdm: The fiqus data model, providing all the necessary data and parameters for mesh generation. 

82 :type fdm: object 

83 :param verbose: A boolean flag that indicates whether additional information should be printed to the console during operation. Defaults to True. 

84 :type verbose: bool, optional 

85 """ 

86 def __init__(self, fdm, verbose=True): 

87 super().__init__(fdm, verbose) 

88 self.geometry_class : TwistedStrand = None # Class from geometry generation step, used to store all information about the model structure. 

89 

90 def filament_field(self): 

91 """ 

92 Generates the filament mesh size field. 

93 

94 The standard mesh size field is linearly interpolated between the mesh size at the filament boundary and the center.  

95 Amplitude dependent meshing is applied if the field is expected to only penetrate a few elements in the filaments, adjusting the mesh size near the filament edges to capture the field penetration region accurately. 

96 

97 :return: The filament mesh size field tag. 

98 :rtype: int 

99 """ 

100 # Define some constants 

101 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius. 

102 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape. 

103 filament_rad = 0 

104 for layer in self.geometry_class.filaments: 

105 for filament in layer: 

106 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag) 

107 for curve in filament.boundary_curves: 

108 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)]) 

109 if r > filament_rad: 

110 filament_rad = r 

111 meshSize_at_filament_boundary = filament_rad * self.cacdm.mesh.filaments.boundary_mesh_size_ratio 

112 meshSize_at_filament_center = filament_rad * self.cacdm.mesh.filaments.center_mesh_size_ratio 

113 filament_edges = sum([filament.boundary_curves for layer in self.geometry_class.filaments for filament in layer], []) 

114 

115 

116 amplitude_dependent_meshing, field_penetration_distance = self.evaluate_amplitude_dependent_meshing().values() 

117 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region 

118 

119 frequency_dependent_meshing, skin_depth = self.evaluate_frequency_dependent_meshing().values() 

120 if frequency_dependent_meshing: 

121 filaments_in_skindepth = [filament for layer in self.geometry_class.filaments for filament in layer if np.linalg.norm(filament.center_point.pos) > self.geometry_class.matrix[-1].rad - 2*skin_depth] 

122 else: 

123 filaments_in_skindepth = [filament for layer in self.geometry_class.filaments for filament in layer] 

124 

125 if amplitude_dependent_meshing: 

126 # The field only penetrates a few elements. We should thus decrease the mesh size close to the filament edges.  

127 # min_mesh_size will be overwritten, while meshSize_at_filament_center will remain.  

128 

129 circular_filament_fields = [] 

130 nodesPerFilament = 2*np.pi*filament_rad / ((field_penetration_distance)/desired_elements_in_field_penetration_region) # We place nodes for meshing at the filament boundary for more efficient meshing (transfinite meshing) 

131 

132 for filament in filaments_in_skindepth: 

133 circular_field = gmsh.model.mesh.field.add("Ball") 

134 gmsh.model.mesh.field.setNumber(circular_field, "Radius", (filament_rad-field_penetration_distance)*0.85 ) 

135 gmsh.model.mesh.field.setNumber(circular_field, "XCenter", filament.center_point.pos[0]) 

136 gmsh.model.mesh.field.setNumber(circular_field, "YCenter", filament.center_point.pos[1]) 

137 gmsh.model.mesh.field.setNumber(circular_field, "Thickness", (filament_rad-field_penetration_distance)*0.15) 

138 gmsh.model.mesh.field.setNumber(circular_field, "VIn", meshSize_at_filament_center) 

139 gmsh.model.mesh.field.setNumber(circular_field, "VOut", (field_penetration_distance)/desired_elements_in_field_penetration_region ) 

140 circular_filament_fields.append(circular_field) 

141 

142 for circle_arc_tag in [CA.tag for CA in filament.boundary_curves]: 

143 gmsh.model.mesh.setTransfiniteCurve(circle_arc_tag, int(nodesPerFilament/len(filament.boundary_curves)), "Progression", 1) 

144 

145 filament_fields = gmsh.model.mesh.field.add("Max") 

146 gmsh.model.mesh.field.setNumbers(filament_fields, "FieldsList", circular_filament_fields) 

147 

148 self.cacdm.mesh.filaments.boundary_mesh_size_ratio = ((field_penetration_distance)/desired_elements_in_field_penetration_region)/filament_rad 

149 

150 

151 else: 

152 # Generalized to arbitrary filament shapes (not only circular) 

153 f_dist_boundary = gmsh.model.mesh.field.add("Distance") 

154 gmsh.model.mesh.field.setNumbers(f_dist_boundary, "CurvesList", [edge.tag for edge in filament_edges]) 

155 

156 mass_centers = [] 

157 for layer in self.geometry_class.filaments: 

158 for filament in layer: 

159 x,y,z = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag) 

160 mass_center_tag = gmsh.model.occ.addPoint(x, y, z) 

161 mass_centers.append(mass_center_tag) 

162 gmsh.model.occ.synchronize() 

163 

164 f_dist_center = gmsh.model.mesh.field.add("Distance") 

165 gmsh.model.mesh.field.setNumbers(f_dist_center, "PointsList", mass_centers) 

166 

167 # Linearly interpolate between the mesh size at the filament center and the mesh size at the filament boundary. 

168 filament_fields = gmsh.model.mesh.field.add("MathEval") 

169 gmsh.model.mesh.field.setString( 

170 filament_fields, 

171 "F", 

172 f"( {meshSize_at_filament_boundary} * F{f_dist_center} + {meshSize_at_filament_center} * F{f_dist_boundary} ) / ( F{f_dist_center} + F{f_dist_boundary} )" 

173 ) 

174 

175 

176 # Restrict the filament fields to the filaments 

177 filament_mesh_field = gmsh.model.mesh.field.add("Restrict") 

178 gmsh.model.mesh.field.setNumber(filament_mesh_field, "InField", filament_fields) 

179 

180 

181 if frequency_dependent_meshing: 

182 gmsh.model.mesh.field.setNumbers(filament_mesh_field, "SurfacesList", [filament.surface_tag for filament in filaments_in_skindepth]) 

183 else: 

184 gmsh.model.mesh.field.setNumbers(filament_mesh_field, "SurfacesList", [filament.surface_tag for layer in self.geometry_class.filaments for filament in layer] + [hole.surface_tag for layer in self.geometry_class.filament_holes for hole in layer]) 

185 

186 return filament_mesh_field 

187 

188 def matrix_threshold_to_filament_field(self, filaments): 

189 """ 

190 Adjusts the matrix mesh size based on the distance from the filaments to improve mesh quality. 

191 

192 This method creates a mesh size field that gradually refines the mesh in the matrix as it approaches the filaments.  

193 Filaments typically have a finer mesh than the matrix, and a large gradient in mesh size at their boundary can degrade the mesh quality.  

194 The method interpolates the mesh size from a minimum value at the filament boundary to a maximum value at a specified distance from the filaments, based on configuration parameters in the FDM. 

195 

196 :param filaments: A list of filament surface objects for which the boundary mesh size adjustments are to be made. 

197 :type filaments: list[list[object]] 

198 

199 :return: The tag of the created mesh size field, which is restricted to the matrix domain. 

200 :rtype: int 

201 """ 

202 matrix = self.geometry_class.matrix 

203 matrix_mesh_options = self.cacdm.mesh.matrix 

204 

205 filament_edges = sum([[CA.tag for CA in filament.boundary_curves] for filament in filaments], []) # List of tags of all filament boundaries 

206 # Distance to filament boundaries field: 

207 distance_to_filament_field = gmsh.model.mesh.field.add("Distance") 

208 gmsh.model.mesh.field.setNumbers(distance_to_filament_field, "CurvesList", filament_edges) 

209 

210 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius. 

211 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape. 

212 filament_rad = 0 

213 for layer in self.geometry_class.filaments: 

214 for filament in layer: 

215 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag) 

216 for curve in filament.boundary_curves: 

217 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)]) 

218 if r > filament_rad: 

219 filament_rad = r 

220 

221 filament_boundary_meshSize = filament_rad * self.cacdm.mesh.filaments.boundary_mesh_size_ratio 

222 

223 # Linarily interpolate the mesh size from the filament boundaries to a specified distance from filaments. 

224 threshold_to_filament_field = gmsh.model.mesh.field.add("Threshold") 

225 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "InField", distance_to_filament_field) 

226 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "SizeMin", filament_boundary_meshSize) 

227 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "SizeMax", filament_rad*matrix_mesh_options.mesh_size_matrix_ratio_inner) 

228 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "DistMin", 0) 

229 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "DistMax", filament_rad*matrix_mesh_options.interpolation_distance_from_filaments_ratio) 

230 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "StopAtDistMax", 1) 

231 # gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "Sigmoid", 0) 

232 

233 meshSize_by_filament_field = gmsh.model.mesh.field.add("Restrict") # Restrict the mesh to the matrix 

234 gmsh.model.mesh.field.setNumber(meshSize_by_filament_field, "InField", threshold_to_filament_field) 

235 gmsh.model.mesh.field.setNumbers(meshSize_by_filament_field, "SurfacesList", [partition.surface_tag for partition in matrix]) 

236 # matrix_mesh_size_fields.append(meshSize_by_filament_field) 

237 return meshSize_by_filament_field 

238 

239 def evaluate_amplitude_dependent_meshing(self): 

240 """ 

241 Evaluates the necessity for amplitude-dependent meshing based on the applied magnetic field's amplitude. 

242 

243 This method assesses whether the amplitude of the applied magnetic field is low enough that it only partially penetrates the filaments.  

244 In such cases, a finer mesh is required in the field penetration region to accurately capture field variations, while a coarser mesh can be used in the filament center where the field does not penetrate.  

245 It calculates the field penetration distance and determines whether an adjustment in the mesh is necessary to ensure accurate simulation results. 

246 

247 :return: A dictionary containing: 

248 - "amplitude_dependent_meshing" (bool): Indicates whether amplitude-dependent meshing is needed based on the magnetic field's penetration. 

249 - "field_penetration_distance" (float): The estimated distance the magnetic field penetrates into the filaments, used to adjust the mesh size accordingly. 

250 :rtype: Dict[str, Union[bool, float]] 

251 """ 

252 if self.cacdm.mesh.filaments.amplitude_dependent_scaling and self.cacdm.solve.source_parameters.sine.current_amplitude == 0 and not self.cacdm.geometry.io_settings.load.load_from_yaml: 

253 filament_radius = self.geometry_class.filaments[0][0].rad 

254 min_meshSize = filament_radius * self.cacdm.mesh.filaments.boundary_mesh_size_ratio 

255 meshSize_at_filament_center = filament_radius * self.cacdm.mesh.filaments.center_mesh_size_ratio 

256 

257 applied_field_amplitude = self.cacdm.solve.source_parameters.sine.field_amplitude 

258 applied_field_amplitude = max(applied_field_amplitude, 0.01) # Limit the amplitude to 0.01 T to avoid excessively fine mesh (Could be removed?). 

259 jc = 3e9 #self.cacdm.solve.material_properties.superconducting.jc.constant 

260 mu0 = np.pi*4e-7 

261 

262 field_penetration_distance = applied_field_amplitude/(mu0*jc) * self.cacdm.mesh.filaments.field_penetration_depth_scaling_factor # Estimate the field penetration distance 

263 

264 # The mesh size in the filaments can be written as mesh_size(x) = a*x+b, with x being the distance from the filament boundary and a,b defined as: 

265 a = (meshSize_at_filament_center - min_meshSize) / filament_radius 

266 b = min_meshSize 

267 # The number of elements in the range x:[x0, x1] is approximated by the function: 

268 number_of_elements_in_range = lambda a, b, x0, x1 : np.log(a*x1+b)/a - np.log(a*x0+b)/a if a != 0 else (x1-x0)/b 

269 

270 number_of_elements_penetrated_by_field = number_of_elements_in_range(a, b, 0, field_penetration_distance) 

271 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region 

272 

273 if (number_of_elements_penetrated_by_field < desired_elements_in_field_penetration_region) and field_penetration_distance < filament_radius: 

274 logger.info(f"The magnetic field is expected to penetrate over approximately {round(number_of_elements_penetrated_by_field, 3)} elements in the filaments. To ensure that the field penetrates over at least {desired_elements_in_field_penetration_region} elements, the mesh in the filaments has been adjusted.") 

275 return dict(amplitude_dependent_meshing = True, field_penetration_distance = field_penetration_distance) 

276 

277 return dict(amplitude_dependent_meshing = False, field_penetration_distance = 1e10) 

278 

279 def evaluate_frequency_dependent_meshing(self): 

280 """ 

281 Evaluates the need for frequency-dependent meshing based on the skin depth due to the applied magnetic field's frequency. 

282 

283 This method determines if the frequency of the applied magnetic field is high enough to confine the currents to a thin layer near the surface of the strand.  

284 A fine mesh is required in this region to capture current variations accurately, while a coarser mesh can suffice in the strand center.  

285 The method calculates the skin depth using the field frequency and the matrix resistivity, then approximates the number of elements within this depth.  

286 If the calculated number is less than desired, it indicates that frequency-dependent meshing adjustments are needed. 

287 

288 :return: A dictionary containing: 

289 - "frequency_dependent_meshing" (bool): Indicates whether frequency-dependent meshing adjustments are necessary. 

290 - "skin_depth" (float): The calculated skin depth, which informs how the mesh should be adjusted to accurately model current variations. 

291 :rtype: Dict[str, Union[bool, float]] 

292 """ 

293 if self.cacdm.mesh.matrix.rate_dependent_scaling_matrix and not self.cacdm.geometry.io_settings.load.load_from_yaml: 

294 matrix_mesh_options = self.cacdm.mesh.matrix 

295 

296 filament_radius = self.geometry_class.filaments[0][0].rad 

297 

298 mesh_size_matrix_inner = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_inner 

299 mesh_size_matrix_middle = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_middle 

300 mesh_size_matrix_outer = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_outer 

301 

302 matrix = self.geometry_class.matrix 

303 

304 rho_matrix = 1.81e-10 #self.cacdm.solve.material_properties.matrix.resistivity.constant 

305 mu0 = np.pi*4e-7 

306 skin_depth = np.sqrt( rho_matrix / (np.pi*mu0 * self.cacdm.solve.source_parameters.sine.frequency ) ) * self.cacdm.mesh.matrix.skindepth_scaling_factor 

307 # Function which approximates number of elements in a range given a linearly interpolated mesh size. 

308 number_of_elements_in_range = lambda a, b, x0, x1 : np.log(a*x1+b)/a - np.log(a*x0+b)/a if a != 0 else (x1-x0)/b 

309 

310 # If the skindepth is smaller than the outer matrix, we only evaluate the outer matrix field to approximate number of elements in skindepth.  

311 if skin_depth <= matrix[-1].rad -matrix[-2].rad: 

312 a = (mesh_size_matrix_outer - mesh_size_matrix_middle) / (matrix[-1].rad - matrix[-2].rad) 

313 b = mesh_size_matrix_middle 

314 elements_in_skindepth = number_of_elements_in_range(a, b, matrix[-1].rad -matrix[-2].rad - skin_depth, matrix[-1].rad-matrix[-2].rad) 

315 

316 elif skin_depth > matrix[-1].rad -matrix[-2].rad and skin_depth < matrix[-1].rad: 

317 # If the skindepth is larger than the outer matrix, we evaluate the inner and outer matrix field to approximate number of elements in skindepth.  

318 a1 = (mesh_size_matrix_middle - mesh_size_matrix_inner) / (matrix[-2].rad) 

319 b1 = mesh_size_matrix_middle 

320 a2 = (mesh_size_matrix_outer - mesh_size_matrix_middle) / (matrix[-1].rad - matrix[-2].rad) 

321 b2 = mesh_size_matrix_middle 

322 elements_in_skindepth = number_of_elements_in_range(a1, b1, matrix[-1].rad-skin_depth, matrix[-2].rad) + number_of_elements_in_range(a2, b2, 0, matrix[-1].rad-matrix[-2].rad) 

323 

324 else: # If the skindepth is greater than the matrix radius we do not want to act, and set the number of elements in skindepth to a high number. 

325 elements_in_skindepth = 1e10 

326 

327 desired_elements_in_skindepth = self.cacdm.mesh.matrix.desired_elements_in_skindepth 

328 

329 if elements_in_skindepth < desired_elements_in_skindepth: 

330 logger.info(f"The skin depth of the matrix is expected to contain only {round(elements_in_skindepth, 3)} elements. The mesh in the matrix has been adjusted to ensure at least {desired_elements_in_skindepth} elements in the skindepth.") 

331 return dict(frequency_dependent_meshing = True, skin_depth = skin_depth) 

332 

333 return dict(frequency_dependent_meshing = False, skin_depth = 1e10) 

334 

335 def matrix_field(self): 

336 """ 

337 Creates a mesh size field for the matrix as a combination of three fields: 

338 - A distance-based interpolation between the matrix center and partition boundaries. The mesh size is linearly interpolated from the center to the boundary of the second matrix partition, and then to the outer matrix boundary. 

339 - Adjustments near the filaments to ensure a smooth transition from the filament mesh to the matrix mesh. 

340 - Frequency-dependent adjustments to capture skin depth effects. 

341 

342 :return: The tag of the composite mesh size field created for the matrix. 

343 :rtype: int 

344 """ 

345 matrix_mesh_options = self.cacdm.mesh.matrix 

346 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius. 

347 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape. 

348 filament_rad = 0 

349 for layer in self.geometry_class.filaments: 

350 for filament in layer: 

351 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag) 

352 for curve in filament.boundary_curves: 

353 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)]) 

354 if r > filament_rad: 

355 filament_rad = r 

356 

357 mesh_size_matrix_inner = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_inner 

358 mesh_size_matrix_middle = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_middle 

359 mesh_size_matrix_outer = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_outer 

360 

361 filament_edge_points = [] 

362 for layer in self.geometry_class.filaments: 

363 for filament in layer: 

364 for curve in filament.boundary_curves: 

365 filament_edge_points.append(curve.P1) 

366 filament_edge_points.append(curve.P2) 

367 filament_outer_point = max(filament_edge_points, key=lambda p: np.linalg.norm(p.pos)) 

368 # Distance from the outermost filament edge to the center of the matrix 

369 # This is done instead of referencing the middle matrix radius, to be able to mesh arbitrary geometries. 

370 filament_outer_point_dist = np.linalg.norm(filament_outer_point.pos) 

371 

372 matrix = self.geometry_class.matrix 

373 matrix_outer_point_dist = max([np.linalg.norm(p.pos) for curve in matrix[-1].boundary_curves for p in [curve.P1, curve.P2]]) # Distance from the outermost matrix edge to the center of the matrix 

374 

375 matrix_mesh_size_fields = [] 

376 

377 # 'Distance from matrix center'- field: 

378 dist_from_center_field = gmsh.model.mesh.field.add("MathEval") 

379 gmsh.model.mesh.field.setString(dist_from_center_field, "F", "sqrt(x^2 + y^2)") 

380 

381 # Mesh size-field in the middle partition of the matrix. The mesh size is large at the center and decreases linearily towards the edge of the middle matrix. 

382 interpField_innerToMiddle = gmsh.model.mesh.field.add("Threshold") 

383 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "InField", dist_from_center_field) 

384 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "SizeMin", mesh_size_matrix_inner) 

385 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "SizeMax", mesh_size_matrix_middle) 

386 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "DistMin", 0) 

387 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "DistMax", filament_outer_point_dist) 

388 

389 # Restrict the middle field to the inner and middle matrix (only middle matrix in the case where we have a center filament and thus no inner matrix). 

390 middle_matrix_field = gmsh.model.mesh.field.add("Restrict") 

391 gmsh.model.mesh.field.setNumber(middle_matrix_field, "InField", interpField_innerToMiddle) 

392 gmsh.model.mesh.field.setNumbers(middle_matrix_field, "SurfacesList", [partition.surface_tag for partition in matrix[:-1]]) 

393 matrix_mesh_size_fields.append(middle_matrix_field) 

394 

395 # Mesh size-field in the outer partition of the matrix. The mesh size is small at the edge of the middle matrix and increases towards the outer edge. 

396 interpField_middleToOuter = gmsh.model.mesh.field.add("Threshold") 

397 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "InField", dist_from_center_field) 

398 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "SizeMin", mesh_size_matrix_middle) 

399 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "SizeMax", mesh_size_matrix_outer) 

400 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "DistMin", filament_outer_point_dist) 

401 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "DistMax", matrix_outer_point_dist) 

402 

403 # Restric the outer field to the outer matrix partition: 

404 outer_matrix_field = gmsh.model.mesh.field.add("Restrict") 

405 gmsh.model.mesh.field.setNumber(outer_matrix_field, "InField", interpField_middleToOuter) 

406 gmsh.model.mesh.field.setNumbers(outer_matrix_field, "SurfacesList", [matrix[-1].surface_tag]) 

407 matrix_mesh_size_fields.append(outer_matrix_field) 

408 

409 

410 # The last field adjusts the matrix mesh size with respect to the distance from the filaments. The reason is that the filaments typically has a finer mesh than the matrix and  

411 # the large mesh size gradient at the boundary between the matrix and the filaments can lead to a bad quality mesh. We thus want to gradually refine the mesh in the matrix close to the filaments.  

412 # This is done by interpolating the mesh size based on the distance to the filament boundaries. 

413 meshSize_by_filament_field = self.matrix_threshold_to_filament_field([filament for layer in self.geometry_class.filaments for filament in layer]) 

414 matrix_mesh_size_fields.append(meshSize_by_filament_field) 

415 

416 frequency_dependent_meshing, skin_depth = self.evaluate_frequency_dependent_meshing().values() 

417 

418 if frequency_dependent_meshing: 

419 # If the skindepth only consists of a few elements, we add an additional mesh size field in the matrix. This field is large outside of the skindepth, but 

420 # in the skindepth it is defined as the constant field which gives the skindepth the 'desired number of elements'. In between the two regions we interpolate between the two values. 

421 # logger.info(f"The skin depth of the matrix is expected to contain only {round(elements_in_skindepth, 3)} elements. The mesh in the matrix has been adjusted to ensure at least {desired_elements_in_skindepth} elements in the skindepth.") 

422 

423 circular_field = gmsh.model.mesh.field.add("Ball") 

424 gmsh.model.mesh.field.setNumber(circular_field, "Radius", (matrix[-1].rad-skin_depth)*0.85 ) 

425 # gmsh.model.mesh.field.setNumber(circular_field, "Thickness", (matrix[-1].rad-skin_depth)*0.15) 

426 gmsh.model.mesh.field.setNumber(circular_field, "VIn", mesh_size_matrix_inner) 

427 gmsh.model.mesh.field.setNumber(circular_field, "VOut", (skin_depth)/self.cacdm.mesh.matrix.desired_elements_in_skindepth ) 

428 

429 

430 

431 # -- Restrict the mesh size field to the matrix and the filaments not in the skindepth region -- # 

432 filaments_not_in_skindepth = [filament for layer in self.geometry_class.filaments for filament in layer if np.linalg.norm(filament.center_point.pos) < self.geometry_class.matrix[-1].rad - 2*skin_depth] 

433 

434 circular_field_restricted = gmsh.model.mesh.field.add("Restrict") 

435 gmsh.model.mesh.field.setNumber(circular_field_restricted, "InField", circular_field) 

436 gmsh.model.mesh.field.setNumbers(circular_field_restricted, "SurfacesList", [partition.surface_tag for partition in matrix] + [filament.surface_tag for filament in filaments_not_in_skindepth] ) 

437 # -- -- # 

438 

439 self.cacdm.mesh.matrix.mesh_size_matrix_ratio_outer = skin_depth/self.cacdm.mesh.matrix.desired_elements_in_skindepth * 1/filament_rad 

440 

441 filaments_in_skindepth = [filament for layer in self.geometry_class.filaments for filament in layer if np.linalg.norm(filament.center_point.pos) > self.geometry_class.matrix[-1].rad - 2*skin_depth] 

442 meshSize_by_filament_field = self.matrix_threshold_to_filament_field(filaments_in_skindepth) # Redefine the interpolation field to apply only to the filaments in the skindepth region. 

443 

444 matrix_meshSize_field = gmsh.model.mesh.field.add("Min") 

445 gmsh.model.mesh.field.setNumbers(matrix_meshSize_field, "FieldsList", [circular_field_restricted, meshSize_by_filament_field]) 

446 return matrix_meshSize_field 

447 

448 

449 

450 

451 # Combine all the matrix mesh size fields by selecting the minimum value of all fields. 

452 matrix_meshSize_field = gmsh.model.mesh.field.add("Min") 

453 gmsh.model.mesh.field.setNumbers(matrix_meshSize_field, "FieldsList", matrix_mesh_size_fields) 

454 

455 return matrix_meshSize_field 

456 

457 def air_field(self): 

458 """ 

459 Generates a mesh size field for the air domain surrounding the model. 

460 The field linearly interpolates between the mesh size at the outer boundary of the matrix and the air boundary. 

461 

462 :return: The tag of the created mesh size field for the air domain. 

463 :rtype: int 

464 """ 

465 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius. 

466 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape. 

467 filament_rad = 0 

468 for layer in self.geometry_class.filaments: 

469 for filament in layer: 

470 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag) 

471 for curve in filament.boundary_curves: 

472 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)]) 

473 if r > filament_rad: 

474 filament_rad = r 

475 

476 mesh_size_outer_matrix = filament_rad * self.cacdm.mesh.matrix.mesh_size_matrix_ratio_outer 

477 mesh_size_air_boundary = filament_rad * self.cacdm.mesh.air.max_mesh_size_ratio 

478 

479 if self.geometry_class.air_composition: 

480 air_outer_boundary_curves = self.geometry_class.air_composition.boundary_curves 

481 air_inner_boundary_curves = sum(self.geometry_class.air_composition.inner_boundary_curves, []) 

482 else: 

483 air_outer_boundary_curves = self.geometry_class.air[0].boundary_curves 

484 air_inner_boundary_curves = sum(self.geometry_class.air[0].inner_boundary_curves, []) 

485 

486 dist_from_outer_air_boundary_field = gmsh.model.mesh.field.add("Distance") 

487 gmsh.model.mesh.field.setNumbers(dist_from_outer_air_boundary_field, "CurvesList", [curve.tag for curve in air_outer_boundary_curves]) 

488 

489 dist_from_matrix_boundary = gmsh.model.mesh.field.add("Distance") 

490 gmsh.model.mesh.field.setNumbers(dist_from_matrix_boundary, "CurvesList", [curve.tag for curve in air_inner_boundary_curves]) 

491 

492 air_field = gmsh.model.mesh.field.add("MathEval") 

493 gmsh.model.mesh.field.setString( 

494 air_field, 

495 "F", 

496 f"( {mesh_size_air_boundary} * F{dist_from_matrix_boundary} + {mesh_size_outer_matrix} * F{dist_from_outer_air_boundary_field} ) / ( F{dist_from_matrix_boundary} + F{dist_from_outer_air_boundary_field} )" 

497 ) 

498 return air_field 

499 

500 def load_geometry_class(self, geom_folder): 

501 """ 

502 Loads the TwistedStrand geometry class from a .pkl file saved within the specified geometry folder. 

503 The geometry class contains all the information about the strand geometry, facilitating mesh generation. 

504 

505 :param geom_folder: The path to the folder where the geometry class .pkl file is saved. 

506 :type geom_folder: str 

507 

508 :return: An instance of the TwistedStrand geometry class. 

509 :rtype: object 

510 """ 

511 geom_save_file = os.path.join(geom_folder, f'{self.magnet_name}.pkl') 

512 

513 with open(geom_save_file, "rb") as geom_save_file: # Unnecessary to return geom instead of setting self.geometry_class 

514 geom = pickle.load(geom_save_file) 

515 return geom 

516 

517 def generate_mesh(self, geom_folder): 

518 """ 

519 Generates the mesh for the entire model geometry based on combined mesh size fields. 

520 

521 The method generates mesh size fields for the filament, matrix, and air domains separately and combines these fields to define a unified mesh size field for the full geometry.  

522 The combined field is applied as the background mesh size field. 

523 

524 :param geom_folder: The path to the folder containing the saved geometry class (.pkl file). This folder is used to load the geometry details required for mesh generation. 

525 :type geom_folder: str 

526 """ 

527 

528 # So far: Adds physical groups and generates an automatic mesh. To be further developed. 

529 self.geometry_class : TwistedStrand = self.load_geometry_class(geom_folder) 

530 self.geometry_class.update_tags() 

531 self.geometry_class.add_physical_groups() 

532 

533 total_filament_boundary_distance_field = self.filament_field() 

534 matrixField = self.matrix_field() 

535 airField = self.air_field() 

536 

537 total_meshSize_field = gmsh.model.mesh.field.add("Min") 

538 gmsh.model.mesh.field.setNumbers(total_meshSize_field, "FieldsList", [total_filament_boundary_distance_field, matrixField, airField]) 

539 

540 gmsh.model.mesh.field.setAsBackgroundMesh(total_meshSize_field) 

541 

542 gmsh.option.setNumber("Mesh.MeshSizeFactor", self.cacdm.mesh.scaling_global) 

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

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

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

546 

547 # enforce center nodes 

548 if self.fdm.magnet.mesh.matrix.force_center_symmetry: 

549 dist = self.geometry_class.filaments[0][0].center_point.pos[0]*5e-2 # 5% of inner filament distance seems sufficient 

550 center = self.geometry_class.matrix[0].center_point.pos 

551 d_tag = gmsh.model.occ.add_point(center[0],center[1]-dist, center[2]) 

552 u_tag = gmsh.model.occ.add_point(center[0],center[1]+dist, center[2]) 

553 l_tag = gmsh.model.occ.add_point(center[0]-dist,center[1], center[2]) 

554 r_tag = gmsh.model.occ.add_point(center[0]+dist,center[1], center[2]) 

555 outDimTags, _ = gmsh.model.occ.fragment([(2, self.geometry_class.matrix[0].surface_tag)],[(0, gmsh.model.occ.add_point(center[0],center[1],center[2])),(0,d_tag),(0,u_tag),(0,l_tag),(0,r_tag)]) #  

556 gmsh.model.occ.synchronize() 

557 # update the physical surface 

558 gmsh.model.removePhysicalGroups([(2, self.geometry_class.matrix[0].physical_surface_tag)]) 

559 self.geometry_class.matrix[0].physical_surface_tag = gmsh.model.add_physical_group(2, [outDimTags[0][1]], tag=self.geometry_class.matrix[0].physical_surface_tag, name='Surface: Matrix partition 0') 

560 

561 gmsh.model.mesh.generate(2) 

562 

563 def generate_cuts(self): 

564 """ 

565 Computes the cuts for imposing global quantities. 

566 """ 

567 # Computes cuts, as needed for imposing global quantities 

568 if self.verbose: 

569 logger.info('(Co)homology computation') 

570 

571 # a) Cut for the total current intensity (through OmegaCC) 

572 if self.fdm.magnet.geometry.type == 'periodic_square': 

573 # only enforce periodic mesh - cut already defined in composite geometry 

574 self.set_periodic_mesh() 

575 else: 

576 # cohomology cut for single air domain 

577 air_region = self.geometry_class.air[0] 

578 gmsh.model.mesh.addHomologyRequest("Homology", domainTags=[air_region.physical_inner_boundary_tags[0]], dims=[1]) 

579 if self.fdm.magnet.geometry.type == 'strand_only': 

580 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[air_region.physical_surface_tag], dims=[1]) 

581 elif self.fdm.magnet.geometry.type == 'coil': 

582 # add semicircle arc as sudomain to force homology cut onto semicircle diameter 

583 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[air_region.physical_surface_tag], subdomainTags=[air_region.physical_cohomology_subdomain], dims=[1]) 

584 

585 cuts = gmsh.model.mesh.computeHomology() 

586 gmsh.model.mesh.clearHomologyRequests() 

587 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains", str(cuts[1][1])) 

588 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains2", str(cuts[0][1])) 

589 gmsh.plugin.run("HomologyPostProcessing") 

590 air_region.strand_bnd_cut_tag = cuts[1][1]+1 # The cut tag for the strand boundary  

591 

592 # b) Cuts for the individual filaments 

593 for layer in self.geometry_class.filaments: 

594 for filament in layer: 

595 filBnd = filament.physical_boundary_tag 

596 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[filBnd],dims=[1]) 

597 cuts = gmsh.model.mesh.computeHomology() 

598 

599 filaments = sum(self.geometry_class.filaments, []) 

600 for i in range(len(cuts)): 

601 filaments[i].cut_tag = cuts[i][1] 

602 

603 def set_periodic_mesh(self): 

604 """ 

605 Enforces symmetric air mesh by mirroring the node positions from bottom to top and left to right. 

606 """ 

607 air = self.geometry_class.air 

608 # link left and right 

609 dx, dy, dz = air[0].cut_curves[1].P2.pos - air[2].cut_curves[0].P2.pos 

610 gmsh.model.mesh.set_periodic(1, [self.geometry_class.air_composition.physical_left_boundary_tag], [self.geometry_class.air_composition.physical_right_boundary_tag], [1,0,0, dx, 0,1,0, dy, 0,0,1, dz, 0,0,0,1]) 

611 # link bottom and top 

612 dx, dy, dz = air[0].cut_curves[0].P2.pos - air[3].cut_curves[1].P2.pos 

613 gmsh.model.mesh.set_periodic(1, [self.geometry_class.air_composition.physical_bottom_boundary_tag], [self.geometry_class.air_composition.physical_top_boundary_tag], [1,0,0, dx, 0,1,0, dy, 0,0,1, dz, 0,0,0,1]) 

614 

615 def generate_regions_file(self): 

616 """ 

617 Generates a .regions file for the GetDP solver, containing all necessary information about the model. 

618 The regions model contains data about physical surfaces, boundaries, points and cuts and is stored in the mesh folder. 

619 """ 

620 regions_model = RegionsModel() 

621 

622 """ -- Initialize data -- """ 

623 regions_model.powered["Filaments"] = RegionsModelFiQuS.Powered() 

624 regions_model.powered["Filaments"].vol.numbers = [] # Filament physical surfaces tags 

625 regions_model.powered["Filaments"].vol.names = [] # Filament physical surfaces names 

626 regions_model.powered["Filaments"].surf.numbers = [] # Filament physical boundaries tags 

627 regions_model.powered["Filaments"].surf.names = [] # Filament physical boundaries names 

628 regions_model.powered["Filaments"].cochain.numbers = [] # Filament boundary cut tags 

629 regions_model.powered["Filaments"].curve.names = [] # Stores physical points at filament boundary (to fix phi=0) 

630 regions_model.powered["Filaments"].curve.numbers = [] # Stores physical points at filament boundary (to fix phi=0) 

631 

632 regions_model.powered["Filaments"].surf_in.numbers = [] # Filament holes physical tags 

633 regions_model.powered["Filaments"].surf_in.names = [] # Filament holes physical names 

634