Coverage for fiqus/geom_generators/GeometryConductorAC_Rutherford.py: 72%

1import os 

2import math 

3import logging 

4import json 

5import shutil 

6import timeit 

7import pickle 

8 

9import numpy as np 

10import gmsh 

11from itertools import zip_longest 

12 

13import fiqus.data.DataFiQuSConductorAC_Strand as geom 

14from fiqus.utils.Utils import GmshUtils 

15from fiqus.utils.Utils import FilesAndFolders 

16from fiqus.parsers.ParserCOND import ParserCOND 

17from fiqus.data import RegionsModelFiQuS as Reg_Mod_FiQ 

18from fiqus.utils.Utils import FilesAndFolders as UFF 

19from abc import ABC, abstractmethod 

20from pydantic import BaseModel 

21from typing import (Dict, List) 

22 

23import matplotlib.pyplot as plt 

24import matplotlib.patches as patches 

25 

26import fiqus.geom_generators.GeometryConductorAC_Strand_RutherfordCopy as StrandGeom 

27 

28logger = logging.getLogger('FiQuS') 

29 

30class RoundedQuadrilateral(StrandGeom.Surface): 

31 def __init__(self, corner_points, corner_arc_rad): 

32 super().__init__() 

33 self.corner_points = corner_points 

34 self.r = corner_arc_rad 

35 self.boundary_curves = self.generate_boundary_curves() 

36 

37 def generate_boundary_curves(self): 

38 if self.r == 0: # If the trapezoid is not rounded, return the edge lines. 

39 return [StrandGeom.Line.create_or_get(self.corner_points[i], self.corner_points[(i+1)%4]) for i in range(4)] 

40 corner_circle_arcs = [] 

41 boundary_curves = [] 

42 

43 for i, p in enumerate(self.corner_points): # Create the corner circle arcs 

44 p_prev = self.corner_points[(i-1)%4] # Previous point in the list of corner points 

45 p_next = self.corner_points[(i+1)%4] # Next point in the list of corner points 

46 v1 = (p_prev - p) / np.linalg.norm(p_prev - p) # Unit vector from p to p_prev 

47 v2 = (p_next - p) / np.linalg.norm(p_next - p) # Unit vector from p to p_next 

48 angle = np.arccos(np.dot(v1, v2)) # Angle between v1 and v2 

49 circle_arc_center = p + self.r / np.sin(angle/2) * (v1 + v2) / np.linalg.norm(v1 + v2) # Center of the circle arc 

50 arc_p1 = p + v1*np.dot(v1, circle_arc_center - p) # First point on the circle arc 

51 arc_p2 = p + v2*np.dot(v2, circle_arc_center - p) # Second point on the circle arc 

52 # The rounded corners should be represented as circle-arcs, but must be ellipse-arcs to allow for snapping close points together. 

53 # Circle-arcs would need the points to be exact. 

54 arc = StrandGeom.EllipseArc(arc_p1, circle_arc_center, arc_p1, arc_p2) # Create the rounded corner ellipse-arc 

55 if np.linalg.norm(arc.P1.pos-arc.C.pos) < np.linalg.norm(arc.P2.pos-arc.C.pos):# Make sure that the major axis is greater than the minor axis 

56 arc.M = arc.P2 

57 arc.points[-1] = arc.M 

58 

59 if i > 0: 

60 if corner_circle_arcs[i-1].P2 is not arc.P1: 

61 # If the subsequent circle-arcs are not connected we want to connect them with a line element 

62 line = StrandGeom.Line.create_or_get(corner_circle_arcs[i-1].P2, arc.P1) # Create the edge line 

63 boundary_curves.append(line) 

64 

65 boundary_curves.append(arc) 

66 corner_circle_arcs.append(arc) 

67 

68 if corner_circle_arcs[-1].P2 is not corner_circle_arcs[0].P1: 

69 line = StrandGeom.Line.create_or_get(corner_circle_arcs[-1].P2, corner_circle_arcs[0].P1) 

70 boundary_curves.append(line) 

71 

72 return boundary_curves 

73 

74 @classmethod 

75 def get_max_arc_radius(cls, corner_points): 

76 """  

77 Returns the maximum radius of the circle arc that can be used to round the corners. 

78 """ 

79 # With increasing radius of the circle arc, the center of a circle arc will move along the bisector of the angle between the two lines that meet at the corner. 

80 # At some maximum radius, two or more circle arc centers will coincide, and at any larger radius, the circle arcs will overlap, rendering the quadrilateral invalid. 

81 # The maximum radius of the circle-arcs is therefore the minimum radius where two or more circle arc centers coincide. 

82 

83 # 1) Find the bisectors of the angles between the lines that meet at the corners. 

84 bisectors = [] 

85 bisector_angles = [] 

86 for i, p in enumerate(corner_points): 

87 p_prev = corner_points[(i-1)%4] # Previous point in the list of corner points 

88 p_next = corner_points[(i+1)%4] # Next point in the list of corner points 

89 v1 = (p_prev - p) / np.linalg.norm(p_prev - p) # Unit vector from p to p_prev 

90 v2 = (p_next - p) / np.linalg.norm(p_next - p) # Unit vector from p to p_next 

91 angle = np.arccos(np.dot(v1, v2)) # Angle between v1 and v2 

92 bisector_angles.append(angle/2) 

93 bisector = (v1 + v2) / np.linalg.norm(v1 + v2) # unit vector from p towards the circle arc center 

94 bisectors.append(bisector) 

95 

96 # 2) Find the radii for which the circle-arc centers will lie at the points where the bisectors intersect 

