Coverage for fiqus/geom_generators/GeometryConductorAC_Strand_RutherfordCopy.py: 46%
791 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-12-19 01:48 +0000
1import os
2import pickle
3import logging
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)
12from scipy.integrate import quad
14logger = logging.getLogger('FiQuS')
18### HELPER FUNCTIONS ###
20def cylindrical_to_cartesian(rad, angle, height):
21 """
22 Convert cylindrical coordinates to Cartesian coordinates.
24 :return: A list of Cartesian coordinates [x, y, z].
25 :rtype: list
26 """
27 return [rad*np.cos(angle), rad*np.sin(angle), height]
29def cartesian_to_cylindrical(x, y, z):
30 """
31 Convert Cartesian coordinates to cylindrical coordinates.
33 :return: A list of cylindrical coordinates [rad, angle, z].
34 :rtype: list
35 """
36 rad = np.sqrt(x**2 + y**2)
37 angle = np.arctan2(y, x)
38 return [rad, angle, z]
40def rotate_vector(vector, angle):
41 """
42 Rotate a 3D vector in the xy-plane by a given angle.
44 :param vector: The 3D vector to rotate. It should be a list or numpy array of three numbers.
45 :type vector: list or numpy.ndarray
46 :param angle: The angle by which to rotate the vector, in radians.
47 :type angle: float
49 :return: The rotated vector.
50 :rtype: numpy.ndarray
51 """
52 rotation_matrix = np.array([[np.cos(angle), -np.sin(angle), 0], [np.sin(angle), np.cos(angle), 0], [0, 0, 1]])
53 return np.matmul(rotation_matrix, vector)
55### END HELPER FUNCTIONS ###
58### GEOMETRY ###
59class Point:
60 """
61 A class to represent a point in 3D space.
63 :cvar points_registry: A list of all points created. This is a class-level attribute.
64 :vartype points_registry: list
65 :cvar point_snap_tolerance: The tolerance for snapping points to existing points. This is a class-level attribute.
66 :vartype point_snap_tolerance: float
68 :ivar pos: A numpy array representing the position of the point in 3D space. This is an instance-level attribute.
69 :vartype pos: numpy.ndarray
70 :ivar tag: A gmsh tag for the point, default is None. This is an instance-level attribute.
71 :vartype tag: int, optional
72 """
73 points_registry = [] # Global registry of points
74 point_snap_tolerance = 1e-10 # Tolerance for snapping points to existing points
75 def __init__(self, pos : List[float]):
76 """
77 Constructs the attributes for the point object.
79 :param pos: A list of three floats representing the position of the point in 3D space.
80 :type pos: List[float]
81 """
82 self.pos = np.array(pos)
83 self.tag = None
85 @classmethod
86 def create_or_get(cls, pos : any):
87 """
88 Creates a new point if no point with the given position exists, otherwise returns the existing one to assure unique points.
90 :param pos: The position of the point.
91 :type pos: list
93 :return: A point object.
94 :rtype: Point
95 """
96 new_point = cls(pos)
97 if new_point in cls.points_registry:
98 return cls.points_registry[cls.points_registry.index(new_point)]
100 cls.points_registry.append(new_point)
101 return new_point
103 def create_gmsh_instance(self, meshSize : float = 0, tag : int = -1):
104 if self.tag == None:
105 self.tag = gmsh.model.occ.addPoint(self.pos[0], self.pos[1], self.pos[2], meshSize, tag)
106 else:
107 pass
109 def __repr__(self) -> str:
110 return f"Point({self.pos})"
112 def __eq__(self, o: object) -> bool:
113 # Two points sharing the same position (within a margin of point_snap_tolerance) should be considered the same point.
114 if isinstance(o, Point):
115 return np.linalg.norm(np.array(self.pos) - np.array(o.pos)) < self.point_snap_tolerance
116 return False
118 def __hash__(self) -> int:
119 return hash(tuple(self.pos))
121### CURVE ELEMENTS ###
123class Curve(ABC):
124 """
125 Abstract base class for curves in 3D space.
127 :cvar curves_registry: A list of all curves created. This is a class-level attribute.
128 :vartype curves_registry: list
130 :ivar P1: The start point of the curve. This is an instance-level attribute.
131 :vartype P1: Point
132 :ivar P2: The end point of the curve. This is an instance-level attribute.
133 :vartype P2: Point
134 :ivar points: A list of points used by the curve. This is an instance-level attribute.
135 :vartype points: list
136 :ivar tag: A tag for the curve. This is an instance-level attribute.
137 :vartype tag: int
138 """
139 curves_registry = []
140 def __init__(self) -> None:
141 """
142 Constructs the attributes for the curve object.
143 """
144 self.P1 : Point = None
145 self.P2 : Point = None
146 self.points = []
147 self.tag = None
149 # Curve.curves_registry.append(self)
151 def set_points(self, *points : Point):
152 """
153 Sets the points used by the curve.
155 :param points: The points to set.
156 :type points: Point
157 """
158 self.points = list(points)
160 @abstractmethod
161 def get_length(self):
162 pass
164 @abstractmethod
165 def create_gmsh_instance(self):
166 pass
168 @classmethod
169 @abstractmethod
170 def create_or_get(cls, *points):
171 pass
173 @classmethod
174 def get_curve_from_tag(cls, tag):
175 """
176 Returns the curve with the given tag.
178 :param tag: The tag of the curve.
179 :type tag: int
181 :return: The curve with the given tag, or None if no such curve exists.
182 :rtype: Curve or None
183 """
184 for curve in cls.curves_registry:
185 if curve.tag == tag:
186 return curve
187 return None
190 @classmethod
191 def get_closed_loops(cls, curves):
192 """
193 Returns a list of lists of curves, where each list of curves defines a closed loop.
194 The function assumes that no two closed loops share a curve and that chains of curves do not split into multiple chains.
196 :param curves: A list of curves.
197 :type curves: list
199 :return: A list of lists of curves, where each list of curves defines a closed loop.
200 :rtype: list
201 """
202 closed_loops = []
204 def get_curve_link(current_curve, curves, curve_link):
205 # Recursive function to find a link of connected curves.
206 # Curves is a list of curves to search through, not containing the current curve.
207 # Curve_link is a list of curves that are connected to the current curve.
208 # Curves are removed from 'curves' as they are added to the curve_link.
209 # The function returns when no more curves are connected to the current link.
210 for curve in curves:
211 for point in [current_curve.P1, current_curve.P2]:
212 if point in [curve.P1, curve.P2]:
213 curve_link.append(curve)
214 curves.remove(curve)
215 return get_curve_link(curve, curves, curve_link)
216 return curves, curve_link # Return the remaining curves and the curve link
218 while len(curves) > 0:
219 curve0 = curves[0]
220 curves.remove(curve0)
221 curves, curve_link = get_curve_link(curve0, curves, [curve0])
223 first_curve_points = set([curve_link[0].P1, curve_link[0].P2])
224 last_curve_points = set([curve_link[-1].P1, curve_link[-1].P2])
226 # TODO: Have to add the case where there is only a single curve in the link. This should not be considered a closed loop.
227 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.
228 closed_loops.append(curve_link)
229 # If the link only contains two curves, the curves must be connected at both ends to form a closed loop.
230 elif len(curve_link) == 2 and len(first_curve_points&last_curve_points) == 2:
231 closed_loops.append(curve_link)
233 return closed_loops
235class CircleArc(Curve):
236 """
237 A class to represent a circular arc, subclass of Curve.
239 :ivar P1: The first point of the arc. This is an instance-level attribute.
240 :vartype P1: Point
241 :ivar P2: The second point of the arc. This is an instance-level attribute.
242 :vartype P2: Point
243 :ivar C: The center of the arc. This is an instance-level attribute.
244 :vartype C: Point
245 """
246 def __init__(self, P1 : any, C : any, P2 : any) -> None:
247 """
248 Constructs the attributes for the arc object.
250 :param P1: The start point of the arc.
251 :type P1: list or Point
252 :param P2: The end point of the arc.
253 :type P2: list or Point
254 :param C: The center of the arc.
255 :type C: list or Point
256 """
257 super().__init__()
258 self.P1 = P1 if isinstance(P1, Point) else Point.create_or_get(P1)
259 self.P2 = P2 if isinstance(P2, Point) else Point.create_or_get(P2)
260 self.C = C if isinstance(C, Point) else Point.create_or_get(C)
262 self.set_points(self.P1, self.C, self.P2)
264 Curve.curves_registry.append(self)
266 def create_gmsh_instance(self, tag : int = -1):
267 if self.tag == None: # If the curve has not been created yet
268 self.P1.create_gmsh_instance()
269 self.P2.create_gmsh_instance()
270 self.C.create_gmsh_instance()
271 self.tag = gmsh.model.occ.addCircleArc(self.P1.tag, self.C.tag, self.P2.tag, tag)
273 def get_length(self):
274 radius = np.linalg.norm(self.P1.pos - self.C.pos)
275 angle = np.arccos(np.dot(self.P1.pos - self.C.pos, self.P2.pos - self.C.pos)/(radius**2))
276 return radius*angle
278 @classmethod
279 def create_or_get(cls, P1 : any, C : any, P2 : any):
280 return cls(P1, C, P2)
282 def __repr__(self) -> str:
283 return f"CircleArc({self.P1.pos}, {self.P2.pos})"
285class Line(Curve):
286 """
287 A class to represent a line, subclass of Curve.
289 :ivar P1: the start point of the line
290 :type P1: list or Point
291 :ivar P2: the end point of the line
292 :type P2: list or Point
294 """
295 def __init__(self, P1 : any, P2 : any) -> None:
296 """
297 Constructs the attributes for the line object.
298 P1 and P2 can be either 'Point' objects or lists of three floats representing the position of the point.
300 :param P1: the start point of the line
301 :type P1: list or Point
302 :param P2: the end point of the line
303 :type P2: list or Point
304 """
305 super().__init__()
306 self.P1 = P1 if isinstance(P1, Point) else Point.create_or_get(P1)
307 self.P2 = P2 if isinstance(P2, Point) else Point.create_or_get(P2)
309 self.set_points(self.P1, self.P2)
311 def get_length(self):
312 return np.linalg.norm(self.P1.pos - self.P2.pos)
314 def get_furthest_point(self):
315 """returns the line point which is further away from the origin"""
316 if np.linalg.norm(self.P1.pos) > np.linalg.norm(self.P2.pos):
317 return self.P1
318 elif np.linalg.norm(self.P1.pos) < np.linalg.norm(self.P2.pos):
319 return self.P2
320 else:
321 ValueError("cant determine furthest point. Dictance is equal")
323 def __repr__(self) -> str:
324 return f"Line from {tuple(self.P1.pos)} to {tuple(self.P2.pos)}"
326 @classmethod
327 def create_or_get(cls, P1 : any, P2 : any):
328 """
329 Creates a new line if it doesn't exist, or returns the existing one.
331 :param P1: the first point of the line
332 :type P1: Point
333 :param P2: the second point of the line
334 :type P2: Point
336 :return: a line object
337 :rtype: Line
338 """
340 new_line = cls(P1, P2)
341 if new_line in cls.curves_registry:
342 return cls.curves_registry[cls.curves_registry.index(new_line)]
344 cls.curves_registry.append(new_line)
345 return new_line
347 @classmethod
348 def is_colinear(cls, line1, line2):
349 """
350 Checks if two lines are colinear.
351 """
352 l1 = line1.P2.pos - line1.P1.pos
353 l2 = line2.P2.pos - line2.P1.pos
354 # Check if the lines are parallel
355 if np.linalg.norm(np.cross(l1, l2)) < 1e-10:
356 # Check if the lines are colinear
357 if np.linalg.norm(np.cross(l1, line2.P1.pos - line1.P1.pos)) < 1e-10:
358 return True
359 return False
361 @classmethod
362 def remove_from_registry(cls, lines):
363 """
364 Removes a list of lines from the registry.
366 :param lines: A list of lines to remove.
367 :type lines: list
368 """
369 for line in lines:
370 if line in cls.curves_registry:
371 cls.curves_registry.remove(line)
373 def create_gmsh_instance(self, tag : int = -1):
374 if self.tag == None:
375 self.P1.create_gmsh_instance()
376 self.P2.create_gmsh_instance()
377 self.tag = gmsh.model.occ.addLine(self.P1.tag, self.P2.tag, tag)
379 def __eq__(self, o: object) -> bool:
380 # Two Line-entities specified by the same two points should be considered equal lines.
381 # Check if the other object is a line
382 if isinstance(o, Line):
383 # Check if the lines have the same points
384 return (self.P1 == o.P1 and self.P2 == o.P2) or (self.P1 == o.P2 and self.P2 == o.P1)
385 return False
387 def __hash__(self) -> int:
388 return hash(frozenset([self.P1, self.P2])) # Frozenset is hashable and order-independent
390class EllipseArc(Curve):
391 """
392 A class to represent an elliptical arc, subclass of Curve.
394 :ivar P1: The start point of the arc. This is an instance-level attribute.
395 :vartype P1: Point
396 :ivar P2: The end point of the arc. This is an instance-level attribute.
397 :vartype P2: Point
398 :ivar C: The center point of the arc. This is an instance-level attribute.
399 :vartype C: Point
400 :ivar M: The major axis point of the arc (point anywhere on the major axis). This is an instance-level attribute.
401 :vartype M: Point
402 """
403 def __init__(self, P_start : any, P_center : any, P_major : any, P_end : any) -> None:
404 """
405 Initializes an elliptical arc object.
407 :param P_start: The start point of the arc. If not a Point instance, attempts to retrieve or create one.
408 :type P_start: any
409 :param P_end: The end point of the arc. If not a Point instance, attempts to retrieve or create one.
410 :type P_end: any
411 :param P_center: The center point of the arc. If not a Point instance, attempts to retrieve or create one.
412 :type P_center: any
413 :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.
414 :type P_major: any
416 :rtype: None
417 """
418 super().__init__()
419 self.P1 = P_start if isinstance(P_start, Point) else Point.create_or_get(P_start) # Start point
420 self.P2 = P_end if isinstance(P_end, Point) else Point.create_or_get(P_end) # End point
421 self.C = P_center if isinstance(P_center, Point) else Point.create_or_get(P_center) # Center point
422 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)
424 self.set_points(self.P1, self.C, self.M, self.P2)
426 Curve.curves_registry.append(self)
427 # self.tag = None
429 def get_length(self):
430 # Approximate the length of the elliptical arc
431 # 1) Center the ellipse on the origin and rotate it so that the major axis is on the x-axis
432 # a) center the ellipse
433 P1 = self.P1.pos - self.C.pos
434 P2 = self.P2.pos - self.C.pos
435 M = self.M.pos - self.C.pos
436 # b) rotate the ellipse
437 angle = np.arctan2(M[1], M[0])
438 P1 = rotate_vector(P1, -angle)
439 P2 = rotate_vector(P2, -angle)
440 # c) calculate the semi-major and semi-minor axes
441 x1, y1 = P1[0], P1[1]
442 x2, y2 = P2[0], P2[1]
444 b = np.sqrt( (x2**2 * y1**2 - x1**2 * y2**2) / (x2**2 - x1**2) ) # semi-minor axis
445 a = np.sqrt( (x2**2 * y1**2 - x1**2 * y2**2) / (y1**2 - y2**2) ) # semi-major axis
447 # 2) Calculate the length of the elliptical arc
448 theta1 = np.arctan2(y1, x1) # angle of the start point
449 theta2 = np.arctan2(y2, x2) # angle of the end point
450 if theta2 < theta1:
451 theta1, theta2 = theta2, theta1
452 t = lambda theta: np.arctan( (a/b) * np.tan(theta) ) # Change of parameter from angle to t
453 arc_length = quad(lambda theta: np.sqrt(a**2 * np.sin(t(theta))**2 + b**2 * np.cos(t(theta))**2), theta1, theta2)[0] # Calculate the arc length
455 return arc_length
458 def create_gmsh_instance(self, tag : int = -1):
459 if self.tag == None:
460 self.P1.create_gmsh_instance()
461 self.P2.create_gmsh_instance()
462 self.C.create_gmsh_instance()
463 self.M.create_gmsh_instance()
464 self.tag = gmsh.model.occ.addEllipseArc(self.P1.tag, self.C.tag, self.M.tag, self.P2.tag, tag)
466 def __repr__(self) -> str:
467 return f"EllipseArc({self.P1.pos}, {self.P2.pos})"
469 @classmethod
470 def create_or_get(cls, P_start : any, P_center : any, P_major : any, P_end : any):
471 return cls(P_start, P_center, P_major, P_end)
474### SURFACE ELEMENTS ###
476class Surface:
477 """
478 A class to represent a surface in 3D space.
480 :cvar surfaces_registry: A class-level attribute that keeps track of all created surfaces.
481 :vartype surfaces_registry: list
482 :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.
483 :vartype curve_loop_registry: list
485 :ivar boundary_curves: A list of Curve objects that define the closed outer boundary of the surface. This is an instance-level attribute.
486 :vartype boundary_curves: list
487 :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.
488 :vartype inner_boundary_curves: list
489 :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.
490 :vartype curve_loop_tag: int
491 :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.
492 :vartype inner_curve_loop_tags: list
493 :ivar surface_tag: A unique identifier for the surface, initially set to None. This is an instance-level attribute.
494 :vartype surface_tag: int
495 :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.
496 :vartype physical_boundary_tag: int
497 :ivar physical_boundary_name: A name for the physical group of the boundary curves, initially set to None. This is an instance-level attribute.
498 :vartype physical_boundary_name: str
499 :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.
500 :vartype physical_inner_boundary_tags: list
501 :ivar physical_inner_boundary_names: A list of names for the physical groups of the inner boundary curves. This is an instance-level attribute.
502 :vartype physical_inner_boundary_names: list
503 :ivar physical_surface_tag: A unique identifier for the physical group of the surface, initially set to None. This is an instance-level attribute.
504 :vartype physical_surface_tag: int
505 :ivar physical_surface_name: A name for the physical group of the surface, initially set to None. This is an instance-level attribute.
506 :vartype physical_surface_name: str
507 :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.
508 :vartype material: any
509 """
511 surfaces_registry = []
512 curve_loop_registry = []
513 def __init__(self, boundary_curves : List[Curve] = [], inner_boundary_curves : List[List[Curve]] = []):
514 """
515 Constructs the attributes for the surface object.
517 :param boundary_curves: A list of Curve objects that define the outer boundary of the surface. The curves must form a closed loop.
518 :type boundary_curves: list[Curve]
519 :param inner_boundary_curves: A list of lists of Curve objects. Each list of curves defines the boundary of a hole in the surface.
520 :type inner_boundary_curves: list[list[Curve]]
521 """
522 self.boundary_curves = boundary_curves
523 self.inner_boundary_curves = inner_boundary_curves
525 self.curve_loop_tag = None
526 self.inner_curve_loop_tags : List[int] = []
527 self.surface_tag = None
529 self.physical_boundary_tag = None
530 self.physical_boundary_name = None
531 self.physical_inner_boundary_tags = []
532 self.physical_inner_boundary_names = []
534 self.physical_surface_tag = None
535 self.physical_surface_name = None
537 self.material = None
539 Surface.surfaces_registry.append(self)
541 def create_gmsh_instance(self):
542 if self.surface_tag == None:
543 for curve in self.boundary_curves + sum(self.inner_boundary_curves, []):
544 curve.create_gmsh_instance()
546 self.curve_loop_tag = self.create_or_get_curve_loop(self.boundary_curves)
547 self.inner_curve_loop_tags = [self.create_or_get_curve_loop(curves) for curves in self.inner_boundary_curves]
548 self.surface_tag = gmsh.model.occ.add_plane_surface([self.curve_loop_tag] + self.inner_curve_loop_tags)
550 def add_physical_boundary(self, name : str = ''):
551 self.physical_boundary_tag = gmsh.model.add_physical_group(1, [curve.tag for curve in self.boundary_curves], name = name)
552 self.physical_boundary_name = name
554 def add_physical_surface(self, name : str = ''):
555 self.physical_surface_tag = gmsh.model.add_physical_group(2, [self.surface_tag], name = name)
556 self.physical_surface_name = name
558 def get_circumference(self):
559 return sum([curve.get_length() for curve in self.boundary_curves])
562 @classmethod
563 def update_tags(cls, surfaces):
564 """
565 Updates the tags of model entities of dimension lower than 2.
567 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.
569 Steps:
570 1) Find the tags of all outer boundary curves of every surface.
571 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.
572 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.
573 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.
575 :param surfaces: A list of Surface objects to update the tags for.
576 :type surfaces: list[Surface]
577 """
579 gmsh.model.occ.synchronize()
580 #1) Find the tags of all outer boundary curves of every surface
581 #1.1) Find the inner surfaces of each surface with holes
582 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.
583 for surface in surfaces:
584 inner_surface_indices = []
585 if surface.inner_boundary_curves: # If the surface has holes
586 for inner_boundary in surface.inner_boundary_curves:
587 for surface in surfaces: # Loop through all surfaces to find the surfaces that are inside the outer surface.
588 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.
589 inner_surface_indices.append(surfaces.index(surface)) # Add the index of the inner surface to the list of inner surfaces.
590 surfaces_inner_surface_indices.append(inner_surface_indices)
592 #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.
593 surfaces_outer_boundary_tags = [set([abs(tag) for dim, tag in gmsh.model.get_boundary([(2, surface.surface_tag)])]) for surface in surfaces]
594 for surface_i in range(len(surfaces)):
595 if surfaces[surface_i].inner_boundary_curves:
596 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
597 surfaces_outer_boundary_tags[surface_i] = surfaces_outer_boundary_tags[surface_i] - surface_inner_surfaces_boundary_tags
599 #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.
600 # We also find the tags of the curves in each group.
601 surface_outer_boundary_curves = [set(surface.boundary_curves) for surface in surfaces]
603 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.
604 curve_group_tags = [] # List of sets of tags. Each set of tags corresponds to the set of curves in curve_groups.
605 for i, si_curves in enumerate(surface_outer_boundary_curves):
606 for j, sj_curves in enumerate(surface_outer_boundary_curves[i+1:]):
607 j += i+1
608 common_curves = si_curves & sj_curves
609 if common_curves:
610 curve_groups.append(common_curves) # Add the common curves to the list of curve groups
611 si_curves -= common_curves # Remove the common curves from the surface boundary curves
612 sj_curves -= common_curves # Remove the common curves from the surface boundary curves
614 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
615 surfaces_outer_boundary_tags[i] -= curve_group_tags[-1] # Remove the tags of the common curves from the surface boundary tags
616 surfaces_outer_boundary_tags[j] -= curve_group_tags[-1] # Remove the tags of the common curves from the surface boundary tags
618 curve_groups.append(si_curves) # Add the remaining curves to the list of curve groups
619 curve_group_tags.append(surfaces_outer_boundary_tags[i]) # Add the tags of the remaining curves to the list of curve group tags
621 #3) Update the tags of the curves in the curve groups
622 # The tags are assigned not to their corresponding curve but to any curve in the group.
623 for group, group_tags in zip(curve_groups, curve_group_tags):
624 for curve, tag in zip(group, group_tags):
625 curve.tag = tag
627 #4) We have now updated the tags of all curves. Next we update the tags of the points.
629 #4.1) Get all points in each group of curves and the tags of the points in each group of curves
630 curve_groups_points = [set([point for curve in group for point in [curve.P1, curve.P2]]) for group in curve_groups]
631 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]
633 #4.2) Points which are in multiple curve-groups are assigned to a new group.
634 # These points will be removed from the groups they are in and assigned to the new group, based on the groups they are in.
635 # Iterate trough all points and check which groups they are in. Points in same groups will be assigned to a new group.
636 all_points = set.union(*curve_groups_points)
637 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.
638 for point in all_points:
639 groups_point_is_in = [i for i, group in enumerate(curve_groups_points) if point in group]
640 # If the point is in multiple groups, remove the point from all groups
641 if len(groups_point_is_in) > 1:
642 for i in groups_point_is_in:
643 curve_groups_points[i].remove(point) # Remove the point from all groups, as it will be assigned to a new group.
644 # 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.
645 point_new_groups[tuple(sorted(groups_point_is_in))] = point_new_groups.get(tuple(sorted(groups_point_is_in)), []) + [point]
647 #4.3) Update the tags of the points in the new groups
648 # Get the tags of all points in each group of points as the boundary of the group of curves
649 for group_indices, points in point_new_groups.items():
650 # 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.
651 point_tags = set.intersection(*[curve_groups_point_tags[i] for i in group_indices])
653 # Update the tags of the points in the group
654 for point, point_tag in zip(points, point_tags):
655 point.tag = point_tag
657 # 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.
658 for i in group_indices:
659 curve_groups_point_tags[i] -= point_tags
661 #4.4) Update the tags of points in the remaining groups 'curve_groups_points'
662 for group, group_tags in zip(curve_groups_points, curve_groups_point_tags):
663 for point, tag in zip(group, group_tags):
664 point.tag = tag
669 @classmethod
670 def remove_from_registry(cls, surfaces):
671 for surface in surfaces:
672 cls.surfaces_registry.remove(surface)
674 @classmethod
675 def create_or_get_curve_loop(cls, curves):
676 """
677 Creates a curve loop if it does not already exist, otherwise returns the existing curve loop.
679 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.
680 If it does, the method returns the tag of the existing curve loop.
681 If it doesn't, the method creates a new curve loop, adds it to the curve_loop_registry, and returns its tag.
683 :param curves: A list of Curve objects that define the curve loop.
684 :type curves: list[Curve]
686 :return: The tag of the curve loop.
687 :rtype: int
688 """
689 for curve_loop in cls.curve_loop_registry:
690 if curve_loop['curves'] == set(curves):
691 return curve_loop['tag']
693 tag = gmsh.model.occ.addCurveLoop([curve.tag for curve in curves])
694 cls.curve_loop_registry.append({'tag': tag, 'curves': set(curves)})
695 return tag
697 @classmethod
698 def replace_overlapping_edges(cls):
699 """
700 Replaces overlapping boundary-curves of type 'Line' across all surfaces.
702 If multiple surface boundaries contain lines which are overlapping, the overlapping lines are replaced by the fraction
703 of the line which is unique to the surface as well as the fraction of the line which is shared with the other surface.
705 Steps:
706 1) Sort all existing lines into groups of lines which are colinear.
707 2) For each group of colinear lines, make a list of all the unique points used, sorted by their location on the colinear line.
708 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.
709 """
710 #1) Sort all existing lines into groups of lines which are colinear
711 lines = [line for line in Curve.curves_registry if type(line) == Line]
712 if len(lines) == 0: # If there are no lines, skip the rest of the function
713 return
714 colinear_groups = [[lines[0]]]
715 for line in lines[1:]:
716 for group in colinear_groups:
717 if Line.is_colinear(line, group[0]):
718 group.append(line)
719 break
720 else:
721 colinear_groups.append([line])
723 #2) For each group of colinear lines, make a list of all the unique points used, sorted by their location on the colinear line
724 colinear_groups_points = []
725 for group in colinear_groups:
726 points = []
727 for line in group:
728 points.append(line.P1)
729 points.append(line.P2)
730 points = list(set(points)) # Remove duplicates
731 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
732 positions = [p.pos - points[0].pos for p in points] # Move the points so that the colinear line passes through the origin
733 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.
734 points = np.array(points)[np.argsort(positions)].tolist() # Sort the points by their position along the colinear line
735 colinear_groups_points.append(points)
737 #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
738 for area in cls.surfaces_registry: # For each rectangle
739 for l, line in enumerate(area.boundary_curves): # For each line
740 for i, group in enumerate(colinear_groups): # Find the group of colinear lines that the line belongs to
741 if line in group:
742 if len(group) == 1: # If the line is not colinear with any other line, skip it
743 break
744 # Find the points that define the line fragments
745 points = colinear_groups_points[i]
746 line_point1_index = points.index(line.P1)
747 line_point2_index = points.index(line.P2)
748 if line_point1_index > line_point2_index:
749 # If the points orientation differs from the line orientation, reverse the list of points
750 # The points are sorted by their position along the colinear line, which is not necessarily the same as the orientation of the line.
751 # The orientation of the fragments must match the orientation of the line.
752 points.reverse()
753 line_point1_index = points.index(line.P1)
754 line_point2_index = points.index(line.P2)
755 line_fragment_points = points[line_point1_index:line_point2_index+1] # The points that define the line fragments
757 # Create the line fragments
758 line_fragments = [Line.create_or_get(line_fragment_points[i], line_fragment_points[i+1]) for i in range(len(line_fragment_points)-1)]
760 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.
761 Line.remove_from_registry([line]) # Remove the old line from the lines-registry.
762 for fragment in line_fragments:
763 area.boundary_curves.insert(area.boundary_curves.index(line), fragment) # Insert the new line fragments into the boundary curves of the surface.
764 area.boundary_curves.remove(line) # Remove the original line from the boundary curves of the surface.
766 @classmethod
767 def set_correct_boundary_orientation(self, surfaces):
768 """
769 When creating surfaces with holes, the boundaries must have a specific orientation to ensure that the surfaces are correctly defined.
770 This method sets the correct orientation of the boundaries of the surfaces.
771 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.
772 We therefore must ensure that the first curve of the curve-loop is oriented correctly.
773 """
774 # It seems to work to just set all the boundaries to have the same orientation. Very simple and does not require any complex logic.
775 for surface in surfaces:
776 outer_boundary_points = sum([[curve.P1, curve.P2] for curve in surface.boundary_curves], [])
777 mean_point = np.mean([point.pos for point in outer_boundary_points], axis=0) # Mean point of the boundary (Center of mass)
778 mean_to_P1 = surface.boundary_curves[0].P1.pos - mean_point # Vector from the mean point to the first point of the boundary curve
779 mean_to_P2 = surface.boundary_curves[0].P2.pos - mean_point # Vector from the mean point to the second point of the boundary curve
781 # 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.
782 if np.linalg.det(np.column_stack((mean_to_P1[:-1], mean_to_P2[:-1]))) < 0:
783 # If the determinant is negative we reverse the orientation so that it is positive (counter-clockwise orientation)
784 surface.boundary_curves[0].P1, surface.boundary_curves[0].P2 = surface.boundary_curves[0].P2, surface.boundary_curves[0].P1
786class Disk(Surface):
787 """
788 A class to represent a disk. Inherits from the Surface class.
790 :ivar rad: The radius of the disk. This is an instance-level attribute.
791 :vartype rad: float
792 :ivar partitions: The number of partitions to divide the boundary of the disk into when generating boundary curves. This is an instance-level attribute.
793 :vartype partitions: int
794 :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.
795 :vartype center_point: Point
796 :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.
797 :vartype physicalEdgePointTag: int
798 """
799 def __init__(self, center_point : any, rad : float, partitions : int = 4):
800 """
801 Constructs the attributes for the disk object.
803 :param center_point: The center point of the disk. Can be a position ([x, y, z]) or a 'Point' object.
804 :type center_point: any
805 :param rad: The radius of the disk.
806 :type rad: float
807 :param partitions: The number of partitions to divide the boundary of the disk into when generating boundary curves, defaults to 4
808 :type partitions: int, optional
809 """
810 super().__init__()
811 self.rad = rad
812 self.partitions = partitions
814 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
815 self.boundary_curves = self.generate_boundary_curves()
816 self.inner_boundary_curves = []
818 self.physicalEdgePointTag = None # Used in pro-template for fixing phi=0 in the outer matrix boundary as well as on the filament boundaries.
820 def generate_boundary_curves(self):
821 """
822 Generates the boundary curves for the disk.
824 This method divides the boundary of the disk into 'partitions' number of segments and creates a CircleArc for each segment.
826 :return: A list of CircleArc objects that define the boundary of the disk.
827 :rtype: list
828 """
829 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)]
830 boundary_curves = [CircleArc(edgePoints[n], self.center_point, edgePoints[(n+1)%self.partitions]) for n in range(self.partitions)]
831 return boundary_curves
833class Rectangle(Surface):
834 """
835 A class to represent a Rectangle surface. Inherits from the Surface class.
837 :ivar width: The width of the rectangle. This is an instance-level attribute.
838 :vartype width: float
839 :ivar height: The height of the rectangle. This is an instance-level attribute.
840 :vartype height: float
841 :ivar center_point: The center point of the Rectangle. Can be a position ([x, y, z]) or a 'Point' object. This is an instance-level attribute.
842 :vartype center_point: Point
843 """
844 def __init__(self, center_point : any, width : float, height : float):
845 """
846 Constructs the attributes for the Rectangle object.
847 """
848 super().__init__()
849 self.width = width
850 self.height = height
852 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
853 self.boundary_curves = self.generate_boundary_lines()
854 self.inner_boundary_curves = []
856 def generate_boundary_lines(self):
857 """
858 Generates the boundary lines for the Rectangle.
860 :return: A list of Line objects that define the closed loop boundary of the Rectangle.
861 :rtype: list
862 """
864 edgePoints = np.array(self.center_point.pos) + (np.array([self.width/2.0, -self.height/2.0, 0]), np.array([self.width/2.0, self.height/2.0, 0]), np.array([-self.width/2.0, self.height/2.0, 0]), np.array([-self.width/2.0, -self.height/2.0, 0]))
865 boundary_curves = [Line(edgePoints[n], edgePoints[(n+1)%4]) for n in range(4)]
866 return boundary_curves
868class Square(Surface):
869 """
870 A class to represent a Square surface. Inherits from the Surface class.
872 :ivar rad: The radius of the biggest circle fitting inside the square(can be interpreted as the smallest boundary distance to center). This is an instance-level attribute.
873 :vartype rad: float
874 :ivar partitions: The number of partitions to divide the boundary of the square into when generating the boundary curves. This is an instance-level attribute.
875 :vartype partitions: int
876 :ivar center_point: The center point of the Square. Can be a position ([x, y, z]) or a 'Point' object. This is an instance-level attribute.
877 :vartype center_point: Point
878 """
879 def __init__(self, center_point : any, rad : float, partitions : int = 4):
880 """
881 Constructs the attributes for the Square object.
883 :param center_point: The center point of the Square. Can be a position ([x, y, z]) or a 'Point' object.
884 :type center_point: any
885 :param rad: The radius of a circle fitting inside the Square.
886 :type rad: float
887 :param partitions: The number of partitions to divide the boundary of the Square into when generating the boundary curves, defaults to 4
888 :type partitions: int, optional
889 """
890 super().__init__()
891 self.rad = rad
892 self.partitions = partitions
894 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
895 self.boundary_curves = self.generate_boundary_lines()
896 self.inner_boundary_curves = []
898 def generate_boundary_lines(self):
899 """
900 Generates the boundary lines for the Square.
902 This method divides the boundary of the square into 'partitions' number of Lines.
904 :return: A list of Line objects that define the closed loop boundary of the Square.
905 :rtype: list
906 """
907 if self.partitions != 4:
908 raise ValueError(
909 f"FiQuS does not support a square air boundary with partition count: {self.partitions}!"
910 )
912 edgePoints = np.array(self.center_point.pos) + (np.array([self.rad, -self.rad, 0]), np.array([self.rad, self.rad, 0]), np.array([-self.rad, self.rad, 0]), np.array([-self.rad, -self.rad, 0]))
913 boundary_curves = [Line(edgePoints[n], edgePoints[(n+1)%(self.partitions)]) for n in range(self.partitions)]
914 return boundary_curves
916class Semicircle(Surface):
917 """
918 A class to represent a Semicircle surface (aligned on the right side of y-axis). Inherits from the Surface class.
920 :ivar offset_x: The offset of the semicircle diameter center in the negative x direction. This is an instance-level attribute.
921 :vartype offset_x: float
922 :ivar rad: The radius of the biggest circle fitting inside the square(can be interpreted as the smallest boundary distance to center). This is an instance-level attribute.
923 :vartype rad: float
924 :ivar partitions: The number of partitions to divide the boundary of the square into when generating the boundary curves. This is an instance-level attribute.
925 :vartype partitions: int
926 :ivar center_point: The center point of the Square. Can be a position ([x, y, z]) or a 'Point' object. This is an instance-level attribute.
927 :vartype center_point: Point
928 """
929 # this is used in the meshing process to force the gmsh cohomology cut on the diameter of the semicircle domain
930 physical_cohomology_subdomain = None
931 def __init__(self, center_point : any, offset_x : float, rad : float, partitions : int = 4):
932 """
933 Constructs the attributes for the Square object.
935 :param center_point: The center point of the Square. Can be a position ([x, y, z]) or a 'Point' object.
936 :type center_point: any
937 :param offset_x: The offset of the semicircle diameter center in the negative x direction. This is an instance-level attribute.
938 :type offset_x: float
939 :param rad: The radius of the Semicircle.
940 :type rad: float
941 :param partitions: The number of partitions to divide the boundary of the Semicircle into when generating the boundary curves, defaults to 4.
942 :type partitions: int, optional
943 """
944 super().__init__()
945 self.offset_x = offset_x
946 self.rad = rad
947 self.partitions = partitions
949 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
950 self.boundary_curves = self.generate_boundary_lines()
951 self.inner_boundary_curves = []
953 def generate_boundary_lines(self):
954 """
955 Generates the boundary lines for the Square.
957 This method divides the boundary of the semicircle into 'partitions' number of Lines.
959 :return: A list of Line objects that define the closed loop boundary of the Semicircle.
960 :rtype: list
961 """
962 if self.partitions != 4:
963 raise ValueError(
964 f"FiQuS does not support a semicircle air boundary with partition count: {self.partitions}!"
965 )
966 if self.offset_x > self.rad:
967 raise ValueError(
968 f"FiQuS does not support a semicircle with center offset bigger than radius!"
969 )
971 edgePoints = np.array(self.center_point.pos) + (np.array([-self.offset_x, -self.rad, 0]), np.array([self.rad-self.offset_x, 0, 0]), np.array([-self.offset_x, self.rad, 0]), np.array([-self.offset_x, 0, 0]))
972 boundary_curves = []
973 for n in range(self.partitions):
974 if n < (self.partitions/2):
975 boundary_curves.append(CircleArc(edgePoints[n], [-self.offset_x, 0, 0], edgePoints[n+1]))
976 else:
977 boundary_curves.append(Line(edgePoints[n], edgePoints[(n+1)%(self.partitions)]))
978 return boundary_curves
980 def add_physical_boundary(self, name : str = ''):
981 """
982 This method extends the ususal procedure defined in 'Surface'.
983 It generates an additional cohomology subdomain physical group and stores its tag as additional attribute """
984 self.physical_boundary_tag = gmsh.model.add_physical_group(1, [curve.tag for curve in self.boundary_curves], name = name)
985 self.physical_boundary_name = name
986 # additional procedure - add the arc of semicircle boundary as cohomology subdomain
987 _, lines = gmsh.model.get_adjacencies(2, self.surface_tag)
988 self.physical_cohomology_subdomain = gmsh.model.add_physical_group(1, lines[0:int(len(self.boundary_curves)/2)], name = name + ' subdomain')
990class SquareSection(Surface):
991 """
992 A class to represent a square surface intersected with the quarter of a circle (used as surface in the 'periodic_square' model geometry). Inherits from the Surface class.
994 :ivar rad: The max radius of a circle fitting inside the initial Square, which is then intersected by intersection_curve. This is an instance-level attribute.
995 :vartype rad: float
996 :ivar intersection_curve: The Curve intersecting the square(defined by exact 2 point on the square boundary). This is an instance-level attribute.
997 :vartype intersection_curve: Curve object
998 :ivar partitions: The number of partitions to divide the resulting boundary of the SquareSection into when generating boundary curves. This is an instance-level attribute.
999 :vartype partitions: int
1000 :ivar center_point: The center point based on which the square section is extruded outwards. Can be a position ([x, y, z]) or a 'Point' object. This is an instance-level attribute.
1001 :vartype center_point: Point
1002 """
1003 def __init__(self, center_point : any, rad : float, intersection_curve : Curve, partitions : int = 5):
1004 """
1005 Constructs the attributes for the SquareSection object.
1007 :param center_point: The center point of extrusion for the SquareSection. Can be a position ([x, y, z]) or a 'Point' object.
1008 :type center_point: any
1009 :param rad: The max radius of a circle fitting inside the Square before intersection.
1010 :type rad: float
1011 :param partitions: The number of partitions to divide the boundary of the SquareSection into when generating boundary curves, should allways be 5
1012 :type partitions: int, optional
1013 """
1014 super().__init__()
1015 self.rad = rad
1016 self.intersection_curve = intersection_curve
1017 self.partitions = partitions
1019 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
1020 self.outer_boundary_curves : List[Curve] = []
1021 self.cut_curves : List[Curve] = []
1022 self.boundary_curves = self.generate_boundary_lines()
1024 def generate_boundary_lines(self):
1025 """
1026 Generates the boundary lines for the SquareSection.
1028 This method divides the boundary of the square into 'partitions' number of segments and creates a line for each segment.
1030 :return: A list of line objects that define the boundary of the IntersectedSquare.
1031 :rtype: list
1032 """
1033 if self.partitions != 5:
1034 raise ValueError(
1035 f"FiQuS does not support a SquareSections with partition count: {self.partitions}!"
1036 )
1037 if not self.center_point.pos[0] == self.center_point.pos[1] == self.center_point.pos[2] == 0:
1038 logger.info(self.center_point.pos)
1039 raise ValueError(
1040 f"FiQuS does not support a SquareSections with center extrusion point: {self.center_point}!"
1041 )
1043 # extrude intersection curves to the outside
1044 line1 = Line.create_or_get(self.intersection_curve.P2, Point.create_or_get(self.intersection_curve.P2.pos/np.linalg.norm(self.intersection_curve.P1.pos)*self.rad))
1045 line4 = Line.create_or_get(self.intersection_curve.P1, Point.create_or_get(self.intersection_curve.P1.pos/np.linalg.norm(self.intersection_curve.P1.pos)*self.rad))
1046 self.cut_curves = [line1, line4]
1047 # close IntersectedSquare of with outer boundary lines
1048 line2 = Line.create_or_get(line1.get_furthest_point(), Point.create_or_get(line1.get_furthest_point().pos+line4.get_furthest_point().pos))
1049 line3 = Line.create_or_get(line4.get_furthest_point(), line2.get_furthest_point())
1050 self.outer_boundary_curves = [line2, line3]
1052 boundary_curves = [self.intersection_curve, line1, line2, line3, line4]
1053 return boundary_curves
1056class Composite():
1057 """
1058 A helper class to handle a composite surface made up by multiple smaller surfaces.
1060 :ivar sections: A List of connected Surfaces which create a Composite regime.
1061 :vartype sections: List of Surfaces
1062 """
1063 def __init__(self, sections : List[Surface]):
1064 """
1065 Constructs the attributes for the Composite object.
1067 :type center_point: any
1068 :param sections: A List of Surface Sections which make up the Composite.
1069 :type sections: List of Surfaces
1071 """
1072 super().__init__()
1074 self.sections = sections
1076 self.physical_surface_tag = None
1077 self.physical_surface_name = None
1078 self.physical_boundary_tag = None
1079 self.physical_boundary_name = None
1080 self.strand_bnd_physicalEdgePointTag = None
1082 self.physical_inner_boundary_tags = []
1083 self.physical_inner_boundary_names = []
1085 self.physical_cuts = []
1086 self.inner_boundary_curves = []
1087 self.boundary_curves = self.generate_boundary_lines()
1089 def generate_boundary_lines(self):
1090 """
1091 Generates the boundary lines for the composite surface, made up by the section surfaces given on initialization.
1093 :return: A list of line objects that define the (outer) boundary of the composition.
1094 :rtype: list
1095 """
1096 boundary_curves = []
1097 for section in self.sections:
1098 boundary_curves.extend(section.outer_boundary_curves)
1099 return boundary_curves
1101 def add_physical_surface(self, name : str = ''):
1102 """
1103 Generates a physical surface group containing all the section surfaces.
1104 """
1105 self.physical_surface_tag = gmsh.model.add_physical_group(2, [section.surface_tag for section in self.sections], name = name)
1106 self.physical_surface_name = name
1108 @abstractmethod
1109 def add_physical_cuts(self, name):
1110 pass
1112 @abstractmethod
1113 def add_physical_boundaries(self, name):
1114 pass
1116class CompositeSquare(Composite):
1117 """
1118 A class to represent the Composite Square structure of the air surface in the 'periodic_square' model. Inherits general functionality from 'Composite'.
1119 """
1120 # special boundary tags
1121 physical_left_boundary_tag = None
1122 physical_right_boundary_tag = None
1123 physical_top_boundary_tag = None
1124 physical_bottom_boundary_tag = None
1126 def add_physical_cuts(self, name : str = ''):
1127 """
1128 Generates the two physical groups for the CompositeSquare air domain in the 'periodic_square' model.
1129 Cuts are used in getDP to imprint a source field in OmegaCC.
1130 """
1131 # generate physical cuts
1132 self.physical_cuts = []
1133 vertical_cuts_tags = [-self.sections[2].cut_curves[0].tag, self.sections[1].cut_curves[1].tag]
1134 self.physical_cuts.append(gmsh.model.add_physical_group(1, vertical_cuts_tags, name = name + " vertical"))
1135 vertical_boundary_tags = [self.sections[3].intersection_curve.tag, self.sections[0].intersection_curve.tag]
1136 self.physical_cuts.append(gmsh.model.add_physical_group(1, vertical_boundary_tags, name = name + " vertical boundary"))
1138 horizontal_cuts_tags = [-self.sections[0].cut_curves[1].tag, self.sections[1].cut_curves[0].tag]
1139 self.physical_cuts.append(gmsh.model.add_physical_group(1, horizontal_cuts_tags, name = name + " horizontal"))
1140 horizontal_boundary_tags = [self.sections[0].intersection_curve.tag, self.sections[1].intersection_curve.tag]
1141 self.physical_cuts.append(gmsh.model.add_physical_group(1, horizontal_boundary_tags, name = name + " horizontal boundary"))
1143 def add_physical_boundaries(self, name : str = ''):
1144 """
1145 Generates direction specific physical boundary groups for the the composite square and stores the group tags as additional attributes.
1146 """
1147 # add complete boundary closed loop
1148 self.physical_boundary_tag = gmsh.model.add_physical_group(1, [curve.tag for curve in self.boundary_curves], name = name)
1149 self.physical_boundary_name = name
1150 # here we actually have to get the line tags based on the sub surfaces because dim1 object-tags are not preserved in gmsh :(
1151 _, lines_sec0 = gmsh.model.get_adjacencies(2, self.sections[0].surface_tag)
1152 _, lines_sec1 = gmsh.model.get_adjacencies(2, self.sections[1].surface_tag)
1153 _, lines_sec2 = gmsh.model.get_adjacencies(2, self.sections[2].surface_tag)
1154 _, lines_sec3 = gmsh.model.get_adjacencies(2, self.sections[3].surface_tag)
1156 self.physical_left_boundary_tag = gmsh.model.add_physical_group(1, [lines_sec3[3], lines_sec0[2]], name = name + " left")
1157 self.physical_bottom_boundary_tag = gmsh.model.add_physical_group(1, [lines_sec0[3], lines_sec1[2]], name = name + " bottom")
1158 self.physical_right_boundary_tag = gmsh.model.add_physical_group(1, [lines_sec1[3], lines_sec2[2]], name = name + " right")
1159 self.physical_top_boundary_tag = gmsh.model.add_physical_group(1, [lines_sec2[3], lines_sec3[2]], name = name + " top")
1161 # set gauge point on boundary
1162 self.strand_bnd_physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.sections[2].outer_boundary_curves[0].points[1].tag])
1164class Hexagon(Surface):
1165 """
1166 A class to represent a hexagon. Inherits from the Surface class.
1168 :ivar rad: The radius of the hexagon.
1169 :vartype rad: float
1170 :ivar center_point: The center point of the hexagon. Can be a position ([x, y, z]) or a 'Point' object.
1171 :vartype center_point: Point
1172 :ivar rotation: The rotation of the hexagon in radians.
1173 :vartype rotation: float
1174 """
1175 def __init__(self, center_point : any, rad : float, rotation : float = 0):
1176 """
1177 Constructs the attributes for the hexagon object.
1179 :param center_point: The center point of the hexagon. Can be a position ([x, y, z]) or a 'Point' object.
1180 :type center_point: any
1181 :param rad: The radius of the hexagon.
1182 :type rad: float
1183 :param rotation: The rotation of the hexagon in radians.
1184 :type rotation: float
1185 """
1186 super().__init__()
1187 self.rad = rad
1188 self.center_point = center_point if isinstance(center_point, Point) else Point.create_or_get(center_point)
1189 self.rotation = rotation
1190 self.boundary_curves = self.generate_boundary_curves()
1191 self.inner_boundary_curves = []
1193 def generate_boundary_curves(self):
1194 """
1195 Generates the boundary curves for the hexagon.
1197 This method creates the boundary of the hexagon by rotating a vector from the center-point to the first edge-point.
1199 :return: A list of Line objects that define the boundary of the hexagon.
1200 :rtype: list
1201 """
1202 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)]
1203 boundary_curves = [Line.create_or_get(edgePoints[n], edgePoints[(n+1)%6]) for n in range(6)]
1204 return boundary_curves
1206### END GEOMETRY ###
1209class TwistedStrand:
1210 """
1211 A class to represent a 2D cross section of a twisted strand.
1213 :ivar filaments: Each list of surfaces represent six filaments in the same layer.
1214 :vartype filaments: list[list[Surface]]
1215 :ivar Matrix: List of surfaces corresponding to the matrix partitions.
1216 :vartype Matrix: list[Surface]
1217 :ivar Air: A surface representing the air region.
1218 :vartype Air: Surface
1219 """
1220 def __init__(self) -> None:
1221 """
1222 Initializes the TwistedStrand object.
1223 """
1224 self.filaments : List[List[Surface]] = []
1225 self.filament_holes : List[List[Surface]] = [] # If the filaments have holes (only for certain geometries obtained via YAML files)
1226 self.matrix : List[Surface] = [] # Inner, middle, outer
1227 self.air : List[Surface] = [] # one or more air surfaces
1228 self.air_composition : Composite = None # Composite structure of multiple air (if more than one air surface is defined)
1229 self.domain_cut : int = None
1231 def create_filament_center_points(self, N_points : int, filament_radius : float, outer_boundary_radius : float, inner_boundary_radius : float, circular_filament_distribution = True):
1232 """
1233 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
1234 6 points satisfying rotational symmetry when rotated by pi/3. The first layer is composed of one group, the second of two groups, etc.
1236 :param N_points: The number of points to create.
1237 :type N_points: int
1238 :param filament_radius: The radius of the filament.
1239 :type filament_radius: float
1240 :param outer_boundary_radius: The radius of the outer boundary.
1241 :type outer_boundary_radius: float
1242 :param inner_boundary_radius: The radius of the inner boundary.
1243 :type inner_boundary_radius: float
1244 :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.
1245 :type circular_filament_distribution: bool, optional
1247 :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.
1248 :rtype: list[list[list[float]]]
1249 """
1250 # The hexagonal grid can be created by transforming the cartesian coordinate grid to a hexagonal grid.
1251 # Point centers are created as unit vectors in the hexagonal grid, and then transformed to cartesian coordinates.
1252 # This function takes as input a vector in the hexagonal grid and returns its distance from the origin in the cartesian coordinate system.
1253 point_relative_magnitude = lambda v : np.sqrt(v[0]**2 + 2*v[0]*v[1]*np.cos(np.pi/3) + v[1]**2)
1254 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))):
1255 """
1256 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.
1257 The function returns the points in the first sextant of the matrix, sorted by distance from the origin.
1259 The points are created in an iterative manner, starting from the points closest to the origin and moving outwards.
1260 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
1261 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.
1262 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.
1263 """
1265 if circular_filament_distribution:
1266 # The next point is the one with the smallest distance from the origin
1267 next_point = min(possible_next_points, key=point_relative_magnitude)
1269 else:
1270 #1) Filter possible next points to include only points in the lowest not full layer.
1271 next_points_in_same_layer = possible_next_points[possible_next_points.sum(axis=1) == possible_next_points.sum(axis=1).min()]
1272 #2) Choose the next point as the one with an angle closest to 30 degrees.
1273 next_point = min(next_points_in_same_layer, key=lambda v : np.abs(v[1]-v[0]))
1275 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
1276 points = np.append(points, np.array([next_point]), axis=0) # Add the selected point to points
1277 points = sorted(points, key=lambda p: point_relative_magnitude(p)) # Sort the points by their distance from the origin
1279 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
1280 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
1283 if len(points) == N_points_per_hexant:
1284 # Check if the outer-inner boundary ratio is satisfied.
1285 # 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.
1286 # The innermost point is then removed and replaced by a point in the outermost layer.
1287 if point_relative_magnitude(points[-1]) / point_relative_magnitude(points[0]) <= outer_inner_boundary_ratio:
1288 return points
1289 else:
1290 # N_points_per_hexant += 1 # Removed: it increases the number of filaments.
1291 points = np.delete(points, 0, axis=0)
1292 return create_first_hexant_points(N_points_per_hexant, outer_inner_boundary_ratio, possible_next_points, points)
1294 else:
1295 return create_first_hexant_points(N_points_per_hexant, outer_inner_boundary_ratio, possible_next_points, points)
1297 outer_boundary_radius -= filament_radius*1.1
1298 if inner_boundary_radius != 0:
1299 inner_boundary_radius += filament_radius*1.1
1301 outer_inner_boundary_ratio = outer_boundary_radius / inner_boundary_radius if inner_boundary_radius != 0 else 1e10
1302 if N_points % 6 == 1 and inner_boundary_radius == 0:
1303 N_points -= 1
1304 points = [[[0,0,0]]]
1305 else:
1306 points = []
1308 if N_points != 0:
1309 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
1311 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
1312 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
1314 # Rotate the points to create the other hexants of points
1315 for point in groups_of_rotational_symmetry_first_points:
1316 transformed_point = np.matmul(cart_to_hex_transformation, point) # Transform the point from the cartesian coordinate system to the hexagonal coordinate system
1317 point_group = []
1318 rotation_angle = 0
1319 for hexagon_side in range(0, 6):
1320 rotation_matrix = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle)], [np.sin(rotation_angle), np.cos(rotation_angle)]])
1322 rotated_point = np.matmul(rotation_matrix, transformed_point)
1323 point_group.append(list(rotated_point)+[0])
1325 rotation_angle += np.pi/3
1326 points.append(point_group)
1327 return points
1329 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, hole_radius : float = 0.0, hexagonal_holes: bool = False):
1330 """
1331 Creates the full geometry of the strand cross-section.
1333 This method generates the geometry by creating a grid of points that represent the center points of the filaments.
1334 It then creates the filaments at these points. The filaments can be either hexagonal or circular, depending on the `hexagonal_filaments` parameter.
1335 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'.
1337 :param filament_radius: The radius of the filaments.
1338 :type filament_radius: float
1339 :param hexagonal_filaments: If True, the filaments are hexagonal. If False, the filaments are circular.
1340 :type hexagonal_filaments: bool
1341 :param N_filaments: The number of filaments to create.
1342 :type N_filaments: int
1343 :param matrix_inner_radius: The radius of the inner area of the matrix.
1344 :type matrix_inner_radius: float
1345 :param matrix_middle_radius: The radius of the middle area of the matrix. This is where the filaments are located.
1346 :type matrix_middle_radius: float
1347 :param matrix_outer_radius: The radius of the outer area of the matrix.
1348 :type matrix_outer_radius: float
1349 :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.
1350 :type circular_filament_distribution: bool, optional
1351 :param hole_radius: The radius of the filaments holes.
1352 :type hole_radius: float, optional
1353 """
1354 # 1) Create the point-grid representing the center points of the filaments
1355 # 1.1) Get the point positions
1356 filament_centers = self.create_filament_center_points(N_filaments, filament_radius, matrix_middle_radius, matrix_inner_radius, circular_filament_distribution)
1357 # 1.2) Create the center points
1358 center_points = []
1359 for layer in filament_centers:
1360 layer_points = []
1361 for pos in layer:
1362 P = Point.create_or_get(pos)
1363 layer_points.append(P)
1364 center_points.append(layer_points)
1366 # 2) Create the filaments
1367 filaments = []
1368 holes = []
1369 for layer_n in range(len(center_points)):
1370 layer_filaments = []
1371 layer_holes = []
1372 for point_i in range(len(center_points[layer_n])):
1373 if hexagonal_filaments:
1374 filament = Hexagon(center_points[layer_n][point_i], filament_radius, np.pi/6)
1375 else:
1376 filament = Disk(center_points[layer_n][point_i], filament_radius)
1377 if hole_radius:
1378 if hexagonal_holes:
1379 hole = Hexagon(center_points[layer_n][point_i], hole_radius, np.pi / 6)
1380 else:
1381 hole = Disk(center_points[layer_n][point_i], hole_radius)
1382 layer_holes.append(hole)
1383 filament.inner_boundary_curves = [hole.boundary_curves]
1384 layer_filaments.append(filament)
1385 filaments.append(layer_filaments)
1386 if hole_radius:
1387 holes.append(layer_holes)
1388 self.filaments = filaments
1389 if hole_radius:
1390 self.filament_holes = holes
1391 # 3) Create the matrix
1392 # 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'.
1393 # No inner section will be made if the matrix inner_radius is 0.
1394 if matrix_inner_radius != 0:
1395 inner_section = Disk([0,0,0], matrix_inner_radius)
1396 middle_section = Disk([0,0,0], matrix_middle_radius) # Middle section
1397 middle_section.inner_boundary_curves.append(inner_section.boundary_curves)
1398 for layer in self.filaments:
1399 for filament in layer:
1400 middle_section.inner_boundary_curves.append(filament.boundary_curves)
1401 else:
1402 middle_section = Disk([0,0,0], matrix_middle_radius) # Middle section
1403 for layer in self.filaments:
1404 for filament in layer:
1405 middle_section.inner_boundary_curves.append(filament.boundary_curves)
1407 outer_section = Disk([0,0,0], matrix_outer_radius)
1408 outer_section.inner_boundary_curves.append(middle_section.boundary_curves)
1410 self.matrix = [middle_section, outer_section] if matrix_inner_radius == 0 else [inner_section, middle_section, outer_section]
1412 def add_air(self, rad : str, coil_rad : str, type : str):
1413 # 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.
1414 def determine_strand_boundary_single_air_domain(matrix):
1415 """
1416 This function finds the combined outer boundary of the strand geometry, which is the inner boundary of the air region.
1417 The outer boundary of the strand geometry is is not necessarily the outer boundary of the matrix, as the outer matrix partition
1418 may not fully contain the full strand (as with a WIRE-IN-CHANNEL geometry).
1419 """
1420 strand_outer_boundary = set(matrix[0].boundary_curves) # Start with the boundary of the inner matrix partition
1421 for i, matrix_partition in enumerate(matrix[:-1]): # Loop over the matrix partitions
1422 next_matrix_partition = matrix[i+1]
1424 inner_partition_boundary = set(matrix_partition.boundary_curves)
1425 next_partition_boundary = set(next_matrix_partition.boundary_curves)
1427 if inner_partition_boundary & next_partition_boundary: # If the inner and outer partition boundaries share some curves
1428 strand_outer_boundary = strand_outer_boundary ^ next_partition_boundary # Get the combined boundary of the inner and outer partition boundaries
1429 else:
1430 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.
1432 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.
1433 return strand_outer_boundary
1435 if type == 'strand_only': # circle w. natural boundary
1436 self.air.append(Disk([0,0,0], rad))
1437 air_inner_boundaries = determine_strand_boundary_single_air_domain(self.matrix)
1438 self.air[0].inner_boundary_curves.extend([air_inner_boundaries])
1439 elif type == 'coil': # offset semicircle w. natural boundary
1440 self.air.append(Semicircle([0,0,0], coil_rad, rad))
1441 air_inner_boundaries = determine_strand_boundary_single_air_domain(self.matrix)
1442 self.air[0].inner_boundary_curves.extend([air_inner_boundaries])
1443 elif type == 'periodic_square':
1444 outer_matrix_curves = self.matrix[-1].boundary_curves # use matrix boundaries to initialize segmented air sections
1445 self.air = [SquareSection([0,0,0], rad, outer_matrix_curves[0]),
1446 SquareSection([0,0,0], rad, outer_matrix_curves[1]),
1447 SquareSection([0,0,0], rad, outer_matrix_curves[2]),
1448 SquareSection([0,0,0], rad, outer_matrix_curves[3])]
1449 self.air_composition = CompositeSquare(self.air)
1450 self.air_composition.inner_boundary_curves.extend([self.matrix[-1].boundary_curves])
1451 else:
1452 raise ValueError(
1453 f"FiQuS does not support type: {type} with coil radius: {coil_rad}!"
1454 )
1456 def update_tags(self):
1457 """
1458 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.
1459 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.
1460 """
1461 surfaces = sum(self.filaments, []) + sum(self.filament_holes, []) + self.matrix + self.air
1463 Surface.update_tags(surfaces)
1465 def create_gmsh_instance(self):
1466 """
1467 Creates the gmsh instances of the geometry.
1468 """
1469 surfaces = sum(self.filaments, []) + sum(self.filament_holes, []) + self.matrix + self.air
1471 Surface.set_correct_boundary_orientation(surfaces)
1473 for surface in surfaces:
1474 surface.create_gmsh_instance()
1476 def add_physical_groups(self):
1477 """
1478 Creates all physical groups.
1479 """
1480 # Filaments: Add physical boundary and surface
1481 for layer_n in range(len(self.filaments)):
1482 for filament_i in range(len(self.filaments[layer_n])):
1483 self.filaments[layer_n][filament_i].add_physical_boundary(f"Boundary: Filament {filament_i} in layer {layer_n}")
1484 self.filaments[layer_n][filament_i].add_physical_surface(f"Surface: Filament {filament_i} in layer {layer_n}")
1486 self.filaments[layer_n][filament_i].physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.filaments[layer_n][filament_i].boundary_curves[0].points[0].tag])
1488 # Add physical surface for the filament holes
1489 for layer_n in range(len(self.filament_holes)):
1490 for filament_i in range(len(self.filament_holes[layer_n])):
1491 self.filament_holes[layer_n][filament_i].add_physical_boundary(f"Boundary: Filament hole {filament_i} in layer {layer_n}")
1492 self.filament_holes[layer_n][filament_i].add_physical_surface(f"Surface: Filament hole {filament_i} in layer {layer_n}")
1494 self.filament_holes[layer_n][filament_i].physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.filament_holes[layer_n][filament_i].boundary_curves[0].points[0].tag])
1496 # Matrix: Add physical boundary and surface for each partition
1497 for i, matrix_partition in enumerate(self.matrix):
1498 matrix_partition.add_physical_boundary(f"Boundary: Matrix partition {i}")
1499 matrix_partition.add_physical_surface(f"Surface: Matrix partition {i}")
1501 # Air: Add physical boundary and surfaces
1502 for i, air_partition in enumerate(self.air):
1503 air_partition.add_physical_boundary(f"Boundary: Air partition {i}")
1504 air_partition.add_physical_surface(f"Surface: Air partition {i}")
1506 if self.air_composition:
1507 # Cut: Add physical cuts, boundaries and the composite surface
1508 self.air_composition.add_physical_boundaries(f"Boundary: Air")
1509 self.air_composition.add_physical_surface(f"Surface: Air")
1511 self.air_composition.add_physical_cuts(f"Cut: Air")
1513 self.air_composition.physical_inner_boundary_tags.append(gmsh.model.addPhysicalGroup(1, [curve.tag for curve in self.air_composition.inner_boundary_curves[0]], name = f"InnerBoundary: Air"))
1514 self.air_composition.physical_inner_boundary_names.append(f"InnerBoundary: Air")
1515 self.air_composition.strand_bnd_physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.air_composition.inner_boundary_curves[0][0].points[0].tag])
1516 else:
1517 # Add inner boundary
1518 self.air[0].physical_inner_boundary_tags.append(gmsh.model.addPhysicalGroup(1, [curve.tag for curve in self.air[0].inner_boundary_curves[0]], name = f"InnerBoundary: Air"))
1519 self.air[0].physical_inner_boundary_names.append(f"InnerBoundary: Air")
1520 self.air[0].strand_bnd_physicalEdgePointTag = gmsh.model.addPhysicalGroup(0, [self.air[0].inner_boundary_curves[0][0].points[0].tag])
1522 # TEST add a physical group CUT through whole domain
1523 # tag = Line.create_or_get(self.air[0].center_point, self.air[0].boundary_curves[0].P1)
1524 # self.domain_cut = gmsh.model.add_physical_group(1, [tag], name = "Domain cut")
1526 def save(self, save_file):
1527 # This function saves the geometry class to a pickle file.
1528 # The geometry class is used again during meshing.
1529 with open(save_file, "wb") as geom_save_file:
1530 pickle.dump(self, geom_save_file)
1531 logger.info(f"Geometry saved to file {save_file}")
1533 def write_geom_to_yaml(self, file_path):
1534 # This function writes the geometry to a yaml file.
1535 # 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.
1536 # Note: Only the strands are written to the yaml file. The air region is not included.
1537 Conductor = geom.Conductor() # Create a data model for the conductor data
1539 # Add Materials to the 'Solution' section of the data model
1540 # 1) Find all unique materials used in the geometry
1541 # materials = set([surface.material for surface in Surface.surfaces_registry if surface.material is not None])
1542 # # Sort the materials into two groups by their type. All materials with the same type are grouped together.
1543 # 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])}
1544 # # 2) Add all unique materials to the data model, represented by a string with the material type and index in the material group
1545 # for material_type, material_group in material_groups.items():
1546 # for i, material in enumerate(material_group):
1547 # material_name = f"{material_type}_{i}"
1548 # Conductor.Solution.Materials[material_name] = material
1550 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.
1551 curves = list(set(sum([surface.boundary_curves for surface in surfaces], []))) # Extract all the boundary curves from the surfaces
1552 points = list(set(sum([curve.points for curve in curves], []))) # Extract all the points from the curves
1555 # Populate the points dictionary with coordinates of each point
1556 # for p, point in enumerate(points):
1557 # Conductor.Geometry.Points[p] = geom.Point(Coordinates=point.pos.tolist()) #{"Coordinates": point.pos.tolist()}
1558 for p, point in enumerate(Point.points_registry):
1559 if point in points:
1560 Conductor.Geometry.Points[p] = geom.Point(Coordinates=point.pos.tolist())
1562 # Populate the curves dictionary with type of each curve and indices of its points
1563 # for c, curve in enumerate(curves):
1564 # curve_points = [Point.points_registry.index(point) for point in curve.points]
1565 # Conductor.Geometry.Curves[c] = geom.Curve(
1566 # Type=curve.__class__.__name__,
1567 # Points=curve_points
1568 # )
1569 for c, curve in enumerate(Curve.curves_registry):
1570 if curve in curves:
1571 curve_points = [Point.points_registry.index(point) for point in curve.points]
1572 Conductor.Geometry.Curves[c] = geom.Curve(
1573 Type=curve.__class__.__name__,
1574 Points=curve_points
1575 )
1578 # Populate the surfaces dictionary with material, boundary curves and inner boundary curves of each surface
1579 for a, surface in enumerate(Surface.surfaces_registry):
1580 if surface in surfaces:
1581 surface_boundary_curves = [Curve.curves_registry.index(curve) for curve in surface.boundary_curves]
1582 surface_inner_boundary_curves = [[Curve.curves_registry.index(curve) for curve in inner_boundary] for inner_boundary in surface.inner_boundary_curves]
1584 if surface in self.matrix: # Add dummy values for writing the matrix surfaces to the data model
1585 surface.layer = None
1586 surface.layer_index = None
1588 elif surface in sum(self.filaments, []): # Add dummy values for writing the filament surfaces to the data model
1589 for l, layer in enumerate(self.filaments):
1590 if surface in layer:
1591 surface.layer = l
1592 surface.layer_index = layer.index(surface)
1593 break
1595 # Name the material based on its type and index in the material groups
1596 # if surface.material is None:
1597 # material_name = None
1598 # else:
1599 # material_type = surface.material.Type
1600 # material_index = material_groups[material_type].index(surface.material)
1601 # material_name = f"{material_type}_{material_index}"
1602 material_name = surface.material
1604 Conductor.Geometry.Areas[a] = geom.Area(
1605 Material=material_name,
1606 Boundary=surface_boundary_curves,
1607 InnerBoundaries=surface_inner_boundary_curves,
1608 Layer=surface.layer,
1609 LayerIndex=surface.layer_index
1610 )
1613 # Write the data model to a yaml file
1614 UFF.write_data_to_yaml(file_path, Conductor.model_dump())
1616 @classmethod
1617 def read_geom_from_yaml(cls, file_path):
1618 """
1619 This function loads a geometry from a yaml file and returns a TwistedStrand object.
1620 The yaml file contains all points, curves and surfaces which define the geometry.
1621 - Points are defined by their position vector and can be referenced by an integer ID.
1622 : Position [x, y, z]
1623 - 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.
1624 : Type ('Line', 'CircleArc', etc.)
1625 : Points ([1,2,3]), list of point-ID defining the curve.
1626 - Surfaces are defined by material, outer boundary, inner boundaries and layer and layer-index if the surface is a strand.
1627 : Material ('Cu', 'NbTi', 'Nb3Sn', etc.)
1628 : Outer boundary ([2,3,4,5...]), list of curve-ID defining the outer boundary closed loop.
1629 : Inner boundaries ([[1,2,3], [4,5,6], ... ]), list of list of curve-IDs defining closed loops.
1630 : Layer (0, 1, 2, ...), the layer of a filament. None if the surface is part of the matrix.
1631 : LayerIndex (0, 1, 2, ...), the index of the filament in the layer. None if the surface is part of the matrix.
1633 :param file_path: The full path to the yaml file.
1634 :type file_path: str
1635 :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.
1637 """
1638 Conductor = UFF.read_data_from_yaml(file_path, geom.Conductor)
1640 # 1) Create the points
1641 for point in Conductor.Geometry.Points.values():
1642 Point.create_or_get(point.Coordinates)
1644 # 2) Create the curves
1645 for curve in Conductor.Geometry.Curves.values():
1646 curve_type = curve.Type
1647 points = [Point.points_registry[p] for p in curve.Points]
1650 c = globals()[curve_type].create_or_get(*points) # Create the curve of the specified type
1652 # TODO: To be added.
1653 # c.contact = curve.Contact
1654 # c.thickness = curve.Thickness
1655 # c.material = curve.Material
1658 # 3) Create the surfaces
1659 # TODO: area.boundary_material and boundary_thickness are not yet used.
1660 strand = cls()
1661 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
1662 strand.filaments = [[None for i in range(6)] for j in range(layers)] # Initialize the filaments list
1663 strand.filament_holes = [[None for i in range(6)] for j in range(layers)] # Initialize the filament holes list
1665 for area_index, area_dm in Conductor.Geometry.Areas.items():
1666 boundary_curves = [Curve.curves_registry[c] for c in area_dm.Boundary]
1667 inner_boundary_curves = [[Curve.curves_registry[c] for c in inner_boundary] for inner_boundary in area_dm.InnerBoundaries]
1668 surface = Surface(boundary_curves, inner_boundary_curves)
1669 if area_dm.Material: # If the material is provided
1670 surface.material = area_dm.Material
1672 if area_dm.Layer is None:
1673 # It is either a matrix partition or it is a hole in a filament.
1674 # 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.
1675 is_hole = False
1676 for other_area in Conductor.Geometry.Areas.values(): # Loop over all areas to check if the area is a hole in a filament
1677 if other_area.Layer is not None: # If the other area is a filament
1678 if area_dm.Boundary in other_area.InnerBoundaries: # If the area is a hole in the other filament
1679 # 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
1680 layer = other_area.Layer
1681 layer_index = other_area.LayerIndex
1682 strand.filament_holes[layer][layer_index] = surface
1683 is_hole = True
1684 break
1686 if not is_hole: # If it is not a hole, it is a matrix partition
1687 strand.matrix.append(surface)
1689 else:
1690 strand.filaments[area_dm.Layer][area_dm.LayerIndex] = surface
1692 # Remove None values from the filaments list
1693 strand.filaments = [[filament for filament in layer if filament is not None] for layer in strand.filaments]
1694 strand.filament_holes = [[hole for hole in layer if hole is not None] for layer in strand.filament_holes]
1697 # Sort the matrix partitions from inner to outer based on the outermost point of the boundary curves of the partitions
1698 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]))
1700 return strand
1702class Geometry:
1703 """
1704 Class to generate the ConductorAC Strand geometry.
1706 This class is responsible for generating the geometry of the twisted strand.
1707 It can either load a geometry from a YAML file or create the model from scratch.
1708 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.
1710 :ivar fdm: The fiqus inputs data model.
1711 :vartype fdm: object
1712 :ivar cacdm: The magnet section of the fiqus inputs data model.
1713 :vartype cacdm: object
1714 :ivar inputs_folder: The full path to the folder with input files, i.e., conductor and STEP files.
1715 :vartype inputs_folder: str
1716 :ivar geom_folder: The full path to the current working directory.
1717 :vartype geom_folder: str
1718 :ivar magnet_name: The name of the magnet.
1719 :vartype magnet_name: str
1720 :ivar geom_file: The full path to the .brep file where the geometry will be saved.
1721 :vartype geom_file: str
1722 :ivar verbose: If True, more information is printed in the Python console.
1723 :vartype verbose: bool
1724 :ivar gu: An instance of the GmshUtils class.
1725 :vartype gu: object
1726 """
1728 def __init__(self, fdm, inputs_folder_path, verbose=True):
1729 """
1730 Initializes the Geometry class.
1732 :param fdm: The fiqus data model.
1733 :type fdm: object
1734 :param inputs_folder_path: The full path to the folder with input files, i.e., conductor and STEP files.
1735 :type inputs_folder_path: str
1736 :param verbose: If True, more information is printed in the Python console. Defaults to True.
1737 :type verbose: bool, optional
1738 """
1739 self.fdm = fdm
1740 self.cacdm = fdm.magnet
1741 self.inputs_folder = inputs_folder_path
1742 self.geom_folder = os.path.join(os.getcwd())
1743 self.magnet_name = fdm.general.magnet_name
1744 self.geom_file = os.path.join(self.geom_folder, f'{self.magnet_name}.brep')
1745 self.verbose = verbose
1746 self.gu = GmshUtils(self.geom_folder, self.verbose)
1747 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh)
1749 # To see the surfaces in a better way in GUI:
1750 gmsh.option.setNumber("Geometry.SurfaceType", 2)
1752 def generate_strand_geometry(self, gui=False):
1753 """
1754 Generates the geometry of a strand based on the settings specified in the input data model.
1755 The geometry can either be loaded from a YAML file or created from scratch. After creation, the geometry
1756 can be saved to a YAML file. The method also creates gmsh instances, adds physical groups to the geometry,
1757 writes the geometry to a .brep file, and saves the geometry-class to a pickle file. If `gui` is True, it
1758 launches an interactive GUI.
1760 :param gui: If True, launches an interactive GUI after generating the geometry. Default is False.
1761 :type gui: bool, optional
1763 :return: None
1764 """
1765 logger.info("Generating geometry")
1766 #0) Clear the registries. Used when generating reference files for tests.
1767 if Point.points_registry: # If the points registry is not empty, clear it.
1768 Point.points_registry.clear()
1769 if Curve.curves_registry: # If the curves registry is not empty, clear it.
1770 Curve.curves_registry.clear()
1771 if Surface.surfaces_registry: # If the surfaces registry is not empty, clear it.
1772 Surface.surfaces_registry.clear()
1773 if Surface.curve_loop_registry: # If the curve loop registry is not empty, clear it.
1774 Surface.curve_loop_registry.clear()
1776 # 1) Either load the geometry from a yaml file or create the model from scratch
1777 if self.cacdm.geometry.io_settings.load.load_from_yaml:
1778 CAC = TwistedStrand.read_geom_from_yaml(os.path.join(self.inputs_folder, self.cacdm.geometry.io_settings.load.filename))
1779 else:
1780 CAC = TwistedStrand()
1781 strand = self.fdm.conductors[self.cacdm.solve.conductor_name].strand
1782 if strand.filament_hole_diameter:
1783 filament_hole_radius = strand.filament_hole_diameter/2
1784 if strand.filament_hole_diameter >= strand.filament_diameter:
1785 raise ValueError(
1786 f"Invalid strand geometry: filament_hole_diameter ({strand.filament_hole_diameter}) "
1787 f"must be smaller than filament_diameter ({strand.filament_diameter})."
1788 )
1789 else:
1790 filament_hole_radius = 0.0
1791 CAC.create_geometry(
1792 strand.filament_diameter/2,
1793 self.cacdm.geometry.hexagonal_filaments,
1794 strand.number_of_filaments,
1795 strand.diameter_core/2,
1796 strand.diameter_filamentary/2,
1797 strand.diameter/2,
1798 self.cacdm.geometry.filament_circular_distribution,
1799 filament_hole_radius,
1800 self.cacdm.geometry.hexagonal_holes
1801 )
1803 CAC.add_air(self.cacdm.geometry.air_radius, self.cacdm.geometry.coil_radius, self.cacdm.geometry.type)
1804 CAC.create_gmsh_instance()
1805 # 2) Save the geometry to a yaml file if specified
1806 if self.cacdm.geometry.io_settings.save.save_to_yaml:
1807 filename = self.cacdm.geometry.io_settings.save.filename
1808 CAC.write_geom_to_yaml(os.path.join(self.geom_folder, filename))
1810 if self.cacdm.geometry.rotate_angle:
1811 dimTags = gmsh.model.occ.getEntities(dim=-1)
1812 gmsh.model.occ.rotate(dimTags, 0.0, 0.0, 0.0, 0, 0, 1, self.cacdm.geometry.rotate_angle * np.pi / 180)
1814 gmsh.model.occ.synchronize()
1816 # Add physical groups to the geometry
1817 CAC.add_physical_groups()
1819 logger.info("Writing geometry")
1820 gmsh.write(self.geom_file) # Write the geometry to a .brep file
1821 CAC.save(os.path.join(self.geom_folder, f'{self.magnet_name}.pkl')) # Save the geometry-class to a pickle file
1823 if gui:
1824 self.gu.launch_interactive_GUI()
1825 else:
1826 if gmsh.isInitialized():
1827 gmsh.clear()
1828 gmsh.finalize()
1830 def load_conductor_geometry(self, gui=False):
1831 """
1832 Loads geometry from .brep file.
1833 """
1835 logger.info("Loading geometry")
1837 gmsh.clear()
1838 gmsh.model.occ.importShapes(self.geom_file, format="brep")
1839 gmsh.model.occ.synchronize()
1841 if gui:
1842 self.gu.launch_interactive_GUI()