Coverage for fiqus/mesh_generators/MeshCCT.py: 95%

1import math 

2import os 

3import timeit 

4import json 

5from typing import List, Any 

6from pathlib import Path 

7 

8import gmsh 

9import numpy as np 

10from fiqus.utils.Utils import FilesAndFolders as uff 

11from fiqus.utils.Utils import GmshUtils 

12from fiqus.data.DataWindingsCCT import WindingsInformation # for volume information 

13from fiqus.data.RegionsModelFiQuS import RegionsModel 

14 

15 

16class Mesh: 

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

18 """ 

19 Class to preparing brep files by adding terminals. 

20 :param fdm: FiQuS data model 

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

22 """ 

23 self.cctdm = fdm.magnet 

24 self.model_folder = os.path.join(os.getcwd()) 

25 self.magnet_name = fdm.general.magnet_name 

26 

27 self.geom_folder = Path(self.model_folder).parent 

28 

29 self.verbose = verbose 

30 regions_file = os.path.join(self.geom_folder, f'{self.magnet_name}.regions') 

31 self.cctrm = uff.read_data_from_yaml(regions_file, RegionsModel) 

32 winding_info_file = os.path.join(self.geom_folder, f'{self.magnet_name}.wi') 

33 self.cctwi = uff.read_data_from_yaml(winding_info_file, WindingsInformation) 

34 self.gu = GmshUtils(self.model_folder, self.verbose) 

35 self.gu.initialize() 

36 self.model_file = f"{os.path.join(self.model_folder, self.magnet_name)}.msh" 

37 self.formers = [] 

38 self.powered_vols = [] 

39 self.air_boundary_tags = [] 

40 

41 def _find_surf(self, volume_tag): 

42 for surf in gmsh.model.getBoundary([(3, volume_tag)], oriented=False): 

43 _, _, zmin, _, _, zmax = gmsh.model.occ.getBoundingBox(*surf) 

44 z = (zmin + zmax) / 2 

45 if math.isclose(z, self.cctdm.geometry.air.z_min, rel_tol=1e-5) or math.isclose(z, self.cctdm.geometry.air.z_max, rel_tol=1e-5): 

46 return surf[1] 

47 

48 def generate_physical_groups(self, gui=False): 

49 if self.verbose: 

50 print('Generating Physical Groups Started') 

51 start_time = timeit.default_timer() 

52 r_types = ['w' for _ in self.cctwi.w_names] + ['f' for _ in self.cctwi.f_names] # needed for picking different pow_surf later on 

53 vol_max_loop = 0 

54 

55 for i, (f_name, r_type, r_name, r_tag) in enumerate(zip(self.cctwi.w_names + self.cctwi.f_names, r_types, self.cctrm.powered.vol.names, self.cctrm.powered.vol.numbers)): 

56 vol_dict = json.load(open(f"{os.path.join(self.geom_folder, f_name)}.vi")) 

57 volumes_file = np.array(vol_dict['all']) 

58 vol_max_file = np.max(volumes_file) 

59 vols_to_use = volumes_file + vol_max_loop 

60 v_tags = (list(map(int, vols_to_use))) 

61 vol_max_loop = vol_max_file + vol_max_loop 

62 self.powered_vols.extend(v_tags) # used later for meshing 

63 gmsh.model.addPhysicalGroup(dim=3, tags=v_tags, tag=r_tag) 

64 gmsh.model.setPhysicalName(dim=3, tag=r_tag, name=r_name) 

65 powered_in_surf = self._find_surf(v_tags[0]) 

66 gmsh.model.addPhysicalGroup(dim=2, tags=[powered_in_surf], tag=self.cctrm.powered.surf_in.numbers[i]) 

67 gmsh.model.setPhysicalName(dim=2, tag=self.cctrm.powered.surf_in.numbers[i], name=self.cctrm.powered.surf_in.names[i]) 

68 powered_out_surf = self._find_surf(v_tags[-1]) 

69 gmsh.model.addPhysicalGroup(dim=2, tags=[powered_out_surf], tag=self.cctrm.powered.surf_out.numbers[i]) 

70 gmsh.model.setPhysicalName(dim=2, tag=self.cctrm.powered.surf_out.numbers[i], name=self.cctrm.powered.surf_out.names[i]) 

71 

72 vol_max_loop = v_tags[-1] 