97 # The smallest such radius is the maximum radius of the circle arc that can be used to round the corners of the quadrilateral. 

98 r = [] 

99 for i, v1 in enumerate(bisectors): 

100 for j, v2 in enumerate(bisectors[i+1:]): 

101 j = j+i+1 

102 p1 = corner_points[i] 

103 angle1 = bisector_angles[i] 

104 

105 p2 = corner_points[j] 

106 angle2 = bisector_angles[j] 

107 

108 r_max = np.linalg.norm(p2-p1) / np.linalg.norm(v1/np.sin(angle1) - v2/np.sin(angle2)) # The radius for which two circle arc centers coincide at the intersection of the bisectors 

109 r.append(r_max) 

110 return min(r) # The smallest radius at which two or more circle arc centers coincide 

111 

112 @classmethod 

113 def get_area(cls, corner_points, r): 

114 """ 

115 Returns the area of a rounded quadrilateral with the given corner points and corner arc radius. 

116 """ 

117 corner_points = corner_points[:, :2] # Remove the z-coordinate from the corner points 

118 # Area of the quadrilateral without rounded corners is calculated using the shoelace formula: 

119 quadrilateral_area = 0.5 * np.abs(np.dot(corner_points[:, 0], np.roll(corner_points[:, 1], -1)) - np.dot(corner_points[:, 1], np.roll(corner_points[:, 0], -1))) 

120 

121 # When the corners are rounded, the area of the quadrilateral is reduced by r^2 * k(angle/2) (for each corner),  

122 # where r is the radius of the circle arcs and k(angle/2) is a factor that depends on the angle between the two lines that meet at the corner. 

123 k = lambda angle: np.pi/2 - angle - 1/np.tan(angle) 

124 area = quadrilateral_area 

125 for i, p in enumerate(corner_points): 

126 p_prev = corner_points[(i-1)%4] # Previous point in the list of corner points 

127 p_next = corner_points[(i+1)%4] # Next point in the list of corner points 

128 v1 = (p_prev - p) / np.linalg.norm(p_prev - p) # Unit vector from p to p_prev 

129 v2 = (p_next - p) / np.linalg.norm(p_next - p) # Unit vector from p to p_next 

130 angle = np.arccos(np.dot(v1, v2)) # Angle between v1 and v2 

131 area += r**2 * k(angle/2) 

132 

133 return area 

134 

135 @classmethod 

136 def get_arc_radius(cls, desired_area, corner_points): 

137 """ 

138 Returns the radius of the circle arc that is required to create a rounded quadrilateral with the desired area. 

139 """ 

140 corner_points = corner_points[:, :2] 

141 # Area of the quadrilateral without rounded corners is calculated using the shoelace formula: 

142 quadrilateral_area = 0.5 * np.abs(np.dot(corner_points[:, 0], np.roll(corner_points[:, 1], -1)) - np.dot(corner_points[:, 1], np.roll(corner_points[:, 0], -1))) 

143 

144 # When the corners are rounded, the area of the quadrilateral is reduced by r^2 * k(angle/2) for each corner,  

145 # where r is the radius of the circle arcs and k(angle/2) is a factor that depends on the angle between the two lines that meet at the corner. 

146 k = lambda angle: np.pi/2 - angle - 1/np.tan(angle) 

147 angles = [] 

148 for i, p in enumerate(corner_points): 

149 p_prev = corner_points[(i-1)%4] # Previous point in the list of corner points 

150 p_next = corner_points[(i+1)%4] # Next point in the list of corner points 

151 v1 = (p_prev - p) / np.linalg.norm(p_prev - p) # Unit vector from p to p_prev 

152 v2 = (p_next - p) / np.linalg.norm(p_next - p) # Unit vector from p to p_next 

153 angle = np.arccos(np.dot(v1, v2)) # Angle between v1 and v2 

154 angles.append(angle) 

155 r = np.sqrt( (desired_area - (quadrilateral_area+1e-10) )/ sum([k(angle/2) for angle in angles]) ) # Radius of the circle arc that is required to create a rounded trapezoid with the desired area. 

156 return r 

157 

158class RutherfordCable: 

159 def __init__(self) -> None: 

160 self.strands : List[StrandGeom.Surface] = [] 

161 self.air : List[StrandGeom.Surface] = [] # TODO: This needs to be a list of all the air-surfaces (including the air between the strands). They must all be added to the physical group 'Air'. 

162 self.coating : StrandGeom.Surface = None 

163 

164 self.excitation_coils : List[StrandGeom.Surface] = [] 

165 

166 self.cable_boundary_physical_group_tag = None 

167 self.cable_boundary_physical_group_name = None 

168 self.cable_boundary_curve_tags = [] 

169 self.cable_boundary_curve_loop_tag = None 

170 

171 # Coating 

172 # self.coating_curve_tags = [] 

173 # self.coating_surface_tag = None 

174 # self.coating_curve_loop_tag = None 

175 # self.coating_physical_group_tag = None 

176 

177 def create_geometry(self, N_strands_per_layer, cable_width, cable_height_min, cable_height_max, strand_area, min_roundness_fraction, coating_thickness, coating_corner_arc_radius, air_radius, keep_strand_area, coil_center_points, coil_widths, coil_heights): 

178 """ 

179 Creates the geometry of the Rutherford cable. 

180 """ 

181 

182 strand_height_min = cable_height_min / 2 

183 strand_height_max = cable_height_max / 2 

184 

185 if not keep_strand_area: 

186 avg_height = (strand_height_min + strand_height_max) / 2 

187 strand_width = cable_width / N_strands_per_layer 

