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

1import timeit 

2import json 

3import logging 

4import math 

5from enum import Enum 

6import operator 

7import itertools 

8import os 

9import pickle 

10 

11import gmsh 

12import numpy as np 

13 

14from fiqus.data import RegionsModelFiQuS 

15from fiqus.utils.Utils import GmshUtils, FilesAndFolders 

16from fiqus.data.RegionsModelFiQuS import RegionsModel 

17 

18from fiqus.geom_generators.GeometryConductorAC_Rutherford import RutherfordCable 

19from fiqus.mesh_generators.MeshConductorAC_Strand_RutherfordCopy import Mesh 

20from abc import ABC, abstractmethod 

21 

22occ = gmsh.model.occ 

23 

24class CableMesh(Mesh): 

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

26 super().__init__(fdm, verbose) 

27 self.geometry_class : RutherfordCable = None # Class from geometry generation step, used to store all information about tags corresponding to everything... To be changed later (probably) 

28 

29 def load_geometry(self, geom_folder): 

30 """ Generating the geometry file also saves the geometry class as a .pkl file. This geometry class can be loaded to reference the different parts of the geometry for meshing.""" 

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

32 

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

34 geom = pickle.load(geom_save_file) 

35 return geom 

36 

37 def generate_mesh(self, geom_folder): 

38 self.geometry_class : RutherfordCable = self.load_geometry(geom_folder) 

39 self.geometry_class.update_tags() 

40 self.geometry_class.add_physical_groups() 

41 

42 # Mesh size field for strands: 

43 strand_field = gmsh.model.mesh.field.add("Constant") 

44 gmsh.model.mesh.field.setNumbers(strand_field, "SurfacesList", [strand.surface_tag for strand in self.geometry_class.strands]) 

45 gmsh.model.mesh.field.setNumber(strand_field, "VIn", self.cacdm.mesh.strand_mesh_size_ratio * self.fdm.conductors[self.cacdm.solve.conductor_name].strand.diameter) 

46 

47 # Mesh size field for coating: 

48 coating_field = gmsh.model.mesh.field.add("Constant") 

49 gmsh.model.mesh.field.setNumbers(coating_field, "SurfacesList", [self.geometry_class.coating.surface_tag]) 

50 gmsh.model.mesh.field.setNumber(coating_field, "VIn", self.cacdm.mesh.coating_mesh_size_ratio * self.fdm.conductors[self.cacdm.solve.conductor_name].strand.diameter) 

51 

52 # Mesh size field for air: 

53 air_field = gmsh.model.mesh.field.add("Constant") 

54 gmsh.model.mesh.field.setNumbers(air_field, "SurfacesList", [self.geometry_class.air[0].surface_tag]) 

55 gmsh.model.mesh.field.setNumber(air_field, "VIn", self.cacdm.mesh.air_boundary_mesh_size_ratio * self.fdm.conductors[self.cacdm.solve.conductor_name].strand.diameter) 

56 

57 # Mesh size field for excitation coils: 

58 coil_field = gmsh.model.mesh.field.add("Constant") 

59 gmsh.model.mesh.field.setNumbers(coil_field, "SurfacesList", [coil.surface_tag for coil in self.geometry_class.excitation_coils]) 

60 gmsh.model.mesh.field.setNumber(coil_field, "VIn", 2*self.cacdm.mesh.coating_mesh_size_ratio * self.fdm.conductors[self.cacdm.solve.conductor_name].strand.diameter) 

61 

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

63 gmsh.model.mesh.field.setNumbers(total_meshSize_field, "FieldsList", [strand_field, coating_field, air_field, coil_field]) 

64 

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

66 

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

68 gmsh.model.mesh.generate(2) 

69 

70 def generate_cuts(self): 

71 """ 

72 Computes the cuts for imposing global quantities (massive or stranded). 

73 """ 

74 # We need: 

75 # 1) one cut for the coating, 

76 # 2) one cut per strand, which should be directly at the interface between the strand and the coating matrix, 

77 # 3) one cut per excitation coil region. 

78 

79 # 1) Compute cut for the coating region and oriente it based on the surface orientation 

80 cable_outer_boundary = self.geometry_class.coating.physical_boundary_tag 

81 gmsh.model.mesh.addHomologyRequest("Homology", domainTags=[cable_outer_boundary],dims=[1]) 

82 air = self.geometry_class.air[0].physical_surface_tag 

83 coils = [coil.physical_surface_tag for coil in self.geometry_class.excitation_coils] 

84 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[air]+coils,dims=[1]) 

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

86 gmsh.model.mesh.clearHomologyRequests() 

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

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

89 gmsh.plugin.run("HomologyPostProcessing") 

90 self.geometry_class.coating.cut_tag = cuts[1][1] + 1 

