diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4cd77cac..23ce0a83 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,8 +1,8 @@ name: CI on: push: - branches: - - main + branches-ignore: + - gh-pages tags: ['*'] pull_request: workflow_dispatch: diff --git a/Project.toml b/Project.toml index 0ad46d40..def66611 100644 --- a/Project.toml +++ b/Project.toml @@ -28,6 +28,7 @@ Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" 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" @@ -68,6 +69,7 @@ Measurements = "2.11.0" NLsolve = "4.5.1" Pkg = "1.11.0" Plots = "1.40.9" +PolygonOps = "0.1.2" Printf = "1.11.0" QuadGK = "2.11.2" Reexport = "1.2.2" diff --git a/examples/tutorial2_sector.jl b/examples/tutorial2_sector.jl new file mode 100644 index 00000000..848b5ac8 --- /dev/null +++ b/examples/tutorial2_sector.jl @@ -0,0 +1,172 @@ +#= +# 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 + +#= +## 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 + +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/src/LineCableModels.jl b/src/LineCableModels.jl index 76b55e46..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 +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 diff --git a/src/datamodel/DataModel.jl b/src/datamodel/DataModel.jl index 07ae4793..529670a9 100644 --- a/src/datamodel/DataModel.jl +++ b/src/datamodel/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 @@ -52,7 +52,8 @@ using DataFrames using Colors using Plots using DisplayAs: DisplayAs - +using LinearAlgebra +using Makie: Point, Point2f # otherwise will require adding GeometryBasics as a dependency # Abstract types & interfaces include("types.jl") include("radii.jl") @@ -70,11 +71,13 @@ include("wirearray.jl") include("strip.jl") include("tubular.jl") include("conductorgroup.jl") +include("sector.jl") # Insulators include("insulator.jl") include("semicon.jl") include("insulatorgroup.jl") +include("sectorinsulator.jl") # Groups diff --git a/src/datamodel/preview.jl b/src/datamodel/preview.jl index 7e8d504a..55afd4c8 100644 --- a/src/datamodel/preview.jl +++ b/src/datamodel/preview.jl @@ -369,6 +369,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 diff --git a/src/datamodel/sector.jl b/src/datamodel/sector.jl new file mode 100644 index 00000000..d7719ca3 --- /dev/null +++ b/src/datamodel/sector.jl @@ -0,0 +1,397 @@ +""" +$(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 + "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} + # 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 + +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) + +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}} + "Geometric centroid of the sector shape." + centroid::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] + + # 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) + 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 + + # 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, + centroid, + ) +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) # ϕ + @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_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) # X_A = -X_F + y_base_corner = x_base_corner * tan(phi_rad) + d_offset # Y_A = Y_F + 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 + 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) # 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 = 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 # X_C + y_side_upper = y_side_center * p.r_back / dist_origin_side_center # Y_C + node_C = 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]) + ) + @debug "(Sector) Nodes: $nodes" + centers = ( + 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) +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) + @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 + +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 + + 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) + @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 + 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) + @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) + 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 + + +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::T = zero(T) + for i in 1:n + p1 = vertices[i] + p2 = vertices[mod1(i + 1, n)] + 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 + + + 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 + signed_area /= 2 + + 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 diff --git a/src/datamodel/sectorinsulator.jl b/src/datamodel/sectorinsulator.jl new file mode 100644 index 00000000..d1938af7 --- /dev/null +++ b/src/datamodel/sectorinsulator.jl @@ -0,0 +1,131 @@ +""" +$(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.vertices, thickness) + + # 2. Calculate areas + inner_area = inner_sector.cross_section + + # 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 + # 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) --- +# 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 + + return new_vertices +end \ No newline at end of file diff --git a/src/engine/fem/FEM.jl b/src/engine/fem/FEM.jl index 5759a593..4522c1ee 100644 --- a/src/engine/fem/FEM.jl +++ b/src/engine/fem/FEM.jl @@ -43,7 +43,7 @@ using ...Utils: using Measurements using LinearAlgebra using Colors - +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 29753b81..73ec50e7 100644 --- a/src/engine/fem/cable.jl +++ b/src/engine/fem/cable.jl @@ -448,3 +448,96 @@ function _make_cablepart!(workspace::FEMWorkspace, part::WireArray, 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 = [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, mesh_size) + + # --- 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 = [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 + 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, mesh_size) + + # --- 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 diff --git a/src/engine/fem/drawing.jl b/src/engine/fem/drawing.jl index 669b2f5c..bfab73c4 100644 --- a/src/engine/fem/drawing.jl +++ b/src/engine/fem/drawing.jl @@ -640,6 +640,98 @@ function draw_transition_region( return tags, all_mesh_points, markers end +function draw_polygon(vertices::Vector{<:Point}) + @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] + + # 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 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 + +""" +$(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 _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 + +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] + 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 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) + + # 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) @@ -737,4 +829,35 @@ 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 = _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] + + # 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 + 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 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