188 max_arc_rad = RoundedQuadrilateral.get_max_arc_radius(np.array([[0, 0, 0], [strand_width, 0, 0], [strand_width, avg_height, 0], [0, avg_height, 0]])) 

189 strand_area = RoundedQuadrilateral.get_area(np.array([[0, 0, 0], [strand_width, 0, 0], [strand_width, avg_height, 0], [0, avg_height, 0]]), min_roundness_fraction * max_arc_rad) 

190 

191 # The rutherford cable strands will now be created as a series of trapezoids. 

192 height = lambda x: strand_height_min + (strand_height_max-strand_height_min)/cable_width*x # Height of the trapezoids as a function of the x-coordinate 

193 w1 = cable_width/(N_strands_per_layer) # Strand width (unaltered) 

194 delta_w = 0 

195 

196 #1) Check if the first strand can be created without altering the width.  

197 # If the first strand, with minimum arc-radius, has an area below the specified area, the width of the first strand must be increased. 

198 s1_max_arc_rad = RoundedQuadrilateral.get_max_arc_radius(np.array([[0, 0, 0], [w1, 0, 0], [w1, height(w1), 0], [0, height(0), 0]])) 

199 s1_min_arc_rad = min_roundness_fraction*s1_max_arc_rad 

200 s1_area_squared = RoundedQuadrilateral.get_area(np.array([[0, 0, 0], [w1, 0, 0], [w1, height(w1), 0], [0, height(0), 0]]), s1_min_arc_rad) 

201 

202 if s1_area_squared < strand_area - 1e-10: 

203 # print("The first strand cannot be created without altering the width.") 

204 # If the first strand cannot be created without altering the width, the width must be adjusted 

205 # by finding the width that gives the specified area for the first strand, given that it is fully squared. 

206 w1 = cable_width*(np.sqrt(strand_height_min**2 + 2*(strand_area+s1_min_arc_rad**2*(np.pi-2))*(strand_height_max-strand_height_min)/cable_width) - strand_height_min)/(strand_height_max-strand_height_min) # Width of the first strand 

207 delta_w = 2*(cable_width-N_strands_per_layer*w1)/(N_strands_per_layer*(N_strands_per_layer-1)) # The remaining strand widths will be incremented by delta_w to keep the total width constant. 

208 

209 #2) Check if the area of the first strand, when fully rounded, is greater than the specified area. 

210 # If it is, then the geometry can not be made since the smallest strand is greater than the specified area. 

211 s1_rounded_arc_radius = RoundedQuadrilateral.get_max_arc_radius(np.array([[0, 0, 0], [w1, 0, 0], [w1, height(w1), 0], [0, height(0), 0]])) 

212 s1_area_rounded = RoundedQuadrilateral.get_area(np.array([[0, 0, 0], [w1, 0, 0], [w1, height(w1), 0], [0, height(0), 0]]), s1_rounded_arc_radius) 

213 if s1_area_rounded > strand_area: 

214 raise Exception(f"The area of the smallest possible strand is greater than the specified area by a factor of {s1_area_rounded/strand_area}.") 

215 

216 x = lambda n: n*w1 + n*(n-1)/2*delta_w # Function for the x-coordinate of the leftmost point of strand n (the width of strand n is x(n+1)-x(n)) 

217 

218 #3) Check if the last strand, when fully squared, is smaller than the specified area.  

219 # If it is, then the geometry can not be made since the largest strand is smaller than the specified area. 

220 sN_area_squared = RoundedQuadrilateral.get_area(np.array([[x(N_strands_per_layer-1), 0, 0], [x(N_strands_per_layer), 0, 0], [x(N_strands_per_layer), height(x(N_strands_per_layer)), 0], [x(N_strands_per_layer-1), height(x(N_strands_per_layer-1)), 0]]), 0) 

221 if sN_area_squared < strand_area - 1e-10: 

222 raise Exception(f"The area of the last strand is less than the specified area by a factor of {sN_area_squared/strand_area}.") 

223 

224 #4) Check if the last strand, when fully rounded, has an area greater than the specified area.  

225 # If it does, then the geometry can not be made since the last strand is larger than the specified area. 

226 # This problem should only arise if the last strand is elongated in the y-direction. In which case the cable is 'pulled', not compressed. 

227 sN_rounded_arc_radius = RoundedQuadrilateral.get_max_arc_radius(np.array([[x(N_strands_per_layer-1), 0, 0], [x(N_strands_per_layer), 0, 0], [x(N_strands_per_layer), height(x(N_strands_per_layer)), 0], [x(N_strands_per_layer-1), height(x(N_strands_per_layer-1)), 0]])) 

228 sN_area_rounded = RoundedQuadrilateral.get_area(np.array([[x(N_strands_per_layer-1), 0, 0], [x(N_strands_per_layer), 0, 0], [x(N_strands_per_layer), height(x(N_strands_per_layer)), 0], [x(N_strands_per_layer-1), height(x(N_strands_per_layer-1)), 0]]), sN_rounded_arc_radius) 

229 if sN_area_rounded > strand_area: 

230 raise Exception(f"The minimum area of the last strand is greater than the specified area by a factor of {sN_area_rounded/strand_area}.") 

231 

232 #5) Create the strands 

233 for layer in range(2): 

234 for n in range(N_strands_per_layer): 

235 sign = 1 if layer == 0 else -1 # Sign of the height of the strand (the second layer is inverted) 

236 strand_corner_points = np.array([[x(n), 0, 0], [x(n+1), 0, 0], [x(n+1), sign*height(x(n+1)), 0], [x(n), sign*height(x(n)), 0]]) 

237 if layer == 1: 

238 strand_corner_points = strand_corner_points[::-1] # Reverse the order of the corner points in the second layer  

