Coverage for fiqus/mesh_generators/MeshConductorAC_Strand.py: 89%

1import os 

2import pickle 

3 

4import gmsh 

5import numpy as np 

6 

7from fiqus.data import RegionsModelFiQuS 

8from fiqus.utils.Utils import GmshUtils, FilesAndFolders 

9from fiqus.data.RegionsModelFiQuS import RegionsModel 

10 

11from fiqus.geom_generators.GeometryConductorAC_Strand import TwistedStrand, Line 

12 

13occ = gmsh.model.occ 

14 

15 

16class Mesh: 

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

18 """ 

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

20 

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

22 :vartype fdm: dict 

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

24 :vartype cacdm: object 

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

26 :vartype mesh_folder: str 

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

28 :vartype magnet_name: str 

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

30 :vartype mesh_file: str 

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

32 :vartype regions_file: str 

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

34 :vartype verbose: bool 

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

36 :vartype gu: GmshUtils 

37 """ 

38 self.fdm = fdm 

39 self.cacdm = fdm.magnet 

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

41 self.magnet_name = fdm.general.magnet_name 

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

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

44 self.verbose = verbose 

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

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

47 

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

49 

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

51 """ 

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

53 """ 

54 gmsh.write(self.mesh_file) 

55 if gui: 

56 self.gu.launch_interactive_GUI() 

57 else: 

58 if gmsh.isInitialized(): 

59 gmsh.clear() 

60 gmsh.finalize() 

61 

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

63 """ 

64 Loads a previously generated mesh. 

65 """ 

66 gmsh.clear() 

67 gmsh.open(self.mesh_file) 

68 

69 if gui: 

70 self.gu.launch_interactive_GUI() 

71 

72 

73class StrandMesh(Mesh): 

74 """ 

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

76 

77 :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. 

78 :vartype geometry_class: TwistedStrand 

79 

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

81 :type fdm: object 

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

83 :type verbose: bool, optional 

84 """ 

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 amplitude_dependent_meshing, field_penetration_distance = self.evaluate_amplitude_dependent_meshing().values() 

116 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region 

117 

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

119 if frequency_dependent_meshing: 

120 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] 

121 else: 

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

123 

124 if amplitude_dependent_meshing: 

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

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

127 

128 circular_filament_fields = [] 

129 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) 

130 

131 for filament in filaments_in_skindepth: 

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

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

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

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

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

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

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

139 circular_filament_fields.append(circular_field) 

140 

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

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

143 

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

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

146 

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

148 

149 

150 else: 

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

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

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

154 

155 mass_centers = [] 

156 for layer in self.geometry_class.filaments: 

157 for filament in layer: 

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

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

160 mass_centers.append(mass_center_tag) 

161 gmsh.model.occ.synchronize() 

162 

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

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

165 

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

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

168 gmsh.model.mesh.field.setString( 

169 filament_fields, 

170 "F", 

171 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} )" 

172 ) 

173 

174 # Restrict the filament fields to the filaments 

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

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

177 

178 if frequency_dependent_meshing: 

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

180 else: 

181 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]) 

182 

183 return filament_mesh_field 

184 

185 def matrix_threshold_to_filament_field(self, filaments): 

186 """ 

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

188 

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

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

191 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. 

192 

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

194 :type filaments: list[list[object]] 

195 

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

197 :rtype: int 

198 """ 

199 matrix = self.geometry_class.matrix 

200 matrix_mesh_options = self.cacdm.mesh.matrix 

201 

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

203 # Distance to filament boundaries field: 

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

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

206 

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

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

209 filament_rad = 0 

210 for layer in self.geometry_class.filaments: 

211 for filament in layer: 

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

213 for curve in filament.boundary_curves: 

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

215 if r > filament_rad: 

216 filament_rad = r 

217 

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

219 

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

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

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

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

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

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

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

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

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

229 

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

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

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

233 # matrix_mesh_size_fields.append(meshSize_by_filament_field) 

234 return meshSize_by_filament_field 

235 

236 def evaluate_amplitude_dependent_meshing(self): 

237 """ 

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

239 

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

241 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. 

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

243 

244 :return: A dictionary containing: 

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

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

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

248 """ 

249 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: 

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

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

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

253 

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

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

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

257 mu0 = np.pi * 4e-7 

258 

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

260 

261 # 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: 

262 a = (meshSize_at_filament_center - min_meshSize) / filament_radius 

263 b = min_meshSize 

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

265 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 

266 

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

268 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region 

269 

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

271 print( 

272 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.") 

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

274 

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

276 

277 def evaluate_frequency_dependent_meshing(self): 

278 """ 

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

280 

281 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. 

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

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

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

285 

286 :return: A dictionary containing: 

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

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

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

290 """ 

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

292 matrix_mesh_options = self.cacdm.mesh.matrix 

293 

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

295 

296 mesh_size_matrix_inner = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_inner 

297 mesh_size_matrix_middle = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_middle 

298 mesh_size_matrix_outer = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_outer 

299 

300 matrix = self.geometry_class.matrix 

301 

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

303 mu0 = np.pi * 4e-7 

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

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

306 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 

307 

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

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

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

311 b = mesh_size_matrix_middle 

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

313 

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

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

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

317 b1 = mesh_size_matrix_middle 

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

319 b2 = mesh_size_matrix_middle 

320 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) 

321 

322 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. 

323 elements_in_skindepth = 1e10 

324 

325 desired_elements_in_skindepth = self.cacdm.mesh.matrix.desired_elements_in_skindepth 

326 

327 if elements_in_skindepth < desired_elements_in_skindepth: 

328 print(f"The skindepth 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.") 

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

330 

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

332 

333 def matrix_field(self): 

334 """ 

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

336 - 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. 

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

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

339 

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

341 :rtype: int 

342 """ 

343 matrix_mesh_options = self.cacdm.mesh.matrix 

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

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

346 filament_rad = 0 

347 for layer in self.geometry_class.filaments: 

348 for filament in layer: 

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

350 for curve in filament.boundary_curves: 

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

352 if r > filament_rad: 

353 filament_rad = r 

354 

355 mesh_size_matrix_inner = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_inner 

356 mesh_size_matrix_middle = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_middle 

357 mesh_size_matrix_outer = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_outer 

358 

359 filament_edge_points = [] 

360 for layer in self.geometry_class.filaments: 

361 for filament in layer: 

362 for curve in filament.boundary_curves: 

363 filament_edge_points.append(curve.P1) 

364 filament_edge_points.append(curve.P2) 

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

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

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

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

369 

370 matrix = self.geometry_class.matrix 

371 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 

372 

373 matrix_mesh_size_fields = [] 

374 

375 # 'Distance from matrix center'- field: 

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

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

378 

379 # 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. 

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

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

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

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

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

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

386 

387 # 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). 

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

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

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

391 matrix_mesh_size_fields.append(middle_matrix_field) 

392 

393 # 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. 

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

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

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

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

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

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

400 

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

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

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

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

405 matrix_mesh_size_fields.append(outer_matrix_field) 

406 

407 # 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 

408 # 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. 

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

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

411 matrix_mesh_size_fields.append(meshSize_by_filament_field) 

412 

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

414 

415 if frequency_dependent_meshing: 

416 # 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 

417 # 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. 

418 # print(f"The skindepth 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.") 

419 

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

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

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

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

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

425 

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

427 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] 

428 

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

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

431 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]) 

432 # -- -- # 

433 

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

435 

436 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] 

437 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. 

438 

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

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

441 return matrix_meshSize_field 

442 

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

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

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

446 

447 return matrix_meshSize_field 

448 

449 def air_field(self): 

450 """ 

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

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

453 

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

455 :rtype: int 

456 """ 

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

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

459 filament_rad = 0 

460 for layer in self.geometry_class.filaments: 

461 for filament in layer: 

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

463 for curve in filament.boundary_curves: 

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

465 if r > filament_rad: 

466 filament_rad = r 

467 

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

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

470 

471 if self.geometry_class.air_composition: 

472 air_outer_boundary_curves = self.geometry_class.air_composition.boundary_curves 

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

474 else: 

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

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

477 

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

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

480 

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

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

483 

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

485 gmsh.model.mesh.field.setString( 

486 air_field, 

487 "F", 

488 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} )" 

489 ) 

490 return air_field 

491 

492 def load_geometry_class(self, geom_folder): 

493 """ 

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

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

496 

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

498 :type geom_folder: str 

499 

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

501 :rtype: object 

502 """ 

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

504 

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

506 geom = pickle.load(geom_save_file) 

507 return geom 

508 

509 def generate_mesh(self, geom_folder): 

510 """ 

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

512 

513 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. 

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

515 

516 :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. 

517 :type geom_folder: str 

518 """ 

519 

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

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

522 self.geometry_class.update_tags() 

523 self.geometry_class.add_physical_groups() 

524 

525 total_filament_boundary_distance_field = self.filament_field() 

526 matrixField = self.matrix_field() 

527 airField = self.air_field() 

528 

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

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

531 

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

533 

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

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

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

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

538 

539 # enforce center nodes 

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

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

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

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

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

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

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

547 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)]) # 

548 gmsh.model.occ.synchronize() 

549 # update the physical surface 

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

551 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') 

552 

553 gmsh.model.mesh.generate(2) 

554 

555 def generate_cuts(self): 

556 """ 

557 Computes the cuts for imposing global quantities. 

558 """ 

559 # Computes cuts, as needed for imposing global quantities 

560 if self.verbose: 

561 print('(Co)homology computation') 

562 

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

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

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

566 self.set_periodic_mesh() 

567 else: 

568 # cohomology cut for single air domain 

569 air_region = self.geometry_class.air[0] 

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

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

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

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

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

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

576 

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

578 gmsh.model.mesh.clearHomologyRequests() 

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

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

581 gmsh.plugin.run("HomologyPostProcessing") 

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

583 

584 # b) Cuts for the individual filaments 

585 for layer in self.geometry_class.filaments: 

586 for filament in layer: 

587 filBnd = filament.physical_boundary_tag 

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

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

590 

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

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

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

594 

595 def set_periodic_mesh(self): 

596 """ 

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

598 """ 

599 air = self.geometry_class.air 

600 # link left and right 

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

602 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]) 

603 # link bottom and top 

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

605 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]) 

606 

607 def generate_regions_file(self): 

608 """ 

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

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

611 """ 

612 regions_model = RegionsModel() 

613 

614 """ -- Initialize data -- """ 

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

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

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

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

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

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

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

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

623 

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

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

626 

627 regions_model.insulator.curve.numbers = [] # Filament holes curves physical tags 

628 regions_model.insulator.curve.names = [] # Filament holes curves physical names 

629 

630 regions_model.induced["Matrix"] = RegionsModelFiQuS.Induced()