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