73 

74 for r_name, r_tag in zip(self.cctrm.induced.vol.names, self.cctrm.induced.vol.numbers): 

75 vol_max_loop += 1 

76 self.formers.append(vol_max_loop) 

77 gmsh.model.addPhysicalGroup(dim=3, tags=[vol_max_loop], tag=r_tag) 

78 gmsh.model.setPhysicalName(dim=3, tag=r_tag, name=r_name) 

79 

80 vol_max_loop += 1 

81 gmsh.model.addPhysicalGroup(dim=3, tags=[vol_max_loop], tag=self.cctrm.air.vol.number) 

82 gmsh.model.setPhysicalName(dim=3, tag=self.cctrm.air.vol.number, name=self.cctrm.air.vol.name) 

83 abt = gmsh.model.getEntities(2)[-3:] 

84 self.air_boundary_tags = [surf[1] for surf in abt] 

85 gmsh.model.addPhysicalGroup(dim=2, tags=self.air_boundary_tags, tag=self.cctrm.air.surf.number) 

86 gmsh.model.setPhysicalName(dim=2, tag=self.cctrm.air.surf.number, name=self.cctrm.air.surf.name) 

87 

88 # air_line_tags = [] 

89 # for air_boundary_tag in self.air_boundary_tags: 

90 # air_line_tags.extend(gmsh.model.getBoundary([(2, air_boundary_tag)], oriented=False)[1]) 

91 # self.air_center_line_tags = [int(np.max(air_line_tags) + 1)] # this assumes that the above found the lines of air boundary but not the one in the middle that is just with the subsequent tag 

92 # gmsh.model.addPhysicalGroup(dim=1, tags=self.air_center_line_tags, tag=self.cctrm.air.line.number) 

93 # gmsh.model.setPhysicalName(dim=1, tag=self.cctrm.air.line.number, name=self.cctrm.air.line.name) 

94 

95 gmsh.model.occ.synchronize() 

96 if gui: 

97 self.gu.launch_interactive_GUI() 

98 if self.verbose: 

99 print(f'Generating Physical Groups Took {timeit.default_timer() - start_time:.2f} s') 

100 

101 def generate_mesh(self, gui=False): 

102 if self.verbose: 

103 print('Generating Mesh Started') 

104 start_time = timeit.default_timer() 

105 # gmsh.option.setNumber("Mesh.AngleToleranceFacetOverlap", 0.01) 

106 gmsh.option.setNumber("Mesh.Algorithm", 5) 

107 gmsh.option.setNumber("Mesh.Algorithm3D", 10) 

108 gmsh.option.setNumber("Mesh.AllowSwapAngle", 20) 

109 line_tags_transfinite = [] 

110 num_div_transfinite = [] 

111 line_tags_mesh: List[Any] = [] 

112 point_tags_mesh = [] 

113 for vol_tag in self.powered_vols: 

114 vol_surf_tags = gmsh.model.getAdjacencies(3, vol_tag)[1] 

115 for surf_tag in vol_surf_tags: 

116 line_tags = gmsh.model.getAdjacencies(2, surf_tag)[1] 

117 # line_tags_mesh.extend(line_tags) 

118 line_lengths = [] 

119 for line_tag in line_tags: 

120 point_tags = gmsh.model.getAdjacencies(1, line_tag)[1] 

121 if line_tag not in line_tags_mesh: 

122 line_tags_mesh.append(line_tag) 

123 x = [] 

124 y = [] 

125 z = [] 

126 for p, point_tag in enumerate(point_tags): 

127 xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.occ.getBoundingBox(0, point_tag) 

128 x.append((xmin + xmax) / 2) 

129 y.append((ymin + ymax) / 2) 

130 z.append((zmin + zmax) / 2) 

131 if point_tag not in point_tags_mesh: 

132 point_tags_mesh.append(point_tag) 

133 line_lengths.append(math.sqrt((x[0] - x[1]) ** 2 + (y[0] - y[1]) ** 2 + (z[0] - z[1]) ** 2)) 

134 # num_div = math.ceil(dist / self.cctdm.mesh_generators.MeshSizeWindings) + 1 

135 l_length_min = np.min(line_lengths) 

136 for line_tag, l_length in zip(line_tags, line_lengths): 

137 aspect = l_length / l_length_min 