91 

92 # 2) Compute cuts and get the oriented surface of each strand 

93 # Cuts are computed one by one for each strand. 

94 # iMPORTANT: the relative cohomology is computed. 

95 # We consider the boundaries of each strand, relative to all the other strand boundaries and inner air regions 

96 # so that we ensure to keep the support of the co-chains in interfaces between strands and coating matrix, 

97 # rather than between strands or adjacent with inner air regions. 

98 # NB: cuts could also be defined one by one, possibly going through other strands and inner air regions. (But it does not seem to work correctly now.) 

99 for strand in self.geometry_class.strands: 

100 gmsh.model.mesh.addHomologyRequest("Homology", domainTags=[strand.physical_surface_tag], subdomainTags=[strand.physical_boundary_tag], dims=[2]) 

101 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[strand.physical_boundary_tag], subdomainTags=[other_strand.physical_boundary_tag for other_strand in self.geometry_class.strands if other_strand != strand]+[air], dims=[1]) 

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

103 gmsh.model.mesh.clearHomologyRequests() 

104 

105 homology = cuts[1::2] # The homology results represent the surfaces of the strands 

106 cuts = cuts[::2] # The (randomly oriented) cuts on the boundary of the strands 

107 

108 # Extract the tags from the homology and cuts 

109 homology = [h[1] for h in homology] 

110 cuts = [c[1] for c in cuts] 

111 

112 # Format the homology and cuts for the plugin 

113 homology_formatted = ', '.join(map(str, homology)) 

114 cuts_formatted = ', '.join(map(str, cuts)) 

115 

116 # Create an identity matrix for the transformation matrix 

117 N = len(self.geometry_class.strands) 

118 identity = '; '.join(', '.join('1' if i == j else '0' for j in range(N)) for i in range(N)) 

119 

120 # Next we get the boundaries of the oriented surfaces, which will also be oriented correctly 

121 # The result of this operation is the oriented boundaries of the strands 

122 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains", homology_formatted) 

123 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains2", "") 

124 gmsh.plugin.setString("HomologyPostProcessing", "TransformationMatrix", identity) 

125 gmsh.plugin.setNumber("HomologyPostProcessing", "ApplyBoundaryOperatorToResults", 1) 

126 gmsh.plugin.run("HomologyPostProcessing") 

127 

128 homology_oriented = [homology[-1] + i + 1 for i in range(N)] 

129 homology_oriented_formatted = ', '.join(map(str, homology_oriented)) 

130 

131 # Finally we get the cuts for each strand as some linear combination of the cuts, which allows inducing currents in each strand separately 

132 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains", homology_oriented_formatted) 

133 gmsh.plugin.setString("HomologyPostProcessing", "PhysicalGroupsOfOperatedChains2", cuts_formatted) 

134 gmsh.plugin.setNumber("HomologyPostProcessing", "ApplyBoundaryOperatorToResults", 0) 

135 gmsh.plugin.run("HomologyPostProcessing") 

136 

137 for i, strand in enumerate(self.geometry_class.strands): 

138 # print(f"Strand {strand.physical_surface_name} cut tag: {homology_oriented[-1] + i + 1}") 

139 strand.cut_tag = int(homology_oriented[-1]) + i + 1 

140 

141 # 3) Compute cuts for excitation coils (will be modelled as stranded conductors) 

142 strands = [strand.physical_surface_tag for strand in self.geometry_class.strands] 

143 coating = self.geometry_class.coating.physical_surface_tag 

144 for coil in self.geometry_class.excitation_coils: 

145 coil_outer_boundary = coil.physical_boundary_tag 

146 gmsh.model.mesh.addHomologyRequest("Homology", domainTags=[coil_outer_boundary],dims=[1]) 

147 other_coils = [other_coil.physical_surface_tag for other_coil in self.geometry_class.excitation_coils if other_coil != coil] 

148 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[air]+other_coils+strands+[coating],dims=[1]) 

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

150 gmsh.model.mesh.clearHomologyRequests() 

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

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

153 gmsh.plugin.run("HomologyPostProcessing") 

154 coil.cut_tag = cuts[1][1] + 1 

155 

156 def generate_regions_file(self): 

157 """ 

158 Generates the .regions file for the GetDP solver. 

159 """ 

160 regions_model = RegionsModel() 

161 

162 """ -- Initialize data -- """ 

163 regions_model.powered["Strands"] = RegionsModelFiQuS.Powered() 

164 regions_model.powered["Strands"].vol.numbers = [] 

165 regions_model.powered["Strands"].vol.names = [] 

166 regions_model.powered["Strands"].surf.numbers = [] 

167 regions_model.powered["Strands"].surf.names = [] 

