From 9db19401d7d3eab5ac668d9c70e620e8fc371cfe Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 31 Aug 2025 15:40:12 +0200 Subject: [PATCH 01/36] feat(sector): add Sector and SectorInsulator types with geometric properties --- src/DataModel.jl | 8 +- src/datamodel/preview.jl | 29 +++++ src/datamodel/sector.jl | 214 +++++++++++++++++++++++++++++++ src/datamodel/sectorinsulator.jl | 91 +++++++++++++ 4 files changed, 340 insertions(+), 2 deletions(-) create mode 100644 src/datamodel/sector.jl create mode 100644 src/datamodel/sectorinsulator.jl diff --git a/src/DataModel.jl b/src/DataModel.jl index d2d96cbb..cf7e514a 100644 --- a/src/DataModel.jl +++ b/src/DataModel.jl @@ -23,8 +23,8 @@ module DataModel # Export public API export Thickness, Diameter # Type definitions -export WireArray, Strip, Tubular # Conductor types -export Semicon, Insulator # Insulator types +export WireArray, Strip, Tubular, SectorParams, Sector # Conductor types +export Semicon, Insulator, SectorInsulator # Insulator types export ConductorGroup, InsulatorGroup # Group types export CableComponent, CableDesign # Cable design types export CablePosition, LineCableSystem # System types @@ -44,6 +44,8 @@ using DataFrames using Colors using Plots using DisplayAs: DisplayAs +using GeometryBasics +using PolygonOps # Abstract types & interfaces include("datamodel/types.jl") @@ -62,11 +64,13 @@ include("datamodel/wirearray.jl") include("datamodel/strip.jl") include("datamodel/tubular.jl") include("datamodel/conductorgroup.jl") +include("datamodel/sector.jl") # Insulators include("datamodel/insulator.jl") include("datamodel/semicon.jl") include("datamodel/insulatorgroup.jl") +include("datamodel/insulatorgroup.jl") # Groups diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index 04af2906..e24a5a82 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -275,6 +275,35 @@ function preview( color=color, label=display_legend ? label : "", ) + elseif layer isa Sector + vertices = layer.vertices + material_props = layer.material_props + color = _get_material_color(material_props) + + plot!( + plt, + Shape(vertices), + linecolor=:black, + color=color, + label=display_legend ? label : "", + ) + elseif layer isa SectorInsulator + outer_vertices = layer.outer_vertices + # The inner boundary is the conductor's vertices. It must be reversed for the hole to be drawn correctly. + inner_vertices = reverse(layer.inner_sector.vertices) + material_props = layer.material_props + color = _get_material_color(material_props) + + # Create a shape with a hole by passing a list of polygons + plot!( + plt, + #Shape([outer_vertices, inner_vertices]), + Shape(outer_vertices), + linecolor=:black, + color=color, + label=display_legend ? label : "", + ) + end end diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl new file mode 100644 index 00000000..c5066566 --- /dev/null +++ b/src/datamodel/sector.jl @@ -0,0 +1,214 @@ +""" +$(TYPEDEF) + +Holds the geometric parameters that define the shape of a sector conductor. + +$(TYPEDFIELDS) +""" +struct SectorParams{T<:REALSCALAR} + "Number of sectors in the full cable (e.g., 3 or 4)." + n_sectors::Int + "Back radius of the sector (outermost curve) \\[m\\]." + r_back::T + "Depth of the sector from the back to the base \\[m\\]." + d_sector::T + "Corner radius for rounding sharp edges \\[m\\]." + r_corner::T + "Angular width of the conductor's flat base/sides [degrees]." + theta_cond_deg::T +end + +""" +$(TYPEDEF) + +Represents a single sector-shaped conductor with defined geometric and material properties. + +$(TYPEDFIELDS) +""" +struct Sector{T<:REALSCALAR} <: AbstractConductorPart{T} + "Inner radius (not applicable, typically 0 for the central point) \\[m\\]." + radius_in::T + "Outer radius (equivalent back radius) \\[m\\]." + radius_ext::T + "Geometric parameters defining the sector's shape." + params::SectorParams{T} + "Rotation angle of this specific sector around the cable's center [degrees]." + rotation_angle_deg::T + "Material properties of the conductor." + material_props::Material{T} + "Operating temperature of the conductor \\[°C\\]." + temperature::T + "Cross-sectional area of the sector \\[m²\\]." + cross_section::T + "Electrical resistance of the sector \\[Ω/m\\]." + resistance::T + "Geometric mean radius (GMR) of the sector (approximated) \\[m\\]." + gmr::T + "Calculated vertices defining the polygon shape." + vertices::Vector{Point{2,T}} +end + + +function Sector( + params::SectorParams{T}, + rotation_angle_deg::T, + material_props::Material{T}; + temperature::T=T₀, +) where {T<:REALSCALAR} + + # 1. Calculate the geometry and vertices for a base (unrotated) sector + base_vertices = _calculate_sector_polygon_points(params) + + # 2. Rotate the vertices to the specified angle + rotation_angle_rad = deg2rad(rotation_angle_deg) + rotated_vertices = [_rotate_point(p, rotation_angle_rad) for p in base_vertices] + + # 3. Calculate cross-sectional area using the Shoelace formula + cross_section = PolygonOps.area(rotated_vertices) + + # 4. Calculate DC resistance + rho_eff = calc_temperature_correction(material_props.alpha, temperature, material_props.T0) * material_props.rho + resistance = rho_eff / cross_section + + # 5. Approximate GMR based on a circle of equivalent area + r_equiv = sqrt(cross_section / π) + gmr = r_equiv * exp(-0.25 * material_props.mu_r) # GMR of an equivalent solid round conductor + + return Sector{T}( + zero(T), + params.r_back, + params, + rotation_angle_deg, + material_props, + temperature, + cross_section, + resistance, + gmr, + rotated_vertices, + ) +end + + + +# --- Geometric Helper Functions (internal) --- + +function _calculate_sector_geometry(p::SectorParams) + phi_deg = 90.0 - p.theta_cond_deg / 2.0 + phi_rad = deg2rad(phi_deg) + + if abs(cos(phi_rad)) < 1e-9 + error("theta_cond_deg is too close to 180, leading to division by zero. Check parameters.") + end + + d_base_corner = p.r_corner * (1.0 / cos(phi_rad) - 1.0) + d_offset = p.r_back - p.d_sector - d_base_corner + + x_base_corner = p.r_corner * sin(phi_rad) + y_base_corner = x_base_corner * tan(phi_rad) + d_offset + node_A = Point2f(x_base_corner, y_base_corner) + + k = p.r_corner / cos(phi_rad) + d_offset + qa = 1.0 + tan(phi_rad)^2 + qb = 2.0 * k * tan(phi_rad) + qc = k^2 - (p.r_back - p.r_corner)^2 + + discriminant = qb^2 - 4.0 * qa * qc + if discriminant < 0 + error("Cannot calculate side corner center: negative discriminant ($discriminant). Check parameters.") + end + + x_side_center = (-qb + sqrt(discriminant)) / (2.0 * qa) + y_side_center = x_side_center * tan(phi_rad) + k + + x_side_lower = x_side_center + p.r_corner * sin(phi_rad) + y_side_lower = y_side_center - p.r_corner * cos(phi_rad) + node_B = Point2f(x_side_lower, y_side_lower) + + dist_origin_side_center = sqrt(x_side_center^2 + y_side_center^2) + if dist_origin_side_center < 1e-9 + error("Side corner center is at the origin, cannot determine upper point direction.") + end + x_side_upper = x_side_center * p.r_back / dist_origin_side_center + y_side_upper = y_side_center * p.r_back / dist_origin_side_center + node_C = Point2f(x_side_upper, y_side_upper) + + nodes = ( + A=node_A, B=node_B, C=node_C, D=Point2f(-node_C[1], node_C[2]), + E=Point2f(-node_B[1], node_B[2]), F=Point2f(-node_A[1], node_A[2]) + ) + centers = ( + Back=Point2f(0, 0), Base=Point2f(0, d_offset+ p.r_corner / cos(phi_rad)), + RightSide=Point2f(x_side_center, y_side_center), LeftSide=Point2f(-x_side_center, y_side_center) + ) + return (Nodes=nodes, Centers=centers, Params=p) +end + +function _generate_arc_points(center, radius, start_angle, end_angle, num_points) + while end_angle < start_angle + @debug "(Sector) end_angle < start_angle: $(rad2deg(end_angle)) < $(rad2deg(start_angle))" + end_angle += 2pi + @debug "(Sector) Adjusted end_angle to be greater than start_angle: $(rad2deg(end_angle)) > $(rad2deg(start_angle))" + end + while end_angle - start_angle > pi + @debug "(Sector) overshoot! end_angle - start_angle > π: $(rad2deg(end_angle - start_angle)) > 180" + end_angle -= 2pi + @debug "(Sector) Adjusted end_angle to be less than start_angle: $(rad2deg(end_angle)) < $(rad2deg(start_angle))" + end + angle_range = range(start_angle, stop=end_angle, length=num_points) + return [Point2f(center[1] + radius * cos(a), center[2] + radius * sin(a)) for a in angle_range] +end + +function _calculate_sector_polygon_points(params; num_arc_points=20) # increase `num_arc_points` for higher accuracy + geom = _calculate_sector_geometry(params) + nodes, centers = geom.Nodes, geom.Centers + + poly_points = Point2f[] + get_angle(p1, p2) = atan((p1[2] - p2[2]), (p1[1] - p2[1])) + # Start at F, go to A (Base) + push!(poly_points, nodes.F) + if params.r_corner > 1e-9 + start_angle = get_angle(nodes.F, centers.Base) + end_angle = get_angle(nodes.A, centers.Base) + append!(poly_points, + _generate_arc_points(centers.Base, params.r_corner, start_angle, end_angle, num_arc_points)[2:end]) + else + push!(poly_points, Point2f(0, params.r_back - params.d_sector), nodes.A) + end + + # Line A to B + push!(poly_points, nodes.B) + + # Arc B to C (Right Side) + if params.r_corner > 1e-9 + start_angle = get_angle(nodes.B, centers.RightSide) + @debug "Arc from B to C: start_angle=$(rad2deg(start_angle))" + end_angle = get_angle(nodes.C, centers.RightSide) + @debug "Arc from B to C: end_angle=$(rad2deg(end_angle))" + append!(poly_points, _generate_arc_points(centers.RightSide, params.r_corner, start_angle, end_angle, num_arc_points)[2:end]) + else + push!(poly_points, nodes.C) + end + + # Arc C to D (Back) + start_angle = get_angle(nodes.C, centers.Back) + end_angle = get_angle(nodes.D, centers.Back) + append!(poly_points, _generate_arc_points(centers.Back, params.r_back, start_angle, end_angle, num_arc_points)[2:end]) + + # Arc D to E (Left Side) + if params.r_corner > 1e-9 + start_angle = get_angle(nodes.D, centers.LeftSide) + @debug "Arc from D to E: start_angle=$(rad2deg(start_angle))" + end_angle = get_angle(nodes.E, centers.LeftSide) + @debug "Arc from D to E: end_angle=$(rad2deg(end_angle))" + append!(poly_points, _generate_arc_points(centers.LeftSide, params.r_corner, start_angle, end_angle, num_arc_points)[2:end]) + else + push!(poly_points, nodes.E) + end + + return poly_points +end + +function _rotate_point(p::Point2f, angle_rad::Real) + cos_a, sin_a = cos(angle_rad), sin(angle_rad) + return Point2f(p[1] * cos_a - p[2] * sin_a, p[1] * sin_a + p[2] * cos_a) +end \ No newline at end of file diff --git a/src/datamodel/sectorinsulator.jl b/src/datamodel/sectorinsulator.jl new file mode 100644 index 00000000..083838b1 --- /dev/null +++ b/src/datamodel/sectorinsulator.jl @@ -0,0 +1,91 @@ +""" +$(TYPEDEF) + +Represents an insulating layer surrounding a sector-shaped conductor. + +$(TYPEDFIELDS) +""" +struct SectorInsulator{T<:REALSCALAR} <: AbstractInsulatorPart{T} + "Inner radius (not applicable, defined by inner sector) \\[m\\]." + radius_in::T + "Outer radius (equivalent back radius of outer boundary) \\[m\\]." + radius_ext::T + "The inner sector conductor that this insulator surrounds." + inner_sector::Sector{T} + "The thickness of the insulating layer \\[m\\]." + thickness::T + "Material properties of the insulator." + material_props::Material{T} + "Operating temperature of the insulator \\[°C\\]." + temperature::T + "Cross-sectional area of the insulating layer \\[m²\\]." + cross_section::T + "Shunt capacitance (approximated) \\[F/m\\]." + shunt_capacitance::T + "Shunt conductance (approximated) \\[S·m\\]." + shunt_conductance::T + "Calculated vertices of the outer boundary polygon." + outer_vertices::Vector{Point{2, T}} +end + + +function SectorInsulator( + inner_sector::Sector{T}, + thickness::T, + material_props::Material{T}; + temperature::T=T₀, +) where {T<:REALSCALAR} + + # 1. Calculate the outer vertices by offsetting the inner sector's geometry + outer_vertices = _calculate_offset_polygon(inner_sector.params, thickness, inner_sector.rotation_angle_deg) + + # 2. Calculate areas + inner_area = inner_sector.cross_section + outer_area = PolygonOps.area(outer_vertices) + cross_section = outer_area - inner_area + + # 3. Approximate capacitance and conductance using equivalent coaxial circles + # This is more consistent with the package's approach than a parallel plate model. + r_eq_in = sqrt(inner_area / π) + r_eq_ext = sqrt(outer_area / π) + shunt_capacitance = calc_shunt_capacitance(r_eq_in, r_eq_ext, material_props.eps_r) + shunt_conductance = calc_shunt_conductance(r_eq_in, r_eq_ext, material_props.rho) + + # 4. Determine outer radius from the new params + outer_r_back = inner_sector.params.r_back + thickness + + return SectorInsulator{T}( + inner_sector.radius_ext, + outer_r_back, + inner_sector, + thickness, + material_props, + temperature, + cross_section, + shunt_capacitance, + shunt_conductance, + outer_vertices + ) +end + +# --- Geometric Helper Functions (internal) --- +#NB: This return another full sector not just the insulation - +#---- Parameter calculations are handled by the equivalent area difference (as in the constructor) +#---- Plots.jl handling can be done +#---- FEM handling not done yet + +function _calculate_offset_polygon(inner_params::SectorParams, thickness::Real, rotation_angle_deg::Real) + # Create a new set of parameters for the outer boundary by adding the thickness + outer_params = SectorParams( + inner_params.n_sectors, + inner_params.r_back + thickness, + inner_params.d_sector + thickness, + inner_params.r_corner + thickness, # Offset corner radius as well + inner_params.theta_cond_deg + ) + + base_vertices = _calculate_sector_polygon_points(outer_params) + rotation_angle_rad = deg2rad(rotation_angle_deg) + rotated_vertices = [_rotate_point(p, rotation_angle_rad) for p in base_vertices] + return rotated_vertices +end From 231f585ea1d386831cf890c4feac18d0f444c71b Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 31 Aug 2025 16:17:54 +0200 Subject: [PATCH 02/36] fix(sector): fix typo insulator group inclusion to SectorInsulator --- src/DataModel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DataModel.jl b/src/DataModel.jl index cf7e514a..2ec8517f 100644 --- a/src/DataModel.jl +++ b/src/DataModel.jl @@ -70,7 +70,7 @@ include("datamodel/sector.jl") include("datamodel/insulator.jl") include("datamodel/semicon.jl") include("datamodel/insulatorgroup.jl") -include("datamodel/insulatorgroup.jl") +include("datamodel/SectorInsulator.jl") # Groups From b767c690e279d333bab31fde5efeaddc5ad1c798 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 31 Aug 2025 16:19:45 +0200 Subject: [PATCH 03/36] feat(sector): add GeometryBasics and PolygonOps dependencies with version updates --- Project.toml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Project.toml b/Project.toml index 2172cc6d..da1bc8eb 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" EzXML = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615" ForceImport = "9dda63f9-cce7-5873-89fa-eccbb2fffcde" Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196" +GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" GetDP = "dbfd83e6-7aba-450b-9c2b-93ccd973023a" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" @@ -29,6 +30,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" PkgTemplates = "14b8a8f1-9102-5b29-a752-f990bacb7fe1" PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -52,6 +54,7 @@ DocStringExtensions = "0.9.3" EzXML = "1.2.0" ForceImport = "0.0.3" Functors = "0.5.2" +GeometryBasics = "0.5.10" Gmsh = "0.3.1" JSON3 = "1.14.2" LinearAlgebra = "1.11.0" @@ -63,6 +66,7 @@ Pkg = "1.11.0" PkgTemplates = "0.7.56" PlotlyJS = "0.18.15" Plots = "1.40.9" +PolygonOps = "0.1.2" Printf = "1.11.0" QuadGK = "2.11.2" Reexport = "1.2.2" From 237e51768f41c459ca62171ddeb83e53845ea1d6 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 01:30:22 +0200 Subject: [PATCH 04/36] feat(sector): add LinearAlgebra import --- src/DataModel.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/DataModel.jl b/src/DataModel.jl index 2ec8517f..4b9c805a 100644 --- a/src/DataModel.jl +++ b/src/DataModel.jl @@ -46,6 +46,7 @@ using Plots using DisplayAs: DisplayAs using GeometryBasics using PolygonOps +using LinearAlgebra # Abstract types & interfaces include("datamodel/types.jl") From 251e0d089279719873dd27c4318b6b9144b26f43 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 01:31:00 +0200 Subject: [PATCH 05/36] feat(sectorinsulator): update _calculate_offset_polygon to use vertices and improve offset calculation --- src/datamodel/sectorinsulator.jl | 61 ++++++++++++++++++++++---------- 1 file changed, 42 insertions(+), 19 deletions(-) diff --git a/src/datamodel/sectorinsulator.jl b/src/datamodel/sectorinsulator.jl index 083838b1..bfe324f9 100644 --- a/src/datamodel/sectorinsulator.jl +++ b/src/datamodel/sectorinsulator.jl @@ -37,7 +37,7 @@ function SectorInsulator( ) where {T<:REALSCALAR} # 1. Calculate the outer vertices by offsetting the inner sector's geometry - outer_vertices = _calculate_offset_polygon(inner_sector.params, thickness, inner_sector.rotation_angle_deg) + outer_vertices = _calculate_offset_polygon(inner_sector.vertices, thickness) # 2. Calculate areas inner_area = inner_sector.cross_section @@ -69,23 +69,46 @@ function SectorInsulator( end # --- Geometric Helper Functions (internal) --- -#NB: This return another full sector not just the insulation - -#---- Parameter calculations are handled by the equivalent area difference (as in the constructor) -#---- Plots.jl handling can be done -#---- FEM handling not done yet - -function _calculate_offset_polygon(inner_params::SectorParams, thickness::Real, rotation_angle_deg::Real) - # Create a new set of parameters for the outer boundary by adding the thickness - outer_params = SectorParams( - inner_params.n_sectors, - inner_params.r_back + thickness, - inner_params.d_sector + thickness, - inner_params.r_corner + thickness, # Offset corner radius as well - inner_params.theta_cond_deg - ) +# REVISED: This function now takes vertices and thickness to compute a geometric offset. +function _calculate_offset_polygon(vertices::Vector{Point{2, T}}, thickness::T) where {T<:REALSCALAR} + num_vertices = length(vertices) + if num_vertices < 3 + error("Polygon must have at least 3 vertices.") + end + + new_vertices = similar(vertices) + + for i in 1:num_vertices + p_prev = vertices[i == 1 ? num_vertices : i - 1] + p_curr = vertices[i] + p_next = vertices[i == num_vertices ? 1 : i + 1] + + v1 = p_curr - p_prev + v2 = p_next - p_curr + + # Normalize the vectors + v1_norm = v1 / norm(v1) + v2_norm = v2 / norm(v2) + + # Normal vectors (rotated 90 degrees clockwise for outward direction) + n1 = Point(v1_norm[2], -v1_norm[1]) + n2 = Point(v2_norm[2], -v2_norm[1]) + + # Bisector of the normals + bisector = (n1 + n2) / norm(n1 + n2) + + # Angle between the two vectors to calculate the correct offset distance + angle = acos(clamp(v1_norm ⋅ v2_norm, -1.0, 1.0)) + + # Miter length + offset_distance = thickness / sin((π - angle) / 2) + + if isinf(offset_distance) # The vectors are parallel + new_vertices[i] = p_curr + n1 * thickness + else + new_vertices[i] = p_curr + bisector * offset_distance + end + end - base_vertices = _calculate_sector_polygon_points(outer_params) - rotation_angle_rad = deg2rad(rotation_angle_deg) - rotated_vertices = [_rotate_point(p, rotation_angle_rad) for p in base_vertices] - return rotated_vertices + return new_vertices end From 6c5b81eb2c1f85f379e9b7c4a10198f9f4d782cb Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 01:33:16 +0200 Subject: [PATCH 06/36] feat(setor): update preview of insulator to be transparent (to be handled later) --- src/datamodel/preview.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index e24a5a82..cccedd77 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -300,7 +300,7 @@ function preview( #Shape([outer_vertices, inner_vertices]), Shape(outer_vertices), linecolor=:black, - color=color, + color=:transparent, label=display_legend ? label : "", ) From 4a0cc2cbd6e33171affd9afaf1e52db072d4c12b Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 11:14:10 +0200 Subject: [PATCH 07/36] feat(sector): update preview function to handle Sector layers with TODOs for nominal length and offset adjustments --- src/datamodel/preview.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index cccedd77..9d83185d 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -275,7 +275,9 @@ function preview( color=color, label=display_legend ? label : "", ) - elseif layer isa Sector + elseif layer isa Sector + #TODO: support the nominal length (using lay ratio) + #TODO: Fix the offset to be taken into account vertices = layer.vertices material_props = layer.material_props color = _get_material_color(material_props) From 05d99b72568a2274c50a60435d951ebdcd445d93 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 11:57:56 +0200 Subject: [PATCH 08/36] feat(sector): adjust vertex calculations for Sector and SectorInsulator layers to take the offset into account --- src/datamodel/preview.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index 9d83185d..081b56a5 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -277,8 +277,8 @@ function preview( ) elseif layer isa Sector #TODO: support the nominal length (using lay ratio) - #TODO: Fix the offset to be taken into account vertices = layer.vertices + vertices = [(v[1] + x0, v[2] + y0) for v in vertices] material_props = layer.material_props color = _get_material_color(material_props) @@ -290,9 +290,10 @@ function preview( label=display_legend ? label : "", ) elseif layer isa SectorInsulator - outer_vertices = layer.outer_vertices - # The inner boundary is the conductor's vertices. It must be reversed for the hole to be drawn correctly. - inner_vertices = reverse(layer.inner_sector.vertices) + outer_vertices = [(v[1] + x0, v[2] + y0) for v in layer.outer_vertices] + # (Not used for now) The inner boundary is the conductor's vertices. It must be reversed for the hole to be drawn correctly. + inner_vertices = [(v[1] + x0, v[2] + y0) for v in layer.inner_sector.vertices] + inner_vertices = reverse(inner_vertices) material_props = layer.material_props color = _get_material_color(material_props) From 08ce9691de8f538f2ea626a996fe8cdd2bba7565 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 13:57:13 +0200 Subject: [PATCH 09/36] feat(sector-FEM): add functions to draw polygons (from points) and polygons with holes in Gmsh --- src/engine/fem/drawing.jl | 58 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index 669b2f5c..96018db1 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -640,6 +640,64 @@ function draw_transition_region( return tags, all_mesh_points, markers end +function draw_polygon(vertices::Vector{<:Point}) + # Create points + points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in vertices] + + # Create lines + lines = [gmsh.model.occ.add_line(points[i], points[i % length(points) + 1]) for i in 1:length(points)] + + # Create curve loop and surface + curve_loop = gmsh.model.occ.add_curve_loop(lines) + surface_tag = gmsh.model.occ.add_plane_surface([curve_loop]) + + # Synchronize to make the new entity available for calculations + gmsh.model.occ.synchronize() + + # Calculate a center of mass for placing a marker + com = gmsh.model.get_center_of_mass(2, surface_tag) + marker = [com[1], com[2], com[3]] + + return surface_tag, marker +end + +""" +$(TYPEDSIGNATURES) + +Draw a polygon with a hole. + +# Arguments +- `outer_vertices`: A vector of (x,y) coordinates for the outer boundary. +- `inner_vertices`: A vector of (x,y) coordinates for the inner boundary (the hole). + +# Returns +- A tuple containing the Gmsh surface tag and a marker point `[x, y, z]`. +""" +function draw_polygon_with_hole(outer_vertices::Vector{<:Point}, inner_vertices::Vector{<:Point}) + # Create outer boundary + outer_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in outer_vertices] + outer_lines = [gmsh.model.occ.add_line(outer_points[i], outer_points[i % length(outer_points) + 1]) for i in 1:length(outer_points)] + outer_loop = gmsh.model.occ.add_curve_loop(outer_lines) + + # Create inner boundary (hole) + inner_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in inner_vertices] + inner_lines = [gmsh.model.occ.add_line(inner_points[i], inner_points[i % length(inner_points) + 1]) for i in 1:length(inner_points)] + inner_loop = gmsh.model.occ.add_curve_loop(inner_lines) + + # Create surface with hole + surface_tag = gmsh.model.occ.add_plane_surface([outer_loop, inner_loop]) + + # Synchronize to make the new entity available for calculations + gmsh.model.occ.synchronize() + + # A marker point must be inside the insulator, but outside the conductor. + # A point halfway between the inner and outer boundaries along one of the vertices should work. + marker_point = (outer_vertices[1] + inner_vertices[1]) / 2 + marker = [marker_point[1], marker_point[2], 0.0] + + return surface_tag, marker +end + function get_system_centroid(cable_system::LineCableSystem, cable_idx::Vector{<:Integer}) # Check if cable_idx is empty if isempty(cable_idx) From b75a9a40aa6994aa574087d5a35882e0be6ae905 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 13:58:31 +0200 Subject: [PATCH 10/36] feat(sector-FEM): add specialized `_make_cablepart!` methods for creating Sector and SectorInsulator geometries --- src/engine/fem/cable.jl | 92 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) diff --git a/src/engine/fem/cable.jl b/src/engine/fem/cable.jl index 7e3dbbfd..1ef9155f 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -390,4 +390,96 @@ function _make_cablepart!(workspace::FEMWorkspace, part::WireArray, # Add physical groups to the workspace register_physical_group!(workspace, physical_group_tag_air_gap, air_material) end +end + +""" +$(TYPEDSIGNATURES) + +Specialized method to create the geometry for a `Sector` conductor. +""" +function _make_cablepart!(workspace::FEMWorkspace, part::Sector, + cable_idx::Int, comp_idx::Int, comp_id::String, + phase::Int, layer_idx::Int) + + # Get the cable's center coordinates from the workspace + cable_position = workspace.problem_def.system.cables[cable_idx] + x_center = to_nominal(cable_position.horz) + y_center = to_nominal(cable_position.vert) + + # The vertices in the Sector object are already rotated and relative to the cable's origin. + # We just need to translate them to the cable's position in the system. + translated_vertices = [(v[1] + x_center, v[2] + y_center) for v in part.vertices] + + # Calculate mesh size for this part + num_elements = workspace.formulation.elements_per_length_conductor + mesh_size = _calc_mesh_size(part.radius_in, part.radius_ext, part.material_props, num_elements, workspace) + + # Create the polygon in Gmsh + surface_tag, marker = draw_polygon(translated_vertices) + + # --- The rest of this function is similar to the other _make_cablepart! methods --- + + # Get material group (1 for conductor) + material_group = get_material_group(part) + material_id = get_or_register_material_id(workspace, part.material_props) + + # Create physical tag + physical_group_tag = encode_physical_group_tag(1, cable_idx, comp_idx, material_group, material_id) + + # Create a descriptive name for the entity + elementary_name = create_cable_elementary_name( + cable_idx=cable_idx, component_id=comp_id, group_type=material_group, + part_type="sector", layer_idx=layer_idx, phase=phase + ) + + # Create the entity data and add it to the workspace's unassigned entities + core_data = CoreEntityData(physical_group_tag, elementary_name, mesh_size) + entity_data = CablePartEntity(core_data, part) + workspace.unassigned_entities[marker] = entity_data + + # Register the physical group + register_physical_group!(workspace, physical_group_tag, part.material_props) +end + + +""" +$(TYPEDSIGNATURES) + +Specialized method to create the geometry for a `SectorInsulator`. +""" +function _make_cablepart!(workspace::FEMWorkspace, part::SectorInsulator, + cable_idx::Int, comp_idx::Int, comp_id::String, + phase::Int, layer_idx::Int) + + cable_position = workspace.problem_def.system.cables[cable_idx] + x_center = to_nominal(cable_position.horz) + y_center = to_nominal(cable_position.vert) + + # Translate the vertices for both outer and inner boundaries + outer_vertices_translated = [(v[1] + x_center, v[2] + y_center) for v in part.outer_vertices] + inner_vertices_translated = [(v[1] + x_center, v[2] + y_center) for v in part.inner_sector.vertices] + + # Calculate mesh size for this part + num_elements = workspace.formulation.elements_per_length_insulator + mesh_size = _calc_mesh_size(part.radius_in, part.radius_ext, part.material_props, num_elements, workspace) + + # Create the polygon with a hole using our new drawing primitive + surface_tag, marker = draw_polygon_with_hole(outer_vertices_translated, inner_vertices_translated) + + # --- The rest is similar to the Sector method --- + + material_group = get_material_group(part) + material_id = get_or_register_material_id(workspace, part.material_props) + physical_group_tag = encode_physical_group_tag(1, cable_idx, comp_idx, material_group, material_id) + + elementary_name = create_cable_elementary_name( + cable_idx=cable_idx, component_id=comp_id, group_type=material_group, + part_type="sector_insulator", layer_idx=layer_idx, phase=phase + ) + + core_data = CoreEntityData(physical_group_tag, elementary_name, mesh_size) + entity_data = CablePartEntity(core_data, part) + workspace.unassigned_entities[marker] = entity_data + + register_physical_group!(workspace, physical_group_tag, part.material_props) end \ No newline at end of file From 15d3085f7d85e5e19cefe1388008b0b5829a1fa8 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 14:04:02 +0200 Subject: [PATCH 11/36] feat(sector-FEM): add GeometryBasics dependency for geometric operations (using Point) --- src/engine/FEM.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/engine/FEM.jl b/src/engine/FEM.jl index f3b68aa0..1fee4e3d 100644 --- a/src/engine/FEM.jl +++ b/src/engine/FEM.jl @@ -39,6 +39,7 @@ using ...Utils: display_path, set_logger!, is_headless, to_nominal using Measurements using LinearAlgebra using Colors +using GeometryBasics # FEM specific dependencies using Gmsh From 9492c87b5ceac974b8d4feaf63882997b54cb4c9 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 14:55:28 +0200 Subject: [PATCH 12/36] feat(sector-FEM): update vertex translation to use GeometryBasics.Point for Sector and SectorInsulator --- src/engine/fem/cable.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/engine/fem/cable.jl b/src/engine/fem/cable.jl index 1ef9155f..2bfb2b91 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -408,13 +408,14 @@ function _make_cablepart!(workspace::FEMWorkspace, part::Sector, # The vertices in the Sector object are already rotated and relative to the cable's origin. # We just need to translate them to the cable's position in the system. - translated_vertices = [(v[1] + x_center, v[2] + y_center) for v in part.vertices] + translated_vertices = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.vertices] # Calculate mesh size for this part num_elements = workspace.formulation.elements_per_length_conductor mesh_size = _calc_mesh_size(part.radius_in, part.radius_ext, part.material_props, num_elements, workspace) # Create the polygon in Gmsh + @debug "the translated vertices are $translated_vertices \n they are of type $(typeof(translated_vertices))" surface_tag, marker = draw_polygon(translated_vertices) # --- The rest of this function is similar to the other _make_cablepart! methods --- @@ -456,8 +457,8 @@ function _make_cablepart!(workspace::FEMWorkspace, part::SectorInsulator, y_center = to_nominal(cable_position.vert) # Translate the vertices for both outer and inner boundaries - outer_vertices_translated = [(v[1] + x_center, v[2] + y_center) for v in part.outer_vertices] - inner_vertices_translated = [(v[1] + x_center, v[2] + y_center) for v in part.inner_sector.vertices] + outer_vertices_translated = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.outer_vertices] + inner_vertices_translated = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.inner_sector.vertices] # Calculate mesh size for this part num_elements = workspace.formulation.elements_per_length_insulator From 036282c0d9b9d816d832625a44ea2c14f66edfbb Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 14:55:47 +0200 Subject: [PATCH 13/36] feat(sector-FEM): enhance draw_polygon to calculate centroid for Point vertices manually --- src/engine/fem/drawing.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index 96018db1..d1ee5720 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -641,6 +641,7 @@ function draw_transition_region( end function draw_polygon(vertices::Vector{<:Point}) + @warn "I am drawing a polygon for a Point vertices using your new defined method" # Create points points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in vertices] @@ -654,9 +655,13 @@ function draw_polygon(vertices::Vector{<:Point}) # Synchronize to make the new entity available for calculations gmsh.model.occ.synchronize() - # Calculate a center of mass for placing a marker - com = gmsh.model.get_center_of_mass(2, surface_tag) - marker = [com[1], com[2], com[3]] + # Workaround: calculate the centroid of the vertices as a marker + if isempty(vertices) + error("Cannot calculate centroid of empty vertex list.") + end + avg_x = sum(v[1] for v in vertices) / length(vertices) + avg_y = sum(v[2] for v in vertices) / length(vertices) + marker = [avg_x, avg_y, 0.0] return surface_tag, marker end From 99528a424369548cbe1f6f673ecd9b1aebcf682a Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 1 Sep 2025 17:31:13 +0200 Subject: [PATCH 14/36] fix(helper): use Base.error for GetDP executable error handling --- src/engine/fem/helpers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/engine/fem/helpers.jl b/src/engine/fem/helpers.jl index 73da7501..f2046f12 100644 --- a/src/engine/fem/helpers.jl +++ b/src/engine/fem/helpers.jl @@ -309,6 +309,6 @@ function _resolve_getdp_path(opts::NamedTuple) return fallback end - error("GetDP executable not found or not working (invocation with -info failed). " * + Base.error("GetDP executable not found or not working (invocation with -info failed). " * "Provide :getdp_executable in opts or ensure GetDP.jl is properly deployed.") end \ No newline at end of file From f67e07dc7ef0727ba9d1ca6a3685dfe69ea5fdb9 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 2 Sep 2025 08:40:07 +0200 Subject: [PATCH 15/36] feat(sector-FEM): update logging in draw_polygon function for clarity --- src/engine/fem/drawing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index d1ee5720..b575921b 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -641,7 +641,7 @@ function draw_transition_region( end function draw_polygon(vertices::Vector{<:Point}) - @warn "I am drawing a polygon for a Point vertices using your new defined method" + @debug "Drawing a polygon for a Point vertices" # Create points points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in vertices] @@ -655,7 +655,7 @@ function draw_polygon(vertices::Vector{<:Point}) # Synchronize to make the new entity available for calculations gmsh.model.occ.synchronize() - # Workaround: calculate the centroid of the vertices as a marker + # calculate the centroid of the vertices as a marker if isempty(vertices) error("Cannot calculate centroid of empty vertex list.") end From 7113f19d53f3b6e4f4e82cfcf01bacc9537d6de2 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 2 Sep 2025 09:24:52 +0200 Subject: [PATCH 16/36] feat(sector-FEM): add draw_polygon function with max edge length parameter --- src/engine/fem/cable.jl | 4 ++-- src/engine/fem/drawing.jl | 48 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 2 deletions(-) diff --git a/src/engine/fem/cable.jl b/src/engine/fem/cable.jl index 2bfb2b91..8a5f2de2 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -416,7 +416,7 @@ function _make_cablepart!(workspace::FEMWorkspace, part::Sector, # Create the polygon in Gmsh @debug "the translated vertices are $translated_vertices \n they are of type $(typeof(translated_vertices))" - surface_tag, marker = draw_polygon(translated_vertices) + surface_tag, marker = draw_polygon(translated_vertices, mesh_size) # --- The rest of this function is similar to the other _make_cablepart! methods --- @@ -463,7 +463,7 @@ function _make_cablepart!(workspace::FEMWorkspace, part::SectorInsulator, # Calculate mesh size for this part num_elements = workspace.formulation.elements_per_length_insulator mesh_size = _calc_mesh_size(part.radius_in, part.radius_ext, part.material_props, num_elements, workspace) - + # Create the polygon with a hole using our new drawing primitive surface_tag, marker = draw_polygon_with_hole(outer_vertices_translated, inner_vertices_translated) diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index b575921b..469d2c0e 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -800,4 +800,52 @@ function get_air_gap_markers(num_wires::Int, radius_wire::Number, radius_in::Num return markers end +function draw_polygon(vertices::Vector{<:Point}, max_edge_length::Number) + if isempty(vertices) + error("Cannot draw a polygon with no vertices.") + end + + new_vertices = Point[] + for i in 1:length(vertices) + p1 = vertices[i] + p2 = vertices[i % length(vertices) + 1] + + push!(new_vertices, p1) + + edge_vec = p2 - p1 + edge_len = norm(edge_vec) + + if edge_len > max_edge_length + num_segments = ceil(Int, edge_len / max_edge_length) + for j in 1:(num_segments - 1) + intermediate_point = p1 + (j / num_segments) * edge_vec + push!(new_vertices, intermediate_point) + end + end + end + + # Create points + points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in new_vertices] + + # Create lines + lines = [gmsh.model.occ.add_line(points[i], points[i % length(points) + 1]) for i in 1:length(points)] + + # Create curve loop and surface + curve_loop = gmsh.model.occ.add_curve_loop(lines) + surface_tag = gmsh.model.occ.add_plane_surface([curve_loop]) + + # Synchronize to make the new entity available for calculations + gmsh.model.occ.synchronize() + + # calculate the centroid of the original vertices as a marker + if isempty(vertices) + error("Cannot calculate centroid of empty vertex list.") + end + avg_x = sum(v[1] for v in vertices) / length(vertices) + avg_y = sum(v[2] for v in vertices) / length(vertices) + marker = [avg_x, avg_y, 0.0] + + return surface_tag, marker +end + From 9c4686c0d7d268e4b6baa261938ab2bb00d6d6fe Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 2 Sep 2025 09:46:44 +0200 Subject: [PATCH 17/36] feat(sector-FEM): enhance draw_polygon_with_hole to support max edge length parameter TODO: it is using thickness to divide! --- src/engine/fem/cable.jl | 2 +- src/engine/fem/drawing.jl | 35 ++++++++++++++++++++++++++++++++--- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/src/engine/fem/cable.jl b/src/engine/fem/cable.jl index 8a5f2de2..e7708792 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -465,7 +465,7 @@ function _make_cablepart!(workspace::FEMWorkspace, part::SectorInsulator, mesh_size = _calc_mesh_size(part.radius_in, part.radius_ext, part.material_props, num_elements, workspace) # Create the polygon with a hole using our new drawing primitive - surface_tag, marker = draw_polygon_with_hole(outer_vertices_translated, inner_vertices_translated) + surface_tag, marker = draw_polygon_with_hole(outer_vertices_translated, inner_vertices_translated, mesh_size) # --- The rest is similar to the Sector method --- diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index 469d2c0e..8472d4f2 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -678,14 +678,43 @@ Draw a polygon with a hole. # Returns - A tuple containing the Gmsh surface tag and a marker point `[x, y, z]`. """ -function draw_polygon_with_hole(outer_vertices::Vector{<:Point}, inner_vertices::Vector{<:Point}) +function draw_polygon_with_hole(outer_vertices::Vector{<:Point}, inner_vertices::Vector{<:Point}, max_edge_length::Number) + + function densify_vertices(vertices::Vector{<:Point}, max_len::Number) + new_vertices = Point[] + if isempty(vertices) + return new_vertices + end + for i in 1:length(vertices) + p1 = vertices[i] + p2 = vertices[i % length(vertices) + 1] + + push!(new_vertices, p1) + + edge_vec = p2 - p1 + edge_len = norm(edge_vec) + + if edge_len > max_len + num_segments = ceil(Int, edge_len / max_len) + for j in 1:(num_segments - 1) + intermediate_point = p1 + (j / num_segments) * edge_vec + push!(new_vertices, intermediate_point) + end + end + end + return new_vertices + end + + new_outer_vertices = densify_vertices(outer_vertices, max_edge_length) + new_inner_vertices = densify_vertices(inner_vertices, max_edge_length) + # Create outer boundary - outer_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in outer_vertices] + outer_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in new_outer_vertices] outer_lines = [gmsh.model.occ.add_line(outer_points[i], outer_points[i % length(outer_points) + 1]) for i in 1:length(outer_points)] outer_loop = gmsh.model.occ.add_curve_loop(outer_lines) # Create inner boundary (hole) - inner_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in inner_vertices] + inner_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in new_inner_vertices] inner_lines = [gmsh.model.occ.add_line(inner_points[i], inner_points[i % length(inner_points) + 1]) for i in 1:length(inner_points)] inner_loop = gmsh.model.occ.add_curve_loop(inner_lines) From a009dcfadfb9875bccf8509a13c8f100d9dc6970 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 3 Sep 2025 16:47:05 +0200 Subject: [PATCH 18/36] feat(sector-FEM): refactor _densify_vertices for better meshing of sector shaped conductor and insulator --- src/engine/fem/drawing.jl | 69 +++++++++++++++------------------------ 1 file changed, 26 insertions(+), 43 deletions(-) diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index 8472d4f2..bfab73c4 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -678,35 +678,35 @@ Draw a polygon with a hole. # Returns - A tuple containing the Gmsh surface tag and a marker point `[x, y, z]`. """ -function draw_polygon_with_hole(outer_vertices::Vector{<:Point}, inner_vertices::Vector{<:Point}, max_edge_length::Number) - - function densify_vertices(vertices::Vector{<:Point}, max_len::Number) - new_vertices = Point[] - if isempty(vertices) - return new_vertices - end - for i in 1:length(vertices) - p1 = vertices[i] - p2 = vertices[i % length(vertices) + 1] - - push!(new_vertices, p1) - - edge_vec = p2 - p1 - edge_len = norm(edge_vec) - - if edge_len > max_len - num_segments = ceil(Int, edge_len / max_len) - for j in 1:(num_segments - 1) - intermediate_point = p1 + (j / num_segments) * edge_vec - push!(new_vertices, intermediate_point) - end +function _densify_vertices(vertices::Vector{<:Point}, max_len::Number) + new_vertices = Point[] + if isempty(vertices) + return new_vertices + end + for i in 1:length(vertices) + p1 = vertices[i] + p2 = vertices[i % length(vertices) + 1] + + push!(new_vertices, p1) + + edge_vec = p2 - p1 + edge_len = norm(edge_vec) + + if edge_len > max_len + num_segments = ceil(Int, edge_len / max_len) + for j in 1:(num_segments - 1) + intermediate_point = p1 + (j / num_segments) * edge_vec + push!(new_vertices, intermediate_point) end end - return new_vertices end + return new_vertices +end - new_outer_vertices = densify_vertices(outer_vertices, max_edge_length) - new_inner_vertices = densify_vertices(inner_vertices, max_edge_length) +function draw_polygon_with_hole(outer_vertices::Vector{<:Point}, inner_vertices::Vector{<:Point}, max_edge_length::Number) + + new_outer_vertices = _densify_vertices(outer_vertices, max_edge_length) + new_inner_vertices = _densify_vertices(inner_vertices, max_edge_length) # Create outer boundary outer_points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in new_outer_vertices] @@ -834,24 +834,7 @@ function draw_polygon(vertices::Vector{<:Point}, max_edge_length::Number) error("Cannot draw a polygon with no vertices.") end - new_vertices = Point[] - for i in 1:length(vertices) - p1 = vertices[i] - p2 = vertices[i % length(vertices) + 1] - - push!(new_vertices, p1) - - edge_vec = p2 - p1 - edge_len = norm(edge_vec) - - if edge_len > max_edge_length - num_segments = ceil(Int, edge_len / max_edge_length) - for j in 1:(num_segments - 1) - intermediate_point = p1 + (j / num_segments) * edge_vec - push!(new_vertices, intermediate_point) - end - end - end + new_vertices = _densify_vertices(vertices, max_edge_length) # Create points points = [gmsh.model.occ.add_point(v[1], v[2], 0.0) for v in new_vertices] From 2c65b6e7379a18e5f6c93b54e9d975ed5f086a85 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Thu, 4 Sep 2025 15:54:02 +0200 Subject: [PATCH 19/36] chore(ci): remove main branch trigger from CI workflow to make it test on my feature branch too --- .github/workflows/CI.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4cd77cac..d8130684 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,8 +1,6 @@ name: CI on: push: - branches: - - main tags: ['*'] pull_request: workflow_dispatch: From 1951c24c57a03ba0f16fbc4ac5c9ac6c7a76882d Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Thu, 4 Sep 2025 15:59:09 +0200 Subject: [PATCH 20/36] chore(ci): add gh-pages to branches-ignore in CI workflow but make it run CI on anything else --- .github/workflows/CI.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d8130684..23ce0a83 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,6 +1,8 @@ name: CI on: push: + branches-ignore: + - gh-pages tags: ['*'] pull_request: workflow_dispatch: From 78bf0082c08b5732d52b91b5f36db25987cad64e Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 9 Sep 2025 23:01:29 +0200 Subject: [PATCH 21/36] refactor(fem): add GeometryBasics dependency for enhanced geometric operations - needed for sector shaped cable --- src/engine/fem/FEM.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/engine/fem/FEM.jl b/src/engine/fem/FEM.jl index 5759a593..6533145c 100644 --- a/src/engine/fem/FEM.jl +++ b/src/engine/fem/FEM.jl @@ -43,6 +43,7 @@ using ...Utils: using Measurements using LinearAlgebra using Colors +using GeometryBasics # FEM specific dependencies using Gmsh From 9ef8c465111f12f9c31e9ca1b2a4ac4148c5fa3e Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 14 Sep 2025 03:01:43 +0200 Subject: [PATCH 22/36] feat(sector): add insulation thickness and enhance geometry calculations with debug logging --- src/datamodel/sector.jl | 63 ++++++++++++++++++++++++++++------------- 1 file changed, 44 insertions(+), 19 deletions(-) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index c5066566..a2e32c88 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -16,6 +16,8 @@ struct SectorParams{T<:REALSCALAR} r_corner::T "Angular width of the conductor's flat base/sides [degrees]." theta_cond_deg::T + "Insulation thickness \\[m\\]." + d_insulation::T # needed to correct for the offset end """ @@ -65,7 +67,7 @@ function Sector( # 3. Calculate cross-sectional area using the Shoelace formula cross_section = PolygonOps.area(rotated_vertices) - + @debug "Sector cross-sectional area: $(cross_section*1e6) mm²" # 4. Calculate DC resistance rho_eff = calc_temperature_correction(material_props.alpha, temperature, material_props.T0) * material_props.rho resistance = rho_eff / cross_section @@ -94,18 +96,28 @@ end function _calculate_sector_geometry(p::SectorParams) phi_deg = 90.0 - p.theta_cond_deg / 2.0 - phi_rad = deg2rad(phi_deg) + phi_rad = deg2rad(phi_deg) # ϕ + @debug "(Sector) phi_deg: $phi_deg" + @debug "(Sector) phi_rad: $phi_rad rad" if abs(cos(phi_rad)) < 1e-9 error("theta_cond_deg is too close to 180, leading to division by zero. Check parameters.") end - d_base_corner = p.r_corner * (1.0 / cos(phi_rad) - 1.0) - d_offset = p.r_back - p.d_sector - d_base_corner + d_base_corner = p.r_corner * (1.0 / cos(phi_rad) - 1.0) # D_B + @debug "(Sector) d_base_corner (D_B) : $d_base_corner m" + d_offset = p.r_back - p.d_sector - d_base_corner # D_O + @debug "(Sector) d_offset (D_O) : $d_offset m" + d_insulation_offset = (p.d_insulation / (cos((pi - ((2 * pi) / p.n_sectors)) / 2.0))) - d_offset # D_I + @debug "(Sector) d_insulation_offset (D_I) : $d_insulation_offset m" + + # trick test (works though different that Urquhart's) + #d_offset = (p.d_insulation / (cos((pi - ((2 * pi) / p.n_sectors)) / 2.0))) # D_I + #d_insulation_offset = 0.0 - x_base_corner = p.r_corner * sin(phi_rad) - y_base_corner = x_base_corner * tan(phi_rad) + d_offset - node_A = Point2f(x_base_corner, y_base_corner) + x_base_corner = p.r_corner * sin(phi_rad) # X_A = -X_F + y_base_corner = x_base_corner * tan(phi_rad) + d_offset # Y_A = Y_F + node_A = GeometryBasics.Point2f(x_base_corner, y_base_corner + d_insulation_offset) k = p.r_corner / cos(phi_rad) + d_offset qa = 1.0 + tan(phi_rad)^2 @@ -117,29 +129,37 @@ function _calculate_sector_geometry(p::SectorParams) error("Cannot calculate side corner center: negative discriminant ($discriminant). Check parameters.") end - x_side_center = (-qb + sqrt(discriminant)) / (2.0 * qa) - y_side_center = x_side_center * tan(phi_rad) + k + x_side_center = (-qb + sqrt(discriminant)) / (2.0 * qa) # X_N + y_side_center = x_side_center * tan(phi_rad) + k # Y_N - x_side_lower = x_side_center + p.r_corner * sin(phi_rad) - y_side_lower = y_side_center - p.r_corner * cos(phi_rad) - node_B = Point2f(x_side_lower, y_side_lower) + x_side_lower = x_side_center + p.r_corner * sin(phi_rad) # X_B + y_side_lower = y_side_center - p.r_corner * cos(phi_rad) # Y_B + node_B = GeometryBasics.Point2f(x_side_lower, y_side_lower + d_insulation_offset) dist_origin_side_center = sqrt(x_side_center^2 + y_side_center^2) if dist_origin_side_center < 1e-9 error("Side corner center is at the origin, cannot determine upper point direction.") end - x_side_upper = x_side_center * p.r_back / dist_origin_side_center - y_side_upper = y_side_center * p.r_back / dist_origin_side_center - node_C = Point2f(x_side_upper, y_side_upper) + x_side_upper = x_side_center * p.r_back / dist_origin_side_center # X_C + y_side_upper = y_side_center * p.r_back / dist_origin_side_center # Y_C + node_C = GeometryBasics.Point2f(x_side_upper, y_side_upper + d_insulation_offset) nodes = ( - A=node_A, B=node_B, C=node_C, D=Point2f(-node_C[1], node_C[2]), - E=Point2f(-node_B[1], node_B[2]), F=Point2f(-node_A[1], node_A[2]) + A=node_A, + B=node_B, + C=node_C, + D=GeometryBasics.Point2f(-node_C[1], node_C[2]), + E=GeometryBasics.Point2f(-node_B[1], node_B[2]), + F=GeometryBasics.Point2f(-node_A[1], node_A[2]) ) + @debug "(Sector) Nodes: $nodes" centers = ( - Back=Point2f(0, 0), Base=Point2f(0, d_offset+ p.r_corner / cos(phi_rad)), - RightSide=Point2f(x_side_center, y_side_center), LeftSide=Point2f(-x_side_center, y_side_center) + Back=GeometryBasics.Point2f(0, 0+ d_insulation_offset), + Base=GeometryBasics.Point2f(0, d_offset + p.r_corner / cos(phi_rad) + d_insulation_offset), + RightSide=GeometryBasics.Point2f(x_side_center, y_side_center+ d_insulation_offset), + LeftSide=GeometryBasics.Point2f(-x_side_center, y_side_center+ d_insulation_offset) ) + @debug "(Sector) Centers: $centers" return (Nodes=nodes, Centers=centers, Params=p) end @@ -155,6 +175,7 @@ function _generate_arc_points(center, radius, start_angle, end_angle, num_points @debug "(Sector) Adjusted end_angle to be less than start_angle: $(rad2deg(end_angle)) < $(rad2deg(start_angle))" end angle_range = range(start_angle, stop=end_angle, length=num_points) + @debug "(Sector) ______________________________ ∠ $(rad2deg(start_angle-end_angle))." return [Point2f(center[1] + radius * cos(a), center[2] + radius * sin(a)) for a in angle_range] end @@ -168,7 +189,9 @@ function _calculate_sector_polygon_points(params; num_arc_points=20) # increase push!(poly_points, nodes.F) if params.r_corner > 1e-9 start_angle = get_angle(nodes.F, centers.Base) + @debug "Arc from F to A: start_angle=$(rad2deg(start_angle))" end_angle = get_angle(nodes.A, centers.Base) + @debug "Arc from F to A: end_angle=$(rad2deg(end_angle))" append!(poly_points, _generate_arc_points(centers.Base, params.r_corner, start_angle, end_angle, num_arc_points)[2:end]) else @@ -191,7 +214,9 @@ function _calculate_sector_polygon_points(params; num_arc_points=20) # increase # Arc C to D (Back) start_angle = get_angle(nodes.C, centers.Back) + @debug "Arc from C to D: start_angle=$(rad2deg(start_angle))" end_angle = get_angle(nodes.D, centers.Back) + @debug "Arc from C to D: end_angle=$(rad2deg(end_angle))" append!(poly_points, _generate_arc_points(centers.Back, params.r_back, start_angle, end_angle, num_arc_points)[2:end]) # Arc D to E (Left Side) From 4675bc063d576d0f259f15639af3db366703602b Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 14 Sep 2025 13:57:22 +0200 Subject: [PATCH 23/36] feat(sector): implement Shoelace formula for area calculations and ensure polygon closure in sector geometry --- src/datamodel/sector.jl | 44 ++++++++++++++++++++++++++++++-- src/datamodel/sectorinsulator.jl | 38 ++++++++++++++++++++++++++- 2 files changed, 79 insertions(+), 3 deletions(-) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index a2e32c88..12cc092e 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -65,9 +65,22 @@ function Sector( rotation_angle_rad = deg2rad(rotation_angle_deg) rotated_vertices = [_rotate_point(p, rotation_angle_rad) for p in base_vertices] + # Ensure the polygon is closed: first point == last point (within tolerance) + # tol = 1e-9 + # if !isempty(rotated_vertices) + # firstp = rotated_vertices[1] + # lastp = rotated_vertices[end] + # if !(isapprox(firstp[1], lastp[1]; atol=tol, rtol=0.0) && + # isapprox(firstp[2], lastp[2]; atol=tol, rtol=0.0)) + # push!(rotated_vertices, firstp) + # @debug "(Sector) rotated_vertices not closed — appended first point to close polygon." + # end + # end # 3. Calculate cross-sectional area using the Shoelace formula - cross_section = PolygonOps.area(rotated_vertices) - @debug "Sector cross-sectional area: $(cross_section*1e6) mm²" + #cross_section = PolygonOps.area(rotated_vertices) + cross_section = _shoelace_area(rotated_vertices) + #@debug "Sector cross-sectional area: $(cross_section*1e6) mm²" + @debug "Sector cross-sectional area (Shoelace): $(shoe_lace_area*1e6) mm²" # 4. Calculate DC resistance rho_eff = calc_temperature_correction(material_props.alpha, temperature, material_props.T0) * material_props.rho resistance = rho_eff / cross_section @@ -130,8 +143,15 @@ function _calculate_sector_geometry(p::SectorParams) end x_side_center = (-qb + sqrt(discriminant)) / (2.0 * qa) # X_N + @debug "(Sector) x_side_center (X_N): $(x_side_center*1e3) mm" + @debug "(Sector) r_corner: $(p.r_corner*1e3) mm" y_side_center = x_side_center * tan(phi_rad) + k # Y_N + # Sector width (for checks): + # w = 2 * X_N + 2 * r_corner + w_sector = 2.0 * x_side_center + 2.0 * p.r_corner + @debug "(Sector) sector width (computed): $(w_sector) m ($(w_sector*1e3) mm)" + x_side_lower = x_side_center + p.r_corner * sin(phi_rad) # X_B y_side_lower = y_side_center - p.r_corner * cos(phi_rad) # Y_B node_B = GeometryBasics.Point2f(x_side_lower, y_side_lower + d_insulation_offset) @@ -236,4 +256,24 @@ end function _rotate_point(p::Point2f, angle_rad::Real) cos_a, sin_a = cos(angle_rad), sin(angle_rad) return Point2f(p[1] * cos_a - p[2] * sin_a, p[1] * sin_a + p[2] * cos_a) +end + +""" +Calculates the area of a polygon using the Shoelace formula. +The vertices are given as a vector of points. +""" +function _shoelace_area(vertices::Vector{Point{2,T}}) where {T} + n = length(vertices) + if n < 3 + return zero(T) + end + + area = zero(T) + for i in 1:n + p1 = vertices[i] + p2 = vertices[mod1(i + 1, n)] # Wrap around for the last segment + area += (p1[1] * p2[2] - p2[1] * p1[2]) + end + + return abs(area) / 2.0 end \ No newline at end of file diff --git a/src/datamodel/sectorinsulator.jl b/src/datamodel/sectorinsulator.jl index bfe324f9..3dc7422a 100644 --- a/src/datamodel/sectorinsulator.jl +++ b/src/datamodel/sectorinsulator.jl @@ -41,7 +41,24 @@ function SectorInsulator( # 2. Calculate areas inner_area = inner_sector.cross_section - outer_area = PolygonOps.area(outer_vertices) + + # tol = 1e-9 + # if !isempty(outer_vertices) + # firstp = outer_vertices[1] + # lastp = outer_vertices[end] + # if !(isapprox(firstp[1], lastp[1]; atol=tol, rtol=0.0) && + # isapprox(firstp[2], lastp[2]; atol=tol, rtol=0.0)) + # push!(outer_vertices, firstp) + # @debug "(Sector) outer_vertices not closed — appended first point to close polygon." + # end + # end + + + + #outer_area = PolygonOps.area(outer_vertices) + outer_area = _shoelace_area(outer_vertices) + @debug "SectorInsulator inner area: $(inner_area*1e6) mm²" + @debug "SectorInsulator outer area: $(outer_area*1e6) mm²" cross_section = outer_area - inner_area # 3. Approximate capacitance and conductance using equivalent coaxial circles @@ -112,3 +129,22 @@ function _calculate_offset_polygon(vertices::Vector{Point{2, T}}, thickness::T) return new_vertices end +""" +Calculates the area of a polygon using the Shoelace formula. +The vertices are given as a vector of points. +""" +function _shoelace_area(vertices::Vector{Point{2,T}}) where {T} + n = length(vertices) + if n < 3 + return zero(T) + end + + area = zero(T) + for i in 1:n + p1 = vertices[i] + p2 = vertices[mod1(i + 1, n)] # Wrap around for the last segment + area += (p1[1] * p2[2] - p2[1] * p1[2]) + end + + return abs(area) / 2.0 +end \ No newline at end of file From 5b684fa6ce93b91670bb25783c4ea83d49a77000 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Fri, 19 Sep 2025 14:24:26 +0200 Subject: [PATCH 24/36] refactor(sector): correct debug logging for cross-sectional area calculation and increase arc points for improved accuracy --- src/datamodel/sector.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index 12cc092e..a98bc7a5 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -80,7 +80,7 @@ function Sector( #cross_section = PolygonOps.area(rotated_vertices) cross_section = _shoelace_area(rotated_vertices) #@debug "Sector cross-sectional area: $(cross_section*1e6) mm²" - @debug "Sector cross-sectional area (Shoelace): $(shoe_lace_area*1e6) mm²" + @debug "Sector cross-sectional area (Shoelace): $(cross_section*1e6) mm²" # 4. Calculate DC resistance rho_eff = calc_temperature_correction(material_props.alpha, temperature, material_props.T0) * material_props.rho resistance = rho_eff / cross_section @@ -199,7 +199,7 @@ function _generate_arc_points(center, radius, start_angle, end_angle, num_points return [Point2f(center[1] + radius * cos(a), center[2] + radius * sin(a)) for a in angle_range] end -function _calculate_sector_polygon_points(params; num_arc_points=20) # increase `num_arc_points` for higher accuracy +function _calculate_sector_polygon_points(params; num_arc_points=30) # increase `num_arc_points` for higher accuracy geom = _calculate_sector_geometry(params) nodes, centers = geom.Nodes, geom.Centers From ad5724ae30697d1f9dddafe0aa0c17938b99e533 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 21 Sep 2025 01:27:26 +0200 Subject: [PATCH 25/36] refactor(sector): remove GeometryBasics dependency and use Makie for point definitions --- Project.toml | 2 -- src/datamodel/DataModel.jl | 4 +--- src/datamodel/sector.jl | 38 ++++++++++++++------------------ src/datamodel/sectorinsulator.jl | 19 ---------------- src/engine/fem/FEM.jl | 3 +-- src/engine/fem/cable.jl | 6 ++--- 6 files changed, 22 insertions(+), 50 deletions(-) diff --git a/Project.toml b/Project.toml index 1682382e..97968a1c 100644 --- a/Project.toml +++ b/Project.toml @@ -18,7 +18,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" EzXML = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615" ForceImport = "9dda63f9-cce7-5873-89fa-eccbb2fffcde" Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196" -GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" GetDP = "dbfd83e6-7aba-450b-9c2b-93ccd973023a" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" @@ -63,7 +62,6 @@ DocStringExtensions = "0.9.3" EzXML = "1.2.0" ForceImport = "0.0.3" Functors = "0.5.2" -GeometryBasics = "0.5.10" GLMakie = "0.13" Gmsh = "0.3.1" JSON3 = "1.14.2" diff --git a/src/datamodel/DataModel.jl b/src/datamodel/DataModel.jl index 95778079..68e48b1b 100644 --- a/src/datamodel/DataModel.jl +++ b/src/datamodel/DataModel.jl @@ -50,10 +50,8 @@ using DataFrames using Colors using Plots using DisplayAs: DisplayAs -using GeometryBasics -using PolygonOps using LinearAlgebra - +using Makie: Point, Point2f # otherwise will require adding GeometryBasics as a dependency # Abstract types & interfaces include("types.jl") include("radii.jl") diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index a98bc7a5..f38b25c3 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -130,7 +130,7 @@ function _calculate_sector_geometry(p::SectorParams) x_base_corner = p.r_corner * sin(phi_rad) # X_A = -X_F y_base_corner = x_base_corner * tan(phi_rad) + d_offset # Y_A = Y_F - node_A = GeometryBasics.Point2f(x_base_corner, y_base_corner + d_insulation_offset) + node_A = Point2f(x_base_corner, y_base_corner + d_insulation_offset) k = p.r_corner / cos(phi_rad) + d_offset qa = 1.0 + tan(phi_rad)^2 @@ -154,7 +154,7 @@ function _calculate_sector_geometry(p::SectorParams) x_side_lower = x_side_center + p.r_corner * sin(phi_rad) # X_B y_side_lower = y_side_center - p.r_corner * cos(phi_rad) # Y_B - node_B = GeometryBasics.Point2f(x_side_lower, y_side_lower + d_insulation_offset) + node_B = Point2f(x_side_lower, y_side_lower + d_insulation_offset) dist_origin_side_center = sqrt(x_side_center^2 + y_side_center^2) if dist_origin_side_center < 1e-9 @@ -162,22 +162,22 @@ function _calculate_sector_geometry(p::SectorParams) end x_side_upper = x_side_center * p.r_back / dist_origin_side_center # X_C y_side_upper = y_side_center * p.r_back / dist_origin_side_center # Y_C - node_C = GeometryBasics.Point2f(x_side_upper, y_side_upper + d_insulation_offset) + node_C = Point2f(x_side_upper, y_side_upper + d_insulation_offset) nodes = ( A=node_A, B=node_B, C=node_C, - D=GeometryBasics.Point2f(-node_C[1], node_C[2]), - E=GeometryBasics.Point2f(-node_B[1], node_B[2]), - F=GeometryBasics.Point2f(-node_A[1], node_A[2]) + D=Point2f(-node_C[1], node_C[2]), + E=Point2f(-node_B[1], node_B[2]), + F=Point2f(-node_A[1], node_A[2]) ) @debug "(Sector) Nodes: $nodes" centers = ( - Back=GeometryBasics.Point2f(0, 0+ d_insulation_offset), - Base=GeometryBasics.Point2f(0, d_offset + p.r_corner / cos(phi_rad) + d_insulation_offset), - RightSide=GeometryBasics.Point2f(x_side_center, y_side_center+ d_insulation_offset), - LeftSide=GeometryBasics.Point2f(-x_side_center, y_side_center+ d_insulation_offset) + Back=Point2f(0, 0+ d_insulation_offset), + Base=Point2f(0, d_offset + p.r_corner / cos(phi_rad) + d_insulation_offset), + RightSide=Point2f(x_side_center, y_side_center+ d_insulation_offset), + LeftSide=Point2f(-x_side_center, y_side_center+ d_insulation_offset) ) @debug "(Sector) Centers: $centers" return (Nodes=nodes, Centers=centers, Params=p) @@ -258,22 +258,18 @@ function _rotate_point(p::Point2f, angle_rad::Real) return Point2f(p[1] * cos_a - p[2] * sin_a, p[1] * sin_a + p[2] * cos_a) end -""" -Calculates the area of a polygon using the Shoelace formula. -The vertices are given as a vector of points. -""" -function _shoelace_area(vertices::Vector{Point{2,T}}) where {T} - n = length(vertices) + +function _shoelace_area(vertices::AbstractVector{Point{2,T}}) where {T<:Real} + n::Int = length(vertices) if n < 3 + @warn "Polygon must have at least 3 vertices to compute area. Returning 0 area" return zero(T) end - - area = zero(T) + area::T = zero(T) for i in 1:n p1 = vertices[i] - p2 = vertices[mod1(i + 1, n)] # Wrap around for the last segment + p2 = vertices[mod1(i + 1, n)] area += (p1[1] * p2[2] - p2[1] * p1[2]) end - - return abs(area) / 2.0 + return abs(area) / T(2) end \ No newline at end of file diff --git a/src/datamodel/sectorinsulator.jl b/src/datamodel/sectorinsulator.jl index 3dc7422a..d1938af7 100644 --- a/src/datamodel/sectorinsulator.jl +++ b/src/datamodel/sectorinsulator.jl @@ -128,23 +128,4 @@ function _calculate_offset_polygon(vertices::Vector{Point{2, T}}, thickness::T) end return new_vertices -end -""" -Calculates the area of a polygon using the Shoelace formula. -The vertices are given as a vector of points. -""" -function _shoelace_area(vertices::Vector{Point{2,T}}) where {T} - n = length(vertices) - if n < 3 - return zero(T) - end - - area = zero(T) - for i in 1:n - p1 = vertices[i] - p2 = vertices[mod1(i + 1, n)] # Wrap around for the last segment - area += (p1[1] * p2[2] - p2[1] * p1[2]) - end - - return abs(area) / 2.0 end \ No newline at end of file diff --git a/src/engine/fem/FEM.jl b/src/engine/fem/FEM.jl index 6533145c..4522c1ee 100644 --- a/src/engine/fem/FEM.jl +++ b/src/engine/fem/FEM.jl @@ -43,8 +43,7 @@ using ...Utils: using Measurements using LinearAlgebra using Colors -using GeometryBasics - +using Makie: Point, Point2f # otherwise will require adding GeometryBasics as a dependency # FEM specific dependencies using Gmsh using GetDP diff --git a/src/engine/fem/cable.jl b/src/engine/fem/cable.jl index 48ce5411..73ec50e7 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -465,7 +465,7 @@ function _make_cablepart!(workspace::FEMWorkspace, part::Sector, # The vertices in the Sector object are already rotated and relative to the cable's origin. # We just need to translate them to the cable's position in the system. - translated_vertices = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.vertices] + translated_vertices = [Point(v[1] + x_center, v[2] + y_center) for v in part.vertices] # Calculate mesh size for this part num_elements = workspace.formulation.elements_per_length_conductor @@ -514,8 +514,8 @@ function _make_cablepart!(workspace::FEMWorkspace, part::SectorInsulator, y_center = to_nominal(cable_position.vert) # Translate the vertices for both outer and inner boundaries - outer_vertices_translated = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.outer_vertices] - inner_vertices_translated = [GeometryBasics.Point(v[1] + x_center, v[2] + y_center) for v in part.inner_sector.vertices] + outer_vertices_translated = [Point(v[1] + x_center, v[2] + y_center) for v in part.outer_vertices] + inner_vertices_translated = [Point(v[1] + x_center, v[2] + y_center) for v in part.inner_sector.vertices] # Calculate mesh size for this part num_elements = workspace.formulation.elements_per_length_insulator From 4717edaa6640988280600966d39b3acdb03bb376 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Sun, 21 Sep 2025 01:52:18 +0200 Subject: [PATCH 26/36] feat(sector): add support for rendering Sector and SectorInsulator layers with Makie --- src/datamodel/preview.jl | 68 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index 66fe5047..13342b1c 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -375,6 +375,74 @@ function _plot_layer_makie!(ax, layer, label::String; return plots end + if layer isa Sector + vertices = layer.vertices + # Convert vertices to Makie.Point2f format with offset + makie_points = [Makie.Point2f(v[1] + x0, v[2] + y0) for v in vertices] + # Ensure polygon is closed by adding first point at the end if needed + if length(makie_points) > 0 && makie_points[1] != makie_points[end] + push!(makie_points, makie_points[1]) + end + + color = get_material_color_makie(layer.material_props) + + poly = Makie.poly!(ax, makie_points; + color = color, + strokecolor = :black, + strokewidth = 0.5, + label = display_legend ? label : "") + + if legend_sink !== nothing && display_legend + push!(legend_sink[1], poly) + push!(legend_sink[2], label) + if length(legend_sink) >= 3 + push!(legend_sink[3], [poly]) + end + if length(legend_sink) >= 4 + push!(legend_sink[4], NaN) # not a wirearray + end + end + return (poly,) + end + + if layer isa SectorInsulator + outer_vertices = [(v[1] + x0, v[2] + y0) for v in layer.outer_vertices] + # Convert to Makie.Point2f format + outer_points = [Makie.Point2f(v[1], v[2]) for v in outer_vertices] + # Ensure polygon is closed + if length(outer_points) > 0 && outer_points[1] != outer_points[end] + push!(outer_points, outer_points[1]) + end + + # (Not used for now) The inner boundary is the conductor's vertices. It must be reversed for the hole to be drawn correctly. + inner_vertices = [(v[1] + x0, v[2] + y0) for v in layer.inner_sector.vertices] + inner_points = [Makie.Point2f(v[1], v[2]) for v in inner_vertices] + # Ensure inner polygon is closed + if length(inner_points) > 0 && inner_points[1] != inner_points[end] + push!(inner_points, inner_points[1]) + end + color = get_material_color_makie(layer.material_props) + # Create a shape with a hole by passing the outer boundary and holes as a vector of vectors + polygon_with_hole = Makie.Polygon(outer_points, [inner_points]) + poly = Makie.poly!(ax, polygon_with_hole; + color = color, + strokecolor = :black, + strokewidth = 0.5, + label = display_legend ? label : "") + + if legend_sink !== nothing && display_legend + push!(legend_sink[1], poly) + push!(legend_sink[2], label) + if length(legend_sink) >= 3 + push!(legend_sink[3], [poly]) + end + if length(legend_sink) >= 4 + push!(legend_sink[4], NaN) # not a wirearray + end + end + return (poly,) + end + @warn "Unknown layer type $(typeof(layer)); skipping" return () end From 66d5610d6d04d11b92c4ae52a5ffa673076eb5f4 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Mon, 22 Sep 2025 11:33:23 +0200 Subject: [PATCH 27/36] feat(sector): add comprehensive tutorial for sector-shaped cable design along with tests for CI --- examples/tutorial2_sector.jl | 157 ++++++++++++++++++++++++++++++ test/test_tutorial_2_sector.jl | 168 +++++++++++++++++++++++++++++++++ 2 files changed, 325 insertions(+) create mode 100644 examples/tutorial2_sector.jl create mode 100644 test/test_tutorial_2_sector.jl diff --git a/examples/tutorial2_sector.jl b/examples/tutorial2_sector.jl new file mode 100644 index 00000000..b020c1ef --- /dev/null +++ b/examples/tutorial2_sector.jl @@ -0,0 +1,157 @@ +#= +# Tutorial 2a - Building a sector-shaped cable design + +This tutorial demonstrates how to model a typical low-voltage three-core power cable with sector-shaped conductors +using the [`LineCableModels.jl`](@ref) package. The objective is to build a complete representation of a three-core 1 kV cable with 95 mm² aluminum sector-shaped conductors and a concentric copper neutral. +=# + +#= +**Tutorial outline** +```@contents +Pages = [ + "tutorial2_sector.md", +] +Depth = 2:3 +``` +=# + +#= +## Introduction + +Three-core power cables with sector-shaped conductors are common in low-voltage distribution networks. Their compact design allows for efficient use of space. This tutorial will guide you through creating a detailed [`CableDesign`](@ref) for such a cable. + +This tutorial covers: + +1. Defining materials with corrected resistivity. +2. Creating sector-shaped conductors using [`SectorParams`](@ref) and [`Sector`](@ref). +3. Assembling a multi-core [`CableDesign`](@ref). +4. Modeling a concentric neutral wire array. +5. Previewing the final cable design. +=# + +#= +## Getting started +=# + +# Load the package and set up the environment: +using LineCableModels +using DataFrames +import LineCableModels.BackendHandler: renderfig #hide +fullfile(filename) = joinpath(@__DIR__, filename); #hide +set_verbosity!(0); #hide + +#= +## Cable and Material Data + +We start by defining the materials. We will create a custom aluminum material with a resistivity corrected based on its nominal DC resistance, and a PVC material for insulation. +=# + +# Initialize materials library and add a PVC material +materials = MaterialsLibrary(add_defaults=true) +pvc = Material(Inf, 8.0, 1.0, 20.0, 0.1) # simple PVC +add!(materials, "pvc", pvc) +copper = get(materials, "copper") +aluminum = get(materials, "aluminum") + + +#= +## Sector-Shaped Core Conductors + +The core of the cable consists of three sector-shaped aluminum conductors. We define the geometry using `SectorParams` based on datasheet or standard values. +=# + +# === Sector (core) geometry (table data) === +# Based on Urquhart's paper for a 3-core 95mm^2 cable +n_sectors = 3 +r_back_mm = 10.24 # sector radius b +d_sector_mm = 9.14 # sector depth s +r_corner_mm = 1.02 # corner radius c +theta_cond_deg = 119.0 # sector angle φ +ins_thick = 1.1e-3 # core insulation thickness + +sector_params = SectorParams( + n_sectors, + r_back_mm / 1000, + d_sector_mm / 1000, + r_corner_mm / 1000, + theta_cond_deg, + ins_thick +) + +#= +With the sector parameters defined, we can create the individual `Sector` conductors and their insulation. Each sector is rotated to form the 3-core bundle. +=# + +rot_angles = (0.0, 120.0, 240.0) +sectors = [Sector(sector_params, ang, aluminum) for ang in rot_angles] +insulators = [SectorInsulator(sectors[i], ins_thick, pvc) for i in 1:3] + +components = [ + CableComponent("core1", ConductorGroup(sectors[1]), InsulatorGroup(insulators[1])), + CableComponent("core2", ConductorGroup(sectors[2]), InsulatorGroup(insulators[2])), + CableComponent("core3", ConductorGroup(sectors[3]), InsulatorGroup(insulators[3])) +] + +#= +## Concentric Neutral and Outer Jacket + +The cable includes a concentric neutral conductor made of copper wires and an outer PVC jacket. +=# + +# === Concentric neutral (30 wires) === +n_neutral = 30 +r_strand = 0.79e-3 +R_N = 14.36e-3 # radius to center of neutral wires +R_O = 17.25e-3 # outer radius of the cable + +inner_radius_neutral = R_N - r_strand +outer_jacket_thickness = R_O - (R_N + r_strand) + +neutral_wires = WireArray( + inner_radius_neutral, + Diameter(2*r_strand), + n_neutral, + 0.0, # lay ratio + copper +) + +neutral_jacket = Insulator(neutral_wires, Thickness(outer_jacket_thickness), pvc) +neutral_component = CableComponent("neutral", ConductorGroup(neutral_wires), InsulatorGroup(neutral_jacket)) + +#= +## Assembling the Cable Design + +Now we assemble the complete `CableDesign` by adding all the components. +=# + +design = CableDesign("NAYCWY_O_3x95_30x2_5", components[1]) +add!(design, components[2]) +add!(design, components[3]) +add!(design, neutral_component) + +#= +## Examining the Cable Design + +We can now display a summary of the cable design and preview it graphically. +=# + +println("Cable design summary:") +detailed_df = DataFrame(design, :detailed) +display(detailed_df) + +println("Previewing cable design...") +plt, _ = preview(design) +plt #hide + +#= +## Conclusion + +This tutorial has demonstrated how to model a three-core cable with sector-shaped conductors. Key takeaways include: + +1. Creating custom materials with corrected properties. +2. Defining complex conductor shapes like sectors. +3. Assembling a multi-core cable design component by component. +4. Visualizing the final design for verification. + +This detailed modeling capability allows for accurate analysis of various cable configurations. +=# diff --git a/test/test_tutorial_2_sector.jl b/test/test_tutorial_2_sector.jl new file mode 100644 index 00000000..65e69abf --- /dev/null +++ b/test/test_tutorial_2_sector.jl @@ -0,0 +1,168 @@ +@testitem "examples/tutorial2_sector.jl tests" setup = [defaults] begin + # Replicate the setup from the tutorial + + # === Materials === + materials = MaterialsLibrary(add_defaults=true) + pvc = Material(Inf, 8.0, 1.0, 20.0, 0.1) + add!(materials, "pvc", pvc) + copper = get(materials, "copper") + aluminum = get(materials, "aluminum") + + @testset "Material setup" begin + @test get(materials, "pvc") isa LineCableModels.Materials.Material + @test get(materials, "aluminum") isa LineCableModels.Materials.Material + @test get(materials, "copper") isa LineCableModels.Materials.Material + end + + # === Sector (core) geometry === + @testset "Sector core construction" begin + n_sectors = 3 + r_back_mm = 10.24 + d_sector_mm = 9.14 + r_corner_mm = 1.02 + theta_cond_deg = 119.0 + ins_thick = 1.1e-3 + + sector_params = SectorParams( + n_sectors, + r_back_mm / 1000, + d_sector_mm / 1000, + r_corner_mm / 1000, + theta_cond_deg, + ins_thick + ) + + rot_angles = (0.0, 120.0, 240.0) + sectors = [Sector(sector_params, ang, aluminum) for ang in rot_angles] + insulators = [SectorInsulator(sectors[i], ins_thick, pvc) for i in 1:3] + + @test length(sectors) == 3 + @test all(s -> s isa Sector, sectors) + @test length(insulators) == 3 + @test all(i -> i isa SectorInsulator, insulators) + + components = [ + CableComponent("core1", ConductorGroup(sectors[1]), InsulatorGroup(insulators[1])), + CableComponent("core2", ConductorGroup(sectors[2]), InsulatorGroup(insulators[2])), + CableComponent("core3", ConductorGroup(sectors[3]), InsulatorGroup(insulators[3])) + ] + @test length(components) == 3 + @test components[1].id == "core1" + end + + # === Concentric neutral === + @testset "Concentric neutral construction" begin + n_neutral = 30 + r_strand = 0.79e-3 + R_N = 14.36e-3 + R_O = 17.25e-3 + + inner_radius_neutral = R_N - r_strand + outer_jacket_thickness = R_O - (R_N + r_strand) + + neutral_wires = WireArray( + inner_radius_neutral, + Diameter(2*r_strand), + n_neutral, + 0.0, + copper + ) + @test neutral_wires isa WireArray + + neutral_jacket = Insulator(neutral_wires, Thickness(outer_jacket_thickness), pvc) + @test neutral_jacket isa Insulator + + neutral_component = CableComponent("neutral", ConductorGroup(neutral_wires), InsulatorGroup(neutral_jacket)) + @test neutral_component.id == "neutral" + end + + # === Assemble cable design === + @testset "Full cable design assembly" begin + # Re-create components for this testset to be self-contained + n_sectors = 3 + r_back_mm = 10.24 + d_sector_mm = 9.14 + r_corner_mm = 1.02 + theta_cond_deg = 119.0 + ins_thick = 1.1e-3 + sector_params = SectorParams(n_sectors, r_back_mm/1000, d_sector_mm/1000, r_corner_mm/1000, theta_cond_deg, ins_thick) + rot_angles = (0.0, 120.0, 240.0) + sectors = [Sector(sector_params, ang, aluminum) for ang in rot_angles] + insulators = [SectorInsulator(sectors[i], ins_thick, pvc) for i in 1:3] + components = [ + CableComponent("core1", ConductorGroup(sectors[1]), InsulatorGroup(insulators[1])), + CableComponent("core2", ConductorGroup(sectors[2]), InsulatorGroup(insulators[2])), + CableComponent("core3", ConductorGroup(sectors[3]), InsulatorGroup(insulators[3])) + ] + + n_neutral = 30 + r_strand = 0.79e-3 + R_N = 14.36e-3 + R_O = 17.25e-3 + inner_radius_neutral = R_N - r_strand + outer_jacket_thickness = R_O - (R_N + r_strand) + neutral_wires = WireArray(inner_radius_neutral, Diameter(2*r_strand), n_neutral, 0.0, copper) + neutral_jacket = Insulator(neutral_wires, Thickness(outer_jacket_thickness), pvc) + neutral_component = CableComponent("neutral", ConductorGroup(neutral_wires), InsulatorGroup(neutral_jacket)) + + design = CableDesign("NAYCWY_O_3x95_30x2_5", components[1]) + add!(design, components[2]) + add!(design, components[3]) + add!(design, neutral_component) + + @test length(design.components) == 4 + @test design.cable_id == "NAYCWY_O_3x95_30x2_5" + @test design.components[1].id == "core1" + @test design.components[2].id == "core2" + @test design.components[3].id == "core3" + @test design.components[4].id == "neutral" + end + + @testset "DataFrame and preview" begin + # Re-create the full design + n_sectors = 3 + r_back_mm = 10.24 + d_sector_mm = 9.14 + r_corner_mm = 1.02 + theta_cond_deg = 119.0 + ins_thick = 1.1e-3 + sector_params = SectorParams(n_sectors, r_back_mm/1000, d_sector_mm/1000, r_corner_mm/1000, theta_cond_deg, ins_thick) + rot_angles = (0.0, 120.0, 240.0) + sectors = [Sector(sector_params, ang, aluminum) for ang in rot_angles] + insulators = [SectorInsulator(sectors[i], ins_thick, pvc) for i in 1:3] + components = [ + CableComponent("core1", ConductorGroup(sectors[1]), InsulatorGroup(insulators[1])), + CableComponent("core2", ConductorGroup(sectors[2]), InsulatorGroup(insulators[2])), + CableComponent("core3", ConductorGroup(sectors[3]), InsulatorGroup(insulators[3])) + ] + n_neutral = 30 + r_strand = 0.79e-3 + R_N = 14.36e-3 + R_O = 17.25e-3 + inner_radius_neutral = R_N - r_strand + outer_jacket_thickness = R_O - (R_N + r_strand) + neutral_wires = WireArray(inner_radius_neutral, Diameter(2*r_strand), n_neutral, 0.0, copper) + neutral_jacket = Insulator(neutral_wires, Thickness(outer_jacket_thickness), pvc) + neutral_component = CableComponent("neutral", ConductorGroup(neutral_wires), InsulatorGroup(neutral_jacket)) + design = CableDesign("NAYCWY_O_3x95_30x2_5", components[1]) + add!(design, components[2]) + add!(design, components[3]) + add!(design, neutral_component) + + # Test that DataFrame constructors do not throw errors + @test DataFrame(design, :detailed) isa DataFrame + @test DataFrame(design, :components) isa DataFrame + @test DataFrame(design, :baseparams) isa DataFrame + + # Test that preview functions execute without error + @test preview(design, display_plot=false) isa Any + end + + @testset "Error handling" begin + # Test invalid geometric parameters for Sector + @test_throws ArgumentError SectorParams(3, -10.24/1000, 9.14/1000, 1.02/1000, 119.0, 1.1e-3) + @test_throws ArgumentError SectorParams(3, 10.24/1000, -9.14/1000, 1.02/1000, 119.0, 1.1e-3) + @test_throws ArgumentError SectorParams(3, 10.24/1000, 9.14/1000, -1.02/1000, 119.0, 1.1e-3) + @test_throws ArgumentError SectorParams(3, 10.24/1000, 9.14/1000, 1.02/1000, 119.0, -1.1e-3) + end +end From 92d49111cde56aa10109be55480e0c2ee01269ba Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 30 Sep 2025 12:34:23 +0200 Subject: [PATCH 28/36] tutorial(sector): added save capability --- examples/tutorial2_sector.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/examples/tutorial2_sector.jl b/examples/tutorial2_sector.jl index b020c1ef..848b5ac8 100644 --- a/examples/tutorial2_sector.jl +++ b/examples/tutorial2_sector.jl @@ -143,6 +143,20 @@ println("Previewing cable design...") plt, _ = preview(design) plt #hide +#= +## Storing in a Library + +Finally, we can store the cable design in a `CablesLibrary` for future reference. +=# + +library = CablesLibrary() +add!(library, design) +library_df = DataFrame(library) + +# Save to file for later use: +output_file = fullfile("cables_library.json") +save(library, file_name = output_file); + #= ## Conclusion @@ -155,3 +169,4 @@ This tutorial has demonstrated how to model a three-core cable with sector-shape This detailed modeling capability allows for accurate analysis of various cable configurations. =# + From c22ab4ecafa797f5b0ae1971224bab3bd29c77a7 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Tue, 30 Sep 2025 17:57:19 +0200 Subject: [PATCH 29/36] feat(sector): add centroid calculation for sector shape --- src/datamodel/sector.jl | 43 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index f38b25c3..136b56cf 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -48,6 +48,8 @@ struct Sector{T<:REALSCALAR} <: AbstractConductorPart{T} gmr::T "Calculated vertices defining the polygon shape." vertices::Vector{Point{2,T}} + "Geometric centroid of the sector shape." + centroid::Point{2,T} end @@ -81,6 +83,10 @@ function Sector( cross_section = _shoelace_area(rotated_vertices) #@debug "Sector cross-sectional area: $(cross_section*1e6) mm²" @debug "Sector cross-sectional area (Shoelace): $(cross_section*1e6) mm²" + + # Calculate centroid + centroid = _calculate_polygon_centroid(rotated_vertices) + @debug "Sector numerically calculated centroid point is: $(centroid)" # 4. Calculate DC resistance rho_eff = calc_temperature_correction(material_props.alpha, temperature, material_props.T0) * material_props.rho resistance = rho_eff / cross_section @@ -100,6 +106,7 @@ function Sector( resistance, gmr, rotated_vertices, + centroid, ) end @@ -272,4 +279,40 @@ function _shoelace_area(vertices::AbstractVector{Point{2,T}}) where {T<:Real} area += (p1[1] * p2[2] - p2[1] * p1[2]) end return abs(area) / T(2) +end + +function _calculate_polygon_centroid(vertices::AbstractVector{Point{2,T}}) where {T<:Real} + n = length(vertices) + if n < 3 + return Point2f(0, 0) + end + + area = _shoelace_area(vertices) + if area < 1e-12 + return Point2f(0, 0) + end + + Cx = zero(T) + Cy = zero(T) + + for i in 1:n + p1 = vertices[i] + p2 = vertices[mod1(i + 1, n)] + cross_prod = (p1[1] * p2[2] - p2[1] * p1[2]) + Cx += (p1[1] + p2[1]) * cross_prod + Cy += (p1[2] + p2[2]) * cross_prod + end + + # The centroid calculation needs to account for the fact that the shoelace formula + # gives a signed area. To get the correct centroid, we must use the signed area + # before taking the absolute value for the final area calculation. + signed_area = zero(T) + for i in 1:n + p1 = vertices[i] + p2 = vertices[mod1(i + 1, n)] + signed_area += (p1[1] * p2[2] - p2[1] * p1[2]) + end + signed_area /= 2 + + return Point2f(Cx / (6 * signed_area), Cy / (6 * signed_area)) end \ No newline at end of file From bfb55085e599ce03aa20d71283d112d4a12c1602 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 1 Oct 2025 01:25:32 +0200 Subject: [PATCH 30/36] fix(sector-geometry): improve centroid calculation --- src/datamodel/sector.jl | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index 136b56cf..39847198 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -287,32 +287,27 @@ function _calculate_polygon_centroid(vertices::AbstractVector{Point{2,T}}) where return Point2f(0, 0) end - area = _shoelace_area(vertices) - if area < 1e-12 - return Point2f(0, 0) - end - + Cx = zero(T) Cy = zero(T) + signed_area = zero(T) for i in 1:n p1 = vertices[i] p2 = vertices[mod1(i + 1, n)] cross_prod = (p1[1] * p2[2] - p2[1] * p1[2]) + signed_area += (p1[1] * p2[2] - p2[1] * p1[2]) Cx += (p1[1] + p2[1]) * cross_prod Cy += (p1[2] + p2[2]) * cross_prod end - - # The centroid calculation needs to account for the fact that the shoelace formula - # gives a signed area. To get the correct centroid, we must use the signed area - # before taking the absolute value for the final area calculation. - signed_area = zero(T) - for i in 1:n - p1 = vertices[i] - p2 = vertices[mod1(i + 1, n)] - signed_area += (p1[1] * p2[2] - p2[1] * p1[2]) - end signed_area /= 2 - return Point2f(Cx / (6 * signed_area), Cy / (6 * signed_area)) + if abs(signed_area) < 1e-12 + return Point2f(0, 0) + end + + Cx /= (6 * signed_area) + Cy /= (6 * signed_area) + + return Point2f(Cx, Cy) end \ No newline at end of file From 14110955e94a80160533e1e85630bf4f06109a27 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Fri, 3 Oct 2025 11:13:27 +0200 Subject: [PATCH 31/36] feat(sector-cableslibrary): enhance handling for Sector and SectorInsulator layer types --- src/importexport/cableslibrary.jl | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/src/importexport/cableslibrary.jl b/src/importexport/cableslibrary.jl index 0ffde1a4..44609dfe 100644 --- a/src/importexport/cableslibrary.jl +++ b/src/importexport/cableslibrary.jl @@ -382,8 +382,28 @@ function _reconstruct_partsgroup(layer_data::Dict) ) elseif LayerType == Semicon radius_ext = get_as(deserialized_layer_dict, :radius_ext, missing, BASE_FLOAT) - ismissing(radius_ext) && Base.error("Missing 'radius_ext' for Semicon first layer.") + ismissing(radius_ext) && Base.error( + "Missing 'radius_ext' for first layer type $LayerType in data: $layer_data", + ) return Semicon(radius_in, radius_ext, material_props; temperature=temperature) + elseif LayerType == Sector + params = get_as(deserialized_layer_dict, :params, missing, BASE_FLOAT) + rotation_angle_deg = get_as(deserialized_layer_dict, :rotation_angle_deg, missing, BASE_FLOAT) + + ismissing(params) && Base.error("Missing 'params' for Sector in data: $layer_data") + !(params isa SectorParams) && error("'params' did not deserialize to a SectorParams object. Got: $(typeof(params))") + ismissing(rotation_angle_deg) && Base.error("Missing 'rotation_angle_deg' for Sector in data: $layer_data") + + return Sector(params, rotation_angle_deg, material_props; temperature=temperature) + elseif LayerType == SectorInsulator + inner_sector = get_as(deserialized_layer_dict, :inner_sector, missing, BASE_FLOAT) + thickness = get_as(deserialized_layer_dict, :thickness, missing, BASE_FLOAT) + + ismissing(inner_sector) && Base.error("Missing 'inner_sector' for SectorInsulator in data: $layer_data") + !(inner_sector isa Sector) && error("'inner_sector' did not deserialize to a Sector object. Got: $(typeof(inner_sector))") + ismissing(thickness) && Base.error("Missing 'thickness' for SectorInsulator in data: $layer_data") + + return SectorInsulator(inner_sector, thickness, material_props; temperature=temperature) else Base.error("Unsupported layer type for first layer reconstruction: $LayerType") end From 4b52082fd601d575d9dbba67f457a133a589bb1f Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Fri, 3 Oct 2025 11:20:30 +0200 Subject: [PATCH 32/36] feat(Sector): export Sector and SectorInsulator in the data model --- src/LineCableModels.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LineCableModels.jl b/src/LineCableModels.jl index 76b55e46..0af446fa 100644 --- a/src/LineCableModels.jl +++ b/src/LineCableModels.jl @@ -9,7 +9,7 @@ export add!, set_verbosity!, set_backend! export Material, MaterialsLibrary # Data model (design + system): -export Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator +export Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator, Sector, SectorInsulator export ConductorGroup, InsulatorGroup export CableComponent, CableDesign, NominalData export CablesLibrary From c0eb865f935863c7fb68bba27a782b473c3b1e0c Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 11 Feb 2026 03:52:51 +0100 Subject: [PATCH 33/36] fix(Sector): fix export of Sector Shape functions --- src/LineCableModels.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LineCableModels.jl b/src/LineCableModels.jl index 0af446fa..a1e0e61e 100644 --- a/src/LineCableModels.jl +++ b/src/LineCableModels.jl @@ -9,7 +9,7 @@ export add!, set_verbosity!, set_backend! export Material, MaterialsLibrary # Data model (design + system): -export Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator, Sector, SectorInsulator +export Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator, Sector, SectorParams, SectorInsulator export ConductorGroup, InsulatorGroup export CableComponent, CableDesign, NominalData export CablesLibrary @@ -61,7 +61,7 @@ using .EarthProps: EarthModel # Submodule `DataModel` include("datamodel/DataModel.jl") -using .DataModel: Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator, +using .DataModel: Thickness, Diameter, WireArray, Strip, Tubular, Semicon, Insulator, Sector, SectorParams, SectorInsulator, ConductorGroup, InsulatorGroup, CableComponent, CableDesign, NominalData, CablesLibrary, CablePosition, LineCableSystem, trifoil_formation, flat_formation, preview, equivalent From d736ca107cea0e423d329f979139a7038b7380c3 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 11 Feb 2026 04:07:42 +0100 Subject: [PATCH 34/36] fix(Sector): correct include statement for SectorInsulator file --- src/datamodel/DataModel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/datamodel/DataModel.jl b/src/datamodel/DataModel.jl index db06a782..529670a9 100644 --- a/src/datamodel/DataModel.jl +++ b/src/datamodel/DataModel.jl @@ -77,7 +77,7 @@ include("sector.jl") include("insulator.jl") include("semicon.jl") include("insulatorgroup.jl") -include("SectorInsulator.jl") +include("sectorinsulator.jl") # Groups From a6d4a0e8017fd6382d9087ca7e10e8526ae6d86c Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 11 Feb 2026 04:49:31 +0100 Subject: [PATCH 35/36] feat(SectorParams): basic validation and constructor parameters for SectorParams struct. just testing validation. --- src/datamodel/sector.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index 39847198..6341cdfa 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -18,8 +18,30 @@ struct SectorParams{T<:REALSCALAR} theta_cond_deg::T "Insulation thickness \\[m\\]." d_insulation::T # needed to correct for the offset + + function SectorParams(n_sectors::Int, r_back::T, d_sector::T, r_corner::T, theta_cond_deg::T, d_insulation::T) where {T<:REALSCALAR} + if r_back < 0 + throw(ArgumentError("Sector geometry: r_back must be non-negative.")) + end + if d_sector < 0 + throw(ArgumentError("Sector geometry: d_sector must be non-negative.")) + end + if r_corner < 0 + throw(ArgumentError("Sector geometry: r_corner must be non-negative.")) + end + if d_insulation < 0 + throw(ArgumentError("Sector geometry: d_insulation must be non-negative.")) + end + + new{T}(n_sectors, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) + end end +SectorParams(n_sectors::Int, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) = + let T = promote_type(typeof(r_back), typeof(d_sector), typeof(r_corner), typeof(theta_cond_deg), typeof(d_insulation)) + SectorParams(n_sectors, convert(T, r_back), convert(T, d_sector), convert(T, r_corner), convert(T, theta_cond_deg), convert(T, d_insulation)) + end + """ $(TYPEDEF) From 5293c2a777346810b0744e6f69de8c122f3ecc42 Mon Sep 17 00:00:00 2001 From: MohamedNumair Date: Wed, 11 Feb 2026 05:00:05 +0100 Subject: [PATCH 36/36] feat(SectorParams): ADDED Validation for sector geometry. I followed the rules so the sectors play nice with other cable parts, In hope of being accepted into the main branch. --- src/datamodel/sector.jl | 96 +++++++++++++++++++++++++++++++++-------- 1 file changed, 79 insertions(+), 17 deletions(-) diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl index 6341cdfa..d7719ca3 100644 --- a/src/datamodel/sector.jl +++ b/src/datamodel/sector.jl @@ -20,27 +20,89 @@ struct SectorParams{T<:REALSCALAR} d_insulation::T # needed to correct for the offset function SectorParams(n_sectors::Int, r_back::T, d_sector::T, r_corner::T, theta_cond_deg::T, d_insulation::T) where {T<:REALSCALAR} - if r_back < 0 - throw(ArgumentError("Sector geometry: r_back must be non-negative.")) - end - if d_sector < 0 - throw(ArgumentError("Sector geometry: d_sector must be non-negative.")) - end - if r_corner < 0 - throw(ArgumentError("Sector geometry: r_corner must be non-negative.")) - end - if d_insulation < 0 - throw(ArgumentError("Sector geometry: d_insulation must be non-negative.")) - end - - new{T}(n_sectors, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) + # validation for the sector geometry constraints (angle limits, no overlap, discriminant check) + _validate_sector_params(n_sectors, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) + new{T}(n_sectors, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) end end -SectorParams(n_sectors::Int, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) = - let T = promote_type(typeof(r_back), typeof(d_sector), typeof(r_corner), typeof(theta_cond_deg), typeof(d_insulation)) - SectorParams(n_sectors, convert(T, r_back), convert(T, d_sector), convert(T, r_corner), convert(T, theta_cond_deg), convert(T, d_insulation)) +const _REQ_SECTOR_PARAMS = (:n_sectors, :r_back, :d_sector, :r_corner, :theta_cond_deg, :d_insulation) + +Validation.has_radii(::Type{SectorParams}) = false +Validation.required_fields(::Type{SectorParams}) = _REQ_SECTOR_PARAMS +Validation.coercive_fields(::Type{SectorParams}) = (:r_back, :d_sector, :r_corner, :theta_cond_deg, :d_insulation) + +# --- Geometry Logic Storage --- +function _validate_sector_params(n_sectors, r_back, d_sector, r_corner, theta_cond_deg, d_insulation) + # 1. Basic Non-negativity (redundant with rules but kept for Core safety) + if r_back < 0 || d_sector < 0 || r_corner < 0 || theta_cond_deg < 0 || d_insulation < 0 + throw(ArgumentError("Sector geometry parameters must be non-negative.")) + end + + # 2. Basic constraints + if theta_cond_deg >= 360.0 + throw(ArgumentError("[SectorParams] theta_cond_deg ($theta_cond_deg) must be < 360")) end + allowed_angle = 360.0 / n_sectors + if theta_cond_deg > allowed_angle + 1e-4 + throw(ArgumentError("[SectorParams] theta_cond_deg ($theta_cond_deg) exceeds allowed 360/n ($allowed_angle)")) + end + + if d_sector >= r_back + throw(ArgumentError("[SectorParams] d_sector ($d_sector) must be less than r_back ($r_back)")) + end + + # 3. Geometric feasibility (Discriminant check) + phi_deg = 90.0 - theta_cond_deg / 2.0 + phi_rad = deg2rad(phi_deg) + + if abs(cos(phi_rad)) < 1e-9 + throw(ArgumentError("[SectorParams] theta_cond_deg too close to 180 (phi ~ 90), valid sector cannot be formed.")) + end + + d_base_corner = r_corner * (1.0 / cos(phi_rad) - 1.0) + d_offset = r_back - d_sector - d_base_corner + + k = r_corner / cos(phi_rad) + d_offset + qa = 1.0 + tan(phi_rad)^2 + qb = 2.0 * k * tan(phi_rad) + qc = k^2 - (r_back - r_corner)^2 + + discriminant = qb^2 - 4.0 * qa * qc + if discriminant < 0 + throw(ArgumentError("[SectorParams] Geometric constraints violated (negative discriminant). The combination of radii, depth and angle does not form a closed sector.")) + end +end + +""" +$(TYPEDEF) + +Rule that enforces specific Sector geometry constraints (discriminant, angle limits, no overlap). +""" +struct SectorGeometryValid <: Validation.Rule + name::Symbol +end + +function Validation._apply(r::SectorGeometryValid, nt, ::Type{SectorParams}) + # Delegate to the shared validation function + # Note: Non-neg rules from Validation framework handle the basic checks, + # but _validate_sector_params repeats them safely. + _validate_sector_params(nt.n_sectors, nt.r_back, nt.d_sector, nt.r_corner, nt.theta_cond_deg, nt.d_insulation) +end + +Validation.extra_rules(::Type{SectorParams}) = ( + Validation.IntegerField(:n_sectors), + Validation.Positive(:n_sectors), + Validation.Nonneg(:r_back), + Validation.Nonneg(:d_sector), + Validation.Nonneg(:r_corner), + Validation.Nonneg(:theta_cond_deg), + Validation.Nonneg(:d_insulation), + Validation.Less(:d_sector, :r_back), + SectorGeometryValid(:sector_geometry), +) + +@construct SectorParams _REQ_SECTOR_PARAMS """ $(TYPEDEF)