138 if aspect > self.cctdm.mesh.MaxAspectWindings: 

139 num_div = math.ceil(l_length / (l_length_min * self.cctdm.mesh.MaxAspectWindings)) + 1 

140 if line_tag in line_tags_transfinite: 

141 idx = line_tags_transfinite.index(line_tag) 

142 num_div_set = num_div_transfinite[idx] 

143 if num_div_set < num_div: 

144 num_div_transfinite[idx] = num_div 

145 else: 

146 line_tags_transfinite.append(line_tag) 

147 num_div_transfinite.append(num_div) 

148 for line_tag, num_div in zip(line_tags_transfinite, num_div_transfinite): 

149 gmsh.model.mesh.setTransfiniteCurve(line_tag, num_div) 

150 

151 gmsh.model.setColor([(1, i) for i in line_tags_mesh], 255, 0, 0) # , recursive=True) # Red 

152 

153 # gmsh.model.mesh.setTransfiniteCurve(self.air_center_line_tags[0], 15) 

154 # gmsh.model.setColor([(1, i) for i in self.air_center_line_tags], 0, 0, 0) 

155 

156 sld = gmsh.model.mesh.field.add("Distance") # straight line distance 

157 gmsh.model.mesh.field.setNumbers(sld, "CurvesList", line_tags_mesh) 

158 # gmsh.model.mesh_generators.field.setNumbers(1, "PointsList", point_tags_mesh) 

159 gmsh.model.mesh.field.setNumber(sld, "Sampling", 100) 

160 slt = gmsh.model.mesh.field.add("Threshold") # straight line threshold 

161 gmsh.model.mesh.field.setNumber(slt, "InField", sld) 

162 gmsh.model.mesh.field.setNumber(slt, "SizeMin", self.cctdm.mesh.ThresholdSizeMin) 

163 gmsh.model.mesh.field.setNumber(slt, "SizeMax", self.cctdm.mesh.ThresholdSizeMax) 

164 gmsh.model.mesh.field.setNumber(slt, "DistMin", self.cctdm.mesh.ThresholdDistMin) 

165 gmsh.model.mesh.field.setNumber(slt, "DistMax", self.cctdm.mesh.ThresholdDistMax) 

166 # gmsh.model.mesh_generators.field.add("Min", 7) 

167 # gmsh.model.mesh_generators.field.setNumbers(7, "FieldsList", [slt]) 

168 gmsh.model.mesh.field.setAsBackgroundMesh(slt) 

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

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

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

172 gmsh.option.setNumber("Mesh.OptimizeNetgen", 0) 

173 gmsh.model.mesh.generate(3) 

174 if self.verbose: 

175 print(f'Generating Mesh Took {timeit.default_timer() - start_time:.2f} s') 

176 if gui: 

177 self.gu.launch_interactive_GUI() 

178 

179 def generate_cuts(self, gui=False): 

180 if self.verbose: 

181 print('Generating Cuts Started') 

182 start_time = timeit.default_timer() 

183 for vol, surf_in, surf_out in zip(self.cctrm.powered.vol.numbers, self.cctrm.powered.surf_in.numbers, self.cctrm.powered.surf_out.numbers): 

184 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[vol], subdomainTags=[surf_in, surf_out], dims=[1, 2, 3]) 

185 for vol in self.cctrm.induced.vol.numbers: 

186 gmsh.model.mesh.addHomologyRequest("Cohomology", domainTags=[vol], dims=[1, 2, 3]) 

187 gmsh.model.mesh.computeHomology() 

188 if self.verbose: 

189 print(f'Generating Cuts Took {timeit.default_timer() - start_time:.2f} s') 

190 if gui: 

191 self.gu.launch_interactive_GUI() 

192 

193 def save_mesh(self, gui=False): 

194 if self.verbose: 

195 print('Saving Mesh Started') 

196 start_time = timeit.default_timer() 

197 gmsh.write(self.model_file) 

198 if self.verbose: 

199 print(f'Saving Mesh Took {timeit.default_timer() - start_time:.2f} s') 

200 if gui: 

201 self.gu.launch_interactive_GUI() 

202 else: 

203 gmsh.clear() 

204 gmsh.finalize() 

205 

206 def load_mesh(self, gui=False): 

207 gmsh.open(self.model_file) 

208 if gui: 

209 self.gu.launch_interactive_GUI()