239 r = RoundedQuadrilateral.get_arc_radius(strand_area, strand_corner_points) # Radius of the circle arc that is required to create a rounded trapezoid with the desired area. 

240 if r <= max((x(n+1) - x(n))/1e4, StrandGeom.Point.point_snap_tolerance): 

241 r = 0 # If the radius is very small, it is likely due to numerical errors and should be set to zero. 

242 strand = RoundedQuadrilateral(strand_corner_points, r) 

243 self.strands.append(strand) 

244 

245 # 6) Reverse the second half of the list of strands to make the second layer of strands inverted. This is to have the strands in the correct counter-clockwise order. 

246 self.strands = self.strands[:N_strands_per_layer] + self.strands[N_strands_per_layer:][::-1] 

247 

248 StrandGeom.Surface.replace_overlapping_edges() # Replace overlapping edges of the strands 

249 

250 # 7) Create a region representing a coating around the cable. 

251 # a) Create the corner points of the trapezoid 

252 corner_points = np.array([[-coating_thickness, strand_height_min+coating_thickness, 0], [-coating_thickness, -strand_height_min-coating_thickness, 0], [cable_width+coating_thickness, -strand_height_max-coating_thickness, 0], [cable_width+coating_thickness, strand_height_max+coating_thickness, 0]]) 

253 # b) Create the trapezoid 

254 self.coating = RoundedQuadrilateral(corner_points, coating_corner_arc_radius) 

255 # c) Get the combined boundary curves of the cable 

256 cable_boundary_curves = sum([strand.boundary_curves for strand in self.strands], []) 

257 # Subtract all curves that are shared between multiple strands 

258 shared_curves = set([curve for curve in cable_boundary_curves if cable_boundary_curves.count(curve) > 1]) 

259 cable_boundary_curves = list(set(cable_boundary_curves) - shared_curves) 

260 