168 regions_model.powered["Strands"].cochain.numbers = [] 

169 regions_model.powered["Strands"].cochain.names = [] 

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

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

172 

173 regions_model.powered["Coating"] = RegionsModelFiQuS.Powered() 

174 regions_model.powered["Coating"].vol.numbers = [self.geometry_class.coating.physical_surface_tag] 

175 regions_model.powered["Coating"].vol.names = [self.geometry_class.coating.physical_surface_name] 

176 regions_model.powered["Coating"].surf_out.numbers = [self.geometry_class.coating.physical_boundary_tag] 

177 regions_model.powered["Coating"].surf_out.names = [self.geometry_class.coating.physical_boundary_name] 

178 ## regions_model.powered["Coating"].surf_in.numbers = [ ] 

179 ## regions_model.powered["Coating"].surf_in.names = [] 

180 regions_model.powered["Coating"].cochain.numbers = [self.geometry_class.coating.cut_tag] 

181 regions_model.powered["Coating"].cochain.names = [f"Cut: Coating"] 

182 regions_model.powered["Coating"].curve.names = ["EdgePoint: Coating"] # Stores physical points at filament boundary (to fix phi=0) 

183 regions_model.powered["Coating"].curve.numbers = [self.geometry_class.coating.physical_edge_point_tag] # Stores physical points at filament boundary (to fix phi=0) 

184 

185 regions_model.air.point.names = [] 

186 regions_model.air.point.numbers = [] 

187 

188 regions_model.powered["ExcitationCoils"] = RegionsModelFiQuS.Powered() 

189 regions_model.powered["ExcitationCoils"].vol.numbers = [] 

190 regions_model.powered["ExcitationCoils"].vol.names = [] 

191 regions_model.powered["ExcitationCoils"].surf.numbers = [] 

192 regions_model.powered["ExcitationCoils"].surf.names = [] 

193 regions_model.powered["ExcitationCoils"].cochain.numbers = [] 

194 regions_model.powered["ExcitationCoils"].cochain.names = [] 

195 

196 for i, strand in enumerate(self.geometry_class.strands): 

197 regions_model.powered["Strands"].vol.numbers.append(strand.physical_surface_tag) # Surfaces in powered.vol 

198 regions_model.powered["Strands"].vol.names.append(strand.physical_surface_name) 

199 

200 regions_model.powered["Strands"].surf.numbers.append(strand.physical_boundary_tag) # Boundaries in powered.surf 

201 regions_model.powered["Strands"].surf.names.append(strand.physical_boundary_name) 

202 

203 regions_model.powered["Strands"].cochain.numbers.append(strand.cut_tag) 

204 regions_model.powered["Strands"].cochain.names.append(f"Cut: Strand {i}") 

205 

206 # Add physical point at Strand boundary to fix phi=0 

207 regions_model.powered["Strands"].curve.names.append(f"EdgePoint: Strand {i}") 

208 regions_model.powered["Strands"].curve.numbers.append(strand.physical_edge_point_tag) 

209 

210 

211 regions_model.air.vol.number = self.geometry_class.air[0].physical_surface_tag 

212 regions_model.air.vol.name = self.geometry_class.air[0].physical_surface_name 

213 

214 regions_model.air.surf.number = self.geometry_class.air[0].physical_boundary_tag 

215 regions_model.air.surf.name = self.geometry_class.air[0].physical_boundary_name 

216 

217 for air_region in self.geometry_class.air[1:]: 

218 regions_model.air.point.names.append(air_region.physical_boundary_name) 

219 regions_model.air.point.numbers.append(air_region.physical_boundary_tag) 

220 

221 for i, coil in enumerate(self.geometry_class.excitation_coils): 

222 regions_model.powered["ExcitationCoils"].vol.numbers.append(coil.physical_surface_tag) # Surfaces in powered.vol 

223 regions_model.powered["ExcitationCoils"].vol.names.append(coil.physical_surface_name) 

224 

225 regions_model.powered["ExcitationCoils"].surf.numbers.append(coil.physical_boundary_tag) # Boundaries in powered.surf 

226 regions_model.powered["ExcitationCoils"].surf.names.append(coil.physical_boundary_name) 

227 

228 regions_model.powered["ExcitationCoils"].cochain.numbers.append(coil.cut_tag) 

229 regions_model.powered["ExcitationCoils"].cochain.names.append(f"Cut: Coil {i}") 

230 

231 # Add physical point at matrix boundary to fix phi=0 

232 ## regions_model.air.point.names = ["Point at matrix boundary"] 

233 ## regions_model.air.point.numbers = [self.geometry_class.matrix[-1].physicalEdgePointTag] 

234 

235 FilesAndFolders.write_data_to_yaml(self.regions_file, regions_model.model_dump())