Coverage for fiqus/geom_generators/GeometryConductorAC_Strand.py: 85%
574 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-27 02:39 +0100
1import os
2import pickle
4import numpy as np
5import gmsh
7import fiqus.data.DataFiQuSConductor as geom
8from fiqus.utils.Utils import GmshUtils
9from fiqus.utils.Utils import FilesAndFolders as UFF
10from abc import ABC, abstractmethod
11from typing import (Dict, List)
15### HELPER FUNCTIONS ###
17def cylindrical_to_cartesian(rad, angle, height):
18 """
19 Convert cylindrical coordinates to Cartesian coordinates.
21 :return: A list of Cartesian coordinates [x, y, z].
22 :rtype: list
23 """
24 return [rad*np.cos(angle), rad*np.sin(angle), height]
26def cartesian_to_cylindrical(x, y, z):
27 """
28 Convert Cartesian coordinates to cylindrical coordinates.
30 :return: A list of cylindrical coordinates [rad, angle, z].
31 :rtype: list
32 """
33 rad = np.sqrt(x**2 + y**2)
34 angle = np.arctan2(y, x)
35 return [rad, angle, z]
37def rotate_vector(vector, angle):
38 """
39 Rotate a 3D vector in the xy-plane by a given angle.
41 :param vector: The 3D vector to rotate. It should be a list or numpy array of three numbers.
42 :type vector: list or numpy.ndarray
43 :param angle: The angle by which to rotate the vector, in radians.
44 :type angle: float
46 :return: The rotated vector.
47 :rtype: numpy.ndarray
48 """
49 rotation_matrix = np.array([[np.cos(angle), -np.sin(angle), 0], [np.sin(angle), np.cos(angle), 0], [0, 0, 1]])
50 return np.matmul(rotation_matrix, vector)
52### END HELPER FUNCTIONS ###
55### GEOMETRY ###
56class Point:
57 """
58 A class to represent a point in 3D space.
60 :cvar points_registry: A list of all points created. This is a class-level attribute.
61 :vartype points_registry: list
62 :cvar point_snap_tolerance: The tolerance for snapping points to existing points. This is a class-level attribute.
63 :vartype point_snap_tolerance: float
65 :ivar pos: A numpy array representing the position of the point in 3D space. This is an instance-level attribute.
66 :vartype pos: numpy.ndarray
67 :ivar tag: A gmsh tag for the point, default is None. This is an instance-level attribute.
68 :vartype tag: int, optional
69 """
70 points_registry = [] # Global registry of points
71 point_snap_tolerance = 1e-10 # Tolerance for snapping points to existing points
72 def __init__(self, pos : List[float]):
73 """
74 Constructs the attributes for the point object.
76 :param pos: A list of three floats representing the position of the point in 3D space.
77 :type pos: List[float]
78 """
79 self.pos = np.array(pos)
80 self.tag = None
82 @classmethod
83 def create_or_get(cls, pos : any):
84 """
85 Creates a new point if no point with the given position exists, otherwise returns the existing one to assure unique points.
87 :param pos: The position of the point.
88 :type pos: list
90 :return: A point object.
91 :rtype: Point
92 """
93 new_point = cls(pos)
94 if new_point in cls.points_registry:
95 return cls.points_registry[cls.points_registry.index(new_point)]
97 cls.points_registry.append(new_point)
98 return new_point
100 def create_gmsh_instance(self, meshSize : float = 0, tag : int = -1):
101 if self.tag == None:
102 self.tag = gmsh.model.occ.addPoint(self.pos[0], self.pos[1], self.pos[2], meshSize, tag)
103 else:
104 pass
106 def __repr__(self) -> str:
107 return f"Point({self.pos})"
109 def __eq__(self, o: object) -> bool:
110 # Two points sharing the same position (within a margin of point_snap_tolerance) should be considered the same point.
111 if isinstance(o, Point):
112 return np.linalg.norm(np.array(self.pos) - np.array(o.pos)) < self.point_snap_tolerance
113 return False
115 def __hash__(self) -> int:
116 return hash(tuple(self.pos))
118### CURVE ELEMENTS ###
120class Curve(ABC):
121 """
122 Abstract base class for curves in 3D space.
124 :cvar curves_registry: A list of all curves created. This is a class-level attribute.
125 :vartype curves_registry: list
127 :ivar P1: The start point of the curve. This is an instance-level attribute.
128 :vartype P1: Point
129 :ivar P2: The end point of the curve. This is an instance-level attribute.
130 :vartype P2: Point
131 :ivar points: A list of points used by the curve. This is an instance-level attribute.
132 :vartype points: list
133 :ivar tag: A tag for the curve. This is an instance-level attribute.
134 :vartype tag: int
135 """
136 curves_registry = []
137 def __init__(self) -> None:
138 """
139 Constructs the attributes for the curve object.
140 """
141 self.P1 : Point = None
142 self.P2 : Point = None
143 self.points = []
144 self.tag = None
146 # Curve.curves_registry.append(self)
148 def set_points(self, *points : Point):
149 """
150 Sets the points used by the curve.
152 :param points: The points to set.
153 :type points: Point
154 """
155 self.points = list(points)
157 @abstractmethod
158 def create_gmsh_instance(self):
159 pass
161 @classmethod
162 @abstractmethod
163 def create_or_get(cls, *points):
164 pass
166 @classmethod
167 def get_curve_from_tag(cls, tag):
168 """
169 Returns the curve with the given tag.
171 :param tag: The tag of the curve.
172 :type tag: int
174 :return: The curve with the given tag, or None if no such curve exists.
175 :rtype: Curve or None
176 """
177 for curve in cls.curves_registry:
178 if curve.tag == tag:
179 return curve
180 return None
183 @classmethod
184 def get_closed_loops(cls, curves):
185 """
186 Returns a list of lists of curves, where each list of curves defines a closed loop.
187 The function assumes that no two closed loops share a curve and that chains of curves do not split into multiple chains.
189 :param curves: A list of curves.
190 :type curves: list
192 :return: A list of lists of curves, where each list of curves defines a closed loop.
193 :rtype: list
194 """
195 closed_loops = []
197 def get_curve_link(current_curve, curves, curve_link):
198 # Recursive function to find a link of connected curves.
199 # Curves is a list of curves to search through, not containing the current curve.
200 # Curve_link is a list of curves that are connected to the current curve.
201 # Curves are removed from 'curves' as they are added to the curve_link.
202 # The function returns when no more curves are connected to the current link.
203 for curve in curves:
204 for point in [current_curve.P1, current_curve.P2]:
205 if point in [curve.P1, curve.P2]:
206 curve_link.append(curve)
207 curves.remove(curve)
208 return get_curve_link(curve, curves, curve_link)
209 return curves, curve_link # Return the remaining curves and the curve link
211 while len(curves) > 0:
212 curve0 = curves[0]
213 curves.remove(curve0)
214 curves, curve_link = get_curve_link(curve0, curves, [curve0])
216 first_curve_points = set([curve_link[0].P1, curve_link[0].P2])
217 last_curve_points = set([curve_link[-1].P1, curve_link[-1].P2])
219 # TODO: Have to add the case where there is only a single curve in the link. This should not be considered a closed loop.
220 if len(curve_link)>2 and len(first_curve_points&last_curve_points) == 1: # Check if the last curve in the link is connected to the first curve in the link. If so, the link is a closed loop.
221 closed_loops.append(curve_link)
222 # If the link only contains two curves, the curves must be connected at both ends to form a closed loop.
223 elif len(curve_link) == 2 and len(first_curve_points&last_curve_points) == 2:
224 closed_loops.append(curve_link)
226 return closed_loops
228class CircleArc(Curve):
229 """
230 A class to represent a circular arc, subclass of Curve.
232 :ivar P1: The first point of the arc. This is an instance-level attribute.
233 :vartype P1: Point
234 :ivar P2: The second point of the arc. This is an instance-level attribute.
235 :vartype P2: Point
236 :ivar C: The center of the arc. This is an instance-level attribute.
237 :vartype C: Point
238 """
239 def __init__(self, P1 : any, C : any, P2 : any) -> None:
240 """
241 Constructs the attributes for the arc object.
243 :param P1: The start point of the arc.
244 :type P1: list or Point
245 :param P2: The end point of the arc.
246 :type P2: list or Point
247 :param C: The center of the arc.
248 :type C: list or Point
249 """
250 super().__init__()
251 self.P1 = P1 if isinstance(P1, Point) else Point.create_or_get(P1)
252 self.P2 = P2 if isinstance(P2, Point) else Point.create_or_get(P2)
253 self.C = C if isinstance(C, Point) else Point.create_or_get(C)
255 self.set_points(self.P1, self.C, self.P2)
257 Curve.curves_registry.append(self)
259 def create_gmsh_instance(self, tag : int = -1):
260 if self.tag == None: # If the curve has not been created yet
261 self.P1.create_gmsh_instance()
262 self.P2.create_gmsh_instance()
263 self.C.create_gmsh_instance()
264 self.tag = gmsh.model.occ.addCircleArc(self.P1.tag, self.C.tag, self.P2.tag, tag)
266 @classmethod
267 def create_or_get(cls, P1 : any, C : any, P2 : any):
268 return cls(P1, C, P2)
270 def __repr__(self) -> str:
271 return f"CircleArc({self.P1.pos}, {self.P2.pos})"
273class Line(Curve):
274 """
275 A class to represent a line, subclass of Curve.
277 :ivar P1: the start point of the line
278 :type P1: list or Point
279 :ivar P2: the end point of the line
280 :type P2: list or Point
282 """
283 def __init__(self, P1 : any, P2 : any) -> None:
284 """
285 Constructs the attributes for the line object.
286 P1 and P2 can be either 'Point' objects or lists of three floats representing the position of the point.
288 :param P1: the start point of the line
289 :type P1: list or Point
290 :param P2: the end point of the line
291 :type P2: list or Point
292 """
293 super().__init__()
294 self.P1 = P1 if isinstance(P1, Point) else Point.create_or_get(P1)
295 self.P2 = P2 if isinstance(P2, Point) else Point.create_or_get(P2)
297 self.set_points(self.P1, self.P2)
299 def get_length(self):
300 return np.linalg.norm(self.P1.pos - self.P2.pos)
302 def __repr__(self) -> str:
303 return f"Line from {tuple(self.P1.pos)} to {tuple(self.P2.pos)}"
305 @classmethod
306 def create_or_get(cls, P1 : any, P2 : any):
307 """
308 Creates a new line if it doesn't exist, or returns the existing one.
310 :param P1: the first point of the line
311 :type P1: Point
312 :param P2: the second point of the line
313 :type P2: Point
315 :return: a line object
316 :rtype: Line
317 """
319 new_line = cls(P1, P2)
320 if new_line in cls.curves_registry:
321 return cls.curves_registry[cls.curves_registry.index(new_line)]
323 cls.curves_registry.append(new_line)
324 return new_line
326 @classmethod
327 def is_colinear(cls, line1, line2):
328 """
329 Checks if two lines are colinear.
330 """
331 l1 = line1.P2.pos - line1.P1.pos
332 l2 = line2.P2.pos - line2.P1.pos
333 # Check if the lines are parallel
334 if np.linalg.norm(np.cross(l1, l2)) < 1e-10:
335 # Check if the lines are colinear
336 if np.linalg.norm(np.cross(l1, line2.P1.pos - line1.P1.pos)) < 1e-10:
337 return True
338 return False
340 @classmethod
341 def remove_from_registry(cls, lines):
342 """
343 Removes a list of lines from the registry.
345 :param lines: A list of lines to remove.
346 :type lines: list
347 """
348 for line in lines:
349 if line in cls.curves_registry:
350 cls.curves_registry.remove(line)
352 def create_gmsh_instance(self, tag : int = -1):
353 if self.tag == None:
354 self.P1.create_gmsh_instance()
355 self.P2.create_gmsh_instance()
356 self.tag = gmsh.model.occ.addLine(self.P1.tag, self.P2.tag, tag)
358 def __eq__(self, o: object) -> bool:
359 # Two Line-entities specified by the same two points should be considered equal lines.
360 # Check if the other object is a line
361 if isinstance(o, Line):
362 # Check if the lines have the same points
363 return (self.P1 == o.P1 and self.P2 == o.P2) or (self.P1 == o.P2 and self.P2 == o.P1)
364 return False
366 def __hash__(self) -> int:
367 return hash(frozenset([self.P1, self.P2])) # Frozenset is hashable and order-independent
369class EllipseArc(Curve):
370 """
371 A class to represent an elliptical arc, subclass of Curve.
373 :ivar P1: The start point of the arc. This is an instance-level attribute.
374 :vartype P1: Point
375 :ivar P2: The end point of the arc. This is an instance-level attribute.
376 :vartype P2: Point
377 :ivar C: The center point of the arc. This is an instance-level attribute.
378 :vartype C: Point
379 :ivar M: The major axis point of the arc (point anywhere on the major axis). This is an instance-level attribute.
380 :vartype M: Point
381 """
382 def __init__(self, P_start : any, P_center : any, P_major : any, P_end : any) -> None:
383 """
384 Initializes an elliptical arc object.
386 :param P_start: The start point of the arc. If not a Point instance, attempts to retrieve or create one.
387 :type P_start: any
388 :param P_end: The end point of the arc. If not a Point instance, attempts to retrieve or create one.
389 :type P_end: any
390 :param P_center: The center point of the arc. If not a Point instance, attempts to retrieve or create one.
391 :type P_center: any
392 :param P_major: The major axis point of the arc (point anywhere on the major axis). If not a Point instance, attempts to retrieve or create one.
393 :type P_major: any
395 :rtype: None
396 """
397 super().__init__()
398 self.P1 = P_start if isinstance(P_start, Point) else Point.create_or_get(P_start) # Start point
399 self.P2 = P_end if isinstance(P_end, Point) else Point.create_or_get(P_end) # End point
400 self.C = P_center if isinstance(P_center, Point) else Point.create_or_get(P_center) # Center point
401 self.M = P_major if isinstance(P_major, Point) else Point.create_or_get(P_major) # Major axis point (point anywhere on the major axis)
403 self.set_points(self.P1, self.C, self.M, self.P2)
405 Curve.curves_registry.append(self)
406 # self.tag = None
408 def create_gmsh_instance(self, tag : int = -1):
409 if self.tag == None:
410 self.P1.create_gmsh_instance()
411 self.P2.create_gmsh_instance()
412 self.C.create_gmsh_instance()
413 self.M.create_gmsh_instance()
414 self.tag = gmsh.model.occ.addEllipseArc(self.P1.tag, self.C.tag, self.M.tag, self.P2.tag, tag)
416 def __repr__(self) -> str:
417 return f"EllipseArc({self.P1.pos}, {self.P2.pos})"
419 @classmethod
420 def create_or_get(cls, P_start : any, P_center : any, P_major : any, P_end : any):
421 return cls(P_start, P_center, P_major, P_end)
424### SURFACE ELEMENTS ###
426class Surface:
427 """
428 A class to represent a surface in 3D space.
430 :cvar surfaces_registry: A class-level attribute that keeps track of all created surfaces.
431 :vartype surfaces_registry: list
432 :cvar curve_loop_registry: A class-level attribute that keeps track of all curve loops. Each curve loop is stored as a dict with the keys 'curves' and 'tag'. This list is necessary when creating surfaces with holes. The curve-loops of the holes may already have been created when creating the surface of the hole, in which case we get the tag of the existing curve-loop instead of creating a new one.
433 :vartype curve_loop_registry: list
435 :ivar boundary_curves: A list of Curve objects that define the closed outer boundary of the surface. This is an instance-level attribute.
436 :vartype boundary_curves: list
437 :ivar inner_boundary_curves: A list of lists of Curve objects. Each list of curves defines the closed boundary of a hole in the surface. This is an instance-level attribute.
438 :vartype inner_boundary_curves: list
439 :ivar curve_loop_tag: A unique identifier for the curve loop of the outer boundary, initially set to None. This is an instance-level attribute.
440 :vartype curve_loop_tag: int
441 :ivar inner_curve_loop_tags: A list of unique identifiers for the inner curve loops. Each tag corresponds to the curve-loop of a hole in the surface. This is an instance-level attribute.
442 :vartype inner_curve_loop_tags: list
443 :ivar surface_tag: A unique identifier for the surface, initially set to None. This is an instance-level attribute.
444 :vartype surface_tag: int
445 :ivar physical_boundary_tag: A unique identifier for the physical group of the boundary curves, initially set to None. This is an instance-level attribute.
446 :vartype physical_boundary_tag: int
447 :ivar physical_boundary_name: A name for the physical group of the boundary curves, initially set to None. This is an instance-level attribute.
448 :vartype physical_boundary_name: str
449 :ivar physical_inner_boundary_tags: A list of unique identifiers for the physical groups of the inner boundary curves. This is an instance-level attribute.
450 :vartype physical_inner_boundary_tags: list
451 :ivar physical_inner_boundary_names: A list of names for the physical groups of the inner boundary curves. This is an instance-level attribute.
452 :vartype physical_inner_boundary_names: list
453 :ivar physical_surface_tag: A unique identifier for the physical group of the surface, initially set to None. This is an instance-level attribute.
454 :vartype physical_surface_tag: int
455 :ivar physical_surface_name: A name for the physical group of the surface, initially set to None. This is an instance-level attribute.
456 :vartype physical_surface_name: str
457 :ivar material: The material-ID of the surface. Used to specify material properties in the .regions file, initially set to None. This is an instance-level attribute.
458 :vartype material: any
459 """
461 surfaces_registry = []
462 curve_loop_registry = []
463 def __init__(self, boundary_curves : List[Curve] = [], inner_boundary_curves : List[List[Curve]] = []):
464 """
465 Constructs the attributes for the surface object.
467 :param boundary_curves: A list of Curve objects that define the outer boundary of the surface. The curves must form a closed loop.
468 :type boundary_curves: list[Curve]
469 :param inner_boundary_curves: A list of lists of Curve objects. Each list of curves defines the boundary of a hole in the surface.
470 :type inner_boundary_curves: list[list[Curve]]
471 """
472 self.boundary_curves = boundary_curves
473 self.inner_boundary_curves = inner_boundary_curves
475 self.curve_loop_tag = None
476 self.inner_curve_loop_tags : List[int] = []
477 self.surface_tag = None
479 self.physical_boundary_tag = None
480 self.physical_boundary_name = None
481 self.physical_inner_boundary_tags = []
482 self.physical_inner_boundary_names = []
484 self.physical_surface_tag = None
485 self.physical_surface_name = None
487 self.material = None
489 Surface.surfaces_registry.append(self)
491 def create_gmsh_instance(self):
492 if self.surface_tag == None:
493 for curve in self.boundary_curves + sum(self.inner_boundary_curves, []):
494 curve.create_gmsh_instance()
497 self.curve_loop_tag = self.create_or_get_curve_loop(self.boundary_curves)
498 self.inner_curve_loop_tags = [self.create_or_get_curve_loop(curves) for curves in self.inner_boundary_curves]
499 self.surface_tag = gmsh.model.occ.add_plane_surface([self.curve_loop_tag] + self.inner_curve_loop_tags)
501 def add_physical_boundary(self, name : str = ''):
502 self.physical_boundary_tag = gmsh.model.add_physical_group(1, [curve.tag for curve in self.boundary_curves], name = name)
503 self.physical_boundary_name = name
505 def add_physical_surface(self, name : str = ''):
506 self.physical_surface_tag = gmsh.model.add_physical_group(2, [self.surface_tag], name = name)
507 self.physical_surface_name = name
510 @classmethod
511 def update_tags(cls, surfaces):
512 """
513 Updates the tags of model entities of dimension lower than 2.
515 When saving the geometry as a .brep file, these tags may be (seemingly arbitrarily) changed. This method ensures that the tags are consistent and correctly assigned.
517 Steps:
518 1) Find the tags of all outer boundary curves of every surface.
519 2) Divide the boundary curves into groups. Each group of curves are either on the boundary of a single surface or on the intersection of the same multiple surfaces.
520 3) Update the tags of the curves in the curve groups. The tags are assigned not to their corresponding curve but to any curve in the group.
521 4) Update the tags of the points. Points which are in multiple curve-groups are assigned to a new group. Point tags are assigned based on the groups they are in.
523 :param surfaces: A list of Surface objects to update the tags for.
524 :type surfaces: list[Surface]
525 """
527 gmsh.model.occ.synchronize()
528 #1) Find the tags of all outer boundary curves of every surface
529 #1.1) Find the inner surfaces of each surface with holes
530 surfaces_inner_surface_indices = [] # List of lists of indices of surfaces. Each list of indices corresponds to the surfaces that are inside the surface with the same index.
531 for surface in surfaces:
532 inner_surface_indices = []
533 if surface.inner_boundary_curves: # If the surface has holes
534 for inner_boundary in surface.inner_boundary_curves:
535 for surface in surfaces: # Loop through all surfaces to find the surfaces that are inside the outer surface.
536 if set(inner_boundary) & set(surface.boundary_curves): # If the two sets of curves have any common curves, the inner surface is inside the outer surface.
537 inner_surface_indices.append(surfaces.index(surface)) # Add the index of the inner surface to the list of inner surfaces.
538 surfaces_inner_surface_indices.append(inner_surface_indices)
540 #1.2) Find the tags of the outer boundary curves of each surface by finding the boundary of the surface (inner and outer) and removing the boundary curves of the inner surfaces.
541 surfaces_outer_boundary_tags = [set([abs(tag) for dim, tag in gmsh.model.get_boundary([(2, surface.surface_tag)])]) for surface in surfaces]
542 for surface_i in range(len(surfaces)):
543 if surfaces[surface_i].inner_boundary_curves:
544 surface_inner_surfaces_boundary_tags = set.union(*[surfaces_outer_boundary_tags[i] for i in surfaces_inner_surface_indices[surface_i]]) # The boundary tags of the inner surfaces
545 surfaces_outer_boundary_tags[surface_i] = surfaces_outer_boundary_tags[surface_i] - surface_inner_surfaces_boundary_tags
547 #2) Divide the boundary curves into groups. Each group of curves are either on the boundary of a single surface or on the intersection of multiple surfaces, but not both.
548 # We also find the tags of the curves in each group.
549 surface_outer_boundary_curves = [set(surface.boundary_curves) for surface in surfaces]
551 curve_groups = [] # List of sets of curves. Each set of curves is a group, defined as curves which lie only on the boundary of a single surface or on the intersection of multiple surfaces.
552 curve_group_tags = [] # List of sets of tags. Each set of tags corresponds to the set of curves in curve_groups.
553 for i, si_curves in enumerate(surface_outer_boundary_curves):
554 for j, sj_curves in enumerate(surface_outer_boundary_curves[i+1:]):
555 j += i+1
556 common_curves = si_curves & sj_curves
557 if common_curves:
558 curve_groups.append(common_curves) # Add the common curves to the list of curve groups
559 si_curves -= common_curves # Remove the common curves from the surface boundary curves
560 sj_curves -= common_curves # Remove the common curves from the surface boundary curves
562 curve_group_tags.append(surfaces_outer_boundary_tags[i] & surfaces_outer_boundary_tags[j]) # Add the tags of the common curves to the list of curve group tags
563 surfaces_outer_boundary_tags[i] -= curve_group_tags[-1] # Remove the tags of the common curves from the surface boundary tags
564 surfaces_outer_boundary_tags[j] -= curve_group_tags[-1] # Remove the tags of the common curves from the surface boundary tags
566 curve_groups.append(si_curves) # Add the remaining curves to the list of curve groups
567 curve_group_tags.append(surfaces_outer_boundary_tags[i]) # Add the tags of the remaining curves to the list of curve group tags
569 #3) Update the tags of the curves in the curve groups
570 # The tags are assigned not to their corresponding curve but to any curve in the group.
571 for group, group_tags in zip(curve_groups, curve_group_tags):
572 for curve, tag in zip(group, group_tags):
573 curve.tag = tag
575 #4) We have now updated the tags of all curves. Next we update the tags of the points.
577 #4.1) Get all points in each group of curves and the tags of the points in each group of curves
578 curve_groups_points = [set([point for curve in group for point in [curve.P1, curve.P2]]) for group in curve_groups]
579 curve_groups_point_tags = [set(sum([list(gmsh.model.get_adjacencies(1, tag)[1]) for tag in group_tags], [])) for group_tags in curve_group_tags]
581 #4.2) Points which are in multiple curve-groups are assigned to a new group.
582 # These points will be removed from the groups they are in and assigned to the new group, based on the groups they are in.
583 # Iterate trough all points and check which groups they are in. Points in same groups will be assigned to a new group.
584 all_points = set.union(*curve_groups_points)
585 point_new_groups = {} # Dictionary with keys as tuples of indices of the groups the point is in, and values as lists of points in the same groups.
586 for point in all_points:
587 groups_point_is_in = [i for i, group in enumerate(curve_groups_points) if point in group]
588 # If the point is in multiple groups, remove the point from all groups
589 if len(groups_point_is_in) > 1:
590 for i in groups_point_is_in:
591 curve_groups_points[i].remove(point) # Remove the point from all groups, as it will be assigned to a new group.
592 # Sort the groups the point is in, make it a tuple and use it as a key in a dictionary. The value is a list of all points in the same groups.
593 point_new_groups[tuple(sorted(groups_point_is_in))] = point_new_groups.get(tuple(sorted(groups_point_is_in)), []) + [point]
595 #4.3) Update the tags of the points in the new groups
596 # Get the tags of all points in each group of points as the boundary of the group of curves
597 for group_indices, points in point_new_groups.items():
598 # The tags of the points in the new group is the intersection of the tags of the points in the groups the point is in.
599 point_tags = set.intersection(*[curve_groups_point_tags[i] for i in group_indices])
601 # Update the tags of the points in the group
602 for point, point_tag in zip(points, point_tags):
603 point.tag = point_tag
605 # Remove the tags of the points in the new group from the tags of the points in the groups the point was in before it was assigned to the new group.
606 for i in group_indices:
607 curve_groups_point_tags[i] -= point_tags
609 #4.4) Update the tags of points in the remaining groups 'curve_groups_points'
610 for group, group_tags in zip(curve_groups_points, curve_groups_point_tags):
611 for point, tag in zip(group, group_tags):
612 point.tag = tag
617 @classmethod
618 def remove_from_registry(cls, surfaces):
619 for surface in surfaces:
620 cls.surfaces_registry.remove(surface)
622 @classmethod
623 def create_or_get_curve_loop(cls, curves):
624 """
625 Creates a curve loop if it does not already exist, otherwise returns the existing curve loop.
627 A curve loop is a sequence of curves that form a closed loop. This method checks if a curve loop with the given curves already exists in the curve_loop_registry.
628 If it does, the method returns the tag of the existing curve loop.
629 If it doesn't, the method creates a new curve loop, adds it to the curve_loop_registry, and returns its tag.
631 :param curves: A list of Curve objects that define the curve loop.
632 :type curves: list[Curve]
634 :return: The tag of the curve loop.
635 :rtype: int
636 """
637 for curve_loop in cls.curve_loop_registry:
638 if curve_loop['curves'] == set(curves):
639 return curve_loop['tag']
641 tag = gmsh.model.occ.addCurveLoop([curve.tag for curve in curves])
642 cls.curve_loop_registry.append({'tag': tag, 'curves': set(curves)})
643 return tag
645 @classmethod
646 def replace_overlapping_edges(cls):
647 """
648 Replaces overlapping boundary-curves of type 'Line' across all surfaces.
650 If multiple surface boundaries contain lines which are overlapping, the overlapping lines are replaced by the fraction
651 of the line which is unique to the surface as well as the fraction of the line which is shared with the other surface.
653 Steps:
654 1) Sort all existing lines into groups of lines which are colinear.
655 2) For each group of colinear lines, make a list of all the unique points used, sorted by their location on the colinear line.
656 3) For each line of each surface, replace the line by fragments of the line, defined by the points in the list from step 2.
657 """
658 #1) Sort all existing lines into groups of lines which are colinear
659 lines = [line for line in Curve.curves_registry if type(line) == Line]
660 if len(lines) == 0: # If there are no lines, skip the rest of the function
661 return
662 colinear_groups = [[lines[0]]]
663 for line in lines[1:]:
664 for group in colinear_groups:
665 if Line.is_colinear(line, group[0]):
666 group.append(line)
667 break
668 else:
669 colinear_groups.append([line])
671 #2) For each group of colinear lines, make a list of all the unique points used, sorted by their location on the colinear line
672 colinear_groups_points = []
673 for group in colinear_groups:
674 points = []
675 for line in group:
676 points.append(line.P1)
677 points.append(line.P2)
678 points = list(set(points)) # Remove duplicates
679 angle = np.arctan2((line.P2.pos-line.P1.pos)[1], (line.P2.pos-line.P1.pos)[0]) # Angle of the lines with respect to the x-axis
680 positions = [p.pos - points[0].pos for p in points] # Move the points so that the colinear line passes through the origin
681 positions = [rotate_vector(p, -angle)[0] for p in positions] # Rotate the points so that the lines are parallel to the x-axis, making it a 1D problem.
682 points = np.array(points)[np.argsort(positions)].tolist() # Sort the points by their position along the colinear line
683 colinear_groups_points.append(points)
685 #3) For each line for each surface, replace the line by fragments of the line, defined by the points in the list from step 2
686 for area in cls.surfaces_registry: # For each rectangle
687 for l, line in enumerate(area.boundary_curves): # For each line
688 for i, group in enumerate(colinear_groups): # Find the group of colinear lines that the line belongs to
689 if line in group:
690 if len(group) == 1: # If the line is not colinear with any other line, skip it
691 break
692 # Find the points that define the line fragments
693 points = colinear_groups_points[i]
694 line_point1_index = points.index(line.P1)
695 line_point2_index = points.index(line.P2)
696 if line_point1_index > line_point2_index:
697 # If the points orientation differs from the line orientation, reverse the list of points
698 # The points are sorted by their position along the colinear line, which is not necessarily the same as the orientation of the line.
699 # The orientation of the fragments must match the orientation of the line.
700 points.reverse()
701 line_point1_index = points.index(line.P1)
702 line_point2_index = points.index(line.P2)
703 line_fragment_points = points[line_point1_index:line_point2_index+1] # The points that define the line fragments
705 # Create the line fragments
706 line_fragments = [Line.create_or_get(line_fragment_points[i], line_fragment_points[i+1]) for i in range(len(line_fragment_points)-1)]
708 if len(line_fragments) > 1: # If the line is split into multiple fragments, remove the old line and insert the new fragments into the boundary curves of the surface.
709 Line.remove_from_registry([line]) # Remove the old line from the lines-registry.
710 for fragment in line_fragments:
711 area.boundary_curves.insert(area.boundary_curves.index(line), fragment) # Insert the new line fragments into the boundary curves of the surface.
712 area.boundary_curves.remove(line) # Remove the original line from the boundary curves of the surface.
714 @classmethod
715 def set_correct_boundary_orientation(self, surfaces):
716 """
717 When creating surfaces with holes, the boundaries must have a specific orientation to ensure that the surfaces are correctly defined.
718 This method sets the correct orientation of the boundaries of the surfaces.
719 The orientation of the boundary is determined by the orientation of the curve-loop, which is determined by the orientation of the first curve in the loop.
720 We therefore must ensure that the first curve of the curve-loop is oriented correctly.
721 """
722 # It seems to work to just set all the boundaries to have the same orientation. Very simple and does not require any complex logic.
723 for surface in surfaces:
724 outer_boundary_points = sum([[curve.P1, curve.P2] for curve in surface.boundary_curves], [])
725 mean_point = np.mean([point.pos for point in outer_boundary_points], axis=0) # Mean point of the boundary (Center of mass)
726 mean_to_P1 = surface.boundary_curves[0].P1.pos - mean_point # Vector from the mean point to the first point of the boundary curve
727 mean_to_P2 = surface.boundary_curves[0].P2.pos - mean_point # Vector from the mean point to the second point of the boundary curve
729 # Using the two vectors we can check the orientation of the first curve in the loop, based on the sign of the determinant of the matrix formed by the two vectors.
730 if np.linalg.det(np.column_stack((mean_to_P1[:-1], mean_to_P2[:-1]))) < 0:
731 # If the determinant is negative we reverse the orientation so that it is positive (counter-clockwise orientation)
732 surface.boundary_curves[0].P1, surface.boundary_curves[0].P2 = surface.boundary_curves[0].P2, surface.boundary_curves[0].P1
734class Disk(Surface):
735 """
736 A class to represent a disk. Inherits from the Surface class.
738 :ivar rad: The radius of the disk. This is an instance-level attribute.
739 :vartype rad: float
740 :ivar partitions: The number of partitions to divide the boundary of the disk into when generating boundary curves. This is an instance-level attribute.
741 :vartype partitions: int
742 :ivar center_point: The center point of the disk. Can be a position ([x, y, z]) or a 'Point' object. This is an instance-level attribute.
743 :vartype center_point: Point
744 :ivar physicalEdgePointTag: A unique identifier for the physical group of the edge point. Used in pro-template for fixing phi=0 in the outer matrix boundary as well as on the filament boundaries. This is an instance-level attribute.
745 :vartype physicalEdgePointTag: int
746 """
747 def __init__(self, center_point : any, rad : float, partitions : int = 4):
748 """
749 Constructs the attributes for the disk object.
751 :param center_point: The center point of the disk. Can be a position ([x, y, z]) or a 'Point' object.
752 :type center_point: any
753 :param rad: The radius of the disk.
754 :type rad: float
755 :param partitions: The number of partitions to divide the boundary of the disk into when generating boundary curves, defaults to 4
756 :type partitions: int, optional
757 """
758 super().__init__()
759 self.rad = rad
760 self.partitions = partitions
762 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
763 self.boundary_curves = self.generate_boundary_curves()
764 self.inner_boundary_curves = []
766 self.physicalEdgePointTag = None # Used in pro-template for fixing phi=0 in the outer matrix boundary as well as on the filament boundaries.
768 def generate_boundary_curves(self):
769 """
770 Generates the boundary curves for the disk.
772 This method divides the boundary of the disk into 'partitions' number of segments and creates a CircleArc for each segment.
774 :return: A list of CircleArc objects that define the boundary of the disk.
775 :rtype: list
776 """
777 edgePoints = [np.array(self.center_point.pos) + np.array(cylindrical_to_cartesian(self.rad, 2*np.pi * n/self.partitions + np.pi + cartesian_to_cylindrical(self.center_point.pos[0], self.center_point.pos[1], self.center_point.pos[2])[1], 0)) for n in range(self.partitions)]
778 boundary_curves = [CircleArc(edgePoints[n], self.center_point, edgePoints[(n+1)%self.partitions]) for n in range(self.partitions)]
779 return boundary_curves
781class Hexagon(Surface):
782 """
783 A class to represent a hexagon. Inherits from the Surface class.
785 :ivar rad: The radius of the hexagon.
786 :vartype rad: float
787 :ivar center_point: The center point of the hexagon. Can be a position ([x, y, z]) or a 'Point' object.
788 :vartype center_point: Point
789 :ivar rotation: The rotation of the hexagon in radians.
790 :vartype rotation: float
791 """
792 def __init__(self, center_point : any, rad : float, rotation : float = 0):
793 """
794 Constructs the attributes for the hexagon object.
796 :param center_point: The center point of the hexagon. Can be a position ([x, y, z]) or a 'Point' object.
797 :type center_point: any
798 :param rad: The radius of the hexagon.
799 :type rad: float
800 :param rotation: The rotation of the hexagon in radians.
801 :type rotation: float
802 """
803 super().__init__()
804 self.rad = rad
805 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
806 self.rotation = rotation
807 self.boundary_curves = self.generate_boundary_curves()
808 self.inner_boundary_curves = []
810 def generate_boundary_curves(self):
811 """
812 Generates the boundary curves for the hexagon.
814 This method creates the boundary of the hexagon by rotating a vector from the center-point to the first edge-point.
816 :return: A list of Line objects that define the boundary of the hexagon.
817 :rtype: list
818 """
819 edgePoints = [np.array(self.center_point.pos) + rotate_vector(np.array([self.rad, 0, 0]), 2*np.pi * n/6 + self.rotation) for n in range(6)]
820 boundary_curves = [Line.create_or_get(edgePoints[n], edgePoints[(n+1)%6]) for n in range(6)]
821 return boundary_curves
823### END GEOMETRY ###
826class TwistedStrand:
827 """
828 A class to represent a 2D cross section of a twisted strand.
830 :ivar filaments: Each list of surfaces represent six filaments in the same layer.
831 :vartype filaments: list[list[Surface]]
832 :ivar Matrix: List of surfaces corresponding to the matrix partitions.
833 :vartype Matrix: list[Surface]
834 :ivar Air: A surface representing the air region.
835 :vartype Air: Surface
836 """
837 def __init__(self) -> None:
838 """
839 Initializes the TwistedStrand object.
840 """
841 self.filaments : List[List[Surface]] = []
842 self.filament_holes : List[List[Surface]] = [] # If the filaments have holes (only for certain geometries obtained via YAML files)
843 self.matrix : List[Surface] = [] # Inner, middle, outer
844 self.air : Surface = None
846 def create_filament_center_points(self, N_points : int, filament_radius : float, outer_boundary_radius : float, inner_boundary_radius : float, circular_filament_distribution = True):
847 """
848 Creates the center-points of N_points filaments. The points are distributed in layers---the first layer containing 6 points, the second 12, the third 18, etc---and groups of
849 6 points satisfying rotational symmetry when rotated by pi/3. The first layer is composed of one group, the second of two groups, etc.
851 :param N_points: The number of points to create.
852 :type N_points: int
853 :param filament_radius: The radius of the filament.
854 :type filament_radius: float
855 :param outer_boundary_radius: The radius of the outer boundary.
856 :type outer_boundary_radius: float
857 :param inner_boundary_radius: The radius of the inner boundary.
858 :type inner_boundary_radius: float
859 :param circular_filament_distribution: If True, points are distributed in a circular pattern. If False, points are distributed in a hexagonal pattern. Defaults to True.
860 :type circular_filament_distribution: bool, optional
862 :return: A list of lists of lists of points representing the filament center-points for each filament in each group. Each point is a list of three coordinates [x, y, z] and each list of points represents a group of 6 points satisfying rotational symmetry.
863 :rtype: list[list[list[float]]]
864 """
865 # The hexagonal grid can be created by transforming the cartesian coordinate grid to a hexagonal grid.
866 # Point centers are created as unit vectors in the hexagonal grid, and then transformed to cartesian coordinates.
867 # This function takes as input a vector in the hexagonal grid and returns its distance from the origin in the cartesian coordinate system.
868 point_relative_magnitude = lambda v : np.sqrt(v[0]**2 + 2*v[0]*v[1]*np.cos(np.pi/3) + v[1]**2)
869 def create_first_hexant_points(N_points_per_hexant, outer_inner_boundary_ratio, possible_next_points = np.array([[1,0]]), points = np.empty((0,2))):
870 """
871 Recursive algorithm used to create the first hexant of an hexagonal point grid. These points can subsequently be rotated by n*pi/3 (for n in [1,5]) radians to create the next points.
872 The function returns the points in the first sextant of the matrix, sorted by distance from the origin.
874 The points are created in an iterative manner, starting from the points closest to the origin and moving outwards.
875 Selecting the next point is done either by selecting the point with the smallest distance from the origin (giving a circular distribution), or by prioritizing points in the
876 same layer as the previous point, choosing the point closest to the 30 degree angle (giving an hexagonal distribution). If no points are available in the same layer, the next layer is considered.
877 Choosing points by distance from the origin results in a more circular distribution of points, while choosing points by angle results in a hexagonal distribution of points.
878 """
880 if circular_filament_distribution:
881 # The next point is the one with the smallest distance from the origin
882 next_point = min(possible_next_points, key=point_relative_magnitude)
884 else:
885 #1) Filter possible next points to include only points in the lowest not full layer.
886 next_points_in_same_layer = possible_next_points[possible_next_points.sum(axis=1) == possible_next_points.sum(axis=1).min()]
887 #2) Choose the next point as the one with an angle closest to 30 degrees.
888 next_point = min(next_points_in_same_layer, key=lambda v : np.abs(v[1]-v[0]))
890 possible_next_points = np.delete(possible_next_points, np.where(np.all(possible_next_points == next_point, axis=1)), axis=0) # Remove the selected point from possible_next_points
891 points = np.append(points, np.array([next_point]), axis=0) # Add the selected point to points
892 points = sorted(points, key=lambda p: point_relative_magnitude(p)) # Sort the points by their distance from the origin
894 new_possible_next_points = np.array([next_point + np.array([1,0]), next_point + np.array([0,1])]) # The possible next points from the current point
895 possible_next_points = np.unique( np.concatenate( (new_possible_next_points, possible_next_points) ), axis = 0 ) # Add the new possible next points to the list of possible next points
898 if len(points) == N_points_per_hexant:
899 # Check if the outer-inner boundary ratio is satisfied.
900 # The outermost points are always placed at the outer boundary, if the outer-inner boundary ratio is not satisfied, the innermost points must be within the inner boundary.
901 # The innermost point is then removed and replaced by a point in the outermost layer.
902 if point_relative_magnitude(points[-1]) / point_relative_magnitude(points[0]) <= outer_inner_boundary_ratio:
903 return points
904 else:
905 N_points_per_hexant += 1
906 points = np.delete(points, 0, axis=0)
907 return create_first_hexant_points(N_points_per_hexant, outer_inner_boundary_ratio, possible_next_points, points)
909 else:
910 return create_first_hexant_points(N_points_per_hexant, outer_inner_boundary_ratio, possible_next_points, points)
912 outer_boundary_radius -= filament_radius*1.1
913 if inner_boundary_radius != 0:
914 inner_boundary_radius += filament_radius*1.1
916 outer_inner_boundary_ratio = outer_boundary_radius / inner_boundary_radius if inner_boundary_radius != 0 else 1e10
917 if N_points % 6 == 1 and inner_boundary_radius == 0:
918 N_points -= 1
919 points = [[[0,0,0]]]
920 else:
921 points = []
923 groups_of_rotational_symmetry_first_points = create_first_hexant_points(N_points//6, outer_inner_boundary_ratio) # The first hexant of points, sorted by distance from the origin
925 R = outer_boundary_radius / point_relative_magnitude(groups_of_rotational_symmetry_first_points[-1]) # Scaling factor to place points at the correct distance from the origin
926 cart_to_hex_transformation = np.array([[R, R*np.cos(np.pi/3)], [0, R*np.sin(np.pi/3)]]) # Transformation matrix from cartesian to hexagonal coordinates
928 # Rotate the points to create the other hexants of points
929 for point in groups_of_rotational_symmetry_first_points:
930 transformed_point = np.matmul(cart_to_hex_transformation, point) # Transform the point from the cartesian coordinate system to the hexagonal coordinate system
931 point_group = []
932 rotation_angle = 0
933 for hexagon_side in range(0, 6):
934 rotation_matrix = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)], [np.sin(rotation_angle), np.cos(rotation_angle)]])
936 rotated_point = np.matmul(rotation_matrix, transformed_point)
937 point_group.append(list(rotated_point)+[0])
939 rotation_angle += np.pi/3
940 points.append(point_group)
941 return points
943 def create_geometry(self, filament_radius : float, hexagonal_filaments : bool, N_filaments : int, matrix_inner_radius : float, matrix_middle_radius : float, matrix_outer_radius : float, circular_filament_distribution : bool = False):
944 """
945 Creates the full geometry of the strand cross-section.
947 This method generates the geometry by creating a grid of points that represent the center points of the filaments.
948 It then creates the filaments at these points. The filaments can be either hexagonal or circular, depending on the `hexagonal_filaments` parameter.
949 Finally, it creates the matrix that contains the filaments. The matrix is divided into an inner, middle, and outer area. The middle section contains the filaments, and the inner and outer areas are 'empty'.
951 :param filament_radius: The radius of the filaments.
952 :type filament_radius: float
953 :param hexagonal_filaments: If True, the filaments are hexagonal. If False, the filaments are circular.
954 :type hexagonal_filaments: bool
955 :param N_filaments: The number of filaments to create.
956 :type N_filaments: int
957 :param matrix_inner_radius: The radius of the inner area of the matrix.
958 :type matrix_inner_radius: float
959 :param matrix_middle_radius: The radius of the middle area of the matrix. This is where the filaments are located.
960 :type matrix_middle_radius: float
961 :param matrix_outer_radius: The radius of the outer area of the matrix.
962 :type matrix_outer_radius: float
963 :param circular_filament_distribution: If True, the filaments are distributed in a circular pattern. If False, the filaments are distributed in a hexagonal pattern. Defaults to False.
964 :type circular_filament_distribution: bool, optional
965 """
966 # 1) Create the point-grid representing the center points of the filaments
967 # 1.1) Get the point positions
968 filament_centers = self.create_filament_center_points(N_filaments, filament_radius, matrix_middle_radius, matrix_inner_radius, circular_filament_distribution)
969 # 1.2) Create the center points
970 center_points = []
971 for layer in filament_centers:
972 layer_points = []
973 for pos in layer:
974 P = Point.create_or_get(pos)
975 layer_points.append(P)
976 center_points.append(layer_points)
978 # 2) Create the filaments
979 filaments = []
980 for layer_n in range(len(center_points)):
981 layer_filaments = []
982 for point_i in range(len(center_points[layer_n])):
983 if hexagonal_filaments:
984 filament = Hexagon(center_points[layer_n][point_i], filament_radius, np.pi/6)
985 else:
986 filament = Disk(center_points[layer_n][point_i], filament_radius)
987 layer_filaments.append(filament)
988 filaments.append(layer_filaments)
989 self.filaments = filaments
991 # 3) Create the matrix
992 # The matrix will be divided into an inner-, middle- and outer area. The middle section contains the filaments and the inner and outer areas are 'empty'.
993 # No inner section will be made if the matrix inner_radius is 0.
994 if matrix_inner_radius != 0:
995 inner_section = Disk([0,0,0], matrix_inner_radius)
996 middle_section = Disk([0,0,0], matrix_middle_radius) # Middle section
997 middle_section.inner_boundary_curves.append(inner_section.boundary_curves)
998 for layer in self.filaments:
999 for filament in layer:
1000 middle_section.inner_boundary_curves.append(filament.boundary_curves)
1001 else:
1002 middle_section = Disk([0,0,0], matrix_middle_radius) # Middle section
1003 for layer in self.filaments:
1004 for filament in layer:
1005 middle_section.inner_boundary_curves.append(filament.boundary_curves)
1007 outer_section = Disk([0,0,0], matrix_outer_radius)
1008 outer_section.inner_boundary_curves.append(middle_section.boundary_curves)
1010 self.matrix = [middle_section, outer_section] if matrix_inner_radius == 0 else [inner_section, middle_section, outer_section]
1012 def add_air(self, rad : float):
1013 # The air region is defined as the region from the matrix outer boundary to the radius 'rad'. The air radius must be greater than the matrix radius.
1014 self.air = Disk([0,0,0], rad) # Initialize the air region as a disk with radius 'rad'
1016 # Find the combined outer boundary of the strand geometry, which is the inner boundary of the air region.
1017 # The outer boundary of the strand geometry is is not necessarily the outer boundary of the matrix, as the outer matrix partition may not fully contain the full strand (as with a wire-in-channel geometry).
1018 strand_outer_boundary = set(self.matrix[0].boundary_curves) # Start with the boundary of the inner matrix partition
1019 for i, matrix_partition in enumerate(self.matrix[:-1]): # Loop over the matrix partitions
1020 next_matrix_partition = self.matrix[i+1]
1022 inner_partition_boundary = set(matrix_partition.boundary_curves)
1023 next_partition_boundary = set(next_matrix_partition.boundary_curves)
1025 if inner_partition_boundary & next_partition_boundary: # If the inner and outer partition boundaries share some curves
1026 strand_outer_boundary = strand_outer_boundary ^ next_partition_boundary # Get the combined boundary of the inner and outer partition boundaries
1027 else:
1028 strand_outer_boundary = next_partition_boundary # If the inner and outer partition boundaries do not share any curves, the outer boundary is simply the boundary of the outer partition.
1030 strand_outer_boundary = Curve.get_closed_loops(list(strand_outer_boundary))[0] # Simply used to sort the curves in the outer boundary into a correct order which can be used to create a closed loop.
1032 self.air.inner_boundary_curves.append(strand_outer_boundary)
1034 def update_tags(self):
1035 """
1036 When the geometry is loaded from a .brep file, the tags of entities with dimensions lower than the highest dimension are not preserved and may change unpredictably.
1037 This function updates the tags of the points, curves and surfaces in the geometry to ensure that they are consistent with the current gmsh model.
1038 """
1039 surfaces = sum(self.filaments, []) + sum(self.filament_holes, []) + self.matrix + [self.air]
1040 Surface.update_tags(surfaces)
1042 def create_gmsh_instance(self):
1043 """
1044 Creates the gmsh instances of the geometry.
1045 """
1046 surfaces = sum(self.filaments, []) + sum(self.filament_holes, []) + self.matrix + [self.air]
1048 Surface.set_correct_boundary_orientation(surfaces)
1050 for surface in surfaces:
1051 surface.create_gmsh_instance()
1053 def add_physical_groups(self):
1054 """
1055 Creates all physical groups.
1056 """
1057 # Filaments: Add physical boundary and surface
1058 for layer_n in range(len(self.filaments)):
1059 for filament_i in range(len(self.filaments[layer_n])):
1060 self.filaments[layer_n][filament_i].add_physical_boundary(f"Boundary: Filament {filament_i} in layer {layer_n}")
1061 self.filaments[layer_n][filament_i].add_physical_surface(f"Surface: Filament {filament_i} in layer {layer_n}")
1063 self.filaments[layer_n][filament_i].physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.filaments[layer_n][filament_i].boundary_curves[0].points[0].tag])
1065 # Add physical surface for the filament holes
1066 for layer_n in range(len(self.filament_holes)):
1067 for filament_i in range(len(self.filament_holes[layer_n])):
1068 self.filament_holes[layer_n][filament_i].add_physical_surface(f"Surface: Filament hole {filament_i} in layer {layer_n}")
1070 # Matrix: Add physical boundary and surface for each partition
1071 for i, matrix_partition in enumerate(self.matrix):
1072 matrix_partition.add_physical_boundary(f"Boundary: Matrix partition {i}")
1073 matrix_partition.add_physical_surface(f"Surface: Matrix partition {i}")
1075 # Air:
1076 self.air.add_physical_boundary(f"Boundary: Air")
1077 self.air.physical_inner_boundary_tags.append(gmsh.model.addPhysicalGroup(1, [curve.tag for curve in self.air.inner_boundary_curves[0]], name = f"InnerBoundary: Air"))
1078 self.air.physical_inner_boundary_names.append(f"InnerBoundary: Air")
1079 self.air.add_physical_surface(f"Surface: Air")
1081 self.air.strand_bnd_physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.air.inner_boundary_curves[0][0].points[0].tag])
1083 def save(self, save_file):
1084 # This function saves the geometry class to a pickle file.
1085 # The geometry class is used again during meshing.
1086 with open(save_file, "wb") as geom_save_file:
1087 pickle.dump(self, geom_save_file)
1088 print(f"Geometry saved to file {save_file}")
1090 def write_geom_to_yaml(self, file_path):
1091 # This function writes the geometry to a yaml file.
1092 # 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.
1093 # Note: Only the strands are written to the yaml file. The air region is not included.
1094 Conductor = geom.Conductor() # Create a data model for the conductor data
1096 # Add Materials to the 'Solution' section of the data model
1097 # 1) Find all unique materials used in the geometry
1098 # materials = set([surface.material for surface in Surface.surfaces_registry if surface.material is not None])
1099 # # Sort the materials into two groups by their type. All materials with the same type are grouped together.
1100 # material_groups = {material_type: [material for material in materials if material.type == material_type] for material_type in set([material.type for material in materials])}
1101 # # 2) Add all unique materials to the data model, represented by a string with the material type and index in the material group
1102 # for material_type, material_group in material_groups.items():
1103 # for i, material in enumerate(material_group):
1104 # material_name = f"{material_type}_{i}"
1105 # Conductor.Solution.Materials[material_name] = material
1107 surfaces = list(set(sum(self.filaments, []) + self.matrix)) # Combine all the filaments and the matrix to get all the surfaces which should be written. The air region is not included.
1108 curves = list(set(sum([surface.boundary_curves for surface in surfaces], []))) # Extract all the boundary curves from the surfaces
1109 points = list(set(sum([curve.points for curve in curves], []))) # Extract all the points from the curves
1112 # Populate the points dictionary with coordinates of each point
1113 # for p, point in enumerate(points):
1114 # Conductor.Geometry.Points[p] = geom.Point(Coordinates=point.pos.tolist()) #{"Coordinates": point.pos.tolist()}
1115 for p, point in enumerate(Point.points_registry):
1116 if point in points:
1117 Conductor.Geometry.Points[p] = geom.Point(Coordinates=point.pos.tolist())
1119 # Populate the curves dictionary with type of each curve and indices of its points
1120 # for c, curve in enumerate(curves):
1121 # curve_points = [Point.points_registry.index(point) for point in curve.points]
1122 # Conductor.Geometry.Curves[c] = geom.Curve(
1123 # Type=curve.__class__.__name__,
1124 # Points=curve_points
1125 # )
1126 for c, curve in enumerate(Curve.curves_registry):
1127 if curve in curves:
1128 curve_points = [Point.points_registry.index(point) for point in curve.points]
1129 Conductor.Geometry.Curves[c] = geom.Curve(
1130 Type=curve.__class__.__name__,
1131 Points=curve_points
1132 )
1135 # Populate the surfaces dictionary with material, boundary curves and inner boundary curves of each surface
1136 for a, surface in enumerate(Surface.surfaces_registry):
1137 if surface in surfaces:
1138 surface_boundary_curves = [Curve.curves_registry.index(curve) for curve in surface.boundary_curves]
1139 surface_inner_boundary_curves = [[Curve.curves_registry.index(curve) for curve in inner_boundary] for inner_boundary in surface.inner_boundary_curves]
1141 if surface in self.matrix: # Add dummy values for writing the matrix surfaces to the data model
1142 surface.layer = None
1143 surface.layer_index = None
1145 elif surface in sum(self.filaments, []): # Add dummy values for writing the filament surfaces to the data model
1146 for l, layer in enumerate(self.filaments):
1147 if surface in layer:
1148 surface.layer = l
1149 surface.layer_index = layer.index(surface)
1150 break
1152 # Name the material based on its type and index in the material groups
1153 # if surface.material is None:
1154 # material_name = None
1155 # else:
1156 # material_type = surface.material.Type
1157 # material_index = material_groups[material_type].index(surface.material)
1158 # material_name = f"{material_type}_{material_index}"
1159 material_name = surface.material
1161 Conductor.Geometry.Areas[a] = geom.Area(
1162 Material=material_name,
1163 Boundary=surface_boundary_curves,
1164 InnerBoundaries=surface_inner_boundary_curves,
1165 Layer=surface.layer,
1166 LayerIndex=surface.layer_index
1167 )
1170 # Write the data model to a yaml file
1171 UFF.write_data_to_yaml(file_path, Conductor.dict())
1173 @classmethod
1174 def read_geom_from_yaml(cls, file_path):
1175 """
1176 This function loads a geometry from a yaml file and returns a TwistedStrand object.
1177 The yaml file contains all points, curves and surfaces which define the geometry.
1178 - Points are defined by their position vector and can be referenced by an integer ID.
1179 : Position [x, y, z]
1180 - Curves are defined by their type (e.g. Line, CircleArc, etc.) and the ID of the points that make up the curves. Curves are referenced by an integer ID as well.
1181 : Type ('Line', 'CircleArc', etc.)
1182 : Points ([1,2,3]), list of point-ID defining the curve.
1183 - Surfaces are defined by material, outer boundary, inner boundaries and layer and layer-index if the surface is a strand.
1184 : Material ('Cu', 'NbTi', 'Nb3Sn', etc.)
1185 : Outer boundary ([2,3,4,5...]), list of curve-ID defining the outer boundary closed loop.
1186 : Inner boundaries ([[1,2,3], [4,5,6], ... ]), list of list of curve-IDs defining closed loops.
1187 : Layer (0, 1, 2, ...), the layer of a filament. None if the surface is part of the matrix.
1188 : LayerIndex (0, 1, 2, ...), the index of the filament in the layer. None if the surface is part of the matrix.
1190 :param file_path: The full path to the yaml file.
1191 :type file_path: str
1192 :param gmsh_curve_convention: If True, the curves are created using the gmsh convention for defining curves. Determines the order of the points in the curves. Defaults to False. This is a temporary solution. In the future, the order of the points in the curves will be updated to fit the gmsh convention.
1194 """
1195 Conductor = UFF.read_data_from_yaml(file_path, geom.Conductor)
1197 # 1) Create the points
1198 for point in Conductor.Geometry.Points.values():
1199 Point.create_or_get(point.Coordinates)
1201 # 2) Create the curves
1202 for curve in Conductor.Geometry.Curves.values():
1203 curve_type = curve.Type
1204 points = [Point.points_registry[p] for p in curve.Points]
1207 c = globals()[curve_type].create_or_get(*points) # Create the curve of the specified type
1209 # TODO: To be added.
1210 # c.contact = curve.Contact
1211 # c.thickness = curve.Thickness
1212 # c.material = curve.Material
1215 # 3) Create the surfaces
1216 # TODO: area.boundary_material and boundary_thickness are not yet used.
1217 strand = cls()
1218 layers = max([area.Layer for area in Conductor.Geometry.Areas.values() if area.Layer is not None]) + 1 # The number of layers in the strand
1219 strand.filaments = [[None for i in range(6)] for j in range(layers)] # Initialize the filaments list
1220 strand.filament_holes = [[None for i in range(6)] for j in range(layers)] # Initialize the filament holes list
1222 for area_index, area_dm in Conductor.Geometry.Areas.items():
1223 boundary_curves = [Curve.curves_registry[c] for c in area_dm.Boundary]
1224 inner_boundary_curves = [[Curve.curves_registry[c] for c in inner_boundary] for inner_boundary in area_dm.InnerBoundaries]
1225 surface = Surface(boundary_curves, inner_boundary_curves)
1226 if area_dm.Material:
1227 surface.material = area_dm.Material
1231 if area_dm.Layer is None:
1232 # It is either a matrix partition or it is a hole in a filament.
1233 # We check if it is a hole in a filament by checking if the area outer boundary is in the inner boundary of a filament.
1234 is_hole = False
1235 for other_area in Conductor.Geometry.Areas.values(): # Loop over all areas to check if the area is a hole in a filament
1236 if other_area.Layer is not None: # If the other area is a filament
1237 if area_dm.Boundary in other_area.InnerBoundaries: # If the area is a hole in the other filament
1238 # boundary_curves[0].P1, boundary_curves[0].P2 = boundary_curves[0].P2, boundary_curves[0].P1 # Reverse the order of the boundary curve points to get the correct orientation
1239 layer = other_area.Layer
1240 layer_index = other_area.LayerIndex
1241 strand.filament_holes[layer][layer_index] = surface
1242 is_hole = True
1243 break
1245 if not is_hole: # If it is not a hole, it is a matrix partition
1246 strand.matrix.append(surface)
1248 else:
1249 strand.filaments[area_dm.Layer][area_dm.LayerIndex] = surface
1251 # Remove None values from the filaments list
1252 strand.filaments = [[filament for filament in layer if filament is not None] for layer in strand.filaments]
1253 strand.filament_holes = [[hole for hole in layer if hole is not None] for layer in strand.filament_holes]
1256 # Sort the matrix partitions from inner to outer based on the outermost point of the boundary curves of the partitions
1257 strand.matrix = sorted(strand.matrix, key=lambda surface: max([max([np.linalg.norm(point.pos) for point in curve.points]) for curve in surface.boundary_curves]))
1259 return strand
1265class Geometry:
1266 """
1267 Class to generate the ConductorAC Strand geometry.
1269 This class is responsible for generating the geometry of the twisted strand.
1270 It can either load a geometry from a YAML file or create the model from scratch.
1271 The geometry is saved to a .brep file and the geometry class is saved as a pickle file. If specified, the geometry representation can also be saved to a YAML file.
1273 :ivar fdm: The fiqus inputs data model.
1274 :vartype fdm: object
1275 :ivar cacdm: The magnet section of the fiqus inputs data model.
1276 :vartype cacdm: object
1277 :ivar inputs_folder: The full path to the folder with input files, i.e., conductor and STEP files.
1278 :vartype inputs_folder: str
1279 :ivar geom_folder: The full path to the current working directory.
1280 :vartype geom_folder: str
1281 :ivar magnet_name: The name of the magnet.
1282 :vartype magnet_name: str
1283 :ivar geom_file: The full path to the .brep file where the geometry will be saved.
1284 :vartype geom_file: str
1285 :ivar verbose: If True, more information is printed in the Python console.
1286 :vartype verbose: bool
1287 :ivar gu: An instance of the GmshUtils class.
1288 :vartype gu: object
1289 """
1291 def __init__(self, fdm, inputs_folder_path, verbose=True):
1292 """
1293 Initializes the Geometry class.
1295 :param fdm: The fiqus data model.
1296 :type fdm: object
1297 :param inputs_folder_path: The full path to the folder with input files, i.e., conductor and STEP files.
1298 :type inputs_folder_path: str
1299 :param verbose: If True, more information is printed in the Python console. Defaults to True.
1300 :type verbose: bool, optional
1301 """
1302 self.fdm = fdm
1303 self.cacdm = fdm.magnet
1304 self.inputs_folder = inputs_folder_path
1305 self.geom_folder = os.path.join(os.getcwd())
1306 self.magnet_name = fdm.general.magnet_name
1307 self.geom_file = os.path.join(self.geom_folder, f'{self.magnet_name}.brep')
1308 self.verbose = verbose
1309 self.gu = GmshUtils(self.geom_folder, self.verbose)
1310 self.gu.initialize()
1312 # To see the surfaces in a better way in GUI:
1313 gmsh.option.setNumber("Geometry.SurfaceType", 2)
1315 def generate_strand_geometry(self, gui=False):
1316 """
1317 Generates the geometry of a strand based on the settings specified in the input data model.
1318 The geometry can either be loaded from a YAML file or created from scratch. After creation, the geometry
1319 can be saved to a YAML file. The method also creates gmsh instances, adds physical groups to the geometry,
1320 writes the geometry to a .brep file, and saves the geometry-class to a pickle file. If `gui` is True, it
1321 launches an interactive GUI.
1323 :param gui: If True, launches an interactive GUI after generating the geometry. Default is False.
1324 :type gui: bool, optional
1326 :return: None
1327 """
1328 print("Generating geometry")
1329 #0) Clear the registries. Used when generating reference files for tests.
1330 if Point.points_registry: # If the points registry is not empty, clear it.
1331 Point.points_registry.clear()
1332 if Curve.curves_registry: # If the curves registry is not empty, clear it.
1333 Curve.curves_registry.clear()
1334 if Surface.surfaces_registry: # If the surfaces registry is not empty, clear it.
1335 Surface.surfaces_registry.clear()
1336 if Surface.curve_loop_registry: # If the curve loop registry is not empty, clear it.
1337 Surface.curve_loop_registry.clear()
1339 # 1) Either load the geometry from a yaml file or create the model from scratch
1340 if self.cacdm.geometry.io_settings.load.load_from_yaml:
1341 CAC = TwistedStrand.read_geom_from_yaml(os.path.join(self.inputs_folder, self.cacdm.geometry.io_settings.load.filename))
1342 else:
1343 CAC = TwistedStrand()
1344 conductor_name = self.cacdm.solve.conductor_name
1345 CAC.create_geometry(
1346 self.fdm.conductors[conductor_name].strand.filament_diameter/2,
1347 self.cacdm.geometry.hexagonal_filaments,
1348 self.fdm.conductors[conductor_name].strand.number_of_filaments,
1349 self.fdm.conductors[conductor_name].strand.diameter_core/2,
1350 self.fdm.conductors[conductor_name].strand.diameter_filamentary/2,
1351 self.fdm.conductors[conductor_name].strand.diameter/2,
1352 self.cacdm.geometry.filament_circular_distribution
1353 )
1354 CAC.add_air(self.cacdm.geometry.air_radius)
1355 CAC.create_gmsh_instance()
1357 # 2) Save the geometry to a yaml file if specified
1358 if self.cacdm.geometry.io_settings.save.save_to_yaml:
1359 filename = self.cacdm.geometry.io_settings.save.filename
1360 CAC.write_geom_to_yaml(os.path.join(self.geom_folder, filename))
1363 gmsh.model.occ.synchronize()
1365 CAC.add_physical_groups() # Add physical groups to the geometry
1366 print("Writing geometry")
1367 gmsh.write(self.geom_file) # Write the geometry to a .brep file
1369 CAC.save(os.path.join(self.geom_folder, f'{self.magnet_name}.pkl')) # Save the geometry-class to a pickle file
1371 if gui:
1372 self.gu.launch_interactive_GUI()
1374 def load_conductor_geometry(self, gui=False):
1375 """
1376 Loads geometry from .brep file.
1377 """
1379 print("Loading geometry")
1381 gmsh.clear()
1382 gmsh.model.occ.importShapes(self.geom_file, format="brep")
1383 gmsh.model.occ.synchronize()
1385 if gui:
1386 self.gu.launch_interactive_GUI()