261 cable_boundary_loops = StrandGeom.Curve.get_closed_loops(cable_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

262 outer_cable_boundary = max(cable_boundary_loops, key=len) # The outer cable-air boundary is the longest closed loop. 

263 

264 # d) Set the outer cable boundary as the inner boundary of the coating 

265 # self.coating.inner_boundary_curves.append(outer_cable_boundary) # Do not append! Mysterious error will occur... 

266 self.coating.inner_boundary_curves = [outer_cable_boundary] 

267 

268 # 8) Create optional excitation coils with source current 

269 if coil_center_points: 

270 n_coils = len(coil_center_points) 

271 else: 

272 n_coils = 0 

273 for n in range(n_coils): 

274 excitation_coil = StrandGeom.Rectangle(coil_center_points[n], coil_widths[n], coil_heights[n]) 

275 self.excitation_coils.append(excitation_coil) 

276 

277 # 9) Create the air region 

278 # a) Create a disk that represents the outer air region 

279 outer_air_region = StrandGeom.Disk([cable_width/2,0,0], air_radius) 

280 # b) Set the inner boundary of the air region 

281 outer_air_region.inner_boundary_curves.append(self.coating.boundary_curves) 

282 for n in range(n_coils): 

283 outer_air_region.inner_boundary_curves.append(self.excitation_coils[n].boundary_curves) 

284 self.air.append(outer_air_region) 

285 # c) Create the inner air regions (the air between the strands) 

286 air_boundaries = cable_boundary_loops 

287 air_boundaries.remove(outer_cable_boundary) # Remove the outer cable boundary from the list of air boundaries. 

288 for air_boundary in air_boundaries: 

289 air_region = StrandGeom.Surface(air_boundary) 

290 self.air.append(air_region) 

291 

292 def write_geom_to_yaml(self, file_path): 

293 # This function writes the geometry to a yaml file.  

294 # The yaml file contains the coordinates of the points, the type of the curves and the indices of the points that make up the curves and the indices of the curves that make up the areas. 

295 # Note: Only the strands are written to the yaml file. The air region is not included. 

296 strands = self.strands 

297 strand_curves = [curve for strand in strands for curve in strand.boundary_curves] 

298 strand_points = [point for curve in strand_curves for point in curve.points] 

299 

300 points = {} 

301 curves = {} 

302 areas = {} 

303 

304 for p, point in enumerate(StrandGeom.Point.points_registry): 

305 if point in strand_points: 

306 points[p] = {"Coordinates": point.pos.tolist()} 

307 

308 for c, curve in enumerate(StrandGeom.Curve.curves_registry): 

309 if curve in strand_curves: 

310 curves[c] = {"Type": curve.__class__.__name__, "Points": []} 

311 for point in curve.points: 

312 point_id = StrandGeom.Point.points_registry.index(point) 

313 curves[c]["Points"].append(point_id) 

314 

315 for a, area in enumerate(StrandGeom.Surface.surfaces_registry): 

316 if area in strands: 

317 areas[a] = {"Boundary": []} 

318 for curve in area.boundary_curves: 

319 curve_id = StrandGeom.Curve.curves_registry.index(curve) 

320 areas[a]["Boundary"].append(curve_id) 

321 

322 geom_dict = {"Points": points, "Curves": curves, "Areas": areas} 

323 

324 UFF.write_data_to_yaml(file_path, geom_dict) 

325 

326 @classmethod 

327 def read_geom_from_yaml(cls, file_path): 

328 # This function loads the geometry from a yaml file. 

329 # The yaml file contains the coordinates of the points, the type of the curves and the indices of the points that make up the curves  

330 # and the indices of the curves that make up the strands. 

331 geom_data = UFF.read_data_from_yaml(file_path, StrandGeom.Geom) 

332 point_coords = [p.Coordinates for p in geom_data.Points.values()] 

333 for pos in point_coords: 

334 StrandGeom.Point.create_or_get(pos) 

335 

336 for curve in geom_data.Curves.values(): 

337 curve_type = curve.Type 

338 points = [StrandGeom.Point.points_registry[p] for p in curve.Points] 

339 globals()[curve_type].create_or_get(*points) # Creates a curve-entity of the specified type (which is added to the curve-registry). 

340 

341 rutherford_cable = cls() 

342 for area in geom_data.Areas.values(): 

343 boundary_curves = [StrandGeom.Curve.curves_registry[c] for c in area.Boundary] 

344 a = StrandGeom.Surface(boundary_curves) 

345 rutherford_cable.strands.append(a) 

346 

347 return rutherford_cable 

348 

349 def create_gmsh_instance(self): 

350 

351 # When creating the gmsh instances of the strands, curve loops are created for each strand. 

352 # The curve loops are oriented according to the orientation of the first curve in the loop and determine the orientation of the surface.  

353 # We want all the strands to have the same orientation, so we need to make sure that the curve loops are oriented in the same direction. 

354 # To do this, we will check the orientation of the first curve in each strand boundary.  

355 # If the orientation is not counter-clockwise, we will 'roll' the curves in the boundary until the first curve is counter-clockwise. 

356 

357 surfaces = self.strands + [self.coating] + self.air + self.excitation_coils 

358 # StrandGeom.Surface.set_correct_boundary_orientation(surfaces) 

359 

360 for s in surfaces: 

361 # Determine orientation by: 

362 # 1) Find the centroid of the boundary curves. 

363 # 2) Look at the cross-product of the vectors from the centroid to the start and end points of the first curve. 

364 # 3) If the cross-product is positive, the orientation is counter-clockwise. If it is negative, the orientation is clockwise. 

365 # centroid = np.mean(sum([[curve.P1.pos, curve.P2.pos] for curve in s.boundary_curves], []), axis=0) 

366 # first_curve = s.boundary_curves[0] 

367 # v1 = first_curve.P1.pos - centroid 

368 # v2 = first_curve.P2.pos - centroid 

369 # orientation = np.sign(np.cross(v1, v2)[2]) 

370 

371 # if orientation < 0: 

372 # print("here for surface ", s) 

373 # If the orientation is not counter-clockwise, check the next curve in the boundary. 

374 # for i in range(1, len(s.boundary_curves)): 

375 # next_curve = s.boundary_curves[i] 

376 # v1 = next_curve.P1.pos - centroid 

377 # v2 = next_curve.P2.pos - centroid 

378 # orientation = np.sign(np.cross(v1, v2)[2]) 

379 # if orientation > 0: 

380 # s.boundary_curves = s.boundary_curves[i:] + s.boundary_curves[:i] # Roll the curves in the boundary until the first curve is counter-clockwise. 

381 # break 

382 

383 s.create_gmsh_instance() 

384 

385 def add_coating(self): 

386 """  

387 This function adds a coating around the entire cable. The coating is added as a surface that surrounds the cable. 

388 We can add the coating by getting the combined outer boundary of the cable, scaling it up while keeping the center fixed, and creating a surface from the scaled boundary. 

389 Alternatively we can simply add a trapezoid that surrounds the cable. 

390 

391 """ 

392 gmsh.model.occ.synchronize() # Synchronize the OCC model with the GMSH model. 

393 

394 # 1) Get the curves that make up the boundary of the cable. 

395 strands = [strand.surface_tag for strand in self.strands] 

396 strand_comnbined_boundary_tags = [tag for dim, tag in gmsh.model.get_boundary([(2, strand) for strand in strands], combined=True)] 

397 strand_combined_boundary_curves = [StrandGeom.Curve.get_curve_from_tag(tag) for tag in strand_comnbined_boundary_tags] 

398 

399 air_boundaries = StrandGeom.Curve.get_closed_loops(strand_combined_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

400 outer_cable_air_boundary = max(air_boundaries, key=len) # The outer cable-air boundary is the longest closed loop. 

401 

402 self.cable_boundary_curve_tags = [curve.tag for curve in outer_cable_air_boundary] 

403 self.cable_boundary_curve_loop_tag = gmsh.model.occ.addCurveLoop(self.cable_boundary_curve_tags) 

404 

405 

406 # 1.1) Find the strand which each curve in the boundary belongs to and store this in a dictionary. 

407 curve_strand_dict = {} 

408 for curve in outer_cable_air_boundary: 

409 for strand in self.strands: 

410 if curve in strand.boundary_curves: 

411 curve_strand_dict[curve] = strand 

412 break 

413 

414 # 2) Create copies of each curve in the boundary 

415 cable_outline_curves = gmsh.model.occ.copy([(1, curve.tag) for curve in outer_cable_air_boundary]) 

416 

417 

418 # 3) Scale the curves individually by a scaling factor, around the center of mass of the strand they belong to. 

419 scaling_factor = 1.1 

420 

421 for new_curve, curve in zip(cable_outline_curves, outer_cable_air_boundary): 

422 strand = curve_strand_dict[curve] 

423 center_of_mass = gmsh.model.occ.get_center_of_mass(2, strand.surface_tag) 

424 gmsh.model.occ.dilate([(1, new_curve[1])], center_of_mass[0], center_of_mass[1], center_of_mass[2], scaling_factor, scaling_factor, scaling_factor) 

425 

426 # 4) Fuse the scaled curves  

427 fused_curve = gmsh.model.occ.fuse([curve for curve in cable_outline_curves[:len(cable_outline_curves)//2]], [curve for curve in cable_outline_curves[len(cable_outline_curves)//2:]])[0] 

428 logger.info(f"Fused curve: {fused_curve}") 

429 

430 gmsh.model.occ.synchronize() 

431 # 5) Now we remove curves which do not belong to the fused curve. 

432 # These curves have one point which is shared by three curves and one point which is only on one curve. 

433 # 5.1) Create a dictionary with the points and the curves they belong to. 

434 point_curves_dict = {} 

435 for curve in fused_curve: 

436 logger.info(f"Curve: {curve}, Boundary points: {gmsh.model.getBoundary([curve])}") 

437 for point in gmsh.model.getBoundary([curve]): 

438 if point not in point_curves_dict: 

439 point_curves_dict[point] = [] 

440 point_curves_dict[point].append(curve) 

441 logger.info(f"Point curves dict: {point_curves_dict}") 

442 # 5.2) Find the curves corresponding to points which lie only on one curve. 

443 single_point_curves = [curve for point, curves in point_curves_dict.items() if len(curves) == 1 for curve in curves] 

444 logger.info(f"Single point curves: {single_point_curves}") 

445 # 5.3) Remove the single point curves 

446 gmsh.model.occ.remove([curve for curve in single_point_curves]) 

447 

448 # 6) Create a surface from the fused curve 

449 # Remove the deleted curves from the list of curves 

450 self.coating_curve_tags = [curve[1] for curve in fused_curve if curve not in single_point_curves] 

451 self.coating_curve_loop_tag = gmsh.model.occ.addCurveLoop(self.coating_curve_tags) 

452 self.coating_surface_tag = gmsh.model.occ.add_plane_surface([self.coating_curve_loop_tag, self.cable_boundary_curve_loop_tag]) 

453 logger.info(f"Coating surface: {self.coating_surface_tag}") 

454 

455 def add_coating2(self, cable_width, cable_height_min, cable_height_max, coating_thickness): 

456 gmsh.model.occ.synchronize() # Synchronize the OCC model with the GMSH model. 

457 # Create a trapezoidal surface that surrounds the cable. 

458 strand_height_min = cable_height_min / 2 

459 strand_height_max = cable_height_max / 2 

460 # 1) Create the corner points of the trapezoid 

461 corner_points = np.array([[-coating_thickness, strand_height_min+coating_thickness, 0], [-coating_thickness, -strand_height_min-coating_thickness, 0], [cable_width+coating_thickness, -strand_height_max-coating_thickness, 0], [cable_width+coating_thickness, strand_height_max+coating_thickness, 0]]) 

462 # 2) Create the trapezoid 

463 coating = RoundedQuadrilateral(corner_points, 0) 

464 # for c in coating.boundary_curves: # Add the boundary curves gmsh instances to the coating. 

465 # c.create_gmsh_instance() 

466 

467 # 3) Get the outer boundary of the cable 

468 strands = [strand.surface_tag for strand in self.strands] 

469 strand_comnbined_boundary_tags = [tag for dim, tag in gmsh.model.get_boundary([(2, strand) for strand in strands], combined=True)] 

470 strand_combined_boundary_curves = [StrandGeom.Curve.get_curve_from_tag(tag) for tag in strand_comnbined_boundary_tags] 

471 

472 air_boundaries = StrandGeom.Curve.get_closed_loops(strand_combined_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

473 outer_cable_boundary = max(air_boundaries, key=len) # The outer cable-air boundary is the longest closed loop. 

474 

475 # 4) Add the boundary curves of the cable to the coating inner boundary 

476 # coating.inner_boundary_curves.append([curve for curve in outer_cable_boundary]) 

477 

478 

479 # 5) Create the curve loops of the coating 

480 # coating.curve_loop_tag = StrandGeom.Surface.create_or_get_curve_loop(coating.boundary_curves) #gmsh.model.occ.addCurveLoop([c.tag for c in coating.boundary_curves]) 

481 # coating.inner_curve_loop_tags.append(StrandGeom.Surface.create_or_get_curve_loop(outer_cable_boundary)) #(gmsh.model.occ.addCurveLoop([c.tag for c in outer_cable_boundary])) 

482 coating.create_gmsh_instance() 

483 

484 # 6) Create the surface of the coating 

485 # coating.surface_tag = gmsh.model.occ.addPlaneSurface([coating.curve_loop_tag, coating.inner_curve_loop_tags[0]]) 

486 

487 self.coating = coating 

488 

489 

490 

491 # def add_air(self, cable_width, air_radius): 

492 # # Note: This must be done after the strand gmsh instances have been created. 

493 # gmsh.model.occ.synchronize() # Synchronize the OCC model with the GMSH model. 

494 

495 # # Get the curves that make up the boundary of the 'fused' strands. This is the boundary separating the strands from the air. 

496 # strands = [strand.surface_tag for strand in self.strands] 

497 # strand_comnbined_boundary_tags = [tag for dim, tag in gmsh.model.get_boundary([(2, strand) for strand in strands], combined=True)] 

498 # strand_combined_boundary_curves = [StrandGeom.Curve.get_curve_from_tag(tag) for tag in strand_comnbined_boundary_tags]  

499 

500 # air_boundaries = StrandGeom.Curve.get_closed_loops(strand_combined_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

501 # outer_cable_air_boundary = max(air_boundaries, key=len) # The outer cable-air boundary is the longest closed loop. 

502 # air_boundaries.remove(outer_cable_air_boundary) # Remove the outer air boundary from the list of air boundaries. 

503 

504 # outer_air_region = StrandGeom.Disk([cable_width/2,0,0], air_radius) # Create the body of the outer air region. 

505 # for c in outer_air_region.boundary_curves: # Add the boundary curves gmsh instances to the outer air region. 

506 # c.create_gmsh_instance() 

507 # outer_air_region.inner_boundary_curves.append([curve for curve in outer_cable_air_boundary]) 

508 

509 # outer_air_region.curve_loop_tag = gmsh.model.occ.addCurveLoop([c.tag for c in outer_air_region.boundary_curves]) # Create the outer curve loop of the outer air region. 

510 # outer_air_region.inner_curve_loop_tags.append(gmsh.model.occ.addCurveLoop([c.tag for c in outer_cable_air_boundary])) # Create the inner curve loop of the outer air region. This is the outline of the entire cable. 

511 # outer_air_region.surface_tag = gmsh.model.occ.addPlaneSurface([outer_air_region.curve_loop_tag, outer_air_region.inner_curve_loop_tags[0]]) 

512 

513 # self.air.append(outer_air_region) # Add the outer air region to the list of air regions. 

514 

515 # # Create the inner air regions: 

516 # for air_boundary in air_boundaries: 

517 # air_region = StrandGeom.Surface(air_boundary) 

518 # air_region.create_gmsh_instance() 

519 # self.air.append(air_region) 

520 

521 def add_air2(self, cable_width, air_radius): 

522 # Note: This must be done after the strand gmsh instances have been created. 

523 gmsh.model.occ.synchronize() # Synchronize the OCC model with the GMSH model. 

524 

525 # Get the curves that make up the boundary of the 'fused' strands. This is the boundary separating the strands from the air. 

526 strands = [strand.surface_tag for strand in self.strands] 

527 strand_comnbined_boundary_tags = [tag for dim, tag in gmsh.model.get_boundary([(2, strand) for strand in strands], combined=True)] 

528 strand_combined_boundary_curves = [StrandGeom.Curve.get_curve_from_tag(tag) for tag in strand_comnbined_boundary_tags] 

529 

530 air_boundaries = StrandGeom.Curve.get_closed_loops(strand_combined_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

531 outer_cable_air_boundary = max(air_boundaries, key=len) # The outer cable-air boundary is the longest closed loop. 

532 air_boundaries.remove(outer_cable_air_boundary) # Remove the outer air boundary from the list of air boundaries. 

533 

534 outer_air_region = StrandGeom.Disk([cable_width/2,0,0], air_radius) # Create the body of the outer air region. 

535 for c in outer_air_region.boundary_curves: # Add the boundary curves gmsh instances to the outer air region. 

536 c.create_gmsh_instance() 

537 # outer_air_region.inner_boundary_curves.append([curve for curve in outer_cable_air_boundary]) 

538 outer_air_region.inner_boundary_curves.append(self.coating.boundary_curves) 

539 

540 outer_air_region.curve_loop_tag = gmsh.model.occ.addCurveLoop([c.tag for c in outer_air_region.boundary_curves]) # Create the outer curve loop of the outer air region. 

541 # outer_air_region.inner_curve_loop_tags.append(gmsh.model.occ.addCurveLoop([c.tag for c in outer_cable_air_boundary])) # Create the inner curve loop of the outer air region. This is the outline of the entire cable. 

542 outer_air_region.inner_curve_loop_tags.append(self.coating.curve_loop_tag) 

543 outer_air_region.surface_tag = gmsh.model.occ.addPlaneSurface([outer_air_region.curve_loop_tag, outer_air_region.inner_curve_loop_tags[0]]) 

544 

545 self.air.append(outer_air_region) # Add the outer air region to the list of air regions. 

546 

547 # Create the inner air regions: 

548 for air_boundary in air_boundaries: 

549 air_region = StrandGeom.Surface(air_boundary) 

550 air_region.create_gmsh_instance() 

551 self.air.append(air_region) 

552 

553 

554 # def add_air(self, cable_width, air_radius): 

555 # # Note: This must be done after the strand gmsh instances have been created. 

556 # gmsh.model.occ.synchronize() # Synchronize the OCC model with the GMSH model. 

557 

558 # # # Get the curves that make up the boundary of the 'fused' strands. This is the boundary separating the strands from the air. 

559 # strands = [strand.surface_tag for strand in self.strands] 

560 # strand_comnbined_boundary_tags = [tag for dim, tag in gmsh.model.get_boundary([(2, strand) for strand in strands], combined=True)] 

561 # strand_combined_boundary_curves = [StrandGeom.Curve.get_curve_from_tag(tag) for tag in strand_comnbined_boundary_tags]  

562 

563 # air_boundaries = StrandGeom.Curve.get_closed_loops(strand_combined_boundary_curves) # Get the closed loops that make up the boundaries between the air region and the strands. 

564 # outer_cable_air_boundary = max(air_boundaries, key=len) # The outer cable-air boundary is the longest closed loop. 

565 # air_boundaries.remove(outer_cable_air_boundary) # Remove the outer air boundary from the list of air boundaries. 

566 

567 # outer_air_region = StrandGeom.Disk([cable_width/2,0,0], air_radius) # Create the body of the outer air region. 

568 # for c in outer_air_region.boundary_curves: # Add the boundary curves gmsh instances to the outer air region. 

569 # c.create_gmsh_instance() 

570 # # outer_air_region.inner_boundary_curves.append([curve for curve in outer_cable_air_boundary]) 

571 

572 # outer_air_region.curve_loop_tag = gmsh.model.occ.addCurveLoop([c.tag for c in outer_air_region.boundary_curves]) # Create the outer curve loop of the outer air region. 

573 # outer_air_region.inner_curve_loop_tags.append(self.coating_curve_loop_tag) # Create the inner curve loop of the outer air region. This is the outline of the entire cable. 

574 # outer_air_region.surface_tag = gmsh.model.occ.addPlaneSurface([outer_air_region.curve_loop_tag, outer_air_region.inner_curve_loop_tags[0]]) 

575 

576 # self.air.append(outer_air_region) # Add the outer air region to the list of air regions. 

577 

578 # # Create the inner air regions: 

579 # for air_boundary in air_boundaries: 

580 # air_region = StrandGeom.Surface(air_boundary) 

581 # air_region.create_gmsh_instance() 

582 # self.air.append(air_region) 

583 

584 def update_tags(self): 

585 # for strand in self.strands: 

586 # strand.update_tags() 

587 

588 # for air in self.air: 

589 # air.update_tags() 

590 surfaces = self.strands + self.air + [self.coating] + self.excitation_coils 

591 StrandGeom.Surface.update_tags(surfaces) 

592 

593 def add_physical_groups(self): 

594 # Add physical groups for the strands and the air regions. 

595 # The strands each have a physical boundary and a physical surface. 

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

597 strand.add_physical_boundary(name=f"Boundary: Strand {i}") 

598 strand.add_physical_surface(name=f"Surface: Strand {i}") 

599 

600 strand.physical_edge_point_tag = gmsh.model.addPhysicalGroup(0, [strand.boundary_curves[0].P1.tag]) 

601 

602 # The air region is composed of multiple surfaces: The outer air region and the inner air regions between strands. 

603 # The air has an outer physical boundary and all the air regions share a physical surface. 

604 self.air[0].add_physical_boundary(f"Outer boundary: Air") 

605 self.air[0].physical_surface_tag = gmsh.model.add_physical_group(2, [air.surface_tag for air in self.air], name = "Surface: Air") 

606 self.air[0].physical_surface_name = "Surface: Air" 

607 for i, air in enumerate(self.air[1:]): 

608 # air.add_physical_boundary(f"Boundary: Inner air-region {i}") 

609 air.physical_surface_tag = self.air[0].physical_surface_tag 

610 air.physical_surface_name = "Surface: Air" 

611 air.add_physical_boundary(f"Boundary: Inner air-region {i}") 

612 

613 # Add a physical boundary for the entire cable boundary 

614 # self.cable_boundary_physical_group_tag = gmsh.model.addPhysicalGroup(1, , name = "Boundary: Cable") 

615 # self.cable_boundary_physical_group_name = "Boundary: Cable" 

616 

617 # if self.coating: 

618 # Add a physical surface and boundary for the coating 

619 self.coating.add_physical_surface("Surface: Cable coating") 

620 self.coating.add_physical_boundary("Boundary: Cable coating") 

621 self.coating.physical_edge_point_tag = gmsh.model.addPhysicalGroup(0, [self.coating.boundary_curves[0].P1.tag]) 

622 

623 # self.cable_boundary_physical_group_tag = gmsh.model.addPhysicalGroup(1, self.coating_curve_tags, name = "Boundary: Cable coating") 

624 # self.cable_boundary_physical_group_name = "Boundary: Cable coating" 

625 

626 # # Add a physical surface for the coating 

627 # self.coating_physical_group_tag = gmsh.model.addPhysicalGroup(2, [self.coating_surface_tag], name = "Surface: Cable coating") 

628 

629 # Add physical groups for the excitation coil regions 

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

631 coil.add_physical_boundary(name=f"Boundary: Excitation coil {i}") 

632 coil.add_physical_surface(name=f"Surface: Excitation coil {i}") 

633 

634 

635 def save(self, save_file): 

636 with open(save_file, "wb") as geom_save_file: 

637 pickle.dump(self, geom_save_file) 

638 logger.info(f"Geometry saved to file {save_file}") 

639 

640class Geometry(StrandGeom.Geometry): 

641 def generate_cable_geometry(self, gui=False): 

642 """ 

643 Generates the geometry of a Rutherford cable. 

644 """ 

645 # gmsh.option.setNumber("Geometry.Tolerance", 1e-12) 

646 conductor_name = self.cacdm.solve.conductor_name 

647 cable_width = self.fdm.conductors[conductor_name].cable.bare_cable_width 

648 cable_height_min = self.fdm.conductors[conductor_name].cable.bare_cable_height_low 

649 cable_height_max = self.fdm.conductors[conductor_name].cable.bare_cable_height_high 

650 cable_strands_per_layer = self.fdm.conductors[conductor_name].cable.n_strands_per_layers 

651 strand_diameter = self.fdm.conductors[conductor_name].strand.diameter 

652 strand_area = (strand_diameter/2)**2 * np.pi 

653 air_radius = self.cacdm.geometry.air_radius 

654 coating_thickness = self.cacdm.geometry.coating_thickness 

655 coating_corner_arc_radius = self.cacdm.geometry.coating_corner_arc_radius 

656 keep_strand_area = self.cacdm.geometry.keep_strand_area 

657 coil_center_points = self.cacdm.geometry.excitation_coils.centers 

658 coil_widths = self.cacdm.geometry.excitation_coils.widths 

659 coil_heights = self.cacdm.geometry.excitation_coils.heights 

660 StrandGeom.Point.point_snap_tolerance = strand_diameter*self.cacdm.geometry.point_snap_tolerance_relative_to_strand_diameter 

661 min_roundness_factor = self.cacdm.geometry.min_roundness_factor 

662 

663 

664 if self.cacdm.geometry.io_settings.load.load_from_yaml: 

665 filename = self.cacdm.geometry.io_settings.load.filename