Coverage for fiqus/mesh_generators/MeshConductorAC_Strand.py: 96%
349 statements
« prev ^ index » next coverage.py v6.4.4, created at 2024-09-01 03:24 +0200
1import os
2import pickle
4import gmsh
5import numpy as np
7from fiqus.data import RegionsModelFiQuS
8from fiqus.utils.Utils import GmshUtils, FilesAndFolders
9from fiqus.data.RegionsModelFiQuS import RegionsModel
11from fiqus.geom_generators.GeometryConductorAC_Strand import TwistedStrand
13occ = gmsh.model.occ
15class Mesh:
16 def __init__(self, fdm, verbose=True):
17 """
18 A base-class used to manage the mesh for the CAC Strand and Rutherford models.
20 :ivar fdm: The fiqus data model for input parameters.
21 :vartype fdm: dict
22 :ivar cacdm: The magnet section from the fdm.
23 :vartype cacdm: object
24 :ivar mesh_folder: The path to the folder where the mesh files are stored.
25 :vartype mesh_folder: str
26 :ivar magnet_name: The name of the magnet model.
27 :vartype magnet_name: str
28 :ivar mesh_file: The path to the .msh file for the mesh.
29 :vartype mesh_file: str
30 :ivar regions_file: The path to the .regions file for the mesh.
31 :vartype regions_file: str
32 :ivar verbose: If True, the class will print additional information during execution.
33 :vartype verbose: bool
34 :ivar gu: An instance of the GmshUtils class for managing the gmsh utility.
35 :vartype gu: GmshUtils
36 """
37 self.fdm = fdm
38 self.cacdm = fdm.magnet
39 self.mesh_folder = os.path.join(os.getcwd())
40 self.magnet_name = fdm.general.magnet_name
41 self.mesh_file = os.path.join(self.mesh_folder, f"{self.magnet_name}.msh")
42 self.regions_file = os.path.join(self.mesh_folder, f"{self.magnet_name}.regions")
43 self.verbose = verbose
44 self.gu = GmshUtils(self.mesh_folder, self.verbose)
45 self.gu.initialize()
47 gmsh.option.setNumber("General.Terminal", verbose)
49 def save_mesh(self, gui: bool = False):
50 """
51 Saves the mesh to a .msh file. If gui is True, the mesh is also loaded in the gmsh GUI.
52 """
53 gmsh.write(self.mesh_file)
54 if gui:
55 self.gu.launch_interactive_GUI()
56 else:
57 gmsh.clear()
58 gmsh.finalize()
60 def load_mesh(self, gui : bool = False):
61 """
62 Loads a previously generated mesh.
63 """
64 gmsh.clear()
65 gmsh.open(self.mesh_file)
67 if gui:
68 self.gu.launch_interactive_GUI()
70class StrandMesh(Mesh):
71 """
72 A subclass of Mesh that handles mesh generation for the twisted strand cross-section geometries.
74 :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.
75 :vartype geometry_class: TwistedStrand
77 :param fdm: The fiqus data model, providing all the necessary data and parameters for mesh generation.
78 :type fdm: object
79 :param verbose: A boolean flag that indicates whether additional information should be printed to the console during operation. Defaults to True.
80 :type verbose: bool, optional
81 """
82 def __init__(self, fdm, verbose=True):
83 super().__init__(fdm, verbose)
84 self.geometry_class : TwistedStrand = None # Class from geometry generation step, used to store all information about the model structure.
86 def filament_field(self):
87 """
88 Generates the filament mesh size field.
90 The standard mesh size field is linearly interpolated between the mesh size at the filament boundary and the center.
91 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.
93 :return: The filament mesh size field tag.
94 :rtype: int
95 """
96 # Define some constants
97 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius.
98 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape.
99 filament_rad = 0
100 for layer in self.geometry_class.filaments:
101 for filament in layer:
102 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag)
103 for curve in filament.boundary_curves:
104 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)])
105 if r > filament_rad:
106 filament_rad = r
107 meshSize_at_filament_boundary = filament_rad * self.cacdm.mesh.filaments.boundary_mesh_size_ratio
108 meshSize_at_filament_center = filament_rad * self.cacdm.mesh.filaments.center_mesh_size_ratio
109 filament_edges = sum([filament.boundary_curves for layer in self.geometry_class.filaments for filament in layer], [])
112 amplitude_dependent_meshing, field_penetration_distance = self.evaluate_amplitude_dependent_meshing().values()
113 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region
115 frequency_dependent_meshing, skin_depth = self.evaluate_frequency_dependent_meshing().values()
116 if frequency_dependent_meshing:
117 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]
118 else:
119 filaments_in_skindepth = [filament for layer in self.geometry_class.filaments for filament in layer]
121 if amplitude_dependent_meshing:
122 # The field only penetrates a few elements. We should thus decrease the mesh size close to the filament edges.
123 # min_mesh_size will be overwritten, while meshSize_at_filament_center will remain.
125 circular_filament_fields = []
126 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)
128 for filament in filaments_in_skindepth:
129 circular_field = gmsh.model.mesh.field.add("Ball")
130 gmsh.model.mesh.field.setNumber(circular_field, "Radius", (filament_rad-field_penetration_distance)*0.85 )
131 gmsh.model.mesh.field.setNumber(circular_field, "XCenter", filament.center_point.pos[0])
132 gmsh.model.mesh.field.setNumber(circular_field, "YCenter", filament.center_point.pos[1])
133 gmsh.model.mesh.field.setNumber(circular_field, "Thickness", (filament_rad-field_penetration_distance)*0.15)
134 gmsh.model.mesh.field.setNumber(circular_field, "VIn", meshSize_at_filament_center)
135 gmsh.model.mesh.field.setNumber(circular_field, "VOut", (field_penetration_distance)/desired_elements_in_field_penetration_region )
136 circular_filament_fields.append(circular_field)
138 for circle_arc_tag in [CA.tag for CA in filament.boundary_curves]:
139 gmsh.model.mesh.setTransfiniteCurve(circle_arc_tag, int(nodesPerFilament/len(filament.boundary_curves)), "Progression", 1)
141 filament_fields = gmsh.model.mesh.field.add("Max")
142 gmsh.model.mesh.field.setNumbers(filament_fields, "FieldsList", circular_filament_fields)
144 self.cacdm.mesh.filaments.boundary_mesh_size_ratio = ((field_penetration_distance)/desired_elements_in_field_penetration_region)/filament_rad
147 else:
148 # Generalized to arbitrary filament shapes (not only circular)
149 f_dist_boundary = gmsh.model.mesh.field.add("Distance")
150 gmsh.model.mesh.field.setNumbers(f_dist_boundary, "CurvesList", [edge.tag for edge in filament_edges])
152 mass_centers = []
153 for layer in self.geometry_class.filaments:
154 for filament in layer:
155 x,y,z = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag)
156 mass_center_tag = gmsh.model.occ.addPoint(x, y, z)
157 mass_centers.append(mass_center_tag)
158 gmsh.model.occ.synchronize()
160 f_dist_center = gmsh.model.mesh.field.add("Distance")
161 gmsh.model.mesh.field.setNumbers(f_dist_center, "PointsList", mass_centers)
163 # Linearly interpolate between the mesh size at the filament center and the mesh size at the filament boundary.
164 filament_fields = gmsh.model.mesh.field.add("MathEval")
165 gmsh.model.mesh.field.setString(
166 filament_fields,
167 "F",
168 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} )"
169 )
172 # Restrict the filament fields to the filaments
173 filament_mesh_field = gmsh.model.mesh.field.add("Restrict")
174 gmsh.model.mesh.field.setNumber(filament_mesh_field, "InField", filament_fields)
177 if frequency_dependent_meshing:
178 gmsh.model.mesh.field.setNumbers(filament_mesh_field, "SurfacesList", [filament.surface_tag for filament in filaments_in_skindepth])
179 else:
180 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 return filament_mesh_field
184 def matrix_threshold_to_filament_field(self, filaments):
185 """
186 Adjusts the matrix mesh size based on the distance from the filaments to improve mesh quality.
188 This method creates a mesh size field that gradually refines the mesh in the matrix as it approaches the filaments.
189 Filaments typically have a finer mesh than the matrix, and a large gradient in mesh size at their boundary can degrade the mesh quality.
190 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 :param filaments: A list of filament surface objects for which the boundary mesh size adjustments are to be made.
193 :type filaments: list[list[object]]
195 :return: The tag of the created mesh size field, which is restricted to the matrix domain.
196 :rtype: int
197 """
198 matrix = self.geometry_class.matrix
199 matrix_mesh_options = self.cacdm.mesh.matrix
201 filament_edges = sum([[CA.tag for CA in filament.boundary_curves] for filament in filaments], []) # List of tags of all filament boundaries
202 # Distance to filament boundaries field:
203 distance_to_filament_field = gmsh.model.mesh.field.add("Distance")
204 gmsh.model.mesh.field.setNumbers(distance_to_filament_field, "CurvesList", filament_edges)
206 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius.
207 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape.
208 filament_rad = 0
209 for layer in self.geometry_class.filaments:
210 for filament in layer:
211 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag)
212 for curve in filament.boundary_curves:
213 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)])
214 if r > filament_rad:
215 filament_rad = r
217 filament_boundary_meshSize = filament_rad * self.cacdm.mesh.filaments.boundary_mesh_size_ratio
219 # Linarily interpolate the mesh size from the filament boundaries to a specified distance from filaments.
220 threshold_to_filament_field = gmsh.model.mesh.field.add("Threshold")
221 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "InField", distance_to_filament_field)
222 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "SizeMin", filament_boundary_meshSize)
223 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "SizeMax", filament_rad*matrix_mesh_options.mesh_size_matrix_ratio_inner)
224 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "DistMin", 0)
225 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "DistMax", filament_rad*matrix_mesh_options.interpolation_distance_from_filaments_ratio)
226 gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "StopAtDistMax", 1)
227 # gmsh.model.mesh.field.setNumber(threshold_to_filament_field, "Sigmoid", 0)
229 meshSize_by_filament_field = gmsh.model.mesh.field.add("Restrict") # Restrict the mesh to the matrix
230 gmsh.model.mesh.field.setNumber(meshSize_by_filament_field, "InField", threshold_to_filament_field)
231 gmsh.model.mesh.field.setNumbers(meshSize_by_filament_field, "SurfacesList", [partition.surface_tag for partition in matrix])
232 # matrix_mesh_size_fields.append(meshSize_by_filament_field)
233 return meshSize_by_filament_field
235 def evaluate_amplitude_dependent_meshing(self):
236 """
237 Evaluates the necessity for amplitude-dependent meshing based on the applied magnetic field's amplitude.
239 This method assesses whether the amplitude of the applied magnetic field is low enough that it only partially penetrates the filaments.
240 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.
241 It calculates the field penetration distance and determines whether an adjustment in the mesh is necessary to ensure accurate simulation results.
243 :return: A dictionary containing:
244 - "amplitude_dependent_meshing" (bool): Indicates whether amplitude-dependent meshing is needed based on the magnetic field's penetration.
245 - "field_penetration_distance" (float): The estimated distance the magnetic field penetrates into the filaments, used to adjust the mesh size accordingly.
246 :rtype: Dict[str, Union[bool, float]]
247 """
248 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:
249 filament_radius = self.geometry_class.filaments[0][0].rad
250 min_meshSize = filament_radius * self.cacdm.mesh.filaments.boundary_mesh_size_ratio
251 meshSize_at_filament_center = filament_radius * self.cacdm.mesh.filaments.center_mesh_size_ratio
253 applied_field_amplitude = self.cacdm.solve.source_parameters.sine.field_amplitude
254 applied_field_amplitude = max(applied_field_amplitude, 0.01) # Limit the amplitude to 0.01 T to avoid excessively fine mesh (Could be removed?).
255 jc = 3e9 #self.cacdm.solve.material_properties.superconducting.jc.constant
256 mu0 = np.pi*4e-7
258 field_penetration_distance = applied_field_amplitude/(mu0*jc) * self.cacdm.mesh.filaments.field_penetration_depth_scaling_factor # Estimate the field penetration distance
260 # 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:
261 a = (meshSize_at_filament_center - min_meshSize) / filament_radius
262 b = min_meshSize
263 # The number of elements in the range x:[x0, x1] is approximated by the function:
264 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 number_of_elements_penetrated_by_field = number_of_elements_in_range(a, b, 0, field_penetration_distance)
267 desired_elements_in_field_penetration_region = self.cacdm.mesh.filaments.desired_elements_in_field_penetration_region
269 if (number_of_elements_penetrated_by_field < desired_elements_in_field_penetration_region) and field_penetration_distance < filament_radius:
270 print(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.")
271 return dict(amplitude_dependent_meshing = True, field_penetration_distance = field_penetration_distance)
273 return dict(amplitude_dependent_meshing = False, field_penetration_distance = 1e10)
275 def evaluate_frequency_dependent_meshing(self):
276 """
277 Evaluates the need for frequency-dependent meshing based on the skin depth due to the applied magnetic field's frequency.
279 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.
280 A fine mesh is required in this region to capture current variations accurately, while a coarser mesh can suffice in the strand center.
281 The method calculates the skin depth using the field frequency and the matrix resistivity, then approximates the number of elements within this depth.
282 If the calculated number is less than desired, it indicates that frequency-dependent meshing adjustments are needed.
284 :return: A dictionary containing:
285 - "frequency_dependent_meshing" (bool): Indicates whether frequency-dependent meshing adjustments are necessary.
286 - "skin_depth" (float): The calculated skin depth, which informs how the mesh should be adjusted to accurately model current variations.
287 :rtype: Dict[str, Union[bool, float]]
288 """
289 if self.cacdm.mesh.matrix.rate_dependent_scaling_matrix and not self.cacdm.geometry.io_settings.load.load_from_yaml:
290 matrix_mesh_options = self.cacdm.mesh.matrix
292 filament_radius = self.geometry_class.filaments[0][0].rad
294 mesh_size_matrix_inner = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_inner
295 mesh_size_matrix_middle = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_middle
296 mesh_size_matrix_outer = filament_radius * matrix_mesh_options.mesh_size_matrix_ratio_outer
298 matrix = self.geometry_class.matrix
300 rho_matrix = 1.81e-10 #self.cacdm.solve.material_properties.matrix.resistivity.constant
301 mu0 = np.pi*4e-7
302 skin_depth = np.sqrt( rho_matrix / (np.pi*mu0 * self.cacdm.solve.source_parameters.sine.frequency ) ) * self.cacdm.mesh.matrix.skindepth_scaling_factor
303 # Function which approximates number of elements in a range given a linearly interpolated mesh size.
304 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
306 # If the skindepth is smaller than the outer matrix, we only evaluate the outer matrix field to approximate number of elements in skindepth.
307 if skin_depth <= matrix[-1].rad -matrix[-2].rad:
308 a = (mesh_size_matrix_outer - mesh_size_matrix_middle) / (matrix[-1].rad - matrix[-2].rad)
309 b = mesh_size_matrix_middle
310 elements_in_skindepth = number_of_elements_in_range(a, b, matrix[-1].rad -matrix[-2].rad - skin_depth, matrix[-1].rad-matrix[-2].rad)
312 elif skin_depth > matrix[-1].rad -matrix[-2].rad and skin_depth < matrix[-1].rad:
313 # If the skindepth is larger than the outer matrix, we evaluate the inner and outer matrix field to approximate number of elements in skindepth.
314 a1 = (mesh_size_matrix_middle - mesh_size_matrix_inner) / (matrix[-2].rad)
315 b1 = mesh_size_matrix_middle
316 a2 = (mesh_size_matrix_outer - mesh_size_matrix_middle) / (matrix[-1].rad - matrix[-2].rad)
317 b2 = mesh_size_matrix_middle
318 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)
320 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.
321 elements_in_skindepth = 1e10
323 desired_elements_in_skindepth = self.cacdm.mesh.matrix.desired_elements_in_skindepth
325 if elements_in_skindepth < desired_elements_in_skindepth:
326 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.")
327 return dict(frequency_dependent_meshing = True, skin_depth = skin_depth)
329 return dict(frequency_dependent_meshing = False, skin_depth = 1e10)
331 def matrix_field(self):
332 """
333 Creates a mesh size field for the matrix as a combination of three fields:
334 - 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.
335 - Adjustments near the filaments to ensure a smooth transition from the filament mesh to the matrix mesh.
336 - Frequency-dependent adjustments to capture skin depth effects.
338 :return: The tag of the composite mesh size field created for the matrix.
339 :rtype: int
340 """
341 matrix_mesh_options = self.cacdm.mesh.matrix
342 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius.
343 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape.
344 filament_rad = 0
345 for layer in self.geometry_class.filaments:
346 for filament in layer:
347 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag)
348 for curve in filament.boundary_curves:
349 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)])
350 if r > filament_rad:
351 filament_rad = r
353 mesh_size_matrix_inner = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_inner
354 mesh_size_matrix_middle = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_middle
355 mesh_size_matrix_outer = filament_rad * matrix_mesh_options.mesh_size_matrix_ratio_outer
357 filament_edge_points = []
358 for layer in self.geometry_class.filaments:
359 for filament in layer:
360 for curve in filament.boundary_curves:
361 filament_edge_points.append(curve.P1)
362 filament_edge_points.append(curve.P2)
363 filament_outer_point = max(filament_edge_points, key=lambda p: np.linalg.norm(p.pos))
364 # Distance from the outermost filament edge to the center of the matrix
365 # This is done instead of referencing the middle matrix radius, to be able to mesh arbitrary geometries.
366 filament_outer_point_dist = np.linalg.norm(filament_outer_point.pos)
368 matrix = self.geometry_class.matrix
369 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
371 matrix_mesh_size_fields = []
373 # 'Distance from matrix center'- field:
374 dist_from_center_field = gmsh.model.mesh.field.add("MathEval")
375 gmsh.model.mesh.field.setString(dist_from_center_field, "F", "sqrt(x^2 + y^2)")
377 # 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.
378 interpField_innerToMiddle = gmsh.model.mesh.field.add("Threshold")
379 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "InField", dist_from_center_field)
380 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "SizeMin", mesh_size_matrix_inner)
381 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "SizeMax", mesh_size_matrix_middle)
382 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "DistMin", 0)
383 gmsh.model.mesh.field.setNumber(interpField_innerToMiddle, "DistMax", filament_outer_point_dist)
385 # 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).
386 middle_matrix_field = gmsh.model.mesh.field.add("Restrict")
387 gmsh.model.mesh.field.setNumber(middle_matrix_field, "InField", interpField_innerToMiddle)
388 gmsh.model.mesh.field.setNumbers(middle_matrix_field, "SurfacesList", [partition.surface_tag for partition in matrix[:-1]])
389 matrix_mesh_size_fields.append(middle_matrix_field)
391 # 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.
392 interpField_middleToOuter = gmsh.model.mesh.field.add("Threshold")
393 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "InField", dist_from_center_field)
394 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "SizeMin", mesh_size_matrix_middle)
395 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "SizeMax", mesh_size_matrix_outer)
396 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "DistMin", filament_outer_point_dist)
397 gmsh.model.mesh.field.setNumber(interpField_middleToOuter, "DistMax", matrix_outer_point_dist)
399 # Restric the outer field to the outer matrix partition:
400 outer_matrix_field = gmsh.model.mesh.field.add("Restrict")
401 gmsh.model.mesh.field.setNumber(outer_matrix_field, "InField", interpField_middleToOuter)
402 gmsh.model.mesh.field.setNumbers(outer_matrix_field, "SurfacesList", [matrix[-1].surface_tag])
403 matrix_mesh_size_fields.append(outer_matrix_field)
406 # 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
407 # 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.
408 # This is done by interpolating the mesh size based on the distance to the filament boundaries.
409 meshSize_by_filament_field = self.matrix_threshold_to_filament_field([filament for layer in self.geometry_class.filaments for filament in layer])
410 matrix_mesh_size_fields.append(meshSize_by_filament_field)
412 frequency_dependent_meshing, skin_depth = self.evaluate_frequency_dependent_meshing().values()
414 if frequency_dependent_meshing:
415 # 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
416 # 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.
417 # 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 circular_field = gmsh.model.mesh.field.add("Ball")
420 gmsh.model.mesh.field.setNumber(circular_field, "Radius", (matrix[-1].rad-skin_depth)*0.85 )
421 # gmsh.model.mesh.field.setNumber(circular_field, "Thickness", (matrix[-1].rad-skin_depth)*0.15)
422 gmsh.model.mesh.field.setNumber(circular_field, "VIn", mesh_size_matrix_inner)
423 gmsh.model.mesh.field.setNumber(circular_field, "VOut", (skin_depth)/self.cacdm.mesh.matrix.desired_elements_in_skindepth )
427 # -- Restrict the mesh size field to the matrix and the filaments not in the skindepth region -- #
428 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]
430 circular_field_restricted = gmsh.model.mesh.field.add("Restrict")
431 gmsh.model.mesh.field.setNumber(circular_field_restricted, "InField", circular_field)
432 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] )
433 # -- -- #
435 self.cacdm.mesh.matrix.mesh_size_matrix_ratio_outer = skin_depth/self.cacdm.mesh.matrix.desired_elements_in_skindepth * 1/filament_rad
437 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]
438 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.
440 matrix_meshSize_field = gmsh.model.mesh.field.add("Min")
441 gmsh.model.mesh.field.setNumbers(matrix_meshSize_field, "FieldsList", [circular_field_restricted, meshSize_by_filament_field])
442 return matrix_meshSize_field
447 # Combine all the matrix mesh size fields by selecting the minimum value of all fields.
448 matrix_meshSize_field = gmsh.model.mesh.field.add("Min")
449 gmsh.model.mesh.field.setNumbers(matrix_meshSize_field, "FieldsList", matrix_mesh_size_fields)
451 return matrix_meshSize_field
453 def air_field(self):
454 """
455 Generates a mesh size field for the air domain surrounding the model.
456 The field linearly interpolates between the mesh size at the outer boundary of the matrix and the air boundary.
458 :return: The tag of the created mesh size field for the air domain.
459 :rtype: int
460 """
461 # Choose the largest distance from the center of mass of the filaments to the filament boundary as the filament radius.
462 # This approximation is done to account for the fact that the filaments may not be circular, but can have any shape.
463 filament_rad = 0
464 for layer in self.geometry_class.filaments:
465 for filament in layer:
466 center_of_mass = gmsh.model.occ.get_center_of_mass(2, filament.surface_tag)
467 for curve in filament.boundary_curves:
468 r = max([np.linalg.norm(np.array(center_of_mass) - curve.P1.pos), np.linalg.norm(np.array(center_of_mass) - curve.P2.pos)])
469 if r > filament_rad:
470 filament_rad = r
472 mesh_size_outer_matrix = filament_rad * self.cacdm.mesh.matrix.mesh_size_matrix_ratio_outer
473 mesh_size_air_boundary = filament_rad * self.cacdm.mesh.air.max_mesh_size_ratio
475 air_outer_boundary_curves = self.geometry_class.air.boundary_curves
476 air_inner_boundary_curves = sum(self.geometry_class.air.inner_boundary_curves, [])
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])
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])
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
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.
497 :param geom_folder: The path to the folder where the geometry class .pkl file is saved.
498 :type geom_folder: str
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')
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
509 def generate_mesh(self, geom_folder):
510 """
511 Generates the mesh for the entire model geometry based on combined mesh size fields.
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.
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 """
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()
525 total_filament_boundary_distance_field = self.filament_field()
526 matrixField = self.matrix_field()
527 airField = self.air_field()
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])
533 gmsh.model.mesh.field.setAsBackgroundMesh(total_meshSize_field)
535 gmsh.option.setNumber("Mesh.MeshSizeFactor", self.cacdm.mesh.scaling_global)
536 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
537 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
538 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
540 gmsh.model.mesh.generate(2)
542 def generate_cuts(self):
543 """
544 Computes the cuts for imposing global quantities.
545 """
546 # Computes cuts, as needed for imposing global quantities
547 if self.verbose:
548 print('(Co)homology computation')
549 # a) Cut for the total current intensity, associated with the strand boundary
550 strand_bnd = self.geometry_class.air.physical_inner_boundary_tags[0]
551 gmsh.model.mesh.addHomologyRequest("Homology", domainTags=[strand_bnd],dims=[1])
552 air = self.geometry_class.air.physical_surface_tag
553 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[air],dims=[1])
554 cuts = gmsh.model.mesh.computeHomology()
555 gmsh.model.mesh.clearHomologyRequests()
556 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains", str(cuts[0][1]))
557 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains2", str(cuts[1][1]))
558 gmsh.plugin.run("HomologyPostProcessing")
560 self.geometry_class.air.strand_bnd_cut_tag = cuts[1][1]+1 # The cut tag for the strand boundary
562 # b) Cuts for the individual filaments
563 for layer in self.geometry_class.filaments:
564 for filament in layer:
565 filBnd = filament.physical_boundary_tag
566 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[filBnd],dims=[1])
567 cuts = gmsh.model.mesh.computeHomology()
569 filaments = sum(self.geometry_class.filaments, [])
570 for i in range(len(cuts)):
571 filaments[i].cut_tag = cuts[i][1]
573 def generate_regions_file(self):
574 """
575 Generates a .regions file for the GetDP solver, containing all necessary information about the model.
576 The regions model contains data about physical surfaces, boundaries, points and cuts and is stored in the mesh folder.
577 """
578 regions_model = RegionsModel()
580 """ -- Initialize data -- """
581 regions_model.powered["Filaments"] = RegionsModelFiQuS.Powered()
582 regions_model.powered["Filaments"].vol.numbers = [] # Filament physical surfaces tags
583 regions_model.powered["Filaments"].vol.names = [] # Filament physical surfaces names
584 regions_model.powered["Filaments"].surf.numbers = [] # Filament physical boundaries tags
585 regions_model.powered["Filaments"].surf.names = [] # Filament physical boundaries names
586 regions_model.powered["Filaments"].cochain.numbers = [] # Filament boundary cut tags
587 regions_model.powered["Filaments"].curve.names = [] # Stores physical points at filament boundary (to fix phi=0)
588 regions_model.powered["Filaments"].curve.numbers = [] # Stores physical points at filament boundary (to fix phi=0)
590 regions_model.powered["Filaments"].surf_in.numbers = [] # Filament holes physical tags
591 regions_model.powered["Filaments"].surf_in.names = [] # Filament holes physical names
594 regions_model.induced["Matrix"] = RegionsModelFiQuS.Induced()
595 regions_model.induced["Matrix"].vol.numbers = [] # Matrix partition physical surfaces tags
596 regions_model.induced["Matrix"].vol.names = [] # Matrix partition physical surfaces names
597 regions_model.induced["Matrix"].surf_in.numbers = [] # Matrix partition physical boundaries tags
598 regions_model.induced["Matrix"].surf_in.names = [] # Matrix partition physical boundaries names
599 regions_model.induced["Matrix"].surf_out.numbers = [] # Strand physical outer boundary tag
600 regions_model.induced["Matrix"].surf_out.names = [] # Strand physical outer boundary name
601 regions_model.induced["Matrix"].cochain.numbers = [] # Strand outer boundary cut tag
602 regions_model.induced["Matrix"].cochain.names = [] # Strand outer boundary cut name
603 """ -- -- """
605 for layer_i, layer in enumerate(self.geometry_class.filaments):
606 for filament_j, filament in enumerate(layer):
607 regions_model.powered["Filaments"].vol.numbers.append(filament.physical_surface_tag) # Surfaces in powered.vol
608 if self.cacdm.geometry.io_settings.load.load_from_yaml:
609 regions_model.powered["Filaments"].vol.names.append(filament.material)
610 else:
611 regions_model.powered["Filaments"].vol.names.append(filament.physical_surface_name)
613 regions_model.powered["Filaments"].surf.numbers.append(filament.physical_boundary_tag) # Boundaries in powered.surf
614 regions_model.powered["Filaments"].surf.names.append(filament.physical_boundary_name)
616 regions_model.powered["Filaments"].cochain.numbers.append(filament.cut_tag)
618 # Add physical point at filament boundary to fix phi=0
619 regions_model.powered["Filaments"].curve.names.append(f"EdgePoint: filament_{layer_i+1}_{filament_j+1}")
620 regions_model.powered["Filaments"].curve.numbers.append(filament.physicalEdgePointTag)
622 for layer in self.geometry_class.filament_holes:
623 for hole in layer:
624 regions_model.powered["Filaments"].surf_in.numbers.append(hole.physical_surface_tag)
625 if self.cacdm.geometry.io_settings.load.load_from_yaml:
626 regions_model.powered["Filaments"].surf_in.names.append(hole.material)
627 else:
628 regions_model.powered["Filaments"].surf_in.names.append(hole.physical_surface_name)
630 for matrixPartition in self.geometry_class.matrix:
631 regions_model.induced["Matrix"].vol.numbers.append(matrixPartition.physical_surface_tag)
632 if self.cacdm.geometry.io_settings.load.load_from_yaml:
633 regions_model.induced["Matrix"].vol.names.append(matrixPartition.material)
634 else:
635 regions_model.induced["Matrix"].vol.names.append(matrixPartition.physical_surface_name)
637 regions_model.induced["Matrix"].surf_in.numbers.append(matrixPartition.physical_boundary_tag)
638 regions_model.induced["Matrix"].surf_in.names.append(matrixPartition.physical_boundary_name)
640 regions_model.induced["Matrix"].surf_out.numbers.append(self.geometry_class.air.physical_inner_boundary_tags[0])
641 regions_model.induced["Matrix"].surf_out.names.append(self.geometry_class.air.physical_inner_boundary_names[0])
643 regions_model.induced["Matrix"].cochain.numbers.append(self.geometry_class.air.strand_bnd_cut_tag)
644 regions_model.induced["Matrix"].cochain.names.append("Strand boundary cut tag")
647 regions_model.air.vol.number = self.geometry_class.air.physical_surface_tag
648 regions_model.air.vol.name = self.geometry_class.air.physical_surface_name
650 regions_model.air.surf.number = self.geometry_class.air.physical_boundary_tag
651 regions_model.air.surf.name = self.geometry_class.air.physical_boundary_name
653 # Add physical point at matrix boundary to fix phi=0
654 regions_model.air.point.names = ["Point at strand boundary"]
655 regions_model.air.point.numbers = [self.geometry_class.air.strand_bnd_physicalEdgePointTag]
657 FilesAndFolders.write_data_to_yaml(self.regions_file, regions_model.